[688] | 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 | { |
---|
[844] | 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; |
---|
[688] | 16 | } |
---|
| 17 | |
---|
| 18 | void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) |
---|
| 19 | { |
---|
[844] | 20 | for (int ne=0; ne<N; ne++) |
---|
| 21 | cptEltGeom(elt[ne], pole); |
---|
[688] | 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 | { |
---|
[844] | 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 | } |
---|
[688] | 40 | } |
---|
| 41 | |
---|
[844] | 42 | |
---|
| 43 | Coord 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 | |
---|
[688] | 68 | Coord gradient(Elt& elt, Elt **neighElts) |
---|
| 69 | { |
---|
[844] | 70 | |
---|
| 71 | Coord grad = ORIGIN; |
---|
| 72 | Coord neighBaryc[3] ; |
---|
[688] | 73 | |
---|
[844] | 74 | double f_i ; |
---|
| 75 | double f_j ; |
---|
| 76 | double f_k ; |
---|
| 77 | |
---|
| 78 | Coord edgeNormal ; |
---|
| 79 | double area=0 ; |
---|
| 80 | int k ; |
---|
[877] | 81 | int count=0 ; |
---|
| 82 | |
---|
[844] | 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) ; |
---|
[877] | 111 | count++ ; |
---|
[844] | 112 | |
---|
| 113 | } |
---|
[877] | 114 | if (count>0) |
---|
| 115 | { |
---|
| 116 | grad=grad*(1./area) ; |
---|
| 117 | return grad - elt.x * scalarprod(elt.x, grad); |
---|
| 118 | } |
---|
| 119 | else return grad ; |
---|
[688] | 120 | } |
---|
| 121 | |
---|
[844] | 122 | |
---|
| 123 | |
---|
| 124 | |
---|
[688] | 125 | void computeGradients(Elt **elts, int N) |
---|
| 126 | { |
---|
[844] | 127 | |
---|
| 128 | for (int j = 0; j < N; j++) |
---|
| 129 | elts[j]->val = 0; |
---|
[688] | 130 | |
---|
[844] | 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]]; |
---|
[688] | 137 | |
---|
[844] | 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 | } |
---|
[688] | 170 | } |
---|
| 171 | |
---|
| 172 | } |
---|