Changeset 2500

Timestamp:

05/12/23 11:34:43 (14 months ago)

Author:

ymipsl

Message:

Modification in remaper to fix problems with non convex cells.

YM

Location:

XIOS2/trunk/extern/remap/src

Files:

• meshutil.cpp (modified) (4 diffs)
• polyg.cpp (modified) (4 diffs)
• polyg.hpp (modified) (2 diffs)

XIOS2/trunk/extern/remap/src/meshutil.cpp

@@ -11 +11 @@
 using namespace std;
 
-double computePolygoneArea(Elt& a, const Coord &pole)
+void computePolygonGeometry(Elt& a, const Coord &pole, double& area, Coord& bary)
 {
   using N = uint32_t;
@@ -49 +49 @@
   vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
   
-  double area_a_gno=0 ;
+  area=0 ;
   for(int i=0;i<indices_a_gno.size()/3;++i)
     {
@@ -55 +55 @@
       Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
       Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
-      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
-    }
-
+      area+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
+    }
+
+
+  bary = exact_barycentre(dstPolygon.data(),na) ;
+  
+  // check signed area
+  double signedArea = 0 ;
+  for(int n=0; n<na;n++) signedArea+= a_gno[n].x*a_gno[(n+1)%na].y-a_gno[(n+1)%na].x*a_gno[n].y ;
+  if (signedArea<0) 
+  {
+    bary = bary * (-1.) ;
+    switchOrientation(a.n, a.vertex,a.edge,a.d) ;
+  }
+  
   vect_points.clear();
   polyline.clear();
   indices_a_gno.clear();
-  return area_a_gno ;
+
 }
 
@@ -67 +79 @@
 void cptEltGeom(Elt& elt, const Coord &pole)
 {
-  orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
-  normals(elt, pole);
-  Coord gg;
-  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
-  elt.x = gg;
+   orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
+   normals(elt, pole);
+//  Coord gg;
+//  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
+//  elt.x = gg;
 // overwrite area computation 
 
-  elt.area =  computePolygoneArea(elt, pole) ;
+   computePolygonGeometry(elt, pole, elt.area, elt.x) ;
+   normals(elt, pole);
+
 }
 

XIOS2/trunk/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;
@@ -99 +104 @@
   else if (n == 1) return x[0];
 }
+
+/* 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)

XIOS2/trunk/extern/remap/src/polyg.hpp

@@ -7 +7 @@
 
 void orient(int n, Coord *vertex, Coord *edge, double *d, const Coord &g);
-
+void switchOrientation(int N, Coord *vertex, Coord *edge, double *d) ;
 void normals(Coord *x, int n, Coord *a);
 
@@ -16 +16 @@
 double polygonarea(Coord *x, int n);
 Coord exact_barycentre(const Coord *x, int n) ;
+Coord new_barycentre(const Coord *x, int n) ;
 
 int packedPolygonSize(const Elt& e);