I am trying to extract the Matrices of “a” and “L” from the assembled system and my “force-vector” is applied to a quarter of the upper boundary of my cube. I saved “L” and also my mesh to a .txt and i thought that the force would be somewhere similar to where the mesh points are, where the force is applied, but that doesn’t seem to be the case, because my force values are all over the place.
So my question is, if someone knows, how i can assemble “a” and “L” but stay true to the mesh or have atleast some consistence.
I had a similar problem with the displacement “u”, but that was pretty easily fixed by this:
u_displacement = np.array([u(Point(x, y, z)) for x, y, z in mesh_points])
but i didn’t find a similar solution for “L”. All i can do at the moment is:
A_num, B_num = assemble_system(a, L, bcs)
but when i save it to .txt:
np.savetxt(“force.txt”, B_num)
and also save my mesh to .txt:
np.savetxt(“meshpoints.txt”, mesh_points)
there seems to be no correlation between the order of these two. I would love some help, cheers.
Hi Jannick, I do not exactly understand the issue but I would suggest trying the following. You can assemble the components (a, L, bcs) more explicitly than assemble_system(a, L, bcs), as shown below.
A = assemble(a)
b = assemble(L)
bcs.apply(A, b)
u = Function(V)
U = u.vector()
solve(A, U, b)
You can dump these into numpy arrays by calling the array function (for example A.array()). Now take note that A is a matrix, while b is a vector and also u.vector() is a vector. You can directly print those arrays onto the screen using the print function or dump those arrays into the MATLAB files using the following commands.
Now you can explore these in MATLAB or Octave directly which I think is more convenient and accurate than in a .txt file. But note that these numpy arrays and MATLAB marices are in dense format.
References and details are in Section 5.3 of the The FEniCS Tutorial 1 - Solving PDEs in Python.
Sadly i can’t really explain the problem any better. The assembling part works for me without problems, but when i look at the vector and Matrix in Matlab it makes absolutely no sense to me in which “order” fenics assembles my system, because the vector should only be applied at the top of my cube but the values of the assembled “L” vector are not where they should be.
An example would be, that my force at point 1, 1, 1 should be the last one in the L array, but thats not the case. Instead it is at position L[87] which to me makes no sense, but maybe i don’t see the greater picture.
When i export u like this:
u_nodal_values = np.array(u.vector())
it has the same “order” as L, but exporting it like this:
u_displacement = np.array([u(Point(x, y, z)) for x, y, z in mesh_points])
it gives me u in a nice order. Maybe that is also possible for L, because i would like to import it into paraview and to do that, i need it in the right order or atleast in reference to an order paraview can understand.
Thank you for your help!
Order of coefficients in vector due to assembly of a linear form over function space V corresponds to the degree of freedom (dof) ordering of that space. For matrix coming from the bilinear form over V\times W, where V, W are respectively trial and test function spaces,
the rows are ordered according to dofs of W while column ordering matches dof order of V. So instead of the mesh points you should be storing the spatial coordinates associated with the degrees of freedom
from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, 'CG', 1)
# Dof coordinates
dofs_x = V.tabulate_dof_coordinates().reshape((V.dim(), -1))
Hey,
thanks for your help first of all, that was actually very insightful. Gladly i found out, that i can change the order of my dof map with the command:
parameters["reorder_dofs_serial"] = False
That still doesn’t give me the exact order paraview wants, but i could change that pretty easily by changing the order of the array in a loop.