EB: First batch of Marching Cubes (#4778)

This adds the marching cubes algorithm of "Efficient implementation of
Marching Cubes' cases with topological guarantees" by Lewiner, Lopes,
Vieira & Tavares, Journal of Graphics Tools 8(2): pp. 1-15 (2003). The
implementation is adapted from the source code available at
http://thomas.lewiner.org/publication_page.php%EF%B9%96pubkey=marching_cubes_jgt.html.

Given a signed distance function (e.g., from our STL tools), this
generates a list of topologically consistent triangles. Every cut cell's
EB surfaces are decomposed into up to 9 triangles. There could be
multiple cuts and/or volumes in a cell. The cuts are shared faces and
edges between neighbor cells are consistent.

For debugging purpose, we have also added a function that can save these
triangles into a STL file.

Future work involves the computation of EB information such as volume
fraction, area faction, etc. based on the triangles.

Author: WeiqunZhang

Src/Base/AMReX.H

@@ -315,6 +315,8 @@ namespace amrex
         // be moved to the top.
         static void push (AMReX* pamrex);
 
+        static void push (std::unique_ptr<AMReX> pamrex);
+
         // This erases `pamrex` from the stack.
         static void erase (AMReX* pamrex);
 

Src/Base/AMReX.cpp

@@ -772,7 +772,7 @@ amrex::Initialize (int& argc, char**& argv, bool build_parm_parse,
 
     BL_TINY_PROFILE_INITIALIZE();
 
-    AMReX::push(new AMReX()); // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
+    AMReX::push(std::make_unique<AMReX>());
     return AMReX::top(); // NOLINT
 }
 
@@ -981,6 +981,12 @@ AMReX::push (AMReX* pamrex)
     }
 }
 
+void
+AMReX::push (std::unique_ptr<AMReX> pamrex)
+{
+    m_instance.push_back(std::move(pamrex));
+}
+
 void
 AMReX::erase (AMReX* pamrex)
 {

Src/EB/AMReX_EB2.H

@@ -41,6 +41,8 @@ public:
     // moved to the top.
     static void push (IndexSpace* ispace);
 
+    static void push (std::unique_ptr<IndexSpace> ispace);
+
     // This erases `ispace` from the stack.
     static void erase (IndexSpace* ispace);
 
@@ -133,13 +135,11 @@ Build (const G& gshop, const Geometry& geom,
        int num_coarsen_opt = NumCoarsenOpt())
 {
     BL_PROFILE("EB2::Initialize()");
-    IndexSpace::push(new IndexSpaceImp<G>(gshop, geom,
-                                          required_coarsening_level,
-                                          max_coarsening_level,
-                                          ngrow, build_coarse_level_by_coarsening,
-                                          extend_domain_face,
-                                          num_coarsen_opt));
-} // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
+    IndexSpace::push(std::make_unique<IndexSpaceImp<G>>
+                     (gshop, geom, required_coarsening_level, max_coarsening_level,
+                      ngrow, build_coarse_level_by_coarsening,extend_domain_face,
+                      num_coarsen_opt));
+}
 
 /// \ingroup amrex_eb
 template <typename G>
@@ -151,11 +151,9 @@ Build (const G& gshop, Vector<Geometry> geom,
 {
     BL_PROFILE("EB2::Initialize()");
     std::sort(geom.begin(), geom.end(), [] (Geometry const& a, Geometry const& b) { return a.Domain().numPts() > b.Domain().numPts(); });
-    IndexSpace::push(new IndexSpaceImp<G>(gshop, geom,
-                                          ngrow,
-                                          extend_domain_face,
-                                          num_coarsen_opt));
-} // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
+    IndexSpace::push(std::make_unique<IndexSpaceImp<G>>
+                     (gshop, geom, ngrow, extend_domain_face, num_coarsen_opt));
+}
 
 /// \ingroup amrex_eb
 void Build (const Geometry& geom,
@@ -164,7 +162,16 @@ void Build (const Geometry& geom,
             int ngrow = 4,
             bool build_coarse_level_by_coarsening = true,
             bool extend_domain_face = ExtendDomainFace(),
-            int num_coarsen_opt = NumCoarsenOpt());
+            int num_coarsen_opt = NumCoarsenOpt(),
+            bool support_mvmc = false);
+
+void BuildMultiValuedMultiCut (const Geometry& geom,
+                               int required_coarsening_level,
+                               int max_coarsening_level,
+                               int ngrow = 4,
+                               bool build_coarse_level_by_coarsening = true,
+                               bool extend_domain_face = ExtendDomainFace(),
+                               int num_coarsen_opt = NumCoarsenOpt());
 
 void BuildFromChkptFile (std::string const& fname,
                          const Geometry& geom,

Src/EB/AMReX_EB2.cpp

@@ -12,6 +12,9 @@
 #include <AMReX_EB2.H>
 #include <AMReX_EB2_IndexSpace_STL.H>
 #include <AMReX_EB2_IndexSpace_chkpt_file.H>
+#if (AMREX_SPACEDIM == 3)
+#  include <AMReX_MarchingCubes.H>
+#endif
 #include <AMReX_ParmParse.H>
 #include <AMReX.H>
 #include <algorithm>
@@ -42,6 +45,9 @@ void Initialize ()
 void Finalize ()
 {
     IndexSpace::clear();
+#if (AMREX_SPACEDIM == 3)
+    amrex::MC::Finalize();
+#endif
 }
 
 bool ExtendDomainFace ()
@@ -67,6 +73,12 @@ IndexSpace::push (IndexSpace* ispace)
     }
 }
 
