Print out the Stiffness Matrix for a vectorial variational problem

General Discussion
1

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 = Enu/((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)(
    fvv
    )
    -int1d(Th, 2) (
    shear    vv
    )
  • 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

2

Could you try to replace varf bid ( u, v, uu, vv ) by varf bid ( [u, v], [uu, vv] ), please?

3

Thanks for your reply.

I did try your suggestion, but I got the same error.

It says to “Check Bilinear Operator”, but it is the same used to solve the problem, modulo a boundary condition (I am interested on stiffness matrix with Dirichlet only boundary conditions).

I would have attached the a txt version of the modified input fil or the .edp file, but new users cannot attach files. My script is anyhow following very closely the elasticity example in the documentation, https://doc.freefem.org/tutorials/elasticity.html

Thanks again a lot.

4

Oh, with respect to that example, you should also replace fespace Vh(Th, P2); by fespace Vh(Th, [P2, P2]);.

1 Like
5

At the moment the FE space is defined, as in the example, as

// Fespace
fespace Vh(Th, P2);
Vh u, v;
Vh uu, vv;

If I follow your suggestion, and I issue

fspace Vh(Th, [P2,P2]);

I get the error

“Vh u, the array size must be 2 not 1
Invalid array size for vectorial fespace function”

I tried then to modify the line after the one you hinted at, issuing

fspace Vh(Th, [P2,P2]);
Vh [u, v];
Vh [uu, vv];

and I do get something that looks like a stiffness matrix.

Thanks a lot, I will now try to understand the difference between
fespace Vh(Th, P2)

and

fespace Vh(Th, [P2, P2])

Both formulation yield the same results, until the stiffness matrix is required.
I assume in the first case 4 coupled scalar fields are considered, instead of a vectorial formulation.

Well, thank you again a lot, your support was much appreciated, have a good evening.