Access to edges (2D faces) of neighboring elements on different MPI ranks with an element iterator

[FH-9]

Hi,

For a given element, I need to evaluate minimum entropy as the minimum of all entropies evaluated at solution, volume and surface quadrature nodes of the element itself and its neighbors. I iterate over the elements. Getting the solution nodes (DOFs) (and volume quadrature nodes) in neighbors on a different MPI rank is straightforward with the following two lines:

  vfes->GetFaceNbrElementVDofs(shifted_indx, vdof_indices);
  sol.FaceNbrData().GetSubVector(vdof_indices, el_vdofs);

el_vdofs now has the DOFs that I seek and with an appropriate quadrature rule I can obtain the values at volume quadrature nodes as follows:

DenseMatrix vdof_mat(el_vdofs.GetData(), ndofs, num_equations);
for (int k = 0; k < ir->GetNPoints(); k++)
 {
    const IntegrationPoint &ip = ir->IntPoint(k);
    fe->CalcShape(ip, shape);
    vdof_mat.MultTranspose(shape, state);
   .
   .
   .
}

However, I haven't found an easy way of extracting faces and face quadrature values of neighbors on a different rank with an iterator over elements. GetFaceNbrElementFaces() would probably work in 3D but is there an analogous method for 2D?

Thank you very much for your time.

Farhad

[v-dobrev]

Why not just use IntegrationRule ir that includes points on the boundary? Then you do not need a separate loop over the faces.

[FH-9]

Yes, using Gauss-Lobatto nodes would include some boundary nodes. But then there is no guarantee that the minimum entropy does not occur at a face node which that volume GL rule does not cover. I would like to check at all the nodes that are used in volume and face integrals in DG semi-discretization.

[v-dobrev]

You can create a custom IntegrationRule that combines the exact volume and face points from volume and face quadrature rules that you want.

[FH-9]

Thanks for pointing that out. I am not sure how that would be done in practice though. Would it be possible to somehow take the union of several integration rules? Or do I have to manually add every node?

[v-dobrev]

I was thinking about just manually merging the volume rule and the transformed face rules (for each of the faces of the element), however, you can just add another loop over the faces and their (mapped) quadrature points after your loop over the volume quadrature points. This will look something like this:

   int face_nbr_id = 0; // loop over this index
   const Element *el = pmesh.face_nbr_elements[face_nbr_id];
   const Element::Type etype = el->GetType();
   const Geometry::Type egeom = el->GetGeometryType();
   IntegrationPointTransformation ip_tr;
   for (int f = 0; f < Geometry::NumBdrArray[egeom]; f++)
   {
      Element::Type ftype;
      Geometry::Type fgeom;
      if (Geometry::Dimension[egeom] == 2)
      {
         ftype = Element::SEGMENT;
         fgeom = Geometry::SEGMENT;
      }
      else // dim == 3
      {
         const int face_nv = el->GetNFaceVertices(f);
         ftype = (face_nv == 3) ? Element::TRIANGLE : Element::QUADRILATERAL;
         fgeom = (face_nv == 3) ? Geometry::TRIANGLE : Geometry::SQUARE;
      }
      pmesh.GetLocalFaceTransformation(ftype, etype, ip_tr.Transf, 64*f);
      const int face_ir_order = 1; // set this appropriately
      const IntegrationRule &face_ir = IntRules.Get(fgeom, face_ir_order);
      IntegrationRule mapped_face_ir(face_ir.GetNPoints());
      ip_tr.Transform(face_ir, mapped_face_ir);
      for (int j = 0; j < mapped_face_ir.GetNPoints(); j++)
      {
         const IntegrationPoint &eip = mapped_face_ir.IntPoint(j);

         // 'eip' is a qudrature point in the volume reference space -- use it
         // as in the loop over the volume quadrature points ...
      }
   }

With some modifications, this can be turned into a code that constructs a new composite IntegrationRule that combines the volume and face quadrature points.

[FH-9]

This is very useful. Thank you so much!
On a related note, I also need to include the effect of the boundary state on the minimum entropy of boundary elements. I was thinking of adopting the following snippet from NonlinearForm::Mult():

for (int i = 0; i < fes -> GetNBE(); i++)
{
    const int bdr_attr = mesh->GetBdrAttribute(i);
    if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }

    tr = mesh->GetBdrFaceTransformations (i);
    if (tr != NULL)
    {
        fes->GetElementVDofs(tr->Elem1No, vdofs);
        px.GetSubVector(vdofs, el_x);

        fe1 = fes->GetFE(tr->Elem1No);
        fe2 = fe1;
        for (int k = 0; k < bfnfi.Size(); k++)
        {
            if (bfnfi_marker[k] &&
                (*bfnfi_marker[k])[bdr_attr-1] == 0) { continue; }

         // Compute the proper outer state based on custom boundary integrators 
         // and boundary markers that mark where the boundary integrators should be applied

        }
    }
}

A couple of issues with this:

1. The iterator i is over the boundary elements while my loop is over all elements in the rank. Would it be possible to convert between the two so that I can use mesh->GetBdrAttribute() and mesh->GetBdrFaceTransformations() methods? This is not such a critical issue - it is just I would have to loop over the boundary elements separately again and update the GridFunction for entropy.
2. A more serious issue is that I do not have access to bfnfi_marker since it is a protected member of the NonlinearForm and this calculation is happening outside (in TimeDependentOperator).

I would appreciate any pointers on this.

[v-dobrev]

Regarding (1): if I understand you correctly, you want to find if an element i has any neighbor boundary elements? To do that you can construct the boundary-element to element mapping (using Mesh::GetBdrElementAdjacentElement) as a Table and the transpose it with mfem::Transpose(const Table &A, Table &At, int ncols_A_ = -1) to get element to boundary-element mapping.

Regarding (2): I see that we don't have a method to access bfnfi_marker -- we should probably add one -- do you want to propose a PR for that? In the mean time, you can derive a custom class that adds such a method (derived classes have access to the inherited protected members) and use that instead.