How to sort an array according to the global degree of freedom indexing.
8 months ago by
I have setup a problem that involves Nedelec Edge elements. Have defined the bilinear form, the linear form, applied the boundary conditions, and assembled the matrix. I can now solve the system. Each row of the assembled matrix corresponds to a degree of freedom. I want (x,y,z) arrays of the coordinates of all such dofs in my system, and I need them sorted (re-ordered) in the same indexing that appears in the assembled matrix.
However, I am not sure (and want to know) if, when I use the statements
n = VN.dim() # VN is my function space (Edge elements: Nedelec 1st kind H(curl) d = mesh.geometry().dim() dof_coordinates = VN.tabulate_dof_coordinates().reshape(d,n) dof_x = dof_coordinates[:, 0] dof_y = dof_coordinates[:, 1] dof_z = dof_coordinates[:, 2]
will the dof_coordinates already be sorted according to the indexing that is applied to the matrix (by FEniCS), or do they have to be "re-ordered" in some way so that their order matches the dof indexing in the matrix's rows! In the beginning of my code I also use the statement
parameters["reorder_dofs_serial"] = False
I use the following statement to assemble the matrix. In the following, c is a bilinear form, rhs1 is a linear form corresponding to the right hand side, and BCN are suitable (Dirichlet type) boundary conditions.
C = PETScMatrix() bc = PETScVector() assemble_system(c,rhs1,BCN, x0=None, form_compiler_parameters=None, add_values=False, finalize_tensor=True, keep_diagonal=False, A_tensor=C, b_tensor=bc, backend=PETSc)
Thanks for your time.
Hisham bin Zubair Syed
Community: FEniCS Project
8 months ago by
What you described works fine for Lagrange-spaces. However, as far as i know, you cannot extract the x,y,z-dofs on a Nedelec space in a straightforward way.
Your question might be related to:
https://fenicsproject.org/qa/11701/h-curl-extract-the-component-of-u-u-is-in-h-curlhave you considered
Are you solving your system with 1.: one of the FEniCS solvers or 2.: externally?
If 1: You can interpolate your solution on a Lagrange-Space and use exactly the code-snippet you provided to access the dof coordinates in x,y,z direction.
If 2: I think the easiest thing would be to export your Martix - solve externally - import your solution vector into FEniCS again and then do the same as for "If 1:".
Note:Of course each dof in your Nedelec space has x,y,z coordinates you can access via "tabulate_dof_coordinates()", but your solution vector contains only one magnuitude value (of e.g. the E-field). for each dof! For p1 polynomials, each Nedlec dof is located at the center of an edge and directed along the local edge direction of the corresponding element.
I hope that might help,
8 months ago by
Thanks for your time and reply.
Actually, I want to "construct" a vector of 1s and 0s. The length of the vector should be the same as the number of rows in my system matrix, and the vector should contain a 1 (or a 0) corresponding to the z-coordinate of the degree of freedom corresponding to the matrix row. 1 if negative or zero, and 0 if positive. Whatever I do, I cannot seem to be sure that I have the order right. From FEniCS Q&A (now the archive), I gathered that I need to lay out the dofs associated with my N1curl space sorted with respect to "global indexing". But I don't know how to access this so-called "global indexing".
I will give you a related example. In the Documented Demo of curl-curl (eddy current example) by Garth Wells, a particular matrix G can be obtained by using the factory method (that he wrote) called DiscreteOperators.build_gradient(). Before I stumbled on this example, I tried numerous times to "construct" this matrix "G" which is (should be) sort of simple to construct as it just contains two nnz per row, a 1 and a -1 on particular column locations. But I got the order incorrect everytime. However, DiscreteOperators.build_gradient() sets the matrix correctly, and I now use it. It leads me to wonder if this method DiscreteOperators.build_gradient() can be possible coded purely in python dolfin (its a cpp function that is interfaced into python). Because if so, then I can use the same technique to address my issue too.
Please login to add an answer/comment or follow this question.