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

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

bugfix : interpolation order 2 : spurious NaN values occur sometime due to 0 / 0.

YM

File size: 4.1 KB
Line 
1#include <list>
2#include "elt.hpp"
3#include "polyg.hpp"
4
5namespace sphereRemap {
6
7using namespace std;
8
9void cptEltGeom(Elt& elt, const Coord &pole)
10{
11  orient(elt.n, elt.vertex, elt.edge, elt.d, elt.x);
12  normals(elt, pole);
13  Coord gg;
14  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg);
15  elt.x = gg;
16}
17
18void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
19{
20  for (int ne=0; ne<N; ne++)
21    cptEltGeom(elt[ne], pole);
22}
23
24/* for all elements of size-N-array `elt`,
25   make centre areaweighted average centres of super mesh elements */
26void update_baryc(Elt *elt, int N)
27{
28  for (int ne=0; ne<N; ne++)
29  {
30    Elt &e = elt[ne];
31    int ns = e.is.size();  // sous-elements
32    Coord *sx = new Coord[ns];
33    int i=0;
34    for (list<Polyg*>::iterator it = e.is.begin(); it != e.is.end(); i++, it++)
35    {
36      sx[i] = (*it)->x * (*it)->area;
37    }
38    e.x = barycentre(sx, ns);
39  }
40}
41
42
43Coord gradient_old(Elt& elt, Elt **neighElts)
44{
45  Coord grad = ORIGIN;
46  Coord *neighBaryc = new Coord[elt.n];
47  for (int j = 0; j < elt.n; j++)
48  {
49    int k = (j + 1) % elt.n;
50    neighBaryc[j] = neighElts[j]->x;
51    Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
52
53    // use nomenclauture form paper
54    double f_i = elt.val;
55    double f_j = neighElts[j]->val;
56    double f_k = neighElts[k]->val;
57    grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i);
58  }
59  // area of the polygon whoes vertices are the barycentres the neighbours
60  grad = grad * (1./polygonarea(neighBaryc, elt.n));
61  delete[] neighBaryc;
62
63  return grad - elt.x * scalarprod(elt.x, grad);
64}
65
66
67
68Coord gradient(Elt& elt, Elt **neighElts)
69{
70   
71  Coord grad = ORIGIN;
72  Coord neighBaryc[3] ;
73
74  double f_i ;
75  double f_j ;
76  double f_k ;
77 
78  Coord edgeNormal ;
79  double area=0 ;
80  int k ;
81  int count=0 ;
82 
83  for (int j = 0; j < elt.n; j++)
84  {
85    f_i = elt.val;
86    k = (j + 1) % elt.n;
87    if (neighElts[j]==NULL || neighElts[k]==NULL) continue ;
88
89    // use nomenclauture form paper
90    f_j = neighElts[j]->val;
91    f_k = neighElts[k]->val;
92
93   
94   
95    neighBaryc[0] = elt.x;
96    neighBaryc[1] = neighElts[j]->x;
97    neighBaryc[2] = neighElts[k]->x;
98
99    edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x);
100    grad = grad + edgeNormal * (0.5*(f_k + f_j) - f_i);
101
102    edgeNormal = crossprod(neighElts[j]->x, elt.x);
103    grad = grad + edgeNormal * (0.5*(f_j + f_i) - f_i);
104
105    edgeNormal = crossprod(elt.x, neighElts[k]->x);
106    grad = grad + edgeNormal * (0.5*(f_i + f_k) - f_i);
107
108  // area of the polygon whoes vertices are the barycentres the neighbours
109
110    area+=polygonarea(neighBaryc, 3) ;
111    count++ ;
112
113  }
114  if (count>0) 
115  {
116    grad=grad*(1./area) ;
117    return grad - elt.x * scalarprod(elt.x, grad);
118  }
119  else return grad ;
120}
121
122
123
124
125void computeGradients(Elt **elts, int N)
126{
127 
128  for (int j = 0; j < N; j++)
129    elts[j]->val = 0;
130
131  Elt *neighbours[NMAX];
132  for (int j = 0; j < N; j++)
133  {
134    for (int i = 0; i < elts[j]->n; i++) 
135      if ( elts[j]->neighbour[i]== NOT_FOUND) neighbours[i]=NULL ;
136      else  neighbours[i] = elts[elts[j]->neighbour[i]];
137
138    for (int i = 0; i < elts[j]->n; i++)
139      if (neighbours[i]!=NULL) neighbours[i]->val = 0;
140     
141    for (int i = 0; i < elts[j]->n; i++)
142    {
143      if (neighbours[i]!=NULL)
144      {
145        elts[j]->neighId[i] = neighbours[i]->src_id;
146        /* for weight computation all values are always kept zero and only set to one when used .. */
147        neighbours[i]->val = 1;
148        elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours);
149        /* .. and right after zeroed again */
150        neighbours[i]->val = 0;
151      }
152      else
153      {
154        elts[j]->neighId[i].rank = -1; // mark end
155        elts[j]->neighId[i].ind = -1; // mark end
156      }
157    }
158
159    for(int i = elts[j]->n ; i < NMAX; i++)
160    {
161      elts[j]->neighId[i].rank = -1; // mark end
162      elts[j]->neighId[i].ind = -1; // mark end
163    }
164    /* For the most naive algorithm the case where the element itself is one must also be considered.
165       Thomas says this can later be optimized out. */
166    elts[j]->val = 1;
167    elts[j]->grad = gradient(*(elts[j]), neighbours);
168    elts[j]->val = 0;
169  }
170}
171
172}
Note: See TracBrowser for help on using the repository browser.