source: XIOS/dev/dev_trunk_omp/extern/remap/src/polyg.cpp @ 1630

Last change on this file since 1630 was 1630, checked in by yushan, 4 years ago

working branch @1608 with bug fix @1628

File size: 9.4 KB
Line 
1/* utilities related to polygons */
2
3#include <cassert>
4#include <iostream>
5#include "elt.hpp"
6#include "errhandle.hpp"
7
8#include "polyg.hpp"
9
10namespace sphereRemap {
11
12using namespace std;
13
14/* given `N` `vertex`es, `N` `edge`s and `N` `d`s (d for small circles)
15   and `g` the barycentre,
16   this function reverses the order of arrays of `vertex`es, `edge`s and `d`s
17   but only if it is required!
18  (because computing intersection area requires both polygons to have same orientation)
19*/
20void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g)
21{
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  }
36}
37
38
39void normals(Coord *x, int n, Coord *a)
40{
41  for (int i = 0; i<n; i++)
42    a[i] = crossprod(x[(i+1)%n], x[i]);
43}
44
45Coord barycentre(const Coord *x, int n)
46{
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);
58}
59
60/** computes the barycentre of the area which is the difference
61  of a side ABO of the spherical tetrahedron and of the straight tetrahedron */
62static Coord tetrah_side_diff_centre(Coord a, Coord b)
63{
64  Coord n = crossprod(a,b);
65        double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z;
66  assert(sinc2 < 1.0 + EPS);
67
68        // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab
69  // approx:
70        // double u = sinc2/6. + (3./40.)*sinc2*sinc2;
71
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);
76        double u = asin(sinc)/sinc - 1;
77
78        return n*u;
79}
80
81/* compute the barycentre as the negative sum of all the barycentres of sides of the tetrahedron */
82Coord gc_normalintegral(const Coord *x, int n)
83{
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;
89}
90
91Coord exact_barycentre(const Coord *x, int n)
92{
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];
100}
101
102Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole)
103{
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  }
145}
146
147/** Computes area of two two-sided polygon
148   needs to have one small and one great circle, otherwise zero
149   (name origin: lun is moon in french)
150*/
151double alun(double b, double d)
152{
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));
160}
161
162/**
163  This function returns the area of a spherical element
164that can be composed of great and small circle arcs.
165 The caller must ensure this function is not called when `alun` should be called instaed.
166  This function also sets `gg` to the barycentre of the element.
167  "air" stands for area and  "bar" for barycentre.
168*/
169double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg)
170{
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];
184                double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
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;
204}
205
206double polygonarea(Coord *vertices, int N)
207{
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;
216}
217
218int packedPolygonSize(const Elt& e)
219{
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));
222}
223
224void packPolygon(const Elt& e, char *buffer, int& pos) 
225{
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  } 
248
249}
250
251void unpackPolygon(Elt& e, const char *buffer, int& pos) 
252{
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  }
275}
276
277int packIntersectionSize(const Elt& elt) 
278{
279  return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
280}
281
282void packIntersection(const Elt& e, char* buffer,int& pos) 
283{
284  for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
285  {
286    *((int *) &(buffer[0])) += 1;
287
288    *((int *) &(buffer[pos])) = e.id.ind;
289    pos += sizeof(int);
290
291    *((double *) &(buffer[pos])) = e.area;
292    pos += sizeof(double);
293
294    *((GloId *) &(buffer[pos])) = (*it)->id;
295    pos += sizeof(GloId);
296 
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  }
305}
306
307void unpackIntersection(Elt* e, const char* buffer) 
308{
309  int ind;
310  int pos = 0;
311 
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];
320
321    elt.area=*((double *) & (buffer[pos]));
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}
Note: See TracBrowser for help on using the repository browser.