# Changeset 1016

Ignore:
Timestamp:
01/09/17 14:33:53 (7 years ago)
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:
5 edited

Unmodified
Removed
• ## XIOS/trunk/extern/remap/src/elt.hpp

 r923 } 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

 r688 } /** 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; } /** 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

 r688 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

 r849 } // 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

 r688 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);
Note: See TracChangeset for help on using the changeset viewer.