EB Boundary Area: Fix issues for anisotropic cell size (#4131)

* In 2D, the scaling in incorrect. For example, if dx >> dy and the
boundary is parallel to the x-direction, it would produce inf for the
scaled boundary area.

* In 3D, the scaling has been modified so that it's easy to convert from
the dimensionless boundary area to physical units.

Docs/sphinx_documentation/source/EB.rst

@@ -271,6 +271,12 @@ following data:
     // embedded boundary centroid
     const MultiCutFab& getBndryCent () const;
 
+    // embedded boundary normal direction
+    const MultiCutFab& getBndryNormal () const;
+
+    // embedded boundary surface area
+    const MultiCutFab& getBndryArea () const;
+
     // area fractions
     Array<const MultiCutFab*,AMREX_SPACEDIM> getAreaFrac () const;
 
@@ -291,6 +297,17 @@ following data:
   of the data is in the range of :math:`[-0.5,0.5]`, based on each
   cell's local coordinates with respect to the regular cell's center.
 
+- **Boundary normal** is in a :cpp:`MultiCutFab` with ``AMREX_SPACEDIM``
+  components representing the unit vector pointing toward the covered part.
+
+- **Boundary area** is in a :cpp:`MultiCutFab` with a single component
+  representing the dimensionless boundary area. When the cell is isotropic
+  (i.e., :math:`\Delta x = \Delta y = \Delta z`), it's trivial to convert it
+  to physical units. If the cell size is anisotropic, the conversion
+  requires multiplying by a factor of :math:`\sqrt{(n_x \Delta y \Delta
+  z)^2 + (n_y \Delta x \Delta z)^2 + (n_z \Delta x \Delta y)^2}`, where
+  :math:`n` is the boundary normal vector.
+
 - **Face centroid** is in a :cpp:`MultiCutFab` with ``AMREX_SPACEDIM`` components.
   Each component of the data is in the range of :math:`[-0.5,0.5]`, based on
   each cell's local coordinates with respect to the embedded boundary.

Src/EB/AMReX_EB2_2D_C.cpp

@@ -30,8 +30,7 @@ void set_eb_data (const int i, const int j,
     const Real apnorm = std::hypot(daxp,dayp) + 1.e-30_rt*std::sqrt(dx[0]*dx[1]);
     const Real nx = daxp * (1.0_rt/apnorm);
     const Real ny = dayp * (1.0_rt/apnorm);
-    const Real bareascaling = std::sqrt( (nx*dx[0])*(nx*dx[0]) +
-            (ny*dx[1])*(ny*dx[1]) );
+    const Real bareascaling = std::sqrt(Math::powi<2>(nx*dx[1]) + Math::powi<2>(ny*dx[0]));
 
     const Real nxabs = std::abs(nx);
     const Real nyabs = std::abs(ny);

Src/EB/AMReX_EB2_3D_C.cpp

@@ -101,7 +101,9 @@ void set_eb_data (const int i, const int j, const int k,
     bnorm(i,j,k,0) = nx;
     bnorm(i,j,k,1) = ny;
     bnorm(i,j,k,2) = nz;
-    barea(i,j,k) = (nx*dapx/(dx[1]*dx[2]) + ny*dapy/(dx[0]*dx[2]) + nz*dapz/(dx[0]*dx[1]));
+    barea(i,j,k) = (nx*dapx + ny*dapy + nz*dapz) / (Math::powi<2>(nx*dx[1]*dx[2]) +
+                                                    Math::powi<2>(ny*dx[0]*dx[2]) +
+                                                    Math::powi<2>(nz*dx[0]*dx[1]));
 
     Real aax = 0.5_rt*(axm+axp);
     Real aay = 0.5_rt*(aym+ayp);