Add STL geom_type for EB in Source/InitEB.cpp

Docs/sphinx/geometry/geometry_init.rst

@@ -3,40 +3,40 @@
 Geometry initialization
 -----------------------
 
-Creating an EB geometry also requires knowledge of the finest level that will be used so that geometries that 'telescope', 
-i.e., coarser volume fractions are consistent with applying the coarsening operator to the finer volumes, can be created. 
-To that end there is a global geometry creation step, facilitated by the `initialize_EB2` function, as well as a step that 
-happens when a new AMRLevel is created. The latter happens by a call to  `PeleC::initialize_eb2_structs`  through `PeleC::init_eb` 
-called from the PeleC constructor. Following construction of the geometry, the geometric information is 
-copied into the structures described in the previous section and the various interpolation stencils are populated. 
-
-Cartesian grid, embedded boundary (EB) methods are methods where the geometric description is formed by cutting a Cartesian 
+Creating an EB geometry also requires knowledge of the finest level that will be used so that geometries that 'telescope',
+i.e., coarser volume fractions are consistent with applying the coarsening operator to the finer volumes, can be created.
+To that end there is a global geometry creation step, facilitated by the `initialize_EB2` function, as well as a step that
+happens when a new AMRLevel is created. The latter happens by a call to  `PeleC::initialize_eb2_structs`  through `PeleC::init_eb`
+called from the PeleC constructor. Following construction of the geometry, the geometric information is
+copied into the structures described in the previous section and the various interpolation stencils are populated.
+
+Cartesian grid, embedded boundary (EB) methods are methods where the geometric description is formed by cutting a Cartesian
 mesh with surface of the geometry.  AMReX's methods to handle EB geometry information, and PeleC's treatment of the
 EB aware update could use many possible sources for geometric description. The necessary information is, on a per-cell basis:
 
 * Apertures for faces intersected by cut cells,
 * cut cell volumes that 'telescope', that is, volumes at a coarser level are consistent with averaging the volumes from finer levels,
 * connectivity indicating which neighbor cells are connected to a given cell, and
-* coordinates of cell and face centroids. 
+* coordinates of cell and face centroids.
 
-Additionally, the algorithms ultimately require surface normals, but these can be trivially recomputed from the aperture. 
+Additionally, the algorithms ultimately require surface normals, but these can be trivially recomputed from the aperture.
 
 GeometryShop and Implicit Functions
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 
-One of the greatest advantages of EB technology is that grid generation is robust and fast and can be done to any accuracy 
+One of the greatest advantages of EB technology is that grid generation is robust and fast and can be done to any accuracy
 as described by `Schwartz et al. <http://dx.doi.org/10.2140/camcos.2015.10.83>`_ The foundation class that AMReX uses for
-geometry generation is called `GeometryShop`. This class is used to initialize geometric information 
-and associated connectivity graph stored in a distributed database class `EBIndexSpace`. 
-Historically, the `EBIndexSpace` database was developed to be used throughout a calculation. 
-Here, we use it only to populate datastructures that can be accessed efficiently in the patterns 
-representative of the Pele motivating problem space. 
-
- Given an implicit function :math:`I`, ``GeometryShop`` interprets the surface upon which 
- :math:`I(\mathbf{x}) = 0` as the surface with which to cut the grid cells. 
- ``GeometryShop`` interprets the positive regions of the implicit function (:math:`\mathbf{x}: I(\mathbf{x}) > 0`) 
- as covered by the geometry and negative regions (:math:`\mathbf{x}: I(\mathbf{x}) < 0`) as part of  the solution domain.  
+geometry generation is called `GeometryShop`. This class is used to initialize geometric information
+and associated connectivity graph stored in a distributed database class `EBIndexSpace`.
+Historically, the `EBIndexSpace` database was developed to be used throughout a calculation.
+Here, we use it only to populate datastructures that can be accessed efficiently in the patterns
+representative of the Pele motivating problem space.
+
+ Given an implicit function :math:`I`, ``GeometryShop`` interprets the surface upon which
+ :math:`I(\mathbf{x}) = 0` as the surface with which to cut the grid cells.
+ ``GeometryShop`` interprets the positive regions of the implicit function (:math:`\mathbf{x}: I(\mathbf{x}) > 0`)
+ as covered by the geometry and negative regions (:math:`\mathbf{x}: I(\mathbf{x}) < 0`) as part of  the solution domain.
  For example, if one defines her implicit function :math:`S` as
 
 .. math:
