source: XIOS/trunk/extern/remap/src/meshutil.cpp @ 1613

Last change on this file since 1613 was 1613, checked in by ymipsl, 5 years ago

interpolation : increase conservation accurency when interpolate integrated flux on the cell. Same methode si now used to compute cell area and intersected cell area.

YM

File size: 6.0 KB
Line 
1#include <list>
2#include "elt.hpp"
3#include "polyg.hpp"
4#include "intersection_ym.hpp"
5#include "earcut.hpp"
6#include <vector>
7
8namespace sphereRemap {
9
10using namespace std;
11
12double computePolygoneArea(Elt& a, const Coord &pole)
13{
14  using N = uint32_t;
15  using Point = array<double, 2>;
16  vector<Point> vect_points;
17  vector< vector<Point> > polyline;
18 
19  vector<Coord> dstPolygon ;
20  createGreatCirclePolygon(a, pole, dstPolygon) ;
21
22  int na=dstPolygon.size() ;
23  Coord *a_gno   = new Coord[na];
24 
25  Coord OC=barycentre(a.vertex,a.n) ;
26  Coord Oz=OC ;
27  Coord Ox=crossprod(Coord(0,0,1),Oz) ;
28// choose Ox not too small to avoid rounding error
29  if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ;
30  Ox=Ox*(1./norm(Ox)) ;
31  Coord Oy=crossprod(Oz,Ox) ;
32  double cos_alpha;
33
34  for(int n=0; n<na;n++)
35  {
36    cos_alpha=scalarprod(OC,dstPolygon[n]) ;
37    a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ;
38    a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ;
39    a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1
40
41    vect_points.push_back( array<double, 2>() );
42    vect_points[n][0] = a_gno[n].x;
43    vect_points[n][1] = a_gno[n].y;
44
45  }
46
47  polyline.push_back(vect_points);
48  vector<N> indices_a_gno = mapbox::earcut<N>(polyline);
49 
50  double area_a_gno=0 ;
51  for(int i=0;i<indices_a_gno.size()/3;++i)
52    {
53      Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ;
54      Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ;
55      Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ;
56      area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ;
57    }
58
59  vect_points.clear();
60  polyline.clear();
61  indices_a_gno.clear();
62  return area_a_gno ;
63}
64
65
66void cptEltGeom(Elt& elt, const Coord &pole)
67{
68  orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
69  normals(elt, pole);
70  Coord gg;
71  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
72  elt.x = gg;
73// overwrite area computation
74
75  elt.area =  computePolygoneArea(elt, pole) ;
76}
77
78
79void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
80{
81  for (int ne=0; ne<N; ne++)
82    cptEltGeom(elt[ne], pole);
83}
84
85/* for all elements of size-N-array `elt`,
86   make centre areaweighted average centres of super mesh elements */
87void update_baryc(Elt *elt, int N)
88{
89  for (int ne=0; ne<N; ne++)
90  {
91    Elt &e = elt[ne];
92    int ns = e.is.size();  // sous-elements
93    Coord *sx = new Coord[ns];
94    int i=0;
95    for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++)
96    {
97      sx[i] = (*it)->x * (*it)->area;
98    }
99    e.x = barycentre(sx, ns);
100  }
101}
102
103
104Coord gradient_old(Elt& elt, Elt **neighElts)
105{
106  Coord grad = ORIGIN;
107  Coord *neighBaryc = new Coord[elt.n];
108  for (int j = 0; j < elt.n; j++)
109  {
110    int k = (j + 1) % elt.n;
111    neighBaryc[j] = neighElts[j]->x;
112    Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
113
114    // use nomenclauture form paper
115    double f_i = elt.val;
116    double f_j = neighElts[j]->val;
117    double f_k = neighElts[k]->val;
118    grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i);
119  }
120  // area of the polygon whoes vertices are the barycentres the neighbours
121  grad = grad * (1./polygonarea(neighBaryc, elt.n));
122  delete[] neighBaryc;
123
124  return grad - elt.x * scalarprod(elt.x, grad);
125}
126
127
128
129Coord gradient(Elt& elt, Elt **neighElts)
130{
131   
132  Coord grad = ORIGIN;
133  Coord neighBaryc[3] ;
134
135  double f_i ;
136  double f_j ;
137  double f_k ;
138 
139  Coord edgeNormal ;
140  double area=0 ;
141  int k ;
142  int count=0 ;
143 
144  for (int j = 0; j < elt.n; j++)
145  {
146    f_i = elt.val;
147    k = (j + 1) % elt.n;
148    if (neighElts[j]==NULL || neighElts[k]==NULL) continue ;
149
150    // use nomenclauture form paper
151    f_j = neighElts[j]->val;
152    f_k = neighElts[k]->val;
153
154   
155   
156    neighBaryc[0] = elt.x;
157    neighBaryc[1] = neighElts[j]->x;
158    neighBaryc[2] = neighElts[k]->x;
159
160    edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
161    grad = grad + edgeNormal * (0.5*(f_k + f_j) - f_i);
162
163    edgeNormal = crossprod(neighElts[j]->x, elt.x);
164    grad = grad + edgeNormal * (0.5*(f_j + f_i) - f_i);
165
166    edgeNormal = crossprod(elt.x, neighElts[k]->x);
167    grad = grad + edgeNormal * (0.5*(f_i + f_k) - f_i);
168
169  // area of the polygon whoes vertices are the barycentres the neighbours
170
171    area+=polygonarea(neighBaryc, 3) ;
172    count++ ;
173
174  }
175  if (count>0) 
176  {
177    grad=grad*(1./area) ;
178    return grad - elt.x * scalarprod(elt.x, grad);
179  }
180  else return grad ;
181}
182
183
184
185
186void computeGradients(Elt **elts, int N)
187{
188 
189  for (int j = 0; j < N; j++)
190    elts[j]->val = 0;
191
192  Elt *neighbours[NMAX];
193  for (int j = 0; j < N; j++)
194  {
195    for (int i = 0; i < elts[j]->n; i++) 
196      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ; // no neighbour
197      else if (elts[elts[j]->neighbour[i]]->is.size() == 0) neighbours[i]=NULL ; // neighbour has none supermesh cell
198      else  neighbours[i] = elts[elts[j]->neighbour[i]];
199
200    for (int i = 0; i < elts[j]->n; i++)
201      if (neighbours[i]!=NULL) neighbours[i]->val = 0;
202     
203    for (int i = 0; i < elts[j]->n; i++)
204    {
205      if (neighbours[i]!=NULL)
206      {
207        elts[j]->neighId[i] = neighbours[i]->src_id;
208        /* for weight computation all values are always kept zero and only set to one when used .. */
209        neighbours[i]->val = 1;
210        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours);
211        /* .. and right after zeroed again */
212        neighbours[i]->val = 0;
213      }
214      else
215      {
216        elts[j]->neighId[i].rank = -1; // mark end
217        elts[j]->neighId[i].ind = -1; // mark end
218      }
219    }
220
221    for(int i = elts[j]->n ; i < NMAX; i++)
222    {
223      elts[j]->neighId[i].rank = -1; // mark end
224      elts[j]->neighId[i].ind = -1; // mark end
225    }
226    /* For the most naive algorithm the case where the element itself is one must also be considered.
227       Thomas says this can later be optimized out. */
228    elts[j]->val = 1;
229    elts[j]->grad = gradient(*(elts[j]), neighbours);
230    elts[j]->val = 0;
231  }
232}
233
234}
Note: See TracBrowser for help on using the repository browser.