How do I extract coordinates of DOF points in high-order FE Space?

[anmelander]

Hello,

This might be a very simple question, but I have been unable to find an answer.

I have a very simple mesh (a 2D square meshed with triangles). I define a finite element space on the mesh using high-order Lagrange basis functions

H1_FECollection fec(order, dim);
FiniteElementSpace fespace(&mesh, &fec);

From here, how do I get the coordinates of the interpolation (DOF) points? GetNode/GetNodes only seem to return the vertices, which I’m assuming is because the mesh itself is affine.

Best,
Anders

[j-signorelli]

Hi,

This is explained in #63, but it can be a bit confusing because by default Mesh::GetNodes() (for low-order straight-edged meshes) just returns the vertices. For higher-order FE spaces, you can use Mesh::GetNodes(GridFunction &gf), where the grid function is defined on the FE space of higher-order.

Alternatively you can set the FiniteElementSpace of the nodes internal to Mesh using Mesh::SetNodalFESpace to one of higher-order and then call Mesh::GetNodes().

[anmelander]

Thanks for the quick answer! A follow-up question:
Using Mesh::SetNodalFESpace and then trying to extract the coordinates by

GridFunction &nodes = *mesh.GetNodes();
cout << nodes << endl;

I seem to only get the x-coordinates.

If I set vdim in FiniteElementSpace fespace(&mesh, &fec,vdim); to 2, it now also returns the y-coordinates. However, as far as I understand, this means I assume I have two DOF per node in my solution, which is not the case (I'm solving a very simple Poisson problem). Am I misunderstanding something here? Do I have to set up two finite element spaces, one for the solution, and one to get the grid coordinates?

Best,
Anders

[j-signorelli]

You are correct that you would then need to define another FiniteElementSpace to get the coordinates - the GridFunction for the node coordinates "live" on an FE space with vdim = physical dim. Just make sure you don't construct the bilinear forms / linear forms / nonlinear forms with that FE space, otherwise you will have 2 solution DOFs / node.

[anmelander]

Thanks!

[v-dobrev]

Another approach that computes the coordinates on the fly, element by element, is:

void print_fes_coords_by_element(FiniteElementSpace &fes)
{
   DenseMatrix coords;
   for (int i = 0; i < fes.GetNE(); i++)
   {
      const FiniteElement *fe = fes.GetFE(i);
      ElementTransformation *T = fes.GetElementTransformation(i);
      T->Transform(fe->GetNodes(), coords);

      cout << "element " << i << ":\n";
      for (int j = 0; j < coords.Width(); j++)
      {
         cout << "   dof " << j << ':';
         for (int d = 0; d < coords.Height(); d++)
         {
            cout << ' ' << coords(d,j);
         }
         cout << '\n';
      }
   }
}