Cantilever beam vibration modes

General Discussion
1

Hello,

I am trying to compute the vibration modes of a cantilever beam in 3D for different material densities.
It can be summarized as solving an eigenvalue problem using stiffness and mass matrices.
Here is the code I use:

load "msh3"
load "medit"
include "MeshSurface.idp"
load "TetGen"
load "iovtk"

// Parameters
real E0 = 1e-9; //void Youngs modulus
real E1 = 1; //material Youngs modulus
real nu = 0.3; //Poisson's modulus
real W0 = 1e-5; //void mass
real W1 = 1; //material mass
real lambda = nu/((1+nu)*(1-2*nu)); //Lamé coefficient
real mu = 1/(2*(1+nu)); //Lamé coefficient
real p = 3; // penalization exponent

// Boundary conditions
int free = 0, loaded = 1, fixed = 2;
// Simple full beam
int[int] BeamNel = [30,3,3];
real [int,int] BeamBounds = [[-5, 5], [-0.5, 0.5], [-0.5, 0.5]];
int [int,int] BeamLabels = [[2, 1], [0, 0], [0, 0]];
meshS ThSBeam = SurfaceHex(BeamNel, BeamBounds, BeamLabels, 1);
mesh3 Th = tetg(ThSBeam, switch="pqaAAYYQ");

// Finite element spaces
fespace Vh (Th, [P1, P1, P1]);
fespace Vh0 (Th, P0);

Vh [ux, uy, uz], [vx, vy, vz], [ax, ay, az];
Vh0 theta, E, W, V;

// Element densities
theta=1;

// Element volumes
V = volume;

// 3D Element stiffness matrix macro
macro D() ([[2*mu + lambda, lambda, lambda, 0, 0, 0],
            [lambda, 2*mu + lambda, lambda, 0, 0, 0],
            [lambda, lambda, 2*mu + lambda, 0, 0, 0],
            [0, 0, 0, mu, 0, 0],
            [0, 0, 0, 0, mu, 0],
            [0, 0, 0, 0, 0, mu]]) // EOM
// Strain macro
macro e (x1,x2,x3) [dx(x1),dy(x2),dz(x3),(dz(x2)+dy(x3)),(dz(x1)+dx(x3)),(dy(x1)+dx(x2))] // EOM

/* Stiffness */
// Variational form of the beam stiffness matrix
varf varfK([ux,uy,uz], [vx,vy,vz]) = int3d(Th) (E*((e(ux,uy,uz)'*D*e(vx,vy,vz))))
                                + on(fixed, ux=0, uy=0, uz=0);
// Elastic modulus penalization formula
E = E0 + (E1-E0)*theta^p;
// Stiffness matrix
matrix K = varfK(Vh, Vh);

/* Mass */
// Variational form of the beam mass matrix
varf varfM([ax, ay, az], [vx, vy, vz]) = int3d(Th) (V*W*[ax, ay, az]'*[vx, vy, vz]);
// Mass penalization formula
W = W0 + (W1-W0)*theta^p;
// Mass matrix
matrix M = varfM(Vh, Vh);

// Eigenvalue setup
int neval = 5;
real[int] eval(neval);
Vh0[int] evec(neval);

// Eigenvalue solving
int k=EigenValue(K, M,nev=neval, sigma=20, value=eval, vector=evec, sym=1);

An error comes from the last line:

current line = 72
Assertion fail : (SameShape(y, *xx))
        line :626, in file lgfem.hpp
Assertion fail : (SameShape(y, *xx))
        line :626, in file lgfem.hpp
 err code 6 ,  mpirank 0

I was not able to find an example in the documentation that was close enough to my problem for me to understand where the error comes from.
Are the stiffness and mass matrices correctly defined ?
Is the EigenValue function call wrong ?
Thank you in advance for your time and answer.
Hugo

2

You are mixing objects of different dimensions (Vh.ndof vs. Vh0.ndof).

3

Thank you for your answer,

However I’m not sure how to modify the code to take your comment into account.
I have tried the following :

  • Setting Vh0 to P1 elements (This way all objects are P1)
  • Removing E, W and V from varfK and varfM (Only Vh P1 objects in the varf definition)
    What could I change to make it work ?
    Regards,
4

Change the definition of evec.

5

Does that mean evec should be an array of size neval, for which each entry is itself an array of 3 vector fields ?
Like the following :

[ [Vh0, Vh0, Vh0], [Vh0, Vh0, Vh0], ...]

Anyhow, removing vector=evec from the last line removes the error.
Thank you !

6

I don’t understand what you are saying, again, you should just change the definition of evec, and stick to Vh.