@@ -56,13 +56,13 @@ There are several basic geometries that are available in AMReX that can be easil
 * *Cylinder* - needs center (cylinder_center), radius (cylinder_radius), height (cylinder_height), direction (cylinder_direction) and fluid inside/outside flag (cylinder_has_fluid_inside).
 * *Box*      - needs the lower corner (box_lo), upper corner (box_hi) and fluid inside/outside flag (box_has_fluid_inside). The box is aligned along coordinate directions.
 * *Spline*   - needs a vector of points to create a 2D function that is a combination of spline and line elements. Currently, this geometry does not have a user interface
-  from the inputs file, but can be used within Pelec_init_eb.cpp with hard coded points. see example in section :ref:`Complicated geometries`<complexGeom>`/ 
+  from the inputs file, but can be used within Pelec_init_eb.cpp with hard coded points. see example in section :ref:`Complicated geometries`<complexGeom>`/
 
 .. code::
 
     eb2.geom_type = box
     eb2.box_lo =  -2.0  -2.0 -2.0
-    eb2.box_hi =   2.0   2.0  2.0 
+    eb2.box_hi =   2.0   2.0  2.0
     eb2.box_has_fluid_inside = 0
 
 
@@ -81,7 +81,8 @@ Adding complicated geometries
 .. _complexGeom:
 
 Geometries beyond the set described above can be built using a combination of basic geometries and EB transformation functions in AMReX.
-It should be noted that building a generic geometry from a user-defined discretized surface (like STL files)  is currently being developed, nonetheless 
+It should be noted that building a generic geometry from a user-defined discretized surface (like STL files)  is currently being developed. An example of this capability is available in ``Exec/Regtests/EB-C9``, but it should be considered to be in beta and
+potentially unstable. Nonetheless,
 engineering relevant geometries can be achieved with the fundamental geometries and transformations.
 
 Some of the relevant transformation handles in AMReX are:
@@ -93,7 +94,7 @@ Some of the relevant transformation handles in AMReX are:
 * *Lathe*       - creates a 3D implicit function from a 2D function by revolving about the z axis (see AMReX_EB2_IF_Lathe.cpp)
 * *Extrusion*   - creates a 3D implicit function from a 2D function by translating along the z axis (see AMReX_EB2_IF_Extrusion.cpp)
 
-The user can copy the file "PeleC_init_eb.cpp" from the Source and add it to his/her test case after which a new geometry can be added in initialize_EB2 
+The user can copy the file "PeleC_init_eb.cpp" from the Source and add it to his/her test case after which a new geometry can be added in initialize_EB2
 function. An example of adding a piston-bowl geometry
 that uses splines, cylinder, lathe and union transform, is shown below.
 
@@ -105,11 +106,11 @@ that uses splines, cylinder, lathe and union transform, is shown below.
     //spline IF object
     EB2::SplineIF Piston;
 
-    // array of points 
+    // array of points
     std::vector<amrex::RealVect> splpts;
 
     amrex::RealVect p;
-    // fill array of points 
+    // fill array of points
     p = amrex::RealVect(D_DECL(36.193*0.1, 7.8583*0.1, 0.0));
     spltpts.push_back(p);
     p = amrex::RealVect(D_DECL(35.924*0.1, 7.7881*0.1, 0.0));
@@ -132,11 +133,11 @@ that uses splines, cylinder, lathe and union transform, is shown below.
     .
     .
     .
-    
+
     //add to straight line elements in splineIF
     Piston.addLineElement(lnpts);
 
-    //create a cylinder 
+    //create a cylinder
     EB2::CylinderIF cylinder(48.0*0.1, 70.0*0.1, 2, {0.0, 0.0, -10.0*0.1}, true);
 
     //revolve the spline IF
