# Solve U=K^-1.F using stiffness matrix for vector fields

**URL:** https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604
**Category:** General Discussion
**Created:** [October 6, 2020, 6:39pm UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604 "2020-10-06T18:39:00Z")
**Posts on this page:** 13
**Page:** 1

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 6, 2020, 6:39pm UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/1 "2020-10-06T18:39:01Z")

</div>

Dear all,

I am trying to solve problem of 3-dimensional elasticity. I have extracted the stiffness matrix (K) and force vector (F) using the variational formulation. Now, I have solved the equation for displacement U, using U=K^-1.F.  
Now, displacement is a vector field, so I need to again distribute U into its respective forms, U1, U2, U3. I did it but the plot is not what I should get for a tip-loaded beam.  
Please find my code below. Thanks in advance.

load “mshmet”

load “msh3”

load “medit”

load “TetGen”

include “cube.idp”

load “lapack”

real l = 10; // lenght

real w = 1; // breadth

real h = 1; // height

real I = w_h_h\*h/12; // area moment of inertia

int[int] Nxyz=[20,2,2];

real [int,int] Bxyz=[[0.,l],[0.,w],[0.,h]];

int [int,int] Lxyz=[[1,2],[3,4],[5,6]];

mesh3 Th=Cube(Nxyz,Bxyz,Lxyz);

real E = 210.0e9;

real nu = 0.3;

real rho = 7850;

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

real lambda = E_nu/((1+nu)_(1-2\*nu));

// Fespace

fespace Vh(Th,[P1,P1,P1]);

Vh [u1,u2,u3], [v1,v2,v3];

real sqrt2=sqrt(2.);

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

macro div(u1,u2,u3) ( dx(u1)+dy(u2)+dz(u3) ) // EOM

varf a([u1,u2,u3],[v1,v2,v3],tgv=1e50)=

```
int3d(Th)(  

    lambda*div(u1,u2,u3)*div(v1,v2,v3)  

    +2.*mu*( epsilon(u1,u2,u3)'*epsilon(v1,v2,v3) )  

         )

+ on(1,u1=0,u2=0,u3=0)

```

;

varf f([u1,u2,u3],[v1,v2,v3])=

```
int2d(Th,2)(60000*v3);

```

matrix A= a(Vh,Vh,solver=“SPARSESOLVER”); //stiffness matrix

real[int] F = f(0,Vh);

real[int,int] K(A.m,A.m);

real c = 0;

for (int i=0; i\<A.m;i++){

```
for(int j=0; j<A.m; j++){

    K(i,j)=A(i,j);

}

```

}

real[int,int] a1realinv(A.m,A.n);

a1realinv = K^-1;

matrix Kinv(A.m,A.n); /_inverse of matrix a1_/

for (int e=0; e\<A.m; e++){

```
for (int ee=0;ee<A.n; ee++){

    Kinv(e,ee) = a1realinv(e,ee);

}

```

}

real[int] U(A.n); /_Displacement_/

U=Kinv\*F; [//U=K^-1.F](https://U=K%5E-1.F)

real[int] U1(A.m/3),U2(A.m/3),U3(A.m/3);

for(int e=0; e\<A.m/3; e++){

```
U1(e) = U(3*e);

U2(e) = U(3*e+1);

U3(e) = U(3*e+2);

```

}

real coef = 1000;

int[int] ref2 = [1, 0, 2, 0];

Vh [uu,vv,ww];

uu[] = U1;

vv[] = U2;

ww[] = U3;

mesh3 Thm = movemesh3(Th, transfo=[x+uu_coef, y+vv_coef, z+ww\*coef], label=ref2);

Thm = change(Thm, label=ref2);

plot(Thm, wait=true, cmm="coef amplification = "+coef);

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/ea0db24ba546597e38b0eb8973ac21f6953f147b.png)

---

<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: [October 7, 2020, 9:14am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/2 "2020-10-07T09:14:43Z")

</div>

Hello,  
What are you trying to do exactly? Computing `K^-1` explicitly should be avoided at all cost.

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 7, 2020, 9:30am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/3 "2020-10-07T09:30:42Z")

</div>

Hi prj,

