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 #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     { std::cout << " error in input " << a << " radius: " << Radius << " error:" << err1 << "\n"; }
00773     CartVect normalOAB = a * b;
00774     CartVect normalOCB = c * b;
00775     return angle( normalOAB, normalOCB );
00776 }
00777 
00778 // could be bigger than M_PI;
00779 // angle at B could be bigger than M_PI, if the orientation is such that ABC points toward the
00780 // interior
00781 double IntxUtils::oriented_spherical_angle( double* A, double* B, double* C )
00782 {
00783     // assume the same radius, sphere at origin
00784     CartVect a( A ), b( B ), c( C );
00785     CartVect normalOAB = a * b;
00786     CartVect normalOCB = c * b;
00787     CartVect orient    = ( c - b ) * ( a - b );
00788     double ang         = angle( normalOAB, normalOCB );  // this is between 0 and M_PI
00789     if( ang != ang )
00790     {
00791         // signal of a nan
00792         std::cout << a << " " << b << " " << c << "\n";
00793         std::cout << ang << "\n";
00794     }
00795     if( orient % b < 0 ) return ( 2 * M_PI - ang );  // the other angle, supplement
00796 
00797     return ang;
00798 }
00799 
00800 double IntxAreaUtils::area_spherical_triangle( double* A, double* B, double* C, double Radius )
00801 {
00802     switch( m_eAreaMethod )
00803     {
00804         case Girard:
00805             return area_spherical_triangle_girard( A, B, C, Radius );
00806 #ifdef MOAB_HAVE_TEMPESTREMAP
00807         case GaussQuadrature:
00808             return area_spherical_triangle_GQ( A, B, C, Radius );
00809 #endif
00810         case lHuiller:
00811         default:
00812             return area_spherical_triangle_lHuiller( A, B, C, Radius );
00813     }
00814 }
00815 
00816 double IntxAreaUtils::area_spherical_polygon( double* A, int N, double Radius, int* sign )
00817 {
00818     switch( m_eAreaMethod )
00819     {
00820         case Girard:
00821             return area_spherical_polygon_girard( A, N, Radius );
00822 #ifdef MOAB_HAVE_TEMPESTREMAP
00823         case GaussQuadrature:
00824             return area_spherical_polygon_GQ( A, N, Radius );
00825 #endif
00826         case lHuiller:
00827         default:
00828             return area_spherical_polygon_lHuiller( A, N, Radius, sign );
00829     }
00830 }
00831 
00832 double IntxAreaUtils::area_spherical_triangle_girard( double* A, double* B, double* C, double Radius )
00833 {
00834     double correction = spherical_angle( A, B, C, Radius ) + spherical_angle( B, C, A, Radius ) +
00835                         spherical_angle( C, A, B, Radius ) - M_PI;
00836     double area = Radius * Radius * correction;
00837     // now, is it negative or positive? is it pointing toward the center or outward?
00838     CartVect a( A ), b( B ), c( C );
00839     CartVect abc = ( b - a ) * ( c - a );
00840     if( abc % a > 0 )  // dot product positive, means ABC points out
00841         return area;
00842     else
00843         return -area;
00844 }
00845 
00846 double IntxAreaUtils::area_spherical_polygon_girard( double* A, int N, double Radius )
00847 {
00848     // this should work for non-convex polygons too
00849     // assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00850     //
00851     if( N <= 2 ) return 0.;
00852     double sum_angles = 0.;
00853     for( int i = 0; i < N; i++ )
00854     {
00855         int i1 = ( i + 1 ) % N;
00856         int i2 = ( i + 2 ) % N;
00857         sum_angles += IntxUtils::oriented_spherical_angle( A + 3 * i, A + 3 * i1, A + 3 * i2 );
00858     }
00859     double correction = sum_angles - ( N - 2 ) * M_PI;
00860     return Radius * Radius * correction;
00861 }
00862 
00863 double IntxAreaUtils::area_spherical_polygon_lHuiller( double* A, int N, double Radius, int* sign )
00864 {
00865     // This should work for non-convex polygons too
00866     // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00867     // We also assume that the orientation is positive;
00868     // If negative orientation, the area will be negative
00869     if( N <= 2 ) return 0.;
00870 
00871     int lsign   = 1;  // assume positive orientain
00872     double area = 0.;
00873     for( int i = 1; i < N - 1; i++ )
00874     {
00875         int i1              = i + 1;
00876         double areaTriangle = area_spherical_triangle_lHuiller( A, A + 3 * i, A + 3 * i1, Radius );
00877         if( areaTriangle < 0 )
00878             lsign = -1;  // signal that we have at least one triangle with negative orientation ;
00879                          // possible nonconvex polygon
00880         area += areaTriangle;
00881     }
00882     if( sign ) *sign = lsign;
00883 
00884     return area;
00885 }
00886 
00887 #ifdef MOAB_HAVE_TEMPESTREMAP
00888 double IntxAreaUtils::area_spherical_polygon_GQ( double* A, int N, double Radius )
00889 {
00890     // this should work for non-convex polygons too
00891     // In the input vector A, assume that the A, A+3, ..., A+3*(N-1) are the coordinates
00892     // We also assume that the orientation is positive;
00893     // If negative orientation, the area can be negative
00894     if( N <= 2 ) return 0.;
00895 
00896     // assume positive orientain
00897     double area = 0.;
00898     for( int i = 1; i < N - 1; i++ )
00899     {
00900         int i1 = i + 1;
00901         area += area_spherical_triangle_GQ( A, A + 3 * i, A + 3 * i1, Radius );
00902     }
00903     return area;
00904 }
00905 
00906 /* compute the area by using Gauss-Quadratures; use TR interfaces directly */
00907 double IntxAreaUtils::area_spherical_triangle_GQ( double* ptA, double* ptB, double* ptC, double )
00908 {
00909     Face face( 3 );
00910     NodeVector nodes( 3 );
00911     nodes[0] = Node( ptA[0], ptA[1], ptA[2] );
00912     nodes[1] = Node( ptB[0], ptB[1], ptB[2] );
00913     nodes[2] = Node( ptC[0], ptC[1], ptC[2] );
00914     face.SetNode( 0, 0 );
00915     face.SetNode( 1, 1 );
00916     face.SetNode( 2, 2 );
00917     return CalculateFaceArea( face, nodes );
00918 }
00919 #endif
00920 
00921 /*
00922  *  l'Huiller's formula for spherical triangle
00923  *  http://williams.best.vwh.net/avform.htm
00924  *  a, b, c are arc measures in radians, too
00925  *  A, B, C are angles on the sphere, for which we already have formula
00926  *               c
00927  *         A -------B
00928  *          \       |
00929  *           \      |
00930  *            \b    |a
00931  *             \    |
00932  *              \   |
00933  *               \  |
00934  *                \C|
00935  *                 \|
00936  *
00937  *  (The angle at B is not necessarily a right angle)
00938  *
00939  *    sin(a)  sin(b)   sin(c)
00940  *    ----- = ------ = ------
00941  *    sin(A)  sin(B)   sin(C)
00942  *
00943  * In terms of the sides (this is excess, as before, but numerically stable)
00944  *
00945  *  E = 4*atan(sqrt(tan(s/2)*tan((s-a)/2)*tan((s-b)/2)*tan((s-c)/2)))
00946  */
00947 double IntxAreaUtils::area_spherical_triangle_lHuiller( double* ptA, double* ptB, double* ptC, double Radius )
00948 {
00949 
00950     // now, a is angle BOC, O is origin
00951     CartVect vA( ptA ), vB( ptB ), vC( ptC );
00952     double a = angle( vB, vC );
00953     double b = angle( vC, vA );
00954     double c = angle( vA, vB );
00955     int sign = 1;
00956     if( ( vA * vB ) % vC < 0 ) sign = -1;
00957     double s   = ( a + b + c ) / 2;
00958     double tmp = tan( s / 2 ) * tan( ( s - a ) / 2 ) * tan( ( s - b ) / 2 ) * tan( ( s - c ) / 2 );
00959     if( tmp < 0. ) tmp = 0.;
00960 
00961     double E = 4 * atan( sqrt( tmp ) );
00962     if( E != E ) std::cout << " NaN at spherical triangle area \n";
00963 
00964     double area = sign * E * Radius * Radius;
00965 
00966 #ifdef CHECKNEGATIVEAREA
00967     if( area < 0 )
00968     {
00969         std::cout << "negative area: " << area << "\n";
00970         std::cout << std::setprecision( 15 );
00971         std::cout << "vA: " << vA << "\n";
00972         std::cout << "vB: " << vB << "\n";
00973         std::cout << "vC: " << vC << "\n";
00974         std::cout << "sign: " << sign << "\n";
00975         std::cout << " a: " << a << "\n";
00976         std::cout << " b: " << b << "\n";
00977         std::cout << " c: " << c << "\n";
00978     }
00979 #endif
00980 
00981     return area;
00982 }
00983 #undef CHECKNEGATIVEAREA
00984 
00985 double IntxAreaUtils::area_on_sphere( Interface* mb, EntityHandle set, double R )
00986 {
00987     // Get all entities of dimension 2
00988     Range inputRange;
00989     ErrorCode rval = mb->get_entities_by_dimension( set, 2, inputRange );MB_CHK_ERR_RET_VAL( rval, -1.0 );
00990 
00991     // Filter by elements that are owned by current process
00992     std::vector< int > ownerinfo( inputRange.size(), -1 );
00993     Tag intxOwnerTag;
00994     rval = mb->tag_get_handle( "ORIG_PROC", intxOwnerTag );
00995     if( MB_SUCCESS == rval )
00996     {
00997         rval = mb->tag_get_data( intxOwnerTag, inputRange, &ownerinfo[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
00998     }
00999 
01000     // compare total area with 4*M_PI * R^2
01001     int ie            = 0;
01002     double total_area = 0.;
01003     for( Range::iterator eit = inputRange.begin(); eit != inputRange.end(); ++eit )
01004     {
01005 
01006         // All zero/positive owner data represents ghosted elems
01007         if( ownerinfo[ie++] >= 0 ) continue;
01008 
01009         EntityHandle eh        = *eit;
01010         const double elem_area = this->area_spherical_element( mb, eh, R );
01011 
01012         // check whether the area of the spherical element is positive.
01013         assert( elem_area > 0 );
01014 
01015         // sum up the contribution
01016         total_area += elem_area;
01017     }
01018 
01019     // return total mesh area
01020     return total_area;
01021 }
01022 
01023 double IntxAreaUtils::area_spherical_element( Interface* mb, EntityHandle elem, double R )
01024 {
01025     // get the nodes, then the coordinates
01026     const EntityHandle* verts;
01027     int nsides;
01028     ErrorCode rval = mb->get_connectivity( elem, verts, nsides );MB_CHK_ERR_RET_VAL( rval, -1.0 );
01029 
01030     // account for possible padded polygons
01031     while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
01032         nsides--;
01033 
01034     // get coordinates
01035     std::vector< double > coords( 3 * nsides );
01036     rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR_RET_VAL( rval, -1.0 );
01037 
01038     // compute and return the area of the polygonal element
01039     return area_spherical_polygon( &coords[0], nsides, R );
01040 }
01041 
01042 double IntxUtils::distance_on_great_circle( CartVect& p1, CartVect& p2 )
01043 {
01044     SphereCoords sph1 = cart_to_spherical( p1 );
01045     SphereCoords sph2 = cart_to_spherical( p2 );
01046     // radius should be the same
01047     return sph1.R *
01048            acos( sin( sph1.lon ) * sin( sph2.lon ) + cos( sph1.lat ) * cos( sph2.lat ) * cos( sph2.lon - sph2.lon ) );
01049 }
01050 
01051 // break the nonconvex quads into triangles; remove the quad from the set? yes.
01052 // maybe radius is not needed;
01053 //
01054 ErrorCode IntxUtils::enforce_convexity( Interface* mb, EntityHandle lset, int my_rank )
01055 {
01056     // look at each polygon; compute all angles; if one is reflex, break that angle with
01057     // the next triangle; put the 2 new polys in the set;
01058     // still look at the next poly
01059     // replace it with 2 triangles, and remove from set;
01060     // it should work for all polygons / tested first for case 1, with dt 0.5 (too much deformation)
01061     // get all entities of dimension 2
01062     // then get the connectivity, etc
01063 
01064     Range inputRange;
01065     ErrorCode rval = mb->get_entities_by_dimension( lset, 2, inputRange );MB_CHK_ERR( rval );
01066 
01067     Tag corrTag       = 0;
01068     EntityHandle dumH = 0;
01069     rval              = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE, &dumH );
01070     if( rval == MB_TAG_NOT_FOUND ) corrTag = 0;
01071 
01072     Tag gidTag;
01073     rval = mb->tag_get_handle( "GLOBAL_ID", 1, MB_TYPE_INTEGER, gidTag, MB_TAG_DENSE );MB_CHK_ERR( rval );
01074 
01075     std::vector< double > coords;
01076     coords.resize( 3 * MAXEDGES );  // at most 10 vertices per polygon
01077     // we should create a queue with new polygons that need processing for reflex angles
01078     //  (obtuse)
01079     std::queue< EntityHandle > newPolys;
01080     int brokenPolys     = 0;
01081     Range::iterator eit = inputRange.begin();
01082     while( eit != inputRange.end() || !newPolys.empty() )
01083     {
01084         EntityHandle eh;
01085         if( eit != inputRange.end() )
01086         {
01087             eh = *eit;
01088             ++eit;
01089         }
01090         else
01091         {
01092             eh = newPolys.front();
01093             newPolys.pop();
01094         }
01095         // get the nodes, then the coordinates
01096         const EntityHandle* verts;
01097         int num_nodes;
01098         rval = mb->get_connectivity( eh, verts, num_nodes );MB_CHK_ERR( rval );
01099         int nsides = num_nodes;
01100         // account for possible padded polygons
01101         while( verts[nsides - 2] == verts[nsides - 1] && nsides > 3 )
01102             nsides--;
01103         EntityHandle corrHandle = 0;
01104         if( corrTag )
01105         {
01106             rval = mb->tag_get_data( corrTag, &eh, 1, &corrHandle );MB_CHK_ERR( rval );
01107         }
01108         int gid = 0;
01109         rval    = mb->tag_get_data( gidTag, &eh, 1, &gid );MB_CHK_ERR( rval );
01110         coords.resize( 3 * nsides );
01111         if( nsides < 4 ) continue;  // if already triangles, don't bother
01112         // get coordinates
01113         rval = mb->get_coords( verts, nsides, &coords[0] );MB_CHK_ERR( rval );
01114         // compute each angle
01115         bool alreadyBroken = false;
01116 
01117         for( int i = 0; i < nsides; i++ )
01118         {
01119             double* A    = &coords[3 * i];
01120             double* B    = &coords[3 * ( ( i + 1 ) % nsides )];
01121             double* C    = &coords[3 * ( ( i + 2 ) % nsides )];
01122             double angle = IntxUtils::oriented_spherical_angle( A, B, C );
01123             if( angle - M_PI > 0. )  // even almost reflex is bad; break it!
01124             {
01125                 if( alreadyBroken )
01126                 {
01127                     mb->list_entities( &eh, 1 );
01128                     mb->list_entities( verts, nsides );
01129                     double* D = &coords[3 * ( ( i + 3 ) % nsides )];
01130                     std::cout << "ABC: " << angle << " \n";
01131                     std::cout << "BCD: " << IntxUtils::oriented_spherical_angle( B, C, D ) << " \n";
01132                     std::cout << "CDA: " << IntxUtils::oriented_spherical_angle( C, D, A ) << " \n";
01133                     std::cout << "DAB: " << IntxUtils::oriented_spherical_angle( D, A, B ) << " \n";
01134                     std::cout << " this cell has at least 2 angles > 180, it has serious issues\n";
01135 
01136                     return MB_FAILURE;
01137                 }
01138                 // the bad angle is at i+1;
01139                 // create 1 triangle and one polygon; add the polygon to the input range, so
01140                 // it will be processed too
01141                 // also, add both to the set :) and remove the original polygon from the set
01142                 // break the next triangle, even though not optimal
01143                 // so create the triangle i+1, i+2, i+3; remove i+2 from original list
01144                 // even though not optimal in general, it is good enough.
01145                 EntityHandle conn3[3] = { verts[( i + 1 ) % nsides], verts[( i + 2 ) % nsides],
01146                                           verts[( i + 3 ) % nsides] };
01147                 // create a polygon with num_nodes-1 vertices, and connectivity
01148                 // verts[i+1], verts[i+3], (all except i+2)
01149                 std::vector< EntityHandle > conn( nsides - 1 );
01150                 for( int j = 1; j < nsides; j++ )
01151                 {
01152                     conn[j - 1] = verts[( i + j + 2 ) % nsides];
01153                 }
01154                 EntityHandle newElement;
01155                 rval = mb->create_element( MBTRI, conn3, 3, newElement );MB_CHK_ERR( rval );
01156 
01157                 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
01158                 if( corrTag )
01159                 {
01160                     rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
01161                 }
01162                 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
01163                 if( nsides == 4 )
01164                 {
01165                     // create another triangle
01166                     rval = mb->create_element( MBTRI, &conn[0], 3, newElement );MB_CHK_ERR( rval );
01167                 }
01168                 else
01169                 {
01170                     // create another polygon, and add it to the inputRange
01171                     rval = mb->create_element( MBPOLYGON, &conn[0], nsides - 1, newElement );MB_CHK_ERR( rval );
01172                     newPolys.push( newElement );  // because it has less number of edges, the
01173                     // reverse should work to find it.
01174                 }
01175                 rval = mb->add_entities( lset, &newElement, 1 );MB_CHK_ERR( rval );
01176                 if( corrTag )
01177                 {
01178                     rval = mb->tag_set_data( corrTag, &newElement, 1, &corrHandle );MB_CHK_ERR( rval );
01179                 }
01180                 rval = mb->tag_set_data( gidTag, &newElement, 1, &gid );MB_CHK_ERR( rval );
01181                 rval = mb->remove_entities( lset, &eh, 1 );MB_CHK_ERR( rval );
01182                 brokenPolys++;
01183                 alreadyBroken = true;  // get out of the loop, element is broken
01184             }
01185         }
01186     }
01187     if (brokenPolys > 0)
01188     {
01189         std::cout << "on local process " << my_rank << ", " << brokenPolys
01190                   << " concave polygons were decomposed in convex ones \n";
01191 #ifdef VERBOSE
01192         std::stringstream fff;
01193         fff << "file_set" <<  mb->id_from_handle(lset) << "rk_"<< my_rank << ".h5m";
01194         rval = mb->write_file(fff.str().c_str(), 0, 0, &lset, 1);MB_CHK_ERR( rval );
01195         std::cout << "wrote new file set: " << fff.str() << "\n";
01196 #endif
01197     }
01198     return MB_SUCCESS;
01199 }
01200 
01201 // looking at quad connectivity, collapse to triangle if 2 nodes equal
01202 // then delete the old quad
01203 ErrorCode IntxUtils::fix_degenerate_quads( Interface* mb, EntityHandle set )
01204 {
01205     Range quads;
01206     ErrorCode rval = mb->get_entities_by_type( set, MBQUAD, quads );MB_CHK_ERR( rval );
01207     Tag gid;
01208     gid = mb->globalId_tag();
01209     for( Range::iterator qit = quads.begin(); qit != quads.end(); ++qit )
01210     {
01211         EntityHandle quad         = *qit;
01212         const EntityHandle* conn4 = NULL;
01213         int num_nodes             = 0;
01214         rval                      = mb->get_connectivity( quad, conn4, num_nodes );MB_CHK_ERR( rval );
01215         for( int i = 0; i < num_nodes; i++ )
01216         {
01217             int next_node_index = ( i + 1 ) % num_nodes;
01218             if( conn4[i] == conn4[next_node_index] )
01219             {
01220                 // form a triangle and delete the quad
01221                 // first get the global id, to set it on triangle later
01222                 int global_id = 0;
01223                 rval          = mb->tag_get_data( gid, &quad, 1, &global_id );MB_CHK_ERR( rval );
01224                 int i2                = ( i + 2 ) % num_nodes;
01225                 int i3                = ( i + 3 ) % num_nodes;
01226                 EntityHandle conn3[3] = { conn4[i], conn4[i2], conn4[i3] };
01227                 EntityHandle tri;
01228                 rval = mb->create_element( MBTRI, conn3, 3, tri );MB_CHK_ERR( rval );
01229                 mb->add_entities( set, &tri, 1 );
01230                 mb->remove_entities( set, &quad, 1 );
01231                 mb->delete_entities( &quad, 1 );
01232                 rval = mb->tag_set_data( gid, &tri, 1, &global_id );MB_CHK_ERR( rval );
01233             }
01234         }
01235     }
01236     return MB_SUCCESS;
01237 }
01238 
01239 ErrorCode IntxAreaUtils::positive_orientation( Interface* mb, EntityHandle set, double R )
01240 {
01241     Range cells2d;
01242     ErrorCode rval = mb->get_entities_by_dimension( set, 2, cells2d );
01243     if( MB_SUCCESS != rval ) return rval;
01244     for( Range::iterator qit = cells2d.begin(); qit != cells2d.end(); ++qit )
01245     {
01246         EntityHandle cell        = *qit;
01247         const EntityHandle* conn = NULL;
01248         int num_nodes            = 0;
01249         rval                     = mb->get_connectivity( cell, conn, num_nodes );
01250         if( MB_SUCCESS != rval ) return rval;
01251         if( num_nodes < 3 ) return MB_FAILURE;
01252 
01253         double coords[9];
01254         rval = mb->get_coords( conn, 3, coords );
01255         if( MB_SUCCESS != rval ) return rval;
01256 
01257         double area;
01258         if( R > 0 )
01259             area = area_spherical_triangle_lHuiller( coords, coords + 3, coords + 6, R );
01260         else
01261             area = IntxUtils::area2D( coords, coords + 3, coords + 6 );
01262         if( area < 0 )
01263         {
01264             // compute all area, do not revert if total area is positive
01265             std::vector< double > coords2( 3 * num_nodes );
01266             // get coordinates
01267             rval = mb->get_coords( conn, num_nodes, &coords2[0] );
01268             if( MB_SUCCESS != rval ) return MB_FAILURE;
01269             double totArea = area_spherical_polygon_lHuiller( &coords2[0], num_nodes, R );
01270             if( totArea < 0 )
01271             {
01272                 std::vector< EntityHandle > newconn( num_nodes );
01273                 for( int i = 0; i < num_nodes; i++ )
01274                 {
01275                     newconn[num_nodes - 1 - i] = conn[i];
01276                 }
01277                 rval = mb->set_connectivity( cell, &newconn[0], num_nodes );
01278                 if( MB_SUCCESS != rval ) return rval;
01279             }
01280             else
01281             {
01282                 std::cout << " nonconvex problem first area:" << area << " total area: " << totArea << std::endl;
01283             }
01284         }
01285     }
01286     return MB_SUCCESS;
01287 }
01288 
01289 // distance along a great circle on a sphere of radius 1
01290 // page 4
01291 double IntxUtils::distance_on_sphere( double la1, double te1, double la2, double te2 )
01292 {
01293     return acos( sin( te1 ) * sin( te2 ) + cos( te1 ) * cos( te2 ) * cos( la1 - la2 ) );
01294 }
01295 
01296 /*
01297  * given 2 great circle arcs, AB and CD, compute the unique intersection point, if it exists
01298  *  in between
01299  */
01300 ErrorCode IntxUtils::intersect_great_circle_arcs( double* A, double* B, double* C, double* D, double R, double* E )
01301 {
01302     // first verify A, B, C, D are on the same sphere
01303     double R2              = R * R;
01304     const double Tolerance = 1.e-12 * R2;
01305 
01306     CartVect a( A ), b( B ), c( C ), d( D );
01307 
01308     if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
01309             fabs( d.length_squared() - R2 ) >
01310         10 * Tolerance )
01311         return MB_FAILURE;
01312 
01313     CartVect n1 = a * b;
01314     if( n1.length_squared() < Tolerance ) return MB_FAILURE;
01315 
01316     CartVect n2 = c * d;
01317     if( n2.length_squared() < Tolerance ) return MB_FAILURE;
01318     CartVect n3 = n1 * n2;
01319     n3.normalize();
01320 
01321     n3 = R * n3;
01322     // the intersection is either n3 or -n3
01323     CartVect n4 = a * n3, n5 = n3 * b;
01324     if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
01325     {
01326         // n3 is good for ab, see if it is good for cd
01327         n4 = c * n3;
01328         n5 = n3 * d;
01329         if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
01330         {
01331             E[0] = n3[0];
01332             E[1] = n3[1];
01333             E[2] = n3[2];
01334         }
01335         else
01336             return MB_FAILURE;
01337     }
01338     else
01339     {
01340         // try -n3
01341         n3 = -n3;
01342         n4 = a * n3, n5 = n3 * b;
01343         if( n1 % n4 >= -Tolerance && n1 % n5 >= -Tolerance )
01344         {
01345             // n3 is good for ab, see if it is good for cd
01346             n4 = c * n3;
01347             n5 = n3 * d;
01348             if( n2 % n4 >= -Tolerance && n2 % n5 >= -Tolerance )
01349             {
01350                 E[0] = n3[0];
01351                 E[1] = n3[1];
01352                 E[2] = n3[2];
01353             }
01354             else
01355                 return MB_FAILURE;
01356         }
01357         else
01358             return MB_FAILURE;
01359     }
01360 
01361     return MB_SUCCESS;
01362 }
01363 
01364 // verify that result is in between a and b on a great circle arc, and between c and d on a constant
01365 // latitude arc
01366 static bool verify( CartVect a, CartVect b, CartVect c, CartVect d, double x, double y, double z )
01367 {
01368     // to check, the point has to be between a and b on a great arc, and between c and d on a const
01369     // lat circle
01370     CartVect s( x, y, z );
01371     CartVect n1 = a * b;
01372     CartVect n2 = a * s;
01373     CartVect n3 = s * b;
01374     if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
01375 
01376     // do the same for c, d, s, in plane z=0
01377     c[2] = d[2] = s[2] = 0.;  // bring everything in the same plane, z=0;
01378 
01379     n1 = c * d;
01380     n2 = c * s;
01381     n3 = s * d;
01382     if( n1 % n2 < 0 || n1 % n3 < 0 ) return false;
01383 
01384     return true;
01385 }
01386 
01387 ErrorCode IntxUtils::intersect_great_circle_arc_with_clat_arc( double* A, double* B, double* C, double* D, double R,
01388                                                                double* E, int& np )
01389 {
01390     const double distTol   = R * 1.e-6;
01391     const double Tolerance = R * R * 1.e-12;  // radius should be 1, usually
01392     np                     = 0;               // number of points in intersection
01393     CartVect a( A ), b( B ), c( C ), d( D );
01394     // check input first
01395     double R2 = R * R;
01396     if( fabs( a.length_squared() - R2 ) + fabs( b.length_squared() - R2 ) + fabs( c.length_squared() - R2 ) +
01397             fabs( d.length_squared() - R2 ) >
01398         10 * Tolerance )
01399         return MB_FAILURE;
01400 
01401     if( ( a - b ).length_squared() < Tolerance ) return MB_FAILURE;
01402     if( ( c - d ).length_squared() < Tolerance )  // edges are too short
01403         return MB_FAILURE;
01404 
01405     // CD is the const latitude arc
01406     if( fabs( C[2] - D[2] ) > distTol )  // cd is not on the same z (constant latitude)
01407         return MB_FAILURE;
01408 
01409     if( fabs( R - C[2] ) < distTol || fabs( R + C[2] ) < distTol ) return MB_FAILURE;  // too close to the poles
01410 
01411     // find the points on the circle P(teta) = (r*sin(teta), r*cos(teta), C[2]) that are on the
01412     // great circle arc AB normal to the AB circle:
01413     CartVect n1 = a * b;  // the normal to the great circle arc (circle)
01414     // solve the system of equations:
01415     /*
01416      *    n1%(x, y, z) = 0  // on the great circle
01417      *     z = C[2];
01418      *    x^2+y^2+z^2 = R^2
01419      */
01420     double z = C[2];
01421     if( fabs( n1[0] ) + fabs( n1[1] ) < 2 * Tolerance )
01422     {
01423         // it is the Equator; check if the const lat edge is Equator too
01424         if( fabs( C[2] ) > distTol )
01425         {
01426             return MB_FAILURE;  // no intx, too far from Eq
01427         }
01428         else
01429         {
01430             // all points are on the equator
01431             //
01432             CartVect cd = c * d;
01433             // by convention,  c<d, positive is from c to d
01434             // is a or b between c , d?
01435             CartVect ca = c * a;
01436             CartVect ad = a * d;
01437             CartVect cb = c * b;
01438             CartVect bd = b * d;
01439             bool agtc   = ( ca % cd >= -Tolerance );  // a>c?
01440             bool dgta   = ( ad % cd >= -Tolerance );  // d>a?
01441             bool bgtc   = ( cb % cd >= -Tolerance );  // b>c?
01442             bool dgtb   = ( bd % cd >= -Tolerance );  // d>b?
01443             if( agtc )
01444             {
01445                 if( dgta )
01446                 {
01447                     // a is for sure a point
01448                     E[0] = a[0];
01449                     E[1] = a[1];
01450                     E[2] = a[2];
01451                     np++;
01452                     if( bgtc )
01453                     {
01454                         if( dgtb )
01455                         {
01456                             // b is also in between c and d
01457                             E[3] = b[0];
01458                             E[4] = b[1];
01459                             E[5] = b[2];
01460                             np++;
01461                         }
01462                         else
01463                         {
01464                             // then order is c a d b, intx is ad
01465                             E[3] = d[0];
01466                             E[4] = d[1];
01467                             E[5] = d[2];
01468                             np++;
01469                         }
01470                     }
01471                     else
01472                     {
01473                         // b is less than c, so b c a d, intx is ac
01474                         E[3] = c[0];
01475                         E[4] = c[1];
01476                         E[5] = c[2];
01477                         np++;  // what if E[0] is E[3]?
01478                     }
01479                 }
01480                 else  // c < d < a
01481                 {
01482                     if( dgtb )  // d is for sure in
01483                     {
01484                         E[0] = d[0];
01485                         E[1] = d[1];
01486                         E[2] = d[2];
01487                         np++;
01488                         if( bgtc )  // c<b<d<a
01489                         {
01490                             // another point is b
01491                             E[3] = b[0];
01492                             E[4] = b[1];
01493                             E[5] = b[2];
01494                             np++;
01495                         }
01496                         else  // b<c<d<a
01497                         {
01498                             // another point is c
01499                             E[3] = c[0];
01500                             E[4] = c[1];
01501                             E[5] = c[2];
01502                             np++;
01503                         }
01504                     }
01505                     else
01506                     {
01507                         // nothing, order is c, d < a, b
01508                     }
01509                 }
01510             }
01511             else  // a < c < d
01512             {
01513                 if( bgtc )
01514                 {
01515                     // c is for sure in
01516                     E[0] = c[0];
01517                     E[1] = c[1];
01518                     E[2] = c[2];
01519                     np++;
01520                     if( dgtb )
01521                     {
01522                         // a < c < b < d; second point is b
01523                         E[3] = b[0];
01524                         E[4] = b[1];
01525                         E[5] = b[2];
01526                         np++;
01527                     }
01528                     else
01529                     {
01530                         // a < c < d < b; second point is d
01531                         E[3] = d[0];
01532                         E[4] = d[1];
01533                         E[5] = d[2];
01534                         np++;
01535                     }
01536                 }
01537                 else  // a, b < c < d
01538                 {
01539                     // nothing
01540                 }
01541             }
01542         }
01543         // for the 2 points selected, see if it is only one?
01544         // no problem, maybe it will be collapsed later anyway
01545         if( np > 0 ) return MB_SUCCESS;
01546         return MB_FAILURE;  // no intersection
01547     }
01548     {
01549         if( fabs( n1[0] ) <= fabs( n1[1] ) )
01550         {
01551             // resolve eq in x:  n0 * x + n1 * y +n2*z = 0; y = -n2/n1*z -n0/n1*x
01552             //  (u+v*x)^2+x^2=R2-z^2
01553             //  (v^2+1)*x^2 + 2*u*v *x + u^2+z^2-R^2 = 0
01554             //  delta = 4*u^2*v^2 - 4*(v^2-1)(u^2+z^2-R^2)
01555             // x1,2 =
01556             double u = -n1[2] / n1[1] * z, v = -n1[0] / n1[1];
01557             double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
01558             double delta = b1 * b1 - 4 * a1 * c1;
01559             if( delta < -Tolerance ) return MB_FAILURE;  // no intersection
01560             if( delta > Tolerance )                      // 2 solutions possible
01561             {
01562                 double x1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
01563                 double x2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
01564                 double y1 = u + v * x1;
01565                 double y2 = u + v * x2;
01566                 if( verify( a, b, c, d, x1, y1, z ) )
01567                 {
01568                     E[0] = x1;
01569                     E[1] = y1;
01570                     E[2] = z;
01571                     np++;
01572                 }
01573                 if( verify( a, b, c, d, x2, y2, z ) )
01574                 {
01575                     E[3 * np + 0] = x2;
01576                     E[3 * np + 1] = y2;
01577                     E[3 * np + 2] = z;
01578                     np++;
01579                 }
01580             }
01581             else
01582             {
01583                 // one solution
01584                 double x1 = -b1 / 2 / a1;
01585                 double y1 = u + v * x1;
01586                 if( verify( a, b, c, d, x1, y1, z ) )
01587                 {
01588                     E[0] = x1;
01589                     E[1] = y1;
01590                     E[2] = z;
01591                     np++;
01592                 }
01593             }
01594         }
01595         else
01596         {
01597             // resolve eq in y, reverse
01598             //   n0 * x + n1 * y +n2*z = 0; x = -n2/n0*z -n1/n0*y = u+v*y
01599             //  (u+v*y)^2+y^2 -R2+z^2 =0
01600             //  (v^2+1)*y^2 + 2*u*v *y + u^2+z^2-R^2 = 0
01601             //
01602             // x1,2 =
01603             double u = -n1[2] / n1[0] * z, v = -n1[1] / n1[0];
01604             double a1 = v * v + 1, b1 = 2 * u * v, c1 = u * u + z * z - R2;
01605             double delta = b1 * b1 - 4 * a1 * c1;
01606             if( delta < -Tolerance ) return MB_FAILURE;  // no intersection
01607             if( delta > Tolerance )                      // 2 solutions possible
01608             {
01609                 double y1 = ( -b1 + sqrt( delta ) ) / 2 / a1;
01610                 double y2 = ( -b1 - sqrt( delta ) ) / 2 / a1;
01611                 double x1 = u + v * y1;
01612                 double x2 = u + v * y2;
01613                 if( verify( a, b, c, d, x1, y1, z ) )
01614                 {
01615                     E[0] = x1;
01616                     E[1] = y1;
01617                     E[2] = z;
01618                     np++;
01619                 }
01620                 if( verify( a, b, c, d, x2, y2, z ) )
01621                 {
01622                     E[3 * np + 0] = x2;
01623                     E[3 * np + 1] = y2;
01624                     E[3 * np + 2] = z;
01625                     np++;
01626                 }
01627             }
01628             else
01629             {
01630                 // one solution
01631                 double y1 = -b1 / 2 / a1;
01632                 double x1 = u + v * y1;
01633                 if( verify( a, b, c, d, x1, y1, z ) )
01634                 {
01635                     E[0] = x1;
01636                     E[1] = y1;
01637                     E[2] = z;
01638                     np++;
01639                 }
01640             }
01641         }
01642     }
01643 
01644     if( np <= 0 ) return MB_FAILURE;
01645     return MB_SUCCESS;
01646 }
01647 
01648 #if 0
01649 ErrorCode set_edge_type_flag(Interface * mb, EntityHandle sf1)
01650 {
01651   Range cells;
01652   ErrorCode rval = mb->get_entities_by_dimension(sf1, 2, cells);
01653   if (MB_SUCCESS!= rval)
01654   return rval;
01655   Range edges;
01656   rval = mb->get_adjacencies(cells, 1, true, edges, Interface::UNION);
01657   if (MB_SUCCESS!= rval)
01658   return rval;
01659 
01660   Tag edgeTypeTag;
01661   int default_int=0;
01662   rval = mb->tag_get_handle("edge_type", 1, MB_TYPE_INTEGER, edgeTypeTag,
01663       MB_TAG_DENSE | MB_TAG_CREAT, &default_int);
01664   if (MB_SUCCESS!= rval)
01665   return rval;
01666   // add edges to the set? not yet, maybe later
01667   // if edge horizontal, set value to 1
01668   int type_constant_lat=1;
01669   for (Range::iterator eit=edges.begin(); eit!=edges.end(); ++eit)
01670   {
01671     EntityHandle edge = *eit;
01672     const EntityHandle *conn=0;
01673     int num_n=0;
01674     rval = mb->get_connectivity(edge, conn, num_n );
01675     if (MB_SUCCESS!= rval)
01676     return rval;
01677     double coords[6];
01678     rval = mb->get_coords(conn, 2, coords);
01679     if (MB_SUCCESS!= rval)
01680     return rval;
01681     if (fabs( coords[2]-coords[5] )< 1.e-6 )
01682     {
01683       rval = mb->tag_set_data(edgeTypeTag, &edge, 1, &type_constant_lat);
01684       if (MB_SUCCESS!= rval)
01685       return rval;
01686     }
01687   }
01688 
01689   return MB_SUCCESS;
01690 }
01691 #endif
01692 
01693 // decide in a different metric if the corners of CS quad are
01694 // in the interior of an RLL quad
01695 int IntxUtils::borderPointsOfCSinRLL( CartVect* redc, double* red2dc, int nsRed, CartVect* bluec, int nsBlue,
01696                                       int* blueEdgeType, double* P, int* side, double epsil )
01697 {
01698     int extraPoints = 0;
01699     // first decide the blue z coordinates
01700     CartVect A( 0. ), B( 0. ), C( 0. ), D( 0. );
01701     for( int i = 0; i < nsBlue; i++ )
01702     {
01703         if( blueEdgeType[i] == 0 )
01704         {
01705             int iP1 = ( i + 1 ) % nsBlue;
01706             if( bluec[i][2] > bluec[iP1][2] )
01707             {
01708                 A = bluec[i];
01709                 B = bluec[iP1];
01710                 C = bluec[( i + 2 ) % nsBlue];
01711                 D = bluec[( i + 3 ) % nsBlue];  // it could be back to A, if triangle on top
01712                 break;
01713             }
01714         }
01715     }
01716     if( nsBlue == 3 && B[2] < 0 )
01717     {
01718         // select D to be C
01719         D = C;
01720         C = B;  // B is the south pole then
01721     }
01722     // so we have A, B, C, D, with A at the top, b going down, then C, D, D > C, A > B
01723     // AB is const longitude, BC and DA constant latitude
01724     // check now each of the red points if they are inside this rectangle
01725     for( int i = 0; i < nsRed; i++ )
01726     {
01727         CartVect& X = redc[i];
01728         if( X[2] > A[2] || X[2] < B[2] ) continue;  // it is above or below the rectangle
01729         // now decide if it is between the planes OAB and OCD
01730         if( ( ( A * B ) % X >= -epsil ) && ( ( C * D ) % X >= -epsil ) )
01731         {
01732             side[i] = 1;  //
01733             // it means point X is in the rectangle that we want , on the sphere
01734             // pass the coords 2d
01735             P[extraPoints * 2]     = red2dc[2 * i];
01736             P[extraPoints * 2 + 1] = red2dc[2 * i + 1];
01737             extraPoints++;
01738         }
01739     }
01740     return extraPoints;
01741 }
01742 
01743 ErrorCode IntxUtils::deep_copy_set_with_quads( Interface* mb, EntityHandle source_set, EntityHandle dest_set )
01744 {
01745     ReadUtilIface* read_iface;
01746     ErrorCode rval = mb->query_interface( read_iface );MB_CHK_ERR( rval );
01747     // create the handle tag for the corresponding element / vertex
01748 
01749     EntityHandle dum = 0;
01750     Tag corrTag      = 0;  // it will be created here
01751     rval             = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );MB_CHK_ERR( rval );
01752 
01753     // give the same global id to new verts and cells created in the lagr(departure) mesh
01754     Tag gid = mb->globalId_tag();
01755 
01756     Range quads;
01757     rval = mb->get_entities_by_type( source_set, MBQUAD, quads );MB_CHK_ERR( rval );
01758 
01759     Range connecVerts;
01760     rval = mb->get_connectivity( quads, connecVerts );MB_CHK_ERR( rval );
01761 
01762     std::map< EntityHandle, EntityHandle > newNodes;
01763 
01764     std::vector< double* > coords;
01765     EntityHandle start_vert, start_elem, *connect;
01766     int num_verts = connecVerts.size();
01767     rval          = read_iface->get_node_coords( 3, num_verts, 0, start_vert, coords );
01768     if( MB_SUCCESS != rval ) return rval;
01769     // fill it up
01770     int i = 0;
01771     for( Range::iterator vit = connecVerts.begin(); vit != connecVerts.end(); ++vit, i++ )
01772     {
01773         EntityHandle oldV = *vit;
01774         CartVect posi;
01775         rval = mb->get_coords( &oldV, 1, &( posi[0] ) );MB_CHK_ERR( rval );
01776 
01777         int global_id;
01778         rval = mb->tag_get_data( gid, &oldV, 1, &global_id );MB_CHK_ERR( rval );
01779         EntityHandle new_vert = start_vert + i;
01780         // Cppcheck warning (false positive): variable coords is assigned a value that is never used
01781         coords[0][i] = posi[0];
01782         coords[1][i] = posi[1];
01783         coords[2][i] = posi[2];
01784 
01785         newNodes[oldV] = new_vert;
01786         // set also the correspondent tag :)
01787         rval = mb->tag_set_data( corrTag, &oldV, 1, &new_vert );MB_CHK_ERR( rval );
01788 
01789         // also the other side
01790         // need to check if we really need this; the new vertex will never need the old vertex
01791         // we have the global id which is the same
01792         rval = mb->tag_set_data( corrTag, &new_vert, 1, &oldV );MB_CHK_ERR( rval );
01793         // set the global id on the corresponding vertex the same as the initial vertex
01794         rval = mb->tag_set_data( gid, &new_vert, 1, &global_id );MB_CHK_ERR( rval );
01795     }
01796     // now create new quads in order (in a sequence)
01797 
01798     rval = read_iface->get_element_connect( quads.size(), 4, MBQUAD, 0, start_elem, connect );
01799     if( MB_SUCCESS != rval ) return rval;
01800     int ie = 0;
01801     for( Range::iterator it = quads.begin(); it != quads.end(); ++it, ie++ )
01802     {
01803         EntityHandle q = *it;
01804         int nnodes;
01805         const EntityHandle* conn;
01806         rval = mb->get_connectivity( q, conn, nnodes );MB_CHK_ERR( rval );
01807         int global_id;
01808         rval = mb->tag_get_data( gid, &q, 1, &global_id );MB_CHK_ERR( rval );
01809 
01810         for( int ii = 0; ii < nnodes; ii++ )
01811         {
01812             EntityHandle v1      = conn[ii];
01813             connect[4 * ie + ii] = newNodes[v1];
01814         }
01815         EntityHandle newElement = start_elem + ie;
01816 
01817         // set the corresponding tag; not sure we need this one, from old to new
01818         rval = mb->tag_set_data( corrTag, &q, 1, &newElement );MB_CHK_ERR( rval );
01819         rval = mb->tag_set_data( corrTag, &newElement, 1, &q );MB_CHK_ERR( rval );
01820 
01821         // set the global id
01822         rval = mb->tag_set_data( gid, &newElement, 1, &global_id );MB_CHK_ERR( rval );
01823 
01824         rval = mb->add_entities( dest_set, &newElement, 1 );MB_CHK_ERR( rval );
01825     }
01826 
01827     rval = read_iface->update_adjacencies( start_elem, quads.size(), 4, connect );MB_CHK_ERR( rval );
01828 
01829     return MB_SUCCESS;
01830 }
01831 
01832 ErrorCode IntxUtils::remove_duplicate_vertices( Interface* mb, EntityHandle file_set, double merge_tol,
01833                                                 std::vector< Tag >& tagList )
01834 {
01835     Range verts;
01836     ErrorCode rval = mb->get_entities_by_dimension( file_set, 0, verts );MB_CHK_ERR( rval );
01837     rval = mb->remove_entities( file_set, verts );MB_CHK_ERR( rval );
01838 
01839     MergeMesh mm( mb );
01840 
01841     // remove the vertices from the set, before merging
01842 
01843     rval = mm.merge_all( file_set, merge_tol );MB_CHK_ERR( rval );
01844 
01845     // now correct vertices that are repeated in polygons
01846     Range cells;
01847     rval = mb->get_entities_by_dimension( file_set, 2, cells );MB_CHK_ERR( rval );
01848 
01849     verts.clear();
01850     rval = mb->get_connectivity( cells, verts );MB_CHK_ERR( rval );
01851 
01852     Range modifiedCells;  // will be deleted at the end; keep the gid
01853     Range newCells;
01854 
01855     for( Range::iterator cit = cells.begin(); cit != cells.end(); ++cit )
01856     {
01857         EntityHandle cell          = *cit;
01858         const EntityHandle* connec = NULL;
01859         int num_verts              = 0;
01860         rval                       = mb->get_connectivity( cell, connec, num_verts );MB_CHK_SET_ERR( rval, "Failed to get connectivity" );
01861 
01862         std::vector< EntityHandle > newConnec;
01863         newConnec.push_back( connec[0] );  // at least one vertex
01864         int index    = 0;
01865         int new_size = 1;
01866         while( index < num_verts - 2 )
01867         {
01868             int next_index = ( index + 1 );
01869             if( connec[next_index] != newConnec[new_size - 1] )
01870             {
01871                 newConnec.push_back( connec[next_index] );
01872                 new_size++;
01873             }
01874             index++;
01875         }
01876         // add the last one only if different from previous and first node
01877         if( ( connec[num_verts - 1] != connec[num_verts - 2] ) && ( connec[num_verts - 1] != connec[0] ) )
01878         {
01879             newConnec.push_back( connec[num_verts - 1] );
01880             new_size++;
01881         }
01882         if( new_size < num_verts )
01883         {
01884             // cout << "new cell from " << cell << " has only " << new_size << " vertices \n";
01885             modifiedCells.insert( cell );
01886             // create a new cell with type triangle, quad or polygon
01887             EntityType type = MBTRI;
01888             if( new_size == 3 )
01889                 type = MBTRI;
01890             else if( new_size == 4 )
01891                 type = MBQUAD;
01892             else if( new_size > 4 )
01893                 type = MBPOLYGON;
01894 
01895             // create new cell
01896             EntityHandle newCell;
01897             rval = mb->create_element( type, &newConnec[0], new_size, newCell );MB_CHK_SET_ERR( rval, "Failed to create new cell" );
01898             // set the old id to the new element
01899             newCells.insert( newCell );
01900             double value;  // use the same value to reset the tags, even if the tags are int (like Global ID)
01901             for( size_t i = 0; i < tagList.size(); i++ )
01902             {
01903                 rval = mb->tag_get_data( tagList[i], &cell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to get tag value" );
01904                 rval = mb->tag_set_data( tagList[i], &newCell, 1, (void*)( &value ) );MB_CHK_SET_ERR( rval, "Failed to set tag value on new cell" );
01905             }
01906         }
01907     }
01908 
01909     rval = mb->remove_entities( file_set, modifiedCells );MB_CHK_SET_ERR( rval, "Failed to remove old cells from file set" );
01910     rval = mb->delete_entities( modifiedCells );MB_CHK_SET_ERR( rval, "Failed to delete old cells" );
01911     rval = mb->add_entities( file_set, newCells );MB_CHK_SET_ERR( rval, "Failed to add new cells to file set" );
01912     rval = mb->add_entities( file_set, verts );MB_CHK_SET_ERR( rval, "Failed to add verts to the file set" );
01913 
01914     return MB_SUCCESS;
01915 }
01916 
01917 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines