### Explicit assembly and external solve

272
views
1
8 months ago by
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

Community: FEniCS Project
I think a good place to start is to verify that A and b are identical between software.
written 8 months ago by pf4d
Hi pf4d, Thanks for the comment.

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
written 8 months ago by AJ
The "export" process is what I refer to; have you verified post export within MATLAB that these tensor elements are the same?

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?
written 8 months ago by pf4d
Thanks for the comment. I get the same solution vector if I copied values from Python tenors "A.array()" and "b.array()" and solve in Matlab.

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?

written 8 months ago by AJ
Oh, I see.  I believe that the internal representation of tensor components has been optimized for speed, and when exporting they change the ordering to something more convenient for post-processing tools...  perhaps.  This is indeed a different question entirely.

https://github.com/mikaem/fenicstools/wiki/DofMap-plotter
written 8 months ago by pf4d
That link can surely help to figure out the reordering. However, if you're aware, is there a way to get the assembly matrix and vector without node reordering, even if it is slow?

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.

-A

written 8 months ago by AJ
Adding "parameters['reorder_dofs_serial'] = False" seems to do the tric.

-A
written 8 months ago by AJ

0
8 months ago by
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
written 8 months ago by AJ
.vtu files are used for viewing with visit or paraview; save your solution with numpy.savetxt() or scipy.savemat().
written 8 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: