source: XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp @ 1220

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

test_remap_omp tested on ADA except two fields

File size: 9.1 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#include <stdio.h>
10
11namespace sphereRemap {
12
13using namespace std;
14
15/* given `N` `vertex`es, `N` `edge`s and `N` `d`s (d for small circles)
16   and `g` the barycentre,
17   this function reverses the order of arrays of `vertex`es, `edge`s and `d`s
18   but only if it is required!
19  (because computing intersection area requires both polygons to have same orientation)
20*/
21void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g)
22{
23        Coord ga = vertex[0] - g;
24        Coord gb = vertex[1] - g;
25        Coord vertical = crossprod(ga, gb);
26        if (N > 2 && scalarprod(g, vertical) < 0)  // (GAxGB).G
27        {
28                for (int i = 0; i < N/2; i++)
29                        swap(vertex[i], vertex[N-1-i]);
30
31                for (int i = 0; i < (N-1)/2; i++)
32                {
33                        swap(edge[N-2-i], edge[i]);
34                        swap(d[i], d[N-2-i]);
35                }
36        }
37}
38
39
40void normals(Coord *x, int n, Coord *a)
41{
42        for (int i = 0; i<n; i++)
43                a[i] = crossprod(x[(i+1)%n], x[i]);
44}
45
46Coord barycentre(const Coord *x, int n)
47{
48
49        if (n == 0) return ORIGIN;
50        Coord bc = ORIGIN;
51        for (int i = 0; i < n; i++)
52        {
53                bc = bc + x[i];
54        }       
55        /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon
56           which can occur when weighted with tiny area */
57 
58       
59  //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));     
60
61        return proj(bc);
62}
63
64/** computes the barycentre of the area which is the difference
65  of a side ABO of the spherical tetrahedron and of the straight tetrahedron */
66static Coord tetrah_side_diff_centre(Coord a, Coord b)
67{
68        Coord n = crossprod(a,b);
69        double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z;
70        assert(sinc2 < 1.0 + EPS);
71
72        // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab
73        // approx:
74        // double u = sinc2/6. + (3./40.)*sinc2*sinc2;
75
76        // exact 
77        if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */
78                return n * (M_PI_2 - 1);
79        double sinc = sqrt(sinc2);
80        double u = asin(sinc)/sinc - 1;
81
82        return n*u;
83}
84
85/* compute the barycentre as the negative sum of all the barycentres of sides of the tetrahedron */
86Coord gc_normalintegral(const Coord *x, int n)
87{
88        Coord m = barycentre(x, n);
89        Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]);
90        for (int i = 1; i < n; i++)
91                bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]);
92        return bc*0.5;
93}
94
95Coord exact_barycentre(const Coord *x, int n)
96{
97        if (n >= 3)
98        {
99                return  proj(gc_normalintegral(x, n));
100        }
101        else if (n == 0) return ORIGIN;
102        else if (n == 2) return midpoint(x[0], x[1]);
103        else if (n == 1) return x[0];
104}
105
106Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole)
107{
108        double hemisphere = (a.z > 0) ? 1: -1;
109
110        double lat = hemisphere * (M_PI_2 - acos(a.z));
111        double lon1 = atan2(a.y, a.x);
112        double lon2 = atan2(b.y, b.x);
113        double lon_diff = lon2 - lon1;
114
115        // wraparound at lon=-pi=pi
116        if (lon_diff < -M_PI) lon_diff += 2.0*M_PI;
117        else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI;
118
119        // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b
120        Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
121                                        0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)),
122                                        hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0));
123        Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1)
124        Coord t[] = {a, b, p};
125        if (hemisphere < 0) swap(t[0], t[1]);
126        return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere;
127}
128
129
130double triarea(Coord& A, Coord& B, Coord& C)
131{
132        double a = ds(B, C);
133        double b = ds(C, A);
134        double c = ds(A, B);
135        double s = 0.5 * (a + b + c);
136        double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c));
137//      assert(t >= 0);
138  if (t<1e-20) return 0. ;
139        return 4 * atan(sqrt(t));
140}
141
142/** Computes area of two two-sided polygon
143   needs to have one small and one great circle, otherwise zero
144   (name origin: lun is moon in french)
145*/
146double alun(double b, double d)
147{
148        double a  = acos(d);
149        assert(b <= 2 * a);
150        double s  = a + 0.5 * b;
151        double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b));
152        double r  = sqrt(1 - d*d);
153        double p  = 2 * asin(sin(0.5*b) / r);
154        return p*(1 - d) - 4*atan(sqrt(t));
155}
156
157/**
158  This function returns the area of a spherical element
159that can be composed of great and small circle arcs.
160 The caller must ensure this function is not called when `alun` should be called instaed.
161  This function also sets `gg` to the barycentre of the element.
162  "air" stands for area and  "bar" for barycentre.
163*/
164double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg)
165{
166        if (N < 3)
167                return 0; /* polygons with less than three vertices have zero area */
168        Coord t[3];
169        t[0] = barycentre(x, N);
170        Coord *g = new Coord[N];
171        double area = 0;
172        Coord gg_exact = gc_normalintegral(x, N);
173        for (int i = 0; i < N; i++)
174        {
175                /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */
176                int ii = (i + 1) % N;
177                t[1] = x[i];
178                t[2] = x[ii];
179    double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ;
180                assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation)
181                double area_gc = triarea(t[0], t[1], t[2]);
182                if(area_gc<=0) printf("area_gc = %e\n", area_gc);
183                double area_sc_gc_moon = 0;
184                if (d[i]) /* handle small circle case */
185                {
186                        Coord m = midpoint(t[1], t[2]);
187                        double mext = scalarprod(m, c[i]) - d[i];
188                        char sgl = (mext > 0) ? -1 : 1;
189                        area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole)));
190                        //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon);
191                        gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole);
192                }
193                area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */
194                g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon);
195                //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc,  area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i]));
196        }
197        gg = barycentre(g, N);
198        gg_exact = proj(gg_exact);
199        delete[] g;
200        gg = gg_exact;
201        return area;
202}
203
204double polygonarea(Coord *vertices, int N)
205{
206        assert(N >= 3);
207
208        /* compute polygon area as sum of triangles */
209        Coord centre = barycentre(vertices, N);
210        double area = 0;
211        for (int i = 0; i < N; i++)
212                area += triarea(centre, vertices[i], vertices[(i+1)%N]);
213        return area;
214}
215
216int packedPolygonSize(const Elt& e)
217{
218        return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) +
219               sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord));
220}
221
222void packPolygon(const Elt& e, char *buffer, int& pos) 
223{
224        *((GloId *) &(buffer[pos])) = e.id;
225        pos += sizeof(e.id);
226        *((GloId *) &(buffer[pos])) = e.src_id;
227        pos += sizeof(e.src_id);
228
229        *((Coord *) &(buffer[pos])) = e.x;
230        pos += sizeof(e.x);
231
232        *((double*) &(buffer[pos])) = e.val;
233        pos += sizeof(e.val);
234
235        *((int *) & (buffer[pos])) = e.n;
236        pos += sizeof(e.n);
237
238        for (int i = 0; i < e.n; i++)
239        {
240                *((double *) & (buffer[pos])) = e.d[i];
241                pos += sizeof(e.d[i]);
242
243                *((Coord *) &(buffer[pos])) = e.vertex[i];
244                pos += sizeof(e.vertex[i]);
245        } 
246
247}
248
249void unpackPolygon(Elt& e, const char *buffer, int& pos) 
250{
251        e.id = *((GloId *) &(buffer[pos]));
252        pos += sizeof(e.id);
253        e.src_id = *((GloId *) &(buffer[pos]));
254        pos += sizeof(e.src_id);
255
256        e.x = *((Coord *) & (buffer[pos]) );
257        pos += sizeof(e.x);
258
259        e.val = *((double *) & (buffer[pos]));
260        pos += sizeof(double);
261
262        e.n = *((int *) & (buffer[pos]));
263        pos += sizeof(int);
264
265        for (int i = 0; i < e.n; i++)
266        {
267                e.d[i] = *((double *) & (buffer[pos]));
268                pos += sizeof(double);
269
270                e.vertex[i] = *((Coord *) & (buffer[pos]));
271                pos += sizeof(Coord);
272        }
273}
274
275int packIntersectionSize(const Elt& elt) 
276{
277        return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double));
278}
279
280void packIntersection(const Elt& e, char* buffer,int& pos) 
281{
282  for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it)
283        {
284                *((int *) &(buffer[0])) += 1;
285
286                *((int *) &(buffer[pos])) = e.id.ind;
287                pos += sizeof(int);
288
289    *((double *) &(buffer[pos])) = e.area;
290    pos += sizeof(double);
291
292                *((GloId *) &(buffer[pos])) = (*it)->id;
293                pos += sizeof(GloId);
294 
295                *((int *) &(buffer[pos])) = (*it)->n;
296                pos += sizeof(int);
297                *((double *) &(buffer[pos])) = (*it)->area;
298                pos += sizeof(double);
299
300                *((Coord *) &(buffer[pos])) = (*it)->x;
301                pos += sizeof(Coord);
302        }
303}
304
305void unpackIntersection(Elt* e, const char* buffer) 
306{
307        int ind;
308        int pos = 0;
309 
310        int n = *((int *) & (buffer[pos]));
311        pos += sizeof(int);
312        for (int i = 0; i < n; i++)
313        {
314                ind = *((int*) &(buffer[pos]));
315                pos+=sizeof(int);
316
317                Elt& elt= e[ind];
318
319    elt.area=*((double *) & (buffer[pos]));
320                pos += sizeof(double);
321
322                Polyg *polygon = new Polyg;
323
324                polygon->id =  *((GloId *) & (buffer[pos]));
325                pos += sizeof(GloId);
326
327                polygon->n =  *((int *) & (buffer[pos]));
328                pos += sizeof(int);
329
330                polygon->area =  *((double *) & (buffer[pos]));
331                pos += sizeof(double);
332
333                polygon->x = *((Coord *) & (buffer[pos]));
334                pos += sizeof(Coord);
335
336                elt.is.push_back(polygon);
337        }
338}
339
340}
Note: See TracBrowser for help on using the repository browser.