1 | #include <list> |
---|
2 | #include "elt.hpp" |
---|
3 | #include "polyg.hpp" |
---|
4 | |
---|
5 | namespace sphereRemap { |
---|
6 | |
---|
7 | using namespace std; |
---|
8 | |
---|
9 | void 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 | |
---|
18 | void 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 */ |
---|
26 | void 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 | Coord gradient(Elt& elt, Elt **neighElts) |
---|
43 | { |
---|
44 | Coord grad = ORIGIN; |
---|
45 | Coord *neighBaryc = new Coord[elt.n]; |
---|
46 | for (int j = 0; j < elt.n; j++) |
---|
47 | { |
---|
48 | int k = (j + 1) % elt.n; |
---|
49 | neighBaryc[j] = neighElts[j]->x; |
---|
50 | Coord edgeNormal = crossprod(neighElts[k]->x, neighElts[j]->x); |
---|
51 | |
---|
52 | /* use nomenclauture form paper */ |
---|
53 | double f_i = elt.val; |
---|
54 | double f_j = neighElts[j]->val; |
---|
55 | double f_k = neighElts[k]->val; |
---|
56 | grad = grad + edgeNormal * (0.5*(f_j + f_k) - f_i); |
---|
57 | } |
---|
58 | /* area of the polygon whoes vertices are the barycentres the neighbours */ |
---|
59 | grad = grad * (1./polygonarea(neighBaryc, elt.n)); |
---|
60 | delete[] neighBaryc; |
---|
61 | |
---|
62 | return grad - elt.x * scalarprod(elt.x, grad); |
---|
63 | } |
---|
64 | |
---|
65 | void computeGradients(Elt **elts, int N) |
---|
66 | { |
---|
67 | |
---|
68 | for (int j = 0; j < N; j++) |
---|
69 | elts[j]->val = 0; |
---|
70 | |
---|
71 | Elt *neighbours[NMAX]; |
---|
72 | for (int j = 0; j < N; j++) |
---|
73 | { |
---|
74 | for (int i = 0; i < elts[j]->n; i++) |
---|
75 | neighbours[i] = elts[elts[j]->neighbour[i]]; |
---|
76 | |
---|
77 | for (int i = 0; i < elts[j]->n; i++) |
---|
78 | neighbours[i]->val = 0; |
---|
79 | for (int i = 0; i < elts[j]->n; i++) |
---|
80 | { |
---|
81 | elts[j]->neighId[i] = neighbours[i]->src_id; |
---|
82 | /* for weight computation all values are always kept zero and only set to one when used .. */ |
---|
83 | neighbours[i]->val = 1; |
---|
84 | elts[j]->gradNeigh[i] = gradient(*(elts[j]), neighbours); |
---|
85 | /* .. and right after zeroed again */ |
---|
86 | neighbours[i]->val = 0; |
---|
87 | } |
---|
88 | elts[j]->neighId[elts[j]->n].rank = -1; // mark end |
---|
89 | elts[j]->neighId[elts[j]->n].ind = -1; // mark end |
---|
90 | /* For the most naive algorithm the case where the element itself is one must also be considered. |
---|
91 | Thomas says this can later be optimized out. */ |
---|
92 | elts[j]->val = 1; |
---|
93 | elts[j]->grad = gradient(*(elts[j]), neighbours); |
---|
94 | elts[j]->val = 0; |
---|
95 | } |
---|
96 | } |
---|
97 | |
---|
98 | } |
---|