MOAB: Mesh Oriented datABase  (version 5.4.0)
Intx2Mesh.cpp
Go to the documentation of this file.
00001 /*
00002  * Intx2Mesh.cpp
00003  *
00004  *  Created on: Oct 2, 2012
00005  */
00006 
00007 #include "moab/IntxMesh/Intx2Mesh.hpp"
00008 #ifdef MOAB_HAVE_MPI
00009 #include "moab/ParallelComm.hpp"
00010 #include "MBParallelConventions.h"
00011 #include "moab/ParallelMergeMesh.hpp"
00012 #endif /* MOAB_HAVE_MPI */
00013 #include "MBTagConventions.hpp"
00014 // this is for DBL_MAX
00015 #include <cfloat>
00016 #include <queue>
00017 #include <sstream>
00018 #include "moab/GeomUtil.hpp"
00019 #include "moab/AdaptiveKDTree.hpp"
00020 
00021 namespace moab
00022 {
00023 
00024 #ifdef ENABLE_DEBUG
00025 int Intx2Mesh::dbg_1 = 0;
00026 #endif
00027 
00028 Intx2Mesh::Intx2Mesh( Interface* mbimpl )
00029     : mb( mbimpl ), mbs1( 0 ), mbs2( 0 ), outSet( 0 ), gid( 0 ), TgtFlagTag( 0 ), tgtParentTag( 0 ), srcParentTag( 0 ),
00030       countTag( 0 ), srcNeighTag( 0 ), tgtNeighTag( 0 ), neighTgtEdgeTag( 0 ), orgSendProcTag( 0 ), tgtConn( NULL ),
00031       srcConn( NULL ), epsilon_1( 0.0 ), epsilon_area( 0.0 ), box_error( 0.0 ), localRoot( 0 ), my_rank( 0 )
00032 #ifdef MOAB_HAVE_MPI
00033       ,
00034       parcomm( NULL ), remote_cells( NULL ), remote_cells_with_tracers( NULL )
00035 #endif
00036       ,
00037       max_edges_1( 0 ), max_edges_2( 0 ), counting( 0 )
00038 {
00039     gid = mbimpl->globalId_tag();
00040 }
00041 
00042 Intx2Mesh::~Intx2Mesh()
00043 {
00044     // TODO Auto-generated destructor stub
00045 #ifdef MOAB_HAVE_MPI
00046     if( remote_cells )
00047     {
00048         delete remote_cells;
00049         remote_cells = NULL;
00050     }
00051 #endif
00052 }
00053 ErrorCode Intx2Mesh::FindMaxEdgesInSet( EntityHandle eset, int& max_edges )
00054 {
00055     Range cells;
00056     ErrorCode rval = mb->get_entities_by_dimension( eset, 2, cells );MB_CHK_ERR( rval );
00057 
00058     max_edges = 0;  // can be 0 for point clouds
00059     for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
00060     {
00061         EntityHandle cell = *cit;
00062         const EntityHandle* conn4;
00063         int nnodes = 3;
00064         rval       = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
00065         if( nnodes > max_edges ) max_edges = nnodes;
00066     }
00067     // if in parallel, communicate the actual max_edges; it is not needed for tgt mesh (to be
00068     // global) but it is better to be consistent
00069 #ifdef MOAB_HAVE_MPI
00070     if( parcomm )
00071     {
00072         int local_max_edges = max_edges;
00073         // now reduce max_edges over all processors
00074         int mpi_err =
00075             MPI_Allreduce( &local_max_edges, &max_edges, 1, MPI_INT, MPI_MAX, parcomm->proc_config().proc_comm() );
00076         if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00077     }
00078 #endif
00079 
00080     return MB_SUCCESS;
00081 }
00082 ErrorCode Intx2Mesh::FindMaxEdges( EntityHandle set1, EntityHandle set2 )
00083 {
00084     ErrorCode rval = FindMaxEdgesInSet( set1, max_edges_1 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 1" );
00085     rval = FindMaxEdgesInSet( set2, max_edges_2 );MB_CHK_SET_ERR( rval, "can't determine max_edges in set 2" );
00086 
00087     return MB_SUCCESS;
00088 }
00089 
00090 ErrorCode Intx2Mesh::createTags()
00091 {
00092     if( tgtParentTag ) mb->tag_delete( tgtParentTag );
00093     if( srcParentTag ) mb->tag_delete( srcParentTag );
00094     if( countTag ) mb->tag_delete( countTag );
00095 
00096     unsigned char def_data_bit = 0;  // unused by default
00097     // maybe the tgt tag is better to be deleted every time, and recreated;
00098     // or is it easy to set all values to something again? like 0?
00099     ErrorCode rval = mb->tag_get_handle( "tgtFlag", 1, MB_TYPE_BIT, TgtFlagTag, MB_TAG_CREAT, &def_data_bit );MB_CHK_SET_ERR( rval, "can't get tgt flag tag" );
00100 
00101     // create tgt edges if they do not exist yet; so when they are looked upon, they are found
00102     // this is the only call that is potentially NlogN, in the whole method
00103     rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
00104 
00105     // now, create a map from each edge to a list of potential new nodes on a tgt edge
00106     // this memory has to be cleaned up
00107     // change it to a vector, and use the index in range of tgt edges
00108     int indx = 0;
00109     extraNodesVec.resize( TgtEdges.size() );
00110     for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
00111     {
00112         std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
00113         extraNodesVec[indx]             = nv;
00114     }
00115 
00116     int defaultInt = -1;
00117 
00118     rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
00119                                &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
00120 
00121     rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
00122                                &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
00123 
00124     rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
00125 
00126     // for each cell in set 1, determine its neigh in set 1 (could be null too)
00127     // for each cell in set 2, determine its neigh in set 2 (if on boundary, could be 0)
00128     rval = DetermineOrderedNeighbors( mbs1, max_edges_1, srcNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 1" );
00129     rval = DetermineOrderedNeighbors( mbs2, max_edges_2, tgtNeighTag );MB_CHK_SET_ERR( rval, "can't determine neighbors for set 2" );
00130 
00131     // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
00132     // them each time edges were for sure created before (tgtEdges)
00133     std::vector< EntityHandle > zeroh( max_edges_2, 0 );
00134     // if we have a tag with this name, it could be of a different size, so delete it
00135     rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
00136     if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
00137     rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
00138                                MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
00139     for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
00140     {
00141         EntityHandle tgtCell = *rit;
00142         int num_nodes        = 0;
00143         rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
00144         // account for padded polygons
00145         while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
00146             num_nodes--;
00147 
00148         int i = 0;
00149         for( i = 0; i < num_nodes; i++ )
00150         {
00151             EntityHandle v[2] = { tgtConn[i],
00152                                   tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
00153             std::vector< EntityHandle > adj_entities;
00154             rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
00155             if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval;  // get out , big error
00156             zeroh[i] = adj_entities[0];                                       // should be only one edge between 2 nodes
00157             // also, even if number of edges is less than max_edges_2, they will be ignored, even if
00158             // the tag is dense
00159         }
00160         // now set the value of the tag
00161         rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
00162     }
00163     return MB_SUCCESS;
00164 }
00165 
00166 ErrorCode Intx2Mesh::DetermineOrderedNeighbors( EntityHandle inputSet, int max_edges, Tag& neighTag )
00167 {
00168     Range cells;
00169     ErrorCode rval = mb->get_entities_by_dimension( inputSet, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells in set" );
00170 
00171     std::vector< EntityHandle > neighbors( max_edges );
00172     std::vector< EntityHandle > zeroh( max_edges, 0 );
00173     // nameless tag, as the name is not important; we will have 2 related tags, but one on tgt mesh,
00174     // one on src mesh
00175     rval = mb->tag_get_handle( "", max_edges, MB_TYPE_HANDLE, neighTag, MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create neighbors tag" );
00176 
00177     for( Range::iterator cit = cells.begin(); cit != cells.end(); cit++ )
00178     {
00179         EntityHandle cell = *cit;
00180         int nnodes        = 3;
00181         // will get the nnodes ordered neighbors;
00182         // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
00183         const EntityHandle* conn4;
00184         rval = mb->get_connectivity( cell, conn4, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity of a cell" );
00185         int nsides = nnodes;
00186         // account for possible padded polygons
00187         while( conn4[nsides - 2] == conn4[nsides - 1] && nsides > 3 )
00188             nsides--;
00189 
00190         for( int i = 0; i < nsides; i++ )
00191         {
00192             EntityHandle v[2];
00193             v[0] = conn4[i];
00194             v[1] = conn4[( i + 1 ) % nsides];
00195             // get all cells adjacent to these 2 vertices on the edge
00196             std::vector< EntityHandle > adjcells;
00197             std::vector< EntityHandle > cellsInSet;
00198             rval = mb->get_adjacencies( v, 2, 2, false, adjcells, Interface::INTERSECT );MB_CHK_SET_ERR( rval, "can't adjacency to 2 verts" );
00199             // now look for the cells contained in the input set;
00200             // the input set should be a correct mesh, not overlapping cells, and manifold
00201             size_t siz = adjcells.size();
00202             for( size_t j = 0; j < siz; j++ )
00203                 if( mb->contains_entities( inputSet, &( adjcells[j] ), 1 ) ) cellsInSet.push_back( adjcells[j] );
00204             siz = cellsInSet.size();
00205 
00206             if( siz > 2 )
00207             {
00208                 std::cout << "non manifold mesh, error" << mb->list_entities( &( cellsInSet[0] ), cellsInSet.size() )
00209                           << "\n";MB_CHK_SET_ERR( MB_FAILURE, "non-manifold input mesh set" );  // non-manifold
00210             }
00211             if( siz == 1 )
00212             {
00213                 // it must be the border of the input mesh;
00214                 neighbors[i] = 0;  // we are guaranteed that ids are !=0; this is marking a border
00215                 // borders do not appear for a sphere in serial, but they do appear for
00216                 // parallel processing anyway
00217                 continue;
00218             }
00219             // here siz ==2, it is either the first or second
00220             if( cell == cellsInSet[0] )
00221                 neighbors[i] = cellsInSet[1];
00222             else
00223                 neighbors[i] = cellsInSet[0];
00224         }
00225         // fill the rest with 0
00226         for( int i = nsides; i < max_edges; i++ )
00227             neighbors[i] = 0;
00228         // now simply set the neighbors tag; the last few positions will not be used, but for
00229         // simplicity will keep them all (MAXEDGES)
00230         rval = mb->tag_set_data( neighTag, &cell, 1, &neighbors[0] );MB_CHK_SET_ERR( rval, "can't set neigh tag" );
00231     }
00232     return MB_SUCCESS;
00233 }
00234 
00235 // slow interface; this will not do the advancing front trick
00236 // some are triangles, some are quads, some are polygons ...
00237 ErrorCode Intx2Mesh::intersect_meshes_kdtree( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
00238 {
00239     ErrorCode rval;
00240     mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
00241     mbs2   = mbset2;
00242     outSet = outputSet;
00243     rval   = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
00244     rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
00245     // from create tags, copy relevant ones
00246     if( tgtParentTag ) mb->tag_delete( tgtParentTag );
00247     if( srcParentTag ) mb->tag_delete( srcParentTag );
00248     if( countTag ) mb->tag_delete( countTag );
00249 
00250     // create tgt edges if they do not exist yet; so when they are looked upon, they are found
00251     // this is the only call that is potentially NlogN, in the whole method
00252     rval = mb->get_adjacencies( rs2, 1, true, TgtEdges, Interface::UNION );MB_CHK_SET_ERR( rval, "can't get adjacent tgt edges" );
00253 
00254     int indx = 0;
00255     extraNodesVec.resize( TgtEdges.size() );
00256     for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
00257     {
00258         std::vector< EntityHandle >* nv = new std::vector< EntityHandle >;
00259         extraNodesVec[indx]             = nv;
00260     }
00261 
00262     int defaultInt = -1;
00263 
00264     rval = mb->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
00265                                &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
00266 
00267     rval = mb->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag, MB_TAG_DENSE | MB_TAG_CREAT,
00268                                &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
00269 
00270     rval = mb->tag_get_handle( "Counting", 1, MB_TYPE_INTEGER, countTag, MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create Counting tag" );
00271 
00272     // for tgt cells, save a dense tag with the bordering edges, so we do not have to search for
00273     // them each time edges were for sure created before (tgtEdges)
00274     // if we have a tag with this name, it could be of a different size, so delete it
00275     rval = mb->tag_get_handle( "__tgtEdgeNeighbors", neighTgtEdgeTag );
00276     if( rval == MB_SUCCESS && neighTgtEdgeTag ) mb->tag_delete( neighTgtEdgeTag );
00277     std::vector< EntityHandle > zeroh( max_edges_2, 0 );
00278     rval = mb->tag_get_handle( "__tgtEdgeNeighbors", max_edges_2, MB_TYPE_HANDLE, neighTgtEdgeTag,
00279                                MB_TAG_DENSE | MB_TAG_CREAT, &zeroh[0] );MB_CHK_SET_ERR( rval, "can't create tgt edge neighbors tag" );
00280     for( Range::iterator rit = rs2.begin(); rit != rs2.end(); rit++ )
00281     {
00282         EntityHandle tgtCell = *rit;
00283         int num_nodes        = 0;
00284         rval                 = mb->get_connectivity( tgtCell, tgtConn, num_nodes );MB_CHK_SET_ERR( rval, "can't get  tgt conn" );
00285         // account for padded polygons
00286         while( tgtConn[num_nodes - 2] == tgtConn[num_nodes - 1] && num_nodes > 3 )
00287             num_nodes--;
00288 
00289         int i = 0;
00290         for( i = 0; i < num_nodes; i++ )
00291         {
00292             EntityHandle v[2] = { tgtConn[i],
00293                                   tgtConn[( i + 1 ) % num_nodes] };  // this is fine even for padded polygons
00294             std::vector< EntityHandle > adj_entities;
00295             rval = mb->get_adjacencies( v, 2, 1, false, adj_entities, Interface::INTERSECT );
00296             if( rval != MB_SUCCESS || adj_entities.size() < 1 ) return rval;  // get out , big error
00297             zeroh[i] = adj_entities[0];                                       // should be only one edge between 2 nodes
00298             // also, even if number of edges is less than max_edges_2, they will be ignored, even if
00299             // the tag is dense
00300         }
00301         // now set the value of the tag
00302         rval = mb->tag_set_data( neighTgtEdgeTag, &tgtCell, 1, &( zeroh[0] ) );MB_CHK_SET_ERR( rval, "can't set edge tgt tag" );
00303     }
00304 
00305     // create the kd tree on source cells, and intersect all targets in an expensive loop
00306     // build a kd tree with the rs1 (source) cells
00307     AdaptiveKDTree kd( mb );
00308     EntityHandle tree_root = 0;
00309     rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
00310 
00311     // find out max edge on source mesh;
00312     double max_length = 0;
00313     {
00314         std::vector< double > coords;
00315         coords.resize( 3 * max_edges_1 );
00316         for( Range::iterator it = rs1.begin(); it != rs1.end(); it++ )
00317         {
00318             const EntityHandle* conn = NULL;
00319             int nnodes;
00320             rval = mb->get_connectivity( *it, conn, nnodes );MB_CHK_SET_ERR( rval, "can't get connectivity" );
00321             while( conn[nnodes - 2] == conn[nnodes - 1] && nnodes > 3 )
00322                 nnodes--;
00323             rval = mb->get_coords( conn, nnodes, &coords[0] );MB_CHK_SET_ERR( rval, "can't get coordinates" );
00324             for( int j = 0; j < nnodes; j++ )
00325             {
00326                 int next = ( j + 1 ) % nnodes;
00327                 double leng;
00328                 leng = ( coords[3 * j] - coords[3 * next] ) * ( coords[3 * j] - coords[3 * next] ) +
00329                        ( coords[3 * j + 1] - coords[3 * next + 1] ) * ( coords[3 * j + 1] - coords[3 * next + 1] ) +
00330                        ( coords[3 * j + 2] - coords[3 * next + 2] ) * ( coords[3 * j + 2] - coords[3 * next + 2] );
00331                 leng = sqrt( leng );
00332                 if( leng > max_length ) max_length = leng;
00333             }
00334         }
00335     }
00336     // maximum sag on a spherical mesh make sense only for intx on a sphere, with radius 1 :(
00337     double tolerance = 1.e-15;
00338     if( max_length < 1. )
00339     {
00340         // basically, the sag for an arc of length max_length on a circle of radius 1
00341         tolerance = 1. - sqrt( 1 - max_length * max_length / 4 );
00342         if( box_error < tolerance ) box_error = tolerance;
00343         tolerance = 3 * tolerance;  // we use it for gnomonic plane too, projected sag could be =* sqrt(2.)
00344         // be more generous, use 1.5 ~= sqrt(2.)
00345 
00346         if( !my_rank )
00347         {
00348             std::cout << " max edge length: " << max_length << "  tolerance for kd tree: " << tolerance << "\n";
00349             std::cout << " box overlap tolerance: " << box_error << "\n";
00350         }
00351     }
00352     for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
00353     {
00354         EntityHandle tcell = *it;
00355         // find vertex positions
00356         const EntityHandle* conn = NULL;
00357         int nnodes               = 0;
00358         rval                     = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
00359         // find leaves close to those positions
00360         double areaTgtCell   = setup_tgt_cell( tcell, nnodes );  // this is the area in the gnomonic plane
00361         double recoveredArea = 0;
00362         std::vector< double > positions;
00363         positions.resize( nnodes * 3 );
00364         rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
00365 
00366         // distance to search will be based on average edge length
00367         double av_len = 0;
00368         for( int k = 0; k < nnodes; k++ )
00369         {
00370             int ik      = ( k + 1 ) % nnodes;
00371             double len1 = 0;
00372             for( int j = 0; j < 3; j++ )
00373             {
00374                 double len2 = positions[3 * k + j] - positions[3 * ik + j];
00375                 len1 += len2 * len2;
00376             }
00377             av_len += sqrt( len1 );
00378         }
00379         if( nnodes > 0 ) av_len /= nnodes;
00380         // find leaves within a distance from each vertex of target
00381         // in those leaves, collect all cells; we will try for an intx in there
00382         Range close_source_cells;
00383         std::vector< EntityHandle > leaves;
00384         for( int i = 0; i < nnodes; i++ )
00385         {
00386 
00387             leaves.clear();
00388             rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
00389 
00390             for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
00391             {
00392                 Range tmp;
00393                 rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
00394 
00395                 close_source_cells.merge( tmp.begin(), tmp.end() );
00396             }
00397         }
00398 #ifdef VERBOSE
00399         if( close_source_cells.empty() )
00400         {
00401             std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
00402                       << mb->id_from_handle( tcell ) << "\n";
00403         }
00404 #endif
00405         for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
00406         {
00407             EntityHandle startSrc = *it2;
00408             double area           = 0;
00409             // if area is > 0 , we have intersections
00410             double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
00411             //
00412             int nP = 0;
00413             int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
00414             int nsTgt, nsSrc;
00415             rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
00416             if( area > 0 )
00417             {
00418                 if( nP > 1 )
00419                 {  // this will also construct triangles/polygons in the new mesh, if needed
00420                     rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
00421                 }
00422                 recoveredArea += area;
00423             }
00424         }
00425         recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fract
00426     }
00427     // before cleaning up , we need to settle the position of the intersection points
00428     // on the boundary edges
00429     // this needs to be collective, so we should maybe wait something
00430 #ifdef MOAB_HAVE_MPI
00431     rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
00432 #endif
00433 
00434     this->clean();
00435     return MB_SUCCESS;
00436 }
00437 // main interface; this will do the advancing front trick
00438 // some are triangles, some are quads, some are polygons ...
00439 ErrorCode Intx2Mesh::intersect_meshes( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
00440 {
00441     ErrorCode rval;
00442     mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
00443     mbs2   = mbset2;
00444     outSet = outputSet;
00445 #ifdef VERBOSE
00446     std::stringstream ffs, fft;
00447     ffs << "source_rank0" << my_rank << ".vtk";
00448     rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
00449     fft << "target_rank0" << my_rank << ".vtk";
00450     rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
00451 
00452 #endif
00453     // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
00454 
00455     EntityHandle startSrc = 0, startTgt = 0;
00456 
00457     rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
00458     rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
00459     // std::cout << "rs1.size() = " << rs1.size() << " and rs2.size() = "  << rs2.size() << "\n";
00460     // std::cout.flush();
00461 
00462     createTags();  // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
00463 
00464     Range rs22 = rs2;  // a copy of the initial range; we will remove from it elements as we
00465                        // advance ; rs2 is needed for marking the polygon to the tgt parent
00466 
00467     // create the local kdd tree with source elements; will use it to search
00468     // more efficiently for the seeds in advancing front;
00469     // some of the target cells will not be covered by source cells, and they need to be eliminated
00470     // early from contention
00471 
00472     // build a kd tree with the rs1 (source) cells
00473     AdaptiveKDTree kd( mb );
00474     EntityHandle tree_root = 0;
00475     rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
00476 
00477     while( !rs22.empty() )
00478     {
00479 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
00480         if( rs22.size() < rs2.size() )
00481         {
00482             std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
00483                       << "\n";
00484             std::stringstream ffo;
00485             ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
00486             rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
00487         }
00488 #endif
00489         bool seedFound = false;
00490         for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
00491         {
00492             startTgt  = *it;
00493             int found = 0;
00494             // find vertex positions
00495             const EntityHandle* conn = NULL;
00496             int nnodes               = 0;
00497             rval                     = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
00498             // find leaves close to those positions
00499             std::vector< double > positions;
00500             positions.resize( nnodes * 3 );
00501             rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
00502             // find leaves within a distance from each vertex of target
00503             // in those leaves, collect all cells; we will try for an intx in there, instead of
00504             // looping over all rs1 cells, as before
00505             Range close_source_cells;
00506             std::vector< EntityHandle > leaves;
00507             for( int i = 0; i < nnodes; i++ )
00508             {
00509 
00510                 leaves.clear();
00511                 rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
00512 
00513                 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
00514                 {
00515                     Range tmp;
00516                     rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
00517 
00518                     close_source_cells.merge( tmp.begin(), tmp.end() );
00519                 }
00520             }
00521 
00522             for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
00523             {
00524                 startSrc    = *it2;
00525                 double area = 0;
00526                 // if area is > 0 , we have intersections
00527                 double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
00528                 //
00529                 int nP = 0;
00530                 int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
00531                 int nsTgt, nsSrc;
00532                 rval =
00533                     computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
00534                 if( area > 0 )
00535                 {
00536                     found     = 1;
00537                     seedFound = true;
00538                     break;  // found 2 elements that intersect; these will be the seeds
00539                 }
00540             }
00541             if( found )
00542                 break;
00543             else
00544             {
00545 #if defined( VERBOSE )
00546                 std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
00547                           << " not intx with any source\n";
00548 #endif
00549                 rs22.erase( startTgt );
00550             }
00551         }
00552         if( !seedFound ) continue;  // continue while(!rs22.empty())
00553 
00554         std::queue< EntityHandle > srcQueue;  // these are corresponding to Ta,
00555         srcQueue.push( startSrc );
00556         std::queue< EntityHandle > tgtQueue;
00557         tgtQueue.push( startTgt );
00558 
00559         Range toResetSrcs;  // will be used to reset src flags for every tgt quad
00560         // processed
00561 
00562         /*if (my_rank==0)
00563           dbg_1 = 1;*/
00564         unsigned char used = 1;
00565         // mark the start tgt quad as used, so it will not come back again
00566         rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
00567         while( !tgtQueue.empty() )
00568         {
00569             // flags for the side : 0 means a src cell not found on side
00570             // a paired src not found yet for the neighbors of tgt
00571             Range nextSrc[MAXEDGES];  // there are new ranges of possible next src cells for
00572                                       // seeding the side j of tgt cell
00573 
00574             EntityHandle currentTgt = tgtQueue.front();
00575             tgtQueue.pop();
00576             int nsidesTgt;                                                   // will be initialized now
00577             double areaTgtCell   = setup_tgt_cell( currentTgt, nsidesTgt );  // this is the area in the gnomonic plane
00578             double recoveredArea = 0;
00579             // get the neighbors of tgt, and if they are solved already, do not bother with that
00580             // side of tgt
00581             EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
00582             rval                                = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
00583 #ifdef ENABLE_DEBUG
00584             if( dbg_1 )
00585             {
00586                 std::cout << "Next: neighbors for current tgt ";
00587                 for( int kk = 0; kk < nsidesTgt; kk++ )
00588                 {
00589                     if( tgtNeighbors[kk] > 0 )
00590                         std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
00591                     else
00592                         std::cout << 0 << " ";
00593                 }
00594                 std::cout << std::endl;
00595             }
00596 #endif
00597             // now get the status of neighbors; if already solved, make them 0, so not to bother
00598             // anymore on that side of tgt
00599             for( int j = 0; j < nsidesTgt; j++ )
00600             {
00601                 EntityHandle tgtNeigh = tgtNeighbors[j];
00602                 unsigned char status  = 1;
00603                 if( tgtNeigh == 0 ) continue;
00604                 rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval );                     // status 0 is unused
00605                 if( 1 == status ) tgtNeighbors[j] = 0;  // so will not look anymore on this side of tgt
00606             }
00607 
00608 #ifdef ENABLE_DEBUG
00609             if( dbg_1 )
00610             {
00611                 std::cout << "reset sources: ";
00612                 for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
00613                     std::cout << mb->id_from_handle( *itr ) << " ";
00614                 std::cout << std::endl;
00615             }
00616 #endif
00617             EntityHandle currentSrc = srcQueue.front();
00618             // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
00619             // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
00620             // pairs;
00621             //  at every moment, the queue contains pairs of cells that intersect, and they form the
00622             //  "advancing front"
00623             srcQueue.pop();
00624             toResetSrcs.clear();  // empty the range of used srcs, will have to be set unused again,
00625             // at the end of tgt element processing
00626             toResetSrcs.insert( currentSrc );
00627             // mb2->set_tag_data
00628             std::queue< EntityHandle > localSrc;
00629             localSrc.push( currentSrc );
00630 #ifdef VERBOSE
00631             int countingStart = counting;
00632 #endif
00633             // will advance-front search in the neighborhood of tgt cell, until we finish processing
00634             // all
00635             //   possible src cells; localSrc queue will contain all possible src cells that cover
00636             //   the current tgt cell
00637             while( !localSrc.empty() )
00638             {
00639                 //
00640                 EntityHandle srcT = localSrc.front();
00641                 localSrc.pop();
00642                 double P[10 * MAXEDGES], area;  //
00643                 int nP           = 0;
00644                 int nb[MAXEDGES] = { 0 };
00645                 int nr[MAXEDGES] = { 0 };
00646 
00647                 int nsidesSrc;  ///
00648                 // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
00649                 // intersection points could include the vertices of initial elements
00650                 // nb [j] = 0 means no intersection on the side j for element src (markers)
00651                 // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
00652                 // tgt poly.  A potential next poly in the tgt queue is the tgt poly that is
00653                 // adjacent to this side
00654                 rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
00655                                                             nsidesTgt );MB_CHK_ERR( rval );
00656                 if( nP > 0 )
00657                 {
00658 #ifdef ENABLE_DEBUG
00659                     if( dbg_1 )
00660                     {
00661                         for( int k = 0; k < 3; k++ )
00662                         {
00663                             std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
00664                         }
00665                     }
00666 #endif
00667 
00668                     // intersection found: output P and original triangles if nP > 2
00669                     EntityHandle neighbors[MAXEDGES] = { 0 };
00670                     rval                             = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
00671                     if( rval != MB_SUCCESS )
00672                     {
00673                         std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
00674                         return MB_FAILURE;
00675                     }
00676 
00677                     // add neighbors to the localSrc queue, if they are not marked
00678                     for( int nn = 0; nn < nsidesSrc; nn++ )
00679                     {
00680                         EntityHandle neighbor = neighbors[nn];
00681                         if( neighbor > 0 && nb[nn] > 0 )  // advance across src boundary nn
00682                         {
00683                             if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
00684                             {
00685                                 localSrc.push( neighbor );
00686 #ifdef ENABLE_DEBUG
00687                                 if( dbg_1 )
00688                                 {
00689                                     std::cout << " local src elem " << mb->id_from_handle( neighbor )
00690                                               << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
00691                                     mb->list_entities( &neighbor, 1 );
00692                                 }
00693 #endif
00694                                 toResetSrcs.insert( neighbor );
00695                             }
00696                         }
00697                     }
00698                     // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
00699                     for( int nn = 0; nn < nsidesTgt; nn++ )
00700                     {
00701                         if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
00702                             nextSrc[nn].insert( srcT );  // potential src cell that can intersect
00703                                                          // the tgt neighbor nn
00704                     }
00705                     if( nP > 1 )
00706                     {  // this will also construct triangles/polygons in the new mesh, if needed
00707                         rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
00708                     }
00709 
00710                     recoveredArea += area;
00711                 }
00712 #ifdef ENABLE_DEBUG
00713                 else if( dbg_1 )
00714                 {
00715                     std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
00716                               << mb->id_from_handle( srcT ) << "\n";
00717                 }
00718 #endif
00719             }                                                               // end while (!localSrc.empty())
00720             recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fraction
00721 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
00722             if( fabs( recoveredArea ) > epsilon_1 )
00723             {
00724 #ifdef VERBOSE
00725                 std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
00726                           << " fraction error recovery:" << recoveredArea
00727                           << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
00728                           << "\n";
00729 #endif
00730             }
00731 #endif
00732             // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
00733             // mesh)
00734             rs22.erase( currentTgt );
00735             // also, look at its neighbors, and add to the seeds a next one
00736 
00737             for( int j = 0; j < nsidesTgt; j++ )
00738             {
00739                 EntityHandle tgtNeigh = tgtNeighbors[j];
00740                 if( tgtNeigh == 0 || nextSrc[j].size() == 0 )  // if tgt is bigger than src, there could be no src
00741                                                                // to advance on that side
00742                     continue;
00743                 int nsidesTgt2 = 0;
00744                 setup_tgt_cell( tgtNeigh,
00745                                 nsidesTgt2 );  // find possible intersection with src cell from nextSrc
00746                 for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
00747                 {
00748                     EntityHandle nextB = *nit;
00749                     // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
00750                     // src quad
00751                     double P[10 * MAXEDGES], area;  //
00752                     int nP           = 0;
00753                     int nb[MAXEDGES] = { 0 };
00754                     int nr[MAXEDGES] = { 0 };
00755 
00756                     int nsidesSrc;  ///
00757                     rval = computeIntersectionBetweenTgtAndSrc(
00758                         /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
00759                     if( area > 0 )
00760                     {
00761                         tgtQueue.push( tgtNeigh );
00762                         srcQueue.push( nextB );
00763 #ifdef ENABLE_DEBUG
00764                         if( dbg_1 )
00765                             std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
00766                                       << mb->id_from_handle( nextB ) << std::endl;
00767 #endif
00768                         rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
00769                         break;  // so we are done with this side of tgt, we have found a proper next
00770                                 // seed
00771                     }
00772                 }
00773             }
00774 
00775         }  // end while (!tgtQueue.empty())
00776     }
00777 #ifdef ENABLE_DEBUG
00778     if( dbg_1 )
00779     {
00780         for( int k = 0; k < 6; k++ )
00781             mout_1[k].close();
00782     }
00783 #endif
00784     // before cleaning up , we need to settle the position of the intersection points
00785     // on the boundary edges
00786     // this needs to be collective, so we should maybe wait something
00787 #ifdef MOAB_HAVE_MPI
00788     rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
00789 #endif
00790 
00791     this->clean();
00792     return MB_SUCCESS;
00793 }
00794 
00795 // clean some memory allocated
00796 void Intx2Mesh::clean()
00797 {
00798     //
00799     int indx = 0;
00800     for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
00801     {
00802         delete extraNodesVec[indx];
00803     }
00804     // extraNodesMap.clear();
00805     extraNodesVec.clear();
00806     // also, delete some bit tags, used to mark processed tgts and srcs
00807     mb->tag_delete( TgtFlagTag );
00808     counting = 0;  // reset counting to original value
00809 }
00810 // this method will reduce number of nodes, collapse edges that are of length 0
00811 // so a polygon like 428 431 431 will become a line 428 431
00812 // or something like 428 431 431 531 -> 428 431 531
00813 void Intx2Mesh::correct_polygon( EntityHandle* nodes, int& nP )
00814 {
00815     int i = 0;
00816     while( i < nP )
00817     {
00818         int nextIndex = ( i + 1 ) % nP;
00819         if( nodes[i] == nodes[nextIndex] )
00820         {
00821 #ifdef ENABLE_DEBUG
00822             // we need to reduce nP, and collapse nodes
00823             if( dbg_1 )
00824             {
00825                 std::cout << " nodes duplicated in list: ";
00826                 for( int j = 0; j < nP; j++ )
00827                     std::cout << nodes[j] << " ";
00828                 std::cout << "\n";
00829                 std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
00830                           << "\n";
00831             }
00832 #endif
00833             // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
00834             // then next thing does nothing
00835             //  (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
00836             for( int k = i; k < nP - 1; k++ )
00837                 nodes[k] = nodes[k + 1];
00838             nP--;  // decrease the number of nodes; also, decrease i, just if we may need to check
00839                    // again
00840             i--;
00841         }
00842         i++;
00843     }
00844     return;
00845 }
00846 #ifdef MOAB_HAVE_MPI
00847 ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
00848 {
00849     localEnts.clear();
00850     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
00851 
00852     rval                = mb->get_connectivity( localEnts, local_verts );
00853     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
00854 
00855     assert( parcomm != NULL );
00856 
00857     // get the position of local vertices, and decide local boxes (allBoxes...)
00858     double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX };
00859     double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX };
00860 
00861     std::vector< double > coords( 3 * num_local_verts );
00862     rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " );
00863 
00864     for( int i = 0; i < num_local_verts; i++ )
00865     {
00866         for( int k = 0; k < 3; k++ )
00867         {
00868             double val = coords[3 * i + k];
00869             if( val < bmin[k] ) bmin[k] = val;
00870             if( val > bmax[k] ) bmax[k] = val;
00871         }
00872     }
00873     int numprocs = parcomm->proc_config().proc_size();
00874     allBoxes.resize( 6 * numprocs );
00875 
00876     my_rank = parcomm->proc_config().proc_rank();
00877     for( int k = 0; k < 3; k++ )
00878     {
00879         allBoxes[6 * my_rank + k]     = bmin[k];
00880         allBoxes[6 * my_rank + 3 + k] = bmax[k];
00881     }
00882 
00883     // now communicate to get all boxes
00884     int mpi_err;
00885 #if( MPI_VERSION >= 2 )
00886     // use "in place" option
00887     mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
00888                              parcomm->proc_config().proc_comm() );
00889 #else
00890     {
00891         std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
00892         mpi_err  = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
00893                                   parcomm->proc_config().proc_comm() );
00894         allBoxes = allBoxes_tmp;
00895     }
00896 #endif
00897     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00898 
00899     if( my_rank == 0 )
00900     {
00901         std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
00902                   << " on second mesh \n";
00903         for( int i = 0; i < numprocs; i++ )
00904         {
00905             std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " "
00906                       << allBoxes[6 * i + 2] << " \n";
00907             std::cout << "        box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " "
00908                       << allBoxes[6 * i + 5] << " \n";
00909         }
00910     }
00911 
00912     return MB_SUCCESS;
00913 }
00914 ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set )
00915 {
00916     // compute the bounding box on each proc
00917     assert( parcomm != NULL );
00918 
00919     localEnts.clear();
00920     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
00921 
00922     Tag dpTag = 0;
00923     std::string tag_name( "DP" );
00924     rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" );
00925 
00926     EntityHandle dum = 0;
00927     Tag corrTag;
00928     rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );ERRORR( rval, "can't get CORR tag" );
00929     // get all local verts
00930     Range local_verts;
00931     rval                = mb->get_connectivity( localEnts, local_verts );
00932     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
00933 
00934     rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" );
00935 
00936     std::vector< int > gids( num_local_verts );
00937     rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
00938 
00939     // now see the departure points; to what boxes should we send them?
00940     std::vector< double > dep_points( 3 * num_local_verts );
00941     rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" );
00942     // ranges to send to each processor; will hold vertices and elements (quads?)
00943     // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
00944     std::map< int, Range > Rto;
00945     int numprocs = parcomm->proc_config().proc_size();
00946 
00947     for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit )
00948     {
00949         EntityHandle q = *eit;
00950         const EntityHandle* conn4;
00951         int num_nodes;
00952         rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
00953         CartVect qbmin( DBL_MAX );
00954         CartVect qbmax( -DBL_MAX );
00955         for( int i = 0; i < num_nodes; i++ )
00956         {
00957             EntityHandle v = conn4[i];
00958             size_t index   = local_verts.find( v ) - local_verts.begin();
00959             CartVect dp( &dep_points[3 * index] );  // will use constructor
00960             for( int j = 0; j < 3; j++ )
00961             {
00962                 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
00963                 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
00964             }
00965         }
00966         for( int p = 0; p < numprocs; p++ )
00967         {
00968             CartVect bbmin( &allBoxes[6 * p] );
00969             CartVect bbmax( &allBoxes[6 * p + 3] );
00970             if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
00971             {
00972                 Rto[p].insert( q );
00973             }
00974         }
00975     }
00976 
00977     // now, build TLv and TLq, for each p
00978     size_t numq = 0;
00979     size_t numv = 0;
00980     for( int p = 0; p < numprocs; p++ )
00981     {
00982         if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
00983         Range& range_to_P = Rto[p];
00984         // add the vertices to it
00985         if( range_to_P.empty() ) continue;  // nothing to send to proc p
00986         Range vertsToP;
00987         rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
00988         numq = numq + range_to_P.size();
00989         numv = numv + vertsToP.size();
00990         range_to_P.merge( vertsToP );
00991     }
00992     TupleList TLv;
00993     TupleList TLq;
00994     TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
00995     TLv.enableWriteAccess();
00996 
00997     int sizeTuple = 2 + max_edges_1;  // determined earlier, for src, first mesh
00998     TLq.initialize( 2 + max_edges_1, 0, 1, 0,
00999                     numq );  // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh
01000     TLq.enableWriteAccess();
01001 #ifdef VERBOSE
01002     std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
01003 #endif
01004     for( int to_proc = 0; to_proc < numprocs; to_proc++ )
01005     {
01006         if( to_proc == (int)my_rank ) continue;
01007         Range& range_to_P = Rto[to_proc];
01008         Range V           = range_to_P.subset_by_type( MBVERTEX );
01009 
01010         for( Range::iterator it = V.begin(); it != V.end(); ++it )
01011         {
01012             EntityHandle v       = *it;
01013             unsigned int index   = local_verts.find( v ) - local_verts.begin();
01014             int n                = TLv.get_n();
01015             TLv.vi_wr[2 * n]     = to_proc;                // send to processor
01016             TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
01017             TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
01018             TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
01019             TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
01020             TLv.inc_n();
01021         }
01022         // also, prep the quad for sending ...
01023         Range Q = range_to_P.subset_by_dimension( 2 );
01024         for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
01025         {
01026             EntityHandle q = *it;
01027             int global_id;
01028             rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
01029             int n                        = TLq.get_n();
01030             TLq.vi_wr[sizeTuple * n]     = to_proc;    //
01031             TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
01032             const EntityHandle* conn4;
01033             int num_nodes;
01034             rval = mb->get_connectivity( q, conn4,
01035                                          num_nodes );  // could be up to MAXEDGES, but it is limited by max_edges_1
01036             ERRORR( rval, "can't get connectivity for cell" );
01037             if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
01038             for( int i = 0; i < num_nodes; i++ )
01039             {
01040                 EntityHandle v                   = conn4[i];
01041                 unsigned int index               = local_verts.find( v ) - local_verts.begin();
01042                 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
01043             }
01044             for( int k = num_nodes; k < max_edges_1; k++ )
01045             {
01046                 TLq.vi_wr[sizeTuple * n + 2 + k] =
01047                     0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
01048             }
01049             TLq.vul_wr[n] = q;  // save here the entity handle, it will be communicated back
01050             // maybe we should forget about global ID
01051             TLq.inc_n();
01052         }
01053     }
01054 
01055     // now we are done populating the tuples; route them to the appropriate processors
01056     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01057     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
01058     // the elements are already in localEnts;
01059 
01060     // maps from global ids to new vertex and quad handles, that are added
01061     std::map< int, EntityHandle > globalID_to_handle;
01062     /*std::map<int, EntityHandle> globalID_to_eh;*/
01063     globalID_to_eh.clear();  // need for next iteration
01064     // now, look at every TLv, and see if we have to create a vertex there or not
01065     int n = TLv.get_n();  // the size of the points received
01066     for( int i = 0; i < n; i++ )
01067     {
01068         int globalId = TLv.vi_rd[2 * i + 1];
01069         if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01070         {
01071             EntityHandle new_vert;
01072             double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
01073             rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01074             globalID_to_handle[globalId] = new_vert;
01075         }
01076     }
01077 
01078     // now, all dep points should be at their place
01079     // look in the local list of q for this proc, and create all those quads and vertices if needed
01080     // it may be an overkill, but because it does not involve communication, we do it anyway
01081     Range& local  = Rto[my_rank];
01082     Range local_q = local.subset_by_dimension( 2 );
01083     // the local should have all the vertices in local_verts
01084     for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
01085     {
01086         EntityHandle q = *it;
01087         int nnodes;
01088         const EntityHandle* conn4;
01089         rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " );
01090         EntityHandle new_conn[MAXEDGES];
01091         for( int i = 0; i < nnodes; i++ )
01092         {
01093             EntityHandle v1    = conn4[i];
01094             unsigned int index = local_verts.find( v1 ) - local_verts.begin();
01095             int globalId       = gids[index];
01096             if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01097             {
01098                 // we need to create that vertex, at this position dep_points
01099                 double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
01100                 EntityHandle new_vert;
01101                 rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01102                 globalID_to_handle[globalId] = new_vert;
01103             }
01104             new_conn[i] = globalID_to_handle[gids[index]];
01105         }
01106         EntityHandle new_element;
01107         //
01108         EntityType entType = MBQUAD;
01109         if( nnodes > 4 ) entType = MBPOLYGON;
01110         if( nnodes < 4 ) entType = MBTRI;
01111 
01112         rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " );
01113         rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
01114         int gid_el;
01115         // get the global ID of the initial quad
01116         rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
01117         globalID_to_eh[gid_el] = new_element;
01118         // is this redundant or not?
01119         rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" );
01120         // set the global id on new elem
01121         rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" );
01122     }
01123     // now look at all elements received through; we do not want to duplicate them
01124     n = TLq.get_n();  // number of elements received by this processor
01125     // form the remote cells, that will be used to send the tracer info back to the originating proc
01126     remote_cells = new TupleList();
01127     remote_cells->initialize( 2, 0, 1, 0, n );  // will not have tracer data anymore
01128     remote_cells->enableWriteAccess();
01129     for( int i = 0; i < n; i++ )
01130     {
01131         int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
01132         int from_proc  = TLq.vi_wr[sizeTuple * i];
01133         // do we already have a quad with this global ID, represented?
01134         if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
01135         {
01136             // construct the conn quad
01137             EntityHandle new_conn[MAXEDGES];
01138             int nnodes = -1;
01139             for( int j = 0; j < max_edges_1; j++ )
01140             {
01141                 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
01142                 if( vgid == 0 )
01143                     new_conn[j] = 0;
01144                 else
01145                 {
01146                     assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
01147                     new_conn[j] = globalID_to_handle[vgid];
01148                     nnodes      = j + 1;  // nodes are at the beginning, and are variable number
01149                 }
01150             }
01151             EntityHandle new_element;
01152             //
01153             EntityType entType = MBQUAD;
01154             if( nnodes > 4 ) entType = MBPOLYGON;
01155             if( nnodes < 4 ) entType = MBTRI;
01156             rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
01157             globalID_to_eh[globalIdEl] = new_element;
01158             rval                       = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
01159             /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/
01160             remote_cells->vi_wr[2 * i]     = from_proc;
01161             remote_cells->vi_wr[2 * i + 1] = globalIdEl;
01162             //     remote_cells->vr_wr[i] = 0.; // no contribution yet sent back
01163             remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
01164             remote_cells->inc_n();
01165             // set the global id on new elem
01166             rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" );
01167         }
01168     }
01169     // order the remote cells tuple list, with the global id, because we will search in it
01170     // remote_cells->print("remote_cells before sorting");
01171     moab::TupleList::buffer sort_buffer;
01172     sort_buffer.buffer_init( n );
01173     remote_cells->sort( 1, &sort_buffer );
01174     sort_buffer.reset();
01175     return MB_SUCCESS;
01176 }
01177 
01178 // this algorithm assumes lagr set is already created, and some elements will be coming from
01179 // other procs, and populate the covering_set
01180 // we need to keep in a tuple list the remote cells from other procs, because we need to send back
01181 // the intersection info (like area of the intx polygon, and the current concentration) maybe total
01182 // mass in that intx
01183 ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set )
01184 {
01185     EntityHandle dum = 0;
01186 
01187     Tag corrTag;
01188     ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
01189     // start copy from 2nd alg
01190     // compute the bounding box on each proc
01191     assert( parcomm != NULL );
01192     if( 1 == parcomm->proc_config().proc_size() )
01193     {
01194         covering_set = lagr_set;  // nothing to communicate, it must be serial
01195         return MB_SUCCESS;
01196     }
01197 
01198     // get all local verts
01199     Range local_verts;
01200     rval                = mb->get_connectivity( localEnts, local_verts );
01201     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
01202 
01203     std::vector< int > gids( num_local_verts );
01204     rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
01205 
01206     Range localDepCells;
01207     rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" );
01208 
01209     // get all lagr verts (departure vertices)
01210     Range lagr_verts;
01211     rval = mb->get_connectivity( localDepCells, lagr_verts );  // they should be created in
01212     // the same order as the euler vertices
01213     int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" );
01214 
01215     // now see the departure points position; to what boxes should we send them?
01216     std::vector< double > dep_points( 3 * num_lagr_verts );
01217     rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" );
01218     // ranges to send to each processor; will hold vertices and elements (quads?)
01219     // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
01220     std::map< int, Range > Rto;
01221     int numprocs = parcomm->proc_config().proc_size();
01222 
01223     for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
01224     {
01225         EntityHandle q = *eit;
01226         const EntityHandle* conn4;
01227         int num_nodes;
01228         rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
01229         CartVect qbmin( DBL_MAX );
01230         CartVect qbmax( -DBL_MAX );
01231         for( int i = 0; i < num_nodes; i++ )
01232         {
01233             EntityHandle v = conn4[i];
01234             int index      = lagr_verts.index( v );
01235             assert( -1 != index );
01236             CartVect dp( &dep_points[3 * index] );  // will use constructor
01237             for( int j = 0; j < 3; j++ )
01238             {
01239                 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
01240                 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
01241             }
01242         }
01243         for( int p = 0; p < numprocs; p++ )
01244         {
01245             CartVect bbmin( &allBoxes[6 * p] );
01246             CartVect bbmax( &allBoxes[6 * p + 3] );
01247             if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) )
01248             {
01249                 Rto[p].insert( q );
01250             }
01251         }
01252     }
01253 
01254     // now, build TLv and TLq, for each p
01255     size_t numq = 0;
01256     size_t numv = 0;
01257     for( int p = 0; p < numprocs; p++ )
01258     {
01259         if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
01260         Range& range_to_P = Rto[p];
01261         // add the vertices to it
01262         if( range_to_P.empty() ) continue;  // nothing to send to proc p
01263         Range vertsToP;
01264         rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
01265         numq = numq + range_to_P.size();
01266         numv = numv + vertsToP.size();
01267         range_to_P.merge( vertsToP );
01268     }
01269     TupleList TLv;
01270     TupleList TLq;
01271     TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
01272     TLv.enableWriteAccess();
01273 
01274     int sizeTuple = 2 + max_edges_1;  // max edges could be up to MAXEDGES :) for polygons
01275     TLq.initialize( 2 + max_edges_1, 0, 1, 0,
01276                     numq );  // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v)
01277     // send also the corresponding tgt cell it will come to
01278     TLq.enableWriteAccess();
01279 #ifdef VERBOSE
01280     std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
01281 #endif
01282 
01283     for( int to_proc = 0; to_proc < numprocs; to_proc++ )
01284     {
01285         if( to_proc == (int)my_rank ) continue;
01286         Range& range_to_P = Rto[to_proc];
01287         Range V           = range_to_P.subset_by_type( MBVERTEX );
01288 
01289         for( Range::iterator it = V.begin(); it != V.end(); ++it )
01290         {
01291             EntityHandle v = *it;
01292             int index = lagr_verts.index( v );  // will be the same index as the corresponding vertex in euler verts
01293             assert( -1 != index );
01294             int n                = TLv.get_n();
01295             TLv.vi_wr[2 * n]     = to_proc;                // send to processor
01296             TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
01297             TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
01298             TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
01299             TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
01300             TLv.inc_n();
01301         }
01302         // also, prep the 2d cells for sending ...
01303         Range Q = range_to_P.subset_by_dimension( 2 );
01304         for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
01305         {
01306             EntityHandle q = *it;  // this is a src cell
01307             int global_id;
01308             rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
01309             int n                        = TLq.get_n();
01310             TLq.vi_wr[sizeTuple * n]     = to_proc;    //
01311             TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
01312             const EntityHandle* conn4;
01313             int num_nodes;
01314             rval = mb->get_connectivity(
01315                 q, conn4, num_nodes );  // could be up to 10;ERRORR( rval, "can't get connectivity for quad" );
01316             if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
01317             for( int i = 0; i < num_nodes; i++ )
01318             {
01319                 EntityHandle v = conn4[i];
01320                 int index      = lagr_verts.index( v );
01321                 assert( -1 != index );
01322                 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
01323             }
01324             for( int k = num_nodes; k < max_edges_1; k++ )
01325             {
01326                 TLq.vi_wr[sizeTuple * n + 2 + k] =
01327                     0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
01328             }
01329             EntityHandle tgtCell;
01330             rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" );
01331             TLq.vul_wr[n] = tgtCell;  // this will be sent to remote_cells, to be able to come back
01332             TLq.inc_n();
01333         }
01334     }
01335     // now we can route them to each processor
01336     // now we are done populating the tuples; route them to the appropriate processors
01337     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01338     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
01339     // the elements are already in localEnts;
01340 
01341     // maps from global ids to new vertex and quad handles, that are added
01342     std::map< int, EntityHandle > globalID_to_handle;
01343     // we already have vertices from lagr set; they are already in the processor, even before
01344     // receiving other verts from neighbors
01345     int k = 0;
01346     for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
01347     {
01348         globalID_to_handle[gids[k]] = *vit;  // a little bit of overkill
01349         // we do know that the global ids between euler and lagr verts are parallel
01350     }
01351     /*std::map<int, EntityHandle> globalID_to_eh;*/  // do we need this one?
01352     globalID_to_eh.clear();
01353     // now, look at every TLv, and see if we have to create a vertex there or not
01354     int n = TLv.get_n();  // the size of the points received
01355     for( int i = 0; i < n; i++ )
01356     {
01357         int globalId = TLv.vi_rd[2 * i + 1];
01358         if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01359         {
01360             EntityHandle new_vert;
01361             double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
01362             rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01363             globalID_to_handle[globalId] = new_vert;
01364         }
01365     }
01366 
01367     // now, all dep points should be at their place
01368     // look in the local list of 2d cells for this proc, and create all those cells if needed
01369     // it may be an overkill, but because it does not involve communication, we do it anyway
01370     Range& local  = Rto[my_rank];
01371     Range local_q = local.subset_by_dimension( 2 );
01372     // the local should have all the vertices in lagr_verts
01373     for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
01374     {
01375         EntityHandle q = *it;  // these are from lagr cells, local
01376         int gid_el;
01377         rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
01378         globalID_to_eh[gid_el] = q;  // do we need this? maybe to just mark the ones on this processor
01379         // maybe a range of global cell ids is fine?
01380     }
01381     // now look at all elements received through; we do not want to duplicate them
01382     n = TLq.get_n();  // number of elements received by this processor
01383     // a cell should be received from one proc only; so why are we so worried about duplicated
01384     // elements? a vertex can be received from multiple sources, that is fine
01385 
01386     remote_cells = new TupleList();
01387     remote_cells->initialize( 2, 0, 1, 0, n );  // no tracers anymore in these tuples
01388     remote_cells->enableWriteAccess();
01389     for( int i = 0; i < n; i++ )
01390     {
01391         int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
01392         int from_proc  = TLq.vi_rd[sizeTuple * i];
01393         // do we already have a quad with this global ID, represented?
01394         if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
01395         {
01396             // construct the conn quad
01397             EntityHandle new_conn[MAXEDGES];
01398             int nnodes = -1;
01399             for( int j = 0; j < max_edges_1; j++ )
01400             {
01401                 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
01402                 if( vgid == 0 )
01403                     new_conn[j] = 0;
01404                 else
01405                 {
01406                     assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
01407                     new_conn[j] = globalID_to_handle[vgid];
01408                     nnodes      = j + 1;  // nodes are at the beginning, and are variable number
01409                 }
01410             }
01411             EntityHandle new_element;
01412             //
01413             EntityType entType = MBQUAD;
01414             if( nnodes > 4 ) entType = MBPOLYGON;
01415             if( nnodes < 4 ) entType = MBTRI;
01416             rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
01417             globalID_to_eh[globalIdEl] = new_element;
01418             local_q.insert( new_element );
01419             rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " );
01420         }
01421         remote_cells->vi_wr[2 * i]     = from_proc;
01422         remote_cells->vi_wr[2 * i + 1] = globalIdEl;
01423         //    remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
01424         remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
01425         remote_cells->inc_n();
01426     }
01427     // now, create a new set, covering_set
01428     rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " );
01429     rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " );
01430     // order the remote cells tuple list, with the global id, because we will search in it
01431     // remote_cells->print("remote_cells before sorting");
01432     moab::TupleList::buffer sort_buffer;
01433     sort_buffer.buffer_init( n );
01434     remote_cells->sort( 1, &sort_buffer );
01435     sort_buffer.reset();
01436     return MB_SUCCESS;
01437     // end copy
01438 }
01439 
01440 ErrorCode Intx2Mesh::resolve_intersection_sharing()
01441 {
01442     if( parcomm && parcomm->size() > 1 )
01443     {
01444         /*
01445             moab::ParallelMergeMesh pm(parcomm, epsilon_1);
01446             ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip
01447            local merge, use dim 2 ERRORR(rval, "can't merge intersection ");
01448         */
01449         // look at non-owned shared vertices, that could be part of original source set
01450         // they should be removed from intx set reference, because they might not have a
01451         // correspondent on the other task
01452         Range nonOwnedVerts;
01453         Range vertsInIntx;
01454         Range intxCells;
01455         ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval );
01456         rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval );
01457 
01458         rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval );
01459 
01460         // some of these vertices can be in original set 1, which was covered, transported;
01461         // but they should not be "shared" from the intx point of view, because they are not shared
01462         // with another task they might have come from coverage as a plain vertex, so losing the
01463         // sharing property ?
01464 
01465         Range coverVerts;
01466         rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval );
01467         // find out those that are on the interface
01468         Range vertsCovInterface;
01469         rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval );
01470         // how many of these are in
01471         Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts );
01472         // first, get all cells connected to these vertices, from intxCells
01473 
01474         Range connectedCells;
01475         rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval );
01476         // only those in intx set:
01477         connectedCells = intersect( connectedCells, intxCells );
01478         // first duplicate vertices in question:
01479         std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
01480         for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
01481         {
01482             EntityHandle vertex = *vit;
01483             double coords[3];
01484             rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval );
01485             EntityHandle newVertex;
01486             rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval );
01487             duplicatedVerticesMap[vertex] = newVertex;
01488         }
01489 
01490         // look now at connectedCells, and change their connectivities:
01491         for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
01492         {
01493             EntityHandle intxCell = *eit;
01494             // replace connectivity
01495             std::vector< EntityHandle > connectivity;
01496             rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval );
01497             for( size_t i = 0; i < connectivity.size(); i++ )
01498             {
01499                 EntityHandle currentVertex                           = connectivity[i];
01500                 std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
01501                 if( mit != duplicatedVerticesMap.end() )
01502                 {
01503                     connectivity[i] = mit->second;  // replace connectivity directly
01504                 }
01505             }
01506             int nnodes = (int)connectivity.size();
01507             rval       = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval );
01508         }
01509     }
01510     return MB_SUCCESS;
01511 }
01512 #endif /* MOAB_HAVE_MPI */
01513 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines