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