Changeset 1220 for XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
- Timestamp:
- 07/20/17 09:18:34 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp
r1205 r1220 15 15 namespace sphereRemap { 16 16 17 extern CRemapGrid srcGrid; 18 #pragma omp threadprivate(srcGrid) 19 20 extern CRemapGrid tgtGrid; 21 #pragma omp threadprivate(tgtGrid) 22 17 23 using namespace std; 18 24 … … 21 27 int neighbour_idx(const Elt& a, const Elt& b) 22 28 { 23 24 25 26 27 28 29 30 31 32 33 34 35 36 29 for (int i = 0; i < a.n; i++) 30 { 31 for (int j = 0; j < b.n; j++) 32 { 33 assert(squaredist(a.vertex[ i ], b.vertex[ j ]) > EPS*EPS || 34 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 35 if ( squaredist(a.vertex[ i ], b.vertex[ j ]) < 1e-13*1e-13 && 36 squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 37 { 38 return i; 39 } 40 } 41 } 42 return NOT_FOUND; 37 43 } 38 44 … … 205 211 bool isNeighbour(Elt& a, const Elt& b) 206 212 { 207 213 // return neighbour_idx(a, b) != NOT_FOUND; 208 214 return insertNeighbour(a,b,false) ; 209 215 } … … 213 219 void intersect(Elt *a, Elt *b) 214 220 { 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 // 221 int na = a->n; /* vertices of a */ 222 int nb = b->n; /* vertices of b */ 223 Coord *c = new Coord[na+nb]; 224 Coord *c2 = new Coord[na+nb]; 225 Coord *xc = new Coord[na+nb]; 226 Coord *xc2 = new Coord[na+nb]; 227 Coord gc, gc2; 228 double *d = new double[na+nb]; 229 double *d2 = new double[na+nb]; 230 double are, are2; 231 Ipt ipt[NMAX*NMAX]; 232 Ipt ipt2[NMAX*NMAX]; 233 ptsec(a, b, ipt); 234 /* make ipt2 transpose of ipt */ 235 for (int ii = 0; ii < na; ii++) 236 for (int jj = 0; jj < nb; jj++) 237 ipt2[jj*na+ii] = ipt[ii*nb+jj]; 238 list<Sgm> iscot; 239 recense(a, b, ipt, iscot, 0); 240 recense(b, a, ipt2, iscot, 1); 241 242 int nseg = iscot.size(); 243 int nc = 0; 244 int nc2 = 0; 245 while (iscot.size() && nc < 2) 246 nc = assemble(iscot, c, d, xc); 247 while (iscot.size() && nc2 < 2) 248 nc2 = assemble(iscot, c2, d2, xc2); 249 // assert(nseg == nc + nc2 || nseg == 1); // unused segment 244 250 245 251 if (!(nseg == nc + nc2 || nseg == 1)) … … 260 266 261 267 // intersect_ym(a,b) ; 262 263 264 265 266 267 268 269 270 271 272 273 268 if (nc == 1) nc = 0; 269 if (nc2 == 1) nc2 = 0; 270 gc = barycentre(xc, nc); 271 gc2 = barycentre(xc2, nc2); 272 orient(nc, xc, c, d, gc); 273 274 Coord pole = srcGrid.pole; 275 if (pole == ORIGIN) pole = tgtGrid.pole; 276 const double MINBASE = 1e-11; 277 if (nc == 2) /* nc is the number of vertices of super mesh element */ 278 { 279 double base = arcdist(xc[0], xc[1]); 274 280 cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 275 276 277 278 279 280 281 282 283 284 285 281 gc = midpoint(gc, midpointSC(xc[0], xc[1])); 282 /* intersection area `are` must be zero here unless we have one great and one small circle */ 283 are = alun(base, fabs(scalarprod(xc[0], pole))); 284 } 285 else 286 { 287 are = airbar(nc, xc, c, d, pole, gc); 288 } 289 if (nc2 == 2) 290 { 291 double base = arcdist(xc2[0], xc2[1]); 286 292 cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 287 288 289 290 291 292 293 294 293 assert(base > MINBASE); 294 gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 295 are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 296 } 297 else 298 { 299 are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 300 } 295 301 296 302 // double ym_area=intersect_ym(a,b) ; 297 303 298 304 if (nc > 1) 299 300 301 302 303 304 305 306 307 308 305 { 306 /* create one super mesh polygon that src and dest point to */ 307 Polyg *is = new Polyg; 308 is->x = gc; 309 is->area = are; 310 is->id = b->id; 311 is->src_id = b->src_id; 312 is->n = nc; 313 (a->is).push_back(is); 314 (b->is).push_back(is); 309 315 /* 310 if (2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)316 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 311 317 { 312 318 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 314 320 } 315 321 */ 316 // 317 318 319 320 321 322 323 324 325 326 327 322 // cout<<"intersection : "<<are<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 323 } 324 if (nc2 > 1) 325 { 326 Polyg *is = new Polyg; 327 is->x = gc2; 328 is->area = are2; 329 is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 330 is->src_id = b->src_id; 331 is->n = nc2; 332 (a->is).push_back(is); 333 (b->is).push_back(is); 328 334 /* 329 if ( 335 if ( 2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 330 336 { 331 337 cout<<"Big area difference : "<<are<<" "<<ym_area<<endl ; … … 333 339 } 334 340 */ 335 // 336 341 // cout<<"intersection : "<<are2<<" "<< ym_area<<" diff : "<<fabs(are-ym_area)<<" ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 342 } 337 343 /* 338 344 if (nc<=1 && nc2<=1) … … 344 350 } 345 351 */ 346 347 348 349 350 351 352 } 353 354 } 352 delete [] c; 353 delete [] c2; 354 delete [] xc; 355 delete [] xc2; 356 delete [] d; 357 delete [] d2; 358 } 359 360 }
Note: See TracChangeset
for help on using the changeset viewer.