Merge pull request #260 from awslabs/sjg/mesh-simplex-dev

Mesh preprocessing operations

CHANGELOG.md

@@ -22,6 +22,9 @@ The format of this changelog is based on
     `"WavePortPEC"`).
   - Fixed a bug in divergence-free projection for problems without essential or mixed
     boundary conditions.
+  - Added `"MakeSimplex"` and `"MakeHexahedral"` mesh options to convert an input mesh to
+    all tetrahedra or all hexahedra. Also adds `"SerialUniformLevels"` option to
+    `config["Model"]["Refinement"]` for testing or debugging.
 
 ## [0.13.0] - 2024-05-20
 

docs/src/config/model.md

@@ -22,11 +22,14 @@ with
 
 `"Mesh" [None]` :  Input mesh file path, an absolute path is recommended.
 
-`"L0" [1.0e-6]` :  Mesh vertex coordinate length unit, m.
+`"L0" [1.0e-6]` :  Unit, relative to m, for mesh vertex coordinates. For example, a value
+of `1.0e-6` implies the mesh coordinates are in μm.
 
 `"Lc" [0.0]` :  Characteristic length scale used for nondimensionalization, specified in
-mesh length units. A value less than or equal to zero uses an internally calculated length
-scale based on the bounding box of the computational domain.
+mesh length units. This keyword should typically not be specified by the user. A value less
+than or equal to zero uses an internally calculated length scale based on the bounding box
+of the computational domain. A value of 1.0 will disable nondimensionalization of lengths
+in the model and all computations will take place in the same units as the mesh.
 
 `"Refinement"` : Top-level object for configuring mesh refinement.
 
@@ -110,10 +113,16 @@ mesh length units.
 
 ### Advanced model options
 
-  - `"Partition" [""]`
-  - `"ReorientTetMesh" [false]`
   - `"RemoveCurvature" [false]`
+  - `"MakeSimplex" [false]`
+  - `"MakeHexahedral" [false]`
+  - `"ReorderElements" [false]`
+  - `"CleanUnusedElements" [true]`
+  - `"AddInterfaceBoundaryElements" [true]`
+  - `"ReorientTetMesh" [false]`
+  - `"Partitioning" [""]`
   - `"MaxNCLevels" [1]`
   - `"MaximumImbalance" [1.1]`
   - `"SaveAdaptIterations" [true]`
   - `"SaveAdaptMesh" [false]`
+  - `"SerialUniformLevels" [0]`

examples/rings/mesh/mesh.jl

@@ -54,7 +54,7 @@ function generate_ring_mesh(;
     end
     gmsh.model.add("rings")
 
-    # Geometry parameters (in um)
+    # Geometry parameters (in μm)
     farfield_radius = 10.0 * outer_radius
 
     # Mesh parameters

examples/rings/rings.json

@@ -8,7 +8,7 @@
   "Model":
   {
     "Mesh": "mesh/rings.msh",
-    "L0": 1.0e-6  // um
+    "L0": 1.0e-6  // μm
   },
   "Domains":
   {

palace/fem/interpolator.cpp

@@ -11,6 +11,14 @@
 namespace palace
 {
 
+namespace
+{
+
+constexpr auto GSLIB_BB_TOL = 0.01;  // MFEM defaults, slightly reduced bounding box
+constexpr auto GSLIB_NEWTON_TOL = 1.0e-12;
+
+}  // namespace
+
 #if defined(MFEM_USE_GSLIB)
 InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh &mesh)
   : op(mesh.GetComm())
@@ -24,27 +32,24 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh
   {
     return;
   }
-  const double bb_t = 0.1;  // MFEM defaults
-  const double newton_tol = 1.0e-12;
-  const int npts = static_cast<int>(iodata.domains.postpro.probe.size());
+  const int dim = mesh.SpaceDimension();
   MFEM_VERIFY(
-      mesh.Dimension() == mesh.SpaceDimension(),
+      mesh.Dimension() == dim,
       "Probe postprocessing functionality requires mesh dimension == space dimension!");
-  mfem::Vector xyz(npts * mesh.SpaceDimension());
+  const int npts = static_cast<int>(iodata.domains.postpro.probe.size());
+  mfem::Vector xyz(npts * dim);
   op_idx.resize(npts);
   int i = 0;
   for (const auto &[idx, data] : iodata.domains.postpro.probe)
   {
-    // Use default ordering byNODES.
-    xyz(i) = data.center[0];
-    xyz(npts + i) = data.center[1];
-    if (mesh.SpaceDimension() == 3)
+    for (int d = 0; d < dim; d++)
     {
-      xyz(2 * npts + i) = data.center[2];
+      // Use default ordering byNODES.
+      xyz(d * npts + i) = data.center[d];
     }
     op_idx[i++] = idx;
   }
-  op.Setup(mesh, bb_t, newton_tol, npts);
+  op.Setup(mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
   op.FindPoints(xyz, mfem::Ordering::byNODES);
   op.SetDefaultInterpolationValue(0.0);
   i = 0;
@@ -69,18 +74,26 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh
 std::vector<double> InterpolationOperator::ProbeField(const mfem::ParGridFunction &U)
 {
 #if defined(MFEM_USE_GSLIB)
-  // Interpolated vector values are returned from GSLIB interpolator byNODES, which we
-  // transform to byVDIM for output.
+  // Interpolated vector values are returned from GSLIB interpolator with the same ordering
+  // as the source grid function, which we transform to byVDIM for output.
   const int npts = op.GetCode().Size();
-  const int dim = U.VectorDim();
-  std::vector<double> vals(npts * dim);
-  mfem::Vector v(npts * dim);
-  op.Interpolate(U, v);
-  for (int d = 0; d < dim; d++)
+  const int vdim = U.VectorDim();
+  std::vector<double> vals(npts * vdim);
+  if (U.FESpace()->GetOrdering() == mfem::Ordering::byVDIM)
   {
-    for (int i = 0; i < npts; i++)
+    mfem::Vector v(vals.data(), npts * vdim);
+    op.Interpolate(U, v);
+  }
+  else
+  {
+    mfem::Vector v(npts * vdim);
+    op.Interpolate(U, v);
+    for (int d = 0; d < vdim; d++)
     {
-      vals[i * dim + d] = v(d * npts + i);
+      for (int i = 0; i < npts; i++)
+      {
+        vals[i * vdim + d] = v(d * npts + i);
+      }
     }
   }
   return vals;
@@ -107,4 +120,184 @@ std::vector<std::complex<double>> InterpolationOperator::ProbeField(const GridFu
   }
 }
 
+namespace fem
+{
+
+void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
+{
+#if defined(MFEM_USE_GSLIB)
+  // Generate list of points where the grid function will be evaluated. If the grid function
+  // to interpolate is an H1 space of the same order as the mesh nodes, we can use the
+  // mesh node points directly. Otherwise, for a different basis order or type, we generate
+  // the interpolation points in the physical space manually.
+  auto &dest_mesh = *V.FESpace()->GetMesh();
+  MFEM_VERIFY(dest_mesh.GetNodes(), "Destination mesh has no nodal FE space!");
+  const int dim = dest_mesh.SpaceDimension();
+  mfem::Vector xyz;
+  mfem::Ordering::Type ordering;
+  const auto *dest_fec_h1 =
+      dynamic_cast<const mfem::H1_FECollection *>(V.FESpace()->FEColl());
+  const auto *dest_nodes_h1 = dynamic_cast<const mfem::H1_FECollection *>(
+      dest_mesh.GetNodes()->FESpace()->FEColl());
+  int dest_fespace_order = V.FESpace()->GetMaxElementOrder();
+  int dest_nodes_order = dest_mesh.GetNodes()->FESpace()->GetMaxElementOrder();
+  dest_mesh.GetNodes()->HostRead();
+  if (dest_fec_h1 && dest_nodes_h1 && dest_fespace_order == dest_nodes_order)
+  {
+    xyz.MakeRef(*dest_mesh.GetNodes(), 0, dest_mesh.GetNodes()->Size());
+    ordering = dest_mesh.GetNodes()->FESpace()->GetOrdering();
+  }
+  else
+  {
+    int npts = 0, offset = 0;
+    for (int i = 0; i < dest_mesh.GetNE(); i++)
+    {
+      npts += V.FESpace()->GetFE(i)->GetNodes().GetNPoints();
+    }
+    xyz.SetSize(npts * dim);
+    mfem::DenseMatrix pointmat;
+    for (int i = 0; i < dest_mesh.GetNE(); i++)
+    {
+      const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
+      mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
+      T.Transform(fe.GetNodes(), pointmat);
+      for (int j = 0; j < pointmat.Width(); j++)
+      {
+        for (int d = 0; d < dim; d++)
+        {
+          // Use default ordering byNODES.
+          xyz(d * npts + offset + j) = pointmat(d, j);
+        }
+      }
+      offset += pointmat.Width();
+    }
+    ordering = mfem::Ordering::byNODES;
+  }
+  const int npts = xyz.Size() / dim;
+
+  // Set up the interpolator.
+  auto &src_mesh = *U.FESpace()->GetMesh();
+  MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!");
+  auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
+  MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
+  mfem::FindPointsGSLIB op(comm);
+  op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
+
+  // Perform the interpolation and fill the target GridFunction (see MFEM's field-interp
+  // miniapp).
+  const int vdim = U.VectorDim();
+  mfem::Vector vals(npts * vdim);
+  op.SetDefaultInterpolationValue(0.0);
+  op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
+  op.Interpolate(xyz, U, vals, ordering);
+  const auto *dest_fec_l2 =
+      dynamic_cast<const mfem::L2_FECollection *>(V.FESpace()->FEColl());
+  if (dest_fec_h1 || dest_fec_l2)
+  {
+    if (dest_fec_h1 && dest_fespace_order != dest_nodes_order)
+    {
+      // H1 with order != mesh order needs to handle duplicated interpolation points.
+      mfem::Vector elem_vals;
+      mfem::Array<int> vdofs;
+      int offset = 0;
+      for (int i = 0; i < dest_mesh.GetNE(); i++)
+      {
+        const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
+        const int elem_npts = fe.GetNodes().GetNPoints();
+        elem_vals.SetSize(elem_npts * vdim);
+        for (int d = 0; d < vdim; d++)
+        {
+          for (int j = 0; j < elem_npts; j++)
+          {
+            // Arrange element values byNODES to align with GetElementVDofs.
+            int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
+                          ? d * npts + offset + j
+                          : (offset + j) * vdim + d;
+            elem_vals(d * elem_npts + j) = vals(idx);
+          }
+        }
+        const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
+        if (dof_trans)
+        {
+          dof_trans->TransformPrimal(elem_vals);
+        }
+        V.SetSubVector(vdofs, elem_vals);
+        offset += elem_npts;
+      }
+    }
+    else
+    {
+      // Otherwise, H1 and L2 copy interpolated values to vdofs.
+      MFEM_ASSERT(V.Size() == vals.Size(),
+                  "Unexpected size mismatch for interpolated values and grid function!");
+      V = vals;
+    }
+  }
+  else
+  {
+    // H(div) or H(curl) use ProjectFromNodes.
+    mfem::Vector elem_vals, v;
+    mfem::Array<int> vdofs;
+    int offset = 0;
+    for (int i = 0; i < dest_mesh.GetNE(); i++)
+    {
+      const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
+      mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
+      const int elem_npts = fe.GetNodes().GetNPoints();
+      elem_vals.SetSize(elem_npts * vdim);
+      for (int d = 0; d < vdim; d++)
+      {
+        for (int j = 0; j < elem_npts; j++)
+        {
+          // Arrange element values byVDIM for ProjectFromNodes.
+          int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
+                        ? d * npts + offset + j
+                        : (offset + j) * vdim + d;
+          elem_vals(j * vdim + d) = vals(idx);
+        }
+      }
+      const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
+      v.SetSize(vdofs.Size());
+      fe.ProjectFromNodes(elem_vals, T, v);
+      if (dof_trans)
+      {
+        dof_trans->TransformPrimal(v);
+      }
+      V.SetSubVector(vdofs, v);
+      offset += elem_npts;
+    }
+  }
+#else
+  MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
+#endif
+}
+
+void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U,
+                         mfem::Vector &vals, mfem::Ordering::Type ordering)
+{
+#if defined(MFEM_USE_GSLIB)
+  // Set up the interpolator.
+  auto &src_mesh = *U.FESpace()->GetMesh();
+  MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!");
+  const int dim = src_mesh.SpaceDimension();
+  const int npts = xyz.Size() / dim;
+  auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
+  MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
+  mfem::FindPointsGSLIB op(comm);
+  op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
+
+  // Perform the interpolation, with the ordering of the returned values matching the
+  // ordering of the source grid function.
+  const int vdim = U.VectorDim();
+  MFEM_VERIFY(vals.Size() == npts * vdim, "Incorrect size for interpolated values vector!");
+  op.SetDefaultInterpolationValue(0.0);
+  op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
+  op.Interpolate(xyz, U, vals, ordering);
+#else
+  MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
+#endif
+}
+
+}  // namespace fem
+
 }  // namespace palace

palace/fem/interpolator.hpp

@@ -35,6 +35,21 @@ class InterpolationOperator
   std::vector<std::complex<double>> ProbeField(const GridFunction &U);
 };
 
+namespace fem
+{
+
+// Interpolate a function on a serial or parallel mesh to a different mesh, using GSLIB.
+// Similar to MFEM's field-interp miniapp.
+void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V);
+
+// Interpolate a function at a specific list of points, specified using the provided
+// ordering. The output vector values are always arranged byVDIM.
+void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U,
+                         mfem::Vector &V,
+                         mfem::Ordering::Type ordering = mfem::Ordering::byNODES);
+
+}  // namespace fem
+
 }  // namespace palace
 
 #endif  // PALACE_FEM_INTERPOLATOR_HPP

palace/main.cpp

@@ -308,7 +308,7 @@ int main(int argc, char *argv[])
   std::vector<std::unique_ptr<Mesh>> mesh;
   {
     std::vector<std::unique_ptr<mfem::ParMesh>> mfem_mesh;
-    mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata, false, true, true, false));
+    mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata));
     iodata.NondimensionalizeInputs(*mfem_mesh[0]);
     mesh::RefineMesh(iodata, mfem_mesh);
     for (auto &m : mfem_mesh)

palace/models/postoperator.cpp

@@ -160,7 +160,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata)
   const bool use_ho = true;
   const int refine_ho = HasE() ? E->ParFESpace()->GetMaxElementOrder()
                                : B->ParFESpace()->GetMaxElementOrder();
-  mesh_Lc0 = iodata.GetLengthScale();
+  mesh_Lc0 = iodata.GetMeshLengthScale();
 
   // Output mesh coordinate units same as input.
   paraview.SetCycle(-1);