EB: First batch of Marching Cubes

Src/EB/AMReX_EB2.H

@@ -159,7 +159,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,7 @@
 #include <AMReX_EB2.H>
 #include <AMReX_EB2_IndexSpace_STL.H>
 #include <AMReX_EB2_IndexSpace_chkpt_file.H>
+#include <AMReX_MarchingCubes.H>
 #include <AMReX_ParmParse.H>
 #include <AMReX.H>
 #include <algorithm>
@@ -31,12 +32,15 @@ void Initialize ()
     pp.queryAdd("extend_domain_face", extend_domain_face);
     pp.queryAdd("num_coarsen_opt", num_coarsen_opt);
 
+    amrex::MC::Initialize();
+
     amrex::ExecOnFinalize(Finalize);
 }
 
 void Finalize ()
 {
     IndexSpace::clear();
+    amrex::MC::Finalize();
 }
 
 bool ExtendDomainFace ()
@@ -83,7 +87,7 @@ 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;
@@ -226,7 +230,7 @@ Build (const Geometry& geom, int required_coarsening_level,
                                            build_coarse_level_by_coarsening,
                                            a_extend_domain_face,
                                            a_num_coarsen_opt,
-                                           stl_use_bvh));
+                                           stl_use_bvh, support_mvmc));
     }
     else
     {
@@ -294,4 +298,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;

Src/EB/AMReX_EB2_Level.H

@@ -11,6 +11,7 @@
 #include <AMReX_Array.H>
 #include <AMReX_EBCellFlag.H>
 #include <AMReX_MultiCutFab.H>
+#include <AMReX_MarchingCubes.H>
 #include <AMReX_EB2_MultiGFab.H>
 #include <AMReX_EB2_C.H>
 #include <AMReX_EB2_IF_AllRegular.H>
@@ -19,10 +20,16 @@
 #include <omp.h>
 #endif
 
-#include <unordered_map>
-#include <limits>
 #include <cmath>
+#include <limits>
+#include <map>
 #include <type_traits>
+#include <unordered_map>
+
+namespace amrex
+{
+    class STLtools;
+}
 
 namespace amrex::EB2 {
 
@@ -83,6 +90,8 @@ protected:
     BoxArray m_grids;
     BoxArray m_covered_grids;
     DistributionMapping m_dmap;
+    MultiFab m_sdf;
+    std::map<int,std::unique_ptr<MC::MCFab>> m_marching_cubes;
     MultiGFab m_mgf;
     MultiFab m_levelset;
     FabArray<EBCellFlagFab> m_cellflag;
@@ -125,13 +134,27 @@ public:
     void define_fine (G const& gshop, const Geometry& geom,
                       int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt);
 
+#if (AMREX_SPACEDIM == 3)
+    template <typename GS = G, std::enable_if_t<std::is_same_v<GS,STLtools>,int> FOO = 0>
+    void define_fine_mvmc (GS const& gshop, const Geometry& geom,
+                           int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt);
+#endif
+
     static GShopLevel<G>
     makeAllRegular(IndexSpace const* is, const Geometry& geom)
     {
         GShopLevel<G> r(is, geom);
         r.setRegularLevel();
         return r;
     }
+
+private:
+
+    Box m_bounding_box;
+
+    void prepare_grids (G const& gshop, Geometry const& geom, int max_grid_size, int ngrow,
+                        bool extend_domain_face, int num_crse_opt);
+
 };
 
 template <typename G>
@@ -155,29 +178,14 @@ GShopLevel<G>::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry&
 
 template <typename G>
 void
-GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
-                            int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
+GShopLevel<G>::prepare_grids (G const& gshop, const Geometry& geom,
+                              int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
 {
     if (amrex::Verbose() > 0 && extend_domain_face == false) {
         amrex::Print() << "AMReX WARNING: extend_domain_face=false is not recommended!\n";
     }
 
-    BL_PROFILE("EB2::GShopLevel()-fine");
-
-#ifdef AMREX_USE_FLOAT
-    Real small_volfrac = 1.e-5_rt;
-#else
-    Real small_volfrac = 1.e-14;
-#endif
-    bool cover_multiple_cuts = false;
-    int maxiter = 32;
-    {
-        ParmParse pp("eb2");
-        pp.queryAdd("small_volfrac", small_volfrac);
-        pp.queryAdd("cover_multiple_cuts", cover_multiple_cuts);
-        pp.queryAdd("maxiter", maxiter);
-    }
-    maxiter = std::min(100000, maxiter);
+    BL_PROFILE("EB2::GShopLevel()-prepare_grids");
 
     // make sure ngrow is multiple of 16
     m_ngrow = IntVect{static_cast<int>(std::ceil(ngrow/16.)) * 16};
@@ -192,8 +200,8 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
         }
     }
     domain_grown.grow(m_ngrow);
-    Box bounding_box = (extend_domain_face) ? domain : domain_grown;
-    bounding_box.surroundingNodes();
+    m_bounding_box = (extend_domain_face) ? domain : domain_grown;
+    m_bounding_box.surroundingNodes();
 
     BoxList cut_boxes;
     BoxList covered_boxes;
@@ -205,7 +213,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
     for (int clev = num_crse_opt; clev >= 0; --clev) {
         IntVect crse_ratio(1 << clev);
         if (domain.coarsenable(crse_ratio)) {
-            Box const& crse_bounding_box = amrex::coarsen(bounding_box, crse_ratio);
+            Box const& crse_bounding_box = amrex::coarsen(m_bounding_box, crse_ratio);
             Geometry const& crse_geom = amrex::coarsen(geom, crse_ratio);
             BoxList test_boxes;
             if (cut_boxes.isEmpty()) {
@@ -277,11 +285,37 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
         m_dmap = DistributionMapping();
         m_allregular = true;
         m_ok = true;
-        return;
+    } else {
+        m_grids = BoxArray(std::move(cut_boxes));
+        m_dmap = DistributionMapping(m_grids);
     }
+}
 
-    m_grids = BoxArray(std::move(cut_boxes));
-    m_dmap = DistributionMapping(m_grids);
+template <typename G>
+void
+GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
+                            int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
+{
+    BL_PROFILE("EB2::GShopLevel()-fine");
+
+#ifdef AMREX_USE_FLOAT
+    Real small_volfrac = 1.e-5_rt;
+#else
+    Real small_volfrac = 1.e-14;
+#endif
+    bool cover_multiple_cuts = false;
+    int maxiter = 32;
+    {
+        ParmParse pp("eb2");
+        pp.queryAdd("small_volfrac", small_volfrac);
+        pp.queryAdd("cover_multiple_cuts", cover_multiple_cuts);
+        pp.queryAdd("maxiter", maxiter);
+    }
+    maxiter = std::min(100000, maxiter);
+
+    prepare_grids(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
+
+    if (m_allregular) { return; }
 
     m_mgf.define(m_grids, m_dmap);
     const int ng = GFab::ng;
@@ -305,6 +339,7 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
     const auto dx = geom.CellSizeArray();
     const auto problo = geom.ProbLoArray();
 
+    auto bounding_box = m_bounding_box;
     for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
         if (!extend_domain_face || geom.isPeriodic(idim)) {
             bounding_box.grow(idim,GFab::ng);
@@ -537,6 +572,41 @@ GShopLevel<G>::define_fine (G const& gshop, const Geometry& geom,
     m_ok = true;
 }
 
+#if (AMREX_SPACEDIM == 3)
+template <typename G>
+template <typename GS, std::enable_if_t<std::is_same_v<GS,STLtools>,int> FOO>
+void
+GShopLevel<G>::define_fine_mvmc (GS const& gshop, const Geometry& geom,
+                                 int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
+{
+    BL_PROFILE("EB2::GShopLevel()-fine-mvmc");
+
+    prepare_grids(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
+
+    if (m_allregular) { return; }
+
+    m_sdf.define(amrex::convert(m_grids,IntVect(1)), m_dmap, 1, IntVect(1));
+    gshop.fillSignedDistance(m_sdf, m_sdf.nGrowVect(), geom);
+
+    for (MFIter mfi(m_sdf,MFItInfo().DisableDeviceSync()); mfi.isValid(); ++mfi) {
+        m_marching_cubes[mfi.index()] = std::make_unique<MC::MCFab>();
+    }
+
+    MFItInfo info{};
+
+#if defined(AMREX_USE_OMP) && !defined(AMREX_USE_GPU)
+    info.SetDynamic(true);
+#pragma omp parallel
+#endif
+    for (MFIter mfi(m_sdf,info); mfi.isValid(); ++mfi) {
+        amrex::MC::marching_cubes(geom, m_sdf[mfi], *m_marching_cubes[mfi.index()]);
+    }
+
+    MC::write_stl("test.stl", m_marching_cubes);
+
+    m_ok = false; // xxxxx TODO
+}
+#endif
 
 template <typename G>
 GShopLevel<G>::GShopLevel (IndexSpace const* is, int /*ilev*/, int max_grid_size, int /*ngrow*/,

Src/EB/AMReX_EB2_Level_STL.H

@@ -13,7 +13,8 @@ class STLLevel
 public:
 
     STLLevel (IndexSpace const* is, STLtools const& stl_tools, const Geometry& geom,
-              int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt);
+              int max_grid_size, int ngrow, bool extend_domain_face,
+              int num_crse_opt, bool support_mvmc);
 
     STLLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow,
               const Geometry& geom, STLLevel& fineLevel);

Src/EB/AMReX_EB2_Level_STL.cpp

@@ -3,12 +3,20 @@
 namespace amrex::EB2 {
 
 STLLevel::STLLevel (IndexSpace const* is, STLtools const& stl_tools, const Geometry& geom,
-                    int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt)
+                    int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt,
+                    bool support_mvmc)
     : GShopLevel<STLtools>(is, geom)
 {
     BL_PROFILE("EB2::STLLevel()-fine");
 
-    define_fine(stl_tools, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
+#if (AMREX_SPACEDIM == 3)
+    if (support_mvmc) {
+        define_fine_mvmc(stl_tools, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
+    } else
+#endif
+    {
+        define_fine(stl_tools, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt);
+    }
 }
 
 STLLevel::STLLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow,