source: XIOS/trunk/extern/remap/src/intersect.cpp @ 688

Last change on this file since 688 was 688, checked in by mhnguyen, 6 years ago

Integrating remap library into XIOS

+) Change name of some files of remap library to be compatible with XIOS
+) Implement function to fill in automatically boundary longitude and latitude

Test
+) On Curie
+) test_remap correct

File size: 5.4 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
40/**
41If `a` and `b` are neighbours (in the sense that they share an edge)
42then this information will be stored in `a` (but not in `b`)
43*/
44void set_neighbour(Elt& a, const Elt& b)
45{
46        if (b.id.ind == a.id.ind) return;
47        int idx = neighbour_idx(a, b);
48        if (idx != NOT_FOUND)
49                a.neighbour[idx] = b.id.ind;
50}
51
52/** return true if `a` and `b` share an edge */
53bool isNeighbour(const Elt& a, const Elt& b)
54{
55        return neighbour_idx(a, b) != NOT_FOUND;
56}
57
58/* computes intersection between elements a and b */
59
60void intersect(Elt *a, Elt *b)
61{
62        int na = a->n; /* vertices of a */
63        int nb = b->n; /* vertices of b */
64        Coord *c   = new Coord[na+nb];
65        Coord *c2  = new Coord[na+nb];
66        Coord *xc  = new Coord[na+nb];
67        Coord *xc2 = new Coord[na+nb];
68        Coord gc, gc2;
69        double *d = new double[na+nb];
70        double *d2 = new double[na+nb];
71        double are, are2;
72        Ipt ipt[NMAX*NMAX];
73        Ipt ipt2[NMAX*NMAX];
74        ptsec(a, b, ipt);
75        /* make ipt2 transpose of ipt */
76        for (int ii = 0; ii < na; ii++)
77                for (int jj = 0; jj < nb; jj++)
78                        ipt2[jj*na+ii] = ipt[ii*nb+jj];
79        list<Sgm> iscot;
80        recense(a, b, ipt, iscot, 0);
81        recense(b, a, ipt2, iscot, 1);
82
83        int nseg = iscot.size();
84        int nc = 0;
85        int nc2 = 0;
86        while (iscot.size() && nc < 2)
87                nc = assemble(iscot, c, d, xc);
88        while (iscot.size() && nc2 < 2)
89                nc2 = assemble(iscot, c2, d2, xc2);
90//      assert(nseg == nc + nc2 || nseg == 1); // unused segment
91
92        if (!(nseg == nc + nc2 || nseg == 1))
93        {
94
95
96//          cout<<a->x.x<<"  "<<a->x.y<<"  "<<a->x.z<<endl ;
97//          cout<<b->x.x<<"  "<<b->x.y<<"  "<<b->x.z<<endl ;
98//          cout<<" List of vertex from a and b"<<endl ;
99//          for(int i=0;i<na;i++) cout <<"polygonPoints.InsertPoint("<<i<<", "<<a->vertex[i].x<<", "<<a->vertex[i].y<<", "<<a->vertex[i].z<<")"<<endl ;
100//          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 ;
101//          cout<<"na : "<<na<<"   nb : "<<nb<<endl;
102//          cout<<"nc :"<<nc<<"  nc2 :"<<nc2<<"   nseg : "<<nseg<<endl ;
103//          abort() ;
104//         cout<<"**********************************************"<<endl ;
105//         intersect_ym(a,b) ;
106        }
107
108//    intersect_ym(a,b) ;
109        if (nc == 1) nc = 0;
110        if (nc2 == 1) nc2 = 0;
111        gc = barycentre(xc, nc);
112        gc2 = barycentre(xc2, nc2);
113        orient(nc, xc, c, d, gc);
114
115        Coord pole = srcGrid.pole;
116        if (pole == ORIGIN) pole = tgtGrid.pole;
117        const double MINBASE = 1e-11;
118        if (nc == 2) /* nc is the number of vertices of super mesh element */
119        {
120                double base = arcdist(xc[0], xc[1]);
121cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl;
122                gc = midpoint(gc, midpointSC(xc[0], xc[1]));
123                /* intersection area `are` must be zero here unless we have one great and one small circle */
124                are = alun(base, fabs(scalarprod(xc[0], pole)));
125        }
126        else
127        {
128                are = airbar(nc, xc, c, d, pole, gc);
129        }
130        if (nc2 == 2)
131        {
132                double base = arcdist(xc2[0], xc2[1]);
133cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl;
134                assert(base > MINBASE);
135                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
136                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
137        }
138        else
139        {
140                are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
141        }
142
143//  double ym_area=intersect_ym(a,b) ;
144
145  if (nc > 1)
146        {
147                /* create one super mesh polygon that src and dest point to */
148                Polyg *is = new Polyg;
149                is->x = gc;
150                is->area = are;
151                is->id = b->id;
152                is->src_id = b->src_id;
153                is->n = nc;
154                (a->is).push_back(is);
155                (b->is).push_back(is);
156/*
157          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
158    {
159      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
160      intersect_ym(a,b) ;
161    }
162*/
163//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
164        }
165        if (nc2 > 1)
166        {
167                Polyg *is = new Polyg;
168                is->x = gc2;
169                is->area = are2;
170                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
171                is->src_id = b->src_id;
172                is->n = nc2;
173                (a->is).push_back(is);
174                (b->is).push_back(is);
175/*
176    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
177    {
178      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
179      intersect_ym(a,b) ;
180    }
181*/
182//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
183        }
184/*
185  if (nc<=1 && nc2<=1)
186  {
187    if (ym_area>1e-12)
188    {
189      cout<<"Big area difference : "<<0<<"  "<<ym_area<<endl ;
190    }
191  }
192*/
193        delete [] c;
194        delete [] c2;
195        delete [] xc;
196        delete [] xc2;
197        delete [] d;
198        delete [] d2;
199}
200
201}
Note: See TracBrowser for help on using the repository browser.