Hi there,

I am an absolute novice to FREEFEM++, and I cannot figure out how to (efficiently) obtain the stiffness matrix.

I found in the documentation reference to the following code

varf a(u,v)

= int2d(Th)(

dx(u)*dx(v)

+ dy(u)*dy(v)

)

+ on(C, u=0)

;

matrix A = a(Vh, Vh); //stiffness matrix

also referenced here https://people.sc.fsu.edu/~jburkardt/freefem_src/stiffnessmatrix/stiffnessmatrix.html

but I cannot make it work for an elasticity bilinear operator (I have started playing with the linear elasticity example provided in the installation).

// Macro

real sqrt2=sqrt(2.);

macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] //

// The sqrt2 is because we want: epsilon(u1,u2)’* epsilon(v1,v2) = epsilon(u): epsilon(v)

macro div(u,v) ( dx(u)+dy(v) ) //

// Problem

real mu= E/(2*(1+nu));

real lambda = E*nu/((1+nu)*(1-2*nu));

solve lame([u, v], [uu, vv])

= int2d(Th)(

lambda * div(u, v) * div(uu, vv)

+ 2.*mu * ( epsilon(u,v)’ * epsilon(uu,vv) )

)

- int2d(Th)(

f*vv*vv

)

-int1d(Th, 2) (

shear

)

- on(4, u=0, v=0)

;

I tried issuing something as

//Stiffness matrix

varf bid ( u, v, uu, vv ) = int2d(Th)(

lambda * div(u, v) * div(uu, vv)

+ 2.*mu * ( epsilon(u,v)’ * epsilon(uu,vv) )

)

- on(4, u=0, v=0)

;

matrix ad = bid ( Vh, Vh);

but i get an error “Check Bilinear Operator”. I suspect something is wrong with the arguments, or the number thereof, but for starters I am not so sure can a function be applied to general finite element functions.

In the scalar example reported above, Vh is the finite element space, wjat does exactly mean to apply a function to an argument as Vh? Is it a shortcut to get the related trial function as an argument?

So how should the vectorial case be treated?

I fear I am missing the obvious, sorry if that is the case.

Any help would be most appreciated, thanks a lot