### Explicit assembly and external solve

326

views

1

Hi,

For the validation purpose, I am trying to export FEniCS assembly (Poisson's equation) matrix and vector to Matlab and solve using (u = A\b). However, when I solve externally the solution vector looks different (same values but different sparsity) compared to the solution vector extracted from internal FEniCS solve(A, U, b).

For example, take a look at the following Snippet:

A = assemble(a)

b = assemble(L)

bc.apply(A, b)

sio.savemat('Ab.mat', {'A':A.array(), 'b':b.array()})

u = Function(V)

U = u.vector()

solve(A, U, b)

U = [0

0

0

0

0

0

0

0

-0.479166666666669

-0.276041666666667

-0.244791666666667

-0.244791666666667

-0.276041666666667]

Matlab solve (u = A\b)

u=[0

0

-0.244791666666667

0

-0.479166666666669

-0.276041666666667

0

-0.244791666666667

0

-0.276041666666667

0

0

0].

Does anyone have a clue what am I missing here?

Thanks,

-A

For the validation purpose, I am trying to export FEniCS assembly (Poisson's equation) matrix and vector to Matlab and solve using (u = A\b). However, when I solve externally the solution vector looks different (same values but different sparsity) compared to the solution vector extracted from internal FEniCS solve(A, U, b).

For example, take a look at the following Snippet:

A = assemble(a)

b = assemble(L)

bc.apply(A, b)

sio.savemat('Ab.mat', {'A':A.array(), 'b':b.array()})

u = Function(V)

U = u.vector()

solve(A, U, b)

U = [0

0

0

0

0

0

0

0

-0.479166666666669

-0.276041666666667

-0.244791666666667

-0.244791666666667

-0.276041666666667]

Matlab solve (u = A\b)

u=[0

0

-0.244791666666667

0

-0.479166666666669

-0.276041666666667

0

-0.244791666666667

0

-0.276041666666667

0

0

0].

Does anyone have a clue what am I missing here?

Thanks,

-A

Community: FEniCS Project

### 1 Answer

0

I'm not sure what MATLAB is doing, but here is a MWE:

```
from fenics import *
mesh = UnitIntervalMesh(3)
Q = FunctionSpace(mesh, 'CG', 1)
u = TrialFunction(Q)
v = TestFunction(Q)
x = Function(Q)
f = Constant(1.0)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(Q, 0.0, boundary)
A = assemble( u.dx(0) * v.dx(0) * dx)
b = assemble( f * v * dx)
bc.apply(A,b)
solve(A, x.vector(), b)
print "A:\n", A.array(), "\n"
print "B:\n", b.array(), "\n"
print "x:\n", x.vector().array(), "\n"
from scipy.linalg import solve as sci_solve
x_sci = sci_solve(A.array(), b.array())
print "x_sci:\n", x_sci, "\n"
print "|| x - x_sci ||_oo:", \
(x.vector() - x_sci).max()
```

It is easy to verify that the tensors and solution are correct for this function space, see for example how to build the tensor from

https://www.amazon.com/J-N-Reddy-Introduction-Element/dp/B008VS0BPC/ref=sr_1_1?s=books&ie=UTF8&qid=1507387554&sr=1-1&keywords=An+Introduction+to+the+Finite+Element+Method+2nd

Hope that helps...

Thank you pf4d! That is helpful.

If I want to visualize the solution after solving externally, let's say using "pdesurf in Matlab". If I use the same mesh and provide "pdesurf(points, elements, U)". I don't get the same plot as I get from FEniCS.

Because the non-zero solution values are associated with different nodal points.

As I posted in the comment:

Printed soution after solving:

U = [ 0. 0. -0.24479167 0. -0.47916667 -0.27604167 0. -0.24479167 0. -0.27604167 0. 0. 0.]

And when I open "poisson000000.vtu" it is:

U = [ 0. 0. 0. 0. 0. 0. 0. 0. -0.47916667 -0.27604167 -0.24479167 -0.24479167 -0.27604167]

-A

If I want to visualize the solution after solving externally, let's say using "pdesurf in Matlab". If I use the same mesh and provide "pdesurf(points, elements, U)". I don't get the same plot as I get from FEniCS.

Because the non-zero solution values are associated with different nodal points.

As I posted in the comment:

Printed soution after solving:

U = [ 0. 0. -0.24479167 0. -0.47916667 -0.27604167 0. -0.24479167 0. -0.27604167 0. 0. 0.]

And when I open "poisson000000.vtu" it is:

U = [ 0. 0. 0. 0. 0. 0. 0. 0. -0.47916667 -0.27604167 -0.24479167 -0.24479167 -0.27604167]

-A

written
10 months ago by
AJ

.vtu files are used for viewing with visit or paraview; save your solution with numpy.savetxt() or scipy.savemat().

written
10 months ago by
pf4d

Or, if you have HDF5 support in Matlab, save in parallel with that format:

https://www.mathworks.com/help/matlab/hdf5-files.html

see my comment here:

https://www.allanswered.com/post/aqox/#pzew

hope that helps!

https://www.mathworks.com/help/matlab/hdf5-files.html

see my comment here:

https://www.allanswered.com/post/aqox/#pzew

hope that helps!

written
10 months ago by
pf4d

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

Both A and b are assembled in FEniCS and after solving for U in FEniCS I export them to Matlab and then solve again. So A and b in both places are one and the same.

I wonder if FEniCS changes something inside during the cal of solve(A, U, b).

-A

That is, if you were to manually create a MATLAB matrix and vector from the result of printing the Python tenors "A.array()" and "b.array()", and then solve, what do you get then?

I posed the question wrongly (my apologies).

My concern is when we print solution vector after solving it is:

U = [ 0. 0. -0.24479167 0. -0.47916667 -0.27604167 0. -0.24479167 0. -0.27604167 0. 0. 0.]

And when I open "poisson000000.vtu" it is:

U = [ 0. 0. 0. 0. 0. 0. 0. 0. -0.47916667 -0.27604167 -0.24479167 -0.24479167 -0.27604167]

How did the non-zero positions change between two?

If I solve the FEniCS assembled system externally how can I do the same so that I can plot the solution vector using mesh loaded externally?

Maybe this will help you:

https://github.com/mikaem/fenicstools/wiki/DofMap-plotter

I am more interested in using assembly routines from FEniCS and then use the in-house DDM solver to solve the system. To do that, I need the Assembly matrix and vector in its original node order.

Thanks for all your comments, I appreciate your help.

-A

-A