EB anisotropic meshes - changes for cut_face_2d routine (#4554)

This PR makes the necessary changes in the variables - `areafrac`, `nx`,
`ny`, `centx`, `centy`, `Sx2`, `Sy2`, `Sxy`, in the `cut_face_2d`
routine, for anisotropic meshes (ie. when mesh spacings are not equal).
The derivation is in the attached pdf in the link below. The changes are
in Eqns. 35-37, 43-44, 46-47, 49-50, 55-56.


[EB_anisotropic.pdf](https://github.com/user-attachments/files/21201504/EB_anisotropic.pdf)

Author: nataraj2

Src/EB/AMReX_EB2_3D_C.cpp

@@ -195,7 +195,8 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
 void cut_face_2d (Real& areafrac, Real& centx, Real& centy,
                   Real& Sx2, Real& Sy2, Real& Sxy,
                   Real axm, Real axp, Real aym, Real ayp,
-                  Real bcx, Real bcy) noexcept
+                  Real bcx, Real bcy,
+                  Real dxval, Real dyval) noexcept
 {
 #ifdef AMREX_USE_FLOAT
     constexpr Real sml = 1.e-5_rt;
@@ -204,9 +205,9 @@ void cut_face_2d (Real& areafrac, Real& centx, Real& centy,
     constexpr Real sml = 1.e-14;
     constexpr Real tiny  = 1.e-15;
 #endif
-    Real apnorm = std::hypot(axm-axp,aym-ayp);
-    Real nx = (axm-axp) * (1.0_rt/apnorm); // pointing to the wall
-    Real ny = (aym-ayp) * (1.0_rt/apnorm);
+    Real apnorm = std::hypot((axm-axp)*dyval,(aym-ayp)*dxval);
+    Real nx = dyval * (axm-axp) * (1.0_rt/apnorm); // pointing to the wall
+    Real ny = dxval * (aym-ayp) * (1.0_rt/apnorm);
 
     Real nxabs = std::abs(nx);
     Real nyabs = std::abs(ny);
@@ -264,9 +265,9 @@ void cut_face_2d (Real& areafrac, Real& centx, Real& centy,
         Real dx2 = dx * (x_ym + x_yp);
         Real dx3 = dx * (x_ym*x_ym + x_ym*x_yp + x_yp*x_yp);
         Real dx4 = dx * (x_ym + x_yp) * (x_ym*x_ym + x_yp*x_yp);
-        Real af1 = 0.5_rt*(axm+axp) + aa*0.5_rt*dx2;
-        centx = 0.125_rt*(axp-axm) + aa*(1.0_rt/6._rt)*dx3;
-        Sx2 = (1.0_rt/24._rt)*(axm+axp) + aa*(1.0_rt/12._rt)*dx4;
+        Real af1 = 0.5_rt*(axm+axp) + aa*0.5_rt*dx2*dxval/dyval;
+        centx = 0.125_rt*(axp-axm) + aa*(1.0_rt/6._rt)*dx3*dxval/dyval;
+        Sx2 = (1.0_rt/24._rt)*(axm+axp) + aa*(1.0_rt/12._rt)*dx4*dxval/dyval;
 
         Real signy = (ny > 0.0_rt) ? 1.0_rt : -1.0_rt;
         Real y_xm = (-0.5_rt + axm)*signy;
@@ -276,13 +277,13 @@ void cut_face_2d (Real& areafrac, Real& centx, Real& centy,
         Real dy2 = dy * (y_xm + y_xp);
         Real dy3 = dy * (y_xm*y_xm + y_xm*y_xp + y_xp*y_xp);
         Real dy4 = dy * (y_xm + y_xp) * (y_xm*y_xm + y_xp*y_xp);
-        Real af2 = 0.5_rt*(aym+ayp) + aa*0.5_rt*dy2;
-        centy = (1.0_rt/8._rt)*(ayp-aym) + aa*(1.0_rt/6._rt)*dy3;
-        Sy2 = (1.0_rt/24._rt)*(aym+ayp) + aa*(1.0_rt/12._rt)*dy4;
+        Real af2 = 0.5_rt*(aym+ayp) + aa*0.5_rt*dy2*dyval/dxval;
+        centy = (1.0_rt/8._rt)*(ayp-aym) + aa*(1.0_rt/6._rt)*dy3*dyval/dxval;
+        Sy2 = (1.0_rt/24._rt)*(aym+ayp) + aa*(1.0_rt/12._rt)*dy4*dyval/dxval;
 
         Real S_b = (nxabs < nyabs)
-            ? (Sx2 - (1.0_rt/24._rt) - signx*(1.0_rt/6._rt)*(x_ym*x_ym*x_ym+x_yp*x_yp*x_yp)) / ny
-            : (Sy2 - (1.0_rt/24._rt) - signy*(1.0_rt/6._rt)*(y_xm*y_xm*y_xm+y_xp*y_xp*y_xp)) / nx;
+            ? (Sx2 - (1.0_rt/24._rt) - signx*(1.0_rt/6._rt)*(x_ym*x_ym*x_ym+x_yp*x_yp*x_yp)) / ny * dxval/dyval
+            : (Sy2 - (1.0_rt/24._rt) - signy*(1.0_rt/6._rt)*(y_xm*y_xm*y_xm+y_xp*y_xp*y_xp)) / nx * dyval/dxval;
         Sxy = (nxabs < nyabs)
             ? -signy*(1.0_rt/16._rt)*dy2 + 0.5_rt*nx*S_b
             : -signx*(1.0_rt/16._rt)*dx2 + 0.5_rt*ny*S_b;
@@ -482,7 +483,7 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 bcz = 0.5_rt*bcz - 0.5_rt;
                 cut_face_2d(apx(i,j,k),fcx(i,j,k,0),fcx(i,j,k,1), // NOLINT(readability-suspicious-call-argument)
                             m2x(i,j,k,0),m2x(i,j,k,1),m2x(i,j,k,2),
-                            lzm,lzp,lym,lyp,bcy,bcz);
+                            lzm,lzp,lym,lyp,bcy,bcz,dx[1],dx[2]);
             }
 
             if (apx(i,j,k) == 0.0_rt) {
@@ -590,7 +591,7 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 bcz = 0.5_rt*bcz - 0.5_rt;
                 cut_face_2d(apy(i,j,k),fcy(i,j,k,0),fcy(i,j,k,1), // NOLINT(readability-suspicious-call-argument)
                             m2y(i,j,k,0),m2y(i,j,k,1),m2y(i,j,k,2),
-                            lzm,lzp,lxm,lxp,bcx,bcz);
+                            lzm,lzp,lxm,lxp,bcx,bcz,dx[0],dx[2]);
             }
 
             if (apy(i,j,k) == 0.0_rt) {
@@ -698,7 +699,7 @@ int build_faces (Box const& bx, Array4<EBCellFlag> const& cell,
                 bcy = 0.5_rt*bcy - 0.5_rt;
                 cut_face_2d(apz(i,j,k),fcz(i,j,k,0),fcz(i,j,k,1), // NOLINT(readability-suspicious-call-argument)
                             m2z(i,j,k,0),m2z(i,j,k,1),m2z(i,j,k,2),
-                            lym,lyp,lxm,lxp,bcx,bcy);
+                            lym,lyp,lxm,lxp,bcx,bcy,dx[0],dx[1]);
             }
 
             if (apz(i,j,k) == 0.0_rt) {