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

2 Answers

8 months ago by
Hi Hisham,

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: 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
Hi Raphael,

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.

Hi again,

I don't know if I understood right, but if you want to fill a vector with 1s or 0s according to the z-coordinates of the dofs, you can just use tabulate_dof_coordinates() in the nedelec space and check if the z coord of each dof is greater or smaller than 0 and fill another FEniCS vector on the same nedelec space accoringly. That should perform very fast and works either in serial or in parallel.

Another noteregarding local and global indexing:
- In FEniCS there is a function to evalualte the "ownership_range" for each process (sry, I don't remember the exact command right now). With determining the ownership_range, you can easily map local to global indices by adding the corresponding range for each process.
written 7 months ago by RR  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »