### How to assemble separately the mass and stiffness matrix

360

views

0

Dear all,

I have been working with some nonlinear equations and wanted to try to implement a predictor corrector scheme, however in order to do so I would need to access both the mass and the stiffness matrix arrays after I assemble them.

I was playing a little bit with FEniCS and I think that I am stuck.

I provide the following MWE to show what I mean:

I expected that I would have the same results regardless the way I set up the problem, however the output of such code is:

All the best, Murilo.

I have been working with some nonlinear equations and wanted to try to implement a predictor corrector scheme, however in order to do so I would need to access both the mass and the stiffness matrix arrays after I assemble them.

I was playing a little bit with FEniCS and I think that I am stuck.

I provide the following MWE to show what I mean:

```
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
import sympy as sy
from sympy import init_printing
init_printing()
%matplotlib inline
dt = 0.1
mesh = IntervalMesh(2, 0, 1)
P1 = FiniteElement('CG', interval, 2)
V = FunctionSpace(mesh, MixedElement([P1, P1]))
v_2, v_1 = TestFunctions(V)
T, P = TrialFunctions(V)
bcl_T = 'near(x[0],0)'
bc = DirichletBC(V.sub(0), 300, bcl_T)
T_0 = Constant('30')
P_0 = Constant('2900')
T_n = interpolate(T_0, V.sub(0).collapse())
P_n = interpolate(P_0, V.sub(1).collapse())
T_en = Constant('30')
P_en = Constant('2900')
F = (P + T)*v_1*dx - (P_n + T_n)*v_1*dx - dt*dot(grad(P), grad(v_1))*dx - \
dt*(P - P_en)*v_1 * ds(1) + \
(P + T)*v_2*dx - (P_n + T_n)*v_2*dx - dt*0.1*dot(grad(T), grad(v_2))*dx - \
dt*inner(grad(T), grad(P))*dx - dt*(T - T_en)*v_2*ds(2)
a_u_v = lhs(F)
A = assemble(a_u_v)
L_v = rhs(F)
b = assemble(L_v)
bc.apply(A)
bc.apply(b)
# Separeted Version
M_wf = (P + T)*v_1*dx + (P + T)*v_2*dx - dt*P*v_1*ds(1) - dt*T*v_2*ds(2)
K_wf = - dt*dot(grad(P), grad(v_1))*dx - 0.1*dt*dot(grad(T), grad(v_2))*dx - \
dt*inner(grad(T), grad(P))*dx - dt*P*v_1*ds(1) - dt*T*v_2*ds(2)
F_wf = - (P_n + T_n)*v_1*dx - (P_n + T_n)*v_2*dx + dt*P_en*v_1*ds(1) + dt*T_en*v_2*ds(2)
M = assemble(lhs(M_wf))
K = assemble(lhs(K_wf))
F_ = assemble(rhs(F_wf))
bc.apply(M+K)
bc.apply(F_)
u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve((M+K), w.vector(), F_)
u.compute_vertex_values(),w.compute_vertex_values()
```

I expected that I would have the same results regardless the way I set up the problem, however the output of such code is:

```
(array([ 300., 300., 300., 2630., 2630., 2630.]),
array([ -4.31058821e+17, -4.31058821e+17, -4.31058821e+17,
4.31058821e+17, 4.31058821e+17, 4.31058821e+17]))
```

I suspect that I am missing something, if anyone could please point me a direction I would appretiate a lot.All the best, Murilo.

Community: FEniCS Project

### 1 Answer

3

Change the last few lines to this

```
M = assemble(lhs(M_wf)+ lhs(K_wf))
F_ = assemble(rhs(F_wf))
bc.apply(M)
bc.apply(F_)
u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve(M, w.vector(), F_)
print u.compute_vertex_values()
print w.compute_vertex_values()
```

and it seems to work

```
fenics@eb194e04ecb5:~/shared$ python test.py
[ 300. 300. 300. 2630. 2630. 2630.]
[ 300. 300. 300. 2630. 2630. 2630.]
```

or like this also works

```
M = assemble(lhs(M_wf))
K = assemble(lhs(K_wf))
F_ = assemble(rhs(F_wf))
A1 = M + K
bc.apply(A1)
bc.apply(F_)
u = Function(V)
w = Function(V)
solve(A, u.vector(), b)
solve(A1, w.vector(), F_)
print u.compute_vertex_values()
print w.compute_vertex_values()
```

1

I added a second way with M and K. Does that not work for you ?

written
8 months ago by
Praveen C

Wow! It works on my simple MWE! Thank you very much. However, for my more complex problem it gives still different results.

I will need to test a little bit more.

Thank you very much for your help!

All the best, Murilo.

I will need to test a little bit more.

Thank you very much for your help!

All the best, Murilo.

written
8 months ago by
Murilo Moreira

And do you have an explanation of why it works this way and not the others?

Thank you.

Thank you.

written
8 months ago by
Murilo Moreira

1

I think if you do

`bc.apply(M+K)`

it creates temporary matrix containing M+K on which bc is applied. Your original matrices are not changed.
written
8 months ago by
Praveen C

It makes sense! Thank you very much!

written
8 months ago by
Murilo Moreira

Please login to add an answer/comment or follow this question.

Thank you very much for you answer.

It works indeed, however, I am sorry for not being so clear about what I would like to have.

I wanted to have the two matrices, for the mass and stiffness, separately as I would need to perform the calculations on the matrix themselves.

I was only solving for the vector to check if I would have the same outcome if I assembled them separately or if I assembled in the more usual way.

With your given answer the two solutions are the same, however I still have only one matrix with both contributions from the mass and stiffness matrices. I would like to have two matrix. Do you think this would be possible?

Anyway thank you very much for your kind help!

All the best, Murilo.