Changeset 2505 for XIOS3/branches/xios-3.0-beta/extern/remap/src/polyg.cpp

Timestamp:

05/23/23 16:20:10 (14 months ago)

Author:

jderouillat

Message:

Backport 2500, 2502 (convex cells) and 2503 (memory_report option) from XIOS2

File:

• XIOS3/branches/xios-3.0-beta/extern/remap/src/polyg.cpp (modified) (3 diffs)

XIOS3/branches/xios-3.0-beta/extern/remap/src/polyg.cpp

@@ -18 +18 @@
   (because computing intersection area requires both polygons to have same orientation)
 */
+
 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g)
 {
@@ -24 +25 @@
   Coord vertical = crossprod(ga, gb);
   if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
-  {
-    for (int i = 0; i < N/2; i++)
-      swap(vertex[i], vertex[N-1-i]);
-
-    for (int i = 0; i < (N-1)/2; i++)
-    {
-      swap(edge[N-2-i], edge[i]);
-      swap(d[i], d[N-2-i]);
-    }
-  }
-}
-
+    switchOrientation(N, vertex, edge, d) ;
+  
+}
+
+void switchOrientation(int N, Coord *vertex, Coord *edge, double *d)
+{
+  for (int i = 0; i < N/2; i++) swap(vertex[i], vertex[N-1-i]);
+
+  for (int i = 0; i < (N-1)/2; i++)
+  {
+    swap(edge[N-2-i], edge[i]);
+    swap(d[i], d[N-2-i]);
+  }
+  
+}
 
 void normals(Coord *x, int n, Coord *a)
@@ -94 +98 @@
   {
     return  proj(gc_normalintegral(x, n));
+    //return new_barycentre(x,n) ;
   }
   else if (n == 0) return ORIGIN;
   else if (n == 2) return midpoint(x[0], x[1]);
   else if (n == 1) return x[0];
-
-  error_exit( "Missing return in : Coord exact_barycentre(const Coord *x, int n)" );
-  return ORIGIN;
-}
+}
+
+/* other methode to compute barycenter of spherical polygon
+   for a spherical polygon, the moment is half the sum of (a x b) / ||a x b|| * (angle between a and b) 
+   for each pair of consecutive vertices a,b.
+*/ 
+
+Coord new_barycentre(const Coord *x, int n)
+{
+  if (n >= 3)
+  {
+    Coord sum=ORIGIN ;
+    for (int i = 0; i < n; i++)
+    { 
+      sum=sum+crossprod(x[i],x[(i+1)%n])*(1./norm(crossprod(x[i],x[(i+1)%n]))*2.*atan2(norm(x[i]-x[(i+1)%n]),norm(x[i]+x[(i+1)%n]))) ; 
+    }
+    return proj(sum) ; 
+  }
+  else if (n == 0) return ORIGIN;
+  else if (n == 2) return midpoint(x[0], x[1]);
+  else if (n == 1) return x[0];
+}
+
+
 
 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole)