Changeset 1642 for XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
- Timestamp:
- 01/23/19 10:31:44 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp
r1328 r1642 20 20 void orient(int N, Coord *vertex, Coord *edge, double *d, const Coord &g) 21 21 { 22 23 24 25 26 27 28 29 30 31 32 33 34 35 22 Coord ga = vertex[0] - g; 23 Coord gb = vertex[1] - g; 24 Coord vertical = crossprod(ga, gb); 25 if (N > 2 && scalarprod(g, vertical) < 0) // (GAxGB).G 26 { 27 for (int i = 0; i < N/2; i++) 28 swap(vertex[i], vertex[N-1-i]); 29 30 for (int i = 0; i < (N-1)/2; i++) 31 { 32 swap(edge[N-2-i], edge[i]); 33 swap(d[i], d[N-2-i]); 34 } 35 } 36 36 } 37 37 … … 39 39 void normals(Coord *x, int n, Coord *a) 40 40 { 41 42 41 for (int i = 0; i<n; i++) 42 a[i] = crossprod(x[(i+1)%n], x[i]); 43 43 } 44 44 45 45 Coord barycentre(const Coord *x, int n) 46 46 { 47 48 49 50 51 52 53 54 55 56 57 47 if (n == 0) return ORIGIN; 48 Coord bc = ORIGIN; 49 for (int i = 0; i < n; i++) 50 bc = bc + x[i]; 51 /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon 52 which can occur when weighted with tiny area */ 53 54 assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 55 //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 56 57 return proj(bc); 58 58 } 59 59 … … 62 62 static Coord tetrah_side_diff_centre(Coord a, Coord b) 63 63 { 64 64 Coord n = crossprod(a,b); 65 65 double sinc2 = n.x*n.x + n.y*n.y + n.z*n.z; 66 66 assert(sinc2 < 1.0 + EPS); 67 67 68 68 // exact: u = asin(sinc)/sinc - 1; asin(sinc) = geodesic length of arc ab 69 69 // approx: 70 70 // double u = sinc2/6. + (3./40.)*sinc2*sinc2; 71 71 72 73 74 75 72 // exact 73 if (sinc2 > 1.0 - EPS) /* if round-off leads to sinc > 1 asin produces NaN */ 74 return n * (M_PI_2 - 1); 75 double sinc = sqrt(sinc2); 76 76 double u = asin(sinc)/sinc - 1; 77 77 … … 82 82 Coord gc_normalintegral(const Coord *x, int n) 83 83 { 84 85 86 87 88 84 Coord m = barycentre(x, n); 85 Coord bc = crossprod(x[n-1]-m, x[0]-m) + tetrah_side_diff_centre(x[n-1], x[0]); 86 for (int i = 1; i < n; i++) 87 bc = bc + crossprod(x[i-1]-m, x[i]-m) + tetrah_side_diff_centre(x[i-1], x[i]); 88 return bc*0.5; 89 89 } 90 90 91 91 Coord exact_barycentre(const Coord *x, int n) 92 92 { 93 94 95 96 97 98 99 93 if (n >= 3) 94 { 95 return proj(gc_normalintegral(x, n)); 96 } 97 else if (n == 0) return ORIGIN; 98 else if (n == 2) return midpoint(x[0], x[1]); 99 else if (n == 1) return x[0]; 100 100 } 101 101 102 102 Coord sc_gc_moon_normalintegral(Coord a, Coord b, Coord pole) 103 103 { 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(Coord& A, Coord& B, Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double s = 0.5 * (a + b + c); 132 double t = tan(0.5*s) * tan(0.5*(s - a)) * tan(0.5*(s - b)) * tan(0.5*(s - c)); 133 // assert(t >= 0); 134 if (t<1e-20) return 0. ; 135 return 4 * atan(sqrt(t)); 104 double hemisphere = (a.z > 0) ? 1: -1; 105 106 double lat = hemisphere * (M_PI_2 - acos(a.z)); 107 double lon1 = atan2(a.y, a.x); 108 double lon2 = atan2(b.y, b.x); 109 double lon_diff = lon2 - lon1; 110 111 // wraparound at lon=-pi=pi 112 if (lon_diff < -M_PI) lon_diff += 2.0*M_PI; 113 else if (lon_diff > M_PI) lon_diff -= 2.0*M_PI; 114 115 // integral of the normal over the surface bound by great arcs a-pole and b-pole and small arc a-b 116 Coord sc_normalintegral = Coord(0.5*(sin(lon2)-sin(lon1))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 117 0.5*(cos(lon1)-cos(lon2))*(M_PI_2 - lat - 0.5*sin(2.0*lat)), 118 hemisphere * lon_diff * 0.25 * (cos(2.0*lat) + 1.0)); 119 Coord p = Coord(0,0,hemisphere); // TODO assumes north pole is (0,0,1) 120 Coord t[] = {a, b, p}; 121 if (hemisphere < 0) swap(t[0], t[1]); 122 return (sc_normalintegral - gc_normalintegral(t, 3)) * hemisphere; 123 } 124 125 126 double triarea(const Coord& A, const Coord& B, const Coord& C) 127 { 128 double a = ds(B, C); 129 double b = ds(C, A); 130 double c = ds(A, B); 131 double tmp ; 132 133 if (a<b) { tmp=a ; a=b ; b=tmp ; } 134 if (c > a) { tmp=a ; a=c ; c=b, b=tmp; } 135 else if (c > b) { tmp=c ; c=b ; b=tmp ; } 136 137 double s = 0.5 * (a + b + c); 138 double t = tan(0.25*(a+(b+c))) * tan(0.25*(c-(a-b))) * tan(0.25*( c + (a-b) )) * tan(0.25*( a + (b - c))); 139 if (t>0) return 4 * atan(sqrt(t)); 140 else 141 { 142 std::cout<<"double triarea(const Coord& A, const Coord& B, const Coord& C) : t < 0 !!! t="<<t<<endl ; 143 return 0 ; 144 } 136 145 } 137 146 … … 142 151 double alun(double b, double d) 143 152 { 144 145 146 147 148 149 150 153 double a = acos(d); 154 assert(b <= 2 * a); 155 double s = a + 0.5 * b; 156 double t = tan(0.5 * s) * tan(0.5 * (s - a)) * tan(0.5 * (s - a)) * tan(0.5 * (s - b)); 157 double r = sqrt(1 - d*d); 158 double p = 2 * asin(sin(0.5*b) / r); 159 return p*(1 - d) - 4*atan(sqrt(t)); 151 160 } 152 161 … … 160 169 double airbar(int N, const Coord *x, const Coord *c, double *d, const Coord& pole, Coord& gg) 161 170 { 162 163 return 0; /* polygons with less than three vertices have zero area */164 165 166 167 168 169 170 171 172 173 174 171 if (N < 3) 172 return 0; /* polygons with less then three vertices have zero area */ 173 Coord t[3]; 174 t[0] = barycentre(x, N); 175 Coord *g = new Coord[N]; 176 double area = 0; 177 Coord gg_exact = gc_normalintegral(x, N); 178 for (int i = 0; i < N; i++) 179 { 180 /* for "spherical circle segment" sum triangular part and "small moon" and => account for small circle */ 181 int ii = (i + 1) % N; 182 t[1] = x[i]; 183 t[2] = x[ii]; 175 184 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 185 assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 186 double area_gc = triarea(t[0], t[1], t[2]); 187 double area_sc_gc_moon = 0; 188 if (d[i]) /* handle small circle case */ 189 { 190 Coord m = midpoint(t[1], t[2]); 191 double mext = scalarprod(m, c[i]) - d[i]; 192 char sgl = (mext > 0) ? -1 : 1; 193 area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 194 gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 195 } 196 area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 197 g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 198 } 199 gg = barycentre(g, N); 200 gg_exact = proj(gg_exact); 201 delete[] g; 202 gg = gg_exact; 203 return area; 195 204 } 196 205 197 206 double polygonarea(Coord *vertices, int N) 198 207 { 199 200 201 202 203 204 205 206 208 assert(N >= 3); 209 210 /* compute polygon area as sum of triangles */ 211 Coord centre = barycentre(vertices, N); 212 double area = 0; 213 for (int i = 0; i < N; i++) 214 area += triarea(centre, vertices[i], vertices[(i+1)%N]); 215 return area; 207 216 } 208 217 209 218 int packedPolygonSize(const Elt& e) 210 219 { 211 return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val)+212 220 return sizeof(e.id) + sizeof(e.src_id) + sizeof(e.x) + sizeof(e.val) + sizeof(e.given_area)+ 221 sizeof(e.n) + e.n*(sizeof(double)+sizeof(Coord)); 213 222 } 214 223 215 224 void packPolygon(const Elt& e, char *buffer, int& pos) 216 225 { 217 *((GloId *) &(buffer[pos])) = e.id; 218 pos += sizeof(e.id); 219 *((GloId *) &(buffer[pos])) = e.src_id; 220 pos += sizeof(e.src_id); 221 222 *((Coord *) &(buffer[pos])) = e.x; 223 pos += sizeof(e.x); 224 225 *((double*) &(buffer[pos])) = e.val; 226 pos += sizeof(e.val); 227 228 *((int *) & (buffer[pos])) = e.n; 229 pos += sizeof(e.n); 230 231 for (int i = 0; i < e.n; i++) 232 { 233 *((double *) & (buffer[pos])) = e.d[i]; 234 pos += sizeof(e.d[i]); 235 236 *((Coord *) &(buffer[pos])) = e.vertex[i]; 237 pos += sizeof(e.vertex[i]); 238 } 226 *((GloId *) &(buffer[pos])) = e.id; 227 pos += sizeof(e.id); 228 *((GloId *) &(buffer[pos])) = e.src_id; 229 pos += sizeof(e.src_id); 230 231 *((Coord *) &(buffer[pos])) = e.x; 232 pos += sizeof(e.x); 233 234 *((double*) &(buffer[pos])) = e.val; 235 pos += sizeof(e.val); 236 237 *((double*) &(buffer[pos])) = e.given_area; 238 pos += sizeof(e.val); 239 240 *((int *) & (buffer[pos])) = e.n; 241 pos += sizeof(e.n); 242 243 for (int i = 0; i < e.n; i++) 244 { 245 *((double *) & (buffer[pos])) = e.d[i]; 246 pos += sizeof(e.d[i]); 247 248 *((Coord *) &(buffer[pos])) = e.vertex[i]; 249 pos += sizeof(e.vertex[i]); 250 } 239 251 240 252 } … … 242 254 void unpackPolygon(Elt& e, const char *buffer, int& pos) 243 255 { 244 e.id = *((GloId *) &(buffer[pos])); 245 pos += sizeof(e.id); 246 e.src_id = *((GloId *) &(buffer[pos])); 247 pos += sizeof(e.src_id); 248 249 e.x = *((Coord *) & (buffer[pos]) ); 250 pos += sizeof(e.x); 251 252 e.val = *((double *) & (buffer[pos])); 253 pos += sizeof(double); 254 255 e.n = *((int *) & (buffer[pos])); 256 pos += sizeof(int); 257 258 for (int i = 0; i < e.n; i++) 259 { 260 e.d[i] = *((double *) & (buffer[pos])); 261 pos += sizeof(double); 262 263 e.vertex[i] = *((Coord *) & (buffer[pos])); 264 pos += sizeof(Coord); 265 } 256 e.id = *((GloId *) &(buffer[pos])); 257 pos += sizeof(e.id); 258 e.src_id = *((GloId *) &(buffer[pos])); 259 pos += sizeof(e.src_id); 260 261 e.x = *((Coord *) & (buffer[pos]) ); 262 pos += sizeof(e.x); 263 264 e.val = *((double *) & (buffer[pos])); 265 pos += sizeof(double); 266 267 e.given_area = *((double *) & (buffer[pos])); 268 pos += sizeof(double); 269 270 e.n = *((int *) & (buffer[pos])); 271 pos += sizeof(int); 272 273 for (int i = 0; i < e.n; i++) 274 { 275 e.d[i] = *((double *) & (buffer[pos])); 276 pos += sizeof(double); 277 278 e.vertex[i] = *((Coord *) & (buffer[pos])); 279 pos += sizeof(Coord); 280 } 266 281 } 267 282 268 283 int packIntersectionSize(const Elt& elt) 269 284 { 270 285 return elt.is.size() * (2*sizeof(int)+ sizeof(GloId) + 5*sizeof(double)); 271 286 } 272 287 … … 274 289 { 275 290 for (list<Polyg *>::const_iterator it = e.is.begin(); it != e.is.end(); ++it) 276 277 278 279 280 291 { 292 *((int *) &(buffer[0])) += 1; 293 294 *((int *) &(buffer[pos])) = e.id.ind; 295 pos += sizeof(int); 281 296 282 297 *((double *) &(buffer[pos])) = e.area; 283 298 pos += sizeof(double); 284 299 285 *((GloId *) &(buffer[pos])) = (*it)->id; 286 pos += sizeof(GloId); 300 301 *((GloId *) &(buffer[pos])) = (*it)->id; 302 pos += sizeof(GloId); 287 303 288 289 290 291 292 293 294 295 304 *((int *) &(buffer[pos])) = (*it)->n; 305 pos += sizeof(int); 306 *((double *) &(buffer[pos])) = (*it)->area; 307 pos += sizeof(double); 308 309 *((Coord *) &(buffer[pos])) = (*it)->x; 310 pos += sizeof(Coord); 311 } 296 312 } 297 313 298 314 void unpackIntersection(Elt* e, const char* buffer) 299 315 { 300 301 316 int ind; 317 int pos = 0; 302 318 303 304 305 306 307 308 309 310 319 int n = *((int *) & (buffer[pos])); 320 pos += sizeof(int); 321 for (int i = 0; i < n; i++) 322 { 323 ind = *((int*) &(buffer[pos])); 324 pos+=sizeof(int); 325 326 Elt& elt= e[ind]; 311 327 312 328 elt.area=*((double *) & (buffer[pos])); 313 pos += sizeof(double); 314 315 Polyg *polygon = new Polyg; 316 317 polygon->id = *((GloId *) & (buffer[pos])); 318 pos += sizeof(GloId); 319 320 polygon->n = *((int *) & (buffer[pos])); 321 pos += sizeof(int); 322 323 polygon->area = *((double *) & (buffer[pos])); 324 pos += sizeof(double); 325 326 polygon->x = *((Coord *) & (buffer[pos])); 327 pos += sizeof(Coord); 328 329 elt.is.push_back(polygon); 330 } 331 } 332 333 } 329 pos += sizeof(double); 330 331 332 Polyg *polygon = new Polyg; 333 334 polygon->id = *((GloId *) & (buffer[pos])); 335 pos += sizeof(GloId); 336 337 polygon->n = *((int *) & (buffer[pos])); 338 pos += sizeof(int); 339 340 polygon->area = *((double *) & (buffer[pos])); 341 pos += sizeof(double); 342 343 polygon->x = *((Coord *) & (buffer[pos])); 344 pos += sizeof(Coord); 345 346 elt.is.push_back(polygon); 347 } 348 } 349 350 }
Note: See TracChangeset
for help on using the changeset viewer.