source: XIOS/dev/branch_yushan_merged/extern/remap/src/intersect.cpp @ 1155

Last change on this file since 1155 was 1155, checked in by yushan, 5 years ago

test_remap OK with openmp

File size: 10.8 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
19extern CRemapGrid srcGrid;
20#pragma omp threadprivate(srcGrid)
21
22extern CRemapGrid tgtGrid;
23#pragma omp threadprivate(tgtGrid)
24
25/** returns index of edge of a that is shared with b,
26    or NOT_FOUND if a and b do not share an edge */
27int neighbour_idx(const Elt& a, const Elt& b)
28{
29        for (int i = 0; i < a.n; i++)
30        {
31                for (int j = 0; j < b.n; j++)
32                {
33                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
34                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
35                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
36                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
37                        {
38                                return i;
39                        }
40                }
41        }
42        return NOT_FOUND;
43}
44
45/** New methods to find an insert a neighbour in a cell of the source mesh.
46 *  return true/false if cell b is a neighbour of a. if "insert" is true, then b will be inserted as a neighbour
47 * in cell a . This is needed for 2 order interpolation that need neighboround for gradient computing.
48 * A cell is a neighbour if :
49 *  - it shares 2 countiguous vertex (ie an edge) with a
50 *  - A vertex of b is located on one of an edge of a.
51 **/
52bool insertNeighbour( Elt& a, const Elt& b, bool insert )
53{
54  // for now suppose pole -> Oz
55  Coord pole(0,0,1) ;
56  Coord O, Oa1, Oa2,Ob1,Ob2,V1,V2 ;
57  double da,db,alpha,alpha1,alpha2,delta ;
58 
59   
60  for (int i = 0; i < a.n; i++)
61  {
62    for (int j = 0; j < b.n; j++)
63    {
64      // share a full edge ? be carefull at the orientation
65     
66      //if(squaredist(a.vertex[i], b.vertex[j]) > 1e-10*1e-10 ||
67      //   squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10 )
68      //{     
69        //printf("A : squaredist(a.vertex[%d], b.vertex[%d]) = %.10e  %d\n",
70        //        i, j, squaredist(a.vertex[i], b.vertex[j]),
71        //        squaredist(a.vertex[i], b.vertex[j]) > 1e-10*1e-10 ? true : false);
72        //printf("B : squaredist(a.vertex[%d], b.vertex[%d]) = %.10e  %d\n",
73        //        (i+1)%a.n, (j+1)%b.n, squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]),
74        //        squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10 ? true : false);
75
76        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > 1e-10*1e-10 ||
77               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > 1e-10*1e-10);
78       
79      //}
80
81      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-10*1e-10 &&
82             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-10*1e-10)
83      {
84        if (insert) a.neighbour[i] = b.id.ind ;
85        return true;
86      }
87     
88
89// 1 or 2 vertices of b is located on an edge of a
90       da=a.d[i] ;
91       if (scalarprod(a.edge[i], pole) < 0) da=-da ;
92       db=b.d[(j+b.n-1)%b.n] ;
93       if (scalarprod(b.edge[(j+b.n-1)%b.n], pole) < 0) db=-db ;
94     
95      if ( fabs(da-db)<1e-10 ) 
96      {
97        O=pole*da ;
98        Oa1=a.vertex[i]-O ;
99        Oa2=a.vertex[(i+1)%a.n]-O ; 
100        Ob1=b.vertex[j]-O ;
101        Ob2=b.vertex[(j+b.n-1)%b.n]-O ;
102        V1=crossprod(Oa1,Oa2) ;
103        V2=crossprod(Ob1,Ob2) ;
104        if (norm(crossprod(V1,V2))/(norm(V1)*norm(V2)) < 1e-10)
105        {
106          alpha = vectAngle(Oa1,Oa2,V1) ;
107          alpha1= vectAngle(Oa1,Ob1,V1) ;
108          alpha2= vectAngle(Oa1,Ob2,V1) ;
109          delta= alpha2-alpha1 ;
110          if (delta >= M_PI) delta=2*M_PI-delta ;
111          else if (delta <= -M_PI) delta=2*M_PI+delta ;
112         
113          if (alpha >= 0)
114          {
115            if (alpha1 > 1e-10 && alpha1 < alpha-1e-10)
116            {
117              if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
118              {
119                assert(delta > 0) ;
120                if (insert)
121                {
122                // insert both
123                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
124                  a.insert_vertex(i,b.vertex[j]);
125                  a.neighbour[i+1] = b.id.ind ;
126                }
127                return true ;
128              }
129              else
130              {
131                assert( delta > 0 ) ;
132                if (insert)
133                {
134                //insert alpha1
135                  a.insert_vertex(i,b.vertex[j]);
136                  a.neighbour[i+1] = b.id.ind ;
137                }
138                return true ;
139              }
140            }
141            else if (alpha2 > 1e-10 && alpha2 < alpha-1e-10)
142            {
143              assert( delta > 0 ) ;
144              if (insert)
145              {
146              // insert alpha2
147                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
148                a.neighbour[i] = b.id.ind ;
149              }
150              return true ;
151            }
152            else
153            {
154              // nothing to do
155            } 
156
157          }
158          else  // alpha < 0
159          {
160            if (alpha1 < -1e-10 && alpha1 > alpha+1e-10)
161            {
162              if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
163              {
164                assert(delta < 0) ;
165                if (insert)
166                {
167                // insert both
168                  a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
169                  a.insert_vertex(i,b.vertex[j]);
170                  a.neighbour[i+1] = b.id.ind ;
171                }
172                return true ;
173              }
174              else
175              {
176                assert(delta < 0) ;
177                if (insert)
178                {
179                //insert alpha1
180                  a.insert_vertex(i,b.vertex[j]);
181                  a.neighbour[i+1] = b.id.ind ;
182                }
183                return true ;
184              }
185            }
186            else if (alpha2 < -1e-10 && alpha2 > alpha+1e-10)
187            {
188              assert(delta < 0) ;
189              if (insert)
190              {
191              // insert alpha2
192                a.insert_vertex(i,b.vertex[(j+b.n-1)%b.n]);
193                a.neighbour[i] = b.id.ind ;
194              }
195              return true ;
196            }
197            else
198            {
199               // nothing to do
200            }
201          }
202        }
203      }
204    }
205  }
206  return false;
207}
208
209/**
210If `a` and `b` are neighbours (in the sense that they share an edge)
211then this information will be stored in `a` (but not in `b`)
212*/
213void set_neighbour(Elt& a, const Elt& b)
214{
215  if (b.id.ind == a.id.ind) return;
216/*
217  int idx = neighbour_idx(a, b);
218  if (idx != NOT_FOUND)
219  a.neighbour[idx] = b.id.ind;
220*/
221  insertNeighbour(a,b,true) ; 
222}
223
224/** return true if `a` and `b` share an edge */
225bool isNeighbour(Elt& a, const Elt& b)
226{
227        // return neighbour_idx(a, b) != NOT_FOUND;
228  return insertNeighbour(a,b,false) ;
229}
230
231/* computes intersection between elements a and b */
232
233void intersect(Elt *a, Elt *b)
234{
235        int na = a->n; /* vertices of a */
236        int nb = b->n; /* vertices of b */
237        Coord *c   = new Coord[na+nb];
238        Coord *c2  = new Coord[na+nb];
239        Coord *xc  = new Coord[na+nb];
240        Coord *xc2 = new Coord[na+nb];
241        Coord gc, gc2;
242        double *d = new double[na+nb];
243        double *d2 = new double[na+nb];
244        double are, are2;
245        Ipt ipt[NMAX*NMAX];
246        Ipt ipt2[NMAX*NMAX];
247        ptsec(a, b, ipt);
248        /* make ipt2 transpose of ipt */
249        for (int ii = 0; ii < na; ii++)
250                for (int jj = 0; jj < nb; jj++)
251                        ipt2[jj*na+ii] = ipt[ii*nb+jj];
252        list<Sgm> iscot;
253        recense(a, b, ipt, iscot, 0);
254        recense(b, a, ipt2, iscot, 1);
255
256        int nseg = iscot.size();
257        int nc = 0;
258        int nc2 = 0;
259        while (iscot.size() && nc < 2)
260                nc = assemble(iscot, c, d, xc);
261        while (iscot.size() && nc2 < 2)
262                nc2 = assemble(iscot, c2, d2, xc2);
263//      assert(nseg == nc + nc2 || nseg == 1); // unused segment
264
265        if (!(nseg == nc + nc2 || nseg == 1))
266        {
267
268
269//          cout<<a->x.x<<"  "<<a->x.y<<"  "<<a->x.z<<endl ;
270//          cout<<b->x.x<<"  "<<b->x.y<<"  "<<b->x.z<<endl ;
271//          cout<<" List of vertex from a and b"<<endl ;
272//          for(int i=0;i<na;i++) cout <<"polygonPoints.InsertPoint("<<i<<", "<<a->vertex[i].x<<", "<<a->vertex[i].y<<", "<<a->vertex[i].z<<")"<<endl ;
273//          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 ;
274//          cout<<"na : "<<na<<"   nb : "<<nb<<endl;
275//          cout<<"nc :"<<nc<<"  nc2 :"<<nc2<<"   nseg : "<<nseg<<endl ;
276//          abort() ;
277//         cout<<"**********************************************"<<endl ;
278//         intersect_ym(a,b) ;
279        }
280
281//    intersect_ym(a,b) ;
282        if (nc == 1) nc = 0;
283        if (nc2 == 1) nc2 = 0;
284        gc = barycentre(xc, nc);
285        gc2 = barycentre(xc2, nc2);
286        orient(nc, xc, c, d, gc);
287
288        Coord pole = srcGrid.pole;
289        if (pole == ORIGIN) pole = tgtGrid.pole;
290        const double MINBASE = 1e-11;
291        if (nc == 2) /* nc is the number of vertices of super mesh element */
292        {
293                double base = arcdist(xc[0], xc[1]);
294cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl;
295                gc = midpoint(gc, midpointSC(xc[0], xc[1]));
296                /* intersection area `are` must be zero here unless we have one great and one small circle */
297                are = alun(base, fabs(scalarprod(xc[0], pole)));
298        }
299        else
300        {
301                are = airbar(nc, xc, c, d, pole, gc);
302        }
303        if (nc2 == 2)
304        {
305                double base = arcdist(xc2[0], xc2[1]);
306cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl;
307                assert(base > MINBASE);
308                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
309                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
310        }
311        else
312        {
313                are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
314        }
315
316//  double ym_area=intersect_ym(a,b) ;
317
318  if (nc > 1)
319        {
320                /* create one super mesh polygon that src and dest point to */
321                Polyg *is = new Polyg;
322                is->x = gc;
323                is->area = are;
324                is->id = b->id;
325                is->src_id = b->src_id;
326                is->n = nc;
327                (a->is).push_back(is);
328                (b->is).push_back(is);
329/*
330          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
331    {
332      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
333      intersect_ym(a,b) ;
334    }
335*/
336//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
337        }
338        if (nc2 > 1)
339        {
340                Polyg *is = new Polyg;
341                is->x = gc2;
342                is->area = are2;
343                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
344                is->src_id = b->src_id;
345                is->n = nc2;
346                (a->is).push_back(is);
347                (b->is).push_back(is);
348/*
349    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
350    {
351      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
352      intersect_ym(a,b) ;
353    }
354*/
355//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
356        }
357/*
358  if (nc<=1 && nc2<=1)
359  {
360    if (ym_area>1e-12)
361    {
362      cout<<"Big area difference : "<<0<<"  "<<ym_area<<endl ;
363    }
364  }
365*/
366        delete [] c;
367        delete [] c2;
368        delete [] xc;
369        delete [] xc2;
370        delete [] d;
371        delete [] d2;
372}
373
374}
Note: See TracBrowser for help on using the repository browser.