Segmentation fault in GetFaceOrder for L2 spaces

[termi-official]

Trying to use the Kelly error indicator on L2 spaces crashes. While the Kelly error indicator is in not necessarily the recommended method for DG problems it should at least not crash. This can be tracked down to a crash in GetFaceOrder. Discovered in while testing https://github.com/mfem/mfem/pull/3693 , but resolving this is a separate issue.

Minimal reproducer:

#include "mfem.hpp"
#include <iostream>

using namespace std;
using namespace mfem;

int main(int argc, char *argv[]) {
    const char *mesh_file = "../data/star-hilbert.mesh";
    int order=1;

    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
    args.AddOption(&order, "-o", "--order",
                  "Finite element order (polynomial degree).");
    args.Parse();
    if (!args.Good())   
    {
        args.PrintUsage(cout);
        return 1;
    }
    args.PrintOptions(cout);

    Mesh mesh(mesh_file);
    auto dim = mesh.Dimension();
    auto sdim = mesh.SpaceDimension();

    L2_FECollection flux_fec(order, dim);
    FiniteElementSpace fespace(&mesh, &flux_fec);
    std::cout << fespace.GetFaceOrder(0) << std::endl; // Crashes

    return 0;
}

We might want to add back these tests after the issue has been fixed to capture a regression in the error indicator, in addition to test coverage for GetFaceOrder on L2 spaces.

TEST_CASE("Kelly Error Estimator on 2D NCMesh L2 Crash Test",
          "[NCMesh]")
{
   // Setup
   const auto order = GENERATE(1, 3, 5);
   const auto element = GENERATE(Element::TRIANGLE, Element::QUADRILATERAL);
   CAPTURE(order, element);
   Mesh mesh = Mesh::MakeCartesian2D(2, 2, element);

   // Make the mesh NC
   mesh.EnsureNCMesh();
   {
      Array<int> elements_to_refine(1);
      elements_to_refine[0] = 1;
      mesh.GeneralRefinement(elements_to_refine, 1, 0);
   }

   L2_FECollection fe_coll(order, mesh.Dimension());
   FiniteElementSpace fespace(&mesh, &fe_coll);

   L2_FECollection flux_fec(order-1, mesh.Dimension());
   FiniteElementSpace flux_fes(&mesh, &flux_fec, mesh.SpaceDimension());

   FunctionCoefficient u_analytic(testhelper::SmoothSolutionX);
   GridFunction u_gf(&fespace);
   u_gf.ProjectCoefficient(u_analytic);

   DiffusionIntegrator di;
   KellyErrorEstimator estimator(di, u_gf, flux_fes);
   estimator.GetLocalErrors();
   CHECK(estimator.GetTotalError() >= 0.0);
}

TEST_CASE("Kelly Error Estimator on 3D NCMesh L2 Crash Test",
          "[NCMesh]")
{
   // Setup
   const auto order = GENERATE(1, 3, 5);
   const auto element = GENERATE(Element::TETRAHEDRON, Element::HEXAHEDRON);
   CAPTURE(order, element);
   Mesh mesh = Mesh::MakeCartesian3D(2, 2, 2, element);

   // Make the mesh NC
   mesh.EnsureNCMesh();
   {
      Array<int> elements_to_refine(1);
      elements_to_refine[0] = 1;
      mesh.GeneralRefinement(elements_to_refine, 1, 0);
   }

   L2_FECollection fe_coll(order, mesh.Dimension());
   FiniteElementSpace fespace(&mesh, &fe_coll);

   L2_FECollection flux_fec(order-1, mesh.Dimension());
   FiniteElementSpace flux_fes(&mesh, &flux_fec, mesh.SpaceDimension());

   FunctionCoefficient u_analytic(testhelper::SmoothSolutionX);
   GridFunction u_gf(&fespace);
   u_gf.ProjectCoefficient(u_analytic);

   DiffusionIntegrator di;
   KellyErrorEstimator estimator(di, u_gf, flux_fes);
   estimator.GetLocalErrors();
   CHECK(estimator.GetTotalError() >= 0.0);
}

[v-dobrev]

L2 spaces have no dofs associated with faces (i.e. dofs that are shared with face-neighbor elements) and consequently they have no FEs associated with faces and no order associated with faces. Thus, producing an error if a user asks for the face order in a DG space is the correct behavior.

When you say "crashes", do you mean the code segfaults? If that is the case, we need to fix this and call one of the error macros, MFEM_ABORT, MFEM_VERIFY, or similar.

[termi-official]

The provided MWE unfortunately just segfaults when compiled with MFEM_USE_LIBUNWIND=YES and MFEM_DEBUG=YES. See:

> ./mwe_l2_faceorder_crash -o 0
Options used:
   --mesh ../data/star-hilbert.mesh
   --order 0
zsh: segmentation fault (core dumped)  ./mwe_l2_faceorder_crash -o 0