EB initialization with Parser (#1980)

If `warpx.eb_implicit_function = ...` is present in input parameters, Parser
will be used to initialize EB.  For example,
```
warpx.eb_implicit_function = "max(max(max(x-0.5,-0.5-x), max(y-0.5,-0.5-y)), max(z-0.5,-0.5-z))"
```
specifies a box with boundaries at x=+-0.5, y=+-0.5 and z=+-0.5 and regular
domain inside the box.
```
warpx.eb_implicit_function = "-(x**2+y**2+z**2-0.2**2)"
```
specifies a solid sphere at (0,0,0) with a radius of 0.2.

1 files changed, 50 insertions, 5 deletions

diff --git a/Source/EmbeddedBoundary/WarpXInitEB.cpp b/Source/EmbeddedBoundary/WarpXInitEB.cpp
index 8214eec9b..b0de40136 100644
--- a/Source/EmbeddedBoundary/WarpXInitEB.cpp
+++ b/Source/EmbeddedBoundary/WarpXInitEB.cpp
@@ -6,6 +6,40 @@
 #   include <AMReX_ParmParse.H>
 #endif
 
+#ifdef AMREX_USE_EB
+namespace {
+    class ParserIF
+        : public amrex::GPUable
+    {
+    public:
+        ParserIF (const ParserWrapper<3>& a_parser)
+            : m_parser(a_parser.getParser())
+            {}
+
+        ParserIF (const ParserIF& rhs) noexcept = default;
+        ParserIF (ParserIF&& rhs) noexcept = default;
+        ParserIF& operator= (const ParserIF& rhs) = delete;
+        ParserIF& operator= (ParserIF&& rhs) = delete;
+
+        AMREX_GPU_HOST_DEVICE inline
+        amrex::Real operator() (AMREX_D_DECL(amrex::Real x, amrex::Real y,
+                                             amrex::Real z)) const noexcept {
+#if (AMREX_SPACEDIM == 2)
+            return m_parser(x,amrex::Real(0.0),y);
+#else
+            return m_parser(x,y,z);
+#endif
+        }
+
+        inline amrex::Real operator() (const amrex::RealArray& p) const noexcept {
+            return this->operator()(AMREX_D_DECL(p[0],p[1],p[2]));
+        }
+
+    private:
+        HostDeviceParser<3> m_parser;
+    };
+}
+#endif
 
 void
 WarpX::InitEB ()
@@ -13,12 +47,23 @@ WarpX::InitEB ()
 #ifdef AMREX_USE_EB
     BL_PROFILE("InitEB");
 
-    amrex::ParmParse pp_eb2("eb2");
-    if (!pp_eb2.contains("geom_type")) {
-        std::string geom_type = "all_regular";
-        pp_eb2.add("geom_type", geom_type); // use all_regular by default
+    amrex::ParmParse pp_warpx("warpx");
+    std::string impf;
+    pp_warpx.query("eb_implicit_function", impf);
+    if (! impf.empty()) {
+        WarpXParser wp = makeParser(impf, {"x", "y", "z"});
+        m_eb_if_parser = std::make_unique<ParserWrapper<3> >(wp);
+        ParserIF pif(*m_eb_if_parser);
+        auto gshop = amrex::EB2::makeShop(pif);
+        amrex::EB2::Build(gshop, Geom(maxLevel()), maxLevel(), maxLevel());
+    } else {
+        amrex::ParmParse pp_eb2("eb2");
+        if (!pp_eb2.contains("geom_type")) {
+            std::string geom_type = "all_regular";
+            pp_eb2.add("geom_type", geom_type); // use all_regular by default
+        }
+        amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel());
     }
-    amrex::EB2::Build(Geom(maxLevel()), maxLevel(), maxLevel());
 
 #endif
 }