Changeset 1016


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

Legend:

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

    r923 r1016  
    105105        } 
    106106 
     107  void insert_vertex(int i, const Coord& v) 
     108  { 
     109    for(int j=n; j > i ; j--) 
     110    { 
     111      vertex[j]=vertex[j-1] ; 
     112      edge[j]=edge[j-1] ; 
     113      d[j]=d[j-1] ; 
     114      neighbour[j]=neighbour[j-1] ; 
     115    } 
     116    vertex[i+1]=v ; 
     117    n++ ; 
     118  } 
     119   
    107120        int neighbour[NMAX]; 
    108121        double d[NMAX]; /**< distance of centre of small circle to origin, zero if great circle */ 
  • XIOS/trunk/extern/remap/src/intersect.cpp

    r688 r1016  
    3737} 
    3838 
     39/** New methods to find an insert a neighbour in a cell of the source mesh. 
     40 *  return true/false if cell b is a neighbour of a. if "insert" is true, then b will be inserted as a neighbour 
     41 * in cell a . This is needed for 2 order interpolation that need neighboround for gradient computing. 
     42 * A cell is a neighbour if : 
     43 *  - it shares 2 countiguous vertex (ie an edge) with a 
     44 *  - A vertex of b is located on one of an edge of a. 
     45 **/ 
     46bool insertNeighbour( Elt& a, const Elt& b, bool insert ) 
     47{ 
     48  // for now suppose pole -> Oz 
     49  Coord pole(0,0,1) ; 
     50  Coord O, Oa1, Oa2,Ob1,Ob2,V1,V2 ; 
     51  double da,db,alpha,alpha1,alpha2,delta ; 
     52   
     53     
     54  for (int i = 0; i < a.n; i++) 
     55  { 
     56    for (int j = 0; j < b.n; j++) 
     57    { 
     58// share a full edge ? be carefull at the orientation 
     59      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > 1e-10*1e-10 || 
     60             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10); 
     61      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-10*1e-10 && 
     62             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-10*1e-10) 
     63      { 
     64        if (insert) a.neighbour[i] = b.id.ind ; 
     65        return true; 
     66      } 
     67       
     68 
     69// 1 or 2 vertices of b is located on an edge of a 
     70       da=a.d[i] ; 
     71       if (scalarprod(a.edge[i], pole) < 0) da=-da ; 
     72       db=b.d[(j+b.n-1)%b.n] ; 
     73       if (scalarprod(b.edge[(j+b.n-1)%b.n], pole) < 0) db=-db ; 
     74       
     75      if ( fabs(da-db)<1e-10 )  
     76      { 
     77        O=pole*da ; 
     78        Oa1=a.vertex[i]-O ; 
     79        Oa2=a.vertex[(i+1)%a.n]-O ;  
     80        Ob1=b.vertex[j]-O ; 
     81        Ob2=b.vertex[(j+b.n-1)%b.n]-O ; 
     82        V1=crossprod(Oa1,Oa2) ; 
     83        V2=crossprod(Ob1,Ob2) ; 
     84        if (norm(crossprod(V1,V2))/(norm(V1)*norm(V2)) < 1e-10) 
     85        { 
     86          alpha = vectAngle(Oa1,Oa2,V1) ; 
     87          alpha1= vectAngle(Oa1,Ob1,V1) ; 
     88          alpha2= vectAngle(Oa1,Ob2,V1) ; 
     89          delta= alpha2-alpha1 ; 
     90          if (delta >= M_PI) delta=2*M_PI-delta ; 
     91          else if (delta <= -M_PI) delta=2*M_PI+delta ; 
     92           
     93          if (alpha >= 0) 
     94          { 
     95            if (alpha1 > 1e-10 && alpha1 < alpha-1e-10) 
     96            { 
     97              if (alpha2 > 1e-10 && alpha2 < alpha-1e-10) 
     98              { 
     99                assert(delta > 0) ; 
     100                if (insert) 
     101                { 
     102                // insert both 
     103                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     104                  a.insert_vertex(i,b.vertex[j]); 
     105                  a.neighbour[i+1] = b.id.ind ; 
     106                } 
     107                return true ; 
     108              } 
     109              else 
     110              { 
     111                assert( delta > 0 ) ; 
     112                if (insert) 
     113                { 
     114                //insert alpha1 
     115                  a.insert_vertex(i,b.vertex[j]); 
     116                  a.neighbour[i+1] = b.id.ind ; 
     117                } 
     118                return true ; 
     119              } 
     120            } 
     121            else if (alpha2 > 1e-10 && alpha2 < alpha-1e-10) 
     122            { 
     123              assert( delta > 0 ) ; 
     124              if (insert) 
     125              { 
     126              // insert alpha2 
     127                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     128                a.neighbour[i] = b.id.ind ; 
     129              } 
     130              return true ; 
     131            } 
     132            else 
     133            { 
     134              // nothing to do 
     135            }  
     136 
     137          } 
     138          else  // alpha < 0 
     139          { 
     140            if (alpha1 < -1e-10 && alpha1 > alpha+1e-10) 
     141            { 
     142              if (alpha2 < -1e-10 && alpha2 > alpha+1e-10) 
     143              { 
     144                assert(delta < 0) ; 
     145                if (insert) 
     146                { 
     147                // insert both 
     148                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     149                  a.insert_vertex(i,b.vertex[j]); 
     150                  a.neighbour[i+1] = b.id.ind ; 
     151                } 
     152                return true ; 
     153              } 
     154              else 
     155              { 
     156                assert(delta < 0) ; 
     157                if (insert) 
     158                { 
     159                //insert alpha1 
     160                  a.insert_vertex(i,b.vertex[j]); 
     161                  a.neighbour[i+1] = b.id.ind ; 
     162                } 
     163                return true ; 
     164              } 
     165            } 
     166            else if (alpha2 < -1e-10 && alpha2 > alpha+1e-10) 
     167            { 
     168              assert(delta < 0) ; 
     169              if (insert) 
     170              { 
     171              // insert alpha2 
     172                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]); 
     173                a.neighbour[i] = b.id.ind ; 
     174              } 
     175              return true ; 
     176            } 
     177            else 
     178            { 
     179               // nothing to do 
     180            } 
     181          } 
     182        } 
     183      } 
     184    } 
     185  } 
     186  return false; 
     187} 
    39188 
    40189/** 
     
    44193void set_neighbour(Elt& a, const Elt& b) 
    45194{ 
    46         if (b.id.ind == a.id.ind) return; 
    47         int idx = neighbour_idx(a, b); 
    48         if (idx != NOT_FOUND) 
    49                 a.neighbour[idx] = b.id.ind; 
     195  if (b.id.ind == a.id.ind) return; 
     196/* 
     197  int idx = neighbour_idx(a, b); 
     198  if (idx != NOT_FOUND) 
     199  a.neighbour[idx] = b.id.ind; 
     200*/ 
     201  insertNeighbour(a,b,true) ;  
    50202} 
    51203 
    52204/** return true if `a` and `b` share an edge */ 
    53 bool isNeighbour(const Elt& a, const Elt& b) 
    54 { 
    55         return neighbour_idx(a, b) != NOT_FOUND; 
     205bool isNeighbour(Elt& a, const Elt& b) 
     206{ 
     207        // return neighbour_idx(a, b) != NOT_FOUND; 
     208  return insertNeighbour(a,b,false) ; 
    56209} 
    57210 
  • XIOS/trunk/extern/remap/src/intersect.hpp

    r688 r1016  
    22 
    33void set_neighbour(Elt& elt, const Elt& elt2); 
    4 bool isNeighbour(const Elt& elt, const Elt& elt2); 
     4bool isNeighbour(Elt& elt, const Elt& elt2); 
    55 
    66void intersect(Elt *a, Elt *b); 
  • XIOS/trunk/extern/remap/src/triple.cpp

    r849 r1016  
    124124} 
    125125 
     126// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b 
     127double vectAngle(const Coord &a, const Coord &b, const Coord &pole) 
     128{  
     129  double nab = 1./(norm(a)*norm(b)) ; 
     130   
     131  Coord a_cross_b=crossprod(a, b)*nab ; 
     132  double sinVect ; 
     133  if (scalarprod(a_cross_b, pole) >= 0) sinVect=norm(a_cross_b) ; 
     134  else sinVect=-norm(a_cross_b) ; 
     135  double cosVect=scalarprod(a,b)*nab ; 
     136 
     137  return atan2(sinVect,cosVect) ; 
    126138} 
     139 
     140} 
  • XIOS/trunk/extern/remap/src/triple.hpp

    r688 r1016  
    121121 
    122122double angle(const Coord &a, const Coord &b, const Coord &pole); 
     123 
     124// return oriented vector angle in range [-pi..pi], pole must be orthogonal to a and b 
     125double vectAngle(const Coord &a, const Coord &b, const Coord &pole) ; 
     126 
    123127void print_Coord(Coord &p); 
    124128 
Note: See TracChangeset for help on using the changeset viewer.