Changeset 1016

Timestamp:

01/09/17 14:33:53 (7 years ago)

Author:

ymipsl

Message:

Second order interpolation : is now able to detect a cell as neighbour if only a part of the edge is common with the target cell. Fix some problem with reduced gaussian grid which cannot detect correct neighboround on the latitude.

YM

Location:

XIOS/trunk/extern/remap/src

Files:

• elt.hpp (modified) (1 diff)
• intersect.cpp (modified) (2 diffs)
• intersect.hpp (modified) (1 diff)
• triple.cpp (modified) (1 diff)
• triple.hpp (modified) (1 diff)

XIOS/trunk/extern/remap/src/elt.hpp

@@ -105 +105 @@
         }
 
+  void insert_vertex(int i, const Coord& v)
+  {
+    for(int j=n; j > i ; j--)
+    {
+      vertex[j]=vertex[j-1] ;
+      edge[j]=edge[j-1] ;
+      d[j]=d[j-1] ;
+      neighbour[j]=neighbour[j-1] ;
+    }
+    vertex[i+1]=v ;
+    n++ ;
+  }
+  
         int neighbour[NMAX];
         double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */

XIOS/trunk/extern/remap/src/intersect.cpp

@@ -37 +37 @@
 }
 
+/** New methods to find an insert a neighbour in a cell of the source mesh.
+ *  return true/false if cell b is a neighbour of a. if "insert" is true, then b will be inserted as a neighbour
+ * in cell a . This is needed for 2 order interpolation that need neighboround for gradient computing.
+ * A cell is a neighbour if :
+ *  - it shares 2 countiguous vertex (ie an edge) with a
+ *  - A vertex of b is located on one of an edge of a.
+ **/
+bool insertNeighbour( Elt& a, const Elt& b, bool insert )
+{
+  // for now suppose pole -> Oz
+  Coord pole(0,0,1) ;
+  Coord O, Oa1, Oa2,Ob1,Ob2,V1,V2 ;
+  double da,db,alpha,alpha1,alpha2,delta ;
+  
+    
+  for (int i = 0; i < a.n; i++)
+  {
+    for (int j = 0; j < b.n; j++)
+    {
+// share a full edge ? be carefull at the orientation
+      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > 1e-10*1e-10 ||
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10);
+      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-10*1e-10 &&
+             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-10*1e-10)
+      {
+        if (insert) a.neighbour[i] = b.id.ind ;
+        return true;
+      }
+      
+
+// 1 or 2 vertices of b is located on an edge of a
+       da=a.d[i] ;
+       if (scalarprod(a.edge[i], pole) < 0) da=-da ;
+       db=b.d[(j+b.n-1)%b.n] ;
+       if (scalarprod(b.edge[(j+b.n-1)%b.n], pole) < 0) db=-db ;
+      
+      if ( fabs(da-db)<1e-10 ) 
+      {
+        O=pole*da ;
+        Oa1=a.vertex[i]-O ;
+        Oa2=a.vertex[(i+1)%a.n]-O ; 
+        Ob1=b.vertex[j]-O ;
+        Ob2=b.vertex[(j+b.n-1)%b.n]-O ;
+        V1=crossprod(Oa1,Oa2) ;
+        V2=crossprod(Ob1,Ob2) ;
+        if (norm(crossprod(V1,V2))/(norm(V1)*norm(V2)) < 1e-10)
+        {
+          alpha = vectAngle(Oa1,Oa2,V1) ;
+          alpha1= vectAngle(Oa1,Ob1,V1) ;
+          alpha2= vectAngle(Oa1,Ob2,V1) ;
+          delta= alpha2-alpha1 ;
+          if (delta >= M_PI) delta=2*M_PI-delta ;
+          else if (delta <= -M_PI) delta=2*M_PI+delta ;
+          
+          if (alpha >= 0)
+          {
+            if (alpha1 > 1e-10 && alpha1 < alpha-1e-10)
+            {
+              if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
+              {
+                assert(delta > 0) ;
+                if (insert)
+                {
+                // insert both
+                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+              else
+              {
+                assert( delta > 0 ) ;
+                if (insert)
+                {
+                //insert alpha1
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+            }
+            else if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
+            {
+              assert( delta > 0 ) ;
+              if (insert)
+              {
+              // insert alpha2
+                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                a.neighbour[i] = b.id.ind ;
+              }
+              return true ;
+            }
+            else
+            {
+              // nothing to do
+            } 
+
+          }
+          else  // alpha < 0
+          {
+            if (alpha1 < -1e-10 && alpha1 > alpha+1e-10)
+            {
+              if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
+              {
+                assert(delta < 0) ;
+                if (insert)
+                {
+                // insert both
+                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+              else
+              {
+                assert(delta < 0) ;
+                if (insert)
+                {
+                //insert alpha1
+                  a.insert_vertex(i,b.vertex[j]);
+                  a.neighbour[i+1] = b.id.ind ;
+                }
+                return true ;
+              }
+            }
+            else if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
+            {
+              assert(delta < 0) ;
+              if (insert)
+              {
+              // insert alpha2
+                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
+                a.neighbour[i] = b.id.ind ;
+              }
+              return true ;
+            }
+            else
+            {
+               // nothing to do
+            }
+          }
+        }
+      }
+    }
+  }
+  return false;
+}
 
 /**
@@ -44 +193 @@
 void set_neighbour(Elt& a, const Elt& b)
 {
-        if (b.id.ind == a.id.ind) return;
-        int idx = neighbour_idx(a, b);
-        if (idx != NOT_FOUND)
-                a.neighbour[idx] = b.id.ind;
+  if (b.id.ind == a.id.ind) return;
+/*
+  int idx = neighbour_idx(a, b);
+  if (idx != NOT_FOUND)
+  a.neighbour[idx] = b.id.ind;
+*/
+  insertNeighbour(a,b,true) ; 
 }
 
 /** return true if `a` and `b` share an edge */
-bool isNeighbour(const Elt& a, const Elt& b)
-{
-        return neighbour_idx(a, b) != NOT_FOUND;
+bool isNeighbour(Elt& a, const Elt& b)
+{
+        // return neighbour_idx(a, b) != NOT_FOUND;
+  return insertNeighbour(a,b,false) ;
 }
 

XIOS/trunk/extern/remap/src/intersect.hpp

@@ -2 +2 @@
 
 void set_neighbour(Elt& elt, const Elt& elt2);
-bool isNeighbour(const Elt& elt, const Elt& elt2);
+bool isNeighbour(Elt& elt, const Elt& elt2);
 
 void intersect(Elt *a, Elt *b);

XIOS/trunk/extern/remap/src/triple.cpp

@@ -124 +124 @@
 }
 
+// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b
+double vectAngle(const Coord &a, const Coord &b, const Coord &pole)
+{ 
+  double nab = 1./(norm(a)*norm(b)) ;
+  
+  Coord a_cross_b=crossprod(a, b)*nab ;
+  double sinVect ;
+  if (scalarprod(a_cross_b, pole) >= 0) sinVect=norm(a_cross_b) ;
+  else sinVect=-norm(a_cross_b) ;
+  double cosVect=scalarprod(a,b)*nab ;
+
+  return atan2(sinVect,cosVect) ;
 }
+
+}

XIOS/trunk/extern/remap/src/triple.hpp

@@ -121 +121 @@
 
 double angle(const Coord &a, const Coord &b, const Coord &pole);
+
+// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b
+double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ;
+
 void print_Coord(Coord &p);