1 | #include <list> |
---|
2 | #include "elt.hpp" |
---|
3 | #include "polyg.hpp" |
---|
4 | #include "intersection_ym.hpp" |
---|
5 | #include "earcut.hpp" |
---|
6 | #include <vector> |
---|
7 | |
---|
8 | namespace sphereRemap { |
---|
9 | |
---|
10 | using namespace std; |
---|
11 | |
---|
12 | double 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 | |
---|
66 | void 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 | |
---|
79 | void 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 */ |
---|
87 | void 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 | |
---|
104 | Coord 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 | |
---|
129 | Coord 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 | |
---|
186 | void 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 | } |
---|