MOAB: Mesh Oriented datABase  (version 5.2.1)
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 <float.h>
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         tolerance = 3 * tolerance;  // why ?
00343         if( !my_rank )
00344             std::cout << " max edge length: " << max_length << "  tolerance for kd tree: " << tolerance << "\n";
00345     }
00346     for( Range::iterator it = rs2.begin(); it != rs2.end(); ++it )
00347     {
00348         EntityHandle tcell = *it;
00349         // find vertex positions
00350         const EntityHandle* conn = NULL;
00351         int nnodes               = 0;
00352         rval                     = mb->get_connectivity( tcell, conn, nnodes );MB_CHK_ERR( rval );
00353         // find leaves close to those positions
00354         double areaTgtCell   = setup_tgt_cell( tcell, nnodes );  // this is the area in the gnomonic plane
00355         double recoveredArea = 0;
00356         std::vector< double > positions;
00357         positions.resize( nnodes * 3 );
00358         rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
00359 
00360         // distance to search will be based on average edge length
00361         double av_len = 0;
00362         for( int k = 0; k < nnodes; k++ )
00363         {
00364             int ik      = ( k + 1 ) % nnodes;
00365             double len1 = 0;
00366             for( int j = 0; j < 3; j++ )
00367             {
00368                 double len2 = positions[3 * k + j] - positions[3 * ik + j];
00369                 len1 += len2 * len2;
00370             }
00371             av_len += sqrt( len1 );
00372         }
00373         if( nnodes > 0 ) av_len /= nnodes;
00374         // find leaves within a distance from each vertex of target
00375         // in those leaves, collect all cells; we will try for an intx in there
00376         Range close_source_cells;
00377         std::vector< EntityHandle > leaves;
00378         for( int i = 0; i < nnodes; i++ )
00379         {
00380 
00381             leaves.clear();
00382             rval = kd.distance_search( &positions[3 * i], av_len, leaves, tolerance, epsilon_1 );MB_CHK_ERR( rval );
00383 
00384             for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
00385             {
00386                 Range tmp;
00387                 rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
00388 
00389                 close_source_cells.merge( tmp.begin(), tmp.end() );
00390             }
00391         }
00392 #ifdef VERBOSE
00393         if( close_source_cells.empty() )
00394         {
00395             std::cout << " there are no close source cells to target cell " << tcell << " id from handle "
00396                       << mb->id_from_handle( tcell ) << "\n";
00397         }
00398 #endif
00399         for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end(); ++it2 )
00400         {
00401             EntityHandle startSrc = *it2;
00402             double area           = 0;
00403             // if area is > 0 , we have intersections
00404             double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
00405             //
00406             int nP = 0;
00407             int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
00408             int nsTgt, nsSrc;
00409             rval = computeIntersectionBetweenTgtAndSrc( tcell, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
00410             if( area > 0 )
00411             {
00412                 if( nP > 1 )
00413                 {  // this will also construct triangles/polygons in the new mesh, if needed
00414                     rval = findNodes( tcell, nnodes, startSrc, nsSrc, P, nP );MB_CHK_ERR( rval );
00415                 }
00416                 recoveredArea += area;
00417             }
00418         }
00419         recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fract
00420     }
00421     // before cleaning up , we need to settle the position of the intersection points
00422     // on the boundary edges
00423     // this needs to be collective, so we should maybe wait something
00424 #ifdef MOAB_HAVE_MPI
00425     rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
00426 #endif
00427 
00428     this->clean();
00429     return MB_SUCCESS;
00430 }
00431 // main interface; this will do the advancing front trick
00432 // some are triangles, some are quads, some are polygons ...
00433 ErrorCode Intx2Mesh::intersect_meshes( EntityHandle mbset1, EntityHandle mbset2, EntityHandle& outputSet )
00434 {
00435     ErrorCode rval;
00436     mbs1   = mbset1;  // set 1 is departure, and it is completely covering the euler set on proc
00437     mbs2   = mbset2;
00438     outSet = outputSet;
00439 #ifdef VERBOSE
00440     std::stringstream ffs, fft;
00441     ffs << "source_rank0" << my_rank << ".vtk";
00442     rval = mb->write_mesh( ffs.str().c_str(), &mbset1, 1 );MB_CHK_ERR( rval );
00443     fft << "target_rank0" << my_rank << ".vtk";
00444     rval = mb->write_mesh( fft.str().c_str(), &mbset2, 1 );MB_CHK_ERR( rval );
00445 
00446 #endif
00447     // really, should be something from t1 and t2; src is 1 (lagrange), tgt is 2 (euler)
00448 
00449     EntityHandle startSrc = 0, startTgt = 0;
00450 
00451     rval = mb->get_entities_by_dimension( mbs1, 2, rs1 );MB_CHK_ERR( rval );
00452     rval = mb->get_entities_by_dimension( mbs2, 2, rs2 );MB_CHK_ERR( rval );
00453     // std::cout << "rs1.size() = " << rs1.size() << " and rs2.size() = "  << rs2.size() << "\n";
00454     // std::cout.flush();
00455 
00456     createTags();  // will also determine max_edges_1, max_edges_2 (for src and tgt meshes)
00457 
00458     Range rs22 = rs2;  // a copy of the initial range; we will remove from it elements as we
00459                        // advance ; rs2 is needed for marking the polygon to the tgt parent
00460 
00461     // create the local kdd tree with source elements; will use it to search
00462     // more efficiently for the seeds in advancing front;
00463     // some of the target cells will not be covered by source cells, and they need to be eliminated
00464     // early from contention
00465 
00466     // build a kd tree with the rs1 (source) cells
00467     AdaptiveKDTree kd( mb );
00468     EntityHandle tree_root = 0;
00469     rval                   = kd.build_tree( rs1, &tree_root );MB_CHK_ERR( rval );
00470 
00471     while( !rs22.empty() )
00472     {
00473 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
00474         if( rs22.size() < rs2.size() )
00475         {
00476             std::cout << " possible not connected arrival mesh; my_rank: " << my_rank << " counting: " << counting
00477                       << "\n";
00478             std::stringstream ffo;
00479             ffo << "file0" << counting << "rank0" << my_rank << ".vtk";
00480             rval = mb->write_mesh( ffo.str().c_str(), &outSet, 1 );MB_CHK_ERR( rval );
00481         }
00482 #endif
00483         bool seedFound = false;
00484         for( Range::iterator it = rs22.begin(); it != rs22.end(); ++it )
00485         {
00486             startTgt  = *it;
00487             int found = 0;
00488             // find vertex positions
00489             const EntityHandle* conn = NULL;
00490             int nnodes               = 0;
00491             rval                     = mb->get_connectivity( startTgt, conn, nnodes );MB_CHK_ERR( rval );
00492             // find leaves close to those positions
00493             std::vector< double > positions;
00494             positions.resize( nnodes * 3 );
00495             rval = mb->get_coords( conn, nnodes, &positions[0] );MB_CHK_ERR( rval );
00496             // find leaves within a distance from each vertex of target
00497             // in those leaves, collect all cells; we will try for an intx in there, instead of
00498             // looping over all rs1 cells, as before
00499             Range close_source_cells;
00500             std::vector< EntityHandle > leaves;
00501             for( int i = 0; i < nnodes; i++ )
00502             {
00503 
00504                 leaves.clear();
00505                 rval = kd.distance_search( &positions[3 * i], epsilon_1, leaves, epsilon_1, epsilon_1 );MB_CHK_ERR( rval );
00506 
00507                 for( std::vector< EntityHandle >::iterator j = leaves.begin(); j != leaves.end(); ++j )
00508                 {
00509                     Range tmp;
00510                     rval = mb->get_entities_by_dimension( *j, 2, tmp );MB_CHK_ERR( rval );
00511 
00512                     close_source_cells.merge( tmp.begin(), tmp.end() );
00513                 }
00514             }
00515 
00516             for( Range::iterator it2 = close_source_cells.begin(); it2 != close_source_cells.end() && !found; ++it2 )
00517             {
00518                 startSrc    = *it2;
00519                 double area = 0;
00520                 // if area is > 0 , we have intersections
00521                 double P[10 * MAXEDGES];  // max 8 intx points + 8 more in the polygon
00522                 //
00523                 int nP = 0;
00524                 int nb[MAXEDGES], nr[MAXEDGES];  // sides 3 or 4? also, check boxes first
00525                 int nsTgt, nsSrc;
00526                 rval =
00527                     computeIntersectionBetweenTgtAndSrc( startTgt, startSrc, P, nP, area, nb, nr, nsSrc, nsTgt, true );MB_CHK_ERR( rval );
00528                 if( area > 0 )
00529                 {
00530                     found     = 1;
00531                     seedFound = true;
00532                     break;  // found 2 elements that intersect; these will be the seeds
00533                 }
00534             }
00535             if( found )
00536                 break;
00537             else
00538             {
00539 #if defined( VERBOSE )
00540                 std::cout << " on rank " << my_rank << " target cell " << ID_FROM_HANDLE( startTgt )
00541                           << " not intx with any source\n";
00542 #endif
00543                 rs22.erase( startTgt );
00544             }
00545         }
00546         if( !seedFound ) continue;  // continue while(!rs22.empty())
00547 
00548         std::queue< EntityHandle > srcQueue;  // these are corresponding to Ta,
00549         srcQueue.push( startSrc );
00550         std::queue< EntityHandle > tgtQueue;
00551         tgtQueue.push( startTgt );
00552 
00553         Range toResetSrcs;  // will be used to reset src flags for every tgt quad
00554         // processed
00555 
00556         /*if (my_rank==0)
00557           dbg_1 = 1;*/
00558         unsigned char used = 1;
00559         // mark the start tgt quad as used, so it will not come back again
00560         rval = mb->tag_set_data( TgtFlagTag, &startTgt, 1, &used );MB_CHK_ERR( rval );
00561         while( !tgtQueue.empty() )
00562         {
00563             // flags for the side : 0 means a src cell not found on side
00564             // a paired src not found yet for the neighbors of tgt
00565             Range nextSrc[MAXEDGES];  // there are new ranges of possible next src cells for
00566                                       // seeding the side j of tgt cell
00567 
00568             EntityHandle currentTgt = tgtQueue.front();
00569             tgtQueue.pop();
00570             int nsidesTgt;                                                   // will be initialized now
00571             double areaTgtCell   = setup_tgt_cell( currentTgt, nsidesTgt );  // this is the area in the gnomonic plane
00572             double recoveredArea = 0;
00573             // get the neighbors of tgt, and if they are solved already, do not bother with that
00574             // side of tgt
00575             EntityHandle tgtNeighbors[MAXEDGES] = { 0 };
00576             rval                                = mb->tag_get_data( tgtNeighTag, &currentTgt, 1, tgtNeighbors );MB_CHK_SET_ERR( rval, "can't get neighbors of current tgt" );
00577 #ifdef ENABLE_DEBUG
00578             if( dbg_1 )
00579             {
00580                 std::cout << "Next: neighbors for current tgt ";
00581                 for( int kk = 0; kk < nsidesTgt; kk++ )
00582                 {
00583                     if( tgtNeighbors[kk] > 0 )
00584                         std::cout << mb->id_from_handle( tgtNeighbors[kk] ) << " ";
00585                     else
00586                         std::cout << 0 << " ";
00587                 }
00588                 std::cout << std::endl;
00589             }
00590 #endif
00591             // now get the status of neighbors; if already solved, make them 0, so not to bother
00592             // anymore on that side of tgt
00593             for( int j = 0; j < nsidesTgt; j++ )
00594             {
00595                 EntityHandle tgtNeigh = tgtNeighbors[j];
00596                 unsigned char status  = 1;
00597                 if( tgtNeigh == 0 ) continue;
00598                 rval = mb->tag_get_data( TgtFlagTag, &tgtNeigh, 1, &status );MB_CHK_ERR( rval );                     // status 0 is unused
00599                 if( 1 == status ) tgtNeighbors[j] = 0;  // so will not look anymore on this side of tgt
00600             }
00601 
00602 #ifdef ENABLE_DEBUG
00603             if( dbg_1 )
00604             {
00605                 std::cout << "reset sources: ";
00606                 for( Range::iterator itr = toResetSrcs.begin(); itr != toResetSrcs.end(); ++itr )
00607                     std::cout << mb->id_from_handle( *itr ) << " ";
00608                 std::cout << std::endl;
00609             }
00610 #endif
00611             EntityHandle currentSrc = srcQueue.front();
00612             // tgt and src queues are parallel; for clarity we should have kept in the queue pairs
00613             // of entity handle std::pair<EntityHandle, EntityHandle>; so just one queue, with
00614             // pairs;
00615             //  at every moment, the queue contains pairs of cells that intersect, and they form the
00616             //  "advancing front"
00617             srcQueue.pop();
00618             toResetSrcs.clear();  // empty the range of used srcs, will have to be set unused again,
00619             // at the end of tgt element processing
00620             toResetSrcs.insert( currentSrc );
00621             // mb2->set_tag_data
00622             std::queue< EntityHandle > localSrc;
00623             localSrc.push( currentSrc );
00624 #ifdef VERBOSE
00625             int countingStart = counting;
00626 #endif
00627             // will advance-front search in the neighborhood of tgt cell, until we finish processing
00628             // all
00629             //   possible src cells; localSrc queue will contain all possible src cells that cover
00630             //   the current tgt cell
00631             while( !localSrc.empty() )
00632             {
00633                 //
00634                 EntityHandle srcT = localSrc.front();
00635                 localSrc.pop();
00636                 double P[10 * MAXEDGES], area;  //
00637                 int nP           = 0;
00638                 int nb[MAXEDGES] = { 0 };
00639                 int nr[MAXEDGES] = { 0 };
00640 
00641                 int nsidesSrc;  ///
00642                 // area is in 2d, points are in 3d (on a sphere), back-projected, or in a plane
00643                 // intersection points could include the vertices of initial elements
00644                 // nb [j] = 0 means no intersection on the side j for element src (markers)
00645                 // nb [j] = 1 means that the side j (from j to j+1) of src poly intersects the
00646                 // tgt poly.  A potential next poly in the tgt queue is the tgt poly that is
00647                 // adjacent to this side
00648                 rval = computeIntersectionBetweenTgtAndSrc( /* tgt */ currentTgt, srcT, P, nP, area, nb, nr, nsidesSrc,
00649                                                             nsidesTgt );MB_CHK_ERR( rval );
00650                 if( nP > 0 )
00651                 {
00652 #ifdef ENABLE_DEBUG
00653                     if( dbg_1 )
00654                     {
00655                         for( int k = 0; k < 3; k++ )
00656                         {
00657                             std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
00658                         }
00659                     }
00660 #endif
00661 
00662                     // intersection found: output P and original triangles if nP > 2
00663                     EntityHandle neighbors[MAXEDGES] = { 0 };
00664                     rval                             = mb->tag_get_data( srcNeighTag, &srcT, 1, neighbors );
00665                     if( rval != MB_SUCCESS )
00666                     {
00667                         std::cout << " can't get the neighbors for src element " << mb->id_from_handle( srcT );
00668                         return MB_FAILURE;
00669                     }
00670 
00671                     // add neighbors to the localSrc queue, if they are not marked
00672                     for( int nn = 0; nn < nsidesSrc; nn++ )
00673                     {
00674                         EntityHandle neighbor = neighbors[nn];
00675                         if( neighbor > 0 && nb[nn] > 0 )  // advance across src boundary nn
00676                         {
00677                             if( toResetSrcs.find( neighbor ) == toResetSrcs.end() )
00678                             {
00679                                 localSrc.push( neighbor );
00680 #ifdef ENABLE_DEBUG
00681                                 if( dbg_1 )
00682                                 {
00683                                     std::cout << " local src elem " << mb->id_from_handle( neighbor )
00684                                               << " for tgt:" << mb->id_from_handle( currentTgt ) << "\n";
00685                                     mb->list_entities( &neighbor, 1 );
00686                                 }
00687 #endif
00688                                 toResetSrcs.insert( neighbor );
00689                             }
00690                         }
00691                     }
00692                     // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
00693                     for( int nn = 0; nn < nsidesTgt; nn++ )
00694                     {
00695                         if( nr[nn] > 0 && tgtNeighbors[nn] > 0 )
00696                             nextSrc[nn].insert( srcT );  // potential src cell that can intersect
00697                                                          // the tgt neighbor nn
00698                     }
00699                     if( nP > 1 )
00700                     {  // this will also construct triangles/polygons in the new mesh, if needed
00701                         rval = findNodes( currentTgt, nsidesTgt, srcT, nsidesSrc, P, nP );MB_CHK_ERR( rval );
00702                     }
00703 
00704                     recoveredArea += area;
00705                 }
00706 #ifdef ENABLE_DEBUG
00707                 else if( dbg_1 )
00708                 {
00709                     std::cout << " tgt, src, do not intersect: " << mb->id_from_handle( currentTgt ) << " "
00710                               << mb->id_from_handle( srcT ) << "\n";
00711                 }
00712 #endif
00713             }                                                               // end while (!localSrc.empty())
00714             recoveredArea = ( recoveredArea - areaTgtCell ) / areaTgtCell;  // replace now with recovery fraction
00715 #if defined( ENABLE_DEBUG ) || defined( VERBOSE )
00716             if( fabs( recoveredArea ) > epsilon_1 )
00717             {
00718 #ifdef VERBOSE
00719                 std::cout << " tgt area: " << areaTgtCell << " recovered :" << recoveredArea * ( 1 + areaTgtCell )
00720                           << " fraction error recovery:" << recoveredArea
00721                           << " tgtID: " << mb->id_from_handle( currentTgt ) << " countingStart:" << countingStart
00722                           << "\n";
00723 #endif
00724             }
00725 #endif
00726             // here, we are finished with tgtCurrent, take it out of the rs22 range (tgt, arrival
00727             // mesh)
00728             rs22.erase( currentTgt );
00729             // also, look at its neighbors, and add to the seeds a next one
00730 
00731             for( int j = 0; j < nsidesTgt; j++ )
00732             {
00733                 EntityHandle tgtNeigh = tgtNeighbors[j];
00734                 if( tgtNeigh == 0 || nextSrc[j].size() == 0 )  // if tgt is bigger than src, there could be no src
00735                                                                // to advance on that side
00736                     continue;
00737                 int nsidesTgt2 = 0;
00738                 setup_tgt_cell( tgtNeigh,
00739                                 nsidesTgt2 );  // find possible intersection with src cell from nextSrc
00740                 for( Range::iterator nit = nextSrc[j].begin(); nit != nextSrc[j].end(); ++nit )
00741                 {
00742                     EntityHandle nextB = *nit;
00743                     // we identified tgt quad n[j] as possibly intersecting with neighbor j of the
00744                     // src quad
00745                     double P[10 * MAXEDGES], area;  //
00746                     int nP           = 0;
00747                     int nb[MAXEDGES] = { 0 };
00748                     int nr[MAXEDGES] = { 0 };
00749 
00750                     int nsidesSrc;  ///
00751                     rval = computeIntersectionBetweenTgtAndSrc(
00752                         /* tgt */ tgtNeigh, nextB, P, nP, area, nb, nr, nsidesSrc, nsidesTgt2 );MB_CHK_ERR( rval );
00753                     if( area > 0 )
00754                     {
00755                         tgtQueue.push( tgtNeigh );
00756                         srcQueue.push( nextB );
00757 #ifdef ENABLE_DEBUG
00758                         if( dbg_1 )
00759                             std::cout << "new polys pushed: src, tgt:" << mb->id_from_handle( tgtNeigh ) << " "
00760                                       << mb->id_from_handle( nextB ) << std::endl;
00761 #endif
00762                         rval = mb->tag_set_data( TgtFlagTag, &tgtNeigh, 1, &used );MB_CHK_ERR( rval );
00763                         break;  // so we are done with this side of tgt, we have found a proper next
00764                                 // seed
00765                     }
00766                 }
00767             }
00768 
00769         }  // end while (!tgtQueue.empty())
00770     }
00771 #ifdef ENABLE_DEBUG
00772     if( dbg_1 )
00773     {
00774         for( int k = 0; k < 6; k++ )
00775             mout_1[k].close();
00776     }
00777 #endif
00778     // before cleaning up , we need to settle the position of the intersection points
00779     // on the boundary edges
00780     // this needs to be collective, so we should maybe wait something
00781 #ifdef MOAB_HAVE_MPI
00782     rval = resolve_intersection_sharing();MB_CHK_SET_ERR( rval, "can't correct position, Intx2Mesh.cpp \n" );
00783 #endif
00784 
00785     this->clean();
00786     return MB_SUCCESS;
00787 }
00788 
00789 // clean some memory allocated
00790 void Intx2Mesh::clean()
00791 {
00792     //
00793     int indx = 0;
00794     for( Range::iterator eit = TgtEdges.begin(); eit != TgtEdges.end(); ++eit, indx++ )
00795     {
00796         delete extraNodesVec[indx];
00797     }
00798     // extraNodesMap.clear();
00799     extraNodesVec.clear();
00800     // also, delete some bit tags, used to mark processed tgts and srcs
00801     mb->tag_delete( TgtFlagTag );
00802     counting = 0;  // reset counting to original value
00803 }
00804 // this method will reduce number of nodes, collapse edges that are of length 0
00805 // so a polygon like 428 431 431 will become a line 428 431
00806 // or something like 428 431 431 531 -> 428 431 531
00807 void Intx2Mesh::correct_polygon( EntityHandle* nodes, int& nP )
00808 {
00809     int i = 0;
00810     while( i < nP )
00811     {
00812         int nextIndex = ( i + 1 ) % nP;
00813         if( nodes[i] == nodes[nextIndex] )
00814         {
00815 #ifdef ENABLE_DEBUG
00816             // we need to reduce nP, and collapse nodes
00817             if( dbg_1 )
00818             {
00819                 std::cout << " nodes duplicated in list: ";
00820                 for( int j = 0; j < nP; j++ )
00821                     std::cout << nodes[j] << " ";
00822                 std::cout << "\n";
00823                 std::cout << " node " << nodes[i] << " at index " << i << " is duplicated"
00824                           << "\n";
00825             }
00826 #endif
00827             // this will work even if we start from 1 2 3 1; when i is 3, we find nextIndex is 0,
00828             // then next thing does nothing
00829             //  (nP-1 is 3, so k is already >= nP-1); it will result in nodes -> 1, 2, 3
00830             for( int k = i; k < nP - 1; k++ )
00831                 nodes[k] = nodes[k + 1];
00832             nP--;  // decrease the number of nodes; also, decrease i, just if we may need to check
00833                    // again
00834             i--;
00835         }
00836         i++;
00837     }
00838     return;
00839 }
00840 #ifdef MOAB_HAVE_MPI
00841 ErrorCode Intx2Mesh::build_processor_euler_boxes( EntityHandle euler_set, Range& local_verts )
00842 {
00843     localEnts.clear();
00844     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
00845 
00846     rval                = mb->get_connectivity( localEnts, local_verts );
00847     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
00848 
00849     assert( parcomm != NULL );
00850 
00851     // get the position of local vertices, and decide local boxes (allBoxes...)
00852     double bmin[3] = { DBL_MAX, DBL_MAX, DBL_MAX };
00853     double bmax[3] = { -DBL_MAX, -DBL_MAX, -DBL_MAX };
00854 
00855     std::vector< double > coords( 3 * num_local_verts );
00856     rval = mb->get_coords( local_verts, &coords[0] );ERRORR( rval, "can't get coords of vertices " );
00857 
00858     for( int i = 0; i < num_local_verts; i++ )
00859     {
00860         for( int k = 0; k < 3; k++ )
00861         {
00862             double val = coords[3 * i + k];
00863             if( val < bmin[k] ) bmin[k] = val;
00864             if( val > bmax[k] ) bmax[k] = val;
00865         }
00866     }
00867     int numprocs = parcomm->proc_config().proc_size();
00868     allBoxes.resize( 6 * numprocs );
00869 
00870     my_rank = parcomm->proc_config().proc_rank();
00871     for( int k = 0; k < 3; k++ )
00872     {
00873         allBoxes[6 * my_rank + k]     = bmin[k];
00874         allBoxes[6 * my_rank + 3 + k] = bmax[k];
00875     }
00876 
00877     // now communicate to get all boxes
00878     int mpi_err;
00879 #if( MPI_VERSION >= 2 )
00880     // use "in place" option
00881     mpi_err = MPI_Allgather( MPI_IN_PLACE, 0, MPI_DATATYPE_NULL, &allBoxes[0], 6, MPI_DOUBLE,
00882                              parcomm->proc_config().proc_comm() );
00883 #else
00884     {
00885         std::vector< double > allBoxes_tmp( 6 * parcomm->proc_config().proc_size() );
00886         mpi_err  = MPI_Allgather( &allBoxes[6 * my_rank], 6, MPI_DOUBLE, &allBoxes_tmp[0], 6, MPI_DOUBLE,
00887                                  parcomm->proc_config().proc_comm() );
00888         allBoxes = allBoxes_tmp;
00889     }
00890 #endif
00891     if( MPI_SUCCESS != mpi_err ) return MB_FAILURE;
00892 
00893     if( my_rank == 0 )
00894     {
00895         std::cout << " maximum number of vertices per cell are " << max_edges_1 << " on first mesh and " << max_edges_2
00896                   << " on second mesh \n";
00897         for( int i = 0; i < numprocs; i++ )
00898         {
00899             std::cout << "proc: " << i << " box min: " << allBoxes[6 * i] << " " << allBoxes[6 * i + 1] << " "
00900                       << allBoxes[6 * i + 2] << " \n";
00901             std::cout << "        box max: " << allBoxes[6 * i + 3] << " " << allBoxes[6 * i + 4] << " "
00902                       << allBoxes[6 * i + 5] << " \n";
00903         }
00904     }
00905 
00906     return MB_SUCCESS;
00907 }
00908 ErrorCode Intx2Mesh::create_departure_mesh_2nd_alg( EntityHandle& euler_set, EntityHandle& covering_lagr_set )
00909 {
00910     // compute the bounding box on each proc
00911     assert( parcomm != NULL );
00912 
00913     localEnts.clear();
00914     ErrorCode rval = mb->get_entities_by_dimension( euler_set, 2, localEnts );ERRORR( rval, "can't get ents by dimension" );
00915 
00916     Tag dpTag = 0;
00917     std::string tag_name( "DP" );
00918     rval = mb->tag_get_handle( tag_name.c_str(), 3, MB_TYPE_DOUBLE, dpTag, MB_TAG_DENSE );ERRORR( rval, "can't get DP tag" );
00919 
00920     EntityHandle dum = 0;
00921     Tag corrTag;
00922     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" );
00923     // get all local verts
00924     Range local_verts;
00925     rval                = mb->get_connectivity( localEnts, local_verts );
00926     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
00927 
00928     rval = Intx2Mesh::build_processor_euler_boxes( euler_set, local_verts );ERRORR( rval, "can't build processor boxes" );
00929 
00930     std::vector< int > gids( num_local_verts );
00931     rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
00932 
00933     // now see the departure points; to what boxes should we send them?
00934     std::vector< double > dep_points( 3 * num_local_verts );
00935     rval = mb->tag_get_data( dpTag, local_verts, (void*)&dep_points[0] );ERRORR( rval, "can't get DP tag values" );
00936     // ranges to send to each processor; will hold vertices and elements (quads?)
00937     // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
00938     std::map< int, Range > Rto;
00939     int numprocs = parcomm->proc_config().proc_size();
00940 
00941     for( Range::iterator eit = localEnts.begin(); eit != localEnts.end(); ++eit )
00942     {
00943         EntityHandle q = *eit;
00944         const EntityHandle* conn4;
00945         int num_nodes;
00946         rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
00947         CartVect qbmin( DBL_MAX );
00948         CartVect qbmax( -DBL_MAX );
00949         for( int i = 0; i < num_nodes; i++ )
00950         {
00951             EntityHandle v = conn4[i];
00952             size_t index   = local_verts.find( v ) - local_verts.begin();
00953             CartVect dp( &dep_points[3 * index] );  // will use constructor
00954             for( int j = 0; j < 3; j++ )
00955             {
00956                 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
00957                 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
00958             }
00959         }
00960         for( int p = 0; p < numprocs; p++ )
00961         {
00962             CartVect bbmin( &allBoxes[6 * p] );
00963             CartVect bbmax( &allBoxes[6 * p + 3] );
00964             if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); }
00965         }
00966     }
00967 
00968     // now, build TLv and TLq, for each p
00969     size_t numq = 0;
00970     size_t numv = 0;
00971     for( int p = 0; p < numprocs; p++ )
00972     {
00973         if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
00974         Range& range_to_P = Rto[p];
00975         // add the vertices to it
00976         if( range_to_P.empty() ) continue;  // nothing to send to proc p
00977         Range vertsToP;
00978         rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
00979         numq = numq + range_to_P.size();
00980         numv = numv + vertsToP.size();
00981         range_to_P.merge( vertsToP );
00982     }
00983     TupleList TLv;
00984     TupleList TLq;
00985     TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
00986     TLv.enableWriteAccess();
00987 
00988     int sizeTuple = 2 + max_edges_1;  // determined earlier, for src, first mesh
00989     TLq.initialize( 2 + max_edges_1, 0, 1, 0,
00990                     numq );  // to proc, elem GLOBAL ID, connectivity[10] (global ID v), local eh
00991     TLq.enableWriteAccess();
00992 #ifdef VERBOSE
00993     std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
00994 #endif
00995     for( int to_proc = 0; to_proc < numprocs; to_proc++ )
00996     {
00997         if( to_proc == (int)my_rank ) continue;
00998         Range& range_to_P = Rto[to_proc];
00999         Range V           = range_to_P.subset_by_type( MBVERTEX );
01000 
01001         for( Range::iterator it = V.begin(); it != V.end(); ++it )
01002         {
01003             EntityHandle v       = *it;
01004             unsigned int index   = local_verts.find( v ) - local_verts.begin();
01005             int n                = TLv.get_n();
01006             TLv.vi_wr[2 * n]     = to_proc;                // send to processor
01007             TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
01008             TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
01009             TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
01010             TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
01011             TLv.inc_n();
01012         }
01013         // also, prep the quad for sending ...
01014         Range Q = range_to_P.subset_by_dimension( 2 );
01015         for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
01016         {
01017             EntityHandle q = *it;
01018             int global_id;
01019             rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
01020             int n                        = TLq.get_n();
01021             TLq.vi_wr[sizeTuple * n]     = to_proc;    //
01022             TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
01023             const EntityHandle* conn4;
01024             int num_nodes;
01025             rval = mb->get_connectivity( q, conn4,
01026                                          num_nodes );  // could be up to MAXEDGES, but it is limited by max_edges_1
01027             ERRORR( rval, "can't get connectivity for cell" );
01028             if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
01029             for( int i = 0; i < num_nodes; i++ )
01030             {
01031                 EntityHandle v                   = conn4[i];
01032                 unsigned int index               = local_verts.find( v ) - local_verts.begin();
01033                 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
01034             }
01035             for( int k = num_nodes; k < max_edges_1; k++ )
01036             {
01037                 TLq.vi_wr[sizeTuple * n + 2 + k] =
01038                     0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
01039             }
01040             TLq.vul_wr[n] = q;  // save here the entity handle, it will be communicated back
01041             // maybe we should forget about global ID
01042             TLq.inc_n();
01043         }
01044     }
01045 
01046     // now we are done populating the tuples; route them to the appropriate processors
01047     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01048     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
01049     // the elements are already in localEnts;
01050 
01051     // maps from global ids to new vertex and quad handles, that are added
01052     std::map< int, EntityHandle > globalID_to_handle;
01053     /*std::map<int, EntityHandle> globalID_to_eh;*/
01054     globalID_to_eh.clear();  // need for next iteration
01055     // now, look at every TLv, and see if we have to create a vertex there or not
01056     int n = TLv.get_n();  // the size of the points received
01057     for( int i = 0; i < n; i++ )
01058     {
01059         int globalId = TLv.vi_rd[2 * i + 1];
01060         if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01061         {
01062             EntityHandle new_vert;
01063             double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
01064             rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01065             globalID_to_handle[globalId] = new_vert;
01066         }
01067     }
01068 
01069     // now, all dep points should be at their place
01070     // look in the local list of q for this proc, and create all those quads and vertices if needed
01071     // it may be an overkill, but because it does not involve communication, we do it anyway
01072     Range& local  = Rto[my_rank];
01073     Range local_q = local.subset_by_dimension( 2 );
01074     // the local should have all the vertices in local_verts
01075     for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
01076     {
01077         EntityHandle q = *it;
01078         int nnodes;
01079         const EntityHandle* conn4;
01080         rval = mb->get_connectivity( q, conn4, nnodes );ERRORR( rval, "can't get connectivity of local q " );
01081         EntityHandle new_conn[MAXEDGES];
01082         for( int i = 0; i < nnodes; i++ )
01083         {
01084             EntityHandle v1    = conn4[i];
01085             unsigned int index = local_verts.find( v1 ) - local_verts.begin();
01086             int globalId       = gids[index];
01087             if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01088             {
01089                 // we need to create that vertex, at this position dep_points
01090                 double dp_pos[3] = { dep_points[3 * index], dep_points[3 * index + 1], dep_points[3 * index + 2] };
01091                 EntityHandle new_vert;
01092                 rval = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01093                 globalID_to_handle[globalId] = new_vert;
01094             }
01095             new_conn[i] = globalID_to_handle[gids[index]];
01096         }
01097         EntityHandle new_element;
01098         //
01099         EntityType entType = MBQUAD;
01100         if( nnodes > 4 ) entType = MBPOLYGON;
01101         if( nnodes < 4 ) entType = MBTRI;
01102 
01103         rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new quad " );
01104         rval = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
01105         int gid_el;
01106         // get the global ID of the initial quad
01107         rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
01108         globalID_to_eh[gid_el] = new_element;
01109         // is this redundant or not?
01110         rval = mb->tag_set_data( corrTag, &new_element, 1, &q );ERRORR( rval, "can't set corr tag on new el" );
01111         // set the global id on new elem
01112         rval = mb->tag_set_data( gid, &new_element, 1, &gid_el );ERRORR( rval, "can't set global id tag on new el" );
01113     }
01114     // now look at all elements received through; we do not want to duplicate them
01115     n = TLq.get_n();  // number of elements received by this processor
01116     // form the remote cells, that will be used to send the tracer info back to the originating proc
01117     remote_cells = new TupleList();
01118     remote_cells->initialize( 2, 0, 1, 0, n );  // will not have tracer data anymore
01119     remote_cells->enableWriteAccess();
01120     for( int i = 0; i < n; i++ )
01121     {
01122         int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
01123         int from_proc  = TLq.vi_wr[sizeTuple * i];
01124         // do we already have a quad with this global ID, represented?
01125         if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
01126         {
01127             // construct the conn quad
01128             EntityHandle new_conn[MAXEDGES];
01129             int nnodes = -1;
01130             for( int j = 0; j < max_edges_1; j++ )
01131             {
01132                 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
01133                 if( vgid == 0 )
01134                     new_conn[j] = 0;
01135                 else
01136                 {
01137                     assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
01138                     new_conn[j] = globalID_to_handle[vgid];
01139                     nnodes      = j + 1;  // nodes are at the beginning, and are variable number
01140                 }
01141             }
01142             EntityHandle new_element;
01143             //
01144             EntityType entType = MBQUAD;
01145             if( nnodes > 4 ) entType = MBPOLYGON;
01146             if( nnodes < 4 ) entType = MBTRI;
01147             rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
01148             globalID_to_eh[globalIdEl] = new_element;
01149             rval                       = mb->add_entities( covering_lagr_set, &new_element, 1 );ERRORR( rval, "can't add new element to dep set" );
01150             /* rval = mb->tag_set_data(corrTag, &new_element, 1, &q);ERRORR(rval, "can't set corr tag on new el");*/
01151             remote_cells->vi_wr[2 * i]     = from_proc;
01152             remote_cells->vi_wr[2 * i + 1] = globalIdEl;
01153             //     remote_cells->vr_wr[i] = 0.; // no contribution yet sent back
01154             remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
01155             remote_cells->inc_n();
01156             // set the global id on new elem
01157             rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set global id tag on new el" );
01158         }
01159     }
01160     // order the remote cells tuple list, with the global id, because we will search in it
01161     // remote_cells->print("remote_cells before sorting");
01162     moab::TupleList::buffer sort_buffer;
01163     sort_buffer.buffer_init( n );
01164     remote_cells->sort( 1, &sort_buffer );
01165     sort_buffer.reset();
01166     return MB_SUCCESS;
01167 }
01168 
01169 // this algorithm assumes lagr set is already created, and some elements will be coming from
01170 // other procs, and populate the covering_set
01171 // we need to keep in a tuple list the remote cells from other procs, because we need to send back
01172 // the intersection info (like area of the intx polygon, and the current concentration) maybe total
01173 // mass in that intx
01174 ErrorCode Intx2Mesh::create_departure_mesh_3rd_alg( EntityHandle& lagr_set, EntityHandle& covering_set )
01175 {
01176     EntityHandle dum = 0;
01177 
01178     Tag corrTag;
01179     ErrorCode rval = mb->tag_get_handle( CORRTAGNAME, 1, MB_TYPE_HANDLE, corrTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
01180     // start copy from 2nd alg
01181     // compute the bounding box on each proc
01182     assert( parcomm != NULL );
01183     if( 1 == parcomm->proc_config().proc_size() )
01184     {
01185         covering_set = lagr_set;  // nothing to communicate, it must be serial
01186         return MB_SUCCESS;
01187     }
01188 
01189     // get all local verts
01190     Range local_verts;
01191     rval                = mb->get_connectivity( localEnts, local_verts );
01192     int num_local_verts = (int)local_verts.size();ERRORR( rval, "can't get local vertices" );
01193 
01194     std::vector< int > gids( num_local_verts );
01195     rval = mb->tag_get_data( gid, local_verts, &gids[0] );ERRORR( rval, "can't get local vertices gids" );
01196 
01197     Range localDepCells;
01198     rval = mb->get_entities_by_dimension( lagr_set, 2, localDepCells );ERRORR( rval, "can't get ents by dimension from lagr set" );
01199 
01200     // get all lagr verts (departure vertices)
01201     Range lagr_verts;
01202     rval = mb->get_connectivity( localDepCells, lagr_verts );  // they should be created in
01203     // the same order as the euler vertices
01204     int num_lagr_verts = (int)lagr_verts.size();ERRORR( rval, "can't get local lagr vertices" );
01205 
01206     // now see the departure points position; to what boxes should we send them?
01207     std::vector< double > dep_points( 3 * num_lagr_verts );
01208     rval = mb->get_coords( lagr_verts, &dep_points[0] );ERRORR( rval, "can't get departure points position" );
01209     // ranges to send to each processor; will hold vertices and elements (quads?)
01210     // will look if the box of the dep quad covers box of euler mesh on proc (with tolerances)
01211     std::map< int, Range > Rto;
01212     int numprocs = parcomm->proc_config().proc_size();
01213 
01214     for( Range::iterator eit = localDepCells.begin(); eit != localDepCells.end(); ++eit )
01215     {
01216         EntityHandle q = *eit;
01217         const EntityHandle* conn4;
01218         int num_nodes;
01219         rval = mb->get_connectivity( q, conn4, num_nodes );ERRORR( rval, "can't get DP tag values" );
01220         CartVect qbmin( DBL_MAX );
01221         CartVect qbmax( -DBL_MAX );
01222         for( int i = 0; i < num_nodes; i++ )
01223         {
01224             EntityHandle v = conn4[i];
01225             int index      = lagr_verts.index( v );
01226             assert( -1 != index );
01227             CartVect dp( &dep_points[3 * index] );  // will use constructor
01228             for( int j = 0; j < 3; j++ )
01229             {
01230                 if( qbmin[j] > dp[j] ) qbmin[j] = dp[j];
01231                 if( qbmax[j] < dp[j] ) qbmax[j] = dp[j];
01232             }
01233         }
01234         for( int p = 0; p < numprocs; p++ )
01235         {
01236             CartVect bbmin( &allBoxes[6 * p] );
01237             CartVect bbmax( &allBoxes[6 * p + 3] );
01238             if( GeomUtil::boxes_overlap( bbmin, bbmax, qbmin, qbmax, box_error ) ) { Rto[p].insert( q ); }
01239         }
01240     }
01241 
01242     // now, build TLv and TLq, for each p
01243     size_t numq = 0;
01244     size_t numv = 0;
01245     for( int p = 0; p < numprocs; p++ )
01246     {
01247         if( p == (int)my_rank ) continue;  // do not "send" it, because it is already here
01248         Range& range_to_P = Rto[p];
01249         // add the vertices to it
01250         if( range_to_P.empty() ) continue;  // nothing to send to proc p
01251         Range vertsToP;
01252         rval = mb->get_connectivity( range_to_P, vertsToP );ERRORR( rval, "can't get connectivity" );
01253         numq = numq + range_to_P.size();
01254         numv = numv + vertsToP.size();
01255         range_to_P.merge( vertsToP );
01256     }
01257     TupleList TLv;
01258     TupleList TLq;
01259     TLv.initialize( 2, 0, 0, 3, numv );  // to proc, GLOBAL ID, DP points
01260     TLv.enableWriteAccess();
01261 
01262     int sizeTuple = 2 + max_edges_1;  // max edges could be up to MAXEDGES :) for polygons
01263     TLq.initialize( 2 + max_edges_1, 0, 1, 0,
01264                     numq );  // to proc, elem GLOBAL ID, connectivity[max_edges] (global ID v)
01265     // send also the corresponding tgt cell it will come to
01266     TLq.enableWriteAccess();
01267 #ifdef VERBOSE
01268     std::cout << "from proc " << my_rank << " send " << numv << " vertices and " << numq << " elements\n";
01269 #endif
01270 
01271     for( int to_proc = 0; to_proc < numprocs; to_proc++ )
01272     {
01273         if( to_proc == (int)my_rank ) continue;
01274         Range& range_to_P = Rto[to_proc];
01275         Range V           = range_to_P.subset_by_type( MBVERTEX );
01276 
01277         for( Range::iterator it = V.begin(); it != V.end(); ++it )
01278         {
01279             EntityHandle v = *it;
01280             int index = lagr_verts.index( v );  // will be the same index as the corresponding vertex in euler verts
01281             assert( -1 != index );
01282             int n                = TLv.get_n();
01283             TLv.vi_wr[2 * n]     = to_proc;                // send to processor
01284             TLv.vi_wr[2 * n + 1] = gids[index];            // global id needs index in the local_verts range
01285             TLv.vr_wr[3 * n]     = dep_points[3 * index];  // departure position, of the node local_verts[i]
01286             TLv.vr_wr[3 * n + 1] = dep_points[3 * index + 1];
01287             TLv.vr_wr[3 * n + 2] = dep_points[3 * index + 2];
01288             TLv.inc_n();
01289         }
01290         // also, prep the 2d cells for sending ...
01291         Range Q = range_to_P.subset_by_dimension( 2 );
01292         for( Range::iterator it = Q.begin(); it != Q.end(); ++it )
01293         {
01294             EntityHandle q = *it;  // this is a src cell
01295             int global_id;
01296             rval = mb->tag_get_data( gid, &q, 1, &global_id );ERRORR( rval, "can't get gid for polygon" );
01297             int n                        = TLq.get_n();
01298             TLq.vi_wr[sizeTuple * n]     = to_proc;    //
01299             TLq.vi_wr[sizeTuple * n + 1] = global_id;  // global id of element, used to identify it ...
01300             const EntityHandle* conn4;
01301             int num_nodes;
01302             rval = mb->get_connectivity(
01303                 q, conn4, num_nodes );  // could be up to 10;ERRORR( rval, "can't get connectivity for quad" );
01304             if( num_nodes > MAXEDGES ) ERRORR( MB_FAILURE, "too many nodes in a polygon" );
01305             for( int i = 0; i < num_nodes; i++ )
01306             {
01307                 EntityHandle v = conn4[i];
01308                 int index      = lagr_verts.index( v );
01309                 assert( -1 != index );
01310                 TLq.vi_wr[sizeTuple * n + 2 + i] = gids[index];
01311             }
01312             for( int k = num_nodes; k < max_edges_1; k++ )
01313             {
01314                 TLq.vi_wr[sizeTuple * n + 2 + k] =
01315                     0;  // fill the rest of node ids with 0; we know that the node ids start from 1!
01316             }
01317             EntityHandle tgtCell;
01318             rval = mb->tag_get_data( corrTag, &q, 1, &tgtCell );ERRORR( rval, "can't get corresponding tgt cell for dep cell" );
01319             TLq.vul_wr[n] = tgtCell;  // this will be sent to remote_cells, to be able to come back
01320             TLq.inc_n();
01321         }
01322     }
01323     // now we can route them to each processor
01324     // now we are done populating the tuples; route them to the appropriate processors
01325     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01326     ( parcomm->proc_config().crystal_router() )->gs_transfer( 1, TLq, 0 );
01327     // the elements are already in localEnts;
01328 
01329     // maps from global ids to new vertex and quad handles, that are added
01330     std::map< int, EntityHandle > globalID_to_handle;
01331     // we already have vertices from lagr set; they are already in the processor, even before
01332     // receiving other verts from neighbors
01333     int k = 0;
01334     for( Range::iterator vit = lagr_verts.begin(); vit != lagr_verts.end(); ++vit, k++ )
01335     {
01336         globalID_to_handle[gids[k]] = *vit;  // a little bit of overkill
01337         // we do know that the global ids between euler and lagr verts are parallel
01338     }
01339     /*std::map<int, EntityHandle> globalID_to_eh;*/  // do we need this one?
01340     globalID_to_eh.clear();
01341     // now, look at every TLv, and see if we have to create a vertex there or not
01342     int n = TLv.get_n();  // the size of the points received
01343     for( int i = 0; i < n; i++ )
01344     {
01345         int globalId = TLv.vi_rd[2 * i + 1];
01346         if( globalID_to_handle.find( globalId ) == globalID_to_handle.end() )
01347         {
01348             EntityHandle new_vert;
01349             double dp_pos[3] = { TLv.vr_wr[3 * i], TLv.vr_wr[3 * i + 1], TLv.vr_wr[3 * i + 2] };
01350             rval             = mb->create_vertex( dp_pos, new_vert );ERRORR( rval, "can't create new vertex " );
01351             globalID_to_handle[globalId] = new_vert;
01352         }
01353     }
01354 
01355     // now, all dep points should be at their place
01356     // look in the local list of 2d cells for this proc, and create all those cells if needed
01357     // it may be an overkill, but because it does not involve communication, we do it anyway
01358     Range& local  = Rto[my_rank];
01359     Range local_q = local.subset_by_dimension( 2 );
01360     // the local should have all the vertices in lagr_verts
01361     for( Range::iterator it = local_q.begin(); it != local_q.end(); ++it )
01362     {
01363         EntityHandle q = *it;  // these are from lagr cells, local
01364         int gid_el;
01365         rval = mb->tag_get_data( gid, &q, 1, &gid_el );ERRORR( rval, "can't get element global ID " );
01366         globalID_to_eh[gid_el] = q;  // do we need this? maybe to just mark the ones on this processor
01367         // maybe a range of global cell ids is fine?
01368     }
01369     // now look at all elements received through; we do not want to duplicate them
01370     n = TLq.get_n();  // number of elements received by this processor
01371     // a cell should be received from one proc only; so why are we so worried about duplicated
01372     // elements? a vertex can be received from multiple sources, that is fine
01373 
01374     remote_cells = new TupleList();
01375     remote_cells->initialize( 2, 0, 1, 0, n );  // no tracers anymore in these tuples
01376     remote_cells->enableWriteAccess();
01377     for( int i = 0; i < n; i++ )
01378     {
01379         int globalIdEl = TLq.vi_rd[sizeTuple * i + 1];
01380         int from_proc  = TLq.vi_rd[sizeTuple * i];
01381         // do we already have a quad with this global ID, represented?
01382         if( globalID_to_eh.find( globalIdEl ) == globalID_to_eh.end() )
01383         {
01384             // construct the conn quad
01385             EntityHandle new_conn[MAXEDGES];
01386             int nnodes = -1;
01387             for( int j = 0; j < max_edges_1; j++ )
01388             {
01389                 int vgid = TLq.vi_rd[sizeTuple * i + 2 + j];  // vertex global ID
01390                 if( vgid == 0 )
01391                     new_conn[j] = 0;
01392                 else
01393                 {
01394                     assert( globalID_to_handle.find( vgid ) != globalID_to_handle.end() );
01395                     new_conn[j] = globalID_to_handle[vgid];
01396                     nnodes      = j + 1;  // nodes are at the beginning, and are variable number
01397                 }
01398             }
01399             EntityHandle new_element;
01400             //
01401             EntityType entType = MBQUAD;
01402             if( nnodes > 4 ) entType = MBPOLYGON;
01403             if( nnodes < 4 ) entType = MBTRI;
01404             rval = mb->create_element( entType, new_conn, nnodes, new_element );ERRORR( rval, "can't create new element " );
01405             globalID_to_eh[globalIdEl] = new_element;
01406             local_q.insert( new_element );
01407             rval = mb->tag_set_data( gid, &new_element, 1, &globalIdEl );ERRORR( rval, "can't set gid on new element " );
01408         }
01409         remote_cells->vi_wr[2 * i]     = from_proc;
01410         remote_cells->vi_wr[2 * i + 1] = globalIdEl;
01411         //    remote_cells->vr_wr[i] = 0.; will have a different tuple for communication
01412         remote_cells->vul_wr[i] = TLq.vul_rd[i];  // this is the corresponding tgt cell (arrival)
01413         remote_cells->inc_n();
01414     }
01415     // now, create a new set, covering_set
01416     rval = mb->create_meshset( MESHSET_SET, covering_set );ERRORR( rval, "can't create new mesh set " );
01417     rval = mb->add_entities( covering_set, local_q );ERRORR( rval, "can't add entities to new mesh set " );
01418     // order the remote cells tuple list, with the global id, because we will search in it
01419     // remote_cells->print("remote_cells before sorting");
01420     moab::TupleList::buffer sort_buffer;
01421     sort_buffer.buffer_init( n );
01422     remote_cells->sort( 1, &sort_buffer );
01423     sort_buffer.reset();
01424     return MB_SUCCESS;
01425     // end copy
01426 }
01427 
01428 ErrorCode Intx2Mesh::resolve_intersection_sharing()
01429 {
01430     if( parcomm && parcomm->size() > 1 )
01431     {
01432         /*
01433             moab::ParallelMergeMesh pm(parcomm, epsilon_1);
01434             ErrorCode rval = pm.merge(outSet, false, 2); // resolve only the output set, do not skip
01435            local merge, use dim 2 ERRORR(rval, "can't merge intersection ");
01436         */
01437         // look at non-owned shared vertices, that could be part of original source set
01438         // they should be removed from intx set reference, because they might not have a
01439         // correspondent on the other task
01440         Range nonOwnedVerts;
01441         Range vertsInIntx;
01442         Range intxCells;
01443         ErrorCode rval = mb->get_entities_by_dimension( outSet, 2, intxCells );MB_CHK_ERR( rval );
01444         rval = mb->get_connectivity( intxCells, vertsInIntx );MB_CHK_ERR( rval );
01445 
01446         rval = parcomm->filter_pstatus( vertsInIntx, PSTATUS_NOT_OWNED, PSTATUS_AND, -1, &nonOwnedVerts );MB_CHK_ERR( rval );
01447 
01448         // some of these vertices can be in original set 1, which was covered, transported;
01449         // but they should not be "shared" from the intx point of view, because they are not shared
01450         // with another task they might have come from coverage as a plain vertex, so losing the
01451         // sharing property ?
01452 
01453         Range coverVerts;
01454         rval = mb->get_connectivity( rs1, coverVerts );MB_CHK_ERR( rval );
01455         // find out those that are on the interface
01456         Range vertsCovInterface;
01457         rval = parcomm->filter_pstatus( coverVerts, PSTATUS_INTERFACE, PSTATUS_AND, -1, &vertsCovInterface );MB_CHK_ERR( rval );
01458         // how many of these are in
01459         Range nodesToDuplicate = intersect( vertsCovInterface, nonOwnedVerts );
01460         // first, get all cells connected to these vertices, from intxCells
01461 
01462         Range connectedCells;
01463         rval = mb->get_adjacencies( nodesToDuplicate, 2, false, connectedCells, Interface::UNION );MB_CHK_ERR( rval );
01464         // only those in intx set:
01465         connectedCells = intersect( connectedCells, intxCells );
01466         // first duplicate vertices in question:
01467         std::map< EntityHandle, EntityHandle > duplicatedVerticesMap;
01468         for( Range::iterator vit = nodesToDuplicate.begin(); vit != nodesToDuplicate.end(); vit++ )
01469         {
01470             EntityHandle vertex = *vit;
01471             double coords[3];
01472             rval = mb->get_coords( &vertex, 1, coords );MB_CHK_ERR( rval );
01473             EntityHandle newVertex;
01474             rval = mb->create_vertex( coords, newVertex );MB_CHK_ERR( rval );
01475             duplicatedVerticesMap[vertex] = newVertex;
01476         }
01477 
01478         // look now at connectedCells, and change their connectivities:
01479         for( Range::iterator eit = connectedCells.begin(); eit != connectedCells.end(); eit++ )
01480         {
01481             EntityHandle intxCell = *eit;
01482             // replace connectivity
01483             std::vector< EntityHandle > connectivity;
01484             rval = mb->get_connectivity( &intxCell, 1, connectivity );MB_CHK_ERR( rval );
01485             for( size_t i = 0; i < connectivity.size(); i++ )
01486             {
01487                 EntityHandle currentVertex                           = connectivity[i];
01488                 std::map< EntityHandle, EntityHandle >::iterator mit = duplicatedVerticesMap.find( currentVertex );
01489                 if( mit != duplicatedVerticesMap.end() )
01490                 {
01491                     connectivity[i] = mit->second;  // replace connectivity directly
01492                 }
01493             }
01494             int nnodes = (int)connectivity.size();
01495             rval       = mb->set_connectivity( intxCell, &connectivity[0], nnodes );MB_CHK_ERR( rval );
01496         }
01497     }
01498     return MB_SUCCESS;
01499 }
01500 #endif /* MOAB_HAVE_MPI */
01501 } /* namespace moab */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines