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