MOAB: Mesh Oriented datABase  (version 5.2.1)
IntxUtils.cpp
Go to the documentation of this file.
00001 /*
00002  * IntxUtils.cpp
00003  *
00004  *  Created on: Oct 3, 2012
00005  */
00006 
00007 #ifdef WIN32               /* windows */
00008 #define _USE_MATH_DEFINES  // For M_PI
00009 #endif
00010 #include <cmath>
00011 #include <cassert>
00012 #include <iostream>
00013 // this is for sstream
00014 #include <sstream>
00015 
00016 #include "moab/IntxMesh/IntxUtils.hpp"
00017 // this is from mbcoupler; maybe it should be moved somewhere in moab src
00018 // right now, add a dependency to mbcoupler
00019 // #include "ElemUtil.hpp"
00020 #include "moab/MergeMesh.hpp"
00021 #include "moab/ReadUtilIface.hpp"
00022 #include "MBTagConventions.hpp"
00023 #define CHECKNEGATIVEAREA
00024 
00025 #ifdef CHECKNEGATIVEAREA
00026 #include <iomanip>
00027 #endif
00028 #include <queue>
00029 #include <map>
00030 
00031 #ifdef MOAB_HAVE_TEMPESTREMAP
00032 #include "GridElements.h"
00033 #endif
00034 
00035 namespace moab
00036 {
00037 
00038 #define CORRTAGNAME "__correspondent"
00039 #define MAXEDGES    10
00040 
00041 int IntxUtils::borderPointsOfXinY2( double* X, int nX, double* Y, int nY, double* P, int* side, double epsilon_area )
00042 {
00043     // 2 triangles, 3 corners, is the corner of X in Y?
00044     // Y must have a positive area
00045     /*
00046      */
00047     int extraPoint = 0;
00048     for( int i = 0; i < nX; i++ )
00049     {
00050         // compute double the area of all nY triangles formed by a side of Y and a corner of X; if
00051         // one is negative, stop (negative means it is outside; X and Y are all oriented such that
00052         // they are positive oriented;
00053         //  if one area is negative, it means it is outside the convex region, for sure)
00054         double* A = X + 2 * i;
00055 
00056         int inside = 1;
00057         for( int j = 0; j < nY; j++ )
00058         {
00059             double* B = Y + 2 * j;
00060             int j1    = ( j + 1 ) % nY;
00061             double* C = Y + 2 * j1;  // no copy of data
00062 
00063             double area2 = ( B[0] - A[0] ) * ( C[1] - A[1] ) - ( C[0] - A[0] ) * ( B[1] - A[1] );
00064             if( area2 < -epsilon_area )
00065             {
00066                 inside = 0;
00067                 break;
00068             }
00069         }
00070         if( inside )
00071         {
00072             side[i] = 1;  // so vertex i of X is inside the convex region formed by Y
00073             // so side has nX dimension (first array)
00074             P[extraPoint * 2]     = A[0];
00075             P[extraPoint * 2 + 1] = A[1];
00076             extraPoint++;
00077         }
00078     }
00079     return extraPoint;
00080 }
00081 
00082 // used to order according to angle, so it can remove double easily
00083 struct angleAndIndex
00084 {
00085     double angle;
00086     int index;
00087 };
00088 
00089 bool angleCompare( angleAndIndex lhs, angleAndIndex rhs )
00090 {
00091     return lhs.angle < rhs.angle;
00092 }
00093 
00094 // nP might be modified too, we will remove duplicates if found
00095 int IntxUtils::SortAndRemoveDoubles2( double* P, int& nP, double epsilon_1 )
00096 {
00097     if( nP < 2 ) return 0;  // nothing to do
00098 
00099     // center of gravity for the points
00100     double c[2] = { 0., 0. };
00101     int k       = 0;
00102     for( k = 0; k < nP; k++ )
00103     {
00104         c[0] += P[2 * k];
00105         c[1] += P[2 * k + 1];
00106     }
00107     c[0] /= nP;
00108     c[1] /= nP;
00109 
00110     // how many? we dimensioned P at MAXEDGES*10; so we imply we could have at most 5*MAXEDGES
00111     // intersection points
00112     struct angleAndIndex pairAngleIndex[5 * MAXEDGES];
00113 
00114     for( k = 0; k < nP; k++ )
00115     {
00116         double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
00117         if( x != 0. || y != 0. ) { pairAngleIndex[k].angle = atan2( y, x ); }
00118         else
00119         {
00120             pairAngleIndex[k].angle = 0;
00121             // it would mean that the cells are touching at a vertex
00122         }
00123         pairAngleIndex[k].index = k;
00124     }
00125 
00126     // this should be faster than the bubble sort we had before
00127     std::sort( pairAngleIndex, pairAngleIndex + nP, angleCompare );
00128     // copy now to a new double array
00129     double PCopy[10 * MAXEDGES];  // the same dimension as P; very conservative, but faster than
00130                                   // reallocate for a vector
00131     for( k = 0; k < nP; k++ )     // index will show where it should go now;
00132     {
00133         int ck           = pairAngleIndex[k].index;
00134         PCopy[2 * k]     = P[2 * ck];
00135         PCopy[2 * k + 1] = P[2 * ck + 1];
00136     }
00137     // now copy from PCopy over original P location
00138     std::copy( PCopy, PCopy + 2 * nP, P );
00139 
00140     // eliminate duplicates, finally
00141 
00142     int i = 0, j = 1;  // the next one; j may advance faster than i
00143     // check the unit
00144     // double epsilon_1 = 1.e-5; // these are cm; 2 points are the same if the distance is less
00145     // than 1.e-5 cm
00146     while( j < nP )
00147     {
00148         double d2 = dist2( &P[2 * i], &P[2 * j] );
00149         if( d2 > epsilon_1 )
00150         {
00151             i++;
00152             P[2 * i]     = P[2 * j];
00153             P[2 * i + 1] = P[2 * j + 1];
00154         }
00155         j++;
00156     }
00157     // test also the last point with the first one (index 0)
00158     // the first one could be at -PI; last one could be at +PI, according to atan2 span
00159 
00160     double d2 = dist2( P, &P[2 * i] );  // check the first and last points (ordered from -pi to +pi)
00161     if( d2 > epsilon_1 ) { nP = i + 1; }
00162     else
00163         nP = i;            // effectively delete the last point (that would have been the same with first)
00164     if( nP == 0 ) nP = 1;  // we should be left with at least one point we already tested if nP is 0 originally
00165     return 0;
00166 }
00167 
00168 // the marks will show what edges of blue intersect the red
00169 
00170 ErrorCode IntxUtils::EdgeIntersections2( double* blue, int nsBlue, double* red, int nsRed, int* markb, int* markr,
00171                                          double* points, int& nPoints )
00172 {
00173     /* EDGEINTERSECTIONS computes edge intersections of two elements
00174      [P,n]=EdgeIntersections(X,Y) computes for the two given elements  * red
00175      and blue ( stored column wise )
00176      (point coordinates are stored column-wise, in counter clock
00177      order) the points P where their edges intersect. In addition,
00178      in n the indices of which neighbors of red  are also intersecting
00179      with blue are given.
00180      */
00181 
00182     // points is an array with enough slots   (24 * 2 doubles)
00183     nPoints = 0;
00184     for( int i = 0; i < MAXEDGES; i++ )
00185     {
00186         markb[i] = markr[i] = 0;
00187     }
00188 
00189     for( int i = 0; i < nsBlue; i++ )
00190     {
00191         for( int j = 0; j < nsRed; j++ )
00192         {
00193             double b[2];
00194             double a[2][2];  // 2*2
00195             int iPlus1 = ( i + 1 ) % nsBlue;
00196             int jPlus1 = ( j + 1 ) % nsRed;
00197             for( int k = 0; k < 2; k++ )
00198             {
00199                 b[k] = red[2 * j + k] - blue[2 * i + k];
00200                 // row k of a: a(k, 0), a(k, 1)
00201                 a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
00202                 a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
00203             }
00204             double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
00205             if( fabs( delta ) > 1.e-14 )  // this is close to machine epsilon
00206             {
00207                 // not parallel
00208                 double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
00209                 double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
00210                 if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
00211                 {
00212                     // the intersection is good
00213                     for( int k = 0; k < 2; k++ )
00214                     {
00215                         points[2 * nPoints + k] = blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
00216                     }
00217                     markb[i] = 1;  // so neighbor number i of blue will be considered too.
00218                     markr[j] = 1;  // this will be used in advancing red around blue quad
00219                     nPoints++;
00220                 }
00221             }
00222             // the case delta ~ 0. will be considered by the interior points logic
00223         }
00224     }
00225     return MB_SUCCESS;
00226 }
00227 
00228 // special one, for intersection between rll (constant latitude)  and cs quads
00229 ErrorCode IntxUtils::EdgeIntxRllCs( double* blue, CartVect* bluec, int* blueEdgeType, int nsBlue, double* red,
00230                                     CartVect* redc, int nsRed, int* markb, int* markr, int plane, double R,
00231                                     double* points, int& nPoints )
00232 {
00233     // if blue edge type is 1, intersect in 3d then project to 2d by gnomonic projection
00234     // everything else the same (except if there are 2 points resulting, which is rare)
00235     for( int i = 0; i < 4; i++ )
00236     {  // always at most 4 , so maybe don't bother
00237         markb[i] = markr[i] = 0;
00238     }
00239 
00240     for( int i = 0; i < nsBlue; i++ )
00241     {
00242         int iPlus1 = ( i + 1 ) % nsBlue;
00243         if( blueEdgeType[i] == 0 )  // old style, just 2d
00244         {
00245             for( int j = 0; j < nsRed; j++ )
00246             {
00247                 double b[2];
00248                 double a[2][2];  // 2*2
00249 
00250                 int jPlus1 = ( j + 1 ) % nsRed;
00251                 for( int k = 0; k < 2; k++ )
00252                 {
00253                     b[k] = red[2 * j + k] - blue[2 * i + k];
00254                     // row k of a: a(k, 0), a(k, 1)
00255                     a[k][0] = blue[2 * iPlus1 + k] - blue[2 * i + k];
00256                     a[k][1] = red[2 * j + k] - red[2 * jPlus1 + k];
00257                 }
00258                 double delta = a[0][0] * a[1][1] - a[0][1] * a[1][0];
00259                 if( fabs( delta ) > 1.e-14 )
00260                 {
00261                     // not parallel
00262                     double alfa = ( b[0] * a[1][1] - a[0][1] * b[1] ) / delta;
00263                     double beta = ( -b[0] * a[1][0] + b[1] * a[0][0] ) / delta;
00264                     if( 0 <= alfa && alfa <= 1. && 0 <= beta && beta <= 1. )
00265                     {
00266                         // the intersection is good
00267                         for( int k = 0; k < 2; k++ )
00268                         {
00269                             points[2 * nPoints + k] =
00270                                 blue[2 * i + k] + alfa * ( blue[2 * iPlus1 + k] - blue[2 * i + k] );
00271                         }
00272                         markb[i] = 1;  // so neighbor number i of blue will be considered too.
00273                         markr[j] = 1;  // this will be used in advancing red around blue quad
00274                         nPoints++;
00275                     }
00276                 }  // if the edges are too "parallel", skip them
00277             }
00278         }
00279         else  // edge type is 1, so use 3d intersection
00280         {
00281             CartVect& C = bluec[i];
00282             CartVect& D = bluec[iPlus1];
00283             for( int j = 0; j < nsRed; j++ )
00284             {
00285                 int jPlus1  = ( j + 1 ) % nsRed;  // nsRed is just 4, forget about it, usually
00286                 CartVect& A = redc[j];
00287                 CartVect& B = redc[jPlus1];
00288                 int np      = 0;
00289                 double E[9];
00290                 intersect_great_circle_arc_with_clat_arc( A.array(), B.array(), C.array(), D.array(), R, E, np );
00291                 if( np == 0 ) continue;
00292                 if( np >= 2 ) { std::cout << "intersection with 2 points :" << A << B << C << D << "\n"; }
00293                 for( int k = 0; k < np; k++ )
00294                 {
00295                     gnomonic_projection( CartVect( E + k * 3 ), R, plane, points[2 * nPoints],
00296                                          points[2 * nPoints + 1] );
00297                     nPoints++;
00298                 }
00299                 markb[i] = 1;  // so neighbor number i of blue will be considered too.
00300                 markr[j] = 1;  // this will be used in advancing red around blue quad
00301             }
00302         }
00303     }
00304     return MB_SUCCESS;
00305 }
00306 
00307 // vec utils related to gnomonic projection on a sphere
00308 
00309 // vec utils
00310 
00311 /*
00312  *
00313  * position on a sphere of radius R
00314  * if plane specified, use it; if not, return the plane, and the point in the plane
00315  * there are 6 planes, numbered 1 to 6
00316  * plane 1: x=R, plane 2: y=R, 3: x=-R, 4: y=-R, 5: z=-R, 6: z=R
00317  *
00318  * projection on the plane will preserve the orientation, such that a triangle, quad pointing
00319  * outside the sphere will have a positive orientation in the projection plane
00320  */
00321 void IntxUtils::decide_gnomonic_plane( const CartVect& pos, int& plane )
00322 {
00323     // decide plane, based on max x, y, z
00324     if( fabs( pos[0] ) < fabs( pos[1] ) )
00325     {
00326         if( fabs( pos[2] ) < fabs( pos[1] ) )
00327         {
00328             // pos[1] is biggest
00329             if( pos[1] > 0 )
00330                 plane = 2;
00331             else
00332                 plane = 4;
00333         }
00334         else
00335         {
00336             // pos[2] is biggest
00337             if( pos[2] < 0 )
00338                 plane = 5;
00339             else
00340                 plane = 6;
00341         }
00342     }
00343     else
00344     {
00345         if( fabs( pos[2] ) < fabs( pos[0] ) )
00346         {
00347             // pos[0] is the greatest
00348             if( pos[0] > 0 )
00349                 plane = 1;
00350             else
00351                 plane = 3;
00352         }
00353         else
00354         {
00355             // pos[2] is biggest
00356             if( pos[2] < 0 )
00357                 plane = 5;
00358             else
00359                 plane = 6;
00360         }
00361     }
00362     return;
00363 }
00364 
00365 // point on a sphere is projected on one of six planes, decided earlier
00366 ErrorCode IntxUtils::gnomonic_projection( const CartVect& pos, double R, int plane, double& c1, double& c2 )
00367 {
00368     double alfa = 1.;  // the new point will be on line alfa*pos
00369 
00370     switch( plane )
00371     {
00372         case 1: {
00373             // the plane with x = R; x>y, x>z
00374             // c1->y, c2->z
00375             alfa = R / pos[0];
00376             c1   = alfa * pos[1];
00377             c2   = alfa * pos[2];
00378             break;
00379         }
00380         case 2: {
00381             // y = R -> zx
00382             alfa = R / pos[1];
00383             c1   = alfa * pos[2];
00384             c2   = alfa * pos[0];
00385             break;
00386         }
00387         case 3: {
00388             // x=-R, -> yz
00389             alfa = -R / pos[0];
00390             c1   = -alfa * pos[1];  // the sign is to preserve orientation
00391             c2   = alfa * pos[2];
00392             break;
00393         }
00394         case 4: {
00395             // y = -R
00396             alfa = -R / pos[1];
00397             c1   = -alfa * pos[2];  // the sign is to preserve orientation
00398             c2   = alfa * pos[0];
00399             break;
00400         }
00401         case 5: {
00402             // z = -R
00403             alfa = -R / pos[2];
00404             c1   = -alfa * pos[0];  // the sign is to preserve orientation
00405             c2   = alfa * pos[1];
00406             break;
00407         }
00408         case 6: {
00409             alfa = R / pos[2];
00410             c1   = alfa * pos[0];
00411             c2   = alfa * pos[1];
00412             break;
00413         }
00414         default:
00415             return MB_FAILURE;  // error
00416     }
00417 
00418     return MB_SUCCESS;  // no error
00419 }
00420 
00421 // given the position on plane (one out of 6), find out the position on sphere
00422 ErrorCode IntxUtils::reverse_gnomonic_projection( const double& c1, const double& c2, double R, int plane,
00423                                                   CartVect& pos )
00424 {
00425 
00426     // the new point will be on line beta*pos
00427     double len  = sqrt( c1 * c1 + c2 * c2 + R * R );
00428     double beta = R / len;  // it is less than 1, in general
00429 
00430     switch( plane )
00431     {
00432         case 1: {
00433             // the plane with x = R; x>y, x>z
00434             // c1->y, c2->z
00435             pos[0] = beta * R;
00436             pos[1] = c1 * beta;
00437             pos[2] = c2 * beta;
00438             break;
00439         }
00440         case 2: {
00441             // y = R -> zx
00442             pos[1] = R * beta;
00443             pos[2] = c1 * beta;
00444             pos[0] = c2 * beta;
00445             break;
00446         }
00447         case 3: {
00448             // x=-R, -> yz
00449             pos[0] = -R * beta;
00450             pos[1] = -c1 * beta;  // the sign is to preserve orientation
00451             pos[2] = c2 * beta;
00452             break;
00453         }
00454         case 4: {
00455             // y = -R
00456             pos[1] = -R * beta;
00457             pos[2] = -c1 * beta;  // the sign is to preserve orientation
00458             pos[0] = c2 * beta;
00459             break;
00460         }
00461         case 5: {
00462             // z = -R
00463             pos[2] = -R * beta;
00464             pos[0] = -c1 * beta;  // the sign is to preserve orientation
00465             pos[1] = c2 * beta;
00466             break;
00467         }
00468         case 6: {
00469             pos[2] = R * beta;
00470             pos[0] = c1 * beta;
00471             pos[1] = c2 * beta;
00472             break;
00473         }
00474         default:
00475             return MB_FAILURE;  // error
00476     }
00477 
00478     return MB_SUCCESS;  // no error
00479 }
00480 
00481 void IntxUtils::gnomonic_unroll( double& c1, double& c2, double R, int plane )
00482 {
00483     double tmp;
00484     switch( plane )
00485     {
00486         case 1:
00487             break;  // nothing
00488         case 2:     // rotate + 90
00489             tmp = c1;
00490             c1  = -c2;
00491             c2  = tmp;
00492             c1  = c1 + 2 * R;
00493             break;
00494         case 3:
00495             c1 = c1 + 4 * R;
00496             break;
00497         case 4:  // rotate with -90 x-> -y; y -> x
00498 
00499             tmp = c1;
00500             c1  = c2;
00501             c2  = -tmp;
00502             c1  = c1 - 2 * R;
00503             break;
00504         case 5:  // South Pole
00505             // rotate 180 then move to (-2, -2)
00506             c1 = -c1 - 2. * R;
00507             c2 = -c2 - 2. * R;
00508             break;
00509             ;
00510         case 6:  // North Pole
00511             c1 = c1 - 2. * R;
00512             c2 = c2 + 2. * R;
00513             break;
00514     }
00515     return;
00516 }
00517 // given a mesh on the sphere, project all centers in 6 gnomonic planes, or project mesh too
00518 ErrorCode IntxUtils::global_gnomonic_projection( Interface* mb, EntityHandle inSet, double R, bool centers_only,
00519                                                  EntityHandle& outSet )
00520 {
00521     std::string parTagName( "PARALLEL_PARTITION" );
00522     Tag part_tag;
00523     Range partSets;
00524     ErrorCode rval = mb->tag_get_handle( parTagName.c_str(), part_tag );
00525     if( MB_SUCCESS == rval && part_tag != 0 )
00526     {
00527         rval = mb->get_entities_by_type_and_tag( inSet, MBENTITYSET, &part_tag, NULL, 1, partSets, Interface::UNION );MB_CHK_ERR( rval );
00528     }
00529     rval = ScaleToRadius( mb, inSet, 1.0 );MB_CHK_ERR( rval );
00530     // Get all entities of dimension 2
00531     Range inputRange;  // get
00532     rval = mb->get_entities_by_dimension( inSet, 1, inputRange );MB_CHK_ERR( rval );
00533     rval = mb->get_entities_by_dimension( inSet, 2, inputRange );MB_CHK_ERR( rval );
00534 
00535     std::map< EntityHandle, int > partsAssign;
00536     std::map< int, EntityHandle > newPartSets;
00537     if( !partSets.empty() )
00538     {
00539         // get all cells, and assign parts
00540         for( Range::iterator setIt = partSets.begin(); setIt != partSets.end(); ++setIt )
00541         {
00542             EntityHandle pSet = *setIt;
00543             Range ents;
00544             rval = mb->get_entities_by_handle( pSet, ents );MB_CHK_ERR( rval );
00545             int val;
00546             rval = mb->tag_get_data( part_tag, &pSet, 1, &val );MB_CHK_ERR( rval );
00547             // create a new set with the same part id tag, in the outSet
00548             EntityHandle newPartSet;
00549             rval = mb->create_meshset( MESHSET_SET, newPartSet );MB_CHK_ERR( rval );
00550             rval = mb->tag_set_data( part_tag, &newPartSet, 1, &val );MB_CHK_ERR( rval );
00551             newPartSets[val] = newPartSet;
00552             rval             = mb->add_entities( outSet, &newPartSet, 1 );MB_CHK_ERR( rval );
00553             for( Range::iterator it = ents.begin(); it != ents.end(); ++it )
00554             {
00555                 partsAssign[*it] = val;
00556             }
00557         }
00558     }
00559 
00560     if( centers_only )
00561     {
00562         for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
00563         {
00564             CartVect center;
00565             EntityHandle cell = *it;
00566             rval              = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
00567             int plane;
00568             decide_gnomonic_plane( center, plane );
00569             double c[3];
00570             c[2] = 0.;
00571             gnomonic_projection( center, R, plane, c[0], c[1] );
00572 
00573             gnomonic_unroll( c[0], c[1], R, plane );
00574 
00575             EntityHandle vertex;
00576             rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
00577             rval = mb->add_entities( outSet, &vertex, 1 );MB_CHK_ERR( rval );
00578         }
00579     }
00580     else
00581     {
00582         // distribute the cells to 6 planes, based on the center
00583         Range subranges[6];
00584         for( Range::iterator it = inputRange.begin(); it != inputRange.end(); ++it )
00585         {
00586             CartVect center;
00587             EntityHandle cell = *it;
00588             rval              = mb->get_coords( &cell, 1, center.array() );MB_CHK_ERR( rval );
00589             int plane;
00590             decide_gnomonic_plane( center, plane );
00591             subranges[plane - 1].insert( cell );
00592         }
00593         for( int i = 1; i <= 6; i++ )
00594         {
00595             Range verts;
00596             rval = mb->get_connectivity( subranges[i - 1], verts );MB_CHK_ERR( rval );
00597             std::map< EntityHandle, EntityHandle > corr;
00598             for( Range::iterator vt = verts.begin(); vt != verts.end(); ++vt )
00599             {
00600                 CartVect vect;
00601                 EntityHandle v = *vt;
00602                 rval           = mb->get_coords( &v, 1, vect.array() );MB_CHK_ERR( rval );
00603                 double c[3];
00604                 c[2] = 0.;
00605                 gnomonic_projection( vect, R, i, c[0], c[1] );
00606                 gnomonic_unroll( c[0], c[1], R, i );
00607                 EntityHandle vertex;
00608                 rval = mb->create_vertex( c, vertex );MB_CHK_ERR( rval );
00609                 corr[v] = vertex;  // for new connectivity
00610             }
00611             EntityHandle new_conn[20];  // max edges in 2d ?
00612             for( Range::iterator eit = subranges[i - 1].begin(); eit != subranges[i - 1].end(); ++eit )
00613             {
00614                 EntityHandle eh          = *eit;
00615                 const EntityHandle* conn = NULL;
00616                 int num_nodes;
00617                 rval = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
00618                 // build a new vertex array
00619                 for( int i = 0; i < num_nodes; i++ )
00620                     new_conn[i] = corr[conn[i]];
00621                 EntityType type = mb->type_from_handle( eh );
00622                 EntityHandle newCell;
00623                 rval = mb->create_element( type, new_conn, num_nodes, newCell );MB_CHK_ERR( rval );
00624                 rval = mb->add_entities( outSet, &newCell, 1 );MB_CHK_ERR( rval );
00625                 std::map< EntityHandle, int >::iterator mit = partsAssign.find( eh );
00626                 if( mit != partsAssign.end() )
00627                 {
00628                     int val = mit->second;
00629                     rval    = mb->add_entities( newPartSets[val], &newCell, 1 );MB_CHK_ERR( rval );
00630                 }
00631             }
00632         }
00633     }
00634 
00635     return MB_SUCCESS;
00636 }
00637 void IntxUtils::transform_coordinates( double* avg_position, int projection_type )
00638 {
00639     if( projection_type == 1 )
00640     {
00641         double R =
00642             avg_position[0] * avg_position[0] + avg_position[1] * avg_position[1] + avg_position[2] * avg_position[2];
00643         R               = sqrt( R );
00644         double lat      = asin( avg_position[2] / R );
00645         double lon      = atan2( avg_position[1], avg_position[0] );
00646         avg_position[0] = lon;
00647         avg_position[1] = lat;
00648         avg_position[2] = R;
00649     }
00650     else if( projection_type == 2 )  // gnomonic projection
00651     {
00652         CartVect pos( avg_position );
00653         int gplane;
00654         IntxUtils::decide_gnomonic_plane( pos, gplane );
00655 
00656         IntxUtils::gnomonic_projection( pos, 1.0, gplane, avg_position[0], avg_position[1] );
00657         avg_position[2] = 0;
00658         IntxUtils::gnomonic_unroll( avg_position[0], avg_position[1], 1.0, gplane );
00659     }
00660 }
00661 /*
00662  *
00663  use physical_constants, only : dd_pi
00664  type(cartesian3D_t), intent(in) :: cart
00665  type(spherical_polar_t)         :: sphere
00666 
00667  sphere%r=distance(cart)
00668  sphere%lat=ASIN(cart%z/sphere%r)
00669  sphere%lon=0
00670 
00671  ! ==========================================================
00672  ! enforce three facts:
00673  !
00674  ! 1) lon at poles is defined to be zero
00675  !
00676  ! 2) Grid points must be separated by about .01 Meter (on earth)
00677  !    from pole to be considered "not the pole".
00678  !
00679  ! 3) range of lon is { 0<= lon < 2*pi }
00680  !
00681  ! ==========================================================
00682 
00683  if (distance(cart) >= DIST_THRESHOLD) then
00684  sphere%lon=ATAN2(cart%y,cart%x)
00685  if (sphere%lon<0) then
00686  sphere%lon=sphere%lon+2*DD_PI
00687  end if
00688  end if
00689 
00690  end function cart_to_spherical
00691  */
00692 IntxUtils::SphereCoords IntxUtils::cart_to_spherical( CartVect& cart3d )
00693 {
00694     SphereCoords res;
00695     res.R = cart3d.length();
00696     if( res.R < 0 )
00697     {
00698         res.lon = res.lat = 0.;
00699         return res;
00700     }
00701     res.lat = asin( cart3d[2] / res.R );
00702     res.lon = atan2( cart3d[1], cart3d[0] );
00703     if( res.lon < 0 ) res.lon += 2 * M_PI;  // M_PI is defined in math.h? it seems to be true, although
00704     // there are some defines it depends on :(
00705     // #if defined __USE_BSD || defined __USE_XOPEN ???
00706 
00707     return res;
00708 }
00709 
00710 /*
00711  * ! ===================================================================
00712  ! spherical_to_cart:
00713  ! converts spherical polar {lon,lat}  to 3D cartesian {x,y,z}
00714  ! on unit sphere
00715  ! ===================================================================
00716 
00717  function spherical_to_cart(sphere) result (cart)
00718 
00719  type(spherical_polar_t), intent(in) :: sphere
00720  type(cartesian3D_t)                 :: cart
00721 
00722  cart%x=sphere%r*COS(sphere%lat)*COS(sphere%lon)
00723  cart%y=sphere%r*COS(sphere%lat)*SIN(sphere%lon)
00724  cart%z=sphere%r*SIN(sphere%lat)
00725 
00726  end function spherical_to_cart
00727  */
00728 CartVect IntxUtils::spherical_to_cart( IntxUtils::SphereCoords& sc )
00729 {
00730     CartVect res;
00731     res[0] = sc.R * cos( sc.lat ) * cos( sc.lon );  // x coordinate
00732     res[1] = sc.R * cos( sc.lat ) * sin( sc.lon );  // y
00733     res[2] = sc.R * sin( sc.lat );                  // z
00734     return res;
00735 }
00736 
00737 ErrorCode IntxUtils::ScaleToRadius( Interface* mb, EntityHandle set, double R )
00738 {
00739     Range nodes;
00740     ErrorCode rval = mb->get_entities_by_type( set, MBVERTEX, nodes, true );  // recursive
00741     if( rval != moab::MB_SUCCESS ) return rval;
00742 
00743     // one by one, get the node and project it on the sphere, with a radius given
00744     // the center of the sphere is at 0,0,0
00745     for( Range::iterator nit = nodes.begin(); nit != nodes.end(); ++nit )
00746     {
00747         EntityHandle nd = *nit;
00748         CartVect pos;
00749         rval = mb->get_coords( &nd, 1, (double*)&( pos[0] ) );
00750         if( rval != moab::MB_SUCCESS ) return rval;
00751         double len = pos.length();
00752         if( len == 0. ) return MB_FAILURE;
00753         pos  = R / len * pos;
00754         rval = mb->set_coords( &nd, 1, (double*)&( pos[0] ) );
00755         if( rval != moab::MB_SUCCESS ) return rval;
00756     }
00757     return MB_SUCCESS;
00758 }
00759 
00760 // assume they are one the same sphere
00761 double IntxAreaUtils::spherical_angle( double* A, double* B, double* C, double Radius )
00762 {
00763     // the angle by definition is between the planes OAB and OBC
00764     CartVect a( A );
00765     CartVect b( B );
00766     CartVect c( C );
00767     double err1 = a.length_squared() - Radius * Radius;
00768     if( fabs( err1 ) > 0.0001 )
00769     { std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n"; }
00770     CartVect normalOAB = a * b;
00771     CartVect normalOCB = c * b;
00772     return angle( normalOAB, normalOCB );
00773 }
00774 
00775 // could be bigger than M_PI;
00776 // angle at B could be bigger than M_PI, if the orientation is such that ABC points toward the
00777 // interior
00778 double IntxUtils::oriented_spherical_angle( double* A, double* B, double* C )
00779 {
00780     // assume the same radius, sphere at origin
00781     CartVect a( A ), b( B ), c( C );
00782     CartVect normalOAB = a * b;
00783     CartVect normalOCB = c * b;
00784     CartVect orient    = ( c - b ) * ( a - b );
00785     double ang         = angle( normalOAB, normalOCB );  // this is between 0 and M_PI
00786     if( ang != ang )
00787     {
00788         // signal of a nan
00789         std::cout << a << " " << b << " " << c << "\n";
00790         std::cout << ang << "\n";
00791     }
00792     if( orient % b < 0 ) return ( 2 * M_PI - ang );  // the other angle, supplement
00793 
00794     return ang;
00795 }
00796 
00797 double IntxAreaUtils::area_spherical_triangle( double* A, double* B, double* C, double Radius )
00798 {
00799     switch( m_eAreaMethod )
00800     {
00801         case Girard:
00802             return area_spherical_triangle_girard( A, B, C, Radius );
00803 #ifdef MOAB_HAVE_TEMPESTREMAP
00804         case GaussQuadrature:
00805             return area_spherical_triangle_GQ( A, B, C, Radius );
00806 #endif
00807         case lHuiller:
00808         default:
00809             return area_spherical_triangle_lHuiller( A, B, C, Radius );
00810     }
00811 }
00812 
00813 double IntxAreaUtils::area_spherical_polygon( double* A, int N, double Radius, int* sign )
00814 {
00815     switch( m_eAreaMethod )
00816     {
00817         case Girard:
00818             return area_spherical_polygon_girard( A, N, Radius );
00819 #ifdef MOAB_HAVE_TEMPESTREMAP
00820         case GaussQuadrature:
00821             return area_spherical_polygon_GQ( A, N, Radius );
00822 #endif
00823         case lHuiller:
00824         default:
00825             return area_spherical_polygon_lHuiller( A, N, Radius, sign );
00826     }
00827 }
00828 
00829 double IntxAreaUtils::area_spherical_triangle_girard( double* A, double* B, double* C, double Radius )
00830 {
00831     double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
00832                         spherical_angle( C, A, B, Radius ) - M_PI;
00833     double area = Radius * Radius * correction;
00834     // now, is it negative or positive? is it pointing toward the center or outward?
00835     CartVect a( A ), b( B ), c( C );
00836     CartVect abc = ( b - a ) * ( c - a );
00837     if( abc % a > 0 )  // dot product positive, means ABC points out
00838         return area;
00839     else
00840         return -area;
00841 }
00842 
00843 double IntxAreaUtils::area_spherical_polygon_girard( double* A, int N, double Radius )
00844 {
00845     // this should work for non-convex polygons too
00846     // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00847     //
00848     if( N <= 2 ) return 0.;
00849     double sum_angles = 0.;
00850     for( int i = 0; i < N; i++ )
00851     {
00852         int i1 = ( i + 1 ) % N;
00853         int i2 = ( i + 2 ) % N;
00854         sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
00855     }
00856     double correction = sum_angles - ( N - 2 ) * M_PI;
00857     return Radius * Radius * correction;
00858 }
00859 
00860 double IntxAreaUtils::area_spherical_polygon_lHuiller( double* A, int N, double Radius, int* sign )
00861 {
00862     // This should work for non-convex polygons too
00863     // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00864     // We also assume that the orientation is positive;
00865     // If negative orientation, the area will be negative
00866     if( N <= 2 ) return 0.;
00867 
00868     int lsign   = 1;  // assume positive orientain
00869     double area = 0.;
00870     for( int i = 1; i < N - 1; i++ )
00871     {
00872         int i1              = i + 1;
00873         double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
00874         if( areaTriangle < 0 )
00875             lsign = -1;  // signal that we have at least one triangle with negative orientation ;
00876                          // possible nonconvex polygon
00877         area += areaTriangle;
00878     }
00879     if( sign ) *sign = lsign;
00880 
00881     return area;
00882 }
00883 
00884 #ifdef MOAB_HAVE_TEMPESTREMAP
00885 double IntxAreaUtils::area_spherical_polygon_GQ( double* A, int N, double Radius )
00886 {
00887     // this should work for non-convex polygons too
00888     // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00889     // We also assume that the orientation is positive;
00890     // If negative orientation, the area can be negative
00891     if( N <= 2 ) return 0.;
00892 
00893     // assume positive orientain
00894     double area = 0.;
00895     for( int i = 1; i < N - 1; i++ )
00896     {
00897         int i1 = i + 1;
00898         area += area_spherical_triangle_GQ( A, A + 3 * i, A + 3 * i1, Radius );
00899     }
00900     return area;
00901 }
00902 
00903 /* compute the area by using Gauss-Quadratures; use TR interfaces directly */
00904 double IntxAreaUtils::area_spherical_triangle_GQ( double* ptA, double* ptB, double* ptC, double )
00905 {
00906     Face face( 3 );
00907     NodeVector nodes( 3 );
00908     nodes[0] = Node( ptA[0], ptA[1], ptA[2] );
00909     nodes[1] = Node( ptB[0], ptB[1], ptB[2] );
00910     nodes[2] = Node( ptC[0], ptC[1], ptC[2] );
00911     face.SetNode( 0, 0 );
00912     face.SetNode( 1, 1 );
00913     face.SetNode( 2, 2 );
00914     return CalculateFaceArea( face, nodes );
00915 }
00916 #endif
00917 
00918 /*
00919  *  l'Huiller's formula for spherical triangle
00920  *  http://williams.best.vwh.net/avform.htm
00921  *  a, b, c are arc measures in radians, too
00922  *  A, B, C are angles on the sphere, for which we already have formula
00923  *               c
00924  *         A -------B
00925  *          \       |
00926  *           \      |
00927  *            \b    |a
00928  *             \    |
00929  *              \   |
00930  *               \  |
00931  *                \C|
00932  *                 \|
00933  *
00934  *  (The angle at B is not necessarily a right angle)
00935  *
00936  *    sin(a)  sin(b)   sin(c)
00937  *    ----- = ------ = ------
00938  *    sin(A)  sin(B)   sin(C)
00939  *
00940  * In terms of the sides (this is excess, as before, but numerically stable)
00941  *
00942  *  E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
00943  */
00944 double IntxAreaUtils::area_spherical_triangle_lHuiller( double* ptA, double* ptB, double* ptC, double Radius )
00945 {
00946 
00947     // now, a is angle BOC, O is origin
00948     CartVect vA( ptA ), vB( ptB ), vC( ptC );
00949     double a = angle( vB, vC );
00950     double b = angle( vC, vA );
00951     double c = angle( vA, vB );
00952     int sign = 1;
00953     if( ( vA * vB ) % vC < 0 ) sign = -1;
00954     double s   = ( a + b + c ) / 2;
00955     double tmp = tan( s / 2 ) * tan( ( s - a ) / 2 ) * tan( ( s - b ) / 2 ) * tan( ( s - c ) / 2 );
00956     if( tmp < 0. ) tmp = 0.;
00957 
00958     double E = 4 * atan( sqrt( tmp ) );
00959     if( E != E ) std::cout << " NaN at spherical triangle area \n";
00960 
00961     double area = sign * E * Radius * Radius;
00962 
00963 #ifdef CHECKNEGATIVEAREA
00964     if( area < 0 )
00965     {
00966         std::cout << "negative area: " << area << "\n";
00967         std::cout << std::setprecision( 15 );
00968         std::cout << "vA: " << vA << "\n";
00969         std::cout << "vB: " << vB << "\n";
00970         std::cout << "vC: " << vC << "\n";
00971         std::cout << "sign: " << sign << "\n";
00972         std::cout << " a: " << a << "\n";
00973         std::cout << " b: " << b << "\n";
00974         std::cout << " c: " << c << "\n";
00975     }
00976 #endif
00977 
00978     return area;
00979 }
00980 #undef CHECKNEGATIVEAREA
00981 
00982 double IntxAreaUtils::area_on_sphere( Interface* mb, EntityHandle set, double R )
00983 {
00984     // Get all entities of dimension 2
00985     Range inputRange;
00986     ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
00987 
00988     // Filter by elements that are owned by current process
00989     std::vector< int > ownerinfo( inputRange.size(), -1 );
00990     Tag intxOwnerTag;
00991     rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
00992     if( MB_SUCCESS == rval )
00993     {
00994         rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
00995     }
00996 
00997     // compare total area with 4*M_PI * R^2
00998     int ie            = 0;
00999     double total_area = 0.;
01000     for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
01001     {
01002 
01003         // All zero/positive owner data represents ghosted elems
01004         if( ownerinfo[ie++] >= 0 ) continue;
01005 
01006         EntityHandle eh        = *eit;
01007         const double elem_area = this->area_spherical_element( mb, eh, R );
01008 
01009         // check whether the area of the spherical element is positive.
01010         assert( elem_area > 0 );
01011 
01012         // sum up the contribution
01013         total_area += elem_area;
01014     }
01015 
01016     // return total mesh area
01017     return total_area;
01018 }
01019 
01020 double IntxAreaUtils::area_spherical_element( Interface* mb, EntityHandle elem, double R )
01021 {
01022     // get the nodes, then the coordinates
01023     const EntityHandle* verts;
01024     int nsides;
01025     ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
01026 
01027     // account for possible padded polygons
01028     while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
01029         nsides--;
01030 
01031     // get coordinates
01032     std::vector< double > coords( 3 * nsides );
01033     rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
01034 
01035     // compute and return the area of the polygonal element
01036     return area_spherical_polygon( &coords[0], nsides, R );
01037 }
01038 
01039 double IntxUtils::distance_on_great_circle( CartVect& p1, CartVect& p2 )
01040 {
01041     SphereCoords sph1 = cart_to_spherical( p1 );
01042     SphereCoords sph2 = cart_to_spherical( p2 );
01043     // radius should be the same
01044     return sph1.R *
01045            acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
01046 }
01047 
01048 // break the nonconvex quads into triangles; remove the quad from the set? yes.
01049 // maybe radius is not needed;
01050 //
01051 ErrorCode IntxUtils::enforce_convexity( Interface* mb, EntityHandle lset, int my_rank )
01052 {
01053     // look at each polygon; compute all angles; if one is reflex, break that angle with
01054     // the next triangle; put the 2 new polys in the set;
01055     // still look at the next poly
01056     // replace it with 2 triangles, and remove from set;
01057     // it should work for all polygons / tested first for case 1, with dt 0.5 (too much deformation)
01058     // get all entities of dimension 2
01059     // then get the connectivity, etc
01060 
01061     Range inputRange;
01062     ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
01063 
01064     Tag corrTag       = 0;
01065     EntityHandle dumH = 0;
01066     rval              = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
01067     if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
01068 
01069     Tag gidTag;
01070     rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
01071 
01072     std::vector< double > coords;
01073     coords.resize( 3 * MAXEDGES );  // at most 10 vertices per polygon
01074     // we should create a queue with new polygons that need processing for reflex angles
01075     //  (obtuse)
01076     std::queue< EntityHandle > newPolys;
01077     int brokenPolys     = 0;
01078     Range::iterator eit = inputRange.begin();
01079     while( eit != inputRange.end() || !newPolys.empty() )
01080     {
01081         EntityHandle eh;
01082         if( eit != inputRange.end() )
01083         {
01084             eh = *eit;
01085             ++eit;
01086         }
01087         else
01088         {
01089             eh = newPolys.front();
01090             newPolys.pop();
01091         }
01092         // get the nodes, then the coordinates
01093         const EntityHandle* verts;
01094         int num_nodes;
01095         rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
01096         int nsides = num_nodes;
01097         // account for possible padded polygons
01098         while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
01099             nsides--;
01100         EntityHandle corrHandle = 0;
01101         if( corrTag )
01102         {
01103             rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
01104         }
01105         int gid = 0;
01106         rval    = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
01107         coords.resize( 3 * nsides );
01108         if( nsides < 4 ) continue;  // if already triangles, don't bother
01109         // get coordinates
01110         rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
01111         // compute each angle
01112         bool alreadyBroken = false;
01113 
01114         for( int i = 0; i < nsides; i++ )
01115         {
01116             double* A    = &coords[3 * i];
01117             double* B    = &coords[3 * ( ( i + 1 ) % nsides )];
01118             double* C    = &coords[3 * ( ( i + 2 ) % nsides )];
01119             double angle = IntxUtils::oriented_spherical_angle( A, B, C );
01120             if( angle - M_PI > 0. )  // even almost reflex is bad; break it!
01121             {
01122                 if( alreadyBroken )
01123                 {
01124                     mb->list_entities( &eh, 1 );
01125                     mb->list_entities( verts, nsides );
01126                     double* D = &coords[3 * ( ( i + 3 ) % nsides )];
01127                     std::cout << "ABC: " << angle << " \n";
01128                     std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
01129                     std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
01130                     std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
01131                     std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
01132 
01133                     return MB_FAILURE;
01134                 }
01135                 // the bad angle is at i+1;
01136                 // create 1 triangle and one polygon; add the polygon to the input range, so
01137                 // it will be processed too
01138                 // also, add both to the set :) and remove the original polygon from the set
01139                 // break the next triangle, even though not optimal
01140                 // so create the triangle i+1, i+2, i+3; remove i+2 from original list
01141                 // even though not optimal in general, it is good enough.
01142                 EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
01143                                           verts[( i + 3 ) % nsides] };
01144                 // create a polygon with num_nodes-1 vertices, and connectivity
01145                 // verts[i+1], verts[i+3], (all except i+2)
01146                 std::vector< EntityHandle > conn( nsides - 1 );
01147                 for( int j = 1; j < nsides; j++ )
01148                 {
01149                     conn[j - 1] = verts[( i + j + 2 ) % nsides];
01150                 }
01151                 EntityHandle newElement;
01152                 rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
01153 
01154                 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
01155                 if( corrTag )
01156                 {
01157                     rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
01158                 }
01159                 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
01160                 if( nsides == 4 )
01161                 {
01162                     // create another triangle
01163                     rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
01164                 }
01165                 else
01166                 {
01167                     // create another polygon, and add it to the inputRange
01168                     rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
01169                     newPolys.push( newElement );  // because it has less number of edges, the
01170                     // reverse should work to find it.
01171                 }
01172                 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
01173                 if( corrTag )
01174                 {
01175                     rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
01176                 }
01177                 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
01178                 rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
01179                 brokenPolys++;
01180                 alreadyBroken = true;  // get out of the loop, element is broken
01181             }
01182         }
01183     }
01184     if (brokenPolys > 0)
01185     {
01186         std::cout << "on local process " << my_rank << ", " << brokenPolys
01187                   << " concave polygons were decomposed in convex ones \n";
01188 #ifdef VERBOSE
01189         std::stringstream fff;
01190         fff << "file_set" <<  mb->id_from_handle(lset) << "rk_"<< my_rank << ".h5m";
01191         rval = mb->write_file(fff.str().c_str(), 0, 0, &lset, 1);MB_CHK_ERR( rval );
01192         std::cout << "wrote new file set: " << fff.str() << "\n";
01193 #endif
01194     }
01195     return MB_SUCCESS;
01196 }
01197 
01198 // looking at quad connectivity, collapse to triangle if 2 nodes equal
01199 // then delete the old quad
01200 ErrorCode IntxUtils::fix_degenerate_quads( Interface* mb, EntityHandle set )
01201 {
01202     Range quads;
01203     ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
01204     Tag gid;
01205     gid = mb->globalId_tag();
01206     for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
01207     {
01208         EntityHandle quad         = *qit;
01209         const EntityHandle* conn4 = NULL;
01210         int num_nodes             = 0;
01211         rval                      = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
01212         for( int i = 0; i < num_nodes; i++ )
01213         {
01214             int next_node_index = ( i + 1 ) % num_nodes;
01215             if( conn4[i] == conn4[next_node_index] )
01216             {
01217                 // form a triangle and delete the quad
01218                 // first get the global id, to set it on triangle later
01219                 int global_id = 0;
01220                 rval          = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
01221                 int i2                = ( i + 2 ) % num_nodes;
01222                 int i3                = ( i + 3 ) % num_nodes;
01223                 EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
01224                 EntityHandle tri;
01225                 rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
01226                 mb->add_entities( set, &tri, 1 );
01227                 mb->remove_entities( set, &quad, 1 );
01228                 mb->delete_entities( &quad, 1 );
01229                 rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
01230             }
01231         }
01232     }
01233     return MB_SUCCESS;
01234 }
01235 
01236 ErrorCode IntxAreaUtils::positive_orientation( Interface* mb, EntityHandle set, double R )
01237 {
01238     Range cells2d;
01239     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );
01240     if( MB_SUCCESS != rval ) return rval;
01241     for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
01242     {
01243         EntityHandle cell        = *qit;
01244         const EntityHandle* conn = NULL;
01245         int num_nodes            = 0;
01246         rval                     = mb->get_connectivity( cell, conn, num_nodes );
01247         if( MB_SUCCESS != rval ) return rval;
01248         if( num_nodes < 3 ) return MB_FAILURE;
01249 
01250         double coords[9];
01251         rval = mb->get_coords( conn, 3, coords );
01252         if( MB_SUCCESS != rval ) return rval;
01253 
01254         double area;
01255         if( R > 0 )
01256             area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
01257         else
01258             area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
01259         if( area < 0 )
01260         {
01261             // compute all area, do not revert if total area is positive
01262             std::vector< double > coords2( 3 * num_nodes );
01263             // get coordinates
01264             rval = mb->get_coords( conn, num_nodes, &coords2[0] );
01265             if( MB_SUCCESS != rval ) return MB_FAILURE;
01266             double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
01267             if( totArea < 0 )
01268             {
01269                 std::vector< EntityHandle > newconn( num_nodes );
01270                 for( int i = 0; i < num_nodes; i++ )
01271                 {
01272                     newconn[num_nodes - 1 - i] = conn[i];
01273                 }
01274                 rval = mb->set_connectivity( cell, &newconn[0], num_nodes );
01275                 if( MB_SUCCESS != rval ) return rval;
01276             }
01277             else
01278             {
01279                 std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
01280             }
01281         }
01282     }
01283     return MB_SUCCESS;
01284 }
01285 
01286 // distance along a great circle on a sphere of radius 1
01287 // page 4
01288 double IntxUtils::distance_on_sphere( double la1, double te1, double la2, double te2 )
01289 {
01290     return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
01291 }
01292 
01293 /*
01294  * given 2 great circle arcs, AB and CD, compute the unique intersection point, if it exists
01295  *  in between
01296  */
01297 ErrorCode IntxUtils::intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E )
01298 {
01299     // first verify A, B, C, D are on the same sphere
01300     double R2              = R * R;
01301     const double Tolerance = 1.e-12 * R2;
01302 
01303     CartVect a( A ), b( B ), c( C ), d( D );
01304 
01305     if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
01306             fabs( d.length_squared() - R2 ) >
01307         10 * Tolerance )
01308         return MB_FAILURE;
01309 
01310     CartVect n1 = a * b;
01311     if( n1.length_squared() < Tolerance ) return MB_FAILURE;
01312 
01313     CartVect n2 = c * d;
01314     if( n2.length_squared() < Tolerance ) return MB_FAILURE;
01315     CartVect n3 = n1 * n2;
01316     n3.normalize();
01317 
01318     n3 = R * n3;
01319     // the intersection is either n3 or -n3
01320     CartVect n4 = a * n3, n5 = n3 * b;
01321     if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
01322     {
01323         // n3 is good for ab, see if it is good for cd
01324         n4 = c * n3;
01325         n5 = n3 * d;
01326         if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
01327         {
01328             E[0] = n3[0];
01329             E[1] = n3[1];
01330             E[2] = n3[2];
01331         }
01332         else
01333             return MB_FAILURE;
01334     }
01335     else
01336     {
01337         // try -n3
01338         n3 = -n3;
01339         n4 = a * n3, n5 = n3 * b;
01340         if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
01341         {
01342             // n3 is good for ab, see if it is good for cd
01343             n4 = c * n3;
01344             n5 = n3 * d;
01345             if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
01346             {
01347                 E[0] = n3[0];
01348                 E[1] = n3[1];
01349                 E[2] = n3[2];
01350             }
01351             else
01352                 return MB_FAILURE;
01353         }
01354         else
01355             return MB_FAILURE;
01356     }
01357 
01358     return MB_SUCCESS;
01359 }
01360 
01361 // verify that result is in between a and b on a great circle arc, and between c and d on a constant
01362 // latitude arc
01363 static bool verify( CartVect a, CartVect b, CartVect c, CartVect d, double x, double y, double z )
01364 {
01365     // to check, the point has to be between a and b on a great arc, and between c and d on a const
01366     // lat circle
01367     CartVect s( x, y, z );
01368     CartVect n1 = a * b;
01369     CartVect n2 = a * s;
01370     CartVect n3 = s * b;
01371     if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
01372 
01373     // do the same for c, d, s, in plane z=0
01374     c[2] = d[2] = s[2] = 0.;  // bring everything in the same plane, z=0;
01375 
01376     n1 = c * d;
01377     n2 = c * s;
01378     n3 = s * d;
01379     if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
01380 
01381     return true;
01382 }
01383 
01384 ErrorCode IntxUtils::intersect_great_circle_arc_with_clat_arc( double* A, double* B, double* C, double* D, double R,
01385                                                                double* E, int& np )
01386 {
01387     const double distTol   = R * 1.e-6;
01388     const double Tolerance = R * R * 1.e-12;  // radius should be 1, usually
01389     np                     = 0;               // number of points in intersection
01390     CartVect a( A ), b( B ), c( C ), d( D );
01391     // check input first
01392     double R2 = R * R;
01393     if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
01394             fabs( d.length_squared() - R2 ) >
01395         10 * Tolerance )
01396         return MB_FAILURE;
01397 
01398     if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
01399     if( ( c - d ).length_squared() < Tolerance )  // edges are too short
01400         return MB_FAILURE;
01401 
01402     // CD is the const latitude arc
01403     if( fabs( C[2] - D[2] ) > distTol )  // cd is not on the same z (constant latitude)
01404         return MB_FAILURE;
01405 
01406     if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE;  // too close to the poles
01407 
01408     // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
01409     // great circle arc AB normal to the AB circle:
01410     CartVect n1 = a * b;  // the normal to the great circle arc (circle)
01411     // solve the system of equations:
01412     /*
01413      *    n1%(x, y, z) = 0  // on the great circle
01414      *     z = C[2];
01415      *    x^2+y^2+z^2 = R^2
01416      */
01417     double z = C[2];
01418     if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
01419     {
01420         // it is the Equator; check if the const lat edge is Equator too
01421         if( fabs( C[2] ) > distTol )
01422         {
01423             return MB_FAILURE;  // no intx, too far from Eq
01424         }
01425         else
01426         {
01427             // all points are on the equator
01428             //
01429             CartVect cd = c * d;
01430             // by convention,  c<d, positive is from c to d
01431             // is a or b between c , d?
01432             CartVect ca = c * a;
01433             CartVect ad = a * d;
01434             CartVect cb = c * b;
01435             CartVect bd = b * d;
01436             bool agtc   = ( ca % cd >= -Tolerance );  // a>c?
01437             bool dgta   = ( ad % cd >= -Tolerance );  // d>a?
01438             bool bgtc   = ( cb % cd >= -Tolerance );  // b>c?
01439             bool dgtb   = ( bd % cd >= -Tolerance );  // d>b?
01440             if( agtc )
01441             {
01442                 if( dgta )
01443                 {
01444                     // a is for sure a point
01445                     E[0] = a[0];
01446                     E[1] = a[1];
01447                     E[2] = a[2];
01448                     np++;
01449                     if( bgtc )
01450                     {
01451                         if( dgtb )
01452                         {
01453                             // b is also in between c and d
01454                             E[3] = b[0];
01455                             E[4] = b[1];
01456                             E[5] = b[2];
01457                             np++;
01458                         }
01459                         else
01460                         {
01461                             // then order is c a d b, intx is ad
01462                             E[3] = d[0];
01463                             E[4] = d[1];
01464                             E[5] = d[2];
01465                             np++;
01466                         }
01467                     }
01468                     else
01469                     {
01470                         // b is less than c, so b c a d, intx is ac
01471                         E[3] = c[0];
01472                         E[4] = c[1];
01473                         E[5] = c[2];
01474                         np++;  // what if E[0] is E[3]?
01475                     }
01476                 }
01477                 else  // c < d < a
01478                 {
01479                     if( dgtb )  // d is for sure in
01480                     {
01481                         E[0] = d[0];
01482                         E[1] = d[1];
01483                         E[2] = d[2];
01484                         np++;
01485                         if( bgtc )  // c<b<d<a
01486                         {
01487                             // another point is b
01488                             E[3] = b[0];
01489                             E[4] = b[1];
01490                             E[5] = b[2];
01491                             np++;
01492                         }
01493                         else  // b<c<d<a
01494                         {
01495                             // another point is c
01496                             E[3] = c[0];
01497                             E[4] = c[1];
01498                             E[5] = c[2];
01499                             np++;
01500                         }
01501                     }
01502                     else
01503                     {
01504                         // nothing, order is c, d < a, b
01505                     }
01506                 }
01507             }
01508             else  // a < c < d
01509             {
01510                 if( bgtc )
01511                 {
01512                     // c is for sure in
01513                     E[0] = c[0];
01514                     E[1] = c[1];
01515                     E[2] = c[2];
01516                     np++;
01517                     if( dgtb )
01518                     {
01519                         // a < c < b < d; second point is b
01520                         E[3] = b[0];
01521                         E[4] = b[1];
01522                         E[5] = b[2];
01523                         np++;
01524                     }
01525                     else
01526                     {
01527                         // a < c < d < b; second point is d
01528                         E[3] = d[0];
01529                         E[4] = d[1];
01530                         E[5] = d[2];
01531                         np++;
01532                     }
01533                 }
01534                 else  // a, b < c < d
01535                 {
01536                     // nothing
01537                 }
01538             }
01539         }
01540         // for the 2 points selected, see if it is only one?
01541         // no problem, maybe it will be collapsed later anyway
01542         if( np > 0 ) return MB_SUCCESS;
01543         return MB_FAILURE;  // no intersection
01544     }
01545     {
01546         if( fabs( n1[0] ) <= fabs( n1[1] ) )
01547         {
01548             // resolve eq in x:  n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
01549             //  (u+v*x)^2+x^2=R2-z^2
01550             //  (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
01551             //  delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
01552             // x1,2 =
01553             double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
01554             double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
01555             double delta = b1 * b1 - 4 * a1 * c1;
01556             if( delta < -Tolerance ) return MB_FAILURE;  // no intersection
01557             if( delta > Tolerance )                      // 2 solutions possible
01558             {
01559                 double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
01560                 double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
01561                 double y1 = u + v * x1;
01562                 double y2 = u + v * x2;
01563                 if( verify( a, b, c, d, x1, y1, z ) )
01564                 {
01565                     E[0] = x1;
01566                     E[1] = y1;
01567                     E[2] = z;
01568                     np++;
01569                 }
01570                 if( verify( a, b, c, d, x2, y2, z ) )
01571                 {
01572                     E[3 * np + 0] = x2;
01573                     E[3 * np + 1] = y2;
01574                     E[3 * np + 2] = z;
01575                     np++;
01576                 }
01577             }
01578             else
01579             {
01580                 // one solution
01581                 double x1 = -b1 / 2 / a1;
01582                 double y1 = u + v * x1;
01583                 if( verify( a, b, c, d, x1, y1, z ) )
01584                 {
01585                     E[0] = x1;
01586                     E[1] = y1;
01587                     E[2] = z;
01588                     np++;
01589                 }
01590             }
01591         }
01592         else
01593         {
01594             // resolve eq in y, reverse
01595             //   n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
01596             //  (u+v*y)^2+y^2 -R2+z^2 =0
01597             //  (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
01598             //
01599             // x1,2 =
01600             double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
01601             double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
01602             double delta = b1 * b1 - 4 * a1 * c1;
01603             if( delta < -Tolerance ) return MB_FAILURE;  // no intersection
01604             if( delta > Tolerance )                      // 2 solutions possible
01605             {
01606                 double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
01607                 double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
01608                 double x1 = u + v * y1;
01609                 double x2 = u + v * y2;
01610                 if( verify( a, b, c, d, x1, y1, z ) )
01611                 {
01612                     E[0] = x1;
01613                     E[1] = y1;
01614                     E[2] = z;
01615                     np++;
01616                 }
01617                 if( verify( a, b, c, d, x2, y2, z ) )
01618                 {
01619                     E[3 * np + 0] = x2;
01620                     E[3 * np + 1] = y2;
01621                     E[3 * np + 2] = z;
01622                     np++;
01623                 }
01624             }
01625             else
01626             {
01627                 // one solution
01628                 double y1 = -b1 / 2 / a1;
01629                 double x1 = u + v * y1;
01630                 if( verify( a, b, c, d, x1, y1, z ) )
01631                 {
01632                     E[0] = x1;
01633                     E[1] = y1;
01634                     E[2] = z;
01635                     np++;
01636                 }
01637             }
01638         }
01639     }
01640 
01641     if( np <= 0 ) return MB_FAILURE;
01642     return MB_SUCCESS;
01643 }
01644 
01645 #if 0
01646 ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1)
01647 {
01648   Range cells;
01649   ErrorCode rval = mb->get_entities_by_dimension(sf1, 2, cells);
01650   if (MB_SUCCESS!= rval)
01651   return rval;
01652   Range edges;
01653   rval = mb->get_adjacencies(cells, 1, true, edges, Interface::UNION);
01654   if (MB_SUCCESS!= rval)
01655   return rval;
01656 
01657   Tag edgeTypeTag;
01658   int default_int=0;
01659   rval = mb->tag_get_handle("edge_type", 1, MB_TYPE_INTEGER, edgeTypeTag,
01660       MB_TAG_DENSE | MB_TAG_CREAT, &default_int);
01661   if (MB_SUCCESS!= rval)
01662   return rval;
01663   // add edges to the set? not yet, maybe later
01664   // if edge horizontal, set value to 1
01665   int type_constant_lat=1;
01666   for (Range::iterator eit=edges.begin(); eit!=edges.end(); ++eit)
01667   {
01668     EntityHandle edge = *eit;
01669     const EntityHandle *conn=0;
01670     int num_n=0;
01671     rval = mb->get_connectivity(edge, conn, num_n );
01672     if (MB_SUCCESS!= rval)
01673     return rval;
01674     double coords[6];
01675     rval = mb->get_coords(conn, 2, coords);
01676     if (MB_SUCCESS!= rval)
01677     return rval;
01678     if (fabs( coords[2]-coords[5] )< 1.e-6 )
01679     {
01680       rval = mb->tag_set_data(edgeTypeTag, &edge, 1, &type_constant_lat);
01681       if (MB_SUCCESS!= rval)
01682       return rval;
01683     }
01684   }
01685 
01686   return MB_SUCCESS;
01687 }
01688 #endif
01689 
01690 // decide in a different metric if the corners of CS quad are
01691 // in the interior of an RLL quad
01692 int IntxUtils::borderPointsOfCSinRLL( CartVect* redc, double* red2dc, int nsRed, CartVect* bluec, int nsBlue,
01693                                       int* blueEdgeType, double* P, int* side, double epsil )
01694 {
01695     int extraPoints = 0;
01696     // first decide the blue z coordinates
01697     CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
01698     for( int i = 0; i < nsBlue; i++ )
01699     {
01700         if( blueEdgeType[i] == 0 )
01701         {
01702             int iP1 = ( i + 1 ) % nsBlue;
01703             if( bluec[i][2] > bluec[iP1][2] )
01704             {
01705                 A = bluec[i];
01706                 B = bluec[iP1];
01707                 C = bluec[( i + 2 ) % nsBlue];
01708                 D = bluec[( i + 3 ) % nsBlue];  // it could be back to A, if triangle on top
01709                 break;
01710             }
01711         }
01712     }
01713     if( nsBlue == 3 && B[2] < 0 )
01714     {
01715         // select D to be C
01716         D = C;
01717         C = B;  // B is the south pole then
01718     }
01719     // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
01720     // AB is const longitude, BC and DA constant latitude
01721     // check now each of the red points if they are inside this rectangle
01722     for( int i = 0; i < nsRed; i++ )
01723     {
01724         CartVect& X = redc[i];
01725         if( X[2] > A[2] || X[2] < B[2] ) continue;  // it is above or below the rectangle
01726         // now decide if it is between the planes OAB and OCD
01727         if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
01728         {
01729             side[i] = 1;  //
01730             // it means point X is in the rectangle that we want , on the sphere
01731             // pass the coords 2d
01732             P[extraPoints * 2]     = red2dc[2 * i];
01733             P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
01734             extraPoints++;
01735         }
01736     }
01737     return extraPoints;
01738 }
01739 
01740 ErrorCode IntxUtils::deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set )
01741 {
01742     ReadUtilIface* read_iface;
01743     ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
01744     // create the handle tag for the corresponding element / vertex
01745 
01746     EntityHandle dum = 0;
01747     Tag corrTag      = 0;  // it will be created here
01748     rval             = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
01749 
01750     // give the same global id to new verts and cells created in the lagr(departure) mesh
01751     Tag gid = mb->globalId_tag();
01752 
01753     Range quads;
01754     rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
01755 
01756     Range connecVerts;
01757     rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
01758 
01759     std::map< EntityHandle, EntityHandle > newNodes;
01760 
01761     std::vector< double* > coords;
01762     EntityHandle start_vert, start_elem, *connect;
01763     int num_verts = connecVerts.size();
01764     rval          = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
01765     if( MB_SUCCESS != rval ) return rval;
01766     // fill it up
01767     int i = 0;
01768     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
01769     {
01770         EntityHandle oldV = *vit;
01771         CartVect posi;
01772         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
01773 
01774         int global_id;
01775         rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
01776         EntityHandle new_vert = start_vert + i;
01777         // Cppcheck warning (false positive): variable coords is assigned a value that is never used
01778         coords[0][i] = posi[0];
01779         coords[1][i] = posi[1];
01780         coords[2][i] = posi[2];
01781 
01782         newNodes[oldV] = new_vert;
01783         // set also the correspondent tag :)
01784         rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
01785 
01786         // also the other side
01787         // need to check if we really need this; the new vertex will never need the old vertex
01788         // we have the global id which is the same
01789         rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
01790         // set the global id on the corresponding vertex the same as the initial vertex
01791         rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
01792     }
01793     // now create new quads in order (in a sequence)
01794 
01795     rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
01796     if( MB_SUCCESS != rval ) return rval;
01797     int ie = 0;
01798     for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
01799     {
01800         EntityHandle q = *it;
01801         int nnodes;
01802         const EntityHandle* conn;
01803         rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
01804         int global_id;
01805         rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
01806 
01807         for( int ii = 0; ii < nnodes; ii++ )
01808         {
01809             EntityHandle v1      = conn[ii];
01810             connect[4 * ie + ii] = newNodes[v1];
01811         }
01812         EntityHandle newElement = start_elem + ie;
01813 
01814         // set the corresponding tag; not sure we need this one, from old to new
01815         rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
01816         rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
01817 
01818         // set the global id
01819         rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
01820 
01821         rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
01822     }
01823 
01824     rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
01825 
01826     return MB_SUCCESS;
01827 }
01828 
01829 ErrorCode IntxUtils::remove_duplicate_vertices( Interface* mb, EntityHandle file_set, double merge_tol,
01830                                                 std::vector< Tag >& tagList )
01831 {
01832     Range verts;
01833     ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
01834     rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
01835 
01836     MergeMesh mm( mb );
01837 
01838     // remove the vertices from the set, before merging
01839 
01840     rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
01841 
01842     // now correct vertices that are repeated in polygons
01843     Range cells;
01844     rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
01845 
01846     verts.clear();
01847     rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
01848 
01849     Range modifiedCells;  // will be deleted at the end; keep the gid
01850     Range newCells;
01851 
01852     for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
01853     {
01854         EntityHandle cell          = *cit;
01855         const EntityHandle* connec = NULL;
01856         int num_verts              = 0;
01857         rval                       = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
01858 
01859         std::vector< EntityHandle > newConnec;
01860         newConnec.push_back( connec[0] );  // at least one vertex
01861         int index    = 0;
01862         int new_size = 1;
01863         while( index < num_verts - 2 )
01864         {
01865             int next_index = ( index + 1 );
01866             if( connec[next_index] != newConnec[new_size - 1] )
01867             {
01868                 newConnec.push_back( connec[next_index] );
01869                 new_size++;
01870             }
01871             index++;
01872         }
01873         // add the last one only if different from previous and first node
01874         if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
01875         {
01876             newConnec.push_back( connec[num_verts - 1] );
01877             new_size++;
01878         }
01879         if( new_size < num_verts )
01880         {
01881             // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
01882             modifiedCells.insert( cell );
01883             // create a new cell with type triangle, quad or polygon
01884             EntityType type = MBTRI;
01885             if( new_size == 3 )
01886                 type = MBTRI;
01887             else if( new_size == 4 )
01888                 type = MBQUAD;
01889             else if( new_size > 4 )
01890                 type = MBPOLYGON;
01891 
01892             // create new cell
01893             EntityHandle newCell;
01894             rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
01895             // set the old id to the new element
01896             newCells.insert( newCell );
01897             double value;  // use the same value to reset the tags, even if the tags are int (like Global ID)
01898             for( size_t i = 0; i < tagList.size(); i++ )
01899             {
01900                 rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
01901                 rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
01902             }
01903         }
01904     }
01905 
01906     rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
01907     rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
01908     rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
01909     rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
01910 
01911     return MB_SUCCESS;
01912 }
01913 
01914 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines