# Print out the Stiffness Matrix for a vectorial variational problem

**URL:** <https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187>\
**Category:** General Discussion\
**Created:** [December 1, 2019, 11:52pm UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187 "2019-12-01T23:52:31Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Smerdjakov](https://avatars.discourse-cdn.com/v4/letter/s/f4b2a3/32.png) [@Smerdjakov](https://community.freefem.org/u/Smerdjakov)\
**Post date:** [December 1, 2019, 11:52pm UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187/1 "2019-12-01T23:52:31Z")

</div>

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](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  
)  
-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

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [December 2, 2019, 7:14am UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187/2 "2019-12-02T07:14:31Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Smerdjakov](https://avatars.discourse-cdn.com/v4/letter/s/f4b2a3/32.png) [@Smerdjakov](https://community.freefem.org/u/Smerdjakov)\
**Post date:** [December 2, 2019, 7:45pm UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187/3 "2019-12-02T19:45:33Z")

</div>

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](https://doc.freefem.org/tutorials/elasticity.html)

Thanks again a lot.

---

<div class="post-metadata">

**Author:** ![prj](https://avatars.discourse-cdn.com/v4/letter/p/ecae2f/32.png) [@prj](https://community.freefem.org/u/prj)\
**Post date:** [December 2, 2019, 8:22pm UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187/4 "2019-12-02T20:22:59Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Smerdjakov](https://avatars.discourse-cdn.com/v4/letter/s/f4b2a3/32.png) [@Smerdjakov](https://community.freefem.org/u/Smerdjakov)\
**Post date:** [December 2, 2019, 8:56pm UTC](https://community.freefem.org/t/print-out-the-stiffness-matrix-for-a-vectorial-variational-problem/187/5 "2019-12-02T20:56:35Z")

</div>

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.