@@ -146,7 +147,7 @@ that uses splines, cylinder, lathe and union transform, is shown below.
     auto PistonCylinder = EB2::makeUnion(revolvePiston, cylinder);
     auto gshop = EB2::makeShop(PistonCylinder);
 
-    
+
 #.. _EB_pistonbowl:
 
 .. figure:: EB_PistonBowl.png

Exec/RegTests/EB-C9/eb-c9.inp

@@ -61,10 +61,16 @@ amr.plot_file       = plt
 amr.plot_int        = 1000
 amr.derive_plot_vars=ALL
 
-eb2.geom_type = "cylinder"
-eb2.cylinder_direction = 0
-eb2.cylinder_center = 0.0 0.0 0.0
-eb2.cylinder_radius = 25.0
-eb2.cylinder_height = 1000.0
-eb2.cylinder_has_fluid_inside = 1
-ebd.boundary_grad_stencil_type = 0
+eb2.geom_type = stl
+eb2.stl_file    = cylinder-r50.stl
+eb2.stl_scale = 1                 # default is 1
+eb2.stl_center = -55  0  0  # default is (0,0,0)
+eb2.stl_reverse_normal = 1   # default is 0
+
+# eb2.geom_type = "cylinder"
+# eb2.cylinder_direction = 0
+# eb2.cylinder_center = 0.0 0.0 0.0
+# eb2.cylinder_radius = 25.0
+# eb2.cylinder_height = 1000.0
+# eb2.cylinder_has_fluid_inside = 1
+# ebd.boundary_grad_stencil_type = 0

Source/InitEB.cpp

@@ -351,7 +351,8 @@ initialize_EB2(
   // Custom types defined here - all_regular, plane, sphere, etc, will get
   // picked up by default (see AMReX_EB2.cpp around L100 )
   amrex::Vector<std::string> amrex_defaults(
-    {"all_regular", "box", "cylinder", "plane", "sphere", "torus", "parser"});
+    {"all_regular", "box", "cylinder", "plane", "sphere", "torus", "parser",
+     "stl"});
   if (!(std::find(amrex_defaults.begin(), amrex_defaults.end(), geom_type) !=
         amrex_defaults.end())) {
     std::unique_ptr<pele::pelec::Geometry> geometry(
@@ -365,8 +366,12 @@ initialize_EB2(
 
   // Add finer level, might be inconsistent with the coarser level created
   // above.
-  if (geom_type != "chkfile") {
+  // EY: This condition is not acceptable in AMReX with stl format
+  if ((geom_type != "chkfile") && (geom_type != "stl")) {
     amrex::EB2::addFineLevels(max_level - eb_max_level);
+  } else {
+    // The AMReX implementation for these does not support addFineLevels
+    AMREX_ALWAYS_ASSERT(max_level == eb_max_level);
   }
 
   bool write_chk_geom = false;

Tests/CMakeLists.txt

@@ -50,10 +50,11 @@ macro(setup_test)
     # Make working directory for test
     file(MAKE_DIRECTORY ${CURRENT_TEST_BINARY_DIR})
     # Gather all files in source directory for test
-    file(GLOB TEST_FILES "${CURRENT_TEST_SOURCE_DIR}/*.dat" "${CURRENT_TEST_SOURCE_DIR}/*.py")
+    file(GLOB TEST_FILES "${CURRENT_TEST_SOURCE_DIR}/*.dat" "${CURRENT_TEST_SOURCE_DIR}/*.py" "${CURRENT_TEST_SOURCE_DIR}/*.stl")
     # Copy files to test working directory
     file(COPY ${CURRENT_TEST_SOURCE_DIR}/${TEST_NAME}.inp DESTINATION "${CURRENT_TEST_BINARY_DIR}/")
     file(COPY ${TEST_FILES} DESTINATION "${CURRENT_TEST_BINARY_DIR}/")
+
     # Set some default runtime options for all tests
     set(RUNTIME_OPTIONS "amr.plot_file=plt amr.checkpoint_files_output=0 amr.plot_files_output=1 amrex.the_arena_is_managed=0 amrex.abort_on_unused_inputs=1")
     if(PELE_ENABLE_FPE_TRAP_FOR_TESTS AND (NOT APPLE))