Branch data Line data Source code
1 : : /*
2 : : * IntxUtils.cpp
3 : : *
4 : : * Created on: Oct 3, 2012
5 : : */
6 : :
7 : : #ifdef WIN32 /* windows */
8 : : #define _USE_MATH_DEFINES // For M_PI
9 : : #endif
10 : : #include <cmath>
11 : : #include <cassert>
12 : : #include <iostream>
13 : : // this is for sstream
14 : : #include <sstream>
15 : :
16 : : #include "moab/IntxMesh/IntxUtils.hpp"
17 : : // this is from mbcoupler; maybe it should be moved somewhere in moab src
18 : : // right now, add a dependency to mbcoupler
19 : : // #include "ElemUtil.hpp"
20 : : #include "moab/MergeMesh.hpp"
21 : : #include "moab/ReadUtilIface.hpp"
22 : : #include "MBTagConventions.hpp"
23 : : #define CHECKNEGATIVEAREA
24 : :
25 : : #ifdef CHECKNEGATIVEAREA
26 : : #include <iomanip>
27 : : #endif
28 : : #include <queue>
29 : : #include <map>
30 : :
31 : : #ifdef MOAB_HAVE_TEMPESTREMAP
32 : : #include "GridElements.h"
33 : : #endif
34 : :
35 : : namespace moab
36 : : {
37 : :
38 : : #define CORRTAGNAME "__correspondent"
39 : : #define MAXEDGES 10
40 : :
41 : 9535 : int IntxUtils::borderPointsOfXinY2( double* X, int nX, double* Y, int nY, double* P, int* side, double epsilon_area )
42 : : {
43 : : // 2 triangles, 3 corners, is the corner of X in Y?
44 : : // Y must have a positive area
45 : : /*
46 : : */
47 : 9535 : int extraPoint = 0;
48 [ + + ]: 46957 : for( int i = 0; i < nX; i++ )
49 : : {
50 : : // compute double the area of all nY triangles formed by a side of Y and a corner of X; if
51 : : // one is negative, stop (negative means it is outside; X and Y are all oriented such that
52 : : // they are positive oriented;
53 : : // if one area is negative, it means it is outside the convex region, for sure)
54 : 37422 : double* A = X + 2 * i;
55 : :
56 : 37422 : int inside = 1;
57 [ + + ]: 107908 : for( int j = 0; j < nY; j++ )
58 : : {
59 : 98911 : double* B = Y + 2 * j;
60 : 98911 : int j1 = ( j + 1 ) % nY;
61 : 98911 : double* C = Y + 2 * j1; // no copy of data
62 : :
63 : 98911 : double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
64 [ + + ]: 98911 : if( area2 < -epsilon_area )
65 : : {
66 : 28425 : inside = 0;
67 : 28425 : break;
68 : : }
69 : : }
70 [ + + ]: 37422 : if( inside )
71 : : {
72 : 8997 : side[i] = 1; // so vertex i of X is inside the convex region formed by Y
73 : : // so side has nX dimension (first array)
74 : 8997 : P[extraPoint * 2] = A[0];
75 : 8997 : P[extraPoint * 2 + 1] = A[1];
76 : 8997 : extraPoint++;
77 : : }
78 : : }
79 : 9535 : return extraPoint;
80 : : }
81 : :
82 : : // used to order according to angle, so it can remove double easily
83 : : struct angleAndIndex
84 : : {
85 : : double angle;
86 : : int index;
87 : : };
88 : :
89 : 57400 : bool angleCompare( angleAndIndex lhs, angleAndIndex rhs )
90 : : {
91 : 57400 : return lhs.angle < rhs.angle;
92 : : }
93 : :
94 : : // nP might be modified too, we will remove duplicates if found
95 : 8562 : int IntxUtils::SortAndRemoveDoubles2( double* P, int& nP, double epsilon_1 )
96 : : {
97 [ + + ]: 8562 : if( nP < 2 ) return 0; // nothing to do
98 : :
99 : : // center of gravity for the points
100 : 8545 : double c[2] = { 0., 0. };
101 : 8545 : int k = 0;
102 [ + + ]: 43700 : for( k = 0; k < nP; k++ )
103 : : {
104 : 35155 : c[0] += P[2 * k];
105 : 35155 : c[1] += P[2 * k + 1];
106 : : }
107 : 8545 : c[0] /= nP;
108 : 8545 : c[1] /= nP;
109 : :
110 : : // how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
111 : : // intersection points
112 : : struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
113 : :
114 [ + + ]: 43700 : for( k = 0; k < nP; k++ )
115 : : {
116 : 35155 : double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
117 [ + + ][ + + ]: 35155 : if( x != 0. || y != 0. ) { pairAngleIndex[k].angle = atan2( y, x ); }
118 : : else
119 : : {
120 : 63 : pairAngleIndex[k].angle = 0;
121 : : // it would mean that the cells are touching at a vertex
122 : : }
123 : 35155 : pairAngleIndex[k].index = k;
124 : : }
125 : :
126 : : // this should be faster than the bubble sort we had before
127 [ + - ]: 8545 : std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
128 : : // copy now to a new double array
129 : : double PCopy[10 * MAXEDGES]; // the same dimension as P; very conservative, but faster than
130 : : // reallocate for a vector
131 [ + + ]: 43700 : for( k = 0; k < nP; k++ ) // index will show where it should go now;
132 : : {
133 : 35155 : int ck = pairAngleIndex[k].index;
134 : 35155 : PCopy[2 * k] = P[2 * ck];
135 : 35155 : PCopy[2 * k + 1] = P[2 * ck + 1];
136 : : }
137 : : // now copy from PCopy over original P location
138 [ + - ]: 8545 : std::copy( PCopy, PCopy + 2 * nP, P );
139 : :
140 : : // eliminate duplicates, finally
141 : :
142 : 8545 : int i = 0, j = 1; // the next one; j may advance faster than i
143 : : // check the unit
144 : : // double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
145 : : // than 1.e-5 cm
146 [ + + ]: 35155 : while( j < nP )
147 : : {
148 [ + - ]: 26610 : double d2 = dist2( &P[2 * i], &P[2 * j] );
149 [ + + ]: 26610 : if( d2 > epsilon_1 )
150 : : {
151 : 22906 : i++;
152 : 22906 : P[2 * i] = P[2 * j];
153 : 22906 : P[2 * i + 1] = P[2 * j + 1];
154 : : }
155 : 26610 : j++;
156 : : }
157 : : // test also the last point with the first one (index 0)
158 : : // the first one could be at -PI; last one could be at +PI, according to atan2 span
159 : :
160 [ + - ]: 8545 : double d2 = dist2( P, &P[2 * i] ); // check the first and last points (ordered from -pi to +pi)
161 [ + + ]: 8545 : if( d2 > epsilon_1 ) { nP = i + 1; }
162 : : else
163 : 326 : nP = i; // effectively delete the last point (that would have been the same with first)
164 [ + + ]: 8545 : if( nP == 0 ) nP = 1; // we should be left with at least one point we already tested if nP is 0 originally
165 : 8562 : return 0;
166 : : }
167 : :
168 : : // the marks will show what edges of blue intersect the red
169 : :
170 : 973 : ErrorCode IntxUtils::EdgeIntersections2( double* blue, int nsBlue, double* red, int nsRed, int* markb, int* markr,
171 : : double* points, int& nPoints )
172 : : {
173 : : /* EDGEINTERSECTIONS computes edge intersections of two elements
174 : : [P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
175 : : and blue ( stored column wise )
176 : : (point coordinates are stored column-wise, in counter clock
177 : : order) the points P where their edges intersect. In addition,
178 : : in n the indices of which neighbors of red are also intersecting
179 : : with blue are given.
180 : : */
181 : :
182 : : // points is an array with enough slots (24 * 2 doubles)
183 : 973 : nPoints = 0;
184 [ + + ]: 10703 : for( int i = 0; i < MAXEDGES; i++ )
185 : : {
186 : 9730 : markb[i] = markr[i] = 0;
187 : : }
188 : :
189 [ + + ]: 4615 : for( int i = 0; i < nsBlue; i++ )
190 : : {
191 [ + + ]: 17460 : for( int j = 0; j < nsRed; j++ )
192 : : {
193 : : double b[2];
194 : : double a[2][2]; // 2*2
195 : 13818 : int iPlus1 = ( i + 1 ) % nsBlue;
196 : 13818 : int jPlus1 = ( j + 1 ) % nsRed;
197 [ + + ]: 41454 : for( int k = 0; k < 2; k++ )
198 : : {
199 : 27636 : b[k] = red[2 * j + k] - blue[2 * i + k];
200 : : // row k of a: a(k, 0), a(k, 1)
201 : 27636 : a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
202 : 27636 : a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
203 : : }
204 : 13818 : double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
205 [ + - ]: 13818 : if( fabs( delta ) > 1.e-14 ) // this is close to machine epsilon
206 : : {
207 : : // not parallel
208 : 13818 : double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
209 : 13818 : double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
210 [ + + ][ + + ]: 13818 : if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
[ + + ][ + + ]
211 : : {
212 : : // the intersection is good
213 [ + + ]: 6048 : for( int k = 0; k < 2; k++ )
214 : : {
215 : 4032 : points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
216 : : }
217 : 2016 : markb[i] = 1; // so neighbor number i of blue will be considered too.
218 : 2016 : markr[j] = 1; // this will be used in advancing red around blue quad
219 : 13818 : nPoints++;
220 : : }
221 : : }
222 : : // the case delta ~ 0. will be considered by the interior points logic
223 : : }
224 : : }
225 : 973 : return MB_SUCCESS;
226 : : }
227 : :
228 : : // special one, for intersection between rll (constant latitude) and cs quads
229 : 7589 : ErrorCode IntxUtils::EdgeIntxRllCs( double* blue, CartVect* bluec, int* blueEdgeType, int nsBlue, double* red,
230 : : CartVect* redc, int nsRed, int* markb, int* markr, int plane, double R,
231 : : double* points, int& nPoints )
232 : : {
233 : : // if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
234 : : // everything else the same (except if there are 2 points resulting, which is rare)
235 [ + + ]: 37945 : for( int i = 0; i < 4; i++ )
236 : : { // always at most 4 , so maybe don't bother
237 : 30356 : markb[i] = markr[i] = 0;
238 : : }
239 : :
240 [ + + ]: 37727 : for( int i = 0; i < nsBlue; i++ )
241 : : {
242 : 30138 : int iPlus1 = ( i + 1 ) % nsBlue;
243 [ + + ]: 30138 : if( blueEdgeType[i] == 0 ) // old style, just 2d
244 : : {
245 [ + + ]: 75890 : for( int j = 0; j < nsRed; j++ )
246 : : {
247 : : double b[2];
248 : : double a[2][2]; // 2*2
249 : :
250 : 60712 : int jPlus1 = ( j + 1 ) % nsRed;
251 [ + + ]: 182136 : for( int k = 0; k < 2; k++ )
252 : : {
253 : 121424 : b[k] = red[2 * j + k] - blue[2 * i + k];
254 : : // row k of a: a(k, 0), a(k, 1)
255 : 121424 : a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
256 : 121424 : a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
257 : : }
258 : 60712 : double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
259 [ + + ]: 60712 : if( fabs( delta ) > 1.e-14 )
260 : : {
261 : : // not parallel
262 : 41362 : double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
263 : 41362 : double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
264 [ + + ][ + + ]: 41362 : if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
[ + + ][ + + ]
265 : : {
266 : : // the intersection is good
267 [ + + ]: 22143 : for( int k = 0; k < 2; k++ )
268 : : {
269 : 14762 : points[2 * nPoints + k] =
270 : 14762 : blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
271 : : }
272 : 7381 : markb[i] = 1; // so neighbor number i of blue will be considered too.
273 : 7381 : markr[j] = 1; // this will be used in advancing red around blue quad
274 : 41362 : nPoints++;
275 : : }
276 : : } // if the edges are too "parallel", skip them
277 : : }
278 : : }
279 : : else // edge type is 1, so use 3d intersection
280 : : {
281 : 14960 : CartVect& C = bluec[i];
282 : 14960 : CartVect& D = bluec[iPlus1];
283 [ + + ]: 74800 : for( int j = 0; j < nsRed; j++ )
284 : : {
285 : 59840 : int jPlus1 = ( j + 1 ) % nsRed; // nsRed is just 4, forget about it, usually
286 : 59840 : CartVect& A = redc[j];
287 : 59840 : CartVect& B = redc[jPlus1];
288 : 59840 : int np = 0;
289 : : double E[9];
290 [ + - ][ + - ]: 59840 : intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
[ + - ][ + - ]
[ + - ]
291 [ + + ]: 59840 : if( np == 0 ) continue;
292 [ - + ][ # # ]: 7734 : if( np >= 2 ) { std::cout << "intersection with 2 points :" << A << B << C << D << "\n"; }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
293 [ + + ]: 15468 : for( int k = 0; k < np; k++ )
294 : : {
295 : 7734 : gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
296 [ + - ][ + - ]: 7734 : points[2 * nPoints + 1] );
297 : 7734 : nPoints++;
298 : : }
299 : 7734 : markb[i] = 1; // so neighbor number i of blue will be considered too.
300 : 7734 : markr[j] = 1; // this will be used in advancing red around blue quad
301 : : }
302 : : }
303 : : }
304 : 7589 : return MB_SUCCESS;
305 : : }
306 : :
307 : : // vec utils related to gnomonic projection on a sphere
308 : :
309 : : // vec utils
310 : :
311 : : /*
312 : : *
313 : : * position on a sphere of radius R
314 : : * if plane specified, use it; if not, return the plane, and the point in the plane
315 : : * there are 6 planes, numbered 1 to 6
316 : : * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
317 : : *
318 : : * projection on the plane will preserve the orientation, such that a triangle, quad pointing
319 : : * outside the sphere will have a positive orientation in the projection plane
320 : : */
321 : 3446 : void IntxUtils::decide_gnomonic_plane( const CartVect& pos, int& plane )
322 : : {
323 : : // decide plane, based on max x, y, z
324 [ + + ]: 3446 : if( fabs( pos[0] ) < fabs( pos[1] ) )
325 : : {
326 [ + + ]: 1645 : if( fabs( pos[2] ) < fabs( pos[1] ) )
327 : : {
328 : : // pos[1] is biggest
329 [ + + ]: 1125 : if( pos[1] > 0 )
330 : 558 : plane = 2;
331 : : else
332 : 1125 : plane = 4;
333 : : }
334 : : else
335 : : {
336 : : // pos[2] is biggest
337 [ + + ]: 520 : if( pos[2] < 0 )
338 : 262 : plane = 5;
339 : : else
340 : 1645 : plane = 6;
341 : : }
342 : : }
343 : : else
344 : : {
345 [ + + ]: 1801 : if( fabs( pos[2] ) < fabs( pos[0] ) )
346 : : {
347 : : // pos[0] is the greatest
348 [ + + ]: 1154 : if( pos[0] > 0 )
349 : 572 : plane = 1;
350 : : else
351 : 1154 : plane = 3;
352 : : }
353 : : else
354 : : {
355 : : // pos[2] is biggest
356 [ + + ]: 647 : if( pos[2] < 0 )
357 : 323 : plane = 5;
358 : : else
359 : 324 : plane = 6;
360 : : }
361 : : }
362 : 3446 : return;
363 : : }
364 : :
365 : : // point on a sphere is projected on one of six planes, decided earlier
366 : 53678 : ErrorCode IntxUtils::gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 )
367 : : {
368 : 53678 : double alfa = 1.; // the new point will be on line alfa*pos
369 : :
370 [ + + + + : 53678 : switch( plane )
+ + - ]
371 : : {
372 : : case 1: {
373 : : // the plane with x = R; x>y, x>z
374 : : // c1->y, c2->z
375 : 8611 : alfa = R / pos[0];
376 : 8611 : c1 = alfa * pos[1];
377 : 8611 : c2 = alfa * pos[2];
378 : 8611 : break;
379 : : }
380 : : case 2: {
381 : : // y = R -> zx
382 : 8749 : alfa = R / pos[1];
383 : 8749 : c1 = alfa * pos[2];
384 : 8749 : c2 = alfa * pos[0];
385 : 8749 : break;
386 : : }
387 : : case 3: {
388 : : // x=-R, -> yz
389 : 8475 : alfa = -R / pos[0];
390 : 8475 : c1 = -alfa * pos[1]; // the sign is to preserve orientation
391 : 8475 : c2 = alfa * pos[2];
392 : 8475 : break;
393 : : }
394 : : case 4: {
395 : : // y = -R
396 : 8369 : alfa = -R / pos[1];
397 : 8369 : c1 = -alfa * pos[2]; // the sign is to preserve orientation
398 : 8369 : c2 = alfa * pos[0];
399 : 8369 : break;
400 : : }
401 : : case 5: {
402 : : // z = -R
403 : 9735 : alfa = -R / pos[2];
404 : 9735 : c1 = -alfa * pos[0]; // the sign is to preserve orientation
405 : 9735 : c2 = alfa * pos[1];
406 : 9735 : break;
407 : : }
408 : : case 6: {
409 : 9739 : alfa = R / pos[2];
410 : 9739 : c1 = alfa * pos[0];
411 : 9739 : c2 = alfa * pos[1];
412 : 9739 : break;
413 : : }
414 : : default:
415 : 0 : return MB_FAILURE; // error
416 : : }
417 : :
418 : 53678 : return MB_SUCCESS; // no error
419 : : }
420 : :
421 : : // given the position on plane (one out of 6), find out the position on sphere
422 : 23466 : ErrorCode IntxUtils::reverse_gnomonic_projection( const double& c1, const double& c2, double R, int plane,
423 : : CartVect& pos )
424 : : {
425 : :
426 : : // the new point will be on line beta*pos
427 : 23466 : double len = sqrt( c1 * c1 + c2 * c2 + R * R );
428 : 23466 : double beta = R / len; // it is less than 1, in general
429 : :
430 [ + + + + : 23466 : switch( plane )
+ + - ]
431 : : {
432 : : case 1: {
433 : : // the plane with x = R; x>y, x>z
434 : : // c1->y, c2->z
435 : 3285 : pos[0] = beta * R;
436 : 3285 : pos[1] = c1 * beta;
437 : 3285 : pos[2] = c2 * beta;
438 : 3285 : break;
439 : : }
440 : : case 2: {
441 : : // y = R -> zx
442 : 3390 : pos[1] = R * beta;
443 : 3390 : pos[2] = c1 * beta;
444 : 3390 : pos[0] = c2 * beta;
445 : 3390 : break;
446 : : }
447 : : case 3: {
448 : : // x=-R, -> yz
449 : 3334 : pos[0] = -R * beta;
450 : 3334 : pos[1] = -c1 * beta; // the sign is to preserve orientation
451 : 3334 : pos[2] = c2 * beta;
452 : 3334 : break;
453 : : }
454 : : case 4: {
455 : : // y = -R
456 : 3284 : pos[1] = -R * beta;
457 : 3284 : pos[2] = -c1 * beta; // the sign is to preserve orientation
458 : 3284 : pos[0] = c2 * beta;
459 : 3284 : break;
460 : : }
461 : : case 5: {
462 : : // z = -R
463 : 5085 : pos[2] = -R * beta;
464 : 5085 : pos[0] = -c1 * beta; // the sign is to preserve orientation
465 : 5085 : pos[1] = c2 * beta;
466 : 5085 : break;
467 : : }
468 : : case 6: {
469 : 5088 : pos[2] = R * beta;
470 : 5088 : pos[0] = c1 * beta;
471 : 5088 : pos[1] = c2 * beta;
472 : 5088 : break;
473 : : }
474 : : default:
475 : 0 : return MB_FAILURE; // error
476 : : }
477 : :
478 : 23466 : return MB_SUCCESS; // no error
479 : : }
480 : :
481 : 294 : void IntxUtils::gnomonic_unroll( double& c1, double& c2, double R, int plane )
482 : : {
483 : : double tmp;
484 [ + + + + : 294 : switch( plane )
+ + - ]
485 : : {
486 : : case 1:
487 : 49 : break; // nothing
488 : : case 2: // rotate + 90
489 : 49 : tmp = c1;
490 : 49 : c1 = -c2;
491 : 49 : c2 = tmp;
492 : 49 : c1 = c1 + 2 * R;
493 : 49 : break;
494 : : case 3:
495 : 49 : c1 = c1 + 4 * R;
496 : 49 : break;
497 : : case 4: // rotate with -90 x-> -y; y -> x
498 : :
499 : 49 : tmp = c1;
500 : 49 : c1 = c2;
501 : 49 : c2 = -tmp;
502 : 49 : c1 = c1 - 2 * R;
503 : 49 : break;
504 : : case 5: // South Pole
505 : : // rotate 180 then move to (-2, -2)
506 : 49 : c1 = -c1 - 2. * R;
507 : 49 : c2 = -c2 - 2. * R;
508 : 49 : break;
509 : : ;
510 : : case 6: // North Pole
511 : 49 : c1 = c1 - 2. * R;
512 : 49 : c2 = c2 + 2. * R;
513 : 49 : break;
514 : : }
515 : 294 : return;
516 : : }
517 : : // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
518 : 1 : ErrorCode IntxUtils::global_gnomonic_projection( Interface* mb, EntityHandle inSet, double R, bool centers_only,
519 : : EntityHandle& outSet )
520 : : {
521 [ + - ]: 1 : std::string parTagName( "PARALLEL_PARTITION" );
522 : : Tag part_tag;
523 [ + - ]: 2 : Range partSets;
524 [ + - ]: 1 : ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
525 [ - + ][ # # ]: 1 : if( MB_SUCCESS == rval && part_tag != 0 )
526 : : {
527 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
[ # # ][ # # ]
528 : : }
529 [ + - ][ - + ]: 1 : rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
[ # # ][ # # ]
530 : : // Get all entities of dimension 2
531 [ + - ]: 2 : Range inputRange; // get
532 [ + - ][ - + ]: 1 : rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
[ # # ][ # # ]
533 [ + - ][ - + ]: 1 : rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
[ # # ][ # # ]
534 : :
535 [ + - ]: 2 : std::map< EntityHandle, int > partsAssign;
536 [ + - ]: 2 : std::map< int, EntityHandle > newPartSets;
537 [ + - ][ - + ]: 1 : if( !partSets.empty() )
538 : : {
539 : : // get all cells, and assign parts
540 [ # # ][ # # ]: 0 : for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
[ # # ][ # # ]
[ # # ]
541 : : {
542 [ # # ]: 0 : EntityHandle pSet = *setIt;
543 [ # # ]: 0 : Range ents;
544 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
[ # # ][ # # ]
545 : : int val;
546 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
[ # # ][ # # ]
547 : : // create a new set with the same part id tag, in the outSet
548 : : EntityHandle newPartSet;
549 [ # # ][ # # ]: 0 : rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
[ # # ][ # # ]
550 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
[ # # ][ # # ]
551 [ # # ]: 0 : newPartSets[val] = newPartSet;
552 [ # # ][ # # ]: 0 : rval = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
553 [ # # ][ # # ]: 0 : for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
554 : : {
555 [ # # ][ # # ]: 0 : partsAssign[*it] = val;
556 : : }
557 : 0 : }
558 : : }
559 : :
560 [ - + ]: 1 : if( centers_only )
561 : : {
562 [ # # ][ # # ]: 0 : for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
[ # # ][ # # ]
[ # # ][ # # ]
563 : : {
564 [ # # ]: 0 : CartVect center;
565 [ # # ]: 0 : EntityHandle cell = *it;
566 [ # # ][ # # ]: 0 : rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
567 : : int plane;
568 [ # # ]: 0 : decide_gnomonic_plane( center, plane );
569 : : double c[3];
570 : 0 : c[2] = 0.;
571 [ # # ]: 0 : gnomonic_projection( center, R, plane, c[0], c[1] );
572 : :
573 [ # # ]: 0 : gnomonic_unroll( c[0], c[1], R, plane );
574 : :
575 : : EntityHandle vertex;
576 [ # # ][ # # ]: 0 : rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
[ # # ][ # # ]
577 [ # # ][ # # ]: 0 : rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
578 : : }
579 : : }
580 : : else
581 : : {
582 : : // distribute the cells to 6 planes, based on the center
583 [ + - ][ + + ]: 14 : Range subranges[6];
[ + - # # ]
584 [ + - ][ + - ]: 289 : for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
[ + - ][ + - ]
[ + + ]
585 : : {
586 [ + - ]: 288 : CartVect center;
587 [ + - ]: 288 : EntityHandle cell = *it;
588 [ + - ][ + - ]: 288 : rval = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
589 : : int plane;
590 [ + - ]: 288 : decide_gnomonic_plane( center, plane );
591 [ + - ]: 288 : subranges[plane - 1].insert( cell );
592 : : }
593 [ + + ][ + + ]: 13 : for( int i = 1; i <= 6; i++ )
594 : : {
595 [ + - ]: 6 : Range verts;
596 [ + - ][ - + ]: 6 : rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
[ # # ][ # # ]
597 [ + - ][ + - ]: 12 : std::map< EntityHandle, EntityHandle > corr;
598 [ + - ][ + - ]: 300 : for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
[ + - ][ + - ]
[ + + ]
599 : : {
600 [ + - ]: 294 : CartVect vect;
601 [ + - ]: 294 : EntityHandle v = *vt;
602 [ + - ][ + - ]: 294 : rval = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
[ - + ][ # # ]
[ # # ]
603 : : double c[3];
604 : 294 : c[2] = 0.;
605 [ + - ]: 294 : gnomonic_projection( vect, R, i, c[0], c[1] );
606 [ + - ]: 294 : gnomonic_unroll( c[0], c[1], R, i );
607 : : EntityHandle vertex;
608 [ + - ][ - + ]: 294 : rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
[ # # ][ # # ]
609 [ + - ]: 294 : corr[v] = vertex; // for new connectivity
610 : : }
611 : : EntityHandle new_conn[20]; // max edges in 2d ?
612 [ + - ][ + - ]: 294 : for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
[ + - ][ + - ]
[ + + ][ + - ]
613 : : {
614 [ + - ]: 288 : EntityHandle eh = *eit;
615 : 288 : const EntityHandle* conn = NULL;
616 : : int num_nodes;
617 [ + - ][ - + ]: 288 : rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
618 : : // build a new vertex array
619 [ + + ]: 1296 : for( int i = 0; i < num_nodes; i++ )
620 [ + - ]: 1008 : new_conn[i] = corr[conn[i]];
621 [ + - ]: 288 : EntityType type = mb->type_from_handle( eh );
622 : : EntityHandle newCell;
623 [ + - ][ - + ]: 288 : rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
[ # # ][ # # ]
624 [ + - ][ - + ]: 288 : rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
625 [ + - ]: 288 : std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
626 [ + - ][ - + ]: 288 : if( mit != partsAssign.end() )
627 : : {
628 [ # # ]: 0 : int val = mit->second;
629 [ # # ][ # # ]: 0 : rval = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
630 : : }
631 : : }
632 [ # # ]: 7 : }
633 : : }
634 : :
635 : 2 : return MB_SUCCESS;
636 : : }
637 : 0 : void IntxUtils::transform_coordinates( double* avg_position, int projection_type )
638 : : {
639 [ # # ]: 0 : if( projection_type == 1 )
640 : : {
641 : : double R =
642 : 0 : avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
643 : 0 : R = sqrt( R );
644 : 0 : double lat = asin( avg_position[2] / R );
645 : 0 : double lon = atan2( avg_position[1], avg_position[0] );
646 : 0 : avg_position[0] = lon;
647 : 0 : avg_position[1] = lat;
648 : 0 : avg_position[2] = R;
649 : : }
650 [ # # ]: 0 : else if( projection_type == 2 ) // gnomonic projection
651 : : {
652 [ # # ]: 0 : CartVect pos( avg_position );
653 : : int gplane;
654 [ # # ]: 0 : IntxUtils::decide_gnomonic_plane( pos, gplane );
655 : :
656 [ # # ]: 0 : IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
657 : 0 : avg_position[2] = 0;
658 [ # # ]: 0 : IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
659 : : }
660 : 0 : }
661 : : /*
662 : : *
663 : : use physical_constants, only : dd_pi
664 : : type(cartesian3D_t), intent(in) :: cart
665 : : type(spherical_polar_t) :: sphere
666 : :
667 : : sphere%r=distance(cart)
668 : : sphere%lat=ASIN(cart%z/sphere%r)
669 : : sphere%lon=0
670 : :
671 : : ! ==========================================================
672 : : ! enforce three facts:
673 : : !
674 : : ! 1) lon at poles is defined to be zero
675 : : !
676 : : ! 2) Grid points must be separated by about .01 Meter (on earth)
677 : : ! from pole to be considered "not the pole".
678 : : !
679 : : ! 3) range of lon is { 0<= lon < 2*pi }
680 : : !
681 : : ! ==========================================================
682 : :
683 : : if (distance(cart) >= DIST_THRESHOLD) then
684 : : sphere%lon=ATAN2(cart%y,cart%x)
685 : : if (sphere%lon<0) then
686 : : sphere%lon=sphere%lon+2*DD_PI
687 : : end if
688 : : end if
689 : :
690 : : end function cart_to_spherical
691 : : */
692 : 0 : IntxUtils::SphereCoords IntxUtils::cart_to_spherical( CartVect& cart3d )
693 : : {
694 : : SphereCoords res;
695 : 0 : res.R = cart3d.length();
696 [ # # ]: 0 : if( res.R < 0 )
697 : : {
698 : 0 : res.lon = res.lat = 0.;
699 : 0 : return res;
700 : : }
701 : 0 : res.lat = asin( cart3d[2] / res.R );
702 : 0 : res.lon = atan2( cart3d[1], cart3d[0] );
703 [ # # ]: 0 : if( res.lon < 0 ) res.lon += 2 * M_PI; // M_PI is defined in math.h? it seems to be true, although
704 : : // there are some defines it depends on :(
705 : : // #if defined __USE_BSD || defined __USE_XOPEN ???
706 : :
707 : 0 : return res;
708 : : }
709 : :
710 : : /*
711 : : * ! ===================================================================
712 : : ! spherical_to_cart:
713 : : ! converts spherical polar {lon,lat} to 3D cartesian {x,y,z}
714 : : ! on unit sphere
715 : : ! ===================================================================
716 : :
717 : : function spherical_to_cart(sphere) result (cart)
718 : :
719 : : type(spherical_polar_t), intent(in) :: sphere
720 : : type(cartesian3D_t) :: cart
721 : :
722 : : cart%x=sphere%r*COS(sphere%lat)*COS(sphere%lon)
723 : : cart%y=sphere%r*COS(sphere%lat)*SIN(sphere%lon)
724 : : cart%z=sphere%r*SIN(sphere%lat)
725 : :
726 : : end function spherical_to_cart
727 : : */
728 : 0 : CartVect IntxUtils::spherical_to_cart( IntxUtils::SphereCoords& sc )
729 : : {
730 : 0 : CartVect res;
731 : 0 : res[0] = sc.R * cos( sc.lat ) * cos( sc.lon ); // x coordinate
732 : 0 : res[1] = sc.R * cos( sc.lat ) * sin( sc.lon ); // y
733 : 0 : res[2] = sc.R * sin( sc.lat ); // z
734 : 0 : return res;
735 : : }
736 : :
737 : 3 : ErrorCode IntxUtils::ScaleToRadius( Interface* mb, EntityHandle set, double R )
738 : : {
739 [ + - ]: 3 : Range nodes;
740 [ + - ]: 3 : ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true ); // recursive
741 [ - + ]: 3 : if( rval != moab::MB_SUCCESS ) return rval;
742 : :
743 : : // one by one, get the node and project it on the sphere, with a radius given
744 : : // the center of the sphere is at 0,0,0
745 [ + - ][ + - ]: 495 : for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
[ + - ][ + - ]
[ + + ]
746 : : {
747 [ + - ]: 492 : EntityHandle nd = *nit;
748 [ + - ]: 492 : CartVect pos;
749 [ + - ][ + - ]: 492 : rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
750 [ - + ]: 492 : if( rval != moab::MB_SUCCESS ) return rval;
751 [ + - ]: 492 : double len = pos.length();
752 [ - + ]: 492 : if( len == 0. ) return MB_FAILURE;
753 [ + - ]: 492 : pos = R / len * pos;
754 [ + - ][ + - ]: 492 : rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
755 [ - + ]: 492 : if( rval != moab::MB_SUCCESS ) return rval;
756 : : }
757 : 3 : return MB_SUCCESS;
758 : : }
759 : :
760 : : // assume they are one the same sphere
761 : 0 : double IntxAreaUtils::spherical_angle( double* A, double* B, double* C, double Radius )
762 : : {
763 : : // the angle by definition is between the planes OAB and OBC
764 [ # # ]: 0 : CartVect a( A );
765 [ # # ]: 0 : CartVect b( B );
766 [ # # ]: 0 : CartVect c( C );
767 [ # # ]: 0 : double err1 = a.length_squared() - Radius * Radius;
768 [ # # ]: 0 : if( fabs( err1 ) > 0.0001 )
769 [ # # ][ # # ]: 0 : { std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n"; }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
770 [ # # ]: 0 : CartVect normalOAB = a * b;
771 [ # # ]: 0 : CartVect normalOCB = c * b;
772 [ # # ]: 0 : return angle( normalOAB, normalOCB );
773 : : }
774 : :
775 : : // could be bigger than M_PI;
776 : : // angle at B could be bigger than M_PI, if the orientation is such that ABC points toward the
777 : : // interior
778 : 864 : double IntxUtils::oriented_spherical_angle( double* A, double* B, double* C )
779 : : {
780 : : // assume the same radius, sphere at origin
781 [ + - ][ + - ]: 864 : CartVect a( A ), b( B ), c( C );
[ + - ]
782 [ + - ]: 864 : CartVect normalOAB = a * b;
783 [ + - ]: 864 : CartVect normalOCB = c * b;
784 [ + - ][ + - ]: 864 : CartVect orient = ( c - b ) * ( a - b );
[ + - ]
785 [ + - ]: 864 : double ang = angle( normalOAB, normalOCB ); // this is between 0 and M_PI
786 [ - + ]: 864 : if( ang != ang )
787 : : {
788 : : // signal of a nan
789 [ # # ][ # # ]: 0 : std::cout << a << " " << b << " " << c << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
790 [ # # ][ # # ]: 0 : std::cout << ang << "\n";
791 : : }
792 [ + - ][ - + ]: 864 : if( orient % b < 0 ) return ( 2 * M_PI - ang ); // the other angle, supplement
793 : :
794 : 864 : return ang;
795 : : }
796 : :
797 : 0 : double IntxAreaUtils::area_spherical_triangle( double* A, double* B, double* C, double Radius )
798 : : {
799 [ # # ]: 0 : switch( m_eAreaMethod )
800 : : {
801 : : case Girard:
802 : 0 : return area_spherical_triangle_girard( A, B, C, Radius );
803 : : #ifdef MOAB_HAVE_TEMPESTREMAP
804 : : case GaussQuadrature:
805 : : return area_spherical_triangle_GQ( A, B, C, Radius );
806 : : #endif
807 : : case lHuiller:
808 : : default:
809 : 0 : return area_spherical_triangle_lHuiller( A, B, C, Radius );
810 : : }
811 : : }
812 : :
813 : 1154 : double IntxAreaUtils::area_spherical_polygon( double* A, int N, double Radius, int* sign )
814 : : {
815 [ + + ]: 1154 : switch( m_eAreaMethod )
816 : : {
817 : : case Girard:
818 : 216 : return area_spherical_polygon_girard( A, N, Radius );
819 : : #ifdef MOAB_HAVE_TEMPESTREMAP
820 : : case GaussQuadrature:
821 : : return area_spherical_polygon_GQ( A, N, Radius );
822 : : #endif
823 : : case lHuiller:
824 : : default:
825 : 938 : return area_spherical_polygon_lHuiller( A, N, Radius, sign );
826 : : }
827 : : }
828 : :
829 : 0 : double IntxAreaUtils::area_spherical_triangle_girard( double* A, double* B, double* C, double Radius )
830 : : {
831 [ # # ][ # # ]: 0 : double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
832 [ # # ]: 0 : spherical_angle( C, A, B, Radius ) - M_PI;
833 : 0 : double area = Radius * Radius * correction;
834 : : // now, is it negative or positive? is it pointing toward the center or outward?
835 [ # # ][ # # ]: 0 : CartVect a( A ), b( B ), c( C );
[ # # ]
836 [ # # ][ # # ]: 0 : CartVect abc = ( b - a ) * ( c - a );
[ # # ]
837 [ # # ][ # # ]: 0 : if( abc % a > 0 ) // dot product positive, means ABC points out
838 : 0 : return area;
839 : : else
840 : 0 : return -area;
841 : : }
842 : :
843 : 216 : double IntxAreaUtils::area_spherical_polygon_girard( double* A, int N, double Radius )
844 : : {
845 : : // this should work for non-convex polygons too
846 : : // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
847 : : //
848 [ - + ]: 216 : if( N <= 2 ) return 0.;
849 : 216 : double sum_angles = 0.;
850 [ + + ]: 1080 : for( int i = 0; i < N; i++ )
851 : : {
852 : 864 : int i1 = ( i + 1 ) % N;
853 : 864 : int i2 = ( i + 2 ) % N;
854 : 864 : sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
855 : : }
856 : 216 : double correction = sum_angles - ( N - 2 ) * M_PI;
857 : 216 : return Radius * Radius * correction;
858 : : }
859 : :
860 : 3488 : double IntxAreaUtils::area_spherical_polygon_lHuiller( double* A, int N, double Radius, int* sign )
861 : : {
862 : : // This should work for non-convex polygons too
863 : : // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
864 : : // We also assume that the orientation is positive;
865 : : // If negative orientation, the area will be negative
866 [ - + ]: 3488 : if( N <= 2 ) return 0.;
867 : :
868 : 3488 : int lsign = 1; // assume positive orientain
869 : 3488 : double area = 0.;
870 [ + + ]: 10392 : for( int i = 1; i < N - 1; i++ )
871 : : {
872 : 6904 : int i1 = i + 1;
873 : 6904 : double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
874 [ + + ]: 6904 : if( areaTriangle < 0 )
875 : 5020 : lsign = -1; // signal that we have at least one triangle with negative orientation ;
876 : : // possible nonconvex polygon
877 : 6904 : area += areaTriangle;
878 : : }
879 [ - + ]: 3488 : if( sign ) *sign = lsign;
880 : :
881 : 3488 : return area;
882 : : }
883 : :
884 : : #ifdef MOAB_HAVE_TEMPESTREMAP
885 : : double IntxAreaUtils::area_spherical_polygon_GQ( double* A, int N, double Radius )
886 : : {
887 : : // this should work for non-convex polygons too
888 : : // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
889 : : // We also assume that the orientation is positive;
890 : : // If negative orientation, the area can be negative
891 : : if( N <= 2 ) return 0.;
892 : :
893 : : // assume positive orientain
894 : : double area = 0.;
895 : : for( int i = 1; i < N - 1; i++ )
896 : : {
897 : : int i1 = i + 1;
898 : : area += area_spherical_triangle_GQ( A, A + 3 * i, A + 3 * i1, Radius );
899 : : }
900 : : return area;
901 : : }
902 : :
903 : : /* compute the area by using Gauss-Quadratures; use TR interfaces directly */
904 : : double IntxAreaUtils::area_spherical_triangle_GQ( double* ptA, double* ptB, double* ptC, double )
905 : : {
906 : : Face face( 3 );
907 : : NodeVector nodes( 3 );
908 : : nodes[0] = Node( ptA[0], ptA[1], ptA[2] );
909 : : nodes[1] = Node( ptB[0], ptB[1], ptB[2] );
910 : : nodes[2] = Node( ptC[0], ptC[1], ptC[2] );
911 : : face.SetNode( 0, 0 );
912 : : face.SetNode( 1, 1 );
913 : : face.SetNode( 2, 2 );
914 : : return CalculateFaceArea( face, nodes );
915 : : }
916 : : #endif
917 : :
918 : : /*
919 : : * l'Huiller's formula for spherical triangle
920 : : * http://williams.best.vwh.net/avform.htm
921 : : * a, b, c are arc measures in radians, too
922 : : * A, B, C are angles on the sphere, for which we already have formula
923 : : * c
924 : : * A -------B
925 : : * \ |
926 : : * \ |
927 : : * \b |a
928 : : * \ |
929 : : * \ |
930 : : * \ |
931 : : * \C|
932 : : * \|
933 : : *
934 : : * (The angle at B is not necessarily a right angle)
935 : : *
936 : : * sin(a) sin(b) sin(c)
937 : : * ----- = ------ = ------
938 : : * sin(A) sin(B) sin(C)
939 : : *
940 : : * In terms of the sides (this is excess, as before, but numerically stable)
941 : : *
942 : : * E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
943 : : */
944 : 9454 : double IntxAreaUtils::area_spherical_triangle_lHuiller( double* ptA, double* ptB, double* ptC, double Radius )
945 : : {
946 : :
947 : : // now, a is angle BOC, O is origin
948 [ + - ][ + - ]: 9454 : CartVect vA( ptA ), vB( ptB ), vC( ptC );
[ + - ]
949 [ + - ]: 9454 : double a = angle( vB, vC );
950 [ + - ]: 9454 : double b = angle( vC, vA );
951 [ + - ]: 9454 : double c = angle( vA, vB );
952 : 9454 : int sign = 1;
953 [ + - ][ + - ]: 9454 : if( ( vA * vB ) % vC < 0 ) sign = -1;
[ + + ]
954 : 9454 : double s = ( a + b + c ) / 2;
955 : 9454 : double tmp = tan( s / 2 ) * tan( ( s - a ) / 2 ) * tan( ( s - b ) / 2 ) * tan( ( s - c ) / 2 );
956 [ - + ]: 9454 : if( tmp < 0. ) tmp = 0.;
957 : :
958 : 9454 : double E = 4 * atan( sqrt( tmp ) );
959 [ - + ][ # # ]: 9454 : if( E != E ) std::cout << " NaN at spherical triangle area \n";
960 : :
961 : 9454 : double area = sign * E * Radius * Radius;
962 : :
963 : : #ifdef CHECKNEGATIVEAREA
964 [ + + ]: 9454 : if( area < 0 )
965 : : {
966 [ + - ][ + - ]: 7570 : std::cout << "negative area: " << area << "\n";
[ + - ]
967 [ + - ][ + - ]: 7570 : std::cout << std::setprecision( 15 );
968 [ + - ][ + - ]: 7570 : std::cout << "vA: " << vA << "\n";
[ + - ]
969 [ + - ][ + - ]: 7570 : std::cout << "vB: " << vB << "\n";
[ + - ]
970 [ + - ][ + - ]: 7570 : std::cout << "vC: " << vC << "\n";
[ + - ]
971 [ + - ][ + - ]: 7570 : std::cout << "sign: " << sign << "\n";
[ + - ]
972 [ + - ][ + - ]: 7570 : std::cout << " a: " << a << "\n";
[ + - ]
973 [ + - ][ + - ]: 7570 : std::cout << " b: " << b << "\n";
[ + - ]
974 [ + - ][ + - ]: 7570 : std::cout << " c: " << c << "\n";
[ + - ]
975 : : }
976 : : #endif
977 : :
978 : 9454 : return area;
979 : : }
980 : : #undef CHECKNEGATIVEAREA
981 : :
982 : 4 : double IntxAreaUtils::area_on_sphere( Interface* mb, EntityHandle set, double R )
983 : : {
984 : : // Get all entities of dimension 2
985 [ + - ]: 4 : Range inputRange;
986 [ + - ][ - + ]: 4 : ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
[ # # ][ # # ]
987 : :
988 : : // Filter by elements that are owned by current process
989 [ + - ][ + - ]: 8 : std::vector< int > ownerinfo( inputRange.size(), -1 );
990 : : Tag intxOwnerTag;
991 [ + - ]: 4 : rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
992 [ - + ]: 4 : if( MB_SUCCESS == rval )
993 : : {
994 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
[ # # ][ # # ]
[ # # ]
995 : : }
996 : :
997 : : // compare total area with 4*M_PI * R^2
998 : 4 : int ie = 0;
999 : 4 : double total_area = 0.;
1000 [ + - ][ + - ]: 1158 : for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
[ + - ][ + - ]
[ + + ]
1001 : : {
1002 : :
1003 : : // All zero/positive owner data represents ghosted elems
1004 [ + - ][ - + ]: 1154 : if( ownerinfo[ie++] >= 0 ) continue;
1005 : :
1006 [ + - ]: 1154 : EntityHandle eh = *eit;
1007 [ + - ]: 1154 : const double elem_area = this->area_spherical_element( mb, eh, R );
1008 : :
1009 : : // check whether the area of the spherical element is positive.
1010 [ - + ]: 1154 : assert( elem_area > 0 );
1011 : :
1012 : : // sum up the contribution
1013 : 1154 : total_area += elem_area;
1014 : : }
1015 : :
1016 : : // return total mesh area
1017 : 8 : return total_area;
1018 : : }
1019 : :
1020 : 1154 : double IntxAreaUtils::area_spherical_element( Interface* mb, EntityHandle elem, double R )
1021 : : {
1022 : : // get the nodes, then the coordinates
1023 : : const EntityHandle* verts;
1024 : : int nsides;
1025 [ + - ][ - + ]: 1154 : ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
[ # # ][ # # ]
1026 : :
1027 : : // account for possible padded polygons
1028 [ - + ][ # # ]: 1154 : while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1029 : 0 : nsides--;
1030 : :
1031 : : // get coordinates
1032 [ + - ]: 1154 : std::vector< double > coords( 3 * nsides );
1033 [ + - ][ + - ]: 1154 : rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
[ - + ][ # # ]
[ # # ]
1034 : :
1035 : : // compute and return the area of the polygonal element
1036 [ + - ][ + - ]: 1154 : return area_spherical_polygon( &coords[0], nsides, R );
1037 : : }
1038 : :
1039 : 0 : double IntxUtils::distance_on_great_circle( CartVect& p1, CartVect& p2 )
1040 : : {
1041 [ # # ]: 0 : SphereCoords sph1 = cart_to_spherical( p1 );
1042 [ # # ]: 0 : SphereCoords sph2 = cart_to_spherical( p2 );
1043 : : // radius should be the same
1044 : 0 : return sph1.R *
1045 : 0 : acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
1046 : : }
1047 : :
1048 : : // break the nonconvex quads into triangles; remove the quad from the set? yes.
1049 : : // maybe radius is not needed;
1050 : : //
1051 : 0 : ErrorCode IntxUtils::enforce_convexity( Interface* mb, EntityHandle lset, int my_rank )
1052 : : {
1053 : : // look at each polygon; compute all angles; if one is reflex, break that angle with
1054 : : // the next triangle; put the 2 new polys in the set;
1055 : : // still look at the next poly
1056 : : // replace it with 2 triangles, and remove from set;
1057 : : // it should work for all polygons / tested first for case 1, with dt 0.5 (too much deformation)
1058 : : // get all entities of dimension 2
1059 : : // then get the connectivity, etc
1060 : :
1061 [ # # ]: 0 : Range inputRange;
1062 [ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );
1063 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1064 : :
1065 : 0 : Tag corrTag = 0;
1066 : 0 : EntityHandle dumH = 0;
1067 [ # # ]: 0 : rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
1068 [ # # ]: 0 : if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
1069 : :
1070 : : Tag gidTag;
1071 [ # # ]: 0 : rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );
1072 : :
1073 [ # # ]: 0 : if( rval != MB_SUCCESS ) return rval;
1074 : :
1075 [ # # ]: 0 : std::vector< double > coords;
1076 [ # # ]: 0 : coords.resize( 3 * MAXEDGES ); // at most 10 vertices per polygon
1077 : : // we should create a queue with new polygons that need processing for reflex angles
1078 : : // (obtuse)
1079 [ # # ][ # # ]: 0 : std::queue< EntityHandle > newPolys;
1080 : 0 : int brokenPolys = 0;
1081 [ # # ]: 0 : Range::iterator eit = inputRange.begin();
1082 [ # # ][ # # ]: 0 : while( eit != inputRange.end() || !newPolys.empty() )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # # # ]
1083 : : {
1084 : : EntityHandle eh;
1085 [ # # ][ # # ]: 0 : if( eit != inputRange.end() )
[ # # ]
1086 : : {
1087 [ # # ]: 0 : eh = *eit;
1088 [ # # ]: 0 : ++eit;
1089 : : }
1090 : : else
1091 : : {
1092 [ # # ]: 0 : eh = newPolys.front();
1093 [ # # ]: 0 : newPolys.pop();
1094 : : }
1095 : : // get the nodes, then the coordinates
1096 : : const EntityHandle* verts;
1097 : : int num_nodes;
1098 [ # # ]: 0 : rval = mb->get_connectivity( eh, verts, num_nodes );
1099 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1100 : 0 : int nsides = num_nodes;
1101 : : // account for possible padded polygons
1102 [ # # ][ # # ]: 0 : while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
1103 : 0 : nsides--;
1104 : 0 : EntityHandle corrHandle = 0;
1105 [ # # ]: 0 : if( corrTag )
1106 : : {
1107 [ # # ]: 0 : rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );
1108 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1109 : : }
1110 : 0 : int gid = 0;
1111 [ # # ]: 0 : rval = mb->tag_get_data( gidTag, &eh, 1, &gid );
1112 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1113 [ # # ]: 0 : coords.resize( 3 * nsides );
1114 [ # # ]: 0 : if( nsides < 4 ) continue; // if already triangles, don't bother
1115 : : // get coordinates
1116 [ # # ][ # # ]: 0 : rval = mb->get_coords( verts, nsides, &coords[0] );
1117 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1118 : : // compute each angle
1119 : 0 : bool alreadyBroken = false;
1120 : :
1121 [ # # ]: 0 : for( int i = 0; i < nsides; i++ )
1122 : : {
1123 [ # # ]: 0 : double* A = &coords[3 * i];
1124 [ # # ]: 0 : double* B = &coords[3 * ( ( i + 1 ) % nsides )];
1125 [ # # ]: 0 : double* C = &coords[3 * ( ( i + 2 ) % nsides )];
1126 [ # # ]: 0 : double angle = IntxUtils::oriented_spherical_angle( A, B, C );
1127 [ # # ]: 0 : if( angle - M_PI > 0. ) // even almost reflex is bad; break it!
1128 : : {
1129 [ # # ]: 0 : if( alreadyBroken )
1130 : : {
1131 [ # # ]: 0 : mb->list_entities( &eh, 1 );
1132 [ # # ]: 0 : mb->list_entities( verts, nsides );
1133 [ # # ]: 0 : double* D = &coords[3 * ( ( i + 3 ) % nsides )];
1134 [ # # ][ # # ]: 0 : std::cout << "ABC: " << angle << " \n";
[ # # ]
1135 [ # # ][ # # ]: 0 : std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
[ # # ][ # # ]
1136 [ # # ][ # # ]: 0 : std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
[ # # ][ # # ]
1137 [ # # ][ # # ]: 0 : std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
[ # # ][ # # ]
1138 [ # # ]: 0 : std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
1139 : :
1140 : 0 : return MB_FAILURE;
1141 : : }
1142 : : // the bad angle is at i+1;
1143 : : // create 1 triangle and one polygon; add the polygon to the input range, so
1144 : : // it will be processed too
1145 : : // also, add both to the set :) and remove the original polygon from the set
1146 : : // break the next triangle, even though not optimal
1147 : : // so create the triangle i+1, i+2, i+3; remove i+2 from original list
1148 : : // even though not optimal in general, it is good enough.
1149 : 0 : EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
1150 : 0 : verts[( i + 3 ) % nsides] };
1151 : : // create a polygon with num_nodes-1 vertices, and connectivity
1152 : : // verts[i+1], verts[i+3], (all except i+2)
1153 [ # # ]: 0 : std::vector< EntityHandle > conn( nsides - 1 );
1154 [ # # ]: 0 : for( int j = 1; j < nsides; j++ )
1155 : : {
1156 [ # # ]: 0 : conn[j - 1] = verts[( i + j + 2 ) % nsides];
1157 : : }
1158 : : EntityHandle newElement;
1159 [ # # ]: 0 : rval = mb->create_element( MBTRI, conn3, 3, newElement );
1160 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1161 : :
1162 [ # # ]: 0 : rval = mb->add_entities( lset, &newElement, 1 );
1163 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1164 [ # # ]: 0 : if( corrTag )
1165 : : {
1166 [ # # ]: 0 : rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );
1167 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1168 : : }
1169 [ # # ]: 0 : rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );
1170 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1171 [ # # ]: 0 : if( nsides == 4 )
1172 : : {
1173 : : // create another triangle
1174 [ # # ][ # # ]: 0 : rval = mb->create_element( MBTRI, &conn[0], 3, newElement );
1175 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1176 : : }
1177 : : else
1178 : : {
1179 : : // create another polygon, and add it to the inputRange
1180 [ # # ][ # # ]: 0 : rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );
1181 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1182 [ # # ]: 0 : newPolys.push( newElement ); // because it has less number of edges, the
1183 : : // reverse should work to find it.
1184 : : }
1185 [ # # ]: 0 : rval = mb->add_entities( lset, &newElement, 1 );
1186 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1187 [ # # ]: 0 : if( corrTag )
1188 : : {
1189 [ # # ]: 0 : rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );
1190 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1191 : : }
1192 [ # # ]: 0 : rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );
1193 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1194 [ # # ]: 0 : mb->remove_entities( lset, &eh, 1 );
1195 : 0 : brokenPolys++;
1196 : : /*std::cout<<"remove: " ;
1197 : : mb->list_entities(&eh, 1);
1198 : :
1199 : : std::stringstream fff;
1200 : : fff << "file0" << brokenQuads<< ".vtk";
1201 : : mb->write_file(fff.str().c_str(), 0, 0, &lset, 1);*/
1202 [ # # ]: 0 : alreadyBroken = true; // get out of the loop, element is broken
1203 : : }
1204 : : }
1205 : : }
1206 [ # # ][ # # ]: 0 : std::cout << "on local process " << my_rank << " " << brokenPolys
[ # # ][ # # ]
1207 [ # # ]: 0 : << " concave polygons were decomposed in convex ones \n";
1208 : 0 : return MB_SUCCESS;
1209 : : }
1210 : :
1211 : : // looking at quad connectivity, collapse to triangle if 2 nodes equal
1212 : : // then delete the old quad
1213 : 1 : ErrorCode IntxUtils::fix_degenerate_quads( Interface* mb, EntityHandle set )
1214 : : {
1215 [ + - ]: 1 : Range quads;
1216 [ + - ][ - + ]: 1 : ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
[ # # ][ # # ]
1217 : : Tag gid;
1218 [ + - ]: 1 : gid = mb->globalId_tag();
1219 [ + - ][ + - ]: 1201 : for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
[ + - ][ + - ]
[ + + ]
1220 : : {
1221 [ + - ]: 1200 : EntityHandle quad = *qit;
1222 : 1200 : const EntityHandle* conn4 = NULL;
1223 : 1200 : int num_nodes = 0;
1224 [ + - ][ - + ]: 1200 : rval = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
1225 [ + + ]: 6000 : for( int i = 0; i < num_nodes; i++ )
1226 : : {
1227 : 4800 : int next_node_index = ( i + 1 ) % num_nodes;
1228 [ + + ]: 4800 : if( conn4[i] == conn4[next_node_index] )
1229 : : {
1230 : : // form a triangle and delete the quad
1231 : : // first get the global id, to set it on triangle later
1232 : 80 : int global_id = 0;
1233 [ + - ][ - + ]: 80 : rval = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1234 : 80 : int i2 = ( i + 2 ) % num_nodes;
1235 : 80 : int i3 = ( i + 3 ) % num_nodes;
1236 : 80 : EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
1237 : : EntityHandle tri;
1238 [ + - ][ - + ]: 80 : rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
[ # # ][ # # ]
1239 [ + - ]: 80 : mb->add_entities( set, &tri, 1 );
1240 [ + - ]: 80 : mb->remove_entities( set, &quad, 1 );
1241 [ + - ]: 80 : mb->delete_entities( &quad, 1 );
1242 [ + - ][ - + ]: 80 : rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1243 : : }
1244 : : }
1245 : : }
1246 : 1 : return MB_SUCCESS;
1247 : : }
1248 : :
1249 : 4 : ErrorCode IntxAreaUtils::positive_orientation( Interface* mb, EntityHandle set, double R )
1250 : : {
1251 [ + - ]: 4 : Range cells2d;
1252 [ + - ]: 4 : ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );
1253 [ - + ]: 4 : if( MB_SUCCESS != rval ) return rval;
1254 [ + - ][ + - ]: 2634 : for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
[ + - ][ + - ]
[ + + ]
1255 : : {
1256 [ + - ]: 2630 : EntityHandle cell = *qit;
1257 : 2630 : const EntityHandle* conn = NULL;
1258 : 2630 : int num_nodes = 0;
1259 [ + - ]: 2630 : rval = mb->get_connectivity( cell, conn, num_nodes );
1260 [ - + ]: 2630 : if( MB_SUCCESS != rval ) return rval;
1261 [ - + ]: 2630 : if( num_nodes < 3 ) return MB_FAILURE;
1262 : :
1263 : : double coords[9];
1264 [ + - ]: 2630 : rval = mb->get_coords( conn, 3, coords );
1265 [ - + ]: 2630 : if( MB_SUCCESS != rval ) return rval;
1266 : :
1267 : : double area;
1268 [ + + ]: 2630 : if( R > 0 )
1269 [ + - ]: 2550 : area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
1270 : : else
1271 [ + - ]: 80 : area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
1272 [ + + ]: 2630 : if( area < 0 )
1273 : : {
1274 : : // compute all area, do not revert if total area is positive
1275 [ + - ]: 2550 : std::vector< double > coords2( 3 * num_nodes );
1276 : : // get coordinates
1277 [ + - ][ + - ]: 2550 : rval = mb->get_coords( conn, num_nodes, &coords2[0] );
1278 [ - + ]: 2550 : if( MB_SUCCESS != rval ) return MB_FAILURE;
1279 [ + - ][ + - ]: 2550 : double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
1280 [ + - ]: 2550 : if( totArea < 0 )
1281 : : {
1282 [ + - ]: 2550 : std::vector< EntityHandle > newconn( num_nodes );
1283 [ + + ]: 12670 : for( int i = 0; i < num_nodes; i++ )
1284 : : {
1285 [ + - ]: 10120 : newconn[num_nodes - 1 - i] = conn[i];
1286 : : }
1287 [ + - ][ + - ]: 2550 : rval = mb->set_connectivity( cell, &newconn[0], num_nodes );
1288 [ - + ][ + - ]: 2550 : if( MB_SUCCESS != rval ) return rval;
1289 : : }
1290 : : else
1291 : : {
1292 [ # # ][ # # ]: 2550 : std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
[ # # ][ # # ]
[ # # ][ + - ]
1293 : 2550 : }
1294 : : }
1295 : : }
1296 : 4 : return MB_SUCCESS;
1297 : : }
1298 : :
1299 : : // distance along a great circle on a sphere of radius 1
1300 : : // page 4
1301 : 0 : double IntxUtils::distance_on_sphere( double la1, double te1, double la2, double te2 )
1302 : : {
1303 : 0 : return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
1304 : : }
1305 : :
1306 : : /*
1307 : : * given 2 great circle arcs, AB and CD, compute the unique intersection point, if it exists
1308 : : * in between
1309 : : */
1310 : 2 : ErrorCode IntxUtils::intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E )
1311 : : {
1312 : : // first verify A, B, C, D are on the same sphere
1313 : 2 : double R2 = R * R;
1314 : 2 : const double Tolerance = 1.e-12 * R2;
1315 : :
1316 [ + - ][ + - ]: 2 : CartVect a( A ), b( B ), c( C ), d( D );
[ + - ][ + - ]
1317 : :
1318 [ + - ][ + - ]: 4 : if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
[ + - ][ - + ]
1319 [ + - ]: 2 : fabs( d.length_squared() - R2 ) >
1320 : 2 : 10 * Tolerance )
1321 : 0 : return MB_FAILURE;
1322 : :
1323 [ + - ]: 2 : CartVect n1 = a * b;
1324 [ + - ][ - + ]: 2 : if( n1.length_squared() < Tolerance ) return MB_FAILURE;
1325 : :
1326 [ + - ]: 2 : CartVect n2 = c * d;
1327 [ + - ][ - + ]: 2 : if( n2.length_squared() < Tolerance ) return MB_FAILURE;
1328 [ + - ]: 2 : CartVect n3 = n1 * n2;
1329 [ + - ]: 2 : n3.normalize();
1330 : :
1331 [ + - ]: 2 : n3 = R * n3;
1332 : : // the intersection is either n3 or -n3
1333 [ + - ][ + - ]: 2 : CartVect n4 = a * n3, n5 = n3 * b;
1334 [ + - ][ - + ]: 2 : if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
[ # # ][ # # ]
[ - + ]
1335 : : {
1336 : : // n3 is good for ab, see if it is good for cd
1337 [ # # ]: 0 : n4 = c * n3;
1338 [ # # ]: 0 : n5 = n3 * d;
1339 [ # # ][ # # ]: 0 : if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
[ # # ][ # # ]
[ # # ]
1340 : : {
1341 [ # # ]: 0 : E[0] = n3[0];
1342 [ # # ]: 0 : E[1] = n3[1];
1343 [ # # ]: 0 : E[2] = n3[2];
1344 : : }
1345 : : else
1346 : 0 : return MB_FAILURE;
1347 : : }
1348 : : else
1349 : : {
1350 : : // try -n3
1351 [ + - ]: 2 : n3 = -n3;
1352 [ + - ][ + - ]: 2 : n4 = a * n3, n5 = n3 * b;
1353 [ + - ][ + - ]: 2 : if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
[ + - ][ + + ]
[ + + ]
1354 : : {
1355 : : // n3 is good for ab, see if it is good for cd
1356 [ + - ]: 1 : n4 = c * n3;
1357 [ + - ]: 1 : n5 = n3 * d;
1358 [ + - ][ + - ]: 1 : if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
[ + - ][ + - ]
[ + - ]
1359 : : {
1360 [ + - ]: 1 : E[0] = n3[0];
1361 [ + - ]: 1 : E[1] = n3[1];
1362 [ + - ]: 1 : E[2] = n3[2];
1363 : : }
1364 : : else
1365 : 0 : return MB_FAILURE;
1366 : : }
1367 : : else
1368 : 1 : return MB_FAILURE;
1369 : : }
1370 : :
1371 : 2 : return MB_SUCCESS;
1372 : : }
1373 : :
1374 : : // verify that result is in between a and b on a great circle arc, and between c and d on a constant
1375 : : // latitude arc
1376 : 96448 : static bool verify( CartVect a, CartVect b, CartVect c, CartVect d, double x, double y, double z )
1377 : : {
1378 : : // to check, the point has to be between a and b on a great arc, and between c and d on a const
1379 : : // lat circle
1380 [ + - ]: 96448 : CartVect s( x, y, z );
1381 [ + - ]: 96448 : CartVect n1 = a * b;
1382 [ + - ]: 96448 : CartVect n2 = a * s;
1383 [ + - ]: 96448 : CartVect n3 = s * b;
1384 [ + - ][ + + ]: 96448 : if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
[ + - ][ + + ]
[ + + ]
1385 : :
1386 : : // do the same for c, d, s, in plane z=0
1387 [ + - ][ + - ]: 15716 : c[2] = d[2] = s[2] = 0.; // bring everything in the same plane, z=0;
[ + - ]
1388 : :
1389 [ + - ]: 15716 : n1 = c * d;
1390 [ + - ]: 15716 : n2 = c * s;
1391 [ + - ]: 15716 : n3 = s * d;
1392 [ + - ][ + + ]: 15716 : if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
[ + - ][ + + ]
[ + + ]
1393 : :
1394 : 96448 : return true;
1395 : : }
1396 : :
1397 : 59842 : ErrorCode IntxUtils::intersect_great_circle_arc_with_clat_arc( double* A, double* B, double* C, double* D, double R,
1398 : : double* E, int& np )
1399 : : {
1400 : 59842 : const double distTol = R * 1.e-6;
1401 : 59842 : const double Tolerance = R * R * 1.e-12; // radius should be 1, usually
1402 : 59842 : np = 0; // number of points in intersection
1403 [ + - ][ + - ]: 59842 : CartVect a( A ), b( B ), c( C ), d( D );
[ + - ][ + - ]
1404 : : // check input first
1405 : 59842 : double R2 = R * R;
1406 [ + - ][ + - ]: 119684 : if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
[ + - ][ - + ]
1407 [ + - ]: 59842 : fabs( d.length_squared() - R2 ) >
1408 : 59842 : 10 * Tolerance )
1409 : 0 : return MB_FAILURE;
1410 : :
1411 [ + - ][ + - ]: 59842 : if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
[ - + ]
1412 [ + - ][ + - ]: 59842 : if( ( c - d ).length_squared() < Tolerance ) // edges are too short
[ - + ]
1413 : 0 : return MB_FAILURE;
1414 : :
1415 : : // CD is the const latitude arc
1416 [ - + ]: 59842 : if( fabs( C[2] - D[2] ) > distTol ) // cd is not on the same z (constant latitude)
1417 : 0 : return MB_FAILURE;
1418 : :
1419 [ + - ][ - + ]: 59842 : if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE; // too close to the poles
1420 : :
1421 : : // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
1422 : : // great circle arc AB normal to the AB circle:
1423 [ + - ]: 59842 : CartVect n1 = a * b; // the normal to the great circle arc (circle)
1424 : : // solve the system of equations:
1425 : : /*
1426 : : * n1%(x, y, z) = 0 // on the great circle
1427 : : * z = C[2];
1428 : : * x^2+y^2+z^2 = R^2
1429 : : */
1430 : 59842 : double z = C[2];
1431 [ + - ][ + - ]: 59842 : if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
[ - + ]
1432 : : {
1433 : : // it is the Equator; check if the const lat edge is Equator too
1434 [ # # ]: 0 : if( fabs( C[2] ) > distTol )
1435 : : {
1436 : 0 : return MB_FAILURE; // no intx, too far from Eq
1437 : : }
1438 : : else
1439 : : {
1440 : : // all points are on the equator
1441 : : //
1442 [ # # ]: 0 : CartVect cd = c * d;
1443 : : // by convention, c<d, positive is from c to d
1444 : : // is a or b between c , d?
1445 [ # # ]: 0 : CartVect ca = c * a;
1446 [ # # ]: 0 : CartVect ad = a * d;
1447 [ # # ]: 0 : CartVect cb = c * b;
1448 [ # # ]: 0 : CartVect bd = b * d;
1449 [ # # ]: 0 : bool agtc = ( ca % cd >= -Tolerance ); // a>c?
1450 [ # # ]: 0 : bool dgta = ( ad % cd >= -Tolerance ); // d>a?
1451 [ # # ]: 0 : bool bgtc = ( cb % cd >= -Tolerance ); // b>c?
1452 [ # # ]: 0 : bool dgtb = ( bd % cd >= -Tolerance ); // d>b?
1453 [ # # ]: 0 : if( agtc )
1454 : : {
1455 [ # # ]: 0 : if( dgta )
1456 : : {
1457 : : // a is for sure a point
1458 [ # # ]: 0 : E[0] = a[0];
1459 [ # # ]: 0 : E[1] = a[1];
1460 [ # # ]: 0 : E[2] = a[2];
1461 : 0 : np++;
1462 [ # # ]: 0 : if( bgtc )
1463 : : {
1464 [ # # ]: 0 : if( dgtb )
1465 : : {
1466 : : // b is also in between c and d
1467 [ # # ]: 0 : E[3] = b[0];
1468 [ # # ]: 0 : E[4] = b[1];
1469 [ # # ]: 0 : E[5] = b[2];
1470 : 0 : np++;
1471 : : }
1472 : : else
1473 : : {
1474 : : // then order is c a d b, intx is ad
1475 [ # # ]: 0 : E[3] = d[0];
1476 [ # # ]: 0 : E[4] = d[1];
1477 [ # # ]: 0 : E[5] = d[2];
1478 : 0 : np++;
1479 : : }
1480 : : }
1481 : : else
1482 : : {
1483 : : // b is less than c, so b c a d, intx is ac
1484 [ # # ]: 0 : E[3] = c[0];
1485 [ # # ]: 0 : E[4] = c[1];
1486 [ # # ]: 0 : E[5] = c[2];
1487 : 0 : np++; // what if E[0] is E[3]?
1488 : : }
1489 : : }
1490 : : else // c < d < a
1491 : : {
1492 [ # # ]: 0 : if( dgtb ) // d is for sure in
1493 : : {
1494 [ # # ]: 0 : E[0] = d[0];
1495 [ # # ]: 0 : E[1] = d[1];
1496 [ # # ]: 0 : E[2] = d[2];
1497 : 0 : np++;
1498 [ # # ]: 0 : if( bgtc ) // c<b<d<a
1499 : : {
1500 : : // another point is b
1501 [ # # ]: 0 : E[3] = b[0];
1502 [ # # ]: 0 : E[4] = b[1];
1503 [ # # ]: 0 : E[5] = b[2];
1504 : 0 : np++;
1505 : : }
1506 : : else // b<c<d<a
1507 : : {
1508 : : // another point is c
1509 [ # # ]: 0 : E[3] = c[0];
1510 [ # # ]: 0 : E[4] = c[1];
1511 [ # # ]: 0 : E[5] = c[2];
1512 : 0 : np++;
1513 : : }
1514 : : }
1515 : : else
1516 : : {
1517 : : // nothing, order is c, d < a, b
1518 : : }
1519 : : }
1520 : : }
1521 : : else // a < c < d
1522 : : {
1523 [ # # ]: 0 : if( bgtc )
1524 : : {
1525 : : // c is for sure in
1526 [ # # ]: 0 : E[0] = c[0];
1527 [ # # ]: 0 : E[1] = c[1];
1528 [ # # ]: 0 : E[2] = c[2];
1529 : 0 : np++;
1530 [ # # ]: 0 : if( dgtb )
1531 : : {
1532 : : // a < c < b < d; second point is b
1533 [ # # ]: 0 : E[3] = b[0];
1534 [ # # ]: 0 : E[4] = b[1];
1535 [ # # ]: 0 : E[5] = b[2];
1536 : 0 : np++;
1537 : : }
1538 : : else
1539 : : {
1540 : : // a < c < d < b; second point is d
1541 [ # # ]: 0 : E[3] = d[0];
1542 [ # # ]: 0 : E[4] = d[1];
1543 [ # # ]: 0 : E[5] = d[2];
1544 : 0 : np++;
1545 : : }
1546 : : }
1547 : : else // a, b < c < d
1548 : : {
1549 : : // nothing
1550 : : }
1551 : : }
1552 : : }
1553 : : // for the 2 points selected, see if it is only one?
1554 : : // no problem, maybe it will be collapsed later anyway
1555 [ # # ]: 0 : if( np > 0 ) return MB_SUCCESS;
1556 : 0 : return MB_FAILURE; // no intersection
1557 : : }
1558 : : {
1559 [ + - ][ + - ]: 59842 : if( fabs( n1[0] ) <= fabs( n1[1] ) )
[ + + ]
1560 : : {
1561 : : // resolve eq in x: n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
1562 : : // (u+v*x)^2+x^2=R2-z^2
1563 : : // (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
1564 : : // delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
1565 : : // x1,2 =
1566 [ + - ][ + - ]: 30516 : double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
[ + - ][ + - ]
1567 : 30516 : double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1568 : 30516 : double delta = b1 * b1 - 4 * a1 * c1;
1569 [ + + ]: 30516 : if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1570 [ + - ]: 24693 : if( delta > Tolerance ) // 2 solutions possible
1571 : : {
1572 : 24693 : double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1573 : 24693 : double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1574 : 24693 : double y1 = u + v * x1;
1575 : 24693 : double y2 = u + v * x2;
1576 [ + - ][ + + ]: 24693 : if( verify( a, b, c, d, x1, y1, z ) )
1577 : : {
1578 : 2004 : E[0] = x1;
1579 : 2004 : E[1] = y1;
1580 : 2004 : E[2] = z;
1581 : 2004 : np++;
1582 : : }
1583 [ + - ][ + + ]: 24693 : if( verify( a, b, c, d, x2, y2, z ) )
1584 : : {
1585 : 2013 : E[3 * np + 0] = x2;
1586 : 2013 : E[3 * np + 1] = y2;
1587 : 2013 : E[3 * np + 2] = z;
1588 : 24693 : np++;
1589 : : }
1590 : : }
1591 : : else
1592 : : {
1593 : : // one solution
1594 : 0 : double x1 = -b1 / 2 / a1;
1595 : 0 : double y1 = u + v * x1;
1596 [ # # ][ # # ]: 0 : if( verify( a, b, c, d, x1, y1, z ) )
1597 : : {
1598 : 0 : E[0] = x1;
1599 : 0 : E[1] = y1;
1600 : 0 : E[2] = z;
1601 : 24693 : np++;
1602 : : }
1603 : : }
1604 : : }
1605 : : else
1606 : : {
1607 : : // resolve eq in y, reverse
1608 : : // n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
1609 : : // (u+v*y)^2+y^2 -R2+z^2 =0
1610 : : // (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
1611 : : //
1612 : : // x1,2 =
1613 [ + - ][ + - ]: 29326 : double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
[ + - ][ + - ]
1614 : 29326 : double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
1615 : 29326 : double delta = b1 * b1 - 4 * a1 * c1;
1616 [ + + ]: 29326 : if( delta < -Tolerance ) return MB_FAILURE; // no intersection
1617 [ + - ]: 23531 : if( delta > Tolerance ) // 2 solutions possible
1618 : : {
1619 : 23531 : double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
1620 : 23531 : double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
1621 : 23531 : double x1 = u + v * y1;
1622 : 23531 : double x2 = u + v * y2;
1623 [ + - ][ + + ]: 23531 : if( verify( a, b, c, d, x1, y1, z ) )
1624 : : {
1625 : 1883 : E[0] = x1;
1626 : 1883 : E[1] = y1;
1627 : 1883 : E[2] = z;
1628 : 1883 : np++;
1629 : : }
1630 [ + - ][ + + ]: 23531 : if( verify( a, b, c, d, x2, y2, z ) )
1631 : : {
1632 : 1835 : E[3 * np + 0] = x2;
1633 : 1835 : E[3 * np + 1] = y2;
1634 : 1835 : E[3 * np + 2] = z;
1635 : 23531 : np++;
1636 : : }
1637 : : }
1638 : : else
1639 : : {
1640 : : // one solution
1641 : 0 : double y1 = -b1 / 2 / a1;
1642 : 0 : double x1 = u + v * y1;
1643 [ # # ][ # # ]: 0 : if( verify( a, b, c, d, x1, y1, z ) )
1644 : : {
1645 : 0 : E[0] = x1;
1646 : 0 : E[1] = y1;
1647 : 0 : E[2] = z;
1648 : 0 : np++;
1649 : : }
1650 : : }
1651 : : }
1652 : : }
1653 : :
1654 [ + + ]: 48224 : if( np <= 0 ) return MB_FAILURE;
1655 : 59842 : return MB_SUCCESS;
1656 : : }
1657 : :
1658 : : #if 0
1659 : : ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1)
1660 : : {
1661 : : Range cells;
1662 : : ErrorCode rval = mb->get_entities_by_dimension(sf1, 2, cells);
1663 : : if (MB_SUCCESS!= rval)
1664 : : return rval;
1665 : : Range edges;
1666 : : rval = mb->get_adjacencies(cells, 1, true, edges, Interface::UNION);
1667 : : if (MB_SUCCESS!= rval)
1668 : : return rval;
1669 : :
1670 : : Tag edgeTypeTag;
1671 : : int default_int=0;
1672 : : rval = mb->tag_get_handle("edge_type", 1, MB_TYPE_INTEGER, edgeTypeTag,
1673 : : MB_TAG_DENSE | MB_TAG_CREAT, &default_int);
1674 : : if (MB_SUCCESS!= rval)
1675 : : return rval;
1676 : : // add edges to the set? not yet, maybe later
1677 : : // if edge horizontal, set value to 1
1678 : : int type_constant_lat=1;
1679 : : for (Range::iterator eit=edges.begin(); eit!=edges.end(); ++eit)
1680 : : {
1681 : : EntityHandle edge = *eit;
1682 : : const EntityHandle *conn=0;
1683 : : int num_n=0;
1684 : : rval = mb->get_connectivity(edge, conn, num_n );
1685 : : if (MB_SUCCESS!= rval)
1686 : : return rval;
1687 : : double coords[6];
1688 : : rval = mb->get_coords(conn, 2, coords);
1689 : : if (MB_SUCCESS!= rval)
1690 : : return rval;
1691 : : if (fabs( coords[2]-coords[5] )< 1.e-6 )
1692 : : {
1693 : : rval = mb->tag_set_data(edgeTypeTag, &edge, 1, &type_constant_lat);
1694 : : if (MB_SUCCESS!= rval)
1695 : : return rval;
1696 : : }
1697 : : }
1698 : :
1699 : : return MB_SUCCESS;
1700 : : }
1701 : : #endif
1702 : :
1703 : : // decide in a different metric if the corners of CS quad are
1704 : : // in the interior of an RLL quad
1705 : 7589 : int IntxUtils::borderPointsOfCSinRLL( CartVect* redc, double* red2dc, int nsRed, CartVect* bluec, int nsBlue,
1706 : : int* blueEdgeType, double* P, int* side, double epsil )
1707 : : {
1708 : 7589 : int extraPoints = 0;
1709 : : // first decide the blue z coordinates
1710 [ + - ][ + - ]: 7589 : CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
[ + - ][ + - ]
1711 [ + - ]: 22659 : for( int i = 0; i < nsBlue; i++ )
1712 : : {
1713 [ + + ]: 22659 : if( blueEdgeType[i] == 0 )
1714 : : {
1715 : 15070 : int iP1 = ( i + 1 ) % nsBlue;
1716 [ + - ][ + - ]: 15070 : if( bluec[i][2] > bluec[iP1][2] )
[ + + ]
1717 : : {
1718 : 7589 : A = bluec[i];
1719 : 7589 : B = bluec[iP1];
1720 : 7589 : C = bluec[( i + 2 ) % nsBlue];
1721 : 7589 : D = bluec[( i + 3 ) % nsBlue]; // it could be back to A, if triangle on top
1722 : 15070 : break;
1723 : : }
1724 : : }
1725 : : }
1726 [ + + ][ + - ]: 7589 : if( nsBlue == 3 && B[2] < 0 )
[ + + ][ + + ]
1727 : : {
1728 : : // select D to be C
1729 : 108 : D = C;
1730 : 108 : C = B; // B is the south pole then
1731 : : }
1732 : : // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
1733 : : // AB is const longitude, BC and DA constant latitude
1734 : : // check now each of the red points if they are inside this rectangle
1735 [ + + ]: 37945 : for( int i = 0; i < nsRed; i++ )
1736 : : {
1737 : 30356 : CartVect& X = redc[i];
1738 [ + - ][ + - ]: 30356 : if( X[2] > A[2] || X[2] < B[2] ) continue; // it is above or below the rectangle
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ]
1739 : : // now decide if it is between the planes OAB and OCD
1740 [ + - ][ + - ]: 15239 : if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
[ + + ][ + - ]
[ + - ][ + + ]
[ + + ][ + - ]
[ + + # #
# # ]
1741 : : {
1742 : 9043 : side[i] = 1; //
1743 : : // it means point X is in the rectangle that we want , on the sphere
1744 : : // pass the coords 2d
1745 : 9043 : P[extraPoints * 2] = red2dc[2 * i];
1746 : 9043 : P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
1747 : 9043 : extraPoints++;
1748 : : }
1749 : : }
1750 : 7589 : return extraPoints;
1751 : : }
1752 : :
1753 : 0 : ErrorCode IntxUtils::deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set )
1754 : : {
1755 : : ReadUtilIface* read_iface;
1756 [ # # ][ # # ]: 0 : ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
[ # # ][ # # ]
1757 : : // create the handle tag for the corresponding element / vertex
1758 : :
1759 : 0 : EntityHandle dum = 0;
1760 : 0 : Tag corrTag = 0; // it will be created here
1761 [ # # ][ # # ]: 0 : rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
[ # # ][ # # ]
1762 : :
1763 : : // give the same global id to new verts and cells created in the lagr(departure) mesh
1764 [ # # ]: 0 : Tag gid = mb->globalId_tag();
1765 : :
1766 [ # # ]: 0 : Range quads;
1767 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
[ # # ][ # # ]
1768 : :
1769 [ # # ]: 0 : Range connecVerts;
1770 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1771 : :
1772 [ # # ]: 0 : std::map< EntityHandle, EntityHandle > newNodes;
1773 : :
1774 [ # # ]: 0 : std::vector< double* > coords;
1775 : : EntityHandle start_vert, start_elem, *connect;
1776 [ # # ]: 0 : int num_verts = connecVerts.size();
1777 [ # # ]: 0 : rval = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
1778 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1779 : : // fill it up
1780 : 0 : int i = 0;
1781 [ # # ][ # # ]: 0 : for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
[ # # ][ # # ]
[ # # ]
1782 : : {
1783 [ # # ]: 0 : EntityHandle oldV = *vit;
1784 [ # # ]: 0 : CartVect posi;
1785 [ # # ][ # # ]: 0 : rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
1786 : :
1787 : : int global_id;
1788 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1789 : 0 : EntityHandle new_vert = start_vert + i;
1790 : : // Cppcheck warning (false positive): variable coords is assigned a value that is never used
1791 [ # # ][ # # ]: 0 : coords[0][i] = posi[0];
1792 [ # # ][ # # ]: 0 : coords[1][i] = posi[1];
1793 [ # # ][ # # ]: 0 : coords[2][i] = posi[2];
1794 : :
1795 [ # # ]: 0 : newNodes[oldV] = new_vert;
1796 : : // set also the correspondent tag :)
1797 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
[ # # ][ # # ]
1798 : :
1799 : : // also the other side
1800 : : // need to check if we really need this; the new vertex will never need the old vertex
1801 : : // we have the global id which is the same
1802 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
[ # # ][ # # ]
1803 : : // set the global id on the corresponding vertex the same as the initial vertex
1804 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1805 : : }
1806 : : // now create new quads in order (in a sequence)
1807 : :
1808 [ # # ][ # # ]: 0 : rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
1809 [ # # ]: 0 : if( MB_SUCCESS != rval ) return rval;
1810 : 0 : int ie = 0;
1811 [ # # ][ # # ]: 0 : for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
[ # # ][ # # ]
[ # # ]
1812 : : {
1813 [ # # ]: 0 : EntityHandle q = *it;
1814 : : int nnodes;
1815 : : const EntityHandle* conn;
1816 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
[ # # ][ # # ]
1817 : : int global_id;
1818 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1819 : :
1820 [ # # ]: 0 : for( int ii = 0; ii < nnodes; ii++ )
1821 : : {
1822 : 0 : EntityHandle v1 = conn[ii];
1823 [ # # ]: 0 : connect[4 * ie + ii] = newNodes[v1];
1824 : : }
1825 : 0 : EntityHandle newElement = start_elem + ie;
1826 : :
1827 : : // set the corresponding tag; not sure we need this one, from old to new
1828 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
[ # # ][ # # ]
1829 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
[ # # ][ # # ]
1830 : :
1831 : : // set the global id
1832 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
[ # # ][ # # ]
1833 : :
1834 [ # # ][ # # ]: 0 : rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
[ # # ][ # # ]
1835 : : }
1836 : :
1837 [ # # ][ # # ]: 0 : rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
[ # # ][ # # ]
[ # # ]
1838 : :
1839 : 0 : return MB_SUCCESS;
1840 : : }
1841 : :
1842 : 0 : ErrorCode IntxUtils::remove_duplicate_vertices( Interface* mb, EntityHandle file_set, double merge_tol,
1843 : : std::vector< Tag >& tagList )
1844 : : {
1845 [ # # ]: 0 : Range verts;
1846 [ # # ][ # # ]: 0 : ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1847 [ # # ][ # # ]: 0 : rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1848 : :
1849 [ # # ]: 0 : MergeMesh mm( mb );
1850 : :
1851 : : // remove the vertices from the set, before merging
1852 : :
1853 [ # # ][ # # ]: 0 : rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
[ # # ][ # # ]
1854 : :
1855 : : // now correct vertices that are repeated in polygons
1856 [ # # ]: 0 : Range cells;
1857 [ # # ][ # # ]: 0 : rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
[ # # ][ # # ]
1858 : :
1859 [ # # ]: 0 : verts.clear();
1860 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
[ # # ][ # # ]
1861 : :
1862 [ # # ]: 0 : Range modifiedCells; // will be deleted at the end; keep the gid
1863 [ # # ]: 0 : Range newCells;
1864 : :
1865 [ # # ][ # # ]: 0 : for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
[ # # ][ # # ]
[ # # ]
1866 : : {
1867 [ # # ]: 0 : EntityHandle cell = *cit;
1868 : 0 : const EntityHandle* connec = NULL;
1869 : 0 : int num_verts = 0;
1870 [ # # ][ # # ]: 0 : rval = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1871 : :
1872 [ # # ]: 0 : std::vector< EntityHandle > newConnec;
1873 [ # # ]: 0 : newConnec.push_back( connec[0] ); // at least one vertex
1874 : 0 : int index = 0;
1875 : 0 : int new_size = 1;
1876 [ # # ]: 0 : while( index < num_verts - 2 )
1877 : : {
1878 : 0 : int next_index = ( index + 1 );
1879 [ # # ][ # # ]: 0 : if( connec[next_index] != newConnec[new_size - 1] )
1880 : : {
1881 [ # # ]: 0 : newConnec.push_back( connec[next_index] );
1882 : 0 : new_size++;
1883 : : }
1884 : 0 : index++;
1885 : : }
1886 : : // add the last one only if different from previous and first node
1887 [ # # ][ # # ]: 0 : if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
1888 : : {
1889 [ # # ]: 0 : newConnec.push_back( connec[num_verts - 1] );
1890 : 0 : new_size++;
1891 : : }
1892 [ # # ]: 0 : if( new_size < num_verts )
1893 : : {
1894 : : // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
1895 [ # # ]: 0 : modifiedCells.insert( cell );
1896 : : // create a new cell with type triangle, quad or polygon
1897 : 0 : EntityType type = MBTRI;
1898 [ # # ]: 0 : if( new_size == 3 )
1899 : 0 : type = MBTRI;
1900 [ # # ]: 0 : else if( new_size == 4 )
1901 : 0 : type = MBQUAD;
1902 [ # # ]: 0 : else if( new_size > 4 )
1903 : 0 : type = MBPOLYGON;
1904 : :
1905 : : // create new cell
1906 : : EntityHandle newCell;
1907 [ # # ][ # # ]: 0 : rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1908 : : // set the old id to the new element
1909 [ # # ]: 0 : newCells.insert( newCell );
1910 : : double value; // use the same value to reset the tags, even if the tags are int (like Global ID)
1911 [ # # ][ # # ]: 0 : for( size_t i = 0; i < tagList.size(); i++ )
1912 : : {
1913 [ # # ][ # # ]: 0 : rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1914 [ # # ][ # # ]: 0 : rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
1915 : : }
1916 : : }
1917 : 0 : }
1918 : :
1919 [ # # ][ # # ]: 0 : rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1920 [ # # ][ # # ]: 0 : rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1921 [ # # ][ # # ]: 0 : rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1922 [ # # ][ # # ]: 0 : rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1923 : :
1924 : 0 : return MB_SUCCESS;
1925 : : }
1926 : :
1927 [ + - ][ + - ]: 228 : } // namespace moab
|