My final goal is to use the Newmark-beta method for transient problems in elasticity.  
So in order to do that, I need the Mass and Stiffness matrices for the body and operate on the same.  
The equations to be solved consist of the inverse of these matrices at final.  
To get the idea, I was first trying with this static deflection problem of beam.  
So how should I do that in case K^-1 should not be used?  
I am attaching the image of one of the equation.  
[m]-mass matrix, [k] - stiffness matrix

![newmark](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/789dbd8cd7c7aac698c376ff97e7849c3dbc4f7d.png)

---

<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: [October 9, 2020, 10:16am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/4 "2020-10-09T10:16:03Z")

</div>

In your code you have:

```auto
U=Kinv*F; //U=K^-1.F

```

What do you do that, and not:

```auto
U=K^-1.F

```

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 9, 2020, 10:50am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/5 "2020-10-09T10:50:27Z")

</div>

In my code,  
**Kinv** is an array real[int, int] which contains the elements of the inverse of stiffness matrix.  
**F** is an array real[int] which contains the force vector for the whole body.  
**U** is an array containing the displacements.  
To carry out the matrix-vector multiplication, which is compatible in the way I use, I have used Kinv. So further, K_U=F, which gives U = K^-1_F = Kinv\*F.

Also, can you let me know what problems can occur when I take an inverse of the matrix obtained from the **varf**.  
The method I follow is:  
I store the elements of obtained matrix in an array real[int,int] and then use _lapack_ to get the inverse.

---

<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: [October 9, 2020, 10:53am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/6 "2020-10-09T10:53:40Z")

</div>

But why do you need the inverse explicitly, and why can’t you just compute `U=A^-1*F`, using `matrix A` assembled from your `varf`?

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 9, 2020, 11:03am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/7 "2020-10-09T11:03:27Z")

</div>

For the Newmark method that I mentioned earlier, I will require the inverse of addition of matrices explicitly, that is:  
[a_ **M** + b_ **C** + **K**]^-1  
where,  
**M** - Mass matrix for the body  
**C** - Damping Matrix  
**K** - Stiffness matrix of the body  
a,b,c - real-constants  
For doing this, I am not able to create a form to be fed to **varf**.  
This is the reason I was searching for an explicit method. It is a time-stepping method and I need to loop it.  
Can you let me know any other way for this method?

---

<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: [October 9, 2020, 11:10am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/8 "2020-10-09T11:10:57Z")

</div>

What is your damping matrix? I don’t see the problem assembling through a `varf` the term `aM + K`.  
You can then compute the sum of both matrix, and use `set(sum, solver = sparsesolver)`.

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 9, 2020, 11:19am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/9 "2020-10-09T11:19:04Z")

</div>

At present, the damping matrix is zero, but I will be considering a model as follows for damping matrix:  
**C** = n **M** + m **K** ;  
where, M,C,K represent the same matrices as earlier and m,n are real constants,  
From your reply, as far as I understand, I should compute the **K** and **M** matrices using **varf** , then store their addition as required in a matrix **sum** = a_ **M** + b_ **C** + **K** , and then use  
set(sum, solver = sparsesolver)  
Am I correct?

---

<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: [October 9, 2020, 2:58pm UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/10 "2020-10-09T14:58:59Z")

</div>

You are indeed correct.

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [October 10, 2020, 4:15pm UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/11 "2020-10-10T16:15:51Z")

</div>

I was able to execute the whole code much faster using the above way.  
Also, I rectified an error regarding the FE array size in my earlier code.  
Now my code seems to work completely as expected.  
Thanks a lot, @prj!

---

<div class="post-metadata">

### Author: ![distractor](https://avatars.discourse-cdn.com/v4/letter/d/57b2e6/32.png) [@distractor](https://community.freefem.org/u/distractor)
#### Post date: [August 20, 2021, 9:35am UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/12 "2021-08-20T09:35:42Z")

</div>

What are `uucoef` and `vvcoef` in the OP?

---

<div class="post-metadata">

### Author: ![pdapte1999](https://avatars.discourse-cdn.com/v4/letter/p/d9b06d/32.png) [@pdapte1999](https://community.freefem.org/u/pdapte1999)
#### Post date: [September 19, 2021, 7:49pm UTC](https://community.freefem.org/t/solve-u-k-1-f-using-stiffness-matrix-for-vector-fields/604/13 "2021-09-19T19:49:53Z")

</div>

Those are the coefficients for amplification in visualization, they do not affect the results but help in visualization by amplifying the displacements.