Changeset 1579


Ignore:
Timestamp:
09/19/18 14:22:53 (6 years ago)
Author:
ymipsl
Message:

XIOS remapper : improvement : Computation the area of spherical triangle : terms are now ordered to improve numerical precision for triangles with small angle

YM

Location:
XIOS/trunk/extern/remap/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • XIOS/trunk/extern/remap/src/polyg.cpp

    r950 r1579  
    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/trunk/extern/remap/src/polyg.hpp

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