From 6f72de283c38ef3e38ed7c5e8760a176434e6be3 Mon Sep 17 00:00:00 2001
From: Weiqun Zhang <WeiqunZhang@lbl.gov>
Date: Wed, 15 Jun 2022 09:37:15 -0700
Subject: [PATCH] Fix a pathological case for 2d EB (#2840)

What could happen is a cell might be cut a tiny bit at a corner such that 3
faces with an area fraction of one and one face with an area fraction of
almost one, and the volume fraction is one.  In that case, the boundary area
and centroid have been set to wrong values based on an incorrect assumption.
---
 Src/EB/AMReX_EB2_2D_C.cpp | 31 +++++++++++--------------------
 1 file changed, 11 insertions(+), 20 deletions(-)

diff --git a/Src/EB/AMReX_EB2_2D_C.cpp b/Src/EB/AMReX_EB2_2D_C.cpp
index 2828a29d9ee..bf17844658c 100644
--- a/Src/EB/AMReX_EB2_2D_C.cpp
+++ b/Src/EB/AMReX_EB2_2D_C.cpp
@@ -31,30 +31,21 @@ void set_eb_data (const int i, const int j,
     const Real nxabs = amrex::Math::abs(nx);
     const Real nyabs = amrex::Math::abs(ny);
 
-    Real x_ym;
-    Real x_yp;
-    Real signx;
-    Real y_xm;
-    Real y_xp;
-    Real signy;
+    Real x_ym, x_yp, y_xm, y_xp;
     if (nx > 0.0_rt) {
         x_ym = -0.5_rt + aym;
         x_yp = -0.5_rt + ayp;
-        signx = 1.0_rt;
     } else {
         x_ym = 0.5_rt - aym;
         x_yp = 0.5_rt - ayp;
-        signx = -1.0_rt;
     }
 
     if (ny > 0.0_rt) {
         y_xm = -0.5_rt + axm;
         y_xp = -0.5_rt + axp;
-        signy = 1.0_rt;
     } else {
         y_xm = 0.5_rt - axm;
         y_xp = 0.5_rt - axp;
-        signy = -1.0_rt;
     }
 
     barea(i,j,0) = nx*(axm-axp) + ny*(aym-ayp);
@@ -64,20 +55,20 @@ void set_eb_data (const int i, const int j,
     bnorm(i,j,0,1) = ny;
 
     if (nxabs < tiny || nyabs > almostone) {
-        barea(i,j,0) = 1.0_rt;
-        bcent(i,j,0,0) = 0.0_rt;
-        bnorm(i,j,0,0) = 0.0_rt;
-        bnorm(i,j,0,1) = signy;
         vfrac(i,j,0) = 0.5_rt*(axm+axp);
         vcent(i,j,0,0) = 0.0_rt;
-        vcent(i,j,0,1) = (0.125_rt*(ayp-aym) + ny*0.5_rt*bcent(i,j,0,1)*bcent(i,j,0,1)) / (vfrac(i,j,0) + 1.e-30_rt);
+        if (vfrac(i,j,0) > almostone) {
+            vcent(i,j,0,1) = 0.0_rt;
+        } else {
+            vcent(i,j,0,1) = (0.125_rt*(ayp-aym) + ny*0.5_rt*bcent(i,j,0,1)*bcent(i,j,0,1)) / (vfrac(i,j,0) + 1.e-30_rt);
+        }
     } else if (nyabs < tiny || nxabs > almostone) {
-        barea(i,j,0) = 1.0_rt;
-        bcent(i,j,0,1) = 0.0_rt;
-        bnorm(i,j,0,0) = signx;
-        bnorm(i,j,0,1) = 0.0_rt;
         vfrac(i,j,0) = 0.5_rt*(aym+ayp);
-        vcent(i,j,0,0) = (0.125_rt*(axp-axm) + nx*0.5_rt*bcent(i,j,0,0)*bcent(i,j,0,0)) / (vfrac(i,j,0) + 1.e-30_rt);
+        if (vfrac(i,j,0) > almostone) {
+            vcent(i,j,0,0) = 0.0_rt;
+        } else {
+            vcent(i,j,0,0) = (0.125_rt*(axp-axm) + nx*0.5_rt*bcent(i,j,0,0)*bcent(i,j,0,0)) / (vfrac(i,j,0) + 1.e-30_rt);
+        }
         vcent(i,j,0,1) = 0.0_rt;
     } else {
         Real aa = nxabs/ny;