Unable to create ParMesh with higher-order Bergot pyramids

I created an mfem::Mesh consisting of a single second-order Bergot-type pyramid. When I tried to use this to create an mfem::ParMesh object I got the following error:

Assertion failed: (dofs.Size() <= elemvect.Size()) is false:
 --> Size mismatch: length of dofs is 45, length of elemvect is 42
 ... in function: void mfem::Vector::SetSubVector(const mfem::Array<int>&, const mfem::Vector&)
 ... in file: /home/cmacmack/Documents/moose-dev/moose/framework/contrib/mfem/linalg/vector.cpp:742

I will work up a minimal-reproducer and post it here.

I've done some digging and worked out the reason this happens. When creating a ParMesh from a Mesh object, there is a call to FiniteElementCollection::New(glob_fes->FEColl()->Name()); (when creating a copy of the mesh nodes). However, the mesh name does not include information on the type of pyramid being used. Therefore, FiniteElementCollection::New invokes the H1_FECollection constructor that uses the default (Fuentes) pyramid type. Fuentes pyramids have a different number of nodes than Bergot pyramids, so we end up with the length-mismatch issue when copy node information.

The solution to this is to incorporate information on the type of pyramid into FE collection name. I suspect this will take a bit of thought, especially if we want to preserve backwards compatibility. Below is a quick-and-dirty fix I made to check I understood the issue properly.

diff --git a/fem/fe_coll.cpp b/fem/fe_coll.cpp
index 7950795ecc..e3c1119d24 100644
--- a/fem/fe_coll.cpp
+++ b/fem/fe_coll.cpp
@@ -245,8 +245,10 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    }
    else if (!strncmp(name, "H1@", 3))
    {
+     const char * final_underscore = strrchr(name, '_');
+     int pyrtype = (final_underscore - name > 7) ? atoi(final_underscore + 1) : 1;
       fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
-                                BasisType::GetType(name[3]));
+                                BasisType::GetType(name[3]), pyrtype);
    }
    else if (!strncmp(name, "L2_T", 4))
       fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
@@ -1724,17 +1726,17 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
    {
       case BasisType::GaussLobatto:
       {
-         snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
+         snprintf(h1_name, 32, "H1_%dD_P%d_%d", dim, p, pyrtype);
          break;
       }
       case BasisType::Positive:
       {
-         snprintf(h1_name, 32, "H1Pos_%dD_P%d", dim, p);
+         snprintf(h1_name, 32, "H1Pos_%dD_P%d_%d", dim, p, pyrtype);
          break;
       }
       case BasisType::Serendipity:
       {
-         snprintf(h1_name, 32, "H1Ser_%dD_P%d", dim, p);
+         snprintf(h1_name, 32, "H1Ser_%dD_P%d_%d", dim, p, pyrtype);
          break;
       }
       default:
@@ -1743,8 +1745,8 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
                      Quadrature1D::Invalid,
                      "unsupported BasisType: " << BasisType::Name(btype));
 
-         snprintf(h1_name, 32, "H1@%c_%dD_P%d",
-                  (int)BasisType::GetChar(btype), dim, p);
+         snprintf(h1_name, 32, "H1@%c_%dD_P%d_%d",
+                  (int)BasisType::GetChar(btype), dim, p, pyrtype);
       }
    }

Here's a minimal reproducer:

#include <mfem.hpp>
#include <mpi.h>

int main(int argc, char *argv[]) {
  mfem::Vector nodes{
      -1.0, -1.0, 0.0,
       1.0, -1.0, 0.0,
       1.0,  1.0, 0.0,
      -1.0,  1.0, 0.0,
       0.0,  0.0, 1.0,
       0.0, -1.0, 0.0,
       1.0,  0.0, 0.0,
       0.0,  1.0, 0.0,
      -1.0,  0.0, 0.0,
      -0.5, -0.5, 0.5,
       0.5, -0.5, 0.5,
       0.5,  0.5, 0.5,
      -0.5,  0.5, 0.5,
       0.0,  0.0, 0.0
    };
  int elem_indices[5] = {0, 1, 2, 3, 4};
  int elem_attrs[1] = {1};
  mfem::Mesh mesh(nodes.GetData(), 5, elem_indices, mfem::Geometry::Type::PYRAMID, elem_attrs, 1,
                  nullptr, mfem::Geometry::Type::SQUARE, nullptr, 0, 3, 3);
  mfem::H1_FECollection fec(2, 3, mfem::BasisType::ClosedUniform, 0);
  mfem::FiniteElementSpace fes(&mesh, &fec, 3);
  mesh.SetNodalFESpace(&fes);

  mesh.SetNodes(nodes);

  mfem::Mpi::Init(argc, argv);
  mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh, nullptr, 1);
}

Related to this, it looks like the H1_FECollection::Clone() method also doesn't handle Bergot pyramids correctly. That should be an easy fix, at least.

Navigation Menu

Repository navigation

Unable to create ParMesh with higher-order Bergot pyramids #5256

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Unable to create ParMesh with higher-order Bergot pyramids #5256

Description

Activity

cmacmackin commented on Mar 4, 2026

cmacmackin commented on Mar 5, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions