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