Marching Cubes for EB Generation

Docs/sphinx_documentation/source/EB.rst

@@ -222,6 +222,100 @@ an STL file using:
 This requires setting :cpp:`ParmParse` parameters ``eb2.geom_type = stl`` and
 ``eb2.stl_file`` to specify the STL file path.
 
+.. _sec:EB:ebinit:mc:
+
+**Marching-Cubes Geometry Generation**
+
+By default the cut-cell geometry (volume and area fractions, centroids,
+boundary normals, and connectivity) is generated by the original EB2
+construction. Setting the :cpp:`ParmParse` parameter
+``eb2.geometry_method = marching_cubes`` selects an alternative generator that
+extracts a marching-cubes (MC33) surface from the geometry and builds
+single-valued cut cells from it. The choice applies to every way of
+describing the geometry: the built-in implicit functions selected by
+``eb2.geom_type``, user-defined implicit functions passed to any of the
+:cpp:`EB2::Build` overloads through :cpp:`EB2::makeShop`, and STL files.
+For example,
+
+.. highlight:: bash
+
+::
+
+   eb2.geometry_method = marching_cubes
+   eb2.geom_type = sphere
+   eb2.sphere_center = 0. 0. 0.
+   eb2.sphere_radius = 0.75
+   eb2.sphere_has_fluid_inside = 0
+
+builds the sphere with the marching-cubes generator. The generator has the
+following restrictions:
+
+- It is available only in 3D on Cartesian grids with cubic cells
+  (``dx == dy == dz``); a non-cubic grid is rejected even when the geometry
+  is all regular.
+- It produces single-valued cut cells and cannot be combined with
+  :cpp:`EB2::BuildMultiValuedMultiCut`.
+- STL input must be watertight and consistently oriented. Shared vertices are
+  matched by exact coordinates; invalid files are rejected with a diagnostic
+  that names the offending edge.
+- It applies to the finest level and to any coarse level that is rebuilt
+  from the geometry description (``eb2.build_coarse_level_by_coarsening = 0``
+  or the multi-level :cpp:`EB2::Build` overload); coarse levels obtained by
+  coarsening are unaffected. It is ignored for ``eb2.geom_type = chkpt_file``.
+
+The generator produces at most one fluid volume per cell and one connected
+fluid aperture per Cartesian face. A cell is kept only when its fluid corners
+form a single group connected along cell edges or across ambiguous faces that
+MC33 resolves as connected; configurations whose fluid regions would only be
+joined through the cell interior (the MC33 "tunnel" tilings) or that split
+into several regions are treated as unsupported topology, so sub-cell tunnels
+are never represented. Unsupported topology is an error by
+default. Setting ``eb2.cover_multiple_cuts = 1`` enables the same monotone
+nodal repair used by the default generator: fluid nodes incident to an
+unsupported face or cell are moved to exact zero and the geometry is rebuilt.
+Cells with volume fraction below ``eb2.small_volfrac`` are repaired in the same
+way. Construction stops when no repair candidates remain or fails after
+``eb2.maxiter`` passes. When ``eb2.extend_domain_face`` is true, the geometry
+outside the domain is extruded straight outward from the domain faces.
+
+Implicit functions are sampled at the nodes, without clamping the sample
+positions to the domain (functions that are only defined inside the domain
+must remain finite up to a few cells outside it), and the crossing of every
+Cartesian edge whose end points differ in sign is found by root finding on the
+function; nodes that lie on the surface within roundoff are snapped to exact
+zero, and an edge whose root finding fails falls back to linear interpolation
+of the nodal values. For STL input, exact
+distances are evaluated only in the band needed by cut cells and their
+vertex-normal stencils; away from that band only the inside/outside sign is
+retained. Sign-changing Cartesian edges are intersected with the original STL
+triangles, and the resulting crossings are used for the surface vertices, face
+geometry, edge centroids and cell moments. The public level set and edge
+centroids describe the final repaired domain: fluid nodes are negative,
+covered nodes are positive, and repaired boundary nodes are zero.
+``eb2.mc_stl_file`` optionally writes the converged, repaired triangulation of
+the finest level to an ASCII STL file.
+
+**Planar EB Surface Output**
+
+The single-valued planar boundary reconstructed from an
+:cpp:`EBFArrayBoxFactory` can be written to an ASCII STL file for diagnostic
+visualization:
+
+.. highlight:: c++
+
+::
+
+   #include <AMReX_WriteEBSurface.H>
+
+   amrex::WriteEBSurfaceSTL(ba, dm, geom, factory.get(), "eb_surface.stl");
+
+Each cut cell contributes its independently reconstructed planar polygon,
+which is triangulated in the output. Neighboring planes are not welded along
+shared cell faces, so this diagnostic STL is generally not watertight.
+Applying this writer to factories produced by different EB generators gives a
+representation-controlled comparison of their projected cut-cell data; it is
+not a closed surface-mesh exporter.
+
 **Managing IndexSpace Objects**
 
 Regardless of which :cpp:`Build` variant is used, the newly built

Docs/sphinx_documentation/source/RuntimeParameters.rst

@@ -800,6 +800,28 @@ Embedded Boundary
    Fixing small and multi-cut cells is an iterative process. This parameter
    specifies the maximum number of iterations for the fix-up process.
 
+.. py:data:: eb2.geometry_method
+   :type: string
+   :value: legacy
+
+   Selects the generator used to build cut-cell geometry from the geometry
+   description (the finest EB level, and any coarse level that is rebuilt
+   rather than coarsened) for every :py:data:`eb2.geom_type` except
+   ``chkpt_file`` and for :cpp:`amrex::EB2::Build` calls with a user-provided
+   implicit function.
+   ``legacy`` is the original EB2 construction. ``marching_cubes`` builds
+   single-valued cut cells from a marching-cubes (MC33) surface; it is
+   available in 3D on Cartesian grids with cubic cells only, and it is not
+   compatible with :cpp:`amrex::EB2::BuildMultiValuedMultiCut`. See
+   :ref:`Marching-Cubes Geometry Generation <sec:EB:ebinit:mc>`.
+
+.. py:data:: eb2.mc_stl_file
+   :type: string
+
+   When ``eb2.geometry_method = marching_cubes``, write the final (repaired)
+   marching-cubes triangulation of the finest level to this ASCII STL file.
+   Mostly useful for visualization and debugging.
+
 Error Handling
 --------------
 

Src/EB/AMReX_EB2.H

@@ -149,6 +149,10 @@ int NumCoarsenOpt ();
 /**
  * \brief Build EB geometry from a GeometryShop \p gshop and a single AMR Geometry.
  *
+ * The fine-level generator is selected by the runtime parameter
+ * `eb2.geometry_method` (`legacy`, the default, or `marching_cubes`); see
+ * EB2::GeometryMethod and EB2::GetGeometryMethod().
+ *
  * \tparam G GeometryShop type (see EB2::GeometryShop helpers).
  * \param gshop  Construct describing the implicit surface.
  * \param geom   Finest-level Geometry.
@@ -208,6 +212,9 @@ Build (const G& gshop, Vector<Geometry> geom,
  * ParmParse namespace (notably `eb2.geom_type` along with its shape-specific parameters)
  * and instantiates the matching implicit function.  Supported values currently include
  * `all_regular`, `box`, `cylinder`, `plane`, `sphere`, `torus`, `parser`, and `stl`.
+ * For every geometry type, runtime `eb2.geometry_method` selects the `legacy`
+ * generator (the default) or the 3-D, cubic-cell, single-valued `marching_cubes`
+ * generator (see EB2::GetGeometryMethod).
  *
  * \param geom   Finest-level Geometry.
  * \param required_coarsening_level Number of coarse levels required by the application.
@@ -218,9 +225,9 @@ Build (const G& gshop, Vector<Geometry> geom,
  * \param extend_domain_face Whether to extrapolate the EB description beyond the domain by extending from the domain face.
  * \param num_coarsen_opt Number of coarse probes to try (defaults to EB2::NumCoarsenOpt /
  *        runtime `eb2.num_coarsen_opt`) before rebuilding at full resolution.
- * \param support_mvmc Whether to build multi-valued multi-cut support. (Only honored when
+ * \param support_mvmc Whether to build experimental multi-valued multi-cut support. (Only honored when
  *        `eb2.geom_type` is `stl`; other geometries ignore it and emit a warning.
- *        Even for STL, this is not supported yet.)
+ *        It is incompatible with `eb2.geometry_method=marching_cubes`.)
  */
 void Build (const Geometry& geom,
             int required_coarsening_level,

Src/EB/AMReX_EB2.cpp

@@ -55,6 +55,22 @@ bool ExtendDomainFace ()
     return extend_domain_face;
 }
 
