PointinPolygon for determing mapped idata in 2D

RevathiJambunathan · Jan 25, 2024 · 66ba468 · 66ba468
1 parent 380ef4b
commit 66ba468
Showing 1 changed file with 60 additions and 12 deletions.
diff --git a/Tests/Particles/Mapped/MappedPC.H b/Tests/Particles/Mapped/MappedPC.H
@@ -82,6 +82,27 @@ public:
         Redistribute(lev_min, lev_max, nGrow, local, remove_neg);
     }
 
+
+    /**
+     * Given the vertices of a polygon (x,y), this function
+     * returns true if the point(xp,yp) is inside the polygon.
+     */
+    bool PointInPolygon (amrex::Real x[4], amrex::Real y[4], amrex::Real xp, amrex::Real yp)
+    {
+        bool inside;
+        int pos = 0;
+        int neg = 0;
+        for (int v = 0; v < 4; ++v) {
+            int w = (v < 3) ? v+1 : 0;
+            amrex::Real CP = (xp - x[v])*(y[w] - y[v]) - ( yp - y[v])*(x[w] - x[v]);
+            if (CP > 0) { pos ++;
+            } else { neg++;
+            }
+        }
+        inside = ( (pos > 0) && (neg > 0) ) ? false : true;
+        return inside;
+    }
+
     template <typename P>
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     void update_mapped_idata (P& p,
@@ -90,18 +111,45 @@ public:
                               const amrex::Array4<amrex::Real const>&  loc_arr)
     {
 #if (AMREX_SPACEDIM == 2)
-        amrex::IntVect iv( int(amrex::Math::floor((p.pos(0)-plo[0])*dxi[0])), p.idata(1) );
-
-        amrex::Real lx = (p.pos(0)-plo[0])*dxi[0] - static_cast<amrex::Real>(iv[0]);
-        auto ylo = loc_arr(iv[0]  ,iv[1]  ,0,1) * (1.0-lx) +
-                   loc_arr(iv[0]+1,iv[1]  ,0,1) *      lx;
-        auto yhi = loc_arr(iv[0]  ,iv[1]+1,0,1) * (1.0-lx) +
-                   loc_arr(iv[0]+1,iv[1]+1,0,1) *      lx;
-
-        if (p.pos(1) > yhi) {
-            p.idata(1) += 1;
-        } else if (p.pos(1) <= ylo) {
-            p.idata(1) -= 1;
+        amrex::IntVect iv( p.idata(0), p.idata(1) );
+        amrex::Real x[5];
+        amrex::Real y[5];
+        int cnt = 0;
+        int k = 0;
+        for (int jj = 0; jj < 2; ++jj) {
+            for (int ii = 0; ii < 2; ++ii) {
+                // the polygon vertices are counterclockwise
+                // i.e., (ii,jj), (ii+1,jj), (ii+1,jj+1), (ii,jj+1)
+                int i = amrex::Math::abs( ii - jj);
+                x[cnt] = loc_arr(iv[0] + i, iv[1] + jj, k, 0);
+                y[cnt] = loc_arr(iv[0] + i, iv[1] + jj, k, 1);
+            }
+        }
+        bool inside = PointInPolygon( x, y, p.pos(0), p.pos(1));
+
+        if (inside) { return;
+        } else {
+            // Looping over neighboring cells
+            for (int lo_j = -1; lo_j < 2; ++lo_j) {
+            for (int lo_i = -1; lo_i < 2; ++lo_i) {
+                if ( (lo_i == 0) && (lo_j == 0) ) continue;
+                amrex::IntVect ivn( p.idata(0) + lo_i, p.idata(1) + lo_j);
+                cnt = 0;
+                // looping over vertices of the polygon to set x,y at each vertex
+                for (int jj = 0; jj < 2; ++jj) {
+                for (int ii = 0; ii < 2; ++ii) {
+                    int i = amrex::Math::abs(ii - jj);
+                    x[cnt] = loc_arr(ivn[0] + i, ivn[1] + jj, k, 0);
+                    y[cnt] = loc_arr(ivn[0] + i, ivn[1] + jj, k, 1);
+                }}
+                inside = PointInPolygon( x, y, p.pos(0), p.pos(1));
+
+                if (inside) {
+                    p.idata(0) = p.idata(0) + lo_i;
+                    p.idata(1) = p.idata(1) + lo_j;
+                    return;
+                }
+            }}
         }
 
 #elif (AMREX_SPACEDIM == 3)