source: XIOS/dev/dev_olga/extern/remap/src/intersect.cpp @ 1158

Last change on this file since 1158 was 1158, checked in by oabramkina, 4 years ago

Two server levels: merging with trunk r1137.
There are bugs.

File size: 10.0 KB
Line 
1#include <cstdlib>
2#include <cmath>
3#include <cassert>
4#include <cstring>
5#include <iostream>
6#include <fstream>
7#include "node.hpp"
8#include "elt.hpp"
9#include "gridRemap.hpp"
10#include "inside.hpp"
11#include "polyg.hpp"
12#include "intersect.hpp"
13#include "intersection_ym.hpp"
14
15namespace sphereRemap {
16
17using namespace std;
18
19/** returns index of edge of a that is shared with b,
20    or NOT_FOUND if a and b do not share an edge */
21int neighbour_idx(const Elt& a, const Elt& b)
22{
23        for (int i = 0; i < a.n; i++)
24        {
25                for (int j = 0; j < b.n; j++)
26                {
27                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
28                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
29                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
30                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
31                        {
32                                return i;
33                        }
34                }
35        }
36        return NOT_FOUND;
37}
38
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}
188
189/**
190If `a` and `b` are neighbours (in the sense that they share an edge)
191then this information will be stored in `a` (but not in `b`)
192*/
193void set_neighbour(Elt& a, const Elt& b)
194{
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) ; 
202}
203
204/** return true if `a` and `b` share an edge */
205bool isNeighbour(Elt& a, const Elt& b)
206{
207        // return neighbour_idx(a, b) != NOT_FOUND;
208  return insertNeighbour(a,b,false) ;
209}
210
211/* computes intersection between elements a and b */
212
213void intersect(Elt *a, Elt *b)
214{
215        int na = a->n; /* vertices of a */
216        int nb = b->n; /* vertices of b */
217        Coord *c   = new Coord[na+nb];
218        Coord *c2  = new Coord[na+nb];
219        Coord *xc  = new Coord[na+nb];
220        Coord *xc2 = new Coord[na+nb];
221        Coord gc, gc2;
222        double *d = new double[na+nb];
223        double *d2 = new double[na+nb];
224        double are, are2;
225        Ipt ipt[NMAX*NMAX];
226        Ipt ipt2[NMAX*NMAX];
227        ptsec(a, b, ipt);
228        /* make ipt2 transpose of ipt */
229        for (int ii = 0; ii < na; ii++)
230                for (int jj = 0; jj < nb; jj++)
231                        ipt2[jj*na+ii] = ipt[ii*nb+jj];
232        list<Sgm> iscot;
233        recense(a, b, ipt, iscot, 0);
234        recense(b, a, ipt2, iscot, 1);
235
236        int nseg = iscot.size();
237        int nc = 0;
238        int nc2 = 0;
239        while (iscot.size() && nc < 2)
240                nc = assemble(iscot, c, d, xc);
241        while (iscot.size() && nc2 < 2)
242                nc2 = assemble(iscot, c2, d2, xc2);
243//      assert(nseg == nc + nc2 || nseg == 1); // unused segment
244
245        if (!(nseg == nc + nc2 || nseg == 1))
246        {
247
248
249//          cout<<a->x.x<<"  "<<a->x.y<<"  "<<a->x.z<<endl ;
250//          cout<<b->x.x<<"  "<<b->x.y<<"  "<<b->x.z<<endl ;
251//          cout<<" List of vertex from a and b"<<endl ;
252//          for(int i=0;i<na;i++) cout <<"polygonPoints.InsertPoint("<<i<<", "<<a->vertex[i].x<<", "<<a->vertex[i].y<<", "<<a->vertex[i].z<<")"<<endl ;
253//          for(int i=0;i<nb;i++) cout <<"polygonPoints.InsertPoint("<<i+6<<", "<<b->vertex[i].x<<", "<<b->vertex[i].y<<", "<<b->vertex[i].z<<")"<<endl ;
254//          cout<<"na : "<<na<<"   nb : "<<nb<<endl;
255//          cout<<"nc :"<<nc<<"  nc2 :"<<nc2<<"   nseg : "<<nseg<<endl ;
256//          abort() ;
257//         cout<<"**********************************************"<<endl ;
258//         intersect_ym(a,b) ;
259        }
260
261//    intersect_ym(a,b) ;
262        if (nc == 1) nc = 0;
263        if (nc2 == 1) nc2 = 0;
264        gc = barycentre(xc, nc);
265        gc2 = barycentre(xc2, nc2);
266        orient(nc, xc, c, d, gc);
267
268        Coord pole = srcGrid.pole;
269        if (pole == ORIGIN) pole = tgtGrid.pole;
270        const double MINBASE = 1e-11;
271        if (nc == 2) /* nc is the number of vertices of super mesh element */
272        {
273                double base = arcdist(xc[0], xc[1]);
274cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl;
275                gc = midpoint(gc, midpointSC(xc[0], xc[1]));
276                /* intersection area `are` must be zero here unless we have one great and one small circle */
277                are = alun(base, fabs(scalarprod(xc[0], pole)));
278        }
279        else
280        {
281                are = airbar(nc, xc, c, d, pole, gc);
282        }
283        if (nc2 == 2)
284        {
285                double base = arcdist(xc2[0], xc2[1]);
286cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl;
287                assert(base > MINBASE);
288                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
289                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
290        }
291        else
292        {
293                are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
294        }
295
296//  double ym_area=intersect_ym(a,b) ;
297
298  if (nc > 1)
299        {
300                /* create one super mesh polygon that src and dest point to */
301                Polyg *is = new Polyg;
302                is->x = gc;
303                is->area = are;
304                is->id = b->id;
305                is->src_id = b->src_id;
306                is->n = nc;
307                (a->is).push_back(is);
308                (b->is).push_back(is);
309/*
310          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
311    {
312      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
313      intersect_ym(a,b) ;
314    }
315*/
316//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
317        }
318        if (nc2 > 1)
319        {
320                Polyg *is = new Polyg;
321                is->x = gc2;
322                is->area = are2;
323                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
324                is->src_id = b->src_id;
325                is->n = nc2;
326                (a->is).push_back(is);
327                (b->is).push_back(is);
328/*
329    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
330    {
331      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
332      intersect_ym(a,b) ;
333    }
334*/
335//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
336        }
337/*
338  if (nc<=1 && nc2<=1)
339  {
340    if (ym_area>1e-12)
341    {
342      cout<<"Big area difference : "<<0<<"  "<<ym_area<<endl ;
343    }
344  }
345*/
346        delete [] c;
347        delete [] c2;
348        delete [] xc;
349        delete [] xc2;
350        delete [] d;
351        delete [] d2;
352}
353
354}
Note: See TracBrowser for help on using the repository browser.