+void
+IndexSpace::push (std::unique_ptr<IndexSpace> ispace)
+{
+    m_instance.push_back(std::move(ispace));
+}
+
 void
 IndexSpace::erase (IndexSpace* ispace)
 {
@@ -88,12 +100,16 @@ const IndexSpace* TopIndexSpaceIfPresent() noexcept {
 void
 Build (const Geometry& geom, int required_coarsening_level,
        int max_coarsening_level, int ngrow, bool build_coarse_level_by_coarsening,
-       bool a_extend_domain_face, int a_num_coarsen_opt)
+       bool a_extend_domain_face, int a_num_coarsen_opt, bool support_mvmc)
 {
     ParmParse pp("eb2");
     std::string geom_type;
     pp.get("geom_type", geom_type);
 
+    if (amrex::Verbose() && support_mvmc && geom_type != "stl") {
+        amrex::Warning("EB2:Build: support_mvmc = true is ignored if eb2.geom_type is not stl.");
+    }
+
     if (geom_type == "all_regular")
     {
         EB2::AllRegularIF rif;
@@ -223,15 +239,12 @@ Build (const Geometry& geom, int required_coarsening_level,
         pp.queryAdd("stl_reverse_normal", stl_reverse_normal);
         bool stl_use_bvh = true;
         pp.queryAdd("stl_use_bvh", stl_use_bvh);
-        IndexSpace::push(new IndexSpaceSTL(stl_file, stl_scale, // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
-                                           {stl_center[0], stl_center[1], stl_center[2]},
-                                           int(stl_reverse_normal),
-                                           geom, required_coarsening_level,
-                                           max_coarsening_level, ngrow,
-                                           build_coarse_level_by_coarsening,
-                                           a_extend_domain_face,
-                                           a_num_coarsen_opt,
-                                           stl_use_bvh));
+        IndexSpace::push(std::make_unique<IndexSpaceSTL>
+                         (stl_file, stl_scale,
+                          Array<Real,3>{stl_center[0], stl_center[1], stl_center[2]},
+                          int(stl_reverse_normal), geom, required_coarsening_level,
+                          max_coarsening_level, ngrow, build_coarse_level_by_coarsening,
+                          a_extend_domain_face, a_num_coarsen_opt, stl_use_bvh, support_mvmc));
     }
     else
     {
@@ -263,12 +276,10 @@ BuildFromChkptFile (std::string const& fname,
                     bool a_extend_domain_face)
 {
     ChkptFile chkpt_file(fname);
-    IndexSpace::push(new IndexSpaceChkptFile(chkpt_file, // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
-                     geom, required_coarsening_level,
-                     max_coarsening_level, ngrow,
-                     build_coarse_level_by_coarsening,
-                     a_extend_domain_face));
-} // NOLINT(clang-analyzer-cplusplus.NewDeleteLeaks)
+    IndexSpace::push(std::make_unique<IndexSpaceChkptFile>
+                     (chkpt_file, geom, required_coarsening_level, max_coarsening_level,
+                      ngrow, build_coarse_level_by_coarsening, a_extend_domain_face));
+}
 
 namespace {
 int comp_max_crse_level (Box cdomain, const Box& domain)
@@ -299,4 +310,15 @@ maxCoarseningLevel (IndexSpace const* ebis, const Geometry& geom)
     return comp_max_crse_level(cdomain,domain);
 }
 
+void
+BuildMultiValuedMultiCut (const Geometry& geom, int required_coarsening_level,
+                          int max_coarsening_level, int ngrow,
+                          bool build_coarse_level_by_coarsening,
+                          bool a_extend_domain_face, int a_num_coarsen_opt)
+{
+    EB2::Build(geom, required_coarsening_level, max_coarsening_level, ngrow,
+               build_coarse_level_by_coarsening, a_extend_domain_face,
+               a_num_coarsen_opt, true);
+}
+
 }

Src/EB/AMReX_EB2_IndexSpace_STL.H

@@ -19,7 +19,8 @@ public:
                   const Geometry& geom, int required_coarsening_level,
                   int max_coarsening_level, int ngrow,
                   bool build_coarse_level_by_coarsening,
-                  bool extend_domain_face, int num_coarsen_opt, bool bvh_optimization);
+                  bool extend_domain_face, int num_coarsen_opt,
+                  bool bvh_optimization, bool support_mvmc);
 
     IndexSpaceSTL (IndexSpaceSTL const&) = delete;
     IndexSpaceSTL (IndexSpaceSTL &&) = delete;

Src/EB/AMReX_EB2_IndexSpace_STL.cpp

@@ -8,7 +8,7 @@ IndexSpaceSTL::IndexSpaceSTL (const std::string& stl_file, Real stl_scale,
                               int max_coarsening_level, int ngrow,
                               bool build_coarse_level_by_coarsening,
                               bool extend_domain_face, int num_coarsen_opt,
-                              bool bvh_optimization)
+                              bool bvh_optimization, bool support_mvmc)
 {
     Gpu::LaunchSafeGuard lsg(true); // Always use GPU
 
@@ -31,7 +31,10 @@ 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);
+                            extend_domain_face, num_coarsen_opt, support_mvmc);
+
+    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.");
 
     for (int ilev = 1; ilev <= max_coarsening_level; ++ilev)
     {
@@ -56,7 +59,7 @@ 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);
+                                            extend_domain_face, num_coarsen_opt-ilev, support_mvmc);
                 }
             } else {
                 break;