Changeset 1590


Ignore:
Timestamp:
11/09/18 09:44:17 (5 years ago)
Author:
oabramkina
Message:

Backporting recent changes to dev before merging it to trunk.

Location:
XIOS/dev/dev_olga
Files:
2 added
12 edited

Legend:

Unmodified
Added
Removed

  • XIOS/dev/dev_olga/extern/remap/src/intersection_ym.cpp

    r849 r1590  
    88#include <stdlib.h> 
    99#include <limits> 
     10#include <array> 
     11#include <cstdint>  
     12#include "earcut.hpp" 
     13#include <fstream>  
     14 
    1015 
    1116#define epsilon 1e-3  // epsilon distance ratio over side lenght for approximate small circle by great circle 
    … …  
    1621using namespace std; 
    1722using namespace ClipperLib ; 
     23         
    1824 
    1925double intersect_ym(Elt *a, Elt *b) 
    2026{ 
     27 
     28  using N = uint32_t; 
     29  using Point = array<double, 2>; 
     30  vector<Point> vect_points; 
     31  vector< vector<Point> > polyline; 
    2132 
    2233// transform small circle into piece of great circle if necessary 
    … …  
    4556  double cos_alpha; 
    4657 
     58  /// vector<p2t::Point*> polyline; 
    4759  for(int n=0; n<na;n++) 
    4860  { 
    … …  
    5163    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 
    5264    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 
    53   } 
     65 
     66    vect_points.push_back( array<double, 2>() ); 
     67    vect_points[n][0] = a_gno[n].x; 
     68    vect_points[n][1] = a_gno[n].y; 
     69 
     70  } 
     71 
     72  polyline.push_back(vect_points); 
     73  vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 
     74   
     75  double area_a_gno=0 ; 
     76  for(int i=0;i<indices_a_gno.size()/3;++i) 
     77    { 
     78      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 
     79      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 
     80      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 
     81      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
     82    } 
     83 
     84  vect_points.clear(); 
     85  polyline.clear(); 
     86  indices_a_gno.clear(); 
     87 
     88  
    5489 
    5590  for(int n=0; n<nb;n++) 
    … …  
    5994    b_gno[n].y=scalarprod(srcPolygon[n],Oy)/cos_alpha ; 
    6095    b_gno[n].z=scalarprod(srcPolygon[n],Oz)/cos_alpha ; // must be equal to 1 
    61   } 
    62  
     96 
     97    vect_points.push_back( array<double, 2>() ); 
     98    vect_points[n][0] = b_gno[n].x; 
     99    vect_points[n][1] = b_gno[n].y; 
     100  } 
     101 
     102 
     103  polyline.push_back(vect_points); 
     104  vector<N> indices_b_gno = mapbox::earcut<N>(polyline); 
     105 
     106  double area_b_gno=0 ; 
     107  for(int i=0;i<indices_b_gno.size()/3;++i) 
     108    { 
     109      Coord x0 = Ox * polyline[0][indices_b_gno[3*i]][0] + Oy* polyline[0][indices_b_gno[3*i]][1] + Oz ; 
     110      Coord x1 = Ox * polyline[0][indices_b_gno[3*i+1]][0] + Oy* polyline[0][indices_b_gno[3*i+1]][1] + Oz ; 
     111      Coord x2 = Ox * polyline[0][indices_b_gno[3*i+2]][0] + Oy* polyline[0][indices_b_gno[3*i+2]][1] + Oz ; 
     112      area_b_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
     113    } 
     114 
     115  vect_points.clear(); 
     116  polyline.clear(); 
     117  indices_b_gno.clear(); 
    63118 
    64119 
    … …  
    109164  clip.AddPaths(dst, ptClip, true); 
    110165  clip.Execute(ctIntersection, intersection); 
    111  
     166  
    112167  double area=0 ; 
    113168 
    114169  for(int ni=0;ni<intersection.size(); ni++) 
    115170  { 
    116     // go back into real coordinate on the sphere 
    117171    Coord* intersectPolygon=new Coord[intersection[ni].size()] ; 
    118172    for(int n=0; n < intersection[ni].size(); n++) 
    119173    { 
    120       double x=intersection[ni][n].X/xscale+xoffset ; 
    121       double y=intersection[ni][n].Y/yscale+yoffset ; 
    122  
    123       intersectPolygon[n]=Ox*x+Oy*y+Oz ; 
    124       intersectPolygon[n]=intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 
    125     } 
    126  
    127 // remove redondants vertex 
    128     int nv=0 ; 
     174      intersectPolygon[n].x=intersection[ni][n].X/xscale+xoffset ; 
     175      intersectPolygon[n].y=intersection[ni][n].Y/yscale+yoffset ; 
     176    } 
     177     
     178 
     179    int nv=0; 
     180 
    129181    for(int n=0; n < intersection[ni].size(); n++) 
    130182    { 
    131       if (norm(intersectPolygon[n]-intersectPolygon[(n+1)%intersection[ni].size()])>fusion_vertex) 
     183       double dx=intersectPolygon[n].x-intersectPolygon[(n+1)%intersection[ni].size()].x ; 
     184       double dy=intersectPolygon[n].y-intersectPolygon[(n+1)%intersection[ni].size()].y ; 
     185      
     186       if (dx*dx+dy*dy>fusion_vertex*fusion_vertex) 
     187       { 
     188          intersectPolygon[nv]=intersectPolygon[n] ; 
     189 
     190          vect_points.push_back( array<double, 2>() ); 
     191          vect_points[nv][0] = intersectPolygon[n].x; 
     192          vect_points[nv][1] = intersectPolygon[n].y; 
     193 
     194          nv++ ; 
     195       } 
     196       
     197 
     198    } 
     199 
     200    polyline.push_back(vect_points); 
     201    vect_points.clear(); 
     202 
     203    if (nv>2) 
     204    { 
     205  
     206      vector<N> indices = mapbox::earcut<N>(polyline); 
     207 
     208      double area2=0 ; 
     209      for(int i=0;i<indices.size()/3;++i) 
    132210      { 
    133         intersectPolygon[nv]=intersectPolygon[n] ; 
    134         nv++ ; 
     211          Coord x0 = Ox * polyline[0][indices[3*i]][0] + Oy* polyline[0][indices[3*i]][1] + Oz ; 
     212          Coord x1 = Ox * polyline[0][indices[3*i+1]][0] + Oy* polyline[0][indices[3*i+1]][1] + Oz ; 
     213          Coord x2 = Ox * polyline[0][indices[3*i+2]][0] + Oy* polyline[0][indices[3*i+2]][1] + Oz ; 
     214          area2+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 
    135215      } 
    136     } 
    137  
    138  
    139     if (nv>2) 
    140     { 
     216 
     217      polyline.clear(); 
     218 
     219      for(int n=0; n < nv; n++) 
     220      { 
     221        intersectPolygon[n] = Ox*intersectPolygon[n].x+Oy*intersectPolygon[n].y+Oz; 
     222        intersectPolygon[n] = intersectPolygon[n]*(1./norm(intersectPolygon[n])) ; 
     223      } 
     224 
     225 
    141226//     assign intersection to source and destination polygons 
    142227       Polyg *is = new Polyg; 
    143228       is->x = exact_barycentre(intersectPolygon,nv); 
    144        is->area = polygonarea(intersectPolygon,nv) ; 
     229//       is->area = polygonarea(intersectPolygon,nv) ; 
     230       is->area = area2 ; 
     231 
    145232//        if (is->area < 1e-12) cout<<"Small intersection : "<<is->area<<endl ; 
    146233       if (is->area==0.) delete is ; 
    … …  
    161248  delete[] b_gno ; 
    162249  return area ; 
     250 
    163251} 
     252 
     253 
    164254 
    165255void createGreatCirclePolygon(const Elt& element, const Coord& pole, vector<Coord>& coordinates) 

  • XIOS/dev/dev_olga/extern/remap/src/meshutil.cpp

    r877 r1590  
    133133  { 
    134134    for (int i = 0; i < elts[j]->n; i++)  
    135       if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; 
     135      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour 
     136      else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell  
    136137      else  neighbours[i] = elts[elts[j]->neighbour[i]]; 
    137138 

  • XIOS/dev/dev_olga/extern/remap/src/polyg.cpp

    r950 r1590  
    2020void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g) 
    2121{ 
    22         Coord ga = vertex[0] - g; 
    23         Coord gb = vertex[1] - g; 
    24         Coord vertical = crossprod(ga, gb); 
    25         if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G 
    26         { 
    27                 for (int i = 0; i < N/2; i++) 
    28                         swap(vertex[i], vertex[N-1-i]); 
    29  
    30                 for (int i = 0; i < (N-1)/2; i++) 
    31                 { 
    32                         swap(edge[N-2-i], edge[i]); 
    33                         swap(d[i], d[N-2-i]); 
    34                 } 
    35         } 
     22  Coord ga = vertex[0] - g; 
     23  Coord gb = vertex[1] - g; 
     24  Coord vertical = crossprod(ga, gb); 
     25  if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G 
     26  { 
     27    for (int i = 0; i < N/2; i++) 
     28      swap(vertex[i], vertex[N-1-i]); 
     29 
     30    for (int i = 0; i < (N-1)/2; i++) 
     31    { 
     32      swap(edge[N-2-i], edge[i]); 
     33      swap(d[i], d[N-2-i]); 
     34    } 
     35  } 
    3636} 
    3737 
    … …  
    3939void normals(Coord *x, int n, Coord *a) 
    4040{ 
    41         for (int i = 0; i<n; i++) 
    42                 a[i] = crossprod(x[(i+1)%n], x[i]); 
     41  for (int i = 0; i<n; i++) 
     42    a[i] = crossprod(x[(i+1)%n], x[i]); 
    4343} 
    4444 
    4545Coord barycentre(const Coord *x, int n) 
    4646{ 
    47         if (n == 0) return ORIGIN; 
    48         Coord bc = ORIGIN; 
    49         for (int i = 0; i < n; i++) 
    50                 bc = bc + x[i]; 
    51         /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
    52            which can occur when weighted with tiny area */ 
    53  
    54         assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
    55         //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
    56  
    57         return proj(bc); 
     47  if (n == 0) return ORIGIN; 
     48  Coord bc = ORIGIN; 
     49  for (int i = 0; i < n; i++) 
     50    bc = bc + x[i]; 
     51  /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
     52     which can occur when weighted with tiny area */ 
     53 
     54  assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
     55  //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
     56 
     57  return proj(bc); 
    5858} 
    5959 
    … …  
    6262static Coord tetrah_side_diff_centre(Coord a, Coord b) 
    6363{ 
    64         Coord n = crossprod(a,b); 
     64  Coord n = crossprod(a,b); 
    6565        double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; 
    66         assert(sinc2 < 1.0 + EPS); 
     66  assert(sinc2 < 1.0 + EPS); 
    6767 
    6868        // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab 
    69         // approx: 
     69  // approx: 
    7070        // double u = sinc2/6. + (3./40.)*sinc2*sinc2; 
    7171 
    72         // exact   
    73         if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ 
    74                 return n * (M_PI_2 - 1); 
    75         double sinc = sqrt(sinc2); 
     72  // exact   
     73  if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ 
     74    return n * (M_PI_2 - 1); 
     75  double sinc = sqrt(sinc2); 
    7676        double u = asin(sinc)/sinc - 1; 
    7777 
    … …  
    8282Coord gc_normalintegral(const Coord *x, int n) 
    8383{ 
    84         Coord m = barycentre(x, n); 
    85         Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); 
    86         for (int i = 1; i < n; i++) 
    87                 bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); 
    88         return bc*0.5; 
     84  Coord m = barycentre(x, n); 
     85  Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); 
     86  for (int i = 1; i < n; i++) 
     87    bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); 
     88  return bc*0.5; 
    8989} 
    9090 
    9191Coord exact_barycentre(const Coord *x, int n) 
    9292{ 
    93         if (n >= 3) 
    94         { 
    95                 return  proj(gc_normalintegral(x, n)); 
    96         } 
    97         else if (n == 0) return ORIGIN; 
    98         else if (n == 2) return midpoint(x[0], x[1]); 
    99         else if (n == 1) return x[0]; 
     93  if (n >= 3) 
     94  { 
     95    return  proj(gc_normalintegral(x, n)); 
     96  } 
     97  else if (n == 0) return ORIGIN; 
     98  else if (n == 2) return midpoint(x[0], x[1]); 
     99  else if (n == 1) return x[0]; 
    100100} 
    101101 
    102102Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) 
    103103{ 
    104         double hemisphere = (a.z > 0) ? 1: -1; 
    105  
    106         double lat = hemisphere * (M_PI_2 - acos(a.z)); 
    107         double lon1 = atan2(a.y, a.x); 
    108         double lon2 = atan2(b.y, b.x); 
    109         double lon_diff = lon2 - lon1; 
    110  
    111         // wraparound at lon=-pi=pi 
    112         if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 
    113         else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 
    114  
    115         // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b  
    116         Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 
    117                                         0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 
    118                                         hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 
    119         Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 
    120         Coord t[] = {a, b, p}; 
    121         if (hemisphere < 0) swap(t[0], t[1]); 
    122         return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 
    123 } 
    124  
    125  
    126 double triarea(Coord& A, Coord& B, Coord& C) 
    127 { 
    128         double a = ds(B, C); 
    129         double b = ds(C, A); 
    130         double c = ds(A, B); 
    131         double s = 0.5 * (a + b + c); 
    132         double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); 
    133 //      assert(t >= 0); 
    134   if (t<1e-20) return 0. ; 
    135         return 4 * atan(sqrt(t)); 
     104  double hemisphere = (a.z > 0) ? 1: -1; 
     105 
     106  double lat = hemisphere * (M_PI_2 - acos(a.z)); 
     107  double lon1 = atan2(a.y, a.x); 
     108  double lon2 = atan2(b.y, b.x); 
     109  double lon_diff = lon2 - lon1; 
     110 
     111  // wraparound at lon=-pi=pi 
     112  if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 
     113  else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 
     114 
     115  // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b  
     116  Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 
     117                                  0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 
     118                                  hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 
     119  Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 
     120  Coord t[] = {a, b, p}; 
     121  if (hemisphere < 0) swap(t[0], t[1]); 
     122  return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 
     123} 
     124 
     125 
     126double triarea(const Coord& A, const Coord& B, const Coord& C) 
     127{ 
     128  double a = ds(B, C); 
     129  double b = ds(C, A); 
     130  double c = ds(A, B); 
     131  double tmp ; 
     132 
     133  if (a<b) { tmp=a ; a=b ; b=tmp ; } 
     134  if (c > a) { tmp=a ; a=c ; c=b, b=tmp;  } 
     135  else if (c > b) { tmp=c ; c=b ; b=tmp ; } 
     136   
     137  double s = 0.5 * (a + b + c); 
     138  double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c))); 
     139  if (t>0) return 4 * atan(sqrt(t)); 
     140  else 
     141  { 
     142    std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; 
     143    return 0 ; 
     144  } 
    136145} 
    137146 
    … …  
    142151double alun(double b, double d) 
    143152{ 
    144         double a  = acos(d); 
    145         assert(b <= 2 * a); 
    146         double s  = a + 0.5 * b; 
    147         double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); 
    148         double r  = sqrt(1 - d*d); 
    149         double p  = 2 * asin(sin(0.5*b) / r); 
    150         return p*(1 - d) - 4*atan(sqrt(t)); 
     153  double a  = acos(d); 
     154  assert(b <= 2 * a); 
     155  double s  = a + 0.5 * b; 
     156  double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); 
     157  double r  = sqrt(1 - d*d); 
     158  double p  = 2 * asin(sin(0.5*b) / r); 
     159  return p*(1 - d) - 4*atan(sqrt(t)); 
    151160} 
    152161 
    … …  
    160169double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) 
    161170{ 
    162         if (N < 3) 
    163                 return 0; /* polygons with less then three vertices have zero area */ 
    164         Coord t[3]; 
    165         t[0] = barycentre(x, N); 
    166         Coord *g = new Coord[N]; 
    167         double area = 0; 
    168         Coord gg_exact = gc_normalintegral(x, N); 
    169         for (int i = 0; i < N; i++) 
    170         { 
    171                 /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ 
    172                 int ii = (i + 1) % N; 
    173                 t[1] = x[i]; 
    174                 t[2] = x[ii]; 
     171  if (N < 3) 
     172    return 0; /* polygons with less then three vertices have zero area */ 
     173  Coord t[3]; 
     174  t[0] = barycentre(x, N); 
     175  Coord *g = new Coord[N]; 
     176  double area = 0; 
     177  Coord gg_exact = gc_normalintegral(x, N); 
     178  for (int i = 0; i < N; i++) 
     179  { 
     180    /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ 
     181    int ii = (i + 1) % N; 
     182    t[1] = x[i]; 
     183    t[2] = x[ii]; 
    175184                double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
    176                 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 
    177                 double area_gc = triarea(t[0], t[1], t[2]); 
    178                 double area_sc_gc_moon = 0; 
    179                 if (d[i]) /* handle small circle case */ 
    180                 { 
    181                         Coord m = midpoint(t[1], t[2]); 
    182                         double mext = scalarprod(m, c[i]) - d[i]; 
    183                         char sgl = (mext > 0) ? -1 : 1; 
    184                         area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 
    185                         gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 
    186                 } 
    187                 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 
    188                 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 
    189         } 
    190         gg = barycentre(g, N); 
    191         gg_exact = proj(gg_exact); 
    192         delete[] g; 
    193         gg = gg_exact; 
    194         return area; 
     185    assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 
     186    double area_gc = triarea(t[0], t[1], t[2]); 
     187    double area_sc_gc_moon = 0; 
     188    if (d[i]) /* handle small circle case */ 
     189    { 
     190      Coord m = midpoint(t[1], t[2]); 
     191      double mext = scalarprod(m, c[i]) - d[i]; 
     192      char sgl = (mext > 0) ? -1 : 1; 
     193      area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 
     194      gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 
     195    } 
     196    area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 
     197    g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 
     198  } 
     199  gg = barycentre(g, N); 
     200  gg_exact = proj(gg_exact); 
     201  delete[] g; 
     202  gg = gg_exact; 
     203  return area; 
    195204} 
    196205 
    197206double polygonarea(Coord *vertices, int N) 
    198207{ 
    199         assert(N >= 3); 
    200  
    201         /* compute polygon area as sum of triangles */ 
    202         Coord centre = barycentre(vertices, N); 
    203         double area = 0; 
    204         for (int i = 0; i < N; i++) 
    205                 area += triarea(centre, vertices[i], vertices[(i+1)%N]); 
    206         return area; 
     208  assert(N >= 3); 
     209 
     210  /* compute polygon area as sum of triangles */ 
     211  Coord centre = barycentre(vertices, N); 
     212  double area = 0; 
     213  for (int i = 0; i < N; i++) 
     214    area += triarea(centre, vertices[i], vertices[(i+1)%N]); 
     215  return area; 
    207216} 
    208217 
    209218int packedPolygonSize(const Elt& e) 
    210219{ 
    211         return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + 
    212                sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 
     220  return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + 
     221         sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 
    213222} 
    214223 
    215224void packPolygon(const Elt& e, char *buffer, int& pos)  
    216225{ 
    217         *((GloId *) &(buffer[pos])) = e.id; 
    218         pos += sizeof(e.id); 
    219         *((GloId *) &(buffer[pos])) = e.src_id; 
    220         pos += sizeof(e.src_id); 
    221  
    222         *((Coord *) &(buffer[pos])) = e.x; 
    223         pos += sizeof(e.x); 
    224  
    225         *((double*) &(buffer[pos])) = e.val; 
    226         pos += sizeof(e.val); 
    227  
    228         *((int *) & (buffer[pos])) = e.n; 
    229         pos += sizeof(e.n); 
    230  
    231         for (int i = 0; i < e.n; i++) 
    232         { 
    233                 *((double *) & (buffer[pos])) = e.d[i]; 
    234                 pos += sizeof(e.d[i]); 
    235  
    236                 *((Coord *) &(buffer[pos])) = e.vertex[i]; 
    237                 pos += sizeof(e.vertex[i]); 
    238         }  
     226  *((GloId *) &(buffer[pos])) = e.id; 
     227  pos += sizeof(e.id); 
     228  *((GloId *) &(buffer[pos])) = e.src_id; 
     229  pos += sizeof(e.src_id); 
     230 
     231  *((Coord *) &(buffer[pos])) = e.x; 
     232  pos += sizeof(e.x); 
     233 
     234  *((double*) &(buffer[pos])) = e.val; 
     235  pos += sizeof(e.val); 
     236 
     237  *((int *) & (buffer[pos])) = e.n; 
     238  pos += sizeof(e.n); 
     239 
     240  for (int i = 0; i < e.n; i++) 
     241  { 
     242    *((double *) & (buffer[pos])) = e.d[i]; 
     243    pos += sizeof(e.d[i]); 
     244 
     245    *((Coord *) &(buffer[pos])) = e.vertex[i]; 
     246    pos += sizeof(e.vertex[i]); 
     247  }  
    239248 
    240249} 
    … …  
    242251void unpackPolygon(Elt& e, const char *buffer, int& pos)  
    243252{ 
    244         e.id = *((GloId *) &(buffer[pos])); 
    245         pos += sizeof(e.id); 
    246         e.src_id = *((GloId *) &(buffer[pos])); 
    247         pos += sizeof(e.src_id); 
    248  
    249         e.x = *((Coord *) & (buffer[pos]) ); 
    250         pos += sizeof(e.x); 
    251  
    252         e.val = *((double *) & (buffer[pos])); 
    253         pos += sizeof(double); 
    254  
    255         e.n = *((int *) & (buffer[pos])); 
    256         pos += sizeof(int); 
    257  
    258         for (int i = 0; i < e.n; i++) 
    259         { 
    260                 e.d[i] = *((double *) & (buffer[pos])); 
    261                 pos += sizeof(double); 
    262  
    263                 e.vertex[i] = *((Coord *) & (buffer[pos])); 
    264                 pos += sizeof(Coord); 
    265         } 
     253  e.id = *((GloId *) &(buffer[pos])); 
     254  pos += sizeof(e.id); 
     255  e.src_id = *((GloId *) &(buffer[pos])); 
     256  pos += sizeof(e.src_id); 
     257 
     258  e.x = *((Coord *) & (buffer[pos]) ); 
     259  pos += sizeof(e.x); 
     260 
     261  e.val = *((double *) & (buffer[pos])); 
     262  pos += sizeof(double); 
     263 
     264  e.n = *((int *) & (buffer[pos])); 
     265  pos += sizeof(int); 
     266 
     267  for (int i = 0; i < e.n; i++) 
     268  { 
     269    e.d[i] = *((double *) & (buffer[pos])); 
     270    pos += sizeof(double); 
     271 
     272    e.vertex[i] = *((Coord *) & (buffer[pos])); 
     273    pos += sizeof(Coord); 
     274  } 
    266275} 
    267276 
    268277int packIntersectionSize(const Elt& elt)  
    269278{ 
    270         return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 
     279  return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 
    271280} 
    272281 
    … …  
    274283{ 
    275284  for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 
    276         { 
    277                 *((int *) &(buffer[0])) += 1; 
    278  
    279                 *((int *) &(buffer[pos])) = e.id.ind; 
    280                 pos += sizeof(int); 
     285  { 
     286    *((int *) &(buffer[0])) += 1; 
     287 
     288    *((int *) &(buffer[pos])) = e.id.ind; 
     289    pos += sizeof(int); 
    281290 
    282291    *((double *) &(buffer[pos])) = e.area; 
    283292    pos += sizeof(double); 
    284293 
    285                 *((GloId *) &(buffer[pos])) = (*it)->id; 
    286                 pos += sizeof(GloId); 
     294    *((GloId *) &(buffer[pos])) = (*it)->id; 
     295    pos += sizeof(GloId); 
    287296   
    288                 *((int *) &(buffer[pos])) = (*it)->n; 
    289                 pos += sizeof(int); 
    290                 *((double *) &(buffer[pos])) = (*it)->area; 
    291                 pos += sizeof(double); 
    292  
    293                 *((Coord *) &(buffer[pos])) = (*it)->x; 
    294                 pos += sizeof(Coord); 
    295         } 
     297    *((int *) &(buffer[pos])) = (*it)->n; 
     298    pos += sizeof(int); 
     299    *((double *) &(buffer[pos])) = (*it)->area; 
     300    pos += sizeof(double); 
     301 
     302    *((Coord *) &(buffer[pos])) = (*it)->x; 
     303    pos += sizeof(Coord); 
     304  } 
    296305} 
    297306 
    298307void unpackIntersection(Elt* e, const char* buffer)  
    299308{ 
    300         int ind; 
    301         int pos = 0; 
     309  int ind; 
     310  int pos = 0; 
    302311   
    303         int n = *((int *) & (buffer[pos])); 
    304         pos += sizeof(int); 
    305         for (int i = 0; i < n; i++) 
    306         { 
    307                 ind = *((int*) &(buffer[pos])); 
    308                 pos+=sizeof(int); 
    309  
    310                 Elt& elt= e[ind]; 
     312  int n = *((int *) & (buffer[pos])); 
     313  pos += sizeof(int); 
     314  for (int i = 0; i < n; i++) 
     315  { 
     316    ind = *((int*) &(buffer[pos])); 
     317    pos+=sizeof(int); 
     318 
     319    Elt& elt= e[ind]; 
    311320 
    312321    elt.area=*((double *) & (buffer[pos])); 
    313                 pos += sizeof(double); 
    314  
    315                 Polyg *polygon = new Polyg; 
    316  
    317                 polygon->id =  *((GloId *) & (buffer[pos])); 
    318                 pos += sizeof(GloId); 
    319  
    320                 polygon->n =  *((int *) & (buffer[pos])); 
    321                 pos += sizeof(int); 
    322  
    323                 polygon->area =  *((double *) & (buffer[pos])); 
    324                 pos += sizeof(double); 
    325  
    326                 polygon->x = *((Coord *) & (buffer[pos])); 
    327                 pos += sizeof(Coord); 
    328  
    329                 elt.is.push_back(polygon); 
    330         } 
    331 } 
    332  
    333 } 
     322    pos += sizeof(double); 
     323 
     324    Polyg *polygon = new Polyg; 
     325 
     326    polygon->id =  *((GloId *) & (buffer[pos])); 
     327    pos += sizeof(GloId); 
     328 
     329    polygon->n =  *((int *) & (buffer[pos])); 
     330    pos += sizeof(int); 
     331 
     332    polygon->area =  *((double *) & (buffer[pos])); 
     333    pos += sizeof(double); 
     334 
     335    polygon->x = *((Coord *) & (buffer[pos])); 
     336    pos += sizeof(Coord); 
     337 
     338    elt.is.push_back(polygon); 
     339  } 
     340} 
     341 
     342} 

  • XIOS/dev/dev_olga/extern/remap/src/polyg.hpp

    r688 r1590  
    2323void unpackIntersection(Elt *e, const char *buffer); 
    2424int packIntersectionSize(const Elt& e); 
     25double triarea( const Coord& A, const Coord& B, const Coord& C) ; 
    2526 
    2627} 

  • XIOS/dev/dev_olga/src/client.cpp

    r1542 r1590  
    1111#include "timer.hpp" 
    1212#include "buffer_client.hpp" 
     13#include "string_tools.hpp" 
    1314 
    1415namespace xios 
    … …  
    2324    StdOFStream CClient::m_infoStream; 
    2425    StdOFStream CClient::m_errorStream; 
    25  
    2626    MPI_Comm& CClient::getInterComm(void)   { return (interComm); } 
    27  
     27      
    2828///--------------------------------------------------------------- 
    2929/*! 
    … …  
    230230      } 
    231231    } 
     232 
     233/*! 
     234 * \fn void CClient::callOasisEnddef(void) 
     235 * \brief Send the order to the servers to call "oasis_enddef". It must be done by each compound of models before calling oasis_enddef on client side 
     236 * Function is only called by client. 
     237 */ 
     238    void CClient::callOasisEnddef(void) 
     239    { 
     240      bool oasisEnddef=CXios::getin<bool>("call_oasis_enddef",true) ; 
     241      if (!oasisEnddef) ERROR("void CClient::callOasisEnddef(void)", <<"Function xios_oasis_enddef called but variable <call_oasis_enddef> is set to false."<<endl 
     242                                                                     <<"Variable <call_oasis_enddef> must be set to true"<<endl) ; 
     243      if (CXios::isServer) 
     244      // Attached mode 
     245      { 
     246        // nothing to do    
     247      } 
     248      else 
     249      { 
     250        int rank ; 
     251        int msg=0 ; 
     252 
     253        MPI_Comm_rank(intraComm,&rank) ; 
     254        if (rank==0)  
     255        { 
     256          MPI_Send(&msg,1,MPI_INT,0,5,interComm) ; // tags oasis_endded = 5 
     257        } 
     258 
     259      } 
     260    } 
     261 
    232262 
    233263    void CClient::finalize(void) 

  • XIOS/dev/dev_olga/src/client.hpp

    r1243 r1590  
    1313        static void finalize(void); 
    1414        static void registerContext(const string& id, MPI_Comm contextComm); 
     15        static void callOasisEnddef(void) ; 
    1516 
    1617        static MPI_Comm intraComm; 

  • XIOS/dev/dev_olga/src/interface/c/icdata.cpp

    r1542 r1590  
    8181     CClient::registerContext(str, comm); 
    8282     CTimer::get("XIOS init context").suspend(); 
     83     CTimer::get("XIOS").suspend(); 
     84   } 
     85 
     86   void cxios_oasis_enddef() 
     87   { 
     88     CTimer::get("XIOS").resume(); 
     89     CClient::callOasisEnddef(); 
    8390     CTimer::get("XIOS").suspend(); 
    8491   } 

  • XIOS/dev/dev_olga/src/interface/fortran/idata.F90

    r1158 r1590  
    4242      END SUBROUTINE cxios_context_finalize 
    4343 
    44  
     44      SUBROUTINE  cxios_oasis_enddef() BIND(C) 
     45         USE ISO_C_BINDING 
     46      END SUBROUTINE cxios_oasis_enddef 
     47       
    4548      SUBROUTINE  cxios_finalize() BIND(C) 
    4649      END SUBROUTINE cxios_finalize 
    … …  
    507510    END SUBROUTINE  xios(finalize) 
    508511 
     512   SUBROUTINE  xios(oasis_enddef) 
     513   IMPLICIT NONE 
     514 
     515      CALL cxios_oasis_enddef 
     516 
     517    END SUBROUTINE  xios(oasis_enddef) 
    509518 
    510519   SUBROUTINE xios(close_context_definition)() 

  • XIOS/dev/dev_olga/src/interface/fortran/ixios.F90

    r981 r1590  
    1313 
    1414USE idata, ONLY : xios(initialize), xios(init_server), xios(finalize), xios(context_initialize), xios(context_is_initialized), & 
    15                   xios(close_context_definition), xios(context_finalize), xios(solve_inheritance) 
     15                  xios(close_context_definition), xios(context_finalize), xios(solve_inheritance), xios(oasis_enddef) 
    1616 
    1717USE idomain, ONLY : txios(domain), txios(domaingroup), xios(is_valid_domain), xios(is_valid_domaingroup) 

  • XIOS/dev/dev_olga/src/io/nc4_data_output.cpp

    r1589 r1590  
    5656      { 
    5757        StdString lonName,latName ; 
    58  
     58          
    5959        domain->computeWrittenIndex(); 
    6060        domain->computeWrittenCompressedIndex(comm_file); 

  • XIOS/dev/dev_olga/src/server.cpp

    r1542 r1590  
    1414#include "timer.hpp" 
    1515#include "event_scheduler.hpp" 
     16#include "string_tools.hpp" 
    1617 
    1718namespace xios 
    … …  
    343344 
    344345        string codesId=CXios::getin<string>("oasis_codes_id") ; 
    345         vector<string> splitted ; 
    346         boost::split( splitted, codesId, boost::is_any_of(","), boost::token_compress_on ) ; 
     346        vector<string> oasisCodeId=splitRegex(codesId,"\\s*,\\s*") ; 
     347  
    347348        vector<string>::iterator it ; 
    348349 
    … …  
    352353 
    353354//      (2) Create interComms with models 
    354         for(it=splitted.begin();it!=splitted.end();it++) 
     355        for(it=oasisCodeId.begin();it!=oasisCodeId.end();it++) 
    355356        { 
    356357          oasis_get_intercomm(newComm,*it) ; 
    … …  
    386387        } 
    387388        if (CXios::usingServer2) delete [] srvGlobalRanks ; 
    388         oasis_enddef() ; 
     389 
     390        bool oasisEnddef=CXios::getin<bool>("call_oasis_enddef",true) ; 
     391        if (!oasisEnddef) oasis_enddef() ; 
    389392      } 
    390393 
    … …  
    442445           listenContext(); 
    443446           listenRootContext(); 
     447           listenOasisEnddef() ; 
     448           listenRootOasisEnddef() ; 
    444449           if (!finished) listenFinalize() ; 
    445450         } 
    … …  
    447452         { 
    448453           listenRootContext(); 
     454           listenRootOasisEnddef() ; 
    449455           if (!finished) listenRootFinalize() ; 
    450456         } 
    … …  
    516522        } 
    517523      } 
     524 
     525 
     526   /*! 
     527    * Root process is listening for an order sent by client to call "oasis_enddef". 
     528    * The root client of a compound send the order (tag 5). It is probed and received. 
     529    * When the order has been received from each coumpound, the server root process ping the order to the root processes of the secondary levels of servers (if any). 
     530    * After, it also inform (asynchronous call) other processes of the communicator that the oasis_enddef call must be done 
     531    */ 
     532     
     533     void CServer::listenOasisEnddef(void) 
     534     { 
     535        int flag ; 
     536        MPI_Status status ; 
     537        list<MPI_Comm>::iterator it; 
     538        int msg ; 
     539        static int nbCompound=0 ; 
     540        int size ; 
     541        static bool sent=false ; 
     542        static MPI_Request* allRequests ; 
     543        static MPI_Status* allStatus ; 
     544 
     545 
     546        if (sent) 
     547        { 
     548          MPI_Comm_size(intraComm,&size) ; 
     549          MPI_Testall(size,allRequests, &flag, allStatus) ; 
     550          if (flag==true) 
     551          { 
     552            delete [] allRequests ; 
     553            delete [] allStatus ; 
     554            sent=false ; 
     555          } 
     556        } 
     557         
     558 
     559        for(it=interCommLeft.begin();it!=interCommLeft.end();it++) 
     560        { 
     561           MPI_Status status ; 
     562           traceOff() ; 
     563           MPI_Iprobe(0,5,*it,&flag,&status) ;  // tags oasis_endded = 5 
     564           traceOn() ; 
     565           if (flag==true) 
     566           { 
     567              MPI_Recv(&msg,1,MPI_INT,0,5,*it,&status) ; // tags oasis_endded = 5 
     568              nbCompound++ ; 
     569              if (nbCompound==interCommLeft.size()) 
     570              { 
     571                for (std::list<MPI_Comm>::iterator it = interCommRight.begin(); it != interCommRight.end(); it++) 
     572                { 
     573                   MPI_Send(&msg,1,MPI_INT,0,5,*it) ; // tags oasis_endded = 5 
     574                } 
     575                MPI_Comm_size(intraComm,&size) ; 
     576                allRequests= new MPI_Request[size] ; 
     577                allStatus= new MPI_Status[size] ; 
     578                for(int i=0;i<size;i++) MPI_Isend(&msg,1,MPI_INT,i,5,intraComm,&allRequests[i]) ; // tags oasis_endded = 5 
     579                sent=true ; 
     580              } 
     581           } 
     582        } 
     583     } 
     584      
     585   /*! 
     586    * Processes probes message from root process if oasis_enddef call must be done. 
     587    * When the order is received it is scheduled to be treated in a synchronized way by all server processes of the communicator 
     588    */ 
     589     void CServer::listenRootOasisEnddef(void) 
     590     { 
     591       int flag ; 
     592       MPI_Status status ; 
     593       const int root=0 ; 
     594       int msg ; 
     595       static bool eventSent=false ; 
     596 
     597       if (eventSent) 
     598       { 
     599         boost::hash<string> hashString; 
     600         size_t hashId = hashString("oasis_enddef"); 
     601         if (eventScheduler->queryEvent(0,hashId)) 
     602         { 
     603           oasis_enddef() ; 
     604           eventSent=false ; 
     605         } 
     606       } 
     607          
     608       traceOff() ; 
     609       MPI_Iprobe(root,5,intraComm, &flag, &status) ; 
     610       traceOn() ; 
     611       if (flag==true) 
     612       { 
     613         MPI_Recv(&msg,1,MPI_INT,root,5,intraComm,&status) ; // tags oasis_endded = 5 
     614         boost::hash<string> hashString; 
     615         size_t hashId = hashString("oasis_enddef"); 
     616         eventScheduler->registerEvent(0,hashId); 
     617         eventSent=true ; 
     618       } 
     619     } 
     620 
     621 
     622 
     623      
    518624 
    519625     void CServer::listenContext(void) 

  • XIOS/dev/dev_olga/src/server.hpp

    r1378 r1590  
    2222        static void listenRootContext(void); 
    2323        static void listenRootFinalize(void); 
     24        static void listenRootOasisEnddef(void); 
     25        static void listenOasisEnddef(void); 
    2426        static void registerContext(void* buff,int count, int leaderRank=0); 
    2527 
    … …  
    7274        static StdOFStream m_infoStream; 
    7375        static StdOFStream m_errorStream; 
    74  
    7576        static void openStream(const StdString& fileName, const StdString& ext, std::filebuf* fb); 
    7677    }; 
Note: See TracChangeset for help on using the changeset viewer.