diff --git a/Docs/sphinx_documentation/source/EB.rst b/Docs/sphinx_documentation/source/EB.rst index 2020aa34d7a..42503ea44db 100644 --- a/Docs/sphinx_documentation/source/EB.rst +++ b/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::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 diff --git a/Docs/sphinx_documentation/source/RuntimeParameters.rst b/Docs/sphinx_documentation/source/RuntimeParameters.rst index 91043b5e3ff..d3ee2ba95af 100644 --- a/Docs/sphinx_documentation/source/RuntimeParameters.rst +++ b/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 `. + +.. 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 -------------- diff --git a/Src/EB/AMReX_EB2.H b/Src/EB/AMReX_EB2.H index 9e6a58867a5..756e21c0083 100644 --- a/Src/EB/AMReX_EB2.H +++ b/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 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 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, diff --git a/Src/EB/AMReX_EB2.cpp b/Src/EB/AMReX_EB2.cpp index 0d4e99a5c3d..80522d42fbb 100644 --- a/Src/EB/AMReX_EB2.cpp +++ b/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; diff --git a/Src/EB/AMReX_EB2_IndexSpaceI.H b/Src/EB/AMReX_EB2_IndexSpaceI.H index ceac19f73a6..43681083266 100644 --- a/Src/EB/AMReX_EB2_IndexSpaceI.H +++ b/Src/EB/AMReX_EB2_IndexSpaceI.H @@ -32,7 +32,7 @@ IndexSpaceImp::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::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::IndexSpaceImp (const G& gshop, const Vector& 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); diff --git a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp b/Src/EB/AMReX_EB2_IndexSpace_STL.cpp index 5811930cb03..bb5089b49c4 100644 --- a/Src/EB/AMReX_EB2_IndexSpace_STL.cpp +++ b/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; diff --git a/Src/EB/AMReX_EB2_Level.H b/Src/EB/AMReX_EB2_Level.H index 7612adfaf07..30d5977cd9f 100644 --- a/Src/EB/AMReX_EB2_Level.H +++ b/Src/EB/AMReX_EB2_Level.H @@ -13,6 +13,7 @@ #include #if (AMREX_SPACEDIM == 3) # include +# include #endif #include #include @@ -38,6 +39,33 @@ namespace amrex::EB2 { class IndexSpace; template class GShopLevel; +/** + * \brief Fine-level EB geometry generator. + * + * `legacy` is the original EB2 cut-cell construction. `marching_cubes` builds + * single-valued cut cells from a marching-cubes (MC33) surface and is available + * for every geometry source (implicit functions, parser expressions, and STL) + * in 3D on Cartesian grids with cubic cells. + */ +enum class GeometryMethod : int { + legacy = 0, + marching_cubes = 1 +}; + +//! Generator selected at runtime by `eb2.geometry_method` (default `legacy`). +GeometryMethod GetGeometryMethod (); + +//! Runtime parameters shared by the small-cell/multi-cut repair loops. +struct RepairParameters +{ + Real small_volfrac; + int maxiter; + bool cover_multiple_cuts; +}; + +//! Read `eb2.small_volfrac`, `eb2.maxiter`, and `eb2.cover_multiple_cuts`. +RepairParameters GetRepairParameters (); + class Level { public: @@ -247,9 +275,6 @@ protected: BoxArray m_covered_grids; DistributionMapping m_dmap; MultiFab m_sdf; -#if (AMREX_SPACEDIM == 3) - std::map> m_marching_cubes; -#endif MultiGFab m_mgf; MultiFab m_levelset; FabArray m_cellflag; @@ -271,6 +296,24 @@ protected: //! Face-centered data storage: [0]=x-face, [1]=y-face, [2]=z-face Array, AMREX_SPACEDIM> m_fc_data; +#if (AMREX_SPACEDIM == 3) + //! Abort unless \p geom is a Cartesian grid with cubic cells. + static void assert_marching_cubes_supported (Geometry const& geom); + /** + * \brief Allocate the per-FAB marching-cubes caches (with edge-crossing + * storage) that the geometry source fills next. Requires m_sdf. + */ + LayoutData define_marching_cubes_caches (); + /** + * \brief Build every EB field of this level from m_sdf and \p mc_fabs. + * + * Runs the marching-cubes extraction and repair loop, classifies cells, + * optionally writes `eb2.mc_stl_file`, and exports m_levelset. + */ + void build_marching_cubes_level (Geometry const& geom, bool extend_domain_face, + LayoutData& mc_fabs, bool write_stl_output); +#endif + private: template friend class GShopLevel; @@ -291,7 +334,8 @@ public: * \brief Build a level directly from a GeometryShop object. */ GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& geom, int max_grid_size, - int ngrow, bool extend_domain_face, int num_crse_opt); + int ngrow, bool extend_domain_face, int num_crse_opt, + bool finest_level); /** * \brief Construct a coarse level by coarsening \p fineLevel. */ @@ -305,13 +349,28 @@ public: #if (AMREX_SPACEDIM == 3) /** - * \brief Define a fine level for multi-valued multi-cut STL geometries (3-D only) + * \brief Define data for the finest level with the marching-cubes generator (3-D only). + * + * Works with any geometry source: STLtools supplies its own MC sampling, + * while GeometryShop implicit functions go through the generic adapters in + * AMReX_EB2_Level_MC.H. */ - template - requires (std::same_as) - void define_fine_mvmc (GS const& gshop, const Geometry& geom, - int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt); + void define_fine_marching_cubes (G const& gshop, const Geometry& geom, + int max_grid_size, int ngrow, + bool extend_domain_face, int num_crse_opt, + bool finest_level); #endif + /** + * \brief Define data for a level directly from the geometry source using the + * generator selected by EB2::GetGeometryMethod(). + * + * \p finest_level is false when a coarse level is rebuilt from the + * geometry source rather than by coarsening; it only controls diagnostics + * such as `eb2.mc_stl_file`, which is written for the finest level only. + */ + void define_fine_with_method (G const& gshop, const Geometry& geom, + int max_grid_size, int ngrow, bool extend_domain_face, + int num_crse_opt, bool finest_level); /** * \brief Build a regular (no EB) level. @@ -343,7 +402,8 @@ GShopLevel::GShopLevel (IndexSpace const* is, const Geometry& geom) template GShopLevel::GShopLevel (IndexSpace const* is, G const& gshop, 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 finest_level) : Level(is, geom) { if (std::is_same_v) { @@ -352,7 +412,27 @@ GShopLevel::GShopLevel (IndexSpace const* is, G const& gshop, const Geometry& return; } - define_fine(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); + define_fine_with_method(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt, + finest_level); +} + +template +void +GShopLevel::define_fine_with_method (G const& gshop, const Geometry& geom, int max_grid_size, + int ngrow, bool extend_domain_face, int num_crse_opt, + bool finest_level) +{ + if (GetGeometryMethod() == GeometryMethod::marching_cubes) { +#if (AMREX_SPACEDIM == 3) + define_fine_marching_cubes(gshop, geom, max_grid_size, ngrow, extend_domain_face, + num_crse_opt, finest_level); +#else + amrex::ignore_unused(finest_level); + amrex::Abort("eb2.geometry_method = marching_cubes is only available in 3D"); +#endif + } else { + define_fine(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); + } } template @@ -477,20 +557,10 @@ GShopLevel::define_fine (G const& gshop, const Geometry& geom, { 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); + auto const repair = GetRepairParameters(); + Real const small_volfrac = repair.small_volfrac; + bool const cover_multiple_cuts = repair.cover_multiple_cuts; + int const maxiter = repair.maxiter; prepare_grids(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); @@ -753,40 +823,36 @@ GShopLevel::define_fine (G const& gshop, const Geometry& geom, #if (AMREX_SPACEDIM == 3) template -template -requires (std::same_as) void -GShopLevel::define_fine_mvmc (GS const& gshop, const Geometry& geom, - int max_grid_size, int ngrow, bool extend_domain_face, int num_crse_opt) +GShopLevel::define_fine_marching_cubes (G const& gshop, const Geometry& geom, + int max_grid_size, int ngrow, + bool extend_domain_face, int num_crse_opt, + bool finest_level) { - BL_PROFILE("EB2::GShopLevel()-fine-mvmc"); + BL_PROFILE("EB2::GShopLevel()-fine-marching-cubes"); - prepare_grids(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); + assert_marching_cubes_supported(geom); - if (m_allregular) { return; } - - MC::Initialize(); - - 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(); + // Keep only the cut-box discovery from the legacy path. Every EB field is + // allocated and populated from the MC working level set by the driver. + prepare_grids(gshop, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); + if (m_allregular) { + return; } - 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()]); + MFInfo mf_info; + mf_info.SetTag("EB2::Level-MC"); + m_sdf.define(amrex::convert(m_grids, IntVect::TheNodeVector()), m_dmap, 1, IntVect(1), + mf_info); + + // Geometry-source specific work: nodal MC level set and exact edge crossings. + detail::fill_marching_cubes_levelset(gshop, m_sdf, geom); + LayoutData mc_fabs = define_marching_cubes_caches(); + for (MFIter mfi(m_sdf, MFItInfo().DisableDeviceSync()); mfi.isValid(); ++mfi) { + detail::fill_marching_cubes_edge_intersections(gshop, m_sdf[mfi], geom, mc_fabs[mfi]); } - MC::write_stl("test.stl", m_marching_cubes); - - m_ok = false; // xxxxx TODO + build_marching_cubes_level(geom, extend_domain_face, mc_fabs, finest_level); } #endif @@ -864,9 +930,11 @@ GShopLevel::buildFCData (G const& gshop, int face_dir, int max_grid_size) // Always create transient refined CC level to avoid grid mismatch // with application-defined BoxArrays (ERF EB-FC development) Geometry refined_geom = amrex::refine(m_geom, 2); + // The transient level is not the finest level of the IndexSpace: its + // diagnostics (e.g. eb2.mc_stl_file) must not replace the real ones. transient_level = std::make_unique>( m_parent, gshop, refined_geom, - max_grid_size, m_ngrow[0] > 0 ? m_ngrow[0] : 4, true, 0); // xxxxx TODO: should we use IndexSpace's member variables extend_domain_face and num_crse_opt? + max_grid_size, m_ngrow[0] > 0 ? m_ngrow[0] : 4, true, 0, false); // xxxxx TODO: should we use IndexSpace's member variables extend_domain_face and num_crse_opt? if (!transient_level->isOK()) { amrex::Abort("GShopLevel::buildFCData: failed to build refined level"); diff --git a/Src/EB/AMReX_EB2_Level.cpp b/Src/EB/AMReX_EB2_Level.cpp index a191d2a8b11..d7609a99105 100644 --- a/Src/EB/AMReX_EB2_Level.cpp +++ b/Src/EB/AMReX_EB2_Level.cpp @@ -6,6 +6,28 @@ namespace amrex::EB2 { +RepairParameters +GetRepairParameters () +{ + RepairParameters p; +#ifdef AMREX_USE_FLOAT + p.small_volfrac = 1.e-5_rt; +#else + p.small_volfrac = 1.e-14; +#endif + p.cover_multiple_cuts = false; + p.maxiter = 32; + { + ParmParse pp("eb2"); + pp.queryAdd("small_volfrac", p.small_volfrac); + pp.queryAdd("cover_multiple_cuts", p.cover_multiple_cuts); + pp.queryAdd("maxiter", p.maxiter); + } + p.maxiter = std::min(100000, p.maxiter); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(p.maxiter > 0, "eb2.maxiter must be positive"); + return p; +} + void Level::prepareForCoarsening (const Level& rhs, int max_grid_size, IntVect const& ngrow) { diff --git a/Src/EB/AMReX_EB2_Level_MC.H b/Src/EB/AMReX_EB2_Level_MC.H new file mode 100644 index 00000000000..39a133988c1 --- /dev/null +++ b/Src/EB/AMReX_EB2_Level_MC.H @@ -0,0 +1,196 @@ +#ifndef AMREX_EB2_LEVEL_MC_H_ +#define AMREX_EB2_LEVEL_MC_H_ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +/** + * \file AMReX_EB2_Level_MC.H + * + * Adapters that produce the two geometry-source inputs of the marching-cubes + * EB builder (see AMReX_MarchingCubes.H, "Geometry-source contract"): + * + * - the nodal MC level set (positive in fluid), and + * - the exact Cartesian-edge crossings cached in MC::MCFab. + * + * A geometry source that provides its own MC sampling (STLtools) is used as is. + * Every other source, i.e. any EB2::GeometryShop implicit function, is served + * by the generic adapters below, which reuse the source's `fillFab` and + * `getIntercept` (Brent root finding on cut edges). + * + * These are implementation details of GShopLevel; they are not part of the + * public EB2 interface. + */ + +namespace amrex::EB2::detail { + +//! Geometry sources with dedicated marching-cubes sampling (e.g. STLtools). +template +concept HasMarchingCubesSampling = requires (G const& g, MultiFab& mf, IntVect const& ng, + Geometry const& geom, FArrayBox const& fab, + MC::MCFab& mc_fab) +{ + g.fillMarchingCubesLevelSet(mf, ng, geom); + g.fillMarchingCubesEdgeIntersections(fab, geom, mc_fab); +}; + +//! Execution target used by the generic marching-cubes adapters for \p gshop. +template +RunOn marching_cubes_run_on (G const& gshop) noexcept +{ + return (Gpu::inLaunchRegion() && gshop.isGPUable()) ? RunOn::Gpu : RunOn::Cpu; +} + +/** + * \brief Fill the nodal marching-cubes level set from geometry source \p gshop. + * + * On return every node of \p sdf (including ghost nodes) holds a value that is + * positive in fluid and negative in the covered region. Every node is a + * deterministic function of the geometry alone, so shared and ghost nodes + * agree across FABs without communication. Implicit functions are sampled + * without clamping to the domain because the builder applies its own + * domain-face extension; as in the legacy generator, ghost nodes across a + * periodic boundary sample the function at their unwrapped coordinates. + */ +template +void fill_marching_cubes_levelset (G const& gshop, MultiFab& sdf, Geometry const& geom) +{ + if constexpr (HasMarchingCubesSampling) { + gshop.fillMarchingCubesLevelSet(sdf, sdf.nGrowVect(), geom); + } else { + BL_PROFILE("EB2::fill_marching_cubes_levelset"); + RunOn const run_on = marching_cubes_run_on(gshop); + + // Roundoff can leave a node that lies on the surface at a tiny nonzero + // value instead of exactly zero, which would create a sliver cut cell + // that the builder cannot represent. Snap such nodes to the exact-zero + // ON state. The threshold is relative to the local magnitude of the + // implicit function, i.e. it bounds the crossing position on the + // incident edges as a fraction of the cell size, so it is independent + // of the function's units; it matches the sliver tolerance below which + // the marching-cubes moments reject a cut. +#ifdef AMREX_USE_FLOAT + constexpr Real snap_factor = 1.e-4_rt; +#else + constexpr Real snap_factor = 1.e-11_rt; +#endif + for (MFIter mfi(sdf); mfi.isValid(); ++mfi) { + Box const nbx = sdf[mfi].box(); + // Sample one more ring of nodes so that the snap stencil is + // complete on every node of the FAB, ghost nodes included. + FArrayBox raw(amrex::grow(nbx,1), 1, The_Arena()); + Elixir raw_eli = raw.elixir(); + gshop.fillFab(raw, geom, run_on, raw.box()); + auto const f = raw.const_array(); + auto const a = sdf.array(mfi); + ParallelFor(nbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const v = f(i,j,k); + Real scale = std::abs(v); + scale = amrex::max(scale, std::abs(f(i-1,j,k)), std::abs(f(i+1,j,k))); + scale = amrex::max(scale, std::abs(f(i,j-1,k)), std::abs(f(i,j+1,k))); + scale = amrex::max(scale, std::abs(f(i,j,k-1)), std::abs(f(i,j,k+1))); + // EB2 implicit functions are positive inside the body; marching + // cubes extracts the positive material, i.e. the fluid. + a(i,j,k) = (std::abs(v) <= snap_factor*scale) ? 0.0_rt : -v; + }); + } + Gpu::streamSynchronize(); + } +} + +/** + * \brief Cache exact edge crossings for one FAB from geometry source \p gshop. + * + * \p sdf is the FAB's nodal marching-cubes level set (positive in fluid). For + * every edge whose end points differ in sign, the crossing position from the + * low node, normalized by the cell size, is stored in + * \p mc_fab.m_edge_intersections; other edges keep MC::invalid_edge_intersection. + */ +template +void fill_marching_cubes_edge_intersections (G const& gshop, FArrayBox const& sdf_fab, + Geometry const& geom, MC::MCFab& mc_fab) +{ + if constexpr (HasMarchingCubesSampling) { + gshop.fillMarchingCubesEdgeIntersections(sdf_fab, geom, mc_fab); + } else { + BL_PROFILE("EB2::fill_marching_cubes_edge_intersections"); + RunOn const run_on = marching_cubes_run_on(gshop); + auto const sdf = sdf_fab.const_array(); + auto const problo = geom.ProbLoArray(); + auto const dx = geom.CellSizeArray(); + + Array,AMREX_SPACEDIM> type_fabs; + Array,AMREX_SPACEDIM> inter_fabs; + Array,AMREX_SPACEDIM> type_arrs; + Array,AMREX_SPACEDIM> inter_arrs; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + Box const ebx = mc_fab.m_edge_intersections[idim].box(); + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(amrex::surroundingNodes(ebx))); + type_fabs[idim].resize(ebx, 1); + inter_fabs[idim].resize(ebx, 1); + IntVect const shift = IntVect::TheDimensionVector(idim); + auto const type = type_fabs[idim].array(); + ParallelFor(ebx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const lo = sdf(i,j,k); + Real const hi = sdf(i+shift[0],j+shift[1],k+shift[2]); + // An exact-zero (ON) end point already is the crossing; the + // builder places it at that node without an exact intercept. + // Root finding on such an edge is not needed and could fail + // when the zero came from snapping a roundoff-level value. + // Away from snapped nodes the field is the negated implicit + // function, so a strict sign change brackets the root for + // getIntercept. + bool const strict_sign_change = (lo > 0.0_rt && hi < 0.0_rt) + || (lo < 0.0_rt && hi > 0.0_rt); + type(i,j,k) = strict_sign_change ? Type::irregular : Type::regular; + }); + type_arrs[idim] = type_fabs[idim].const_array(); + inter_arrs[idim] = inter_fabs[idim].array(); + } + Gpu::streamSynchronize(); + + // Physical crossing coordinates on irregular edges; NaN elsewhere. + // GeometryShop ignores the level-set argument (it is only a hint for + // STL ray casting), so the MC-convention values are passed unchanged. + gshop.getIntercept(inter_arrs, type_arrs, sdf, geom, run_on, sdf_fab.box()); + + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + Box const ebx = mc_fab.m_edge_intersections[idim].box(); + auto const type = type_fabs[idim].const_array(); + auto const inter = inter_fabs[idim].const_array(); + auto const crossing = mc_fab.m_edge_intersections[idim].array(); + Real const problo_d = problo[idim]; + Real const dx_d = dx[idim]; + ParallelFor(ebx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + if (type(i,j,k) == Type::irregular) { + Real const x = inter(i,j,k); + if (amrex::isnan(x)) { + crossing(i,j,k) = MC::invalid_edge_intersection; + } else { + int const iv[3] = {i,j,k}; + Real const lo = problo_d + static_cast(iv[idim])*dx_d; + crossing(i,j,k) = amrex::Clamp((x-lo)/dx_d, 0.0_rt, 1.0_rt); + } + } + }); + } + Gpu::streamSynchronize(); + } +} + +} // namespace amrex::EB2::detail + +#endif diff --git a/Src/EB/AMReX_EB2_Level_MC.cpp b/Src/EB/AMReX_EB2_Level_MC.cpp new file mode 100644 index 00000000000..fa62fd6afca --- /dev/null +++ b/Src/EB/AMReX_EB2_Level_MC.cpp @@ -0,0 +1,386 @@ +#include +#include +#include +#include + +#include +#include +#include +#include + +/** + * \file AMReX_EB2_Level_MC.cpp + * + * Geometry-source independent driver of the marching-cubes EB builder. The + * only geometry-specific work (filling the nodal MC level set and the exact + * edge crossings) happens in GShopLevel::define_fine_marching_cubes; see + * AMReX_EB2_Level_MC.H for the adapters. Everything below consumes just the + * nodal field in m_sdf and the per-FAB MC::MCFab caches. + */ + +namespace amrex::EB2 { + +namespace { + +// Device kernels live in free functions: CUDA does not allow extended device +// lambdas inside protected member functions. + +//! Pre-fill the volume fraction of \p bx from the corner signs: 0 covered, +//! 1 regular, -1 for cut cells that build_cell_fractions overwrites. +void prefill_volume_fractions (Box const& bx, Array4 const& sdf, + Array4 const& vfrac) +{ + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int nfluid = 0; + nfluid += sdf(i, j, k) > 0.0_rt; + nfluid += sdf(i + 1, j, k) > 0.0_rt; + nfluid += sdf(i, j + 1, k) > 0.0_rt; + nfluid += sdf(i + 1, j + 1, k) > 0.0_rt; + nfluid += sdf(i, j, k + 1) > 0.0_rt; + nfluid += sdf(i + 1, j, k + 1) > 0.0_rt; + nfluid += sdf(i, j + 1, k + 1) > 0.0_rt; + nfluid += sdf(i + 1, j + 1, k + 1) > 0.0_rt; + vfrac(i, j, k) = nfluid == 0 ? 0.0_rt : (nfluid == 8 ? 1.0_rt : -1.0_rt); + }); +} + +//! Flip the sign of the nodal field: EB2's public convention is negative in fluid. +void negate_levelset (Box const& bx, Array4 const& phi) +{ + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + phi(i, j, k) = -phi(i, j, k); + }); +} + +//! Per-FAB counter blocks (MC::num_fab_counters ints each) for one level. +class FabCounters +{ +public: + explicit FabCounters (int nfabs) + : m_buffer(static_cast(nfabs) * MC::num_fab_counters) + {} + + //! Zero every block on the host and push to the device. + void reset () + { + std::fill_n(m_buffer.hostData(), m_buffer.size(), 0); + m_buffer.copyToDeviceAsync(); + // The FAB loops that follow may use other streams. + Gpu::streamSynchronize(); + } + + //! Device pointer to the block of the FAB with local index \p local_index. + int* block (int local_index) noexcept + { + return m_buffer.data() + static_cast(local_index) * MC::num_fab_counters; + } + + //! Copy every block to the host and reduce over FABs and MPI ranks. + void reduce () + { + Gpu::streamSynchronizeAll(); + m_buffer.copyToHost(); + m_totals.fill(0); + for (std::size_t n = 0; n < m_buffer.size(); ++n) { + m_totals[n % MC::num_fab_counters] += m_buffer.hostData()[n]; + } + ParallelAllReduce::Sum(m_totals.data(), MC::num_fab_counters, + ParallelContext::CommunicatorSub()); + } + + //! Level-wide total of \p counter after reduce(). + [[nodiscard]] int total (MC::FabCounter counter) const noexcept { return m_totals[counter]; } + +private: + Gpu::Buffer m_buffer; + std::array m_totals{}; +}; + +} // namespace + +void +Level::assert_marching_cubes_supported (Geometry const& geom) +{ + auto const cell_size = geom.CellSizeArray(); + Real const max_cell_size = amrex::max(cell_size[0], amrex::max(cell_size[1], cell_size[2])); + Real const cubic_tolerance = 16.0_rt * std::numeric_limits::epsilon() * max_cell_size; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + geom.Coord() == CoordSys::cartesian && + std::abs(cell_size[0] - cell_size[1]) <= cubic_tolerance && + std::abs(cell_size[0] - cell_size[2]) <= cubic_tolerance, + "Marching-cubes EB construction requires a 3D Cartesian grid with dx == " + "dy == dz"); +} + +LayoutData +Level::define_marching_cubes_caches () +{ + BL_PROFILE("EB2::Level::define_marching_cubes_caches"); + + MC::Initialize(); + + LayoutData mc_fabs(m_grids, m_dmap); + for (MFIter mfi(m_sdf, MFItInfo().DisableDeviceSync()); mfi.isValid(); ++mfi) { + mc_fabs[mfi].defineEdgeIntersections(m_sdf[mfi].box()); + } + return mc_fabs; +} + +void +Level::build_marching_cubes_level (Geometry const& geom, bool extend_domain_face, + LayoutData& mc_fabs, bool write_stl_output) +{ + BL_PROFILE("EB2::Level::build_marching_cubes_level"); + + auto const repair = GetRepairParameters(); + Real const small_volfrac = repair.small_volfrac; + int const maxiter = repair.maxiter; + bool const cover_multiple_cuts = repair.cover_multiple_cuts; + + // The geometry source has filled m_sdf (MC convention: > 0 fluid) and the + // exact edge crossings from the raw geometry. The domain-face extension + // is owned by the builder so that every geometry source behaves + // identically: it is applied to the nodal field and to the crossings. + // It runs on whole FABs, ghost nodes included, so every node stays a + // deterministic function of the geometry and no communication is needed. + Box const domain = geom.Domain(); + GpuArray const is_periodic{geom.isPeriodic(0), geom.isPeriodic(1), geom.isPeriodic(2)}; + // Every MC entry point accumulates its counts into the FAB's device + // counter block; the blocks are copied to the host once per pass. + FabCounters counters(m_sdf.local_size()); + counters.reset(); + if (extend_domain_face) { + for (MFIter mfi(m_sdf, MFItInfo().DisableDeviceSync()); mfi.isValid(); ++mfi) { + MC::extend_domain_face_levelset(m_sdf[mfi].box(), domain, is_periodic, m_sdf[mfi], + counters.block(mfi.LocalIndex())); + MC::extend_domain_face_edge_intersections(domain, is_periodic, mc_fabs[mfi]); + } + } + + // Every EB field is allocated here and populated from the MC working SDF. + int const ng = GFab::ng; + MFInfo mf_info; + mf_info.SetTag("EB2::Level-MC"); + m_cellflag.define(m_grids, m_dmap, 1, ng, mf_info); + m_volfrac.define(m_grids, m_dmap, 1, ng, mf_info); + m_centroid.define(m_grids, m_dmap, AMREX_SPACEDIM, ng, mf_info); + m_bndryarea.define(m_grids, m_dmap, 1, ng, mf_info); + m_bndrycent.define(m_grids, m_dmap, AMREX_SPACEDIM, ng, mf_info); + m_bndrynorm.define(m_grids, m_dmap, AMREX_SPACEDIM, ng, mf_info); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + m_areafrac[idim].define(amrex::convert(m_grids, IntVect::TheDimensionVector(idim)), m_dmap, + 1, ng, mf_info); + m_facecent[idim].define(amrex::convert(m_grids, IntVect::TheDimensionVector(idim)), m_dmap, + AMREX_SPACEDIM - 1, ng, mf_info); + IntVect edge_type{1}; + edge_type[idim] = 0; + m_edgecent[idim].define(amrex::convert(m_grids, edge_type), m_dmap, 1, ng, mf_info); + } + // prefill_volume_fractions rewrites grow(vbx,1) every pass; this only + // provides the regular default beyond it. + m_volfrac.setVal(1.0_rt, 0, 1, m_volfrac.nGrowVect()); + + // Like the moments, the rejection marks are computed on grow(vbx,1), so + // the ghost ring of these arrays holds exactly the marks the neighboring + // FAB computes for its valid cells and faces, and the nodal repair needs + // no communication. Degenerate faces in the outer ring are marked and + // repaired like any other. + iMultiFab rejected_cells(m_grids, m_dmap, 1, IntVect(1), mf_info); + Array rejected_faces; + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + rejected_faces[idim].define(m_areafrac[idim].boxArray(), m_dmap, 1, IntVect(1), mf_info); + } + bool converged = false; + for (int iter = 0; iter < maxiter; ++iter) { + rejected_cells.setVal(0); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + rejected_faces[idim].setVal(0); + } + + counters.reset(); + 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) { + MC::marching_cubes(geom, m_sdf[mfi], mc_fabs[mfi], counters.block(mfi.LocalIndex())); + } + + m_centroid.setVal(0.0_rt, 0, AMREX_SPACEDIM, m_centroid.nGrowVect()); + m_bndryarea.setVal(0.0_rt, 0, 1, m_bndryarea.nGrowVect()); + m_bndrycent.setVal(-1.0_rt, 0, AMREX_SPACEDIM, m_bndrycent.nGrowVect()); + m_bndrynorm.setVal(0.0_rt, 0, AMREX_SPACEDIM, m_bndrynorm.nGrowVect()); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(m_sdf, info); mfi.isValid(); ++mfi) { + Box const vbx = amrex::enclosedCells(mfi.validbox()); + Box const gbx = amrex::grow(vbx, 1); + auto& mc_fab = mc_fabs[mfi]; + int* const fab_counters = counters.block(mfi.LocalIndex()); + + prefill_volume_fractions(gbx, m_sdf.const_array(mfi), m_volfrac.array(mfi)); + + MC::build_face_fractions( + gbx, mc_fab, m_sdf[mfi], m_areafrac[0][mfi], m_areafrac[1][mfi], + m_areafrac[2][mfi], m_facecent[0][mfi], m_facecent[1][mfi], m_facecent[2][mfi], + rejected_faces[0][mfi], rejected_faces[1][mfi], rejected_faces[2][mfi], + fab_counters); + MC::build_edge_centroids(gbx, mc_fab, m_sdf[mfi], m_edgecent[0][mfi], + m_edgecent[1][mfi], m_edgecent[2][mfi]); + // Moments and rejection marks include the ghost ring so that they + // are available locally for the repair and the final topology. + MC::build_cell_fractions( + gbx, geom, mc_fab, m_sdf[mfi], m_areafrac[0][mfi], + m_areafrac[1][mfi], m_areafrac[2][mfi], m_volfrac[mfi], m_centroid[mfi], + m_bndryarea[mfi], m_bndrycent[mfi], m_bndrynorm[mfi], fab_counters); + + MC::mark_faces_for_cleanup( + gbx, mc_fab, m_sdf[mfi], rejected_faces[0][mfi], + rejected_faces[1][mfi], rejected_faces[2][mfi], fab_counters); + MC::mark_cells_for_cleanup(gbx, mc_fab, m_sdf[mfi], m_volfrac[mfi], small_volfrac, + rejected_cells[mfi], fab_counters); + } + + // One device-to-host copy and one reduction per pass. + counters.reduce(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(counters.total(MC::counter_invalid_triangles) == 0, + "Marching Cubes: invalid triangle"); + int const face_decision_errors = counters.total(MC::counter_face_decision_errors); + int const degenerate_faces = counters.total(MC::counter_degenerate_faces); + int const cell_geometry_errors = counters.total(MC::counter_closure_errors) + + counters.total(MC::counter_volume_errors) + + counters.total(MC::counter_centroid_errors) + + counters.total(MC::counter_area_vector_errors); + int face_rejections = counters.total(MC::counter_face_rejections); + int const topology_rejections = counters.total(MC::counter_topology_rejections); + int const small_cell_rejections = counters.total(MC::counter_small_cell_rejections); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(face_decision_errors == 0, + "Marching-cubes EB encountered " + + std::to_string(face_decision_errors) + + " Cartesian faces whose two cells resolved the " + "MC33 ambiguity differently"); + // Degenerate face polygons were marked for nodal repair alongside the + // multi-aperture faces. + face_rejections += degenerate_faces; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + cell_geometry_errors == 0 || topology_rejections + face_rejections > 0, + "Marching-cubes EB could not map an invalid cell-moment record to " + "the nodal repair set"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + cover_multiple_cuts || face_rejections + topology_rejections == 0, + "Marching-cubes EB found " + std::to_string(face_rejections) + + " unsupported faces and " + std::to_string(topology_rejections) + + " unsupported cells; " + "set eb2.cover_multiple_cuts=1 to enable legacy-style nodal " + "repair"); + + if (face_rejections + topology_rejections + small_cell_rejections == 0) { + converged = true; + break; + } + + counters.reset(); +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(m_sdf, info); mfi.isValid(); ++mfi) { + // Ghost nodes are repaired too, like the nodes of the legacy + // generator's grown level set: ghost nodes beyond every grid have + // no owner and must be fixed here for the repair to converge. + MC::zero_nodes_for_cleanup( + m_sdf[mfi].box(), rejected_cells[mfi], rejected_faces[0][mfi], + rejected_faces[1][mfi], rejected_faces[2][mfi], m_sdf[mfi], + counters.block(mfi.LocalIndex())); + } + counters.reduce(); + int const changed_nodes = counters.total(MC::counter_changed_nodes); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(changed_nodes > 0, + "Marching-cubes EB repair found rejected " + "geometry but changed no fluid nodes"); + + // Shared valid nodes were repaired from identical marks on every FAB + // and already agree. Owned ghost nodes also depend on cells outside + // grow(vbx,1), so refresh them from their owners, as the legacy + // generator does after each pass. The domain-face extension follows + // so that it reads up-to-date reference nodes on every FAB. + m_sdf.FillBoundary(geom.periodicity()); + if (extend_domain_face) { +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(m_sdf, info); mfi.isValid(); ++mfi) { + MC::extend_domain_face_levelset(m_sdf[mfi].box(), domain, is_periodic, + m_sdf[mfi], counters.block(mfi.LocalIndex())); + } + } + if (amrex::Verbose() > 0) { + if (small_cell_rejections > 0) { + amrex::Print() << "AMReX MC EB: Iter. " << iter + 1 << " fixed " + << small_cell_rejections << " small cells\n"; + } + if (face_rejections + topology_rejections > 0) { + amrex::Print() << "AMReX MC EB: Iter. " << iter + 1 << " fixed " + << face_rejections << " unsupported faces and " + << topology_rejections << " unsupported cells\n"; + } + if (cell_geometry_errors > 0) { + amrex::Print() << "AMReX MC EB: Iter. " << iter + 1 << " included " + << cell_geometry_errors + << " invalid cell-moment records in the repair set (closure=" + << counters.total(MC::counter_closure_errors) + << ", volume=" << counters.total(MC::counter_volume_errors) + << ", centroid=" << counters.total(MC::counter_centroid_errors) + << ", area-vector=" + << counters.total(MC::counter_area_vector_errors) << ")\n"; + } + } + } + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(converged, + "Marching-cubes EB failed to fix small or unsupported cells"); + + MFItInfo final_info{}; +#if defined(AMREX_USE_OMP) && !defined(AMREX_USE_GPU) + final_info.SetDynamic(true); +#endif + m_cellflag.setVal(EBCellFlag::TheDefaultCell()); + int geometry_errors = 0; +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) reduction(+ : geometry_errors) +#endif + for (MFIter mfi(m_sdf, final_info); mfi.isValid(); ++mfi) { + Box const vbx = amrex::enclosedCells(mfi.validbox()); + geometry_errors += MC::build_cell_topology( + vbx, mc_fabs[mfi], m_sdf[mfi], m_cellflag[mfi], m_volfrac[mfi], + m_areafrac[0][mfi], m_areafrac[1][mfi], m_areafrac[2][mfi]); + m_cellflag[mfi].resetType(GFab::ng); + } + ParallelAllReduce::Sum(geometry_errors, ParallelContext::CommunicatorSub()); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(geometry_errors == 0, + "Final repaired marching-cubes topology is inconsistent"); + + // The converged triangulation and the public EB data now describe the same + // repaired domain. Only the finest level is written; coarse levels that + // are rebuilt from the geometry source would otherwise overwrite it. + std::string stl_output; + ParmParse("eb2").query("mc_stl_file", stl_output); + if (write_stl_output && !stl_output.empty()) { + MC::write_stl(stl_output, mc_fabs); + } + mc_fabs.clear(); + +#ifdef AMREX_USE_OMP +#pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi(m_sdf, final_info); mfi.isValid(); ++mfi) { + negate_levelset(m_sdf[mfi].box(), m_sdf.array(mfi)); + } + m_levelset = std::move(m_sdf); + + m_ok = true; +} + +} diff --git a/Src/EB/AMReX_EB2_Level_STL.H b/Src/EB/AMReX_EB2_Level_STL.H index 6052e886359..f1762dc13f5 100644 --- a/Src/EB/AMReX_EB2_Level_STL.H +++ b/Src/EB/AMReX_EB2_Level_STL.H @@ -24,7 +24,7 @@ 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, bool support_mvmc); + int num_crse_opt, bool support_mvmc, bool finest_level); //! Coarsen from an existing fine-level STL description. STLLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow, diff --git a/Src/EB/AMReX_EB2_Level_STL.cpp b/Src/EB/AMReX_EB2_Level_STL.cpp index 615a3321d19..d4fe056888f 100644 --- a/Src/EB/AMReX_EB2_Level_STL.cpp +++ b/Src/EB/AMReX_EB2_Level_STL.cpp @@ -4,27 +4,21 @@ 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, - bool support_mvmc) + bool support_mvmc, bool finest_level) : GShopLevel(is, geom) { BL_PROFILE("EB2::STLLevel()-fine"); -#if (AMREX_SPACEDIM == 3) - bool test_marching_cubes = false; - { - ParmParse pp("eb2"); - pp.query("test_marching_cubes", test_marching_cubes); - } - if (support_mvmc && test_marching_cubes) { - define_fine_mvmc(stl_tools, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); - } else -#endif - { - if (amrex::Verbose() && support_mvmc) { - amrex::Warning("STLlevel: support_mvmc = true is not supported yet"); - } - define_fine(stl_tools, geom, max_grid_size, ngrow, extend_domain_face, num_crse_opt); + if (GetGeometryMethod() == GeometryMethod::marching_cubes) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !support_mvmc, + "The supported marching-cubes STL generator is single-valued; " + "use ordinary EB2::Build instead of BuildMultiValuedMultiCut"); + } else if (amrex::Verbose() && support_mvmc) { + amrex::Warning("STLlevel: support_mvmc = true is not supported yet"); } + define_fine_with_method(stl_tools,geom,max_grid_size,ngrow,extend_domain_face,num_crse_opt, + finest_level); } STLLevel::STLLevel (IndexSpace const* is, int ilev, int max_grid_size, int ngrow, diff --git a/Src/EB/AMReX_EBToPVD.H b/Src/EB/AMReX_EBToPVD.H index b8eef52cd90..fc68b3b092e 100644 --- a/Src/EB/AMReX_EBToPVD.H +++ b/Src/EB/AMReX_EBToPVD.H @@ -8,6 +8,7 @@ #include #include #include +#include /** * \file AMReX_EBToPVD.H @@ -40,6 +41,8 @@ public: //! Write the accumulated EB polygons for rank \p myID into a VTP file. void WriteEBVTP(int myID) const; + //! Triangulate and write the accumulated planar EB polygons to one ASCII STL file. + void WriteSTL(std::string const& filename) const; //! Write a PVD file referencing per-rank outputs. static void WritePVTP(int nProcs) ; diff --git a/Src/EB/AMReX_EBToPVD.cpp b/Src/EB/AMReX_EBToPVD.cpp index 959c85aa89e..e1b23a1af28 100644 --- a/Src/EB/AMReX_EBToPVD.cpp +++ b/Src/EB/AMReX_EBToPVD.cpp @@ -1,9 +1,11 @@ #include #include #include +#include #include #include +#include #include #include #include @@ -215,6 +217,85 @@ void EBToPVD::WriteEBVTP(const int myID) const } } +void EBToPVD::WriteSTL (std::string const& filename) const +{ + int const myproc = ParallelContext::MyProcSub(); + int const nprocs = ParallelContext::NProcsSub(); + std::ofstream ofs; + + if (myproc == 0) { + ofs.open(filename); + ofs << "solid Created by AMReX planar EB reconstruction\n"; + } + +#ifdef AMREX_USE_MPI + if (myproc > 0) { + int token = 0; + ParallelDescriptor::Recv(&token, 1, myproc-1, 101, ParallelContext::CommunicatorSub()); + } +#endif + + if (!ofs.is_open()) { + ofs.open(filename, std::ios_base::app); + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ofs.good(), "Could not open planar EB STL output " + filename); + ofs << std::setprecision(std::numeric_limits::max_digits10); + + for (auto const& polygon : m_connectivity) { + int const npoints = polygon[0]; + AMREX_ALWAYS_ASSERT(npoints <= 6); + if (npoints < 3) { continue; } // no polygon to triangulate + auto const& v1 = m_points[polygon[1]]; + for (int n = 2; n < npoints; ++n) { + auto const& v2 = m_points[polygon[n]]; + auto const& v3 = m_points[polygon[n+1]]; + std::array const e1{ + v2[0]-v1[0],v2[1]-v1[1],v2[2]-v1[2]}; + std::array const e2{ + v3[0]-v1[0],v3[1]-v1[1],v3[2]-v1[2]}; + std::array normal{ + e1[1]*e2[2]-e1[2]*e2[1], + e1[2]*e2[0]-e1[0]*e2[2], + e1[0]*e2[1]-e1[1]*e2[0]}; + Real const norm = std::sqrt(dot_product(normal,normal)); + Real const edge_scale_sq = std::max( + dot_product(e1,e1),dot_product(e2,e2)); + Real const degenerate_tolerance = 64.0_rt + * std::numeric_limits::epsilon()*edge_scale_sq; + if (norm <= degenerate_tolerance) { + continue; + } + for (Real& value : normal) { + value /= norm; + } + ofs << "facet normal " + << normal[0] << " " << normal[1] << " " << normal[2] << "\n" + << " outer loop\n" + << " vertex " << v1[0] << " " << v1[1] << " " << v1[2] << "\n" + << " vertex " << v2[0] << " " << v2[1] << " " << v2[2] << "\n" + << " vertex " << v3[0] << " " << v3[1] << " " << v3[2] << "\n" + << " endloop\n" + << "endfacet\n"; + } + } + + if (myproc == nprocs-1) { + ofs << "endsolid Created by AMReX planar EB reconstruction\n"; + } + ofs.close(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !ofs.fail(), "Could not complete planar EB STL output " + filename); + +#ifdef AMREX_USE_MPI + if (myproc < nprocs-1) { + int token = 0; + ParallelDescriptor::Send(&token, 1, myproc+1, 101, ParallelContext::CommunicatorSub()); + } + ParallelDescriptor::Barrier(ParallelContext::CommunicatorSub()); +#endif +} + void EBToPVD::WritePVTP(const int nProcs) { std::ofstream myfile("eb.pvtp"); diff --git a/Src/EB/AMReX_EB_STL_utils.H b/Src/EB/AMReX_EB_STL_utils.H index fbd433e8059..89d5ec8f241 100644 --- a/Src/EB/AMReX_EB_STL_utils.H +++ b/Src/EB/AMReX_EB_STL_utils.H @@ -21,6 +21,10 @@ namespace amrex { +namespace MC { +struct MCFab; +} + //! Utility class that loads STL meshes and samples them onto AMReX grids. class STLtools { @@ -75,6 +79,12 @@ public: //! Enable/disable BVH acceleration when sampling. void setBVHOptimization (bool flag) { m_bvh_optimization = flag; } + //! Select the strict STL contract required by the marching-cubes builder. + void setUseMarchingCubes (bool flag) noexcept { m_use_marching_cubes = flag; } + [[nodiscard]] bool usesMarchingCubes () const noexcept { + return m_use_marching_cubes; + } + /** * \brief Read an STL file, apply scaling/translation, and store triangles. * @@ -115,12 +125,34 @@ public: //! Fill \p mf with signed-distance values sampled from the STL mesh. void fillSignedDistance (MultiFab& mf, IntVect const& nghost, Geometry const& geom) const; + /** + * \brief Fill the nodal scalar field used by the marching-cubes builder. + * + * Exact STL distances are evaluated only in the two-node band needed by + * cut-cell interpolation and vertex-normal stencils. Away from that + * band, the field retains unit-magnitude inside/outside values. + */ + void fillMarchingCubesLevelSet (MultiFab& mf, IntVect const& nghost, + Geometry const& geom) const; + + /** + * \brief Cache exact normalized Cartesian-edge crossings for MC geometry. + * + * The cache is evaluated against the original triangles rather than the + * sampled signed-distance field. It remains valid through the monotone + * nodal cleanup used by marching cubes. + */ + void fillMarchingCubesEdgeIntersections ( + FArrayBox const& levelset, Geometry const& geom, + MC::MCFab& mc_fab) const; + //! Upload triangle data to device memory (exposed for CUDA workflows). void prepare (Gpu::PinnedVector a_tri_pts); // public for cuda private: bool m_bvh_optimization = true; + bool m_use_marching_cubes = false; Gpu::DeviceVector m_tri_pts_d; Gpu::DeviceVector m_tri_normals_d; diff --git a/Src/EB/AMReX_EB_STL_utils.cpp b/Src/EB/AMReX_EB_STL_utils.cpp index 98417b2adbb..ecf7477798a 100644 --- a/Src/EB/AMReX_EB_STL_utils.cpp +++ b/Src/EB/AMReX_EB_STL_utils.cpp @@ -2,11 +2,18 @@ #include #include #include +#include #include #include +#include +#include #include +#include +#include #include +#include +#include // Reference for BVH: https://rmrsk.github.io/EBGeometry/Concepts.html#bounding-volume-hierarchies @@ -15,6 +22,145 @@ namespace amrex namespace { + struct STLVertexKey { + Real x; + Real y; + Real z; + + bool operator== (STLVertexKey const& rhs) const noexcept { + return x == rhs.x && y == rhs.y && z == rhs.z; + } + }; + + struct STLVertexKeyHash { + std::size_t operator() (STLVertexKey const& key) const noexcept { + std::hash const h; + std::size_t seed = h(key.x); + seed ^= h(key.y) + std::size_t(0x9e3779b97f4a7c15ULL) + (seed << 6) + (seed >> 2); + seed ^= h(key.z) + std::size_t(0x9e3779b97f4a7c15ULL) + (seed << 6) + (seed >> 2); + return seed; + } + }; + + struct STLEdgeKeyHash { + std::size_t operator() (std::pair const& key) const noexcept { + return (static_cast(static_cast(key.first)) << 32) + ^ static_cast(static_cast(key.second)); + } + }; + + void validate_marching_cubes_stl ( + Gpu::PinnedVector const& triangles) + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + !triangles.empty(), + "Marching-cubes STL contains no triangles"); + Real lo[3] = {std::numeric_limits::max(), + std::numeric_limits::max(), + std::numeric_limits::max()}; + Real hi[3] = {std::numeric_limits::lowest(), + std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + for (std::size_t n = 0; n < triangles.size(); ++n) { + auto const& tri = triangles[n]; + std::array const vertices{tri.v1,tri.v2,tri.v3}; + for (auto const& vertex : vertices) { + Real const values[3] = {vertex.x,vertex.y,vertex.z}; + for (int d = 0; d < 3; ++d) { + if (!std::isfinite(values[d])) { + amrex::Abort( + "Marching-cubes STL contains a non-finite vertex in triangle " + + std::to_string(n)); + } + lo[d] = std::min(lo[d],values[d]); + hi[d] = std::max(hi[d],values[d]); + } + } + } + + Real const scale = std::max( + {hi[0]-lo[0],hi[1]-lo[1],hi[2]-lo[2]}); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + scale > 0.0_rt, + "Marching-cubes STL has a zero-size bounding box"); + Real const epsilon = std::numeric_limits::epsilon(); + // Squaring the coordinate resolution gives the smallest meaningful + // cross-product scale. A larger multiplier rejects valid, very small + // facets after otherwise representable single-precision conversion. + Real const area_tolerance = epsilon*epsilon*scale*scale; + + // Vertices are welded only when their coordinates are bit-identical + // after scaling and centering; a file whose shared vertices differ in + // the last digit is reported as open below, with the coordinates. + std::unordered_map vertex_ids; + std::unordered_map,std::pair,STLEdgeKeyHash> edges; + vertex_ids.reserve(triangles.size()*3); + edges.reserve(triangles.size()*3); + Vector vertex_coordinates; + vertex_coordinates.reserve(triangles.size()*3); + int next_vertex_id = 0; + auto vertex_id = [&] (XDim3 const& vertex) -> int { + STLVertexKey const key{.x = vertex.x, .y = vertex.y, .z = vertex.z}; + auto [it,inserted] = vertex_ids.emplace(key,next_vertex_id); + if (inserted) { + ++next_vertex_id; + vertex_coordinates.push_back(vertex); + } + return it->second; + }; + + for (std::size_t n = 0; n < triangles.size(); ++n) { + auto const& tri = triangles[n]; + Real const e1[3] = {tri.v2.x-tri.v1.x, + tri.v2.y-tri.v1.y, + tri.v2.z-tri.v1.z}; + Real const e2[3] = {tri.v3.x-tri.v1.x, + tri.v3.y-tri.v1.y, + tri.v3.z-tri.v1.z}; + Real const cross[3] = { + e1[1]*e2[2]-e1[2]*e2[1], + e1[2]*e2[0]-e1[0]*e2[2], + e1[0]*e2[1]-e1[1]*e2[0]}; + Real const twice_area = std::sqrt( + cross[0]*cross[0]+cross[1]*cross[1]+cross[2]*cross[2]); + if (!(twice_area > area_tolerance)) { + amrex::Abort( + "Marching-cubes STL contains a degenerate triangle at index " + + std::to_string(n)); + } + + int const ids[3] = { + vertex_id(tri.v1),vertex_id(tri.v2),vertex_id(tri.v3)}; + for (int edge = 0; edge < 3; ++edge) { + int const from = ids[edge]; + int const to = ids[(edge+1)%3]; + std::pair const key{ + std::min(from,to),std::max(from,to)}; + auto& record = edges[key]; + ++record.first; + record.second += from < to ? 1 : -1; + } + } + + for (auto const& [edge,record] : edges) { + if (record.first != 2 || record.second != 0) { + auto const& a = vertex_coordinates[edge.first]; + auto const& b = vertex_coordinates[edge.second]; + std::ostringstream message; + message << std::setprecision(std::numeric_limits::max_digits10) + << "Marching-cubes STL is open, nonmanifold, or " + << "inconsistently oriented at the edge from (" + << a.x << ',' << a.y << ',' << a.z << ") to (" + << b.x << ',' << b.y << ',' << b.z << "): incidence=" + << record.first << " (expected 2), orientation_sum=" + << record.second << " (expected 0). Vertices are matched " + << "exactly; repair or re-export the STL so shared vertices " + << "have identical coordinates."; + amrex::Abort(message.str()); + } + } + } + AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE XDim3 triangle_norm (STLtools::Triangle const& tri) { @@ -352,6 +498,9 @@ STLtools::read_stl_file (std::string const& fname, Real scale, Array con } else { read_ascii_stl_file(fname, scale, center, reverse_normal, tri_pts); } + if (m_use_marching_cubes) { + validate_marching_cubes_stl(tri_pts); + } } prepare(std::move(tri_pts)); @@ -706,8 +855,22 @@ STLtools::build_bvh (Triangle* begin, Triangle* end, Gpu::PinnedVector& bv centmax.max(cent); } int max_dir = (centmax-centmin).maxDir(false); - std::sort(begin, end, [max_dir] (Triangle const& a, Triangle const& b) -> bool - { return a.cent(max_dir) < b.cent(max_dir); }); + { + // Sort by precomputed centroid keys. Recomputing Triangle::cent() + // inside the comparator is not a strict weak ordering under + // fast-math (the sum can be re-associated differently at different + // call sites), which lets std::sort run out of bounds. + Vector> keys(ntri); + for (int i = 0; i < ntri; ++i) { + keys[i] = std::make_pair(begin[i].cent(max_dir), i); + } + std::sort(keys.begin(), keys.end()); + Vector sorted(ntri); + for (int i = 0; i < ntri; ++i) { + sorted[i] = begin[keys[i].second]; + } + std::copy(sorted.begin(), sorted.end(), begin); + } int nsplits = std::min((ntri + (m_bvh_max_size-1)) / m_bvh_max_size, m_bvh_max_splits); int tsize = ntri / nsplits; @@ -1376,4 +1539,283 @@ STLtools::fillSignedDistance (MultiFab& mf, IntVect const& nghost, Geometry cons #endif } +void STLtools::fillMarchingCubesLevelSet (MultiFab& mf, IntVect const& nghost, + Geometry const& geom) const +{ + BL_PROFILE("STLtools::fillMarchingCubesLevelSet"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(AMREX_SPACEDIM == 3, + "STLtools::fillMarchingCubesLevelSet is only available in 3D"); + +#if (AMREX_SPACEDIM != 3) + amrex::ignore_unused(this, mf, nghost, geom); +#else + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + nghost.allGE(IntVect(1)) && mf.nGrowVect().allGE(nghost), + "Marching-cubes STL sampling requires at least one nodal ghost cell"); + + // The sign is sufficient in regular/covered regions. Canonicalize it + // before constructing the band so every FAB makes the same decision at a + // shared node. + this->fill(mf, nghost, geom, 1.0_rt, -1.0_rt); + mf.OverrideSync(geom.periodicity()); + mf.FillBoundary(geom.periodicity()); + + FabArray> exact_band( + mf.boxArray(), mf.DistributionMap(), 1, nghost); + + // A radius-one sign search includes every corner of a mixed cell. A byte + // mask preserves the level-wide kernel batching while using one quarter of + // the old integer-mask storage. + + for (MFIter mfi(mf); mfi.isValid(); ++mfi) { + Box const bx = mf[mfi].box(); + int const ilo = bx.smallEnd(0); + int const jlo = bx.smallEnd(1); + int const klo = bx.smallEnd(2); + int const ihi = bx.bigEnd(0); + int const jhi = bx.bigEnd(1); + int const khi = bx.bigEnd(2); + auto const phi = mf.const_array(mfi); + auto const band = exact_band.array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool const fluid = phi(i, j, k) > 0.0_rt; + bool mixed = false; + for (int kk = amrex::max(k - 1, klo); kk <= amrex::min(k + 1, khi) && !mixed; ++kk) { + for (int jj = amrex::max(j - 1, jlo); jj <= amrex::min(j + 1, jhi) && !mixed; + ++jj) { + for (int ii = amrex::max(i - 1, ilo); ii <= amrex::min(i + 1, ihi); ++ii) { + if ((phi(ii, jj, kk) > 0.0_rt) != fluid) { + mixed = true; + break; + } + } + } + } + band(i, j, k) = static_cast(mixed); + }); + } + + auto const plo = geom.ProbLoArray(); + auto const dx = geom.CellSizeArray(); + auto const ixt = mf.ixType(); + RealVect const offset(AMREX_D_DECL(ixt.cellCentered(0) ? 0.5_rt : 0.0_rt, + ixt.cellCentered(1) ? 0.5_rt : 0.0_rt, + ixt.cellCentered(2) ? 0.5_rt : 0.0_rt)); + + auto const* bvh_root = m_bvh_nodes.data(); + auto const* tri_pts = m_tri_pts_d.data(); + int const num_triangles = m_num_tri; + bool const use_bvh = m_bvh_optimization; + + for (MFIter mfi(mf); mfi.isValid(); ++mfi) { + Box const bx = mf[mfi].box(); + auto const phi = mf.array(mfi); + auto const band = exact_band.const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (band(i, j, k) == 0) { + return; + } + XDim3 const coords{.x = plo[0] + (static_cast(i) + offset[0]) * dx[0], + .y = plo[1] + (static_cast(j) + offset[1]) * dx[1], + .z = plo[2] + (static_cast(k) + offset[2]) * dx[2]}; + Real d2 = std::numeric_limits::max(); + if (use_bvh) { + d2 = bvh_d2(coords, bvh_root); + } else { + for (int tr = 0; tr < num_triangles; ++tr) { + d2 = amrex::min(d2, pt_tri_min_d2(coords, tri_pts[tr])); + } + } + Real const distance = std::sqrt(d2); + Real const coordinate_scale = amrex::max( + 1.0_rt, amrex::max(std::abs(coords.x), + amrex::max(std::abs(coords.y), std::abs(coords.z)))); + Real const on_surface_tolerance = + 16.0_rt * std::numeric_limits::epsilon() * coordinate_scale; + phi(i, j, k) = distance <= on_surface_tolerance + ? 0.0_rt + : phi(i, j, k) * distance; + }); + } + Gpu::streamSynchronize(); + mf.OverrideSync(geom.periodicity()); + mf.FillBoundary(geom.periodicity()); +#endif +} + +void STLtools::fillMarchingCubesEdgeIntersections ( + FArrayBox const& levelset, Geometry const& geom, MC::MCFab& mc_fab) const +{ + BL_PROFILE("STLtools::fillMarchingCubesEdgeIntersections"); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + AMREX_SPACEDIM == 3, + "STL marching-cubes edge intersections are only available in 3D"); + +#if (AMREX_SPACEDIM != 3) + amrex::ignore_unused(this, levelset, geom, mc_fab); +#else + auto const sdf = levelset.const_array(); + auto const plo = geom.ProbLoArray(); + auto const dx = geom.CellSizeArray(); + auto const* tri_pts = m_tri_pts_d.data(); + auto const* tri_norm = m_tri_normals_d.data(); + auto const* bvh_root = m_bvh_nodes.data(); + int const num_triangles = m_num_tri; + bool const use_bvh = m_bvh_optimization; + + auto find_x = [=] AMREX_GPU_DEVICE (int i, int j, int k, + Real& fraction) noexcept { + XDim3 const p1{.x = plo[0]+static_cast(i)*dx[0], + .y = plo[1]+static_cast(j)*dx[1], + .z = plo[2]+static_cast(k)*dx[2]}; + Real const x2 = p1.x + dx[0]; + Real const target = sdf(i,j,k)/(sdf(i,j,k)-sdf(i+1,j,k)); + Real best_distance = std::numeric_limits::max(); + bool found = false; + auto test = [&] AMREX_GPU_DEVICE (int ntri, Triangle const* tris, + XDim3 const* norms) noexcept -> int { + for (int it = 0; it < ntri; ++it) { + auto const hit = edge_tri_intersects( + p1.x, x2, p1.y, p1.z, tris[it].v1, tris[it].v2, + tris[it].v3, norms[it], sdf(i,j,k)-sdf(i+1,j,k)); + if (hit.first) { + Real const candidate = amrex::Clamp( + (hit.second-p1.x)/dx[0], 0.0_rt, 1.0_rt); + Real const distance = std::abs(candidate-target); + if (distance < best_distance) { + fraction = candidate; + best_distance = distance; + found = true; + } + } + } + return 0; + }; + if (use_bvh) { + Real a[3] = {p1.x,p1.y,p1.z}; + Real b[3] = {x2,p1.y,p1.z}; + bvh_line_tri_intersects(a,b,bvh_root,test); + } else { + test(num_triangles,tri_pts,tri_norm); + } + return found; + }; + + auto find_y = [=] AMREX_GPU_DEVICE (int i, int j, int k, + Real& fraction) noexcept { + XDim3 const p1{.x = plo[0]+static_cast(i)*dx[0], + .y = plo[1]+static_cast(j)*dx[1], + .z = plo[2]+static_cast(k)*dx[2]}; + Real const y2 = p1.y + dx[1]; + Real const target = sdf(i,j,k)/(sdf(i,j,k)-sdf(i,j+1,k)); + Real best_distance = std::numeric_limits::max(); + bool found = false; + auto test = [&] AMREX_GPU_DEVICE (int ntri, Triangle const* tris, + XDim3 const* norms) noexcept -> int { + for (int it = 0; it < ntri; ++it) { + auto const& tri = tris[it]; + auto const& norm = norms[it]; + auto const hit = edge_tri_intersects( + p1.y, y2, p1.z, p1.x, + XDim3{.x = tri.v1.y, .y = tri.v1.z, .z = tri.v1.x}, + XDim3{.x = tri.v2.y, .y = tri.v2.z, .z = tri.v2.x}, + XDim3{.x = tri.v3.y, .y = tri.v3.z, .z = tri.v3.x}, + XDim3{.x = norm.y, .y = norm.z, .z = norm.x}, + sdf(i,j,k)-sdf(i,j+1,k)); + if (hit.first) { + Real const candidate = amrex::Clamp( + (hit.second-p1.y)/dx[1], 0.0_rt, 1.0_rt); + Real const distance = std::abs(candidate-target); + if (distance < best_distance) { + fraction = candidate; + best_distance = distance; + found = true; + } + } + } + return 0; + }; + if (use_bvh) { + Real a[3] = {p1.x,p1.y,p1.z}; + Real b[3] = {p1.x,y2,p1.z}; + bvh_line_tri_intersects(a,b,bvh_root,test); + } else { + test(num_triangles,tri_pts,tri_norm); + } + return found; + }; + + auto find_z = [=] AMREX_GPU_DEVICE (int i, int j, int k, + Real& fraction) noexcept { + XDim3 const p1{.x = plo[0]+static_cast(i)*dx[0], + .y = plo[1]+static_cast(j)*dx[1], + .z = plo[2]+static_cast(k)*dx[2]}; + Real const z2 = p1.z + dx[2]; + Real const target = sdf(i,j,k)/(sdf(i,j,k)-sdf(i,j,k+1)); + Real best_distance = std::numeric_limits::max(); + bool found = false; + auto test = [&] AMREX_GPU_DEVICE (int ntri, Triangle const* tris, + XDim3 const* norms) noexcept -> int { + for (int it = 0; it < ntri; ++it) { + auto const& tri = tris[it]; + auto const& norm = norms[it]; + auto const hit = edge_tri_intersects( + p1.z, z2, p1.x, p1.y, + XDim3{.x = tri.v1.z, .y = tri.v1.x, .z = tri.v1.y}, + XDim3{.x = tri.v2.z, .y = tri.v2.x, .z = tri.v2.y}, + XDim3{.x = tri.v3.z, .y = tri.v3.x, .z = tri.v3.y}, + XDim3{.x = norm.z, .y = norm.x, .z = norm.y}, + sdf(i,j,k)-sdf(i,j,k+1)); + if (hit.first) { + Real const candidate = amrex::Clamp( + (hit.second-p1.z)/dx[2], 0.0_rt, 1.0_rt); + Real const distance = std::abs(candidate-target); + if (distance < best_distance) { + fraction = candidate; + best_distance = distance; + found = true; + } + } + } + return 0; + }; + if (use_bvh) { + Real a[3] = {p1.x,p1.y,p1.z}; + Real b[3] = {p1.x,p1.y,z2}; + bvh_line_tri_intersects(a,b,bvh_root,test); + } else { + test(num_triangles,tri_pts,tri_norm); + } + return found; + }; + + auto const exact_x = mc_fab.m_edge_intersections[0].array(); + auto const exact_y = mc_fab.m_edge_intersections[1].array(); + auto const exact_z = mc_fab.m_edge_intersections[2].array(); + ParallelFor(Box{exact_x}, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + if ((sdf(i,j,k) > 0.0_rt) != (sdf(i+1,j,k) > 0.0_rt)) { + Real fraction = MC::invalid_edge_intersection; + find_x(i,j,k,fraction); + exact_x(i,j,k) = fraction; + } + }); + ParallelFor(Box{exact_y}, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + if ((sdf(i,j,k) > 0.0_rt) != (sdf(i,j+1,k) > 0.0_rt)) { + Real fraction = MC::invalid_edge_intersection; + find_y(i,j,k,fraction); + exact_y(i,j,k) = fraction; + } + }); + ParallelFor(Box{exact_z}, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept { + if ((sdf(i,j,k) > 0.0_rt) != (sdf(i,j,k+1) > 0.0_rt)) { + Real fraction = MC::invalid_edge_intersection; + find_z(i,j,k,fraction); + exact_z(i,j,k) = fraction; + } + }); +#endif } + +} // namespace amrex diff --git a/Src/EB/AMReX_MarchingCubes.H b/Src/EB/AMReX_MarchingCubes.H index 6830eb51002..a51b09d1db4 100644 --- a/Src/EB/AMReX_MarchingCubes.H +++ b/Src/EB/AMReX_MarchingCubes.H @@ -2,36 +2,113 @@ #define AMREX_MARCHING_CUBES_H_ #include +#include +#include #include #include +#include +#include -#include -#include #include /** * \file AMReX_MarchingCubes.H * * GPU-enabled marching-cubes routines used to extract EB surfaces. + * + * Geometry-source contract + * ------------------------ + * The marching-cubes EB builder is independent of how the geometry is + * described (STL triangulation, analytic implicit function, parser + * expression, ...). A geometry source feeds it through exactly two inputs: + * + * 1. A nodal scalar field (the "MC level set"). It uses the marching-cubes + * sign convention: values greater than zero are fluid, values less than + * zero are covered, and an exact zero marks a node lying on the boundary + * (it takes the covered sign for the two-state MC33 lookup, and edge + * crossings are placed exactly at that node). Only the sign and local + * continuity are required; magnitudes are used for the linear crossing + * fallback and the MC33 asymptotic decider. Note that this is the + * opposite of the EB2 implicit-function convention (negative in fluid), + * which is why the builder negates the field on ingest and again when it + * exports the public level set. + * + * 2. Optional exact crossings on Cartesian edges, stored in + * MCFab::m_edge_intersections as the normalized position in [0,1] measured + * from the low node. Edges without an exact crossing keep + * invalid_edge_intersection and fall back to linear interpolation of the + * nodal field. + * + * See AMReX_EB2_Level_MC.H for the adapters that produce these inputs from + * EB2::GeometryShop implicit functions and from STLtools. */ namespace amrex::MC { +//! Sentinel used when an original-STL edge crossing is unavailable. +inline constexpr Real invalid_edge_intersection = -1.0_rt; + +enum RejectionReason : int { + invalid_topology = 1, + small_volume = 2 +}; + +/** + * Components stored in MCFab::m_cell_data. Besides compact triangle and + * interior-vertex ranges, these retain the MC33 face decisions consumed by + * downstream face construction. + */ +enum CellDataComponent : int { + triangle_count = 0, + triangle_offset, + interior_vertex_count, + interior_vertex_offset, + face_decision_valid_mask, + face_fluid_connected_mask, + num_cell_data_components +}; + +/** + * \brief Slots of the per-FAB counter block that the entry points below + * accumulate into. + * + * The caller allocates num_fab_counters ints per FAB in device-accessible + * memory (Gpu::Buffer or Gpu::DeviceVector), zeroes them once per repair + * pass, passes each FAB's block to the entry points, and copies all blocks to + * the host once per pass. The entry points themselves never synchronize with + * the host to report their counts. + */ +enum FabCounter : int { + counter_invalid_triangles = 0, //!< marching_cubes: triangle without a crossing vertex + counter_face_decision_errors, //!< build_face_fractions + counter_degenerate_faces, //!< build_face_fractions + counter_closure_errors, //!< build_cell_fractions + counter_volume_errors, //!< build_cell_fractions + counter_centroid_errors, //!< build_cell_fractions + counter_area_vector_errors, //!< build_cell_fractions + counter_face_rejections, //!< mark_faces_for_cleanup + counter_topology_rejections, //!< mark_cells_for_cleanup + counter_small_cell_rejections, //!< mark_cells_for_cleanup + counter_changed_nodes, //!< zero_nodes_for_cleanup + counter_extended_nodes, //!< extend_domain_face_levelset + num_fab_counters +}; + //! Initialize internal lookup tables and device buffers for marching cubes. void Initialize (); //! Destroy any static storage allocated by Initialize(). void Finalize (); -//! Vertex buffer storing coordinates and normals for extracted surfaces. +//! Vertex buffer storing coordinates for extracted surfaces. struct Vertex { - Gpu::DeviceVector x, y, z, nx, ny, nz; + Gpu::DeviceVector x, y, z; //! Resize the vertex buffers to hold \p n entries. void resize (int n); //! Raw device pointers used when launching kernels. - GpuArray dataPtrs () { return {x.data(), y.data(), z.data(), nx.data(), ny.data(), nz.data()}; } + GpuArray dataPtrs () { return {x.data(), y.data(), z.data()}; } }; //! Triangle index buffer referencing Vertex entries. @@ -47,19 +124,175 @@ struct Triangle }; //! Minimal container holding marching-cubes output for one FAB. -class MCFab +struct MCFab { -public: - // We will store more data here in the future + //! Per-cell triangle/interior ranges and resolved MC33 face decisions. + BaseFab m_cell_data; + //! Exact normalized STL crossings on Cartesian edges, when available. + Array,AMREX_SPACEDIM> m_edge_intersections; Triangle m_triangles; Vertex m_vertices; + + //! Allocate edge-intersection caches for the usable nodal SDF region. + void defineEdgeIntersections (Box const& sdf_box); }; -//! Run marching cubes on signed-distance field \p sdf_fab. -void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab); +/** + * \brief Run marching cubes on signed-distance field \p sdf_fab. + * + * \p counters is the FAB's counter block (see FabCounter); a nonzero + * counter_invalid_triangles afterwards indicates a lookup-table invariant + * violation and the caller must abort. + */ +void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab, int* counters); + +/** + * \brief Build Cartesian-face aperture fractions and centroids from the MC + * signed-distance field. + * + * \p bx is cell-centered; all faces surrounding it are updated. Faces whose + * fluid polygon is degenerate (for example a four-crossing face whose + * crossings all sit on repaired zero nodes) receive zero aperture and are + * marked with 1 in \p rejected_x / \p rejected_y / \p rejected_z where those + * arrays contain the face, so the nodal repair pass covers them. + * + * Adds the number of faces whose two cells resolved the MC33 ambiguity + * differently (an invariant violation) to counter_face_decision_errors and + * the number of degenerate faces marked for repair to counter_degenerate_faces + * of \p counters. + */ +void build_face_fractions ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf, + FArrayBox& apx, FArrayBox& apy, FArrayBox& apz, + FArrayBox& fcx, FArrayBox& fcy, FArrayBox& fcz, + IArrayBox& rejected_x, IArrayBox& rejected_y, IArrayBox& rejected_z, + int* counters); + +/** + * \brief Build EB edge-centroid data from the MC signed-distance field. + * + * Regular and covered edges receive the standard EB2 sentinels 1 and -1. + * Cut-edge values are the centroid of the fluid part of the unit edge. + */ +void build_edge_centroids ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf, + FArrayBox& ecx, FArrayBox& ecy, FArrayBox& ecz); + +/** + * \brief Replace single-valued cut-cell geometry with moments computed from + * the marching-cubes triangles. + * + * Coordinate-face apertures previously built from the MC signed-distance + * field and the cell-owned embedded-boundary triangles form the closed + * boundary used to compute volume and first moments. Requires cubic cells + * (dx == dy == dz); the EB2 driver checks this once per level. + * + * Single-valued cells whose triangulation produced invalid geometry are + * counted in counter_closure_errors, counter_volume_errors, + * counter_centroid_errors and counter_area_vector_errors of \p counters. + */ +void build_cell_fractions ( + Box const& bx, Geometry const& geom, MCFab const& mc_fab, + FArrayBox const& sdf, + FArrayBox& apx, FArrayBox& apy, FArrayBox& apz, + FArrayBox& vfrac, FArrayBox& vcent, FArrayBox& barea, + FArrayBox& bcent, FArrayBox& bnorm, int* counters); + +/** + * \brief Classify cells and build all EBCellFlag connectivity from MC SDF and + * face apertures. + * + * Cell types are derived from the converged volume and face apertures, + * including coincident nodal surfaces. MC ownership is checked against the + * repaired SDF before coordinate connectivity is populated. + * + * \return Number of inconsistent cell classifications. + */ +int build_cell_topology ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf, + EBCellFlagFab& cellflag, FArrayBox const& vfrac, + FArrayBox const& apx, FArrayBox const& apy, FArrayBox const& apz); + +/** + * \brief Mark Cartesian faces with more than one resolved fluid aperture. + * + * Marks are OR-ed into the arrays, so faces already rejected by + * build_face_fractions are kept. The number of newly rejected faces across + * all three coordinate directions is added to counter_face_rejections of + * \p counters. + */ +void mark_faces_for_cleanup ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf, + IArrayBox& rejected_x, IArrayBox& rejected_y, + IArrayBox& rejected_z, int* counters); + +/** + * \brief Mark cells that the single-VoF MC builder must remove. + * + * A cell is rejected if its MC boundary is invalid, if its fluid corners form + * more than one group, or if its volume fraction is below \p small_volfrac. + * Corner groups are joined along cube edges and across ambiguous faces whose + * MC33 decision is "fluid connected" (the same decision the surface and the + * face apertures use), but never through the cell interior: MC33 tunnel + * tilings (4.1.2, 10.1.2, ...) and the split tilings 7.3/10.2/12.2/13.x are + * always rejected, so an exported single-valued cell holds exactly one + * face-connected fluid region. Face-aperture rejection is handled separately + * by mark_faces_for_cleanup and is consistent with this rule (a face whose + * decision is "not connected" is a rejected face and also splits the corner + * groups of both cells). + * + * Topology and small-cell rejections are added to counter_topology_rejections + * and counter_small_cell_rejections of \p counters. \p rejected records the + * reason of each rejected cell as one of the RejectionReason values (topology + * takes precedence) and 0 otherwise. + */ +void mark_cells_for_cleanup ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf, + FArrayBox const& vfrac, Real small_volfrac, IArrayBox& rejected, + int* counters); + +/** + * \brief Apply the legacy EB repair transform to the MC working level set. + * + * Positive (fluid) nodes of \p node_box incident to a rejected face or cell + * are moved to exact zero. The operation is monotone and leaves all other + * samples intact. The marks must cover the cells and faces enclosed by + * \p node_box; rejected cells and faces beyond the arrays are ignored, so a + * caller can repair its ghost nodes from the marks it holds. The number of + * nodes changed from positive to zero is added to counter_changed_nodes of + * \p counters. + */ +void zero_nodes_for_cleanup ( + Box const& node_box, IArrayBox const& rejected_cells, + IArrayBox const& rejected_x, IArrayBox const& rejected_y, + IArrayBox const& rejected_z, FArrayBox& sdf, int* counters); + +/** + * \brief Extrude the nodal level set straight outward from physical domain faces. + * + * An exterior node is assigned the value of the corresponding physical-face + * node only when it is outside in exactly one non-periodic direction. This + * preserves the feature's transverse footprint and does not extend diagonally + * from domain edges or corners. Exterior nodes that the nodal repair has + * moved to exact zero are left at zero, so the repair stays monotone. The + * number of exterior nodal values changed is added to counter_extended_nodes + * of \p counters. + */ +void extend_domain_face_levelset ( + Box const& node_box, Box const& domain, GpuArray const& is_periodic, + FArrayBox& sdf, int* counters); + +/** + * \brief Extrude exact transverse-edge crossings with the domain-face level set. + * + * Crossings on edges transverse to the extension normal are copied straight + * outward from the physical face. Edges parallel to the normal are unchanged. + */ +void extend_domain_face_edge_intersections ( + Box const& domain, GpuArray const& is_periodic, MCFab& mc_fab); -//! Write the collected marching-cubes output to an STL file. -void write_stl (std::string const& filename, std::map> const& mc_fabs); +//! Write the converged MC triangles owned by the valid boxes of \p mc_fabs to an STL file. +void write_stl (std::string const& filename, LayoutData const& mc_fabs); } diff --git a/Src/EB/AMReX_MarchingCubes.cpp b/Src/EB/AMReX_MarchingCubes.cpp index a76e037cc98..497bd010656 100644 --- a/Src/EB/AMReX_MarchingCubes.cpp +++ b/Src/EB/AMReX_MarchingCubes.cpp @@ -2,12 +2,17 @@ #include #include #include +#include #include +#include #include +#include #include +#include #include +#include /* * http://thomas.lewiner.org/publication_page.php%EF%B9%96pubkey=marching_cubes_jgt.html @@ -27,11 +32,285 @@ namespace { LookUpTable* h_table = nullptr; LookUpTable* d_table = nullptr; +/** + * Sign of a*c - b*d with a scale-free tie break: the products have the units + * of the level set squared, so a tie is declared relative to their magnitude + * rather than against an absolute epsilon. Returns -1, 0 (tie) or 1. Shared + * by the MC33 face decider and the interior (tunnel) test so that both are + * invariant under scaling of the implicit function. + */ +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +int ambiguous_product_sign (Real a, Real c, Real b, Real d) noexcept +{ + Real const p = a*c; + Real const q = b*d; + Real const r = p - q; + Real const tol = Real(8.0)*std::numeric_limits::epsilon()*(std::abs(p)+std::abs(q)); + if (std::abs(r) <= tol) { return 0; } + return (r > 0.0_rt) ? 1 : -1; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +bool four_crossing_fluid_is_connected (Real const* levelset) noexcept +{ + // The tie break is scale free (see ambiguous_product_sign); |a|+|b| is + // invariant under the index rotations/reversals neighboring cells apply + // to a shared face. + int const sign = ambiguous_product_sign(levelset[0], levelset[2], levelset[1], levelset[3]); + if (sign == 0) { + // MC33's test_face tie break always selects the positive material. + // Expressing the decision in material terms makes it invariant under + // the rotations/reversals used by neighboring cells on a shared face. + return true; + } + return (levelset[0] > 0.0_rt) == (sign > 0); +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +EB2::Type_t face_type (Real area, Real tolerance) noexcept +{ + if (area <= tolerance) { + return EB2::Type::covered; + } else if (area >= 1.0_rt-tolerance) { + return EB2::Type::regular; + } else { + return EB2::Type::irregular; + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +Real edge_intersection_fraction (Real lo, Real hi, Real exact) noexcept +{ + // Cleanup moves fluid nodes to exactly zero. Those repaired crossings + // belong at the node, rather than at the original STL intersection. + if (lo == 0.0_rt) { + return 0.0_rt; + } + if (hi == 0.0_rt) { + return 1.0_rt; + } + return (exact < 0.0_rt || exact > 1.0_rt) ? lo/(lo-hi) : exact; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void consume_face_decision (Array4 const& cell_data, + Box const& cell_box, int i, int j, int k, + int face, int& valid, int& connected) noexcept +{ + if (valid == 0 && cell_box.contains(i,j,k)) { + int const bit = 1 << face; + if ((cell_data(i,j,k,face_decision_valid_mask) & bit) != 0) { + valid = 1; + connected = + (cell_data(i,j,k,face_fluid_connected_mask) & bit) != 0; + } + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +void accumulate_face_decision (Array4 const& cell_data, + Box const& cell_box, int i, int j, int k, + int face, int* error, + int& valid, int& connected) noexcept +{ + if (cell_box.contains(i,j,k)) { + int const bit = 1 << face; + int const valid_mask = + cell_data(i,j,k,face_decision_valid_mask); + if ((valid_mask & bit) != 0) { + int const decision = + (cell_data(i,j,k,face_fluid_connected_mask) & bit) != 0; + if (valid != 0 && connected != decision) { + Gpu::Atomic::AddNoRet(error,1); + } else if (valid == 0) { + connected = decision; + } + valid = 1; + } + } +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +GpuArray resolved_face_decision ( + Array4 const& cell_data, Box const& cell_box, int* error, + int ilo, int jlo, int klo, int flo, + int ihi, int jhi, int khi, int fhi) noexcept +{ + int valid = 0; + int connected = 0; + accumulate_face_decision(cell_data, cell_box, ilo, jlo, klo, flo, + error, valid, connected); + accumulate_face_decision(cell_data, cell_box, ihi, jhi, khi, fhi, + error, valid, connected); + return {valid,connected}; +} + +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +bool face_is_rejected (Real a, Real b, Real c, Real d, + Array4 const& cell_data, + Box const& cell_box, + int ilo, int jlo, int klo, int flo, + int ihi, int jhi, int khi, int fhi) noexcept +{ + bool const f0 = a > 0.0_rt; + bool const f1 = b > 0.0_rt; + bool const f2 = c > 0.0_rt; + bool const f3 = d > 0.0_rt; + if (!(f0 == f2 && f1 == f3 && f0 != f1)) { + return false; + } + + int valid = 0; + int connected = 0; + consume_face_decision(cell_data, cell_box, ilo, jlo, klo, flo, + valid, connected); + consume_face_decision(cell_data, cell_box, ihi, jhi, khi, fhi, + valid, connected); + + Real const values[4] = {a, b, c, d}; + bool const fluid_connected = valid != 0 + ? connected != 0 + : four_crossing_fluid_is_connected(values); + return !fluid_connected; +} + +//! Returns false when the face polygon is degenerate; the caller then routes +//! the face into the nodal repair set. +AMREX_GPU_DEVICE AMREX_FORCE_INLINE +bool cut_face_fraction (Real const* levelset, Real const* intersections, + Real& area, + Real& centroid_x, Real& centroid_y, + bool has_resolved_decision, + bool resolved_fluid_connected) noexcept +{ + // Safe sentinels for every early-error path. The caller will reject the + // geometry, but no downstream kernel may observe uninitialized storage. + area = 0.0_rt; + centroid_x = 0.0_rt; + centroid_y = 0.0_rt; + + constexpr Real vertex_x[4] = {0.0_rt, 1.0_rt, 1.0_rt, 0.0_rt}; + constexpr Real vertex_y[4] = {0.0_rt, 0.0_rt, 1.0_rt, 1.0_rt}; + + Real polygon_x[8]; + Real polygon_y[8]; + int polygon_size = 0; + int crossing_count = 0; + + for (int n = 0; n < 4; ++n) { + int const next = (n+1) % 4; + bool const fluid = levelset[n] > 0.0_rt; + bool const next_fluid = levelset[next] > 0.0_rt; + + if (fluid) { + polygon_x[polygon_size] = vertex_x[n]; + polygon_y[polygon_size] = vertex_y[n]; + ++polygon_size; + } + if (fluid != next_fluid) { + Real const alpha = edge_intersection_fraction( + levelset[n], levelset[next], intersections[n]); + polygon_x[polygon_size] = + vertex_x[n] + alpha*(vertex_x[next]-vertex_x[n]); + polygon_y[polygon_size] = + vertex_y[n] + alpha*(vertex_y[next]-vertex_y[n]); + ++polygon_size; + ++crossing_count; + } + } + + if (crossing_count == 0) { + area = (levelset[0] > 0.0_rt) ? 1.0_rt : 0.0_rt; + centroid_x = 0.0_rt; + centroid_y = 0.0_rt; + return true; + } + + if (crossing_count == 4) { + Real edge_x[4]; + Real edge_y[4]; + for (int n = 0; n < 4; ++n) { + int const next = (n+1) % 4; + Real const alpha = edge_intersection_fraction( + levelset[n], levelset[next], intersections[n]); + edge_x[n] = vertex_x[n] + + alpha*(vertex_x[next]-vertex_x[n]); + edge_y[n] = vertex_y[n] + + alpha*(vertex_y[next]-vertex_y[n]); + } + + bool const fluid_connected = has_resolved_decision + ? resolved_fluid_connected + : four_crossing_fluid_is_connected(levelset); + Real fluid_area = fluid_connected ? 1.0_rt : 0.0_rt; + Real fluid_moment_x = fluid_connected ? 0.5_rt : 0.0_rt; + Real fluid_moment_y = fluid_connected ? 0.5_rt : 0.0_rt; + + // The asymptotic decider selects which diagonal pair is connected. + // With straight MC face segments, the other sign consists of the two + // corner triangles. Add the fluid triangles or subtract the covered + // triangles from the unit square. + for (int n = 0; n < 4; ++n) { + bool const fluid = levelset[n] > 0.0_rt; + if (fluid != fluid_connected) { + int const previous = (n+3) % 4; + Real const ax = edge_x[n] - vertex_x[n]; + Real const ay = edge_y[n] - vertex_y[n]; + Real const bx = edge_x[previous] - vertex_x[n]; + Real const by = edge_y[previous] - vertex_y[n]; + Real const triangle_area = + 0.5_rt*std::abs(ax*by-ay*bx); + Real const triangle_centroid_x = + (vertex_x[n]+edge_x[n]+edge_x[previous])/3.0_rt; + Real const triangle_centroid_y = + (vertex_y[n]+edge_y[n]+edge_y[previous])/3.0_rt; + Real const sign = fluid_connected ? -1.0_rt : 1.0_rt; + fluid_area += sign*triangle_area; + fluid_moment_x += sign*triangle_area*triangle_centroid_x; + fluid_moment_y += sign*triangle_area*triangle_centroid_y; + } + } + + if (fluid_area <= 0.0_rt || fluid_area >= 1.0_rt) { + return false; + } + area = fluid_area; + centroid_x = fluid_moment_x/fluid_area - 0.5_rt; + centroid_y = fluid_moment_y/fluid_area - 0.5_rt; + return true; + } + + if (crossing_count != 2 || polygon_size < 3) { + return false; + } + + Real twice_area = 0.0_rt; + Real centroid_x_numerator = 0.0_rt; + Real centroid_y_numerator = 0.0_rt; + for (int n = 0; n < polygon_size; ++n) { + int const next = (n+1) % polygon_size; + Real const cross = polygon_x[n]*polygon_y[next] + - polygon_x[next]*polygon_y[n]; + twice_area += cross; + centroid_x_numerator += (polygon_x[n]+polygon_x[next])*cross; + centroid_y_numerator += (polygon_y[n]+polygon_y[next])*cross; + } + + if (twice_area <= 0.0_rt) { + return false; + } + + area = 0.5_rt*twice_area; + centroid_x = centroid_x_numerator/(3.0_rt*twice_area) - 0.5_rt; + centroid_y = centroid_y_numerator/(3.0_rt*twice_area) - 0.5_rt; + return true; +} + AMREX_GPU_DEVICE AMREX_FORCE_INLINE void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int k, Array4 const& sdf, Array4 const& ex, Array4 const& ey, Array4 const& ez, - GpuArray const& pvrtx, Array4 const& ntri, + GpuArray const& pvrtx, Array4 const& ntri, GpuArray const& ptri, int* error) { std::uint8_t lut_entry = 0; @@ -46,6 +325,9 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int if (cube[6] > 0) { lut_entry += 64; } if (cube[7] > 0) { lut_entry += 128; } + int face_valid_mask = 0; + int face_connected_mask = 0; + auto add_triangle = [&] (const std::int8_t* trig, std::int8_t n, int v12 = -1) -> int { if (ipass == 0) { return n; } @@ -77,7 +359,7 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int if( t%3 == 2 ) { - auto m = ntri(i,j,k,1) + r; + auto m = ntri(i,j,k,triangle_offset) + r; ptri[0][m] = tv[0]; ptri[1][m] = tv[1]; ptri[2][m] = tv[2]; @@ -103,14 +385,18 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int default : *error = 1 ; A = B = C = D = 0 ; }; - if( std::abs( A*C - B*D ) < std::numeric_limits::epsilon() ) { - return face >= 0 ; - } else { - return Real(face) * A * ( A*C - B*D ) >= 0 ; // face and A invert signs + Real const values[4] = {A,B,C,D}; + bool const fluid_connected = four_crossing_fluid_is_connected(values); + int const bit = 1 << (std::abs(int(face))-1); + face_valid_mask |= bit; + if (fluid_connected) { + face_connected_mask |= bit; } + return fluid_connected == (face > 0); }; - auto test_interior = [&] (std::int8_t _case, std::int8_t _config, std::int8_t _subconfig, std::int8_t s) + auto test_interior_impl = [&] (std::int8_t _case, std::int8_t _config, + std::int8_t _subconfig, std::int8_t s) { Real t, At=0, Bt=0, Ct=0, Dt=0, a, b ; std::int8_t test = 0 ; @@ -136,6 +422,14 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int case 7 : case 12 : case 13 : + // Lewiner's reference-edge shortcut: the interior test is taken on + // the plane through the reference edge with At = 0. With that + // choice the tunnel alternatives 6.1.2, 7.4.2, 12.1.2 and 13.5.2 + // are never selected (Bt >= 0 is equivalent to the face test that + // has already failed), so these cases always take the split + // tiling. That is conservative for the EB builder: split cells + // hold two fluid corner groups and are rejected and repaired by + // mark_cells_for_cleanup, exactly as tunnel cells would be. switch( _case ) // NOLINT(bugprone-switch-missing-default-case) { case 6 : edge = lut->test6 [_config][2] ; break ; @@ -247,12 +541,12 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int case 2 : return s>0 ; case 3 : return s>0 ; case 4 : return s>0 ; - case 5 : if( At * Ct - Bt * Dt < FLT_EPSILON ) { return s>0 ; } break; + case 5 : if( ambiguous_product_sign(At,Ct,Bt,Dt) <= 0 ) { return s>0 ; } break; case 6 : return s>0 ; case 7 : return s<0 ; case 8 : return s>0 ; case 9 : return s>0 ; - case 10 : if( At * Ct - Bt * Dt >= FLT_EPSILON ) { return s>0 ; } break; + case 10 : if( ambiguous_product_sign(At,Ct,Bt,Dt) >= 0 ) { return s>0 ; } break; case 11 : return s<0 ; case 12 : return s>0 ; case 13 : return s<0 ; @@ -263,31 +557,34 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int return s<0 ; }; + auto test_interior = [&] (std::int8_t c, std::int8_t config, + std::int8_t subconfig, std::int8_t s) -> bool + { + return test_interior_impl(c,config,subconfig,s); + }; + auto add_c_vertex = [&] () -> int { - if (ipass == 0) { ntri(i,j,k,2) = 1; return -1; } + if (ipass == 0) { + ntri(i,j,k,interior_vertex_count) = 1; + return -1; + } Real u = 0 ; int vid ; - auto m = ntri(i,j,k,3); + auto m = ntri(i,j,k,interior_vertex_offset); auto& vert_x = pvrtx[0][m]; auto& vert_y = pvrtx[1][m]; auto& vert_z = pvrtx[2][m]; - auto& vert_nx = pvrtx[3][m]; - auto& vert_ny = pvrtx[4][m]; - auto& vert_nz = pvrtx[5][m]; - vert_x = vert_y = vert_z = vert_nx = vert_ny = vert_nz = 0 ; + vert_x = vert_y = vert_z = 0 ; auto update_vertex = [&] () { ++u ; vert_x += pvrtx[0][vid]; vert_y += pvrtx[1][vid]; vert_z += pvrtx[2][vid]; - vert_nx += pvrtx[3][vid]; - vert_ny += pvrtx[4][vid]; - vert_nz += pvrtx[5][vid]; }; @@ -321,14 +618,6 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int vert_y *= Real(1)/u ; vert_z *= Real(1)/u ; - u = std::sqrt( vert_nx * vert_nx + vert_ny * vert_ny +vert_nz * vert_nz ) ; - if( u > 0 ) - { - vert_nx *= Real(1)/u ; - vert_ny *= Real(1)/u ; - vert_nz *= Real(1)/u ; - } - return m; }; @@ -642,7 +931,9 @@ void process_cube (std::int8_t ipass, LookUpTable const* lut, int i, int j, int }; if (ipass == 0) { - ntri(i,j,k,0) = nt; + ntri(i,j,k,triangle_count) = nt; + ntri(i,j,k,face_decision_valid_mask) = face_valid_mask; + ntri(i,j,k,face_fluid_connected_mask) = face_connected_mask; } } @@ -667,6 +958,7 @@ void Finalize () delete h_table; h_table = nullptr; #ifdef AMREX_USE_GPU + Gpu::streamSynchronizeAll(); The_Arena()->free(d_table); #endif d_table = nullptr; @@ -678,9 +970,6 @@ void Vertex::resize (int n) x.resize(n); y.resize(n); z.resize(n); - nx.resize(n); - ny.resize(n); - nz.resize(n); } void Triangle::resize (int n) @@ -690,23 +979,28 @@ void Triangle::resize (int n) v3.resize(n); } -void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) +void MCFab::defineEdgeIntersections (Box const& sdf_box) +{ + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + Box const edge_box = amrex::enclosedCells(sdf_box, idim); + m_edge_intersections[idim].resize(edge_box, 1); + m_edge_intersections[idim].setVal(invalid_edge_intersection); + } +} + +void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab, int* counters) { BL_PROFILE("marching_cubes"); - AMREX_ALWAYS_ASSERT(sdf_fab.numPts() < Long(std::numeric_limits::max())); + // The prefix sums below index vertices (up to 3 per node) and triangles + // (up to 12 per cell) with int, so bound the node count accordingly. + AMREX_ALWAYS_ASSERT(sdf_fab.numPts() < Long(std::numeric_limits::max())/12); - // Remove small numbers. + // Exact zeros belong to the covered side. This lets the cleanup loop move + // nodes to ON without introducing a new positive fluid sample. auto const& sdf = sdf_fab.array(); - ParallelFor(sdf_fab.box(), [=] AMREX_GPU_DEVICE (int i, int j, int k) - { - if (std::abs(sdf(i,j,k)) < std::numeric_limits::epsilon()) { - sdf(i,j,k) = std::numeric_limits::epsilon(); - } - }); - Box nbox = sdf_fab.box(); - nbox.grow(-1); // Shrink the box by 1 so that we can compute gradient of sdf + Box const nbox = sdf_fab.box(); Box cbox = amrex::enclosedCells(nbox); Box exbox = amrex::enclosedCells(nbox, 0); Box eybox = amrex::enclosedCells(nbox, 1); @@ -718,6 +1012,20 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) Array4 ex = ex_fab.array(); Array4 ey = ey_fab.array(); Array4 ez = ez_fab.array(); + if (mc_fab.m_edge_intersections[0].nComp() == 0 + || mc_fab.m_edge_intersections[1].nComp() == 0 + || mc_fab.m_edge_intersections[2].nComp() == 0 + || !mc_fab.m_edge_intersections[0].box().contains(exbox) + || !mc_fab.m_edge_intersections[1].box().contains(eybox) + || !mc_fab.m_edge_intersections[2].box().contains(ezbox)) + { + // Standalone MC callers do not have an original STL to refine + // against. Keep their established sampled-SDF interpolation. + mc_fab.defineEdgeIntersections(sdf_fab.box()); + } + auto const exact_x = mc_fab.m_edge_intersections[0].const_array(); + auto const exact_y = mc_fab.m_edge_intersections[1].const_array(); + auto const exact_z = mc_fab.m_edge_intersections[2].const_array(); BoxIndexer n_bi(nbox); auto nvx = Scan::PrefixSum(int(nbox.numPts()), @@ -725,24 +1033,24 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) auto [i,j,k] = n_bi(m); int vx = 0, vy = 0, vz = 0; if (ex.contains(i,j,k)) { - if (((sdf(i,j,k) < 0) && (sdf(i+1,j,k) > 0)) || - ((sdf(i,j,k) > 0) && (sdf(i+1,j,k) < 0))) + if ((sdf(i,j,k) > 0) + != (sdf(i+1,j,k) > 0)) { vx = 1; } ex(i,j,k,0) = vx; } if (ey.contains(i,j,k)) { - if (((sdf(i,j,k) < 0) && (sdf(i,j+1,k) > 0)) || - ((sdf(i,j,k) > 0) && (sdf(i,j+1,k) < 0))) + if ((sdf(i,j,k) > 0) + != (sdf(i,j+1,k) > 0)) { vy = 1; } ey(i,j,k,0) = vy; } if (ez.contains(i,j,k)) { - if (((sdf(i,j,k) < 0) && (sdf(i,j,k+1) > 0)) || - ((sdf(i,j,k) > 0) && (sdf(i,j,k+1) < 0))) + if ((sdf(i,j,k) > 0) + != (sdf(i,j,k+1) > 0)) { vz = 1; } @@ -752,16 +1060,21 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) }, [=] AMREX_GPU_DEVICE (int m, int ps) { auto [i,j,k] = n_bi(m); + // Component 1 holds the vertex index of the + // crossing on this edge, or -1 when the edge has + // no crossing. The -1 sentinel is what + // add_c_vertex and the triangle assembly test. if (ex.contains(i,j,k)) { - ex(i,j,k,1) = ps; - if (ex(i,j,k,0)) { ++ps; } + if (ex(i,j,k,0)) { ex(i,j,k,1) = ps++; } + else { ex(i,j,k,1) = -1; } } if (ey.contains(i,j,k)) { - ey(i,j,k,1) = ps; - if (ey(i,j,k,0)) { ++ps; } + if (ey(i,j,k,0)) { ey(i,j,k,1) = ps++; } + else { ey(i,j,k,1) = -1; } } if (ez.contains(i,j,k)) { - ez(i,j,k,1) = ps; + if (ez(i,j,k,0)) { ez(i,j,k,1) = ps; } + else { ez(i,j,k,1) = -1; } } }, Scan::Type::exclusive, Scan::retSum); @@ -774,69 +1087,27 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) { if (ex.contains(i,j,k) && ex(i,j,k,0)) { int m = ex(i,j,k,1); - Real u = sdf(i,j,k) / (sdf(i,j,k) - sdf(i+1,j,k)); + Real u = edge_intersection_fraction( + sdf(i,j,k), sdf(i+1,j,k), exact_x(i,j,k)); pvrtx[0][m] = Real(i) + u; pvrtx[1][m] = Real(j); pvrtx[2][m] = Real(k); - Real nx = (Real(1)-u) * (sdf(i+1,j ,k ) - sdf(i-1,j ,k )) - + u * (sdf(i+2,j, k ) - sdf(i ,j ,k )); - Real ny = (Real(1)-u) * (sdf(i ,j+1,k ) - sdf(i ,j-1,k )) - + u * (sdf(i+1,j+1,k ) - sdf(i+1,j-1,k )); - Real nz = (Real(1)-u) * (sdf(i ,j ,k+1) - sdf(i ,j ,k-1)) - + u * (sdf(i+1,j ,k+1) - sdf(i+1,j ,k-1)); - Real norm = std::sqrt(nx*nx + ny*ny + nz*nz); - if (norm > 0) { - nx *= Real(1)/norm; - ny *= Real(1)/norm; - nz *= Real(1)/norm; - } - pvrtx[3][m] = nx; - pvrtx[4][m] = ny; - pvrtx[5][m] = nz; } if (ey.contains(i,j,k) && ey(i,j,k,0)) { int m = ey(i,j,k,1); - Real u = sdf(i,j,k) / (sdf(i,j,k) - sdf(i,j+1,k)); + Real u = edge_intersection_fraction( + sdf(i,j,k), sdf(i,j+1,k), exact_y(i,j,k)); pvrtx[0][m] = Real(i); pvrtx[1][m] = Real(j) + u; pvrtx[2][m] = Real(k); - Real nx = (Real(1)-u) * (sdf(i+1,j ,k ) - sdf(i-1,j ,k )) - + u * (sdf(i+1,j+1,k ) - sdf(i-1,j+1,k )); - Real ny = (Real(1)-u) * (sdf(i ,j+1,k ) - sdf(i ,j-1,k )) - + u * (sdf(i ,j+2,k ) - sdf(i ,j ,k )); - Real nz = (Real(1)-u) * (sdf(i ,j ,k+1) - sdf(i ,j ,k-1)) - + u * (sdf(i ,j+1,k+1) - sdf(i ,j+1,k-1)); - Real norm = std::sqrt(nx*nx + ny*ny + nz*nz); - if (norm > 0) { - nx *= Real(1)/norm; - ny *= Real(1)/norm; - nz *= Real(1)/norm; - } - pvrtx[3][m] = nx; - pvrtx[4][m] = ny; - pvrtx[5][m] = nz; } if (ez.contains(i,j,k) && ez(i,j,k,0)) { int m = ez(i,j,k,1); - Real u = sdf(i,j,k) / (sdf(i,j,k) - sdf(i,j,k+1)); + Real u = edge_intersection_fraction( + sdf(i,j,k), sdf(i,j,k+1), exact_z(i,j,k)); pvrtx[0][m] = Real(i); pvrtx[1][m] = Real(j); pvrtx[2][m] = Real(k) + u; - Real nx = (Real(1)-u) * (sdf(i+1,j ,k ) - sdf(i-1,j ,k )) - + u * (sdf(i+1,j, k+1) - sdf(i-1,j ,k+1)); - Real ny = (Real(1)-u) * (sdf(i ,j+1,k ) - sdf(i ,j-1,k )) - + u * (sdf(i ,j+1,k+1) - sdf(i ,j-1,k+1)); - Real nz = (Real(1)-u) * (sdf(i ,j ,k+1) - sdf(i ,j ,k-1)) - + u * (sdf(i ,j ,k+2) - sdf(i ,j ,k )); - Real norm = std::sqrt(nx*nx + ny*ny + nz*nz); - if (norm > 0) { - nx *= Real(1)/norm; - ny *= Real(1)/norm; - nz *= Real(1)/norm; - } - pvrtx[3][m] = nx; - pvrtx[4][m] = ny; - pvrtx[5][m] = nz; } }); @@ -846,13 +1117,12 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) auto const& ey_c = ey_fab.const_array(); auto const& ez_c = ez_fab.const_array(); - // 0: # tri, 1: partial sum of # tri, 2: # c_vertex, 3: partial sum of # c_vertex - BaseFab ntri_fab(cbox,4); + BaseFab ntri_fab(cbox,num_cell_data_components); + ntri_fab.setVal(0); auto const& ntri = ntri_fab.array(); GpuArray ptri{nullptr,nullptr,nullptr}; - Gpu::Buffer error({0}); - auto* perror = error.data(); + int* const perror = counters + counter_invalid_triangles; BoxIndexer c_bi(cbox); @@ -860,32 +1130,29 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) [=] AMREX_GPU_DEVICE (int m) { auto [i,j,k] = c_bi(m); int err = 0; - ntri(i,j,k,2) = 0; + ntri(i,j,k,interior_vertex_count) = 0; process_cube(0,lut,i,j,k,sdf_c,ex_c,ey_c,ez_c,pvrtx,ntri,ptri,&err); if (err != 0) { Gpu::Atomic::AddNoRet(perror, 1); } - return ntri(i,j,k,0); + return ntri(i,j,k,triangle_count); }, [=] AMREX_GPU_DEVICE (int m, int ps) { auto [i,j,k] = c_bi(m); - ntri(i,j,k,1) = ps; + ntri(i,j,k,triangle_offset) = ps; }, Scan::Type::exclusive, Scan::retSum); - auto* nerror = error.copyToHost(); - if (*nerror > 0) { - amrex::Abort("Marching Cubes: invalid triangle"); - } - + int const edge_vertex_count = nvx; int nvx_c = Scan::PrefixSum(int(cbox.numPts()), [=] AMREX_GPU_DEVICE (int m) { auto [i,j,k] = c_bi(m); - return ntri(i,j,k,2); + return ntri(i,j,k,interior_vertex_count); }, [=] AMREX_GPU_DEVICE (int m, int ps) { auto [i,j,k] = c_bi(m); - ntri(i,j,k,3) = ps; + ntri(i,j,k,interior_vertex_offset) = + edge_vertex_count + ps; }, Scan::Type::exclusive, Scan::retSum); @@ -903,6 +1170,9 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) { int err = 0; process_cube(1,lut,i,j,k,sdf_c,ex_c,ey_c,ez_c,pvrtx,ntri,ptri,&err); + if (err != 0) { + Gpu::Atomic::AddNoRet(perror, 1); + } }); // Shift vertices @@ -915,16 +1185,934 @@ void marching_cubes (Geometry const& geom, FArrayBox& sdf_fab, MCFab& mc_fab) pvrtx[2][m] = problo[2] + dx[2] * pvrtx[2][m]; }); + // The scratch edge fabs are destroyed on return, so wait for the kernels + // that read them. Both passes accumulate into the same counter: an + // unknown MC33 face id in pass 0 or a triangle referencing an edge without + // a crossing in pass 1 are lookup-table invariants that the driver checks. Gpu::streamSynchronize(); + mc_fab.m_cell_data = std::move(ntri_fab); mc_fab.m_triangles = std::move(tri); - mc_fab.m_vertices = std::move(vrtx); // We can probably release the memory used by nx, ny, nz. + mc_fab.m_vertices = std::move(vrtx); +} + +void build_face_fractions ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf_fab, + FArrayBox& apx_fab, FArrayBox& apy_fab, FArrayBox& apz_fab, + FArrayBox& fcx_fab, FArrayBox& fcy_fab, FArrayBox& fcz_fab, + IArrayBox& rejected_x_fab, IArrayBox& rejected_y_fab, + IArrayBox& rejected_z_fab, int* counters) +{ + BL_PROFILE("MC::build_face_fractions"); + + AMREX_ALWAYS_ASSERT(mc_fab.m_cell_data.box().contains(bx)); + auto const rejected_x = rejected_x_fab.array(); + auto const rejected_y = rejected_y_fab.array(); + auto const rejected_z = rejected_z_fab.array(); + Box const rejected_xbox = rejected_x_fab.box(); + Box const rejected_ybox = rejected_y_fab.box(); + Box const rejected_zbox = rejected_z_fab.box(); + + auto const sdf = sdf_fab.const_array(); + auto const cell_data = mc_fab.m_cell_data.const_array(); + auto const exact_x = mc_fab.m_edge_intersections[0].const_array(); + auto const exact_y = mc_fab.m_edge_intersections[1].const_array(); + auto const exact_z = mc_fab.m_edge_intersections[2].const_array(); + Box const cell_box = mc_fab.m_cell_data.box(); + auto const apx = apx_fab.array(); + auto const apy = apy_fab.array(); + auto const apz = apz_fab.array(); + auto const fcx = fcx_fab.array(); + auto const fcy = fcy_fab.array(); + auto const fcz = fcz_fab.array(); + + // The two cells sharing a face resolved its MC33 ambiguity differently + // (an invariant violation, fatal in the driver), and degenerate face + // polygons that were marked in the rejected face arrays for nodal repair. + int* const error = counters + counter_face_decision_errors; + int* const degenerate = counters + counter_degenerate_faces; + + // Chombo's moment construction starts with boundary-face moments. Build + // those apertures from the same signed-distance edge intersections used + // by marching cubes, so the six Cartesian patches and EB triangles close. + // Edges traversed against their storage direction use 1 - exact; the + // invalid_edge_intersection sentinel (-1) maps to 2, which + // edge_intersection_fraction also treats as "no exact crossing". + ParallelFor(amrex::surroundingNodes(bx,0), + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const face_levelset[4] = { + sdf(i,j ,k ), sdf(i,j+1,k ), + sdf(i,j+1,k+1), sdf(i,j ,k+1) + }; + Real const intersections[4] = { + exact_y(i,j,k), exact_z(i,j+1,k), + 1.0_rt-exact_y(i,j,k+1), 1.0_rt-exact_z(i,j,k) + }; + auto const decision = resolved_face_decision( + cell_data, cell_box, error, i-1,j,k,1, i,j,k,3); + bool const ok = cut_face_fraction(face_levelset, intersections, apx(i,j,k), + fcx(i,j,k,0), fcx(i,j,k,1), + decision[0] != 0, decision[1] != 0); + if (!ok && rejected_xbox.contains(i,j,k)) { + rejected_x(i,j,k) = 1; + Gpu::Atomic::AddNoRet(degenerate, 1); + } + }); + + ParallelFor(amrex::surroundingNodes(bx,1), + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const face_levelset[4] = { + sdf(i ,j,k ), sdf(i+1,j,k ), + sdf(i+1,j,k+1), sdf(i ,j,k+1) + }; + Real const intersections[4] = { + exact_x(i,j,k), exact_z(i+1,j,k), + 1.0_rt-exact_x(i,j,k+1), 1.0_rt-exact_z(i,j,k) + }; + auto const decision = resolved_face_decision( + cell_data, cell_box, error, i,j-1,k,2, i,j,k,0); + bool const ok = cut_face_fraction(face_levelset, intersections, apy(i,j,k), + fcy(i,j,k,0), fcy(i,j,k,1), + decision[0] != 0, decision[1] != 0); + if (!ok && rejected_ybox.contains(i,j,k)) { + rejected_y(i,j,k) = 1; + Gpu::Atomic::AddNoRet(degenerate, 1); + } + }); + + ParallelFor(amrex::surroundingNodes(bx,2), + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const face_levelset[4] = { + sdf(i ,j ,k), sdf(i+1,j ,k), + sdf(i+1,j+1,k), sdf(i ,j+1,k) + }; + Real const intersections[4] = { + exact_x(i,j,k), exact_y(i+1,j,k), + 1.0_rt-exact_x(i,j+1,k), 1.0_rt-exact_y(i,j,k) + }; + auto const decision = resolved_face_decision( + cell_data, cell_box, error, i,j,k-1,5, i,j,k,4); + bool const ok = cut_face_fraction(face_levelset, intersections, apz(i,j,k), + fcz(i,j,k,0), fcz(i,j,k,1), + decision[0] != 0, decision[1] != 0); + if (!ok && rejected_zbox.contains(i,j,k)) { + rejected_z(i,j,k) = 1; + Gpu::Atomic::AddNoRet(degenerate, 1); + } + }); + +} + +void build_edge_centroids ( + Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf_fab, + FArrayBox& ecx_fab, FArrayBox& ecy_fab, FArrayBox& ecz_fab) +{ + BL_PROFILE("MC::build_edge_centroids"); + + Box const exbx = amrex::convert(bx,IntVect(0,1,1)); + Box const eybx = amrex::convert(bx,IntVect(1,0,1)); + Box const ezbx = amrex::convert(bx,IntVect(1,1,0)); + AMREX_ALWAYS_ASSERT(ecx_fab.box().contains(exbx)); + AMREX_ALWAYS_ASSERT(ecy_fab.box().contains(eybx)); + AMREX_ALWAYS_ASSERT(ecz_fab.box().contains(ezbx)); + + auto const sdf = sdf_fab.const_array(); + auto const exact_x = mc_fab.m_edge_intersections[0].const_array(); + auto const exact_y = mc_fab.m_edge_intersections[1].const_array(); + auto const exact_z = mc_fab.m_edge_intersections[2].const_array(); + auto const ecx = ecx_fab.array(); + auto const ecy = ecy_fab.array(); + auto const ecz = ecz_fab.array(); + + ParallelFor(exbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const lo = sdf(i,j,k); + Real const hi = sdf(i+1,j,k); + bool const lo_fluid = lo > 0.0_rt; + bool const hi_fluid = hi > 0.0_rt; + if (lo_fluid && hi_fluid) { + ecx(i,j,k) = 1.0_rt; + } else if (!lo_fluid && !hi_fluid) { + ecx(i,j,k) = -1.0_rt; + } else { + Real const cut = edge_intersection_fraction( + lo, hi, exact_x(i,j,k)); + ecx(i,j,k) = lo_fluid + ? 0.5_rt*cut-0.5_rt + : 0.5_rt*cut; + } + }); + + ParallelFor(eybx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const lo = sdf(i,j,k); + Real const hi = sdf(i,j+1,k); + bool const lo_fluid = lo > 0.0_rt; + bool const hi_fluid = hi > 0.0_rt; + if (lo_fluid && hi_fluid) { + ecy(i,j,k) = 1.0_rt; + } else if (!lo_fluid && !hi_fluid) { + ecy(i,j,k) = -1.0_rt; + } else { + Real const cut = edge_intersection_fraction( + lo, hi, exact_y(i,j,k)); + ecy(i,j,k) = lo_fluid + ? 0.5_rt*cut-0.5_rt + : 0.5_rt*cut; + } + }); + + ParallelFor(ezbx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + Real const lo = sdf(i,j,k); + Real const hi = sdf(i,j,k+1); + bool const lo_fluid = lo > 0.0_rt; + bool const hi_fluid = hi > 0.0_rt; + if (lo_fluid && hi_fluid) { + ecz(i,j,k) = 1.0_rt; + } else if (!lo_fluid && !hi_fluid) { + ecz(i,j,k) = -1.0_rt; + } else { + Real const cut = edge_intersection_fraction( + lo, hi, exact_z(i,j,k)); + ecz(i,j,k) = lo_fluid + ? 0.5_rt*cut-0.5_rt + : 0.5_rt*cut; + } + }); +} + +void build_cell_fractions ( + Box const& bx, Geometry const& geom, MCFab const& mc_fab, + FArrayBox const& sdf_fab, + FArrayBox& apx_fab, FArrayBox& apy_fab, FArrayBox& apz_fab, + FArrayBox& vfrac_fab, FArrayBox& vcent_fab, FArrayBox& barea_fab, + FArrayBox& bcent_fab, FArrayBox& bnorm_fab, int* counters) +{ + BL_PROFILE("MC::build_cell_fractions"); + + AMREX_ALWAYS_ASSERT(mc_fab.m_cell_data.box().contains(bx)); + AMREX_ALWAYS_ASSERT( + sdf_fab.box().contains(amrex::surroundingNodes(bx))); + auto const dx = geom.CellSizeArray(); + + auto const cell_data = mc_fab.m_cell_data.const_array(); + auto const sdf = sdf_fab.const_array(); + auto const apx = apx_fab.array(); + auto const apy = apy_fab.array(); + auto const apz = apz_fab.array(); + auto const vfrac = vfrac_fab.array(); + auto const vcent = vcent_fab.array(); + auto const barea = barea_fab.array(); + auto const bcent = bcent_fab.array(); + auto const bnorm = bnorm_fab.array(); + + auto const* tri_v1 = mc_fab.m_triangles.v1.data(); + auto const* tri_v2 = mc_fab.m_triangles.v2.data(); + auto const* tri_v3 = mc_fab.m_triangles.v3.data(); + auto const* vert_x = mc_fab.m_vertices.x.data(); + auto const* vert_y = mc_fab.m_vertices.y.data(); + auto const* vert_z = mc_fab.m_vertices.z.data(); + + auto const problo = geom.ProbLoArray(); + // Area-vector closure, volume, centroid, and boundary-normal errors in + // four consecutive counter slots. Every quantity checked against + // tolerance is dimensionless. + static_assert(counter_volume_errors == counter_closure_errors + 1 + && counter_centroid_errors == counter_closure_errors + 2 + && counter_area_vector_errors == counter_closure_errors + 3); + int* const errors = counters + counter_closure_errors; + +#ifdef AMREX_USE_FLOAT + constexpr Real tolerance = 2.e-5_rt; +#else + constexpr Real tolerance = 2.e-12_rt; +#endif + + ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + int const triangle_count = + cell_data(i,j,k,CellDataComponent::triangle_count); + int const triangle_offset = + cell_data(i,j,k,CellDataComponent::triangle_offset); + if (triangle_count <= 0) { + return; + } + + Real const cell_lo[3] = { + problo[0] + Real(i)*dx[0], + problo[1] + Real(j)*dx[1], + problo[2] + Real(k)*dx[2] + }; + + // The coordinate-face patches and the EB triangles form a closed + // surface. Work in cell-local coordinates so volume is already a + // volume fraction and first moments normalize directly to AMReX + // centroid coordinates. + Real eb_area_vector[3] = {0.0_rt, 0.0_rt, 0.0_rt}; + Real eb_area = 0.0_rt; + Real eb_centroid_numerator[3] = {0.0_rt, 0.0_rt, 0.0_rt}; + + for (int n = 0; n < triangle_count; ++n) { + int const m = triangle_offset + n; + int const iv1 = tri_v1[m]; + int const iv2 = tri_v2[m]; + int const iv3 = tri_v3[m]; + + Real const p1[3] = {vert_x[iv1], vert_y[iv1], vert_z[iv1]}; + Real const p2[3] = {vert_x[iv2], vert_y[iv2], vert_z[iv2]}; + Real const p3[3] = {vert_x[iv3], vert_y[iv3], vert_z[iv3]}; + Real q1[3]; + Real q2[3]; + Real q3[3]; + for (int d = 0; d < 3; ++d) { + q1[d] = (p1[d]-cell_lo[d])/dx[d] - 0.5_rt; + q2[d] = (p2[d]-cell_lo[d])/dx[d] - 0.5_rt; + q3[d] = (p3[d]-cell_lo[d])/dx[d] - 0.5_rt; + } + + Real const q21[3] = {q2[0]-q1[0], q2[1]-q1[1], q2[2]-q1[2]}; + Real const q31[3] = {q3[0]-q1[0], q3[1]-q1[1], q3[2]-q1[2]}; + Real const local_av[3] = { + 0.5_rt*(q21[1]*q31[2]-q21[2]*q31[1]), + 0.5_rt*(q21[2]*q31[0]-q21[0]*q31[2]), + 0.5_rt*(q21[0]*q31[1]-q21[1]*q31[0]) + }; + for (int d = 0; d < 3; ++d) { + eb_area_vector[d] += local_av[d]; + } + Real const triangle_area = std::sqrt( + local_av[0]*local_av[0] + + local_av[1]*local_av[1] + + local_av[2]*local_av[2]); + eb_area += triangle_area; + for (int d = 0; d < 3; ++d) { + eb_centroid_numerator[d] += + triangle_area*(q1[d]+q2[d]+q3[d])/3.0_rt; + } + } + + Real const expected_area_vector[3] = { + apx(i,j,k)-apx(i+1,j,k), + apy(i,j,k)-apy(i,j+1,k), + apz(i,j,k)-apz(i,j,k+1) + }; + Real const orientation_dot = + eb_area_vector[0]*expected_area_vector[0] + + eb_area_vector[1]*expected_area_vector[1] + + eb_area_vector[2]*expected_area_vector[2]; + Real const orientation = (orientation_dot < 0.0_rt) ? -1.0_rt : 1.0_rt; + if (std::abs(orientation*eb_area_vector[0]-expected_area_vector[0]) > tolerance + || std::abs(orientation*eb_area_vector[1]-expected_area_vector[1]) > tolerance + || std::abs(orientation*eb_area_vector[2]-expected_area_vector[2]) > tolerance) + { + Gpu::Atomic::AddNoRet(errors, 1); + return; + } + + Real volume = ( + 0.5_rt*(apx(i,j,k)+apx(i+1,j,k) + +apy(i,j,k)+apy(i,j+1,k) + +apz(i,j,k)+apz(i,j,k+1)))/3.0_rt; + Real first_moment[3] = { + (apx(i+1,j,k)-apx(i,j,k))/8.0_rt, + (apy(i,j+1,k)-apy(i,j,k))/8.0_rt, + (apz(i,j,k+1)-apz(i,j,k))/8.0_rt + }; + + for (int n = 0; n < triangle_count; ++n) { + int const m = triangle_offset + n; + int const iv1 = tri_v1[m]; + int const iv2 = tri_v2[m]; + int const iv3 = tri_v3[m]; + + Real const q1[3] = { + (vert_x[iv1]-cell_lo[0])/dx[0] - 0.5_rt, + (vert_y[iv1]-cell_lo[1])/dx[1] - 0.5_rt, + (vert_z[iv1]-cell_lo[2])/dx[2] - 0.5_rt + }; + Real const q2[3] = { + (vert_x[iv2]-cell_lo[0])/dx[0] - 0.5_rt, + (vert_y[iv2]-cell_lo[1])/dx[1] - 0.5_rt, + (vert_z[iv2]-cell_lo[2])/dx[2] - 0.5_rt + }; + Real const q3[3] = { + (vert_x[iv3]-cell_lo[0])/dx[0] - 0.5_rt, + (vert_y[iv3]-cell_lo[1])/dx[1] - 0.5_rt, + (vert_z[iv3]-cell_lo[2])/dx[2] - 0.5_rt + }; + Real const q21[3] = {q2[0]-q1[0], q2[1]-q1[1], q2[2]-q1[2]}; + Real const q31[3] = {q3[0]-q1[0], q3[1]-q1[1], q3[2]-q1[2]}; + Real const area_vector[3] = { + orientation*0.5_rt*(q21[1]*q31[2]-q21[2]*q31[1]), + orientation*0.5_rt*(q21[2]*q31[0]-q21[0]*q31[2]), + orientation*0.5_rt*(q21[0]*q31[1]-q21[1]*q31[0]) + }; + Real const triangle_centroid[3] = { + (q1[0]+q2[0]+q3[0])/3.0_rt, + (q1[1]+q2[1]+q3[1])/3.0_rt, + (q1[2]+q2[2]+q3[2])/3.0_rt + }; + volume += (triangle_centroid[0]*area_vector[0] + + triangle_centroid[1]*area_vector[1] + + triangle_centroid[2]*area_vector[2])/3.0_rt; + + for (int d = 0; d < 3; ++d) { + Real const quadratic = + q1[d]*q1[d] + q2[d]*q2[d] + q3[d]*q3[d] + + q1[d]*q2[d] + q2[d]*q3[d] + q3[d]*q1[d]; + first_moment[d] += area_vector[d]*quadratic/12.0_rt; + } + } + + bool nodal_plane_owner = true; + bool nodal_plane_covered = true; + bool has_on_node = false; + for (int kk = 0; kk <= 1; ++kk) { + for (int jj = 0; jj <= 1; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + Real const value = sdf(i+ii,j+jj,k+kk); + nodal_plane_owner = nodal_plane_owner + && value >= 0.0_rt; + nodal_plane_covered = nodal_plane_covered + && value <= 0.0_rt; + has_on_node = has_on_node || value == 0.0_rt; + } + } + } + nodal_plane_owner = nodal_plane_owner && has_on_node; + nodal_plane_covered = nodal_plane_covered && has_on_node; + if (nodal_plane_covered) { + // A coincident patch belongs to the fluid-side cell. On the + // solid side, the zero-area MC triangles are only a classification + // artifact and the cell remains fully covered. + vfrac(i,j,k) = 0.0_rt; + return; + } + if (nodal_plane_owner) { + if (eb_area <= tolerance) { + // The surface only touches this cell at a node or edge. It + // has no measure inside the cell, so the fluid-side cell is + // regular rather than a zero-area cut cell. Write the full + // set of "no boundary in this cell" values explicitly in case + // a cut face still makes the cell single-valued downstream. + vfrac(i,j,k) = 1.0_rt; + barea(i,j,k) = 0.0_rt; + for (int d = 0; d < 3; ++d) { + vcent(i,j,k,d) = 0.0_rt; + bcent(i,j,k,d) = 0.0_rt; + bnorm(i,j,k,d) = 0.0_rt; + } + return; + } + // The EB lies on one or more cell faces and this cell owns the + // coincident patch. Its open volume is still the complete cell. + volume = 1.0_rt; + first_moment[0] = 0.0_rt; + first_moment[1] = 0.0_rt; + first_moment[2] = 0.0_rt; + } + + if (volume <= tolerance || volume > 1.0_rt+tolerance + || eb_area <= tolerance) + { + Gpu::Atomic::AddNoRet(errors+1, 1); + return; + } + volume = amrex::min(amrex::max(volume,0.0_rt),1.0_rt); + + Real centroid[3] = { + first_moment[0]/volume, + first_moment[1]/volume, + first_moment[2]/volume + }; + if (centroid[0] < -0.5_rt-tolerance || centroid[0] > 0.5_rt+tolerance + || centroid[1] < -0.5_rt-tolerance || centroid[1] > 0.5_rt+tolerance + || centroid[2] < -0.5_rt-tolerance || centroid[2] > 0.5_rt+tolerance) + { + Gpu::Atomic::AddNoRet(errors+2, 1); + return; + } + + Real boundary_normal[3] = { + orientation*eb_area_vector[0], + orientation*eb_area_vector[1], + orientation*eb_area_vector[2] + }; + Real const area_vector_norm = std::sqrt( + boundary_normal[0]*boundary_normal[0] + + boundary_normal[1]*boundary_normal[1] + + boundary_normal[2]*boundary_normal[2]); + if (area_vector_norm <= tolerance) { + Gpu::Atomic::AddNoRet(errors+3, 1); + return; + } + for (Real& normal_component : boundary_normal) { + normal_component /= area_vector_norm; + } + + Real boundary_centroid[3] = { + eb_centroid_numerator[0]/eb_area, + eb_centroid_numerator[1]/eb_area, + eb_centroid_numerator[2]/eb_area + }; + if (boundary_centroid[0] < -0.5_rt-tolerance + || boundary_centroid[0] > 0.5_rt+tolerance + || boundary_centroid[1] < -0.5_rt-tolerance + || boundary_centroid[1] > 0.5_rt+tolerance + || boundary_centroid[2] < -0.5_rt-tolerance + || boundary_centroid[2] > 0.5_rt+tolerance) + { + Gpu::Atomic::AddNoRet(errors+2, 1); + return; + } + + vfrac(i,j,k) = volume; + for (int d = 0; d < 3; ++d) { + vcent(i,j,k,d) = amrex::min(amrex::max(centroid[d],-0.5_rt),0.5_rt); + bcent(i,j,k,d) = amrex::Clamp( + boundary_centroid[d], -0.5_rt, 0.5_rt); + bnorm(i,j,k,d) = boundary_normal[d]; + } + barea(i,j,k) = eb_area; + }); + +} + +int build_cell_topology (Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf_fab, + EBCellFlagFab& cellflag, FArrayBox const& vfrac_fab, + FArrayBox const& apx_fab, FArrayBox const& apy_fab, + FArrayBox const& apz_fab) +{ + BL_PROFILE("MC::build_cell_topology"); + + Box const bxg1 = amrex::grow(bx, 1); + Box const nbxg1 = amrex::surroundingNodes(bxg1); + IntVect const x_transverse_ghost = + IntVect::TheUnitVector() - IntVect::TheDimensionVector(0); + IntVect const y_transverse_ghost = + IntVect::TheUnitVector() - IntVect::TheDimensionVector(1); + IntVect const z_transverse_ghost = + IntVect::TheUnitVector() - IntVect::TheDimensionVector(2); + Box const valid_fxbx = amrex::surroundingNodes(bxg1, 0); + Box const valid_fybx = amrex::surroundingNodes(bxg1, 1); + Box const valid_fzbx = amrex::surroundingNodes(bxg1, 2); + Box const fxbx = amrex::grow(valid_fxbx, x_transverse_ghost); + Box const fybx = amrex::grow(valid_fybx, y_transverse_ghost); + Box const fzbx = amrex::grow(valid_fzbx, z_transverse_ghost); + + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(nbxg1)); + AMREX_ALWAYS_ASSERT(cellflag.box().contains(bxg1)); + AMREX_ALWAYS_ASSERT(vfrac_fab.box().contains(bxg1)); + AMREX_ALWAYS_ASSERT(mc_fab.m_cell_data.box().contains(bx)); + AMREX_ALWAYS_ASSERT(apx_fab.box().contains(valid_fxbx)); + AMREX_ALWAYS_ASSERT(apy_fab.box().contains(valid_fybx)); + AMREX_ALWAYS_ASSERT(apz_fab.box().contains(valid_fzbx)); + + auto const sdf = sdf_fab.const_array(); + auto const cell = cellflag.array(); + auto const cell_data = mc_fab.m_cell_data.const_array(); + auto const vfrac = vfrac_fab.const_array(); + auto const apx = apx_fab.const_array(); + auto const apy = apy_fab.const_array(); + auto const apz = apz_fab.const_array(); + + BaseFab fx_fab(fxbx); + BaseFab fy_fab(fybx); + BaseFab fz_fab(fzbx); + fx_fab.setVal(EB2::Type::regular); + fy_fab.setVal(EB2::Type::regular); + fz_fab.setVal(EB2::Type::regular); + auto const fx = fx_fab.array(); + auto const fy = fy_fab.array(); + auto const fz = fz_fab.array(); + +#ifdef AMREX_USE_FLOAT + constexpr Real tolerance = 2.e-5_rt; +#else + constexpr Real tolerance = 2.e-12_rt; +#endif + + ParallelFor(bxg1, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool const all_faces_open = + apx(i, j, k) >= 1.0_rt - tolerance && apx(i + 1, j, k) >= 1.0_rt - tolerance && + apy(i, j, k) >= 1.0_rt - tolerance && apy(i, j + 1, k) >= 1.0_rt - tolerance && + apz(i, j, k) >= 1.0_rt - tolerance && apz(i, j, k + 1) >= 1.0_rt - tolerance; + + if (vfrac(i, j, k) <= tolerance) { + cell(i, j, k).setCovered(); + } else if (vfrac(i, j, k) >= 1.0_rt - tolerance && all_faces_open) { + cell(i, j, k).setRegular(); + } else { + cell(i, j, k).setSingleValued(); + } + }); + + ParallelFor(valid_fxbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + fx(i, j, k) = face_type(apx(i, j, k), tolerance); + }); + ParallelFor(valid_fybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + fy(i, j, k) = face_type(apy(i, j, k), tolerance); + }); + ParallelFor(valid_fzbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + fz(i, j, k) = face_type(apz(i, j, k), tolerance); + }); + + Gpu::DeviceScalar error_count(0); + int* const errors = error_count.dataPtr(); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const triangle_count = cell_data(i, j, k, CellDataComponent::triangle_count); + bool const has_triangles = triangle_count > 0; + bool const source_has_fluid = sdf(i, j, k) > 0.0_rt || sdf(i + 1, j, k) > 0.0_rt || + sdf(i, j + 1, k) > 0.0_rt || sdf(i + 1, j + 1, k) > 0.0_rt || + sdf(i, j, k + 1) > 0.0_rt || sdf(i + 1, j, k + 1) > 0.0_rt || + sdf(i, j + 1, k + 1) > 0.0_rt || + sdf(i + 1, j + 1, k + 1) > 0.0_rt; + bool const source_has_covered = + sdf(i, j, k) <= 0.0_rt || sdf(i + 1, j, k) <= 0.0_rt || sdf(i, j + 1, k) <= 0.0_rt || + sdf(i + 1, j + 1, k) <= 0.0_rt || sdf(i, j, k + 1) <= 0.0_rt || + sdf(i + 1, j, k + 1) <= 0.0_rt || sdf(i, j + 1, k + 1) <= 0.0_rt || + sdf(i + 1, j + 1, k + 1) <= 0.0_rt; + bool const source_is_cut = source_has_fluid && source_has_covered; + if (has_triangles != source_is_cut) { + Gpu::Atomic::AddNoRet(errors, 1); + } + }); + + EBCellFlagFab cellflagtmp(cellflag.box(), 1); + Elixir cellflagtmp_eli = cellflagtmp.elixir(); + EB2::set_connection_flags(bx, bxg1, cell, cellflagtmp.array(), fx, fy, fz); + + Gpu::streamSynchronize(); + return error_count.dataValue(); +} + +void mark_faces_for_cleanup (Box const& bx, MCFab const& mc_fab, FArrayBox const& sdf_fab, + IArrayBox& rejected_x_fab, IArrayBox& rejected_y_fab, + IArrayBox& rejected_z_fab, int* counters) +{ + BL_PROFILE("MC::mark_faces_for_cleanup"); + + Box const xbx = amrex::surroundingNodes(bx, 0); + Box const ybx = amrex::surroundingNodes(bx, 1); + Box const zbx = amrex::surroundingNodes(bx, 2); + AMREX_ALWAYS_ASSERT(mc_fab.m_cell_data.box().contains(bx)); + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(amrex::surroundingNodes(bx))); + AMREX_ALWAYS_ASSERT(rejected_x_fab.box().contains(xbx)); + AMREX_ALWAYS_ASSERT(rejected_y_fab.box().contains(ybx)); + AMREX_ALWAYS_ASSERT(rejected_z_fab.box().contains(zbx)); + + auto const cell_data = mc_fab.m_cell_data.const_array(); + Box const cell_box = mc_fab.m_cell_data.box(); + auto const sdf = sdf_fab.const_array(); + auto const rejected_x = rejected_x_fab.array(); + auto const rejected_y = rejected_y_fab.array(); + auto const rejected_z = rejected_z_fab.array(); + + int* const count = counters + counter_face_rejections; + + ParallelFor(xbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool const rejected = + face_is_rejected(sdf(i, j, k), sdf(i, j + 1, k), sdf(i, j + 1, k + 1), sdf(i, j, k + 1), + cell_data, cell_box, i - 1, j, k, 1, i, j, k, 3); + if (rejected && rejected_x(i, j, k) == 0) { + rejected_x(i, j, k) = 1; + Gpu::Atomic::AddNoRet(count, 1); + } + }); + ParallelFor(ybx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool const rejected = + face_is_rejected(sdf(i, j, k), sdf(i + 1, j, k), sdf(i + 1, j, k + 1), sdf(i, j, k + 1), + cell_data, cell_box, i, j - 1, k, 2, i, j, k, 0); + if (rejected && rejected_y(i, j, k) == 0) { + rejected_y(i, j, k) = 1; + Gpu::Atomic::AddNoRet(count, 1); + } + }); + ParallelFor(zbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool const rejected = + face_is_rejected(sdf(i, j, k), sdf(i + 1, j, k), sdf(i + 1, j + 1, k), sdf(i, j + 1, k), + cell_data, cell_box, i, j, k - 1, 5, i, j, k, 4); + if (rejected && rejected_z(i, j, k) == 0) { + rejected_z(i, j, k) = 1; + Gpu::Atomic::AddNoRet(count, 1); + } + }); + +} + +void zero_nodes_for_cleanup (Box const& node_box, IArrayBox const& rejected_cells_fab, + IArrayBox const& rejected_x_fab, IArrayBox const& rejected_y_fab, + IArrayBox const& rejected_z_fab, FArrayBox& sdf_fab, int* counters) +{ + BL_PROFILE("MC::zero_nodes_for_cleanup"); + + // The marks must cover every cell and face enclosed by node_box. The + // outermost cells and faces incident to its boundary nodes may lie beyond + // the arrays; they are ignored, as in the legacy generator. + Box const cell_box = amrex::enclosedCells(node_box); + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(node_box)); + AMREX_ALWAYS_ASSERT(rejected_cells_fab.box().contains(cell_box)); + AMREX_ALWAYS_ASSERT(rejected_x_fab.box().contains(amrex::surroundingNodes(cell_box, 0))); + AMREX_ALWAYS_ASSERT(rejected_y_fab.box().contains(amrex::surroundingNodes(cell_box, 1))); + AMREX_ALWAYS_ASSERT(rejected_z_fab.box().contains(amrex::surroundingNodes(cell_box, 2))); + auto const rejected_cells = rejected_cells_fab.const_array(); + auto const rejected_x = rejected_x_fab.const_array(); + auto const rejected_y = rejected_y_fab.const_array(); + auto const rejected_z = rejected_z_fab.const_array(); + auto const sdf = sdf_fab.array(); + + int* const changed = counters + counter_changed_nodes; + ParallelFor(node_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (sdf(i, j, k) <= 0.0_rt) { + return; + } + + bool rejected = false; + for (int kk = 0; kk <= 1 && !rejected; ++kk) { + for (int jj = 0; jj <= 1 && !rejected; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + rejected = rejected || (rejected_cells.contains(i - ii, j - jj, k - kk) + && rejected_cells(i - ii, j - jj, k - kk) != 0); + } + } + } + for (int kk = 0; kk <= 1 && !rejected; ++kk) { + for (int jj = 0; jj <= 1; ++jj) { + rejected = rejected || (rejected_x.contains(i, j - jj, k - kk) + && rejected_x(i, j - jj, k - kk) != 0); + } + } + for (int kk = 0; kk <= 1 && !rejected; ++kk) { + for (int ii = 0; ii <= 1; ++ii) { + rejected = rejected || (rejected_y.contains(i - ii, j, k - kk) + && rejected_y(i - ii, j, k - kk) != 0); + } + } + for (int jj = 0; jj <= 1 && !rejected; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + rejected = rejected || (rejected_z.contains(i - ii, j - jj, k) + && rejected_z(i - ii, j - jj, k) != 0); + } + } + if (rejected) { + sdf(i, j, k) = 0.0_rt; + Gpu::Atomic::AddNoRet(changed, 1); + } + }); +} + +void extend_domain_face_levelset (Box const& node_box, Box const& domain, + GpuArray const& is_periodic, FArrayBox& sdf_fab, + int* counters) +{ + BL_PROFILE("MC::extend_domain_face_levelset"); + + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(node_box)); + Box const nodal_domain = amrex::surroundingNodes(domain); + Box reference_box = node_box; + for (int d = 0; d < 3; ++d) { + if (!is_periodic[d]) { + reference_box.setSmall( + d, amrex::Clamp(node_box.smallEnd(d), nodal_domain.smallEnd(d), + nodal_domain.bigEnd(d))); + reference_box.setBig( + d, amrex::Clamp(node_box.bigEnd(d), nodal_domain.smallEnd(d), + nodal_domain.bigEnd(d))); + } + } + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(reference_box)); + auto const sdf = sdf_fab.array(); + + int const domlo_x = nodal_domain.smallEnd(0); + int const domlo_y = nodal_domain.smallEnd(1); + int const domlo_z = nodal_domain.smallEnd(2); + int const domhi_x = nodal_domain.bigEnd(0); + int const domhi_y = nodal_domain.bigEnd(1); + int const domhi_z = nodal_domain.bigEnd(2); + + int* const changed = counters + counter_extended_nodes; + ParallelFor(node_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const ii = is_periodic[0] ? i : amrex::Clamp(i, domlo_x, domhi_x); + int const jj = is_periodic[1] ? j : amrex::Clamp(j, domlo_y, domhi_y); + int const kk = is_periodic[2] ? k : amrex::Clamp(k, domlo_z, domhi_z); + int const outside_directions = (ii != i) + (jj != j) + (kk != k); + if (outside_directions != 1) { + return; + } + Real const extended_value = sdf(ii, jj, kk); + // A repaired (exact-zero) exterior node stays covered: re-extending a + // fluid value onto it would undo the repair and could keep the repair + // loop from converging. + if (sdf(i, j, k) == 0.0_rt && extended_value > 0.0_rt) { + return; + } + if (sdf(i, j, k) != extended_value) { + sdf(i, j, k) = extended_value; + Gpu::Atomic::AddNoRet(changed, 1); + } + }); +} + +void extend_domain_face_edge_intersections (Box const& domain, + GpuArray const& is_periodic, + MCFab& mc_fab) +{ + BL_PROFILE("MC::extend_domain_face_edge_intersections"); + + GpuArray const domain_lo{ + domain.smallEnd(0), domain.smallEnd(1), domain.smallEnd(2)}; + GpuArray const domain_hi{ + domain.bigEnd(0), domain.bigEnd(1), domain.bigEnd(2)}; + for (int edge_direction = 0; edge_direction < 3; ++edge_direction) { + auto const crossing = mc_fab.m_edge_intersections[edge_direction].array(); + Box const edge_box = mc_fab.m_edge_intersections[edge_direction].box(); + ParallelFor(edge_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int index[3] = {i, j, k}; + int reference[3] = {i, j, k}; + int outside_direction = -1; + int outside_directions = 0; + for (int d = 0; d < 3; ++d) { + if (is_periodic[d]) { + continue; + } + int const lo = domain_lo[d]; + int const hi = domain_hi[d] + (d != edge_direction); + reference[d] = amrex::Clamp(index[d], lo, hi); + if (reference[d] != index[d]) { + outside_direction = d; + ++outside_directions; + } + } + if (outside_directions == 1 && outside_direction != edge_direction) { + crossing(i, j, k) = crossing(reference[0], reference[1], reference[2]); + } + }); + } +} + +void mark_cells_for_cleanup (Box const& bx, MCFab const& mc_fab, + FArrayBox const& sdf_fab, FArrayBox const& vfrac_fab, + Real small_volfrac, IArrayBox& rejected_fab, int* counters) +{ + BL_PROFILE("MC::mark_cells_for_cleanup"); + + AMREX_ALWAYS_ASSERT(mc_fab.m_cell_data.box().contains(bx)); + AMREX_ALWAYS_ASSERT(vfrac_fab.box().contains(bx)); + AMREX_ALWAYS_ASSERT(rejected_fab.box().contains(bx)); + AMREX_ALWAYS_ASSERT(sdf_fab.box().contains(amrex::surroundingNodes(bx))); + + auto const cell_data = mc_fab.m_cell_data.const_array(); + auto const sdf = sdf_fab.const_array(); + auto const vfrac = vfrac_fab.const_array(); + auto const rejected = rejected_fab.array(); + + static_assert(counter_small_cell_rejections == counter_topology_rejections + 1); + int* const counts = counters + counter_topology_rejections; + + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + Real const cube[8] = { + sdf(i, j, k), sdf(i + 1, j, k), sdf(i + 1, j + 1, k), sdf(i, j + 1, k), + sdf(i, j, k + 1), sdf(i + 1, j, k + 1), sdf(i + 1, j + 1, k + 1), sdf(i, j + 1, k + 1)}; + bool fluid[8]; + int nfluid = 0; + for (int n = 0; n < 8; ++n) { + fluid[n] = cube[n] > 0.0_rt; + nfluid += fluid[n]; + } + + int const triangle_count = cell_data(i, j, k, CellDataComponent::triangle_count); + bool const is_cut = nfluid > 0 && nfluid < 8; + bool bad_topology = is_cut && triangle_count <= 0; + + // Count the fluid corner groups of the cell. Corners are joined along + // the 12 cube edges and across every ambiguous (four-crossing) face + // whose MC33 decision says the two diagonal fluid corners are + // connected: that decision is what the extracted surface and the face + // apertures already use, so a group counted here is exactly a fluid + // region that touches the cell faces. Corner groups that MC33 joins + // only through the cell interior (the tunnel tilings 4.1.2/10.1.2 and + // the 7.3/10.2/12.2/13.x variants) are deliberately NOT merged: a + // single-valued EB cell may hold one face-connected fluid region only, + // so such cells are rejected and repaired. + int corner_parent[8]; + for (int n = 0; n < 8; ++n) { + corner_parent[n] = n; + } + auto find_root = [&] (int n) { + while (corner_parent[n] != n) { n = corner_parent[n]; } + return n; + }; + auto join = [&] (int a, int b) { + int const ra = find_root(a); + int const rb = find_root(b); + if (ra != rb) { corner_parent[rb] = ra; } + }; + constexpr int edge_lo[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3}; + constexpr int edge_hi[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7}; + for (int n = 0; n < 12; ++n) { + if (fluid[edge_lo[n]] && fluid[edge_hi[n]]) { + join(edge_lo[n], edge_hi[n]); + } + } + // MC33 face ids 1..6 in cube-corner order, matching test_face()'s + // A,B,C,D and the bits stored in m_cell_data. + constexpr int face_corner[6][4] = {{0, 4, 5, 1}, {1, 5, 6, 2}, {2, 6, 7, 3}, + {3, 7, 4, 0}, {0, 3, 2, 1}, {4, 7, 6, 5}}; + int const valid_mask = cell_data(i, j, k, CellDataComponent::face_decision_valid_mask); + int const connected_mask = + cell_data(i, j, k, CellDataComponent::face_fluid_connected_mask); + for (int f = 0; f < 6; ++f) { + int const c0 = face_corner[f][0]; + int const c1 = face_corner[f][1]; + int const c2 = face_corner[f][2]; + int const c3 = face_corner[f][3]; + if (!(fluid[c0] == fluid[c2] && fluid[c1] == fluid[c3] && fluid[c0] != fluid[c1])) { + continue; // not an ambiguous face + } + int const bit = 1 << f; + Real const face_levelset[4] = {cube[c0], cube[c1], cube[c2], cube[c3]}; + bool const connected = ((valid_mask & bit) != 0) + ? ((connected_mask & bit) != 0) + : four_crossing_fluid_is_connected(face_levelset); + if (connected) { + if (fluid[c0]) { join(c0, c2); } else { join(c1, c3); } + } + } + int fluid_components = 0; + for (int n = 0; n < 8; ++n) { + if (fluid[n]) { + fluid_components += find_root(n) == n; + } + } + bad_topology = bad_topology || fluid_components > 1; + + // A negative sentinel means geometry construction rejected the closed + // boundary. Cover it and let the next MC pass rebuild its neighbors. + bad_topology = bad_topology || (is_cut && vfrac(i, j, k) < 0.0_rt); + bool const small_cell = is_cut && !bad_topology && vfrac(i, j, k) < small_volfrac; + + rejected(i, j, k) = bad_topology ? RejectionReason::invalid_topology + : (small_cell ? RejectionReason::small_volume : 0); + if (bad_topology) { + Gpu::Atomic::AddNoRet(counts, 1); + } else if (small_cell) { + Gpu::Atomic::AddNoRet(counts + 1, 1); + } + }); + } -void write_stl (std::string const& filename, std::map> const& mc_fabs) +void write_stl (std::string const& filename, LayoutData const& mc_fabs) { - int myproc = ParallelDescriptor::MyProc(); - int nprocs = ParallelDescriptor::NProcs(); + BoxArray const& grids = mc_fabs.boxArray(); + // The EB may be built inside a ParallelContext sub-frame, so the token + // chain runs over the sub-communicator like the rest of the builder. + int const myproc = ParallelContext::MyProcSub(); + int const nprocs = ParallelContext::NProcsSub(); std::ofstream ofs; @@ -936,16 +2124,23 @@ void write_stl (std::string const& filename, std::map #ifdef AMREX_USE_MPI if (myproc > 0) { int foo = 0; - ParallelDescriptor::Recv(&foo, 1, myproc-1, 100); + ParallelDescriptor::Recv(&foo, 1, myproc - 1, 100, ParallelContext::CommunicatorSub()); } #endif - if (! ofs.is_open()) { + if (!ofs.is_open()) { ofs.open(filename, std::ios_base::app); } - - for (auto const& [k,p] : mc_fabs) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ofs.good(), + "Could not open marching-cubes STL output " + filename); + ofs << std::setprecision(std::numeric_limits::max_digits10); + + for (MFIter mfi(mc_fabs); mfi.isValid(); ++mfi) { + int const k = mfi.index(); + MCFab const* p = &mc_fabs[mfi]; + AMREX_ALWAYS_ASSERT(p->m_cell_data.box().contains(grids[k])); auto ntri = int(p->m_triangles.v1.size()); + amrex::ignore_unused(ntri); #ifdef AMREX_USE_GPU Gpu::PinnedVector tri_v1(ntri); @@ -955,12 +2150,21 @@ void write_stl (std::string const& filename, std::map Gpu::PinnedVector vert_x(nvert); Gpu::PinnedVector vert_y(nvert); Gpu::PinnedVector vert_z(nvert); - Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v1.begin(), p->m_triangles.v1.end(), tri_v1.begin()); - Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v2.begin(), p->m_triangles.v2.end(), tri_v2.begin()); - Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v3.begin(), p->m_triangles.v3.end(), tri_v3.begin()); - Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.x.begin(), p->m_vertices.x.end(), vert_x.begin()); - Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.y.begin(), p->m_vertices.y.end(), vert_y.begin()); - Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.z.begin(), p->m_vertices.z.end(), vert_z.begin()); + BaseFab cell_data(p->m_cell_data.box(), p->m_cell_data.nComp(), The_Pinned_Arena()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v1.begin(), p->m_triangles.v1.end(), + tri_v1.begin()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v2.begin(), p->m_triangles.v2.end(), + tri_v2.begin()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_triangles.v3.begin(), p->m_triangles.v3.end(), + tri_v3.begin()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.x.begin(), p->m_vertices.x.end(), + vert_x.begin()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.y.begin(), p->m_vertices.y.end(), + vert_y.begin()); + Gpu::copyAsync(Gpu::deviceToHost, p->m_vertices.z.begin(), p->m_vertices.z.end(), + vert_z.begin()); + Gpu::dtoh_memcpy_async(cell_data.dataPtr(), p->m_cell_data.dataPtr(), + p->m_cell_data.nBytes(p->m_cell_data.box(), p->m_cell_data.nComp())); Gpu::streamSynchronize(); #else auto const& tri_v1 = p->m_triangles.v1; @@ -969,41 +2173,61 @@ void write_stl (std::string const& filename, std::map auto const& vert_x = p->m_vertices.x; auto const& vert_y = p->m_vertices.y; auto const& vert_z = p->m_vertices.z; + auto const& cell_data = p->m_cell_data; #endif - for (int itri = 0; itri < ntri; ++itri) { - auto iv1 = tri_v1[itri]; - auto iv2 = tri_v2[itri]; - auto iv3 = tri_v3[itri]; - XDim3 v1{.x = vert_x[iv1], .y = vert_y[iv1], .z = vert_z[iv1]}; - XDim3 v2{.x = vert_x[iv2], .y = vert_y[iv2], .z = vert_z[iv2]}; - XDim3 v3{.x = vert_x[iv3], .y = vert_y[iv3], .z = vert_z[iv3]}; - XDim3 vec1{.x = v2.x-v1.x, .y = v2.y-v1.y, .z = v2.z-v1.z}; - XDim3 vec2{.x = v3.x-v2.x, .y = v3.y-v2.y, .z = v3.z-v2.z}; - XDim3 norm{.x = vec1.y*vec2.z-vec1.z*vec2.y, - .y = vec1.z*vec2.x-vec1.x*vec2.z, - .z = vec1.x*vec2.y-vec1.y*vec2.x}; - auto tmp = std::sqrt(norm.x*norm.x + norm.y*norm.y + norm.z*norm.z); - if (tmp != 0) { tmp = Real(1) / tmp; } - ofs << "facet normal " << norm.x*tmp << " " << norm.y*tmp << " " << norm.z*tmp << "\n" - << " outer loop\n" - << " vertex " << v1.x << " " << v1.y << " " << v1.z << "\n" - << " vertex " << v2.x << " " << v2.y << " " << v2.z << "\n" - << " vertex " << v3.x << " " << v3.y << " " << v3.z << "\n" - << " endloop\n" - << "endfacet\n"; - } + auto const cell = cell_data.const_array(); + amrex::LoopOnCpu(grids[k], [&] (int i, int j, int kk) noexcept { + int const count = cell(i, j, kk, 0); + int const offset = cell(i, j, kk, 1); + for (int n = 0; n < count; ++n) { + int const itri = offset + n; + AMREX_ASSERT(itri >= 0 && itri < ntri); + auto iv1 = tri_v1[itri]; + auto iv2 = tri_v2[itri]; + auto iv3 = tri_v3[itri]; + XDim3 v1{.x = vert_x[iv1], .y = vert_y[iv1], .z = vert_z[iv1]}; + XDim3 v2{.x = vert_x[iv2], .y = vert_y[iv2], .z = vert_z[iv2]}; + XDim3 v3{.x = vert_x[iv3], .y = vert_y[iv3], .z = vert_z[iv3]}; + XDim3 vec1{.x = v2.x - v1.x, .y = v2.y - v1.y, .z = v2.z - v1.z}; + XDim3 vec2{.x = v3.x - v2.x, .y = v3.y - v2.y, .z = v3.z - v2.z}; + XDim3 norm{.x = vec1.y * vec2.z - vec1.z * vec2.y, + .y = vec1.z * vec2.x - vec1.x * vec2.z, + .z = vec1.x * vec2.y - vec1.y * vec2.x}; + auto tmp = std::sqrt(norm.x * norm.x + norm.y * norm.y + norm.z * norm.z); + Real const edge_scale_sq = + amrex::max(vec1.x * vec1.x + vec1.y * vec1.y + vec1.z * vec1.z, + vec2.x * vec2.x + vec2.y * vec2.y + vec2.z * vec2.z); + Real const degenerate_tolerance = + 64.0_rt * std::numeric_limits::epsilon() * edge_scale_sq; + if (tmp <= degenerate_tolerance) { + continue; + } + tmp = Real(1) / tmp; + ofs << "facet normal " << norm.x * tmp << " " << norm.y * tmp << " " << norm.z * tmp + << "\n" + << " outer loop\n" + << " vertex " << v1.x << " " << v1.y << " " << v1.z << "\n" + << " vertex " << v2.x << " " << v2.y << " " << v2.z << "\n" + << " vertex " << v3.x << " " << v3.y << " " << v3.z << "\n" + << " endloop\n" + << "endfacet\n"; + } + }); } + if (myproc == nprocs - 1) { + ofs << "endsolid Created by AMReX\n"; + } + ofs.close(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(!ofs.fail(), + "Could not complete marching-cubes STL output " + filename); + #ifdef AMREX_USE_MPI - if (myproc < nprocs-1) { + if (myproc < nprocs - 1) { int foo = 0; - ParallelDescriptor::Send(&foo, 1, myproc+1, 100); + ParallelDescriptor::Send(&foo, 1, myproc + 1, 100, ParallelContext::CommunicatorSub()); } + ParallelDescriptor::Barrier(ParallelContext::CommunicatorSub()); #endif - - if (myproc == nprocs-1) { - ofs << "endsolid Created by AMReX\n"; - } -} - } +} // namespace amrex::MC diff --git a/Src/EB/AMReX_WriteEBSurface.H b/Src/EB/AMReX_WriteEBSurface.H index ac27516e937..2c349176ba9 100644 --- a/Src/EB/AMReX_WriteEBSurface.H +++ b/Src/EB/AMReX_WriteEBSurface.H @@ -5,10 +5,12 @@ #include #include +#include + /** * \file AMReX_WriteEBSurface.H * - * Writes EB surface geometry (in pvd format) for diagnostic visualization. + * Writes EB surface geometry for diagnostic visualization. */ namespace amrex { @@ -21,6 +23,15 @@ class EBFArrayBoxFactory; void WriteEBSurface (const amrex::BoxArray & ba, const amrex::DistributionMapping & dmap, const amrex::Geometry & geom, const amrex::EBFArrayBoxFactory * ebf); +/** + * \brief Triangulate and write the planar EB reconstruction to an ASCII STL file. + */ +void WriteEBSurfaceSTL (const amrex::BoxArray& ba, + const amrex::DistributionMapping& dmap, + const amrex::Geometry& geom, + const amrex::EBFArrayBoxFactory* ebf, + std::string const& filename); + } #endif diff --git a/Src/EB/AMReX_WriteEBSurface.cpp b/Src/EB/AMReX_WriteEBSurface.cpp index 2e94ccb7ffb..079bcc92b9f 100644 --- a/Src/EB/AMReX_WriteEBSurface.cpp +++ b/Src/EB/AMReX_WriteEBSurface.cpp @@ -9,16 +9,17 @@ namespace amrex { -void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, const Geometry & geom, - const EBFArrayBoxFactory * ebf) { - +namespace { +//! Collect the planar EB reconstruction of every cut cell of \p mf_ba, an +//! EB-factory-backed MultiFab. +EBToPVD make_eb_surface (const MultiFab& mf_ba, const Geometry& geom, + const EBFArrayBoxFactory* ebf) +{ EBToPVD eb_to_pvd; const Real* dx = geom.CellSize(); const Real* problo = geom.ProbLo(); - MultiFab mf_ba(ba, dmap, 1, 0, MFInfo(), *ebf); - for (MFIter mfi(mf_ba); mfi.isValid(); ++mfi) { const auto & sfab = static_cast(mf_ba[mfi]); @@ -81,6 +82,16 @@ void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, cons areafrac[2]->const_array()); } + return eb_to_pvd; +} +} + +void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, const Geometry & geom, + const EBFArrayBoxFactory * ebf) { + + MultiFab mf_ba(ba, dmap, 1, 0, MFInfo(), *ebf); + EBToPVD eb_to_pvd = make_eb_surface(mf_ba,geom,ebf); + int cpu = ParallelDescriptor::MyProc(); int nProcs = ParallelDescriptor::NProcs(); @@ -90,6 +101,9 @@ void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, cons EBToPVD::WritePVTP(nProcs); } + const Real* dx = geom.CellSize(); + const Real* problo = geom.ProbLo(); + for (MFIter mfi(mf_ba); mfi.isValid(); ++mfi) { const auto & sfab = static_cast(mf_ba[mfi]); @@ -117,4 +131,13 @@ void WriteEBSurface (const BoxArray & ba, const DistributionMapping & dmap, cons } } +void WriteEBSurfaceSTL (const BoxArray& ba, const DistributionMapping& dmap, + const Geometry& geom, const EBFArrayBoxFactory* ebf, + std::string const& filename) +{ + MultiFab mf_ba(ba, dmap, 1, 0, MFInfo(), *ebf); + EBToPVD eb_to_pvd = make_eb_surface(mf_ba,geom,ebf); + eb_to_pvd.WriteSTL(filename); +} + } diff --git a/Src/EB/CMakeLists.txt b/Src/EB/CMakeLists.txt index 13a8c95fb8b..61653641c41 100644 --- a/Src/EB/CMakeLists.txt +++ b/Src/EB/CMakeLists.txt @@ -113,6 +113,8 @@ foreach(D IN LISTS AMReX_SPACEDIM) AMReX_mc_jgt_table.H AMReX_MarchingCubes.H AMReX_MarchingCubes.cpp + AMReX_EB2_Level_MC.H + AMReX_EB2_Level_MC.cpp ) endif () diff --git a/Src/EB/Make.package b/Src/EB/Make.package index a61923f1976..d6af355b0d0 100644 --- a/Src/EB/Make.package +++ b/Src/EB/Make.package @@ -102,8 +102,8 @@ CEXE_sources += AMReX_EB2_Level_chkpt_file.cpp AMReX_EB2_IndexSpace_chkpt_file.c ifeq ($(DIM),3) CEXE_sources += AMReX_WriteEBSurface.cpp AMReX_EBToPVD.cpp CEXE_headers += AMReX_WriteEBSurface.H AMReX_EBToPVD.H - CEXE_sources += AMReX_MarchingCubes.cpp - CEXE_headers += AMReX_MarchingCubes.H AMReX_mc_jgt_table.H + CEXE_sources += AMReX_MarchingCubes.cpp AMReX_EB2_Level_MC.cpp + CEXE_headers += AMReX_MarchingCubes.H AMReX_mc_jgt_table.H AMReX_EB2_Level_MC.H endif ifneq ($(BL_NO_FORT),TRUE) diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt index 26b98a39198..2ccea884678 100644 --- a/Tests/CMakeLists.txt +++ b/Tests/CMakeLists.txt @@ -3,7 +3,7 @@ # function (setup_test _dim _srcs _inputs) - cmake_parse_arguments( "" "HAS_FORTRAN_MODULES" + cmake_parse_arguments( "" "HAS_FORTRAN_MODULES;COMBINE_MPI_OMP" "BASE_NAME;RUNTIME_SUBDIR;EXTRA_DEFINITIONS;CMDLINE_PARAMS;NTASKS;NTHREADS" "CMDLINE_PARAMS_AFTER_INPUTS" ${ARGN} ) @@ -79,7 +79,30 @@ function (setup_test _dim _srcs _inputs) # # Add the test # - if (AMReX_OMP) + if (_COMBINE_MPI_OMP) + set(_parallel_cmd ${_cmd}) + if (AMReX_MPI) + set(_ntasks 2) + if (_NTASKS) + set(_ntasks ${_NTASKS}) + endif () + set(_parallel_cmd ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} ${_ntasks} + ${MPIEXEC_PREFLAGS} ${_cmd} ${MPIEXEC_POSTFLAGS}) + endif () + add_test( + NAME ${_test_name} + COMMAND ${_parallel_cmd} + WORKING_DIRECTORY ${_exe_dir} + ) + if (AMReX_OMP) + set(_nthreads 2) + if (_NTHREADS) + set(_nthreads ${_NTHREADS}) + endif () + set_tests_properties(${_test_name} PROPERTIES + ENVIRONMENT OMP_NUM_THREADS=${_nthreads}) + endif () + elseif (AMReX_OMP) add_test( NAME ${_test_name} COMMAND ${_cmd} @@ -89,7 +112,8 @@ function (setup_test _dim _srcs _inputs) elseif (AMReX_MPI) add_test( NAME ${_test_name} - COMMAND mpiexec -n 2 ${_cmd} + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + ${MPIEXEC_PREFLAGS} ${_cmd} ${MPIEXEC_POSTFLAGS} WORKING_DIRECTORY ${_exe_dir} ) else () @@ -102,6 +126,35 @@ function (setup_test _dim _srcs _inputs) endfunction () +# Register another inputs-file variant for an executable created by setup_test. +# Keep launcher behavior in one place so variants cannot silently become +# serial-only in MPI or use a different OpenMP thread count than the base test. +function (add_test_variant _test_name _target _inputs_filename) + set(_cmd $ ${_inputs_filename} ${ARGN}) + + if (AMReX_OMP) + add_test( + NAME ${_test_name} + COMMAND ${_cmd} + WORKING_DIRECTORY $ + ) + set_tests_properties(${_test_name} PROPERTIES ENVIRONMENT OMP_NUM_THREADS=2) + elseif (AMReX_MPI) + add_test( + NAME ${_test_name} + COMMAND ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + ${MPIEXEC_PREFLAGS} ${_cmd} ${MPIEXEC_POSTFLAGS} + WORKING_DIRECTORY $ + ) + else () + add_test( + NAME ${_test_name} + COMMAND ${_cmd} + WORKING_DIRECTORY $ + ) + endif () +endfunction () + if (AMReX_TEST_TYPE STREQUAL "Small") add_subdirectory("Base/TrackedVector") @@ -121,6 +174,7 @@ if (AMReX_TEST_TYPE STREQUAL "Small") if (AMReX_EB) add_subdirectory("EB/CNS") + add_subdirectory("EB/MarchingCubes") endif() if (AMReX_LINEAR_SOLVERS) diff --git a/Tests/EB/CNS/CMakeLists.txt b/Tests/EB/CNS/CMakeLists.txt index 972447373ea..c3b452d1487 100644 --- a/Tests/EB/CNS/CMakeLists.txt +++ b/Tests/EB/CNS/CMakeLists.txt @@ -35,7 +35,7 @@ unset(_hydro_sources) # ########################################################################################## set(_sod_sources Exec/Sod/cns_prob.F90 ${_sources}) -set(_input_files Exec/Sod/inputs-ci) +set(_input_files Exec/Sod/inputs-ci ../MarchingCubes/cube.stl) setup_test(3 _sod_sources _input_files HAS_FORTRAN_MODULES @@ -45,6 +45,16 @@ setup_test(3 _sod_sources _input_files target_link_options(Test_CNS_Sod_3d PRIVATE "$<$:-Mnomain>") +add_test_variant(CNS_Sod_MarchingCubes_3d Test_CNS_Sod_3d inputs-ci + max_step=1 + amr.plot_files_output=0 + amr.checkpoint_files_output=0 + eb2.geom_type=stl + eb2.stl_file=cube.stl + eb2.stl_scale=0.125 + "eb2.stl_center=0.7 0.5 0.5" + eb2.geometry_method=marching_cubes) + unset(_sod_sources) diff --git a/Tests/EB/MarchingCubes/CMakeLists.txt b/Tests/EB/MarchingCubes/CMakeLists.txt index 24280207166..c4c91efd0e5 100644 --- a/Tests/EB/MarchingCubes/CMakeLists.txt +++ b/Tests/EB/MarchingCubes/CMakeLists.txt @@ -4,9 +4,128 @@ endif () set(_sources main.cpp) -set(_input_files inputs cube.stl) +set(_input_files + inputs.cube.mc + inputs.cube.legacy + inputs.cubes.mc + inputs.cubes.legacy + inputs.bunny.mc + inputs.bunny.legacy + inputs.multilevel.mc + inputs.multilevel.rebuild.mc + inputs.nodal.mc + inputs.open-stl.mc + inputs.box.mc + inputs.box.nodal.mc + inputs.sphere.mc + inputs.sphere.legacy + inputs.sphere.rebuild.mc + inputs.parser.mc + inputs.api.mc + inputs.api.host.mc + inputs.api.multilevel.mc + inputs.box.periodic.mc + inputs.box.periodic.legacy + inputs.sphere.periodic.mc + inputs.box.crossing.mc + inputs.scale.mc + inputs.gyroid.mc + inputs.gyroid.nocover.mc + inputs.sphere.noncubic.mc + cube.stl + cubes.stl + stanford_bunny.stl + STANFORD_BUNNY_LICENSE.md + open_surface.stl +) -setup_test(3 _sources _input_files) +setup_test(3 _sources _input_files COMBINE_MPI_OMP NTASKS 2 NTHREADS 2) +# Input files stay in the configuration-independent runtime directory even +# when a multi-config generator puts the executable in a configuration subdir. +get_target_property(_mc_test_working_directory Test_EB_MarchingCubes_3d + RUNTIME_OUTPUT_DIRECTORY) + +set(_mc_test_labels EB MarchingCubes) +if (NOT AMReX_GPU_BACKEND STREQUAL NONE) + list(APPEND _mc_test_labels GPU) +endif() +set_tests_properties(EB_MarchingCubes_3d PROPERTIES + LABELS "${_mc_test_labels}") + +function(add_mc_test_variant _test_name _inputs_filename) + set(_cmd $ ${_inputs_filename} ${ARGN}) + if (AMReX_MPI) + set(_cmd ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + ${MPIEXEC_PREFLAGS} ${_cmd} ${MPIEXEC_POSTFLAGS}) + endif() + add_test( + NAME ${_test_name} + COMMAND ${_cmd} + WORKING_DIRECTORY ${_mc_test_working_directory} + ) + if (AMReX_OMP) + set_tests_properties(${_test_name} PROPERTIES ENVIRONMENT OMP_NUM_THREADS=2) + endif() + set_tests_properties(${_test_name} PROPERTIES LABELS "${_mc_test_labels}") +endfunction() + +add_mc_test_variant(EB_MarchingCubesMultilevel_3d inputs.multilevel.mc) +add_mc_test_variant(EB_MarchingCubesMultilevelRebuild_3d inputs.multilevel.rebuild.mc) +add_mc_test_variant(EB_MarchingCubesNodal_3d inputs.nodal.mc) +add_mc_test_variant(EB_MarchingCubesCubeLegacy_3d inputs.cube.legacy) +add_mc_test_variant(EB_MarchingCubesCubes_3d inputs.cubes.mc) +add_mc_test_variant(EB_MarchingCubesCubesLegacy_3d inputs.cubes.legacy) +add_mc_test_variant(EB_MarchingCubesBunny_3d inputs.bunny.mc) +add_mc_test_variant(EB_MarchingCubesBunnyLegacy_3d inputs.bunny.legacy) +add_mc_test_variant(EB_MarchingCubesBox_3d inputs.box.mc) +add_mc_test_variant(EB_MarchingCubesBoxNodal_3d inputs.box.nodal.mc) +add_mc_test_variant(EB_MarchingCubesSphere_3d inputs.sphere.mc) +add_mc_test_variant(EB_MarchingCubesSphereLegacy_3d inputs.sphere.legacy) +add_mc_test_variant(EB_MarchingCubesSphereRebuild_3d inputs.sphere.rebuild.mc) +add_mc_test_variant(EB_MarchingCubesParser_3d inputs.parser.mc) +add_mc_test_variant(EB_MarchingCubesApi_3d inputs.api.mc) +add_mc_test_variant(EB_MarchingCubesApiHost_3d inputs.api.host.mc) +add_mc_test_variant(EB_MarchingCubesApiMultilevel_3d inputs.api.multilevel.mc) +add_mc_test_variant(EB_MarchingCubesBoxPeriodic_3d inputs.box.periodic.mc) +add_mc_test_variant(EB_MarchingCubesBoxPeriodicLegacy_3d inputs.box.periodic.legacy) +add_mc_test_variant(EB_MarchingCubesSpherePeriodic_3d inputs.sphere.periodic.mc) +add_mc_test_variant(EB_MarchingCubesBoxCrossing_3d inputs.box.crossing.mc) +add_mc_test_variant(EB_MarchingCubesScaleInvariance_3d inputs.scale.mc) +add_mc_test_variant(EB_MarchingCubesGyroid_3d inputs.gyroid.mc) + +function(add_expected_failure _test_name _inputs_filename _expected_regex) + set(_test_prefix) + set(_test_suffix) + if (AMReX_MPI) + set(_test_prefix ${MPIEXEC_EXECUTABLE} ${MPIEXEC_NUMPROC_FLAG} 2 + ${MPIEXEC_PREFLAGS}) + set(_test_suffix ${MPIEXEC_POSTFLAGS}) + endif() + add_test( + NAME ${_test_name} + COMMAND ${CMAKE_COMMAND} + "-DTEST_PREFIX=${_test_prefix}" + "-DTEST_SUFFIX=${_test_suffix}" + -DTEST_EXECUTABLE=$ + -DTEST_INPUT=${_inputs_filename} + -DEXPECTED_REGEX=${_expected_regex} + -P ${CMAKE_CURRENT_LIST_DIR}/expect_failure.cmake + WORKING_DIRECTORY ${_mc_test_working_directory} + ) + if (AMReX_OMP) + set_tests_properties(${_test_name} PROPERTIES ENVIRONMENT OMP_NUM_THREADS=2) + endif() + set_tests_properties(${_test_name} PROPERTIES LABELS "${_mc_test_labels}") +endfunction() + +add_expected_failure(EB_MarchingCubesRejectOpenSTL_3d + inputs.open-stl.mc "STL is open, nonmanifold, or inconsistently oriented") +add_expected_failure(EB_MarchingCubesRejectNoCover_3d + inputs.gyroid.nocover.mc "set eb2.cover_multiple_cuts=1") +add_expected_failure(EB_MarchingCubesRejectNonCubic_3d + inputs.sphere.noncubic.mc "requires a 3D Cartesian grid with dx == dy == dz") unset(_sources) unset(_input_files) +unset(_mc_test_labels) +unset(_mc_test_working_directory) diff --git a/Tests/EB/MarchingCubes/GNUmakefile b/Tests/EB/MarchingCubes/GNUmakefile index e82a072d116..08572a4640a 100644 --- a/Tests/EB/MarchingCubes/GNUmakefile +++ b/Tests/EB/MarchingCubes/GNUmakefile @@ -4,7 +4,7 @@ USE_ASSERTION = TRUE USE_EB = TRUE -USE_MPI = FALSE +USE_MPI = TRUE USE_OMP = FALSE COMP = gnu @@ -26,4 +26,3 @@ Ppack += $(foreach dir, $(Pdirs), $(AMREX_HOME)/Src/$(dir)/Make.package) include $(Ppack) include $(AMREX_HOME)/Tools/GNUMake/Make.rules - diff --git a/Tests/EB/MarchingCubes/STANFORD_BUNNY_LICENSE.md b/Tests/EB/MarchingCubes/STANFORD_BUNNY_LICENSE.md new file mode 100644 index 00000000000..e4eda0a930b --- /dev/null +++ b/Tests/EB/MarchingCubes/STANFORD_BUNNY_LICENSE.md @@ -0,0 +1,20 @@ +# Stanford Bunny STL attribution + +`stanford_bunny.stl` is the Stanford Bunny model attributed to MakerBot. + +- Source: https://commons.wikimedia.org/wiki/File:Stanford_Bunny.stl +- Original source listed by Wikimedia Commons: + https://www.thingiverse.com/thing:88208/files +- License: Creative Commons Attribution 3.0 Unported + (https://creativecommons.org/licenses/by/3.0/) +- Retrieved from Wikimedia Commons on 2026-07-28. +- SHA-256: + `e1ff1293a49eb066de3c02cde6ccd260835e9da8544d43b1411f83f2d55c2eba` +- The checked-in STL is unmodified. The AMReX test applies a runtime scale and + translation without changing the source file. +- `stanford_bunny_mc.stl` and `stanford_bunny_legacy.stl` are the comparable + planar EB reconstructions produced from the MC and legacy factory data, + respectively. Both carry this attribution to the source model. +- `stanford_bunny_mc_triangles.stl` is the converged MC-owned triangle surface + after legacy-style nodal repair and carries this attribution to the source + model. diff --git a/Tests/EB/MarchingCubes/cubes.stl b/Tests/EB/MarchingCubes/cubes.stl new file mode 100644 index 00000000000..06231584202 Binary files /dev/null and b/Tests/EB/MarchingCubes/cubes.stl differ diff --git a/Tests/EB/MarchingCubes/expect_failure.cmake b/Tests/EB/MarchingCubes/expect_failure.cmake new file mode 100644 index 00000000000..ce0d565e049 --- /dev/null +++ b/Tests/EB/MarchingCubes/expect_failure.cmake @@ -0,0 +1,24 @@ +if (NOT DEFINED TEST_EXECUTABLE OR NOT DEFINED TEST_INPUT + OR NOT DEFINED EXPECTED_REGEX) + message(FATAL_ERROR "Expected-failure test is missing a required argument") +endif () + +set(test_command ${TEST_PREFIX} "${TEST_EXECUTABLE}" "${TEST_INPUT}" ${TEST_SUFFIX}) + +execute_process( + COMMAND ${test_command} + RESULT_VARIABLE test_result + OUTPUT_VARIABLE test_stdout + ERROR_VARIABLE test_stderr +) + +if (test_result STREQUAL "0") + message(FATAL_ERROR "Command unexpectedly succeeded: ${test_command}") +endif () + +set(test_output "${test_stdout}\n${test_stderr}") +if (NOT test_output MATCHES "${EXPECTED_REGEX}") + message(FATAL_ERROR + "Command failed for the wrong reason. Expected '${EXPECTED_REGEX}'.\n" + "Result: ${test_result}\nOutput:\n${test_output}") +endif () diff --git a/Tests/EB/MarchingCubes/inputs b/Tests/EB/MarchingCubes/inputs deleted file mode 100644 index 9cedd90c588..00000000000 --- a/Tests/EB/MarchingCubes/inputs +++ /dev/null @@ -1,6 +0,0 @@ - -nx = 64 -ny = 64 -nz = 64 - -max_grid_size = 32 diff --git a/Tests/EB/MarchingCubes/inputs.api.host.mc b/Tests/EB/MarchingCubes/inputs.api.host.mc new file mode 100644 index 00000000000..6ef32fd22a0 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.api.host.mc @@ -0,0 +1,14 @@ +# EB2::Build(makeShop(host-only functor), ...) through the C++ API with the +# marching-cubes generator selected by eb2.geometry_method. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 +api_build = host_sphere + +eb2.geometry_method = marching_cubes +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.7 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.api.mc b/Tests/EB/MarchingCubes/inputs.api.mc new file mode 100644 index 00000000000..9f8a0601e68 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.api.mc @@ -0,0 +1,14 @@ +# EB2::Build(makeShop(SphereIF), ...) through the C++ API with the +# marching-cubes generator selected by eb2.geometry_method. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 +api_build = sphere + +eb2.geometry_method = marching_cubes +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.7 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.api.multilevel.mc b/Tests/EB/MarchingCubes/inputs.api.multilevel.mc new file mode 100644 index 00000000000..1576364cdbe --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.api.multilevel.mc @@ -0,0 +1,21 @@ +# EB2::Build(makeShop(SphereIF), Vector, ...) builds every coarse +# level directly with the marching-cubes generator instead of coarsening. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 +api_build = sphere +api_all_levels = 1 +required_coarsening_level = 2 +max_coarsening_level = 2 + +eb2.geometry_method = marching_cubes +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.7 +eb2.sphere_has_fluid_inside = 0 +# The STL output must hold the finest level (same facet count as the +# single-level build in inputs.api.mc), not the last coarse level built. +eb2.mc_stl_file = api_multilevel_mc.stl +expected_mc_stl_facets = 13112 diff --git a/Tests/EB/MarchingCubes/inputs.box.crossing.mc b/Tests/EB/MarchingCubes/inputs.box.crossing.mc new file mode 100644 index 00000000000..44241007c00 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.box.crossing.mc @@ -0,0 +1,16 @@ +# Box that crosses the -y and +z domain faces with the default domain-face +# extension enabled: the geometry is extruded straight outward from the +# physical faces into the ghost region, and the in-domain volume is unchanged. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.extend_domain_face = 1 +eb2.geom_type = box +eb2.box_lo = -0.55 -1.5 -0.55 +eb2.box_hi = 0.55 0.47 1.5 +eb2.box_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.box.mc b/Tests/EB/MarchingCubes/inputs.box.mc new file mode 100644 index 00000000000..2c60647ea20 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.box.mc @@ -0,0 +1,14 @@ +# Analytic box implicit function through the marching-cubes generator. +# Same body as cube.stl, so the expected volume matches inputs.cube.mc. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = box +eb2.box_lo = -1.0 -1.0 -1.0 +eb2.box_hi = 1.0 1.0 1.0 +eb2.box_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.box.nodal.mc b/Tests/EB/MarchingCubes/inputs.box.nodal.mc new file mode 100644 index 00000000000..bba676beb9d --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.box.nodal.mc @@ -0,0 +1,24 @@ +# Box faces lying exactly on grid nodes: exercises the exact-zero (ON) node +# handling of the marching-cubes generator with an implicit function. +nx = 64 +ny = 64 +nz = 64 + +xmin = -1.0 +xmax = 1.0 +ymin = -1.0 +ymax = 1.0 +zmin = -1.0 +zmax = 1.0 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = box +eb2.box_lo = -0.5 -0.5 -0.5 +eb2.box_hi = 0.5 0.5 0.5 +eb2.box_has_fluid_inside = 0 +# The box faces lie on grid nodes, so cells owning a node-coincident EB patch +# legitimately have unit volume and a nonzero boundary area. +allow_full_volume_cut_cells = 1 diff --git a/Tests/EB/MarchingCubes/inputs.box.periodic.legacy b/Tests/EB/MarchingCubes/inputs.box.periodic.legacy new file mode 100644 index 00000000000..339b0aa9282 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.box.periodic.legacy @@ -0,0 +1,20 @@ +# Legacy counterpart of inputs.box.periodic.mc. +# Box that crosses the +x domain face, in a domain that is periodic in z, with +# the domain-face extension disabled. Exercises periodic ghost handling +# and eb2.extend_domain_face = 0; only the part of the box +# inside the domain counts toward the expected volume. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +geometry.is_periodic = 0 0 1 + +eb2.geometry_method = legacy +eb2.extend_domain_face = 0 +eb2.geom_type = box +eb2.box_lo = 0.4 -0.7 -0.7 +eb2.box_hi = 1.6 0.7 0.7 +eb2.box_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.box.periodic.mc b/Tests/EB/MarchingCubes/inputs.box.periodic.mc new file mode 100644 index 00000000000..19cf0879945 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.box.periodic.mc @@ -0,0 +1,20 @@ +# Box that crosses the +x domain face, in a domain that is periodic in z, with +# the domain-face extension disabled. Exercises periodic ghost handling +# and eb2.extend_domain_face = 0; only the part of the box +# inside the domain counts toward the expected volume. (The body stays clear +# of the periodic faces so the geometry is consistent with the periodicity.) +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +geometry.is_periodic = 0 0 1 + +eb2.geometry_method = marching_cubes +eb2.extend_domain_face = 0 +eb2.geom_type = box +eb2.box_lo = 0.4 -0.7 -0.7 +eb2.box_hi = 1.6 0.7 0.7 +eb2.box_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.bunny.legacy b/Tests/EB/MarchingCubes/inputs.bunny.legacy new file mode 100644 index 00000000000..2f436ccfdd3 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.bunny.legacy @@ -0,0 +1,17 @@ +nx = 128 +ny = 128 +nz = 128 +algorithm_tests = 0 + +max_grid_size = 32 + +custom_stl_test = 1 +eb2.geometry_method = legacy +eb2.stl_file = stanford_bunny.stl +eb2.stl_scale = 0.016 +eb2.stl_center = -0.48516235 -0.03025528 -0.94753541 + +eb2.small_volfrac = 1.e-6 +eb2.maxiter = 64 +eb2.cover_multiple_cuts = 1 +# eb2.eb_surface_stl_file = stanford_bunny_legacy.stl diff --git a/Tests/EB/MarchingCubes/inputs.bunny.mc b/Tests/EB/MarchingCubes/inputs.bunny.mc new file mode 100644 index 00000000000..463dc83d153 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.bunny.mc @@ -0,0 +1,17 @@ +nx = 128 +ny = 128 +nz = 128 +algorithm_tests = 0 + +max_grid_size = 32 + +custom_stl_test = 1 +eb2.geometry_method = marching_cubes +eb2.stl_file = stanford_bunny.stl +eb2.stl_scale = 0.016 +eb2.stl_center = -0.48516235 -0.03025528 -0.94753541 + +eb2.small_volfrac = 1.e-6 +eb2.maxiter = 64 +eb2.cover_multiple_cuts = 1 +# eb2.eb_surface_stl_file = stanford_bunny_mc.stl diff --git a/Tests/EB/MarchingCubes/inputs.cube.legacy b/Tests/EB/MarchingCubes/inputs.cube.legacy new file mode 100644 index 00000000000..6a278663454 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.cube.legacy @@ -0,0 +1,11 @@ +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +narrow_band_test = 1 +algorithm_tests = 0 + +eb2.geometry_method = legacy +eb2.stl_file = cube.stl +# eb2.eb_surface_stl_file = cube_legacy.stl diff --git a/Tests/EB/MarchingCubes/inputs.cube.mc b/Tests/EB/MarchingCubes/inputs.cube.mc new file mode 100644 index 00000000000..f4afbcc3d8f --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.cube.mc @@ -0,0 +1,11 @@ + +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +narrow_band_test = 1 + +eb2.geometry_method = marching_cubes +eb2.stl_file = cube.stl +# eb2.eb_surface_stl_file = cube_mc.stl diff --git a/Tests/EB/MarchingCubes/inputs.cubes.legacy b/Tests/EB/MarchingCubes/inputs.cubes.legacy new file mode 100644 index 00000000000..407256f76b9 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.cubes.legacy @@ -0,0 +1,19 @@ +nx = 256 +ny = 128 +nz = 64 +algorithm_tests = 0 + +max_grid_size = 32 + +custom_stl_test = 1 + +geometry.prob_lo = -2.4 -1.2 0.00 +geometry.prob_hi = 2.4 1.2 1.20 + +eb2.geometry_method = legacy +eb2.stl_file = cubes.stl +eb2.small_volfrac = 1.e-6 +eb2.stl_scale = 1e-3 +eb2.stl_center = 0.0 0.0 0.0 +eb2.cover_multiple_cuts = 1 +# eb2.eb_surface_stl_file = cubes_legacy.stl diff --git a/Tests/EB/MarchingCubes/inputs.cubes.mc b/Tests/EB/MarchingCubes/inputs.cubes.mc new file mode 100644 index 00000000000..3980bcd93f4 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.cubes.mc @@ -0,0 +1,19 @@ +nx = 256 +ny = 128 +nz = 64 +max_grid_size = 32 + +algorithm_tests = 0 +custom_stl_test = 1 + +geometry.prob_lo = -2.4 -1.2 0.00 +geometry.prob_hi = 2.4 1.2 1.20 + +eb2.geometry_method = marching_cubes +eb2.stl_file = cubes.stl +eb2.small_volfrac = 1.e-6 +eb2.stl_scale = 1e-3 +eb2.stl_center = 0.0 0.0 0.0 +eb2.cover_multiple_cuts = 1 +# eb2.mc_stl_file = cubes_mc_triangles.stl +# eb2.eb_surface_stl_file = cubes_mc.stl diff --git a/Tests/EB/MarchingCubes/inputs.gyroid.mc b/Tests/EB/MarchingCubes/inputs.gyroid.mc new file mode 100644 index 00000000000..a4413a89f22 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.gyroid.mc @@ -0,0 +1,15 @@ +# Coarsely resolved gyroid (parser implicit function) with the nodal repair +# enabled: its saddles produce ambiguous MC33 cells and multi-component fluid +# cells that the repair loop must cover. +nx = 32 +ny = 32 +nz = 32 + +max_grid_size = 16 +algorithm_tests = 0 +custom_stl_test = 1 + +eb2.geometry_method = marching_cubes +eb2.geom_type = parser +eb2.parser_function = "sin(6*x)*cos(6*y)+sin(6*y)*cos(6*z)+sin(6*z)*cos(6*x)-0.9" +eb2.cover_multiple_cuts = 1 diff --git a/Tests/EB/MarchingCubes/inputs.gyroid.nocover.mc b/Tests/EB/MarchingCubes/inputs.gyroid.nocover.mc new file mode 100644 index 00000000000..2f575620150 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.gyroid.nocover.mc @@ -0,0 +1,17 @@ +# Coarsely resolved gyroid (parser implicit function). Its saddles produce +# ambiguous MC33 cells whose fluid splits into more than one component, which +# the generator can only handle through the legacy-style nodal repair; with +# eb2.cover_multiple_cuts = 0 the build must abort with a message pointing at +# that parameter. +nx = 32 +ny = 32 +nz = 32 + +max_grid_size = 16 +algorithm_tests = 0 +custom_stl_test = 1 + +eb2.geometry_method = marching_cubes +eb2.geom_type = parser +eb2.parser_function = "sin(6*x)*cos(6*y)+sin(6*y)*cos(6*z)+sin(6*z)*cos(6*x)-0.9" +eb2.cover_multiple_cuts = 0 diff --git a/Tests/EB/MarchingCubes/inputs.multilevel.mc b/Tests/EB/MarchingCubes/inputs.multilevel.mc new file mode 100644 index 00000000000..793c497d0e1 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.multilevel.mc @@ -0,0 +1,11 @@ +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 + +algorithm_tests = 0 + +required_coarsening_level = 2 +max_coarsening_level = 2 +build_coarse_level_by_coarsening = 1 diff --git a/Tests/EB/MarchingCubes/inputs.multilevel.rebuild.mc b/Tests/EB/MarchingCubes/inputs.multilevel.rebuild.mc new file mode 100644 index 00000000000..a6683005555 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.multilevel.rebuild.mc @@ -0,0 +1,11 @@ +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 + +algorithm_tests = 0 + +required_coarsening_level = 2 +max_coarsening_level = 2 +build_coarse_level_by_coarsening = 0 diff --git a/Tests/EB/MarchingCubes/inputs.nodal.mc b/Tests/EB/MarchingCubes/inputs.nodal.mc new file mode 100644 index 00000000000..f94dd461bc1 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.nodal.mc @@ -0,0 +1,18 @@ +nx = 64 +ny = 64 +nz = 64 + +xmin = -1.0 +xmax = 1.0 +ymin = -1.0 +ymax = 1.0 +zmin = -1.0 +zmax = 1.0 + +max_grid_size = 32 +eb2.stl_scale = 0.5 + +algorithm_tests = 0 +# The STL faces lie on grid nodes, so cells owning a node-coincident EB patch +# legitimately have unit volume and a nonzero boundary area. +allow_full_volume_cut_cells = 1 diff --git a/Tests/EB/MarchingCubes/inputs.open-stl.mc b/Tests/EB/MarchingCubes/inputs.open-stl.mc new file mode 100644 index 00000000000..055b6643011 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.open-stl.mc @@ -0,0 +1,9 @@ +nx = 8 +ny = 8 +nz = 8 + +max_grid_size = 8 + +algorithm_tests = 0 + +eb2.stl_file = open_surface.stl diff --git a/Tests/EB/MarchingCubes/inputs.parser.mc b/Tests/EB/MarchingCubes/inputs.parser.mc new file mode 100644 index 00000000000..774f416128e --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.parser.mc @@ -0,0 +1,11 @@ +# Parser implicit function through the marching-cubes generator: a paraboloid. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = parser +eb2.parser_function = "0.36 - (x*x + y*y) - 0.5*z" diff --git a/Tests/EB/MarchingCubes/inputs.scale.mc b/Tests/EB/MarchingCubes/inputs.scale.mc new file mode 100644 index 00000000000..16a8d6bfb47 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.scale.mc @@ -0,0 +1,18 @@ +# Scale and decomposition invariance of the marching-cubes generator: the same two-sphere +# implicit function is built with the function multiplied by 1e-6, 1 and 1e6 +# and with two different eb2.max_grid_size values, and must give the same fluid +# volume and nodal repair. The main build is a +# plain sphere. +nx = 32 +ny = 32 +nz = 32 + +max_grid_size = 16 +algorithm_tests = 0 +scale_invariance_test = 1 + +eb2.geometry_method = marching_cubes +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.75 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.sphere.legacy b/Tests/EB/MarchingCubes/inputs.sphere.legacy new file mode 100644 index 00000000000..a0e208d6805 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.sphere.legacy @@ -0,0 +1,14 @@ +# Sphere implicit function through the legacy generator. The radius +# places the six poles exactly on grid nodes. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = legacy +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.75 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.sphere.mc b/Tests/EB/MarchingCubes/inputs.sphere.mc new file mode 100644 index 00000000000..7e21e279218 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.sphere.mc @@ -0,0 +1,15 @@ +# Sphere implicit function through the marching-cubes generator. The radius +# places the six poles exactly on grid nodes. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.75 +eb2.sphere_has_fluid_inside = 0 +eb2.mc_stl_file = sphere_mc.stl diff --git a/Tests/EB/MarchingCubes/inputs.sphere.noncubic.mc b/Tests/EB/MarchingCubes/inputs.sphere.noncubic.mc new file mode 100644 index 00000000000..cd091093acf --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.sphere.noncubic.mc @@ -0,0 +1,14 @@ +# Non-cubic cells (dz != dx): the marching-cubes generator must abort with a +# clear message instead of building an inconsistent geometry. +nx = 64 +ny = 64 +nz = 32 + +max_grid_size = 32 +algorithm_tests = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.75 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/inputs.sphere.periodic.mc b/Tests/EB/MarchingCubes/inputs.sphere.periodic.mc new file mode 100644 index 00000000000..8f2baf60da2 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.sphere.periodic.mc @@ -0,0 +1,19 @@ +# Sphere of radius 0.75 centered at x = 1.0 in the fully periodic domain +# [-1.2,1.2]^3, together with its periodic image at x = -1.4, so that the body +# is consistent with the periodicity while crossing the x faces. Ghost nodes +# sample the implicit function at their unwrapped coordinates, as in the +# legacy generator, which is only meaningful for a periodic body; the exported +# level set and cell flags must therefore agree on both sides of the seam. +# Fluid volume: 2.4^3 - 4/3*pi*0.75^3 = 12.0567. +nx = 32 +ny = 32 +nz = 32 + +max_grid_size = 16 +algorithm_tests = 0 + +geometry.is_periodic = 1 1 1 + +eb2.geometry_method = marching_cubes +eb2.geom_type = parser +eb2.parser_function = "0.5625 - min((x-1.0)*(x-1.0) + y*y + z*z, (x+1.4)*(x+1.4) + y*y + z*z)" diff --git a/Tests/EB/MarchingCubes/inputs.sphere.rebuild.mc b/Tests/EB/MarchingCubes/inputs.sphere.rebuild.mc new file mode 100644 index 00000000000..4ff0e5fb105 --- /dev/null +++ b/Tests/EB/MarchingCubes/inputs.sphere.rebuild.mc @@ -0,0 +1,18 @@ +# Coarse levels rebuilt from the implicit function (not coarsened) so that the +# marching-cubes generator runs through the GeometryShop rebuild path. +nx = 64 +ny = 64 +nz = 64 + +max_grid_size = 32 +algorithm_tests = 0 + +required_coarsening_level = 2 +max_coarsening_level = 2 +build_coarse_level_by_coarsening = 0 + +eb2.geometry_method = marching_cubes +eb2.geom_type = sphere +eb2.sphere_center = 0.0 0.0 0.0 +eb2.sphere_radius = 0.7 +eb2.sphere_has_fluid_inside = 0 diff --git a/Tests/EB/MarchingCubes/main.cpp b/Tests/EB/MarchingCubes/main.cpp index 3fc9b5756b3..ef799f51ff8 100644 --- a/Tests/EB/MarchingCubes/main.cpp +++ b/Tests/EB/MarchingCubes/main.cpp @@ -1,59 +1,1193 @@ #include #include +#include +#include +#include +#include #include -#include +#include + +#include +#include +#include +#include +#include +#include +#include using namespace amrex; +namespace { + +/** + * Host-only sphere implicit function. It is deliberately not GPU-callable so + * that the generic marching-cubes adapters exercise their RunOn::Cpu path + * (pinned staging of fillFab/getIntercept), the same way a user-defined + * host functor would. + */ +struct HostSphereIF +{ + Real m_radius; + RealArray m_center; + bool m_fluid_inside; + + Real operator() (RealArray const& p) const noexcept + { + Real d2 = 0.0_rt; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + d2 += (p[d]-m_center[d])*(p[d]-m_center[d]); + } + Real const v = std::sqrt(d2) - m_radius; + return m_fluid_inside ? v : -v; + } +}; + +/** + * Union of two overlapping spheres (fluid outside) with the implicit function + * multiplied by an arbitrary scale. The marching-cubes reconstruction, + * including its ambiguous-face decisions and the nodal repair, must not depend + * on the units of the implicit function. + */ +struct ScaledTwoSphereIF : amrex::GPUable +{ + Real m_scale; + + AMREX_GPU_HOST_DEVICE Real operator() (AMREX_D_DECL(Real x, Real y, Real z)) const noexcept + { + constexpr Real radius = 0.6_rt; + constexpr Real offset = 0.35_rt; + Real const d1 = std::sqrt((x-offset)*(x-offset) + y*y + z*z); + Real const d2 = std::sqrt((x+offset)*(x+offset) + y*y + z*z); + // Positive inside the body (EB2 convention), fluid outside. + return m_scale*amrex::max(radius-d1, radius-d2); + } + + Real operator() (RealArray const& p) const noexcept + { + return this->operator()(AMREX_D_DECL(p[0], p[1], p[2])); + } +}; + +/** + * Coarsely resolved high-frequency gyroid (fluid outside). Its saddles produce + * MC33 case-4/10 cells, so the interior (tunnel) test and the nodal repair of + * tunnel cells are exercised; both must be independent of the scale. + */ +struct ScaledGyroidIF : amrex::GPUable +{ + Real m_scale; + + AMREX_GPU_HOST_DEVICE Real operator() (AMREX_D_DECL(Real x, Real y, Real z)) const noexcept + { + constexpr Real k = 21.0_rt; + return m_scale*(std::sin(k*x)*std::cos(k*y) + std::sin(k*y)*std::cos(k*z) + + std::sin(k*z)*std::cos(k*x) - 0.1_rt); + } + + Real operator() (RealArray const& p) const noexcept + { + return this->operator()(AMREX_D_DECL(p[0], p[1], p[2])); + } +}; + +//! Implicit-function factories for check_scale_invariance(). These are named +//! types rather than function-local lambdas because nvcc rejects an extended +//! __device__ lambda inside a function template instantiated with a type local +//! to a function. +struct MakeTwoSphereShop +{ + auto operator() (Real scale) const { return EB2::makeShop(ScaledTwoSphereIF{{}, scale}); } +}; + +struct MakeGyroidShop +{ + auto operator() (Real scale) const { return EB2::makeShop(ScaledGyroidIF{{}, scale}); } +}; + +//! Build \p make_shop(scale) with the implicit function scaled by 1e-6, 1 and +//! 1e6 (and once with a different EB2 box size) and require the same fluid +//! volume and repaired-node count. +template +void check_scale_invariance (std::string const& name, int ncell, MakeShop const& make_shop) +{ + Box const domain(IntVect(0), IntVect(ncell - 1)); + Geometry const geom(domain, RealBox({-1.2_rt, -1.2_rt, -1.2_rt}, {1.2_rt, 1.2_rt, 1.2_rt}), + 0, {0, 0, 0}); + BoxArray ba(domain); + ba.maxSize(16); + DistributionMapping const dm(ba); + + // {scale of the implicit function, EB2 box size used for the build}: the + // reconstruction must not depend on either. + struct Variant { Real scale; int eb_max_grid_size; }; + Variant const variants[] = {{.scale = 1.e-6_rt, .eb_max_grid_size = 16}, + {.scale = 1.0_rt, .eb_max_grid_size = 16}, + {.scale = 1.e6_rt, .eb_max_grid_size = 16}, + {.scale = 1.0_rt, .eb_max_grid_size = 8}}; + int const saved_max_grid_size = EB2::max_grid_size; + Real reference_volume = -1.0_rt; + Long reference_zero_nodes = -1; + for (auto const& [scale, eb_max_grid_size] : variants) { + EB2::max_grid_size = eb_max_grid_size; + EB2::Build(make_shop(scale), geom, 0, 0, 4); + auto factory = makeEBFabFactory(geom, ba, dm, {1, 1, 1}, EBSupport::full); + Real const volume = factory->getVolFrac().sum() * geom.CellSize(0) * geom.CellSize(1) + * geom.CellSize(2); + auto const& levelset = factory->getLevelSet(); + Gpu::DeviceScalar zero_node_count(static_cast(0)); + Long* const zero_nodes = zero_node_count.dataPtr(); + for (MFIter mfi(levelset); mfi.isValid(); ++mfi) { + auto const phi = levelset.const_array(mfi); + ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (phi(i, j, k) == 0.0_rt) { + Gpu::Atomic::AddNoRet(zero_nodes, static_cast(1)); + } + }); + } + Long global_zero_nodes = zero_node_count.dataValue(); + ParallelAllReduce::Sum(global_zero_nodes, ParallelContext::CommunicatorSub()); + amrex::Print() << name << ": scale " << scale << ", eb2.max_grid_size " + << eb_max_grid_size << ": fluid volume " << volume + << ", zero level-set nodes " << global_zero_nodes << "\n"; + // The variants share every sign decision, so only the root finding + // (scale) and the summation order (decomposition) may differ at + // roundoff level. + if (reference_volume < 0.0_rt) { + reference_volume = volume; + reference_zero_nodes = global_zero_nodes; + } else { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + std::abs(volume - reference_volume) <= 1.e-8_rt * reference_volume, + "Marching-cubes fluid volume depends on the scale of the implicit " + "function or on the box decomposition"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_zero_nodes == reference_zero_nodes, + "Marching-cubes nodal repair depends on the scale of the implicit " + "function or on the box decomposition"); + } + EB2::IndexSpace::pop(); + } + EB2::max_grid_size = saved_max_grid_size; +} + +void validate_scale_invariance () +{ + check_scale_invariance("two spheres", 48, MakeTwoSphereShop{}); + // Tunnel-rich, badly under-resolved geometry: needs the nodal repair and + // takes a long cascade of repair passes. + ParmParse ppeb2("eb2"); + int cover_multiple_cuts = 0; + int maxiter = 32; + ppeb2.queryAdd("cover_multiple_cuts", cover_multiple_cuts); + ppeb2.queryAdd("maxiter", maxiter); + ppeb2.add("cover_multiple_cuts", 1); + ppeb2.add("maxiter", 200); + check_scale_invariance("gyroid", 24, MakeGyroidShop{}); + ppeb2.add("cover_multiple_cuts", cover_multiple_cuts); + ppeb2.add("maxiter", maxiter); +} + +//! A zeroed marching-cubes counter block for one FAB. +Gpu::Buffer make_fab_counters () +{ + Gpu::Buffer counters(MC::num_fab_counters); + std::fill_n(counters.hostData(), counters.size(), 0); + counters.copyToDeviceAsync(); + return counters; +} + +//! Level-wide sum of one counter after copying the block to the host. +int fab_counter (Gpu::Buffer& counters, MC::FabCounter which) +{ + counters.copyToHost(); + return counters.hostData()[which]; +} + +bool resolved_bottom_face_is_fluid_connected (GpuArray const& values) +{ + Box const cell_box(IntVect(0), IntVect(0)); + Box const node_box = amrex::surroundingNodes(cell_box); + FArrayBox sdf(node_box, 1); + sdf.setVal(-1.0_rt); + auto const phi = sdf.array(); + Box const cube_nodes = amrex::surroundingNodes(cell_box); + ParallelFor(cube_nodes, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const n = 4 * k + (j == 0 ? i : 3 - i); + phi(i, j, k) = values[n]; + }); + + Geometry geom(cell_box, RealBox({0.0_rt, 0.0_rt, 0.0_rt}, + {1.0_rt, 1.0_rt, 1.0_rt}), 0, {0, 0, 0}); + MC::MCFab result; + auto counters = make_fab_counters(); + MC::marching_cubes(geom, sdf, result, counters.data()); + AMREX_ALWAYS_ASSERT(fab_counter(counters, MC::counter_invalid_triangles) == 0); + + GpuArray host{}; + Gpu::dtoh_memcpy(host.data(), result.m_cell_data.dataPtr(), sizeof(host)); + Gpu::streamSynchronize(); + int const bit = 1 << 4; // MC33 face 5: the low-z face. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE((host[MC::face_decision_valid_mask] & bit) != 0, + "MC33 did not retain its resolved low-z face decision"); + return (host[MC::face_fluid_connected_mask] & bit) != 0; +} + +void validate_mc33_face_decisions () +{ + MC::Initialize(); + AMREX_ALWAYS_ASSERT(resolved_bottom_face_is_fluid_connected( + {2.0_rt, -1.0_rt, 2.0_rt, -1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt})); + AMREX_ALWAYS_ASSERT(!resolved_bottom_face_is_fluid_connected( + {1.0_rt, -2.0_rt, 1.0_rt, -2.0_rt, -1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt})); + AMREX_ALWAYS_ASSERT(resolved_bottom_face_is_fluid_connected( + {1.0_rt, -1.0_rt, 1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt, -1.0_rt})); +} + +void validate_mc33_vertex_indices () +{ + Box const cell_box(IntVect(0), IntVect(0)); + Box const node_box = amrex::surroundingNodes(cell_box); + Box const cube_nodes = amrex::surroundingNodes(cell_box); + Geometry geom(cell_box, RealBox({0.0_rt, 0.0_rt, 0.0_rt}, + {1.0_rt, 1.0_rt, 1.0_rt}), 0, {0, 0, 0}); + for (int mask = 0; mask < 256; ++mask) { + FArrayBox sdf(node_box, 1); + sdf.setVal(-1.0_rt); + auto const phi = sdf.array(); + ParallelFor(cube_nodes, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const n = 4 * k + (j == 0 ? i : 3 - i); + Real const magnitude = 1.0_rt + 0.125_rt * static_cast(n); + phi(i, j, k) = (mask & (1 << n)) != 0 ? magnitude : -magnitude; + }); + + MC::MCFab result; + auto counters = make_fab_counters(); + MC::marching_cubes(geom, sdf, result, counters.data()); + AMREX_ALWAYS_ASSERT(fab_counter(counters, MC::counter_invalid_triangles) == 0); + + GpuArray cell_data{}; + Gpu::dtoh_memcpy(cell_data.data(), result.m_cell_data.dataPtr(), sizeof(cell_data)); + + int const ntri = cell_data[MC::triangle_count]; + int const interior_count = cell_data[MC::interior_vertex_count]; + int const interior_offset = cell_data[MC::interior_vertex_offset]; + int const nvertices = static_cast(result.m_vertices.x.size()); + int const edge_vertex_count = nvertices - interior_count; + AMREX_ALWAYS_ASSERT(std::cmp_equal(ntri, result.m_triangles.v1.size())); + AMREX_ALWAYS_ASSERT(interior_count == 0 || interior_count == 1); + if (interior_count != 0) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + interior_offset == edge_vertex_count, + "MC33 interior vertex overlaps the Cartesian-edge vertex block"); + // The interior vertex is the average of the crossings on the 12 + // cube edges, so it must be their centroid and lie inside the cell. + Gpu::PinnedVector vx(nvertices), vy(nvertices), vz(nvertices); + Gpu::dtoh_memcpy(vx.data(), result.m_vertices.x.data(), nvertices * sizeof(Real)); + Gpu::dtoh_memcpy(vy.data(), result.m_vertices.y.data(), nvertices * sizeof(Real)); + Gpu::dtoh_memcpy(vz.data(), result.m_vertices.z.data(), nvertices * sizeof(Real)); + Gpu::streamSynchronize(); + Real cx = 0.0_rt, cy = 0.0_rt, cz = 0.0_rt; + for (int n = 0; n < edge_vertex_count; ++n) { + cx += vx[n]; cy += vy[n]; cz += vz[n]; + } + cx /= static_cast(edge_vertex_count); + cy /= static_cast(edge_vertex_count); + cz /= static_cast(edge_vertex_count); + Real const tol = 16.0_rt * std::numeric_limits::epsilon(); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + std::abs(vx[interior_offset] - cx) <= tol && + std::abs(vy[interior_offset] - cy) <= tol && + std::abs(vz[interior_offset] - cz) <= tol, + "MC33 interior vertex is not the centroid of the edge crossings"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + vx[interior_offset] >= 0.0_rt && vx[interior_offset] <= 1.0_rt && + vy[interior_offset] >= 0.0_rt && vy[interior_offset] <= 1.0_rt && + vz[interior_offset] >= 0.0_rt && vz[interior_offset] <= 1.0_rt, + "MC33 interior vertex lies outside the cell"); + } + + Gpu::PinnedVector v1(ntri); + Gpu::PinnedVector v2(ntri); + Gpu::PinnedVector v3(ntri); + if (ntri != 0) { + Gpu::dtoh_memcpy(v1.data(), result.m_triangles.v1.data(), ntri * sizeof(int)); + Gpu::dtoh_memcpy(v2.data(), result.m_triangles.v2.data(), ntri * sizeof(int)); + Gpu::dtoh_memcpy(v3.data(), result.m_triangles.v3.data(), ntri * sizeof(int)); + } + Gpu::streamSynchronize(); + for (int n = 0; n < ntri; ++n) { + AMREX_ALWAYS_ASSERT(v1[n] >= 0 && v1[n] < nvertices); + AMREX_ALWAYS_ASSERT(v2[n] >= 0 && v2[n] < nvertices); + AMREX_ALWAYS_ASSERT(v3[n] >= 0 && v3[n] < nvertices); + } + } +} + +/** + * A single-valued cell may hold exactly one face-connected fluid region. For + * every corner-sign mask and several magnitude patterns (chosen so that the + * MC33 face and interior tests take both branches, including the 4.1.2 / + * 10.1.2 tunnel tilings), run the extraction and the cell rejection rule and + * check that a cell is rejected exactly when its fluid corners form more than + * one group under cube-edge adjacency plus MC33-resolved connected ambiguous + * faces. In particular every tunnel (one connected surface patch joining two + * corner groups through the interior) must be rejected. + */ +void validate_cell_topology_rejection () +{ + Box const cell_box(IntVect(0), IntVect(0)); + Box const node_box = amrex::surroundingNodes(cell_box); + Geometry const geom(cell_box, RealBox({0.0_rt, 0.0_rt, 0.0_rt}, {1.0_rt, 1.0_rt, 1.0_rt}), + 0, {0, 0, 0}); + // MC33 face ids 1..6 in cube-corner order (test_face's A,B,C,D). + constexpr int face_corner[6][4] = {{0, 4, 5, 1}, {1, 5, 6, 2}, {2, 6, 7, 3}, + {3, 7, 4, 0}, {0, 3, 2, 1}, {4, 7, 6, 5}}; + constexpr int edge_lo[12] = {0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3}; + constexpr int edge_hi[12] = {1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7}; + + unsigned int seed = 12345U; + auto next_magnitude = [&seed] () { + seed = seed * 1664525U + 1013904223U; // LCG, deterministic + return 0.05_rt + 3.0_rt * static_cast(seed >> 8) / 16777216.0_rt; + }; + + // Every mask a few times, plus many draws of the MC33 case-4 masks (two + // body-diagonal fluid corners), whose interior test selects the 4.1.2 + // tunnel for roughly one draw in a hundred. + constexpr int case4_masks[8] = {65, 130, 20, 40, 190, 125, 235, 215}; + auto configurations = [&] (int trial, int slot) { + return trial < 24 ? slot + 1 : case4_masks[slot % 8]; + }; + Long checked = 0, rejected_cells = 0, tunnels = 0, mismatches = 0; + for (int trial = 0; trial < 24 + 40; ++trial) { + int const nslots = trial < 24 ? 254 : 8 * 8; + for (int slot = 0; slot < nslots; ++slot) { + int const mask = configurations(trial, slot); + GpuArray values{}; + for (int n = 0; n < 8; ++n) { + Real const magnitude = next_magnitude(); + values[n] = (mask & (1 << n)) != 0 ? magnitude : -magnitude; + } + FArrayBox sdf(node_box, 1); + auto const phi = sdf.array(); + ParallelFor(node_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const n = 4 * k + (j == 0 ? i : 3 - i); + phi(i, j, k) = values[n]; + }); + + MC::MCFab result; + auto counters = make_fab_counters(); + MC::marching_cubes(geom, sdf, result, counters.data()); + FArrayBox vfrac(cell_box, 1); + vfrac.setVal(0.5_rt); // no small-cell or sentinel rejection + IArrayBox rejected(cell_box, 1); + rejected.setVal(0); + MC::mark_cells_for_cleanup(cell_box, result, sdf, vfrac, 0.0_rt, rejected, + counters.data()); + Gpu::streamSynchronize(); + + GpuArray cell_data{}; + Gpu::dtoh_memcpy(cell_data.data(), result.m_cell_data.dataPtr(), sizeof(cell_data)); + int rejected_flag = 0; + Gpu::dtoh_memcpy(&rejected_flag, rejected.dataPtr(), sizeof(int)); + AMREX_ALWAYS_ASSERT(fab_counter(counters, MC::counter_invalid_triangles) == 0); + + // Reference: fluid corner groups joined along edges and across + // ambiguous faces whose stored MC33 decision is "connected". + bool fluid[8]; + for (int n = 0; n < 8; ++n) { fluid[n] = values[n] > 0.0_rt; } + int parent[8]; + for (int n = 0; n < 8; ++n) { parent[n] = n; } + auto root = [&] (int n) { while (parent[n] != n) { n = parent[n]; } return n; }; + auto join = [&] (int a, int b) { int ra = root(a), rb = root(b); if (ra != rb) { parent[rb] = ra; } }; + for (int e = 0; e < 12; ++e) { + if (fluid[edge_lo[e]] && fluid[edge_hi[e]]) { join(edge_lo[e], edge_hi[e]); } + } + for (int f = 0; f < 6; ++f) { + int const c0 = face_corner[f][0], c1 = face_corner[f][1]; + int const c2 = face_corner[f][2], c3 = face_corner[f][3]; + if (!(fluid[c0] == fluid[c2] && fluid[c1] == fluid[c3] && fluid[c0] != fluid[c1])) { + continue; + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + (cell_data[MC::face_decision_valid_mask] & (1 << f)) != 0, + "An ambiguous face has no stored MC33 decision"); + if ((cell_data[MC::face_fluid_connected_mask] & (1 << f)) != 0) { + if (fluid[c0]) { join(c0, c2); } else { join(c1, c3); } + } + } + int groups = 0; + for (int n = 0; n < 8; ++n) { groups += fluid[n] && root(n) == n; } + + // Tunnel: several corner groups but a single connected surface patch. + int const ntri = cell_data[MC::triangle_count]; + bool single_patch = ntri > 0; + if (groups > 1 && ntri > 0) { + Gpu::PinnedVector v1(ntri), v2(ntri), v3(ntri); + Gpu::dtoh_memcpy(v1.data(), result.m_triangles.v1.data(), ntri * sizeof(int)); + Gpu::dtoh_memcpy(v2.data(), result.m_triangles.v2.data(), ntri * sizeof(int)); + Gpu::dtoh_memcpy(v3.data(), result.m_triangles.v3.data(), ntri * sizeof(int)); + Gpu::streamSynchronize(); + std::vector tparent(ntri); + for (int t = 0; t < ntri; ++t) { tparent[t] = t; } + auto troot = [&] (int t) { while (tparent[t] != t) { t = tparent[t]; } return t; }; + for (int a = 0; a < ntri; ++a) { + for (int b = a + 1; b < ntri; ++b) { + int const va[3] = {v1[a], v2[a], v3[a]}; + int const vb[3] = {v1[b], v2[b], v3[b]}; + bool share = false; + for (int const x : va) { for (int const y : vb) { share = share || x == y; } } + if (share) { int ra = troot(a), rb = troot(b); if (ra != rb) { tparent[rb] = ra; } } + } + } + int patches = 0; + for (int t = 0; t < ntri; ++t) { patches += troot(t) == t; } + single_patch = patches == 1; + if (single_patch) { ++tunnels; } + } + + bool const expect_rejected = groups > 1; + ++checked; + rejected_cells += rejected_flag != 0; + if ((rejected_flag != 0) != expect_rejected) { + ++mismatches; + amrex::Print() << " mask " << mask << " trial " << trial << ": groups=" << groups + << " rejected=" << rejected_flag << " triangles=" << ntri << "\n"; + } + if (groups > 1 && single_patch) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rejected_flag != 0, + "An MC33 tunnel cell was not rejected"); + } + } + } + amrex::Print() << "Cell topology rule: " << checked << " configurations, " << rejected_cells + << " rejected, " << tunnels << " tunnel tilings (all rejected)\n"; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(mismatches == 0, + "Cell rejection disagrees with the face-connected corner groups"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(tunnels > 0, + "The topology test did not exercise any tunnel tiling"); +} + +void validate_exact_ambiguous_face_fraction () +{ + Box const cell_box(IntVect(0), IntVect(0)); + Box const node_box = amrex::surroundingNodes(cell_box); + FArrayBox sdf(node_box, 1); + sdf.setVal(1.0_rt); + auto const phi = sdf.array(); + ParallelFor(amrex::surroundingNodes(cell_box), + [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (k == 0) { + phi(i, j, k) = (i == j) ? 1.0_rt : -1.0_rt; + } + }); + + MC::MCFab mc_fab; + mc_fab.m_cell_data.resize(cell_box, MC::num_cell_data_components); + mc_fab.m_cell_data.setVal(0); + auto const cell_data = mc_fab.m_cell_data.array(); + ParallelFor(cell_box, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + cell_data(i, j, k, MC::face_decision_valid_mask) = 1 << 4; + }); + + mc_fab.defineEdgeIntersections(sdf.box()); + auto const exact_x = mc_fab.m_edge_intersections[0].array(); + auto const exact_y = mc_fab.m_edge_intersections[1].array(); + ParallelFor(cell_box, [=] AMREX_GPU_DEVICE(int, int, int) noexcept { + exact_x(0, 0, 0) = 0.2_rt; + exact_y(1, 0, 0) = 0.3_rt; + exact_x(0, 1, 0) = 0.6_rt; + exact_y(0, 0, 0) = 0.9_rt; + }); + + FArrayBox apx(amrex::surroundingNodes(cell_box, 0), 1); + FArrayBox apy(amrex::surroundingNodes(cell_box, 1), 1); + FArrayBox apz(amrex::surroundingNodes(cell_box, 2), 1); + FArrayBox fcx(amrex::surroundingNodes(cell_box, 0), 2); + FArrayBox fcy(amrex::surroundingNodes(cell_box, 1), 2); + FArrayBox fcz(amrex::surroundingNodes(cell_box, 2), 2); + IArrayBox rejected_x(amrex::surroundingNodes(cell_box, 0), 1); + IArrayBox rejected_y(amrex::surroundingNodes(cell_box, 1), 1); + IArrayBox rejected_z(amrex::surroundingNodes(cell_box, 2), 1); + rejected_x.setVal(0); + rejected_y.setVal(0); + rejected_z.setVal(0); + auto counters = make_fab_counters(); + MC::build_face_fractions(cell_box, mc_fab, sdf, apx, apy, apz, fcx, fcy, fcz, rejected_x, + rejected_y, rejected_z, counters.data()); + int const errors = fab_counter(counters, MC::counter_face_decision_errors) + + fab_counter(counters, MC::counter_degenerate_faces); + + Gpu::Buffer area_buffer(1); + Box const sample_box = amrex::makeSingleCellBox(IntVect(0), apz.box().ixType()); + apz.copyToMem(sample_box, 0, 1, area_buffer.data()); + area_buffer.copyToHost(); + Real const area = area_buffer.hostData()[0]; + AMREX_ALWAYS_ASSERT(errors == 0); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + std::abs(area - 0.23_rt) < 64.0_rt * std::numeric_limits::epsilon(), + "Ambiguous MC face fraction did not use exact STL edge crossings"); +} + +void validate_nodal_stl_face_levelset () +{ + Box const domain(IntVect(0), IntVect(7)); + Geometry const geom( + domain, + RealBox({-0.4875_rt, -0.3625_rt, -0.075_rt}, + {-0.3375_rt, -0.2125_rt, 0.075_rt}), + 0, {0, 0, 0}); + BoxArray const nodal_grid(amrex::surroundingNodes(domain)); + DistributionMapping const dm(nodal_grid); + MultiFab levelset(nodal_grid, dm, 1, 1); + + STLtools stl; + stl.setUseMarchingCubes(true); + stl.read_stl_file("cube.stl", 0.05_rt, {-0.4_rt, -0.275_rt, 0.0_rt}, 0); + stl.fillMarchingCubesLevelSet(levelset, levelset.nGrowVect(), geom); + + Gpu::DeviceScalar error_count(static_cast(0)); + Long* const errors = error_count.dataPtr(); + for (MFIter mfi(levelset); mfi.isValid(); ++mfi) { + auto const phi = levelset.const_array(mfi); + Box const face_nodes(IntVect(2, 2, 2), IntVect(2, 7, 6)); + ParallelFor(face_nodes, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (phi(i, j, k) != 0.0_rt) { + Gpu::Atomic::AddNoRet(errors, static_cast(1)); + } + }); + } + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + error_count.dataValue() == 0, + "STL samples on a grid-aligned face were not consistently classified as on-surface"); +} + +void validate_narrow_band_levelset () +{ + Box const domain(IntVect(0), IntVect(15)); + Geometry const geom(domain, RealBox({-1.0_rt, -1.0_rt, -1.0_rt}, {1.0_rt, 1.0_rt, 1.0_rt}), 0, + {0, 0, 0}); + BoxArray ba(domain); + ba.maxSize(8); + DistributionMapping dm(ba); + BoxArray const nba = amrex::convert(ba, IntVect::TheNodeVector()); + MultiFab exact(nba, dm, 1, 1); + MultiFab narrow(nba, dm, 1, 1); + + STLtools stl; + stl.setUseMarchingCubes(true); + stl.read_stl_file("cube.stl", 0.43_rt, {0.037_rt, -0.021_rt, 0.015_rt}, 0); + stl.fillSignedDistance(exact, exact.nGrowVect(), geom); + exact.OverrideSync(geom.periodicity()); + exact.FillBoundary(geom.periodicity()); + stl.fillMarchingCubesLevelSet(narrow, narrow.nGrowVect(), geom); + + Gpu::Buffer errors(5); + std::fill_n(errors.hostData(), errors.size(), static_cast(0)); + errors.copyToDeviceAsync(); + Long* const error = errors.data(); + for (MFIter mfi(exact); mfi.isValid(); ++mfi) { + auto const full = exact.const_array(mfi); + auto const band = narrow.const_array(mfi); + MC::MCFab mc_fab; + mc_fab.defineEdgeIntersections(narrow[mfi].box()); + stl.fillMarchingCubesEdgeIntersections(narrow[mfi], geom, mc_fab); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + auto const crossing = mc_fab.m_edge_intersections[idim].const_array(); + ParallelFor(Box{crossing}, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const di = idim == 0; + int const dj = idim == 1; + int const dk = idim == 2; + Real const lo = band(i, j, k); + Real const hi = band(i + di, j + dj, k + dk); + Real const refined = crossing(i, j, k); + // Edges without an exact STL crossing keep the sentinel; a + // sign-changing edge inside the band must always have one. + if ((lo > 0.0_rt) != (hi > 0.0_rt)) { + if (refined == MC::invalid_edge_intersection) { + Gpu::Atomic::AddNoRet(error + 4, static_cast(1)); + return; + } + Gpu::Atomic::AddNoRet(error + 2, static_cast(1)); + Real const linear = lo / (lo - hi); + if (std::abs(refined - linear) > + 64.0_rt * std::numeric_limits::epsilon()) { + Gpu::Atomic::AddNoRet(error + 3, static_cast(1)); + } + } + }); + } + ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if ((full(i, j, k) > 0.0_rt) != (band(i, j, k) > 0.0_rt)) { + Gpu::Atomic::AddNoRet(error, static_cast(1)); + } + }); + + Box const cbx = amrex::enclosedCells(mfi.validbox()); + ParallelFor(cbx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + bool has_fluid = false; + bool has_covered = false; + for (int kk = 0; kk <= 1; ++kk) { + for (int jj = 0; jj <= 1; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + bool const fluid = full(i + ii, j + jj, k + kk) > 0.0_rt; + has_fluid = has_fluid || fluid; + has_covered = has_covered || !fluid; + } + } + } + if (has_fluid && has_covered) { + for (int kk = 0; kk <= 1; ++kk) { + for (int jj = 0; jj <= 1; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + Real const a = full(i + ii, j + jj, k + kk); + Real const b = band(i + ii, j + jj, k + kk); + Real const tolerance = 16.0_rt * std::numeric_limits::epsilon() * + amrex::max(1.0_rt, std::abs(a)); + if (std::abs(a - b) > tolerance) { + Gpu::Atomic::AddNoRet(error + 1, static_cast(1)); + } + } + } + } + } + }); + } + errors.copyToHost(); + GpuArray global_errors{errors.hostData()[0], errors.hostData()[1], + errors.hostData()[2], errors.hostData()[3], + errors.hostData()[4]}; + ParallelAllReduce::Sum(global_errors.data(), static_cast(global_errors.size()), + ParallelContext::CommunicatorSub()); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_errors[0] == 0, "Narrow-band MC field changed the global STL sign classification"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(global_errors[1] == 0, + "Narrow-band MC field is not exact on a " + "cut-cell interpolation stencil"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(global_errors[2] > 0, + "Exact STL refinement found no sign-changing Cartesian edges"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_errors[3] > 0, "Exact STL refinement did not move any sampled-SDF edge crossing"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_errors[4] == 0, "A sign-changing edge in the band has no exact STL crossing"); +} + +char const* execution_backend () noexcept +{ +#if defined(AMREX_USE_CUDA) + return "CUDA"; +#elif defined(AMREX_USE_HIP) + return "HIP"; +#elif defined(AMREX_USE_SYCL) + return "SYCL"; +#else + return "CPU"; +#endif +} + +} // namespace + void main_main () { +#ifdef AMREX_USE_GPU + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + Gpu::inLaunchRegion(), "The marching-cubes GPU test must execute in a GPU launch region"); +#endif + bool algorithm_tests = true; + ParmParse pp; + pp.query("algorithm_tests", algorithm_tests); + if (algorithm_tests) { + validate_mc33_face_decisions(); + validate_mc33_vertex_indices(); + validate_cell_topology_rejection(); + validate_exact_ambiguous_face_fraction(); + validate_nodal_stl_face_levelset(); + } + bool scale_invariance_test = false; + pp.query("scale_invariance_test", scale_invariance_test); + int nx = 64; int ny = 64; int nz = 64; int max_grid_size = 32; + int required_coarsening_level = 0; + int max_coarsening_level = 0; + bool build_coarse_level_by_coarsening = true; Real xmin = -1.2; - Real xmax = 1.2; + Real xmax = 1.2; Real ymin = -1.2; - Real ymax = 1.2; + Real ymax = 1.2; Real zmin = -1.2; - Real zmax = 1.2; + Real zmax = 1.2; + bool custom_stl_test = false; + bool narrow_band_test = false; + std::string eb_method("marching_cubes"); + std::string geom_type("stl"); + std::string api_build; + bool api_all_levels = false; + bool allow_full_volume_cut_cells = false; + std::vector is_periodic{0, 0, 0}; + std::string stl_file("cube.stl"); + Real stl_scale = 1.0; + std::vector stl_center{0.0, 0.0, 0.0}; + Real body_volume = -1.0_rt; // analytic solid volume, < 0 when unknown + bool fluid_inside = false; { - ParmParse pp; pp.query("nx", nx); pp.query("ny", ny); pp.query("nz", nz); pp.query("max_grid_size", max_grid_size); + pp.query("xmin", xmin); + pp.query("xmax", xmax); + pp.query("ymin", ymin); + pp.query("ymax", ymax); + pp.query("zmin", zmin); + pp.query("zmax", zmax); + pp.query("required_coarsening_level", required_coarsening_level); + pp.query("max_coarsening_level", max_coarsening_level); + pp.query("build_coarse_level_by_coarsening", build_coarse_level_by_coarsening); + pp.query("custom_stl_test", custom_stl_test); + pp.query("narrow_band_test", narrow_band_test); + // Build through the C++ API with an explicit GeometryShop instead of + // the ParmParse-driven EB2::Build: "sphere" (GPU-callable SphereIF) or + // "host_sphere" (host-only functor). + pp.query("api_build", api_build); + // With api_build, also build every coarse level directly from the + // GeometryShop (EB2::Build overload taking a Vector) instead + // of coarsening the fine level. + pp.query("api_all_levels", api_all_levels); + std::vector prob_lo{xmin, ymin, zmin}; + std::vector prob_hi{xmax, ymax, zmax}; + ParmParse ppgeom("geometry"); + if (ppgeom.queryarr("prob_lo", prob_lo)) { + xmin = prob_lo[0]; + ymin = prob_lo[1]; + zmin = prob_lo[2]; + } + if (ppgeom.queryarr("prob_hi", prob_hi)) { + xmax = prob_hi[0]; + ymax = prob_hi[1]; + zmax = prob_hi[2]; + } + ppgeom.queryarr("is_periodic", is_periodic); + pp.query("allow_full_volume_cut_cells", allow_full_volume_cut_cells); ParmParse ppeb2("eb2"); - std::string geom_type("stl"); - ppeb2.add("geom_type", geom_type); - ppeb2.add("test_marching_cubes", 1); - } + ppeb2.queryAdd("geom_type", geom_type); + ppeb2.queryAdd("geometry_method", eb_method); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + eb_method == "legacy" || eb_method == "marching_cubes", + "eb2.geometry_method must be legacy or marching_cubes"); - { - std::string stl_file("cube.stl"); - Real stl_scale = 1.0; - std::vector stl_center{0.0, 0.0, 0.0}; + // Analytic solid volume for the geometry types the test knows. + if (!api_build.empty() || geom_type == "sphere") { + Real radius = 0.5_rt; + ppeb2.queryAdd("sphere_radius", radius); + std::vector center{0.0_rt, 0.0_rt, 0.0_rt}; + ppeb2.queryAdd("sphere_center", center); + int has_fluid_inside = 0; + ppeb2.queryAdd("sphere_has_fluid_inside", has_fluid_inside); + fluid_inside = has_fluid_inside != 0; + body_volume = 4.0_rt/3.0_rt*std::acos(-1.0_rt)*radius*radius*radius; + } else if (geom_type == "stl") { + ppeb2.queryAdd("stl_file", stl_file); + ppeb2.queryAdd("stl_scale", stl_scale); + ppeb2.queryAdd("stl_center", stl_center); + // cube.stl is the [-1,1]^3 cube; other STL files use custom_stl_test. + body_volume = 8.0_rt*stl_scale*stl_scale*stl_scale; + } else if (geom_type == "box") { + std::vector lo, hi; + ppeb2.getarr("box_lo", lo); + ppeb2.getarr("box_hi", hi); + int has_fluid_inside = 0; + ppeb2.queryAdd("box_has_fluid_inside", has_fluid_inside); + fluid_inside = has_fluid_inside != 0; + // Only the part of the box inside the domain contributes. + Real const clipped[3] = { + amrex::max(0.0_rt, amrex::min(hi[0], xmax) - amrex::max(lo[0], xmin)), + amrex::max(0.0_rt, amrex::min(hi[1], ymax) - amrex::max(lo[1], ymin)), + amrex::max(0.0_rt, amrex::min(hi[2], zmax) - amrex::max(lo[2], zmin))}; + body_volume = clipped[0]*clipped[1]*clipped[2]; + } else if (geom_type == "cylinder") { + Real radius = 0.5_rt; + Real height = -1.0_rt; + ppeb2.get("cylinder_radius", radius); + ppeb2.queryAdd("cylinder_height", height); + int has_fluid_inside = 0; + ppeb2.queryAdd("cylinder_has_fluid_inside", has_fluid_inside); + fluid_inside = has_fluid_inside != 0; + if (height > 0.0_rt) { + body_volume = std::acos(-1.0_rt)*radius*radius*height; + } + } + if (body_volume < 0.0_rt) { + // parser, torus, plane, ...: only range and consistency checks. + custom_stl_test = true; + } + } - ParmParse pp("eb2"); - pp.add("stl_file", stl_file); - pp.add("stl_scale", stl_scale); - pp.addarr("stl_center", stl_center); + if (narrow_band_test) { + validate_narrow_band_levelset(); + } + if (scale_invariance_test) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(eb_method == "marching_cubes", + "scale_invariance_test requires marching_cubes"); + validate_scale_invariance(); } - Geometry geom(Box(IntVect(0),IntVect(AMREX_D_DECL(nx,ny,nz))), - RealBox({AMREX_D_DECL(xmin,ymin,zmin)}, - {AMREX_D_DECL(xmax,ymax,zmax)}), - 0, {AMREX_D_DECL(0,0,0)}); + Geometry geom(Box(IntVect(0), IntVect(AMREX_D_DECL(nx - 1, ny - 1, nz - 1))), + RealBox({AMREX_D_DECL(xmin, ymin, zmin)}, {AMREX_D_DECL(xmax, ymax, zmax)}), 0, + {AMREX_D_DECL(is_periodic[0], is_periodic[1], is_periodic[2])}); BoxArray ba(geom.Domain()); ba.maxSize(max_grid_size); DistributionMapping dm(ba); double t0 = amrex::second(); - EB2::BuildMultiValuedMultiCut(geom,0,0); + if (api_build.empty()) { + EB2::Build(geom, required_coarsening_level, max_coarsening_level, 4, + build_coarse_level_by_coarsening); + } else { + Real radius = 0.5_rt; + std::vector center{0.0_rt, 0.0_rt, 0.0_rt}; + ParmParse ppeb2("eb2"); + ppeb2.queryAdd("sphere_radius", radius); + ppeb2.queryAdd("sphere_center", center); + RealArray const c{AMREX_D_DECL(center[0], center[1], center[2])}; + Vector all_geoms; + if (api_all_levels) { + all_geoms.push_back(geom); + for (int ilev = 1; ilev <= max_coarsening_level; ++ilev) { + all_geoms.push_back(amrex::coarsen(all_geoms.back(), 2)); + } + } + auto build_with = [&] (auto const& gshop) { + if (api_all_levels) { + EB2::Build(gshop, all_geoms, 4); + } else { + EB2::Build(gshop, geom, required_coarsening_level, max_coarsening_level, 4, + build_coarse_level_by_coarsening); + } + }; + if (api_build == "sphere") { + build_with(EB2::makeShop(EB2::SphereIF(radius, c, fluid_inside))); + } else if (api_build == "host_sphere") { + static_assert(!EB2::IsGPUable::value, + "HostSphereIF must exercise the host-only adapter path"); + build_with(EB2::makeShop(HostSphereIF{.m_radius = radius, .m_center = c, + .m_fluid_inside = fluid_inside})); + } else { + amrex::Abort("api_build must be sphere or host_sphere"); + } + } double t1 = amrex::second(); - amrex::Print() << "Build time: " << t1-t0 << "\n"; + amrex::Print() << "EB method: " << eb_method << ", build time: " << t1 - t0 << "\n"; + auto factory = makeEBFabFactory(geom, ba, dm, {1, 1, 1}, EBSupport::full); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(factory->maxCoarseningLevel() >= required_coarsening_level, + "Marching-cubes EB did not build every required coarse level"); + std::string eb_surface_stl_file; + if (ParmParse("eb2").query("eb_surface_stl_file", eb_surface_stl_file)) { + WriteEBSurfaceSTL(ba, dm, geom, factory.get(), eb_surface_stl_file); + } + std::string mc_stl_file; + if (ParmParse("eb2").query("mc_stl_file", mc_stl_file) && eb_method == "marching_cubes") { + // The file holds the converged triangulation of the finest level. + Long facets = 0; + if (ParallelDescriptor::IOProcessor()) { + std::ifstream ifs(mc_stl_file); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(ifs.good(), "eb2.mc_stl_file was not written"); + std::string line; + while (std::getline(ifs, line)) { + if (line.find("facet normal") != std::string::npos) { ++facets; } + } + } + ParallelDescriptor::Bcast(&facets, 1); + amrex::Print() << "Marching-cubes STL facets: " << facets << "\n"; + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(facets > 0, "eb2.mc_stl_file contains no facets"); + Long expected_facets = -1; + pp.query("expected_mc_stl_facets", expected_facets); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(expected_facets < 0 || facets == expected_facets, + "eb2.mc_stl_file facet count differs from the " + "expected finest-level count"); + } + auto const& volfrac = factory->getVolFrac(); + Real const fluid_volume = + volfrac.sum() * geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2); + Real const domain_volume = (xmax - xmin) * (ymax - ymin) * (zmax - zmin); + Real const expected_volume = fluid_inside ? body_volume : domain_volume - body_volume; + Real const error = std::abs(fluid_volume - expected_volume); + Real const max_dx = + amrex::max(geom.CellSize(0), amrex::max(geom.CellSize(1), geom.CellSize(2))); + Real const roundoff_tolerance = 64.0_rt * std::numeric_limits::epsilon(); + if (custom_stl_test) { + amrex::Print() << "Fluid volume: " << fluid_volume << "\n"; + } else { + amrex::Print() << "Fluid volume: " << fluid_volume << ", expected: " << expected_volume + << ", error: " << error << "\n"; + } + AMREX_ALWAYS_ASSERT(volfrac.min(0) >= -roundoff_tolerance); + AMREX_ALWAYS_ASSERT(volfrac.max(0) <= 1.0_rt + roundoff_tolerance); + if (custom_stl_test) { + AMREX_ALWAYS_ASSERT(fluid_volume > 0.0_rt); + AMREX_ALWAYS_ASSERT(fluid_volume < domain_volume); + } + + auto volume_centroid = factory->getCentroid().ToMultiFab(0.0_rt, 0.0_rt); + auto boundary_centroid = factory->getBndryCent().ToMultiFab(0.0_rt, 0.0_rt); + auto boundary_normal = factory->getBndryNormal().ToMultiFab(0.0_rt, 0.0_rt); + auto bndry_area = factory->getBndryArea().ToMultiFab(0.0_rt, 0.0_rt); + AMREX_ALWAYS_ASSERT(!volume_centroid.contains_nan()); + AMREX_ALWAYS_ASSERT(!boundary_centroid.contains_nan()); + AMREX_ALWAYS_ASSERT(!boundary_normal.contains_nan()); + AMREX_ALWAYS_ASSERT(!bndry_area.contains_nan()); + AMREX_ALWAYS_ASSERT(!volume_centroid.contains_inf()); + AMREX_ALWAYS_ASSERT(!boundary_centroid.contains_inf()); + AMREX_ALWAYS_ASSERT(!boundary_normal.contains_inf()); + AMREX_ALWAYS_ASSERT(!bndry_area.contains_inf()); + AMREX_ALWAYS_ASSERT(bndry_area.min(0) >= -roundoff_tolerance); + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + AMREX_ALWAYS_ASSERT(volume_centroid.min(d) >= -0.5_rt - roundoff_tolerance); + AMREX_ALWAYS_ASSERT(volume_centroid.max(d) <= 0.5_rt + roundoff_tolerance); + // Non-cut entries in a mixed CutFab retain the EB sentinel -1. + AMREX_ALWAYS_ASSERT(boundary_centroid.min(d) >= -1.0_rt - roundoff_tolerance); + AMREX_ALWAYS_ASSERT(boundary_centroid.max(d) <= 0.5_rt + roundoff_tolerance); + AMREX_ALWAYS_ASSERT(boundary_normal.min(d) >= -1.0_rt - roundoff_tolerance); + AMREX_ALWAYS_ASSERT(boundary_normal.max(d) <= 1.0_rt + roundoff_tolerance); + } + + auto const area_fraction = factory->getAreaFrac(); + auto const face_centroid = factory->getFaceCent(); + Array, AMREX_SPACEDIM> dense_area; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + dense_area[d] = std::make_unique(area_fraction[d]->ToMultiFab(1.0_rt, 0.0_rt)); + auto face_cent = face_centroid[d]->ToMultiFab(0.0_rt, 0.0_rt); + AMREX_ALWAYS_ASSERT(!dense_area[d]->contains_nan()); + AMREX_ALWAYS_ASSERT(!face_cent.contains_nan()); + AMREX_ALWAYS_ASSERT(!dense_area[d]->contains_inf()); + AMREX_ALWAYS_ASSERT(!face_cent.contains_inf()); + AMREX_ALWAYS_ASSERT(dense_area[d]->min(0) >= -roundoff_tolerance); + AMREX_ALWAYS_ASSERT(dense_area[d]->max(0) <= 1.0_rt + roundoff_tolerance); + for (int n = 0; n < AMREX_SPACEDIM - 1; ++n) { + AMREX_ALWAYS_ASSERT(face_cent.min(n) >= -0.5_rt - roundoff_tolerance); + AMREX_ALWAYS_ASSERT(face_cent.max(n) <= 0.5_rt + roundoff_tolerance); + } + } + + auto const& flags = factory->getMultiEBCellFlagFab(); + auto const& levelset = factory->getLevelSet(); + AMREX_ALWAYS_ASSERT(!levelset.contains_nan()); + AMREX_ALWAYS_ASSERT(!levelset.contains_inf()); + + auto const edge_centroid = factory->getEdgeCent(); + Array, AMREX_SPACEDIM> dense_edge; + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + dense_edge[d] = std::make_unique(edge_centroid[d]->ToMultiFab(1.0_rt, -1.0_rt)); + AMREX_ALWAYS_ASSERT(!dense_edge[d]->contains_nan()); + AMREX_ALWAYS_ASSERT(!dense_edge[d]->contains_inf()); + AMREX_ALWAYS_ASSERT(dense_edge[d]->min(0) >= -1.0_rt - roundoff_tolerance); + AMREX_ALWAYS_ASSERT(dense_edge[d]->max(0) <= 1.0_rt + roundoff_tolerance); + } + + Gpu::Buffer levelset_errors(2); + std::fill_n(levelset_errors.hostData(), levelset_errors.size(), 0); + levelset_errors.copyToDeviceAsync(); + int* const level_errors = levelset_errors.data(); + for (MFIter mfi(flags); mfi.isValid(); ++mfi) { + auto const flag = flags.const_array(mfi); + auto const phi = levelset.const_array(mfi); + ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (flag(i, j, k).isCovered()) { + for (int kk = 0; kk <= 1; ++kk) { + for (int jj = 0; jj <= 1; ++jj) { + for (int ii = 0; ii <= 1; ++ii) { + if (phi(i + ii, j + jj, k + kk) < 0.0_rt) { + Gpu::Atomic::AddNoRet(level_errors, 1); + } + } + } + } + } + }); + } + for (int d = 0; d < AMREX_SPACEDIM; ++d) { + for (MFIter mfi(*dense_edge[d]); mfi.isValid(); ++mfi) { + auto const edge = dense_edge[d]->const_array(mfi); + auto const phi = levelset.const_array(mfi); + ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + int const ih = i + (d == 0); + int const jh = j + (d == 1); + int const kh = k + (d == 2); + bool const lo_fluid = phi(i, j, k) < 0.0_rt; + bool const hi_fluid = phi(ih, jh, kh) < 0.0_rt; + bool mismatch = false; + if (lo_fluid && hi_fluid) { + mismatch = edge(i, j, k) != 1.0_rt; + } else if (!lo_fluid && !hi_fluid) { + mismatch = edge(i, j, k) != -1.0_rt; + } else { + mismatch = edge(i, j, k) < -0.5_rt || edge(i, j, k) > 0.5_rt; + } + if (mismatch) { + Gpu::Atomic::AddNoRet(level_errors + 1, 1); + } + }); + } + } + levelset_errors.copyToHost(); + GpuArray global_levelset_errors{levelset_errors.hostData()[0], + levelset_errors.hostData()[1]}; + ParallelAllReduce::Sum(global_levelset_errors.data(), + static_cast(global_levelset_errors.size()), + ParallelContext::CommunicatorSub()); + // The legacy generator repairs the level set independently in every FAB + // and the exported level set takes an arbitrary owner on shared nodes. + // Without deterministic floating point (-ffast-math) two FABs can round + // the same borderline cell differently, so cross-FAB agreement of the + // exported signs is only guaranteed for the marching-cubes generator. +#ifdef __FAST_MATH__ + bool const check_levelset_signs = (eb_method == "marching_cubes"); +#else + bool const check_levelset_signs = true; +#endif + if (check_levelset_signs) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_levelset_errors[0] == 0, + "A covered cell retains a public negative-in-fluid level-set node"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_levelset_errors[1] == 0, + "Edge centroids disagree with repaired public level-set signs"); + } else { + amrex::Print() << "Skipping level-set sign consistency checks for the " + << "legacy generator in a fast-math build (" + << global_levelset_errors[0] << " covered-cell, " + << global_levelset_errors[1] << " edge mismatches)\n"; + } + + Gpu::DeviceScalar zero_node_count(static_cast(0)); + Long* const zero_nodes = zero_node_count.dataPtr(); + for (MFIter mfi(levelset); mfi.isValid(); ++mfi) { + auto const phi = levelset.const_array(mfi); + ParallelFor(mfi.validbox(), [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (phi(i, j, k) == 0.0_rt) { + Gpu::Atomic::AddNoRet(zero_nodes, static_cast(1)); + } + }); + } + Long repaired_nodes = zero_node_count.dataValue(); + ParallelAllReduce::Sum(repaired_nodes, ParallelContext::CommunicatorSub()); + + Real small_volfrac = 0.0_rt; + ParmParse("eb2").query("small_volfrac", small_volfrac); + Gpu::Buffer topology_counts(6); + std::fill_n(topology_counts.hostData(), topology_counts.size(), 0); + topology_counts.copyToDeviceAsync(); + int* const topology = topology_counts.data(); + Box const domain = geom.Domain(); + for (MFIter mfi(flags); mfi.isValid(); ++mfi) { + Box const bx = mfi.validbox(); + auto const flag = flags.const_array(mfi); + auto const apx = dense_area[0]->const_array(mfi); + auto const apy = dense_area[1]->const_array(mfi); + auto const apz = dense_area[2]->const_array(mfi); + auto const vf = volfrac.const_array(mfi); + auto const bndry_area_arr = bndry_area.const_array(mfi); + auto const bn = boundary_normal.const_array(mfi); + ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept { + if (flag(i, j, k).isSingleValued()) { + Gpu::Atomic::AddNoRet(topology + 1, 1); + if (vf(i, j, k) + roundoff_tolerance < small_volfrac) { + Gpu::Atomic::AddNoRet(topology + 2, 1); + } + if (vf(i, j, k) == 1.0_rt && bndry_area_arr(i, j, k) > 0.0_rt) { + Gpu::Atomic::AddNoRet(topology + 4, 1); + } + Real const area_x = apx(i, j, k) - apx(i + 1, j, k); + Real const area_y = apy(i, j, k) - apy(i, j + 1, k); + Real const area_z = apz(i, j, k) - apz(i, j, k + 1); + Real const norm = + std::sqrt(area_x * area_x + area_y * area_y + area_z * area_z); + Real const normal_tolerance = + amrex::max(5.e-12_rt, roundoff_tolerance * amrex::max(1.0_rt, norm)); + if (norm <= 1.e-12_rt || + std::abs(norm * bn(i, j, k, 0) - area_x) > normal_tolerance || + std::abs(norm * bn(i, j, k, 1) - area_y) > normal_tolerance || + std::abs(norm * bn(i, j, k, 2) - area_z) > normal_tolerance) { + Gpu::Atomic::AddNoRet(topology + 5, 1); + } + } + if (flag(i, j, k).isMultiValued()) { + Gpu::Atomic::AddNoRet(topology + 3, 1); + } + + int mismatches = 0; + if (i > domain.smallEnd(0)) { + mismatches += flag(i, j, k).isConnected(-1, 0, 0) != (apx(i, j, k) > 0.0_rt); + } + if (i < domain.bigEnd(0)) { + mismatches += flag(i, j, k).isConnected(1, 0, 0) != (apx(i + 1, j, k) > 0.0_rt); + } + if (j > domain.smallEnd(1)) { + mismatches += flag(i, j, k).isConnected(0, -1, 0) != (apy(i, j, k) > 0.0_rt); + } + if (j < domain.bigEnd(1)) { + mismatches += flag(i, j, k).isConnected(0, 1, 0) != (apy(i, j + 1, k) > 0.0_rt); + } + if (k > domain.smallEnd(2)) { + mismatches += flag(i, j, k).isConnected(0, 0, -1) != (apz(i, j, k) > 0.0_rt); + } + if (k < domain.bigEnd(2)) { + mismatches += flag(i, j, k).isConnected(0, 0, 1) != (apz(i, j, k + 1) > 0.0_rt); + } + if (mismatches != 0) { + Gpu::Atomic::AddNoRet(topology, mismatches); + } + }); + } + topology_counts.copyToHost(); + GpuArray global_topology_counts; + std::copy_n(topology_counts.hostData(), global_topology_counts.size(), + global_topology_counts.begin()); + ParallelAllReduce::Sum(global_topology_counts.data(), + static_cast(global_topology_counts.size()), + ParallelContext::CommunicatorSub()); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_topology_counts[0] == 0, + "EBCellFlag coordinate connectivity disagrees with MC face apertures"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(global_topology_counts[1] > 0, + "Marching-cubes test did not produce any single-valued cells"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_topology_counts[2] == 0, + "Nodal repair left a single-valued cell below eb2.small_volfrac"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(global_topology_counts[3] == 0, + "Single-valued EB regression retained a multi-valued cell"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_topology_counts[5] == 0, + "EB boundary normals disagree with repaired face-aperture closure"); + // A single-valued cell with unit volume and a nonzero boundary area owns + // an EB patch coincident with one of its faces or nodes. This state only + // arises when the surface passes exactly through grid nodes, either + // because the geometry does (allow_full_volume_cut_cells = 1) or because + // the nodal repair moved nodes onto the surface. + if (!allow_full_volume_cut_cells && repaired_nodes == 0) { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + global_topology_counts[4] == 0, + "A single-valued cell has unit volume and a nonzero boundary area"); + } + amrex::Print() << "Single-valued cells: " << global_topology_counts[1] + << ", retained multi-valued cells: " << global_topology_counts[3] + << ", full-volume cut cells: " << global_topology_counts[4] << "\n"; + amrex::Print() << std::setprecision(std::numeric_limits::max_digits10) + << "MC_TEST_SIGNATURE backend=" << execution_backend() + << " precision=" << 8 * sizeof(Real) << " method=" << eb_method + << " fluid_volume=" << fluid_volume + << " single_valued_cells=" << global_topology_counts[1] + << " repaired_nodes=" << repaired_nodes << '\n'; + if (!custom_stl_test) { + Real const domain_scale = amrex::max(xmax - xmin, amrex::max(ymax - ymin, zmax - zmin)); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + error < 6.0_rt * max_dx * max_dx * domain_scale, + "Marching-cubes cut-cell volume error exceeds the expected O(dx^2) bound"); + } } int main (int argc, char* argv[]) diff --git a/Tests/EB/MarchingCubes/open_surface.stl b/Tests/EB/MarchingCubes/open_surface.stl new file mode 100644 index 00000000000..6b1b659203e --- /dev/null +++ b/Tests/EB/MarchingCubes/open_surface.stl @@ -0,0 +1,23 @@ +solid open_tetrahedron + facet normal 0 0 -1 + outer loop + vertex 0 0 0 + vertex 0 1 0 + vertex 1 0 0 + endloop + endfacet + facet normal 0 -1 0 + outer loop + vertex 0 0 0 + vertex 1 0 0 + vertex 0 0 1 + endloop + endfacet + facet normal -1 0 0 + outer loop + vertex 0 0 0 + vertex 0 0 1 + vertex 0 1 0 + endloop + endfacet +endsolid open_tetrahedron diff --git a/Tests/EB/MarchingCubes/stanford_bunny.stl b/Tests/EB/MarchingCubes/stanford_bunny.stl new file mode 100644 index 00000000000..bec02fcb5de Binary files /dev/null and b/Tests/EB/MarchingCubes/stanford_bunny.stl differ