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