+GeometryMethod GetGeometryMethod ()
+{
+    // Queried at build time rather than in Initialize so that parameters set
+    // programmatically after amrex::Initialize are honored.
+    ParmParse pp("eb2");
+    std::string method_name("legacy");
+    pp.queryAdd("geometry_method", method_name);
+    if (method_name == "legacy") {
+        return GeometryMethod::legacy;
+    } else if (method_name == "marching_cubes") {
+        return GeometryMethod::marching_cubes;
+    }
+    amrex::Abort("eb2.geometry_method must be legacy or marching_cubes, not " + method_name);
+    return GeometryMethod::legacy;
+}
+
 int NumCoarsenOpt ()
 {
     return num_coarsen_opt;

Src/EB/AMReX_EB2_IndexSpaceI.H

@@ -32,7 +32,7 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
     m_ngrow.push_back(ngrow_finest);
     m_gslevel.reserve(max_coarsening_level+1);
     m_gslevel.emplace_back(this, gshop, geom, EB2::max_grid_size, ngrow_finest, extend_domain_face,
-                           num_coarsen_opt);
+                           num_coarsen_opt, true);
 
     for (int ilev = 1; ilev <= max_coarsening_level; ++ilev)
     {
@@ -57,7 +57,7 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Geometry& geom,
                     amrex::Abort("Failed to build required coarse EB level "+std::to_string(ilev));
                 } else {
                     m_gslevel.emplace_back(this, gshop, cgeom, EB2::max_grid_size, ng, extend_domain_face,
-                                           num_coarsen_opt-ilev);
+                                           num_coarsen_opt-ilev, false);
                 }
             } else {
                 break;
@@ -89,7 +89,7 @@ IndexSpaceImp<G>::IndexSpaceImp (const G& gshop, const Vector<Geometry>& geom,
         Geometry cgeom = geom[ilev];
         Box cdomain = cgeom.Domain();
         m_gslevel.emplace_back(this, gshop, cgeom, EB2::max_grid_size, ng, extend_domain_face,
-                                           num_coarsen_opt-ilev);
+                                           num_coarsen_opt-ilev, ilev == 0);
         m_geom.push_back(cgeom);
         m_domain.push_back(cdomain);
         m_ngrow.push_back(ng);

Src/EB/AMReX_EB2_IndexSpace_STL.cpp

@@ -13,6 +13,9 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
     Gpu::LaunchSafeGuard lsg(true); // Always use GPU
 
     STLtools stl_tools;
+    // The marching-cubes generator requires a watertight, consistently oriented
+    // STL; the flag enables that validation when the file is read.
+    stl_tools.setUseMarchingCubes(GetGeometryMethod() == GeometryMethod::marching_cubes);
     stl_tools.setBVHOptimization(bvh_optimization);
     stl_tools.read_stl_file(stl_file, stl_scale, stl_center, stl_reverse_normal);
 
@@ -31,7 +34,7 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
     m_ngrow.push_back(ngrow_finest);
     m_stllevel.reserve(max_coarsening_level+1);
     m_stllevel.emplace_back(this, stl_tools, geom, EB2::max_grid_size, ngrow_finest,
-                            extend_domain_face, num_coarsen_opt, support_mvmc);
+                            extend_domain_face, num_coarsen_opt, support_mvmc, true);
 
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(max_coarsening_level == 0 || support_mvmc == false,
                                      "We don't support multiple levels when multi-valued and multi-cut are enabled.");
@@ -59,7 +62,8 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
                     amrex::Abort("Failed to build required coarse EB level "+std::to_string(ilev));
                 } else {
                     m_stllevel.emplace_back(this, stl_tools, cgeom, EB2::max_grid_size, ng,
-                                            extend_domain_face, num_coarsen_opt-ilev, support_mvmc);
+                                            extend_domain_face, num_coarsen_opt-ilev, support_mvmc,
+                                            false);
                 }
             } else {
                 break;