Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
TempestRemapper.cpp
Go to the documentation of this file.
00001 /*
00002  * =====================================================================================
00003  *
00004  *       Filename:  TempestRemapper.hpp
00005  *
00006  *    Description:  Interface to the TempestRemap library to enable intersection and
00007  *                  high-order conservative remapping of climate solution from
00008  *                  arbitrary resolution of source and target grids on the sphere.
00009  *
00010  *         Author:  Vijay S. Mahadevan (vijaysm), [email protected]
00011  *
00012  * =====================================================================================
00013  */
00014 
00015 #include <string>
00016 #include <iostream>
00017 #include <cassert>
00018 #include "DebugOutput.hpp"
00019 #include "moab/Remapping/TempestRemapper.hpp"
00020 #include "moab/ReadUtilIface.hpp"
00021 
00022 // Intersection includes
00023 #include "moab/IntxMesh/Intx2MeshOnSphere.hpp"
00024 #include "moab/IntxMesh/IntxUtils.hpp"
00025 
00026 #include "moab/AdaptiveKDTree.hpp"
00027 #include "moab/SpatialLocator.hpp"
00028 
00029 // skinner for augmenting overlap mesh to complete coverage
00030 #include "moab/Skinner.hpp"
00031 #include "MBParallelConventions.h"
00032 
00033 #ifdef MOAB_HAVE_TEMPESTREMAP
00034 #include "GaussLobattoQuadrature.h"
00035 #endif
00036 
00037 // #define VERBOSE
00038 
00039 namespace moab
00040 {
00041 
00042 ///////////////////////////////////////////////////////////////////////////////////
00043 
00044 ErrorCode TempestRemapper::initialize( bool initialize_fsets )
00045 {
00046     ErrorCode rval;
00047     if( initialize_fsets )
00048     {
00049         rval = m_interface->create_meshset( moab::MESHSET_SET, m_source_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00050         rval = m_interface->create_meshset( moab::MESHSET_SET, m_target_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00051         rval = m_interface->create_meshset( moab::MESHSET_SET, m_overlap_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
00052     }
00053     else
00054     {
00055         m_source_set  = 0;
00056         m_target_set  = 0;
00057         m_overlap_set = 0;
00058     }
00059 
00060     is_parallel = false;
00061     is_root     = true;
00062     rank        = 0;
00063     size        = 1;
00064 #ifdef MOAB_HAVE_MPI
00065     int flagInit;
00066     MPI_Initialized( &flagInit );
00067     if( flagInit )
00068     {
00069         is_parallel = true;
00070         assert( m_pcomm != NULL );
00071         rank    = m_pcomm->rank();
00072         size    = m_pcomm->size();
00073         is_root = ( rank == 0 );
00074     }
00075 #endif
00076 
00077     m_source          = NULL;
00078     m_target          = NULL;
00079     m_overlap         = NULL;
00080     m_covering_source = NULL;
00081 
00082     point_cloud_source = false;
00083     point_cloud_target = false;
00084 
00085     return MB_SUCCESS;
00086 }
00087 
00088 ///////////////////////////////////////////////////////////////////////////////////
00089 
00090 TempestRemapper::~TempestRemapper()
00091 {
00092     this->clear();
00093 }
00094 
00095 ErrorCode TempestRemapper::clear()
00096 {
00097     // destroy all meshes
00098     if( m_source )
00099     {
00100         delete m_source;
00101         m_source = NULL;
00102     }
00103     if( m_target )
00104     {
00105         delete m_target;
00106         m_target = NULL;
00107     }
00108     if( m_overlap )
00109     {
00110         delete m_overlap;
00111         m_overlap = NULL;
00112     }
00113     if( m_covering_source && size > 1 )
00114     {
00115         delete m_covering_source;
00116         m_covering_source = NULL;
00117     }
00118 
00119     point_cloud_source = false;
00120     point_cloud_target = false;
00121 
00122     m_source_entities.clear();
00123     m_source_vertices.clear();
00124     m_covering_source_entities.clear();
00125     m_covering_source_vertices.clear();
00126     m_target_entities.clear();
00127     m_target_vertices.clear();
00128     m_overlap_entities.clear();
00129     gid_to_lid_src.clear();
00130     gid_to_lid_tgt.clear();
00131     gid_to_lid_covsrc.clear();
00132     lid_to_gid_src.clear();
00133     lid_to_gid_tgt.clear();
00134     lid_to_gid_covsrc.clear();
00135 
00136     return MB_SUCCESS;
00137 }
00138 
00139 ///////////////////////////////////////////////////////////////////////////////////
00140 
00141 ErrorCode TempestRemapper::LoadMesh( Remapper::IntersectionContext ctx,
00142                                      std::string inputFilename,
00143                                      TempestMeshType type )
00144 {
00145     if( ctx == Remapper::SourceMesh )
00146     {
00147         m_source_type = type;
00148         return load_tempest_mesh_private( inputFilename, &m_source );
00149     }
00150     else if( ctx == Remapper::TargetMesh )
00151     {
00152         m_target_type = type;
00153         return load_tempest_mesh_private( inputFilename, &m_target );
00154     }
00155     else if( ctx != Remapper::DEFAULT )
00156     {
00157         m_overlap_type = type;
00158         return load_tempest_mesh_private( inputFilename, &m_overlap );
00159     }
00160     else
00161     {
00162         MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
00163     }
00164 }
00165 
00166 ErrorCode TempestRemapper::load_tempest_mesh_private( std::string inputFilename, Mesh** tempest_mesh )
00167 {
00168     const bool outputEnabled = ( TempestRemapper::verbose && is_root );
00169     if( outputEnabled ) std::cout << "\nLoading TempestRemap Mesh object from file = " << inputFilename << " ...\n";
00170 
00171     {
00172         NcError error( NcError::silent_nonfatal );
00173 
00174         try
00175         {
00176             // Load input mesh
00177             if( outputEnabled ) std::cout << "Loading mesh ...\n";
00178             Mesh* mesh = new Mesh( inputFilename );
00179             mesh->RemoveZeroEdges();
00180             if( outputEnabled ) std::cout << "----------------\n";
00181 
00182             // Validate mesh
00183             if( meshValidate )
00184             {
00185                 if( outputEnabled ) std::cout << "Validating mesh ...\n";
00186                 mesh->Validate();
00187                 if( outputEnabled ) std::cout << "-------------------\n";
00188             }
00189 
00190             // Construct the edge map on the mesh
00191             if( constructEdgeMap )
00192             {
00193                 if( outputEnabled ) std::cout << "Constructing edge map on mesh ...\n";
00194                 mesh->ConstructEdgeMap( false );
00195                 if( outputEnabled ) std::cout << "---------------------------------\n";
00196             }
00197 
00198             if( tempest_mesh ) *tempest_mesh = mesh;
00199         }
00200         catch( Exception& e )
00201         {
00202             std::cout << "TempestRemap ERROR: " << e.ToString() << "\n";
00203             return MB_FAILURE;
00204         }
00205         catch( ... )
00206         {
00207             return MB_FAILURE;
00208         }
00209     }
00210     return MB_SUCCESS;
00211 }
00212 
00213 ///////////////////////////////////////////////////////////////////////////////////
00214 
00215 ErrorCode TempestRemapper::ConvertTempestMesh( Remapper::IntersectionContext ctx )
00216 {
00217     const bool outputEnabled = ( TempestRemapper::verbose && is_root );
00218     if( ctx == Remapper::SourceMesh )
00219     {
00220         if( outputEnabled ) std::cout << "Converting (source) TempestRemap Mesh object to MOAB representation ...\n";
00221         return convert_tempest_mesh_private( m_source_type, m_source, m_source_set, m_source_entities,
00222                                              &m_source_vertices );
00223     }
00224     else if( ctx == Remapper::TargetMesh )
00225     {
00226         if( outputEnabled ) std::cout << "Converting (target) TempestRemap Mesh object to MOAB representation ...\n";
00227         return convert_tempest_mesh_private( m_target_type, m_target, m_target_set, m_target_entities,
00228                                              &m_target_vertices );
00229     }
00230     else if( ctx != Remapper::DEFAULT )
00231     {
00232         if( outputEnabled ) std::cout << "Converting (overlap) TempestRemap Mesh object to MOAB representation ...\n";
00233         return convert_tempest_mesh_private( m_overlap_type, m_overlap, m_overlap_set, m_overlap_entities, NULL );
00234     }
00235     else
00236     {
00237         MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
00238     }
00239 }
00240 
00241 ErrorCode TempestRemapper::convert_tempest_mesh_private( TempestMeshType meshType,
00242                                                          Mesh* mesh,
00243                                                          EntityHandle& mesh_set,
00244                                                          Range& entities,
00245                                                          Range* vertices )
00246 {
00247     ErrorCode rval;
00248 
00249     const bool outputEnabled = ( TempestRemapper::verbose && is_root );
00250     const NodeVector& nodes  = mesh->nodes;
00251     const FaceVector& faces  = mesh->faces;
00252 
00253     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
00254     dbgprint.set_prefix( "[TempestToMOAB]: " );
00255 
00256     ReadUtilIface* iface;
00257     rval = m_interface->query_interface( iface );MB_CHK_SET_ERR( rval, "Can't get reader interface" );
00258 
00259     Tag gidTag = m_interface->globalId_tag();
00260 
00261     // Set the data for the vertices
00262     std::vector< double* > arrays;
00263     std::vector< int > gidsv( nodes.size() );
00264     EntityHandle startv;
00265     rval = iface->get_node_coords( 3, nodes.size(), 0, startv, arrays );MB_CHK_SET_ERR( rval, "Can't get node coords" );
00266     for( unsigned iverts = 0; iverts < nodes.size(); ++iverts )
00267     {
00268         const Node& node  = nodes[iverts];
00269         arrays[0][iverts] = node.x;
00270         arrays[1][iverts] = node.y;
00271         arrays[2][iverts] = node.z;
00272         gidsv[iverts]     = iverts + 1;
00273     }
00274     Range mbverts( startv, startv + nodes.size() - 1 );
00275     rval = m_interface->add_entities( mesh_set, mbverts );MB_CHK_SET_ERR( rval, "Can't add entities" );
00276     rval = m_interface->tag_set_data( gidTag, mbverts, &gidsv[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00277 
00278     gidsv.clear();
00279     entities.clear();
00280 
00281     Tag srcParentTag, tgtParentTag;
00282     std::vector< int > srcParent, tgtParent;
00283     bool storeParentInfo = ( mesh->vecSourceFaceIx.size() > 0 );
00284 
00285     if( storeParentInfo )
00286     {
00287         int defaultInt = -1;
00288         rval           = m_interface->tag_get_handle( "TargetParent", 1, MB_TYPE_INTEGER, tgtParentTag,
00289                                             MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create positive tag" );
00290 
00291         rval = m_interface->tag_get_handle( "SourceParent", 1, MB_TYPE_INTEGER, srcParentTag,
00292                                             MB_TAG_DENSE | MB_TAG_CREAT, &defaultInt );MB_CHK_SET_ERR( rval, "can't create negative tag" );
00293     }
00294 
00295     // Let us first perform a full pass assuming arbitrary polygons. This is especially true for
00296     // overlap meshes.
00297     //   1. We do a first pass over faces, decipher edge size and group into categories based on
00298     //   element type
00299     //   2. Next we loop over type, and add blocks of elements into MOAB
00300     //   3. For each block within the loop, also update the connectivity of elements.
00301     {
00302         if( outputEnabled )
00303             dbgprint.printf( 0, "..Mesh size: Nodes [%zu]  Elements [%zu].\n", nodes.size(), faces.size() );
00304         const int NMAXPOLYEDGES = 15;
00305         std::vector< unsigned > nPolys( NMAXPOLYEDGES, 0 );
00306         std::vector< std::vector< int > > typeNSeqs( NMAXPOLYEDGES );
00307         for( unsigned ifaces = 0; ifaces < faces.size(); ++ifaces )
00308         {
00309             const int iType = faces[ifaces].edges.size();
00310             nPolys[iType]++;
00311             typeNSeqs[iType].push_back( ifaces );
00312         }
00313         int iBlock = 0;
00314         for( unsigned iType = 0; iType < NMAXPOLYEDGES; ++iType )
00315         {
00316             if( !nPolys[iType] ) continue;  // Nothing to do
00317 
00318             const unsigned num_v_per_elem = iType;
00319             EntityHandle starte;  // Connectivity
00320             EntityHandle* conn;
00321 
00322             // Allocate the connectivity array, depending on the element type
00323             switch( num_v_per_elem )
00324             {
00325                 case 3:
00326                     if( outputEnabled )
00327                         dbgprint.printf( 0, "....Block %d: Triangular Elements [%u].\n", iBlock++, nPolys[iType] );
00328                     rval = iface->get_element_connect( nPolys[iType], num_v_per_elem, MBTRI, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00329                     break;
00330                 case 4:
00331                     if( outputEnabled )
00332                         dbgprint.printf( 0, "....Block %d: Quadrilateral Elements [%u].\n", iBlock++, nPolys[iType] );
00333                     rval = iface->get_element_connect( nPolys[iType], num_v_per_elem, MBQUAD, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00334                     break;
00335                 default:
00336                     if( outputEnabled )
00337                         dbgprint.printf( 0, "....Block %d: Polygonal [%u] Elements [%u].\n", iBlock++, iType,
00338                                          nPolys[iType] );
00339                     rval = iface->get_element_connect( nPolys[iType], num_v_per_elem, MBPOLYGON, 0, starte, conn );MB_CHK_SET_ERR( rval, "Can't get element connectivity" );
00340                     break;
00341             }
00342 
00343             Range mbcells( starte, starte + nPolys[iType] - 1 );
00344             m_interface->add_entities( mesh_set, mbcells );
00345 
00346             if( storeParentInfo )
00347             {
00348                 srcParent.resize( mbcells.size(), -1 );
00349                 tgtParent.resize( mbcells.size(), -1 );
00350             }
00351 
00352             std::vector< int > gids( typeNSeqs[iType].size() );
00353             for( unsigned ifaces = 0, offset = 0; ifaces < typeNSeqs[iType].size(); ++ifaces )
00354             {
00355                 const int fIndex = typeNSeqs[iType][ifaces];
00356                 const Face& face = faces[fIndex];
00357                 // conn[offset++] = startv + face.edges[0].node[0];
00358                 for( unsigned iedges = 0; iedges < face.edges.size(); ++iedges )
00359                 {
00360                     conn[offset++] = startv + face.edges[iedges].node[0];
00361                 }
00362 
00363                 if( storeParentInfo )
00364                 {
00365                     srcParent[ifaces] = mesh->vecSourceFaceIx[fIndex] + 1;
00366                     tgtParent[ifaces] = mesh->vecTargetFaceIx[fIndex] + 1;
00367                 }
00368 
00369                 gids[ifaces] = typeNSeqs[iType][ifaces] + 1;
00370             }
00371             rval = m_interface->tag_set_data( gidTag, mbcells, &gids[0] );MB_CHK_SET_ERR( rval, "Can't set global_id tag" );
00372 
00373             if( meshType == OVERLAP_FILES )
00374             {
00375                 // Now let us update the adjacency data, because some elements are new
00376                 rval = iface->update_adjacencies( starte, nPolys[iType], num_v_per_elem, conn );MB_CHK_SET_ERR( rval, "Can't update adjacencies" );
00377                 // Generate all adj entities dimension 1 and 2 (edges and faces/ tri or qua)
00378                 Range edges;
00379                 rval = m_interface->get_adjacencies( mbcells, 1, true, edges, Interface::UNION );MB_CHK_SET_ERR( rval, "Can't get edges" );
00380             }
00381 
00382             if( storeParentInfo )
00383             {
00384                 rval = m_interface->tag_set_data( srcParentTag, mbcells, &srcParent[0] );MB_CHK_SET_ERR( rval, "Can't set tag data" );
00385                 rval = m_interface->tag_set_data( tgtParentTag, mbcells, &tgtParent[0] );MB_CHK_SET_ERR( rval, "Can't set tag data" );
00386             }
00387             entities.merge( mbcells );
00388         }
00389     }
00390 
00391     if( vertices ) *vertices = mbverts;
00392 
00393     return MB_SUCCESS;
00394 }
00395 
00396 ///////////////////////////////////////////////////////////////////////////////////
00397 
00398 ErrorCode TempestRemapper::ConvertMeshToTempest( Remapper::IntersectionContext ctx )
00399 {
00400     ErrorCode rval;
00401     const bool outputEnabled = ( TempestRemapper::verbose && is_root );
00402 
00403     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
00404     dbgprint.set_prefix( "[MOABToTempest]: " );
00405 
00406     if( ctx == Remapper::SourceMesh )
00407     {
00408         if( !m_source ) m_source = new Mesh();
00409         if( outputEnabled ) dbgprint.printf( 0, "Converting (source) MOAB to TempestRemap Mesh representation ...\n" );
00410         rval = convert_mesh_to_tempest_private( m_source, m_source_set, m_source_entities, &m_source_vertices );MB_CHK_SET_ERR( rval, "Can't convert source mesh to Tempest" );
00411         if( m_source_entities.size() == 0 && m_source_vertices.size() != 0 )
00412         {
00413             this->point_cloud_source = true;
00414         }
00415     }
00416     else if( ctx == Remapper::TargetMesh )
00417     {
00418         if( !m_target ) m_target = new Mesh();
00419         if( outputEnabled ) dbgprint.printf( 0, "Converting (target) MOAB to TempestRemap Mesh representation ...\n" );
00420         rval = convert_mesh_to_tempest_private( m_target, m_target_set, m_target_entities, &m_target_vertices );MB_CHK_SET_ERR( rval, "Can't convert target mesh to Tempest" );
00421         if( m_target_entities.size() == 0 && m_target_vertices.size() != 0 ) this->point_cloud_target = true;
00422     }
00423     else if( ctx == Remapper::OverlapMesh )  // Overlap mesh
00424     {
00425         if( !m_overlap ) m_overlap = new Mesh();
00426         if( outputEnabled ) dbgprint.printf( 0, "Converting (overlap) MOAB to TempestRemap Mesh representation ...\n" );
00427         rval = convert_overlap_mesh_sorted_by_source();MB_CHK_SET_ERR( rval, "Can't convert overlap mesh to Tempest" );
00428     }
00429     else
00430     {
00431         MB_CHK_SET_ERR( MB_FAILURE, "Invalid IntersectionContext context provided" );
00432     }
00433 
00434     return rval;
00435 }
00436 
00437 ErrorCode TempestRemapper::convert_mesh_to_tempest_private( Mesh* mesh,
00438                                                             EntityHandle mesh_set,
00439                                                             moab::Range& elems,
00440                                                             moab::Range* pverts )
00441 {
00442     ErrorCode rval;
00443     Range verts;
00444 
00445     NodeVector& nodes = mesh->nodes;
00446     FaceVector& faces = mesh->faces;
00447 
00448     elems.clear();
00449     rval = m_interface->get_entities_by_dimension( mesh_set, 2, elems );MB_CHK_ERR( rval );
00450 
00451     // resize the number of elements in Tempest mesh
00452     faces.resize( elems.size() );
00453 
00454     // let us now get the vertices from all the elements
00455     rval = m_interface->get_connectivity( elems, verts );MB_CHK_ERR( rval );
00456     if( verts.size() == 0 )
00457     {
00458         rval = m_interface->get_entities_by_dimension( mesh_set, 0, verts );MB_CHK_ERR( rval );
00459     }
00460     // assert(verts.size() > 0); // If not, this may be an invalid mesh ! possible for unbalanced
00461     // loads
00462 
00463     std::map< EntityHandle, int > indxMap;
00464     bool useRange = true;
00465     if( verts.compactness() > 0.01 )
00466     {
00467         int j = 0;
00468         for( Range::iterator it = verts.begin(); it != verts.end(); it++ )
00469             indxMap[*it] = j++;
00470         useRange = false;
00471     }
00472 
00473     for( unsigned iface = 0; iface < elems.size(); ++iface )
00474     {
00475         Face& face           = faces[iface];
00476         EntityHandle ehandle = elems[iface];
00477 
00478         // get the connectivity for each edge
00479         const EntityHandle* connectface;
00480         int nnodesf;
00481         rval = m_interface->get_connectivity( ehandle, connectface, nnodesf );MB_CHK_ERR( rval );
00482 
00483         face.edges.resize( nnodesf );
00484         for( int iverts = 0; iverts < nnodesf; ++iverts )
00485         {
00486             int indx = ( useRange ? verts.index( connectface[iverts] ) : indxMap[connectface[iverts]] );
00487             assert( indx >= 0 );
00488             face.SetNode( iverts, indx );
00489         }
00490     }
00491 
00492     unsigned nnodes = verts.size();
00493     nodes.resize( nnodes );
00494 
00495     // Set the data for the vertices
00496     std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
00497     rval = m_interface->get_coords( verts, &coordx[0], &coordy[0], &coordz[0] );MB_CHK_ERR( rval );
00498     for( unsigned inode = 0; inode < nnodes; ++inode )
00499     {
00500         Node& node = nodes[inode];
00501         node.x     = coordx[inode];
00502         node.y     = coordy[inode];
00503         node.z     = coordz[inode];
00504     }
00505     coordx.clear();
00506     coordy.clear();
00507     coordz.clear();
00508 
00509     mesh->RemoveZeroEdges();
00510     mesh->RemoveCoincidentNodes();
00511 
00512     // Generate reverse node array and edge map
00513     if( constructEdgeMap ) mesh->ConstructEdgeMap( false );
00514     // mesh->ConstructReverseNodeArray();
00515 
00516     // mesh->Validate();
00517 
00518     if( pverts )
00519     {
00520         pverts->clear();
00521         *pverts = verts;
00522     }
00523     verts.clear();
00524 
00525     return MB_SUCCESS;
00526 }
00527 
00528 ///////////////////////////////////////////////////////////////////////////////////
00529 
00530 bool IntPairComparator( const std::pair< int, int >& a, const std::pair< int, int >& b )
00531 {
00532     if( a.first == b.first )
00533         return a.second < b.second;
00534     else
00535         return a.first < b.first;
00536 }
00537 
00538 moab::ErrorCode moab::TempestRemapper::GetOverlapAugmentedEntities( moab::Range& sharedGhostEntities )
00539 {
00540     sharedGhostEntities.clear();
00541 #ifdef MOAB_HAVE_MPI
00542     moab::ErrorCode rval;
00543 
00544     // Remove entities in the intersection mesh that are part of the ghosted overlap
00545     if( is_parallel && size > 1 )
00546     {
00547         moab::Range allents;
00548         rval = m_interface->get_entities_by_dimension( m_overlap_set, 2, allents );MB_CHK_SET_ERR( rval, "Getting entities dim 2 failed" );
00549 
00550         moab::Range sharedents;
00551         moab::Tag ghostTag;
00552         std::vector< int > ghFlags( allents.size() );
00553         rval = m_interface->tag_get_handle( "ORIG_PROC", ghostTag );MB_CHK_ERR( rval );
00554         rval = m_interface->tag_get_data( ghostTag, allents, &ghFlags[0] );MB_CHK_ERR( rval );
00555         for( unsigned i = 0; i < allents.size(); ++i )
00556             if( ghFlags[i] >= 0 )                 // it means it is a ghost overlap element
00557                 sharedents.insert( allents[i] );  // this should not participate in smat!
00558 
00559         allents = subtract( allents, sharedents );
00560 
00561         // Get connectivity from all ghosted elements and filter out
00562         // the vertices that are not owned
00563         moab::Range ownedverts, sharedverts;
00564         rval = m_interface->get_connectivity( allents, ownedverts );MB_CHK_SET_ERR( rval, "Deleting entities dim 0 failed" );
00565         rval = m_interface->get_connectivity( sharedents, sharedverts );MB_CHK_SET_ERR( rval, "Deleting entities dim 0 failed" );
00566         sharedverts = subtract( sharedverts, ownedverts );
00567         // rval = m_interface->remove_entities(m_overlap_set, sharedents);MB_CHK_SET_ERR(rval,
00568         // "Deleting entities dim 2 failed"); rval = m_interface->remove_entities(m_overlap_set,
00569         // sharedverts);MB_CHK_SET_ERR(rval, "Deleting entities dim 0 failed");
00570 
00571         sharedGhostEntities.merge( sharedents );
00572         // sharedGhostEntities.merge(sharedverts);
00573     }
00574 #endif
00575     return moab::MB_SUCCESS;
00576 }
00577 
00578 ErrorCode TempestRemapper::convert_overlap_mesh_sorted_by_source()
00579 {
00580     ErrorCode rval;
00581 
00582     m_overlap_entities.clear();
00583     rval = m_interface->get_entities_by_dimension( m_overlap_set, 2, m_overlap_entities );MB_CHK_ERR( rval );
00584 
00585     // Allocate for the overlap mesh
00586     if( !m_overlap ) m_overlap = new Mesh();
00587 
00588     size_t n_overlap_entities = m_overlap_entities.size();
00589 
00590     std::vector< std::pair< int, int > > sorted_overlap_order( n_overlap_entities );
00591     {
00592         Tag srcParentTag, tgtParentTag;
00593         rval = m_interface->tag_get_handle( "SourceParent", srcParentTag );MB_CHK_ERR( rval );
00594         rval = m_interface->tag_get_handle( "TargetParent", tgtParentTag );MB_CHK_ERR( rval );
00595         // Overlap mesh: resize the source and target connection arrays
00596         m_overlap->vecTargetFaceIx.resize( n_overlap_entities );
00597         m_overlap->vecSourceFaceIx.resize( n_overlap_entities );
00598 
00599         // Overlap mesh: resize the source and target connection arrays
00600         std::vector< int > rbids_src( n_overlap_entities ), rbids_tgt( n_overlap_entities );
00601         rval = m_interface->tag_get_data( srcParentTag, m_overlap_entities, &rbids_src[0] );MB_CHK_ERR( rval );
00602         rval = m_interface->tag_get_data( tgtParentTag, m_overlap_entities, &rbids_tgt[0] );MB_CHK_ERR( rval );
00603         for( size_t ix = 0; ix < n_overlap_entities; ++ix )
00604         {
00605             sorted_overlap_order[ix].first =
00606                 ( gid_to_lid_covsrc.size() ? gid_to_lid_covsrc[rbids_src[ix]] : rbids_src[ix] );
00607             sorted_overlap_order[ix].second = ix;
00608         }
00609         std::sort( sorted_overlap_order.begin(), sorted_overlap_order.end(), IntPairComparator );
00610         // sorted_overlap_order[ie].second , ie=0,nOverlap-1 is the order such that overlap elems
00611         // are ordered by source parent
00612 
00613         std::vector< int > ghFlags;
00614         if( is_parallel && size > 1 )
00615         {
00616             Tag ghostTag;
00617             ghFlags.resize( n_overlap_entities );
00618             rval = m_interface->tag_get_handle( "ORIG_PROC", ghostTag );MB_CHK_ERR( rval );
00619             rval = m_interface->tag_get_data( ghostTag, m_overlap_entities, &ghFlags[0] );MB_CHK_ERR( rval );
00620         }
00621         for( unsigned ie = 0; ie < n_overlap_entities; ++ie )
00622         {
00623             int ix = sorted_overlap_order[ie].second;  // original index of the element
00624             m_overlap->vecSourceFaceIx[ie] =
00625                 ( gid_to_lid_covsrc.size() ? gid_to_lid_covsrc[rbids_src[ix]] : rbids_src[ix] - 1 );
00626             if( is_parallel && size > 1 && ghFlags[ix] >= 0 )  // it means it is a ghost overlap element
00627                 m_overlap->vecTargetFaceIx[ie] = -1;           // this should not participate in smat!
00628             else
00629                 m_overlap->vecTargetFaceIx[ie] =
00630                     ( gid_to_lid_tgt.size() ? gid_to_lid_tgt[rbids_tgt[ix]] : rbids_tgt[ix] - 1 );
00631         }
00632     }
00633 
00634     FaceVector& faces = m_overlap->faces;
00635     faces.resize( n_overlap_entities );
00636 
00637     Range verts;
00638     // let us now get the vertices from all the elements
00639     rval = m_interface->get_connectivity( m_overlap_entities, verts );MB_CHK_ERR( rval );
00640     // std::cout << "Vertices size = " << verts.size() << " , psize = " << verts.psize() << ",
00641     // compactness = " << verts.compactness() << std::endl;
00642 
00643     std::map< EntityHandle, int > indxMap;
00644     bool useRange = true;
00645     if( verts.compactness() > 0.01 )
00646     {
00647         int j = 0;
00648         for( Range::iterator it = verts.begin(); it != verts.end(); ++it )
00649             indxMap[*it] = j++;
00650         useRange = false;
00651     }
00652 
00653     for( unsigned ifac = 0; ifac < m_overlap_entities.size(); ++ifac )
00654     {
00655         const unsigned iface = sorted_overlap_order[ifac].second;
00656         Face& face           = faces[ifac];
00657         EntityHandle ehandle = m_overlap_entities[iface];
00658 
00659         // get the connectivity for each edge
00660         const EntityHandle* connectface;
00661         int nnodesf;
00662         rval = m_interface->get_connectivity( ehandle, connectface, nnodesf );MB_CHK_ERR( rval );
00663 
00664         face.edges.resize( nnodesf );
00665         for( int iverts = 0; iverts < nnodesf; ++iverts )
00666         {
00667             int indx = ( useRange ? verts.index( connectface[iverts] ) : indxMap[connectface[iverts]] );
00668             assert( indx >= 0 );
00669             face.SetNode( iverts, indx );
00670         }
00671     }
00672 
00673     unsigned nnodes   = verts.size();
00674     NodeVector& nodes = m_overlap->nodes;
00675     nodes.resize( nnodes );
00676 
00677     // Set the data for the vertices
00678     std::vector< double > coordx( nnodes ), coordy( nnodes ), coordz( nnodes );
00679     rval = m_interface->get_coords( verts, &coordx[0], &coordy[0], &coordz[0] );MB_CHK_ERR( rval );
00680     for( unsigned inode = 0; inode < nnodes; ++inode )
00681     {
00682         Node& node = nodes[inode];
00683         node.x     = coordx[inode];
00684         node.y     = coordy[inode];
00685         node.z     = coordz[inode];
00686     }
00687     coordx.clear();
00688     coordy.clear();
00689     coordz.clear();
00690     verts.clear();
00691 
00692     m_overlap->RemoveZeroEdges();
00693     m_overlap->RemoveCoincidentNodes( false );
00694 
00695     // Generate reverse node array and edge map
00696     // if ( constructEdgeMap ) m_overlap->ConstructEdgeMap(false);
00697     // m_overlap->ConstructReverseNodeArray();
00698 
00699     // m_overlap->Validate();
00700     return MB_SUCCESS;
00701 }
00702 
00703 // Should be ordered as Source, Target, Overlap
00704 ErrorCode TempestRemapper::ComputeGlobalLocalMaps()
00705 {
00706     ErrorCode rval;
00707 
00708     if( 0 == m_covering_source )
00709     {
00710         m_covering_source = new Mesh();
00711         rval = convert_mesh_to_tempest_private( m_covering_source, m_covering_source_set, m_covering_source_entities,
00712                                                 &m_covering_source_vertices );MB_CHK_SET_ERR( rval, "Can't convert source Tempest mesh" );
00713         // std::cout << "ComputeGlobalLocalMaps: " << rank << ", " << " covering entities = [" <<
00714         // m_covering_source_vertices.size() << ", " << m_covering_source_entities.size() << "]\n";
00715     }
00716     gid_to_lid_src.clear();
00717     lid_to_gid_src.clear();
00718     gid_to_lid_covsrc.clear();
00719     lid_to_gid_covsrc.clear();
00720     gid_to_lid_tgt.clear();
00721     lid_to_gid_tgt.clear();
00722     {
00723         Tag gidtag = m_interface->globalId_tag();
00724 
00725         std::vector< int > gids;
00726         if( point_cloud_source )
00727         {
00728             gids.resize( m_covering_source_vertices.size(), -1 );
00729             rval = m_interface->tag_get_data( gidtag, m_covering_source_vertices, &gids[0] );MB_CHK_ERR( rval );
00730         }
00731         else
00732         {
00733             gids.resize( m_covering_source_entities.size(), -1 );
00734             rval = m_interface->tag_get_data( gidtag, m_covering_source_entities, &gids[0] );MB_CHK_ERR( rval );
00735         }
00736         for( unsigned ie = 0; ie < gids.size(); ++ie )
00737         {
00738             gid_to_lid_covsrc[gids[ie]] = ie;
00739             lid_to_gid_covsrc[ie]       = gids[ie];
00740         }
00741 
00742         if( point_cloud_source )
00743         {
00744             gids.resize( m_source_vertices.size(), -1 );
00745             rval = m_interface->tag_get_data( gidtag, m_source_vertices, &gids[0] );MB_CHK_ERR( rval );
00746         }
00747         else
00748         {
00749             gids.resize( m_source_entities.size(), -1 );
00750             rval = m_interface->tag_get_data( gidtag, m_source_entities, &gids[0] );MB_CHK_ERR( rval );
00751         }
00752         for( unsigned ie = 0; ie < gids.size(); ++ie )
00753         {
00754             gid_to_lid_src[gids[ie]] = ie;
00755             lid_to_gid_src[ie]       = gids[ie];
00756         }
00757 
00758         if( point_cloud_target )
00759         {
00760             gids.resize( m_target_vertices.size(), -1 );
00761             rval = m_interface->tag_get_data( gidtag, m_target_vertices, &gids[0] );MB_CHK_ERR( rval );
00762         }
00763         else
00764         {
00765             gids.resize( m_target_entities.size(), -1 );
00766             rval = m_interface->tag_get_data( gidtag, m_target_entities, &gids[0] );MB_CHK_ERR( rval );
00767         }
00768         for( unsigned ie = 0; ie < gids.size(); ++ie )
00769         {
00770             gid_to_lid_tgt[gids[ie]] = ie;
00771             lid_to_gid_tgt[ie]       = gids[ie];
00772         }
00773     }
00774 
00775     return MB_SUCCESS;
00776 }
00777 
00778 ///////////////////////////////////////////////////////////////////////////////////
00779 
00780 moab::ErrorCode moab::TempestRemapper::WriteTempestIntersectionMesh( std::string strOutputFileName,
00781                                                                      const bool fAllParallel,
00782                                                                      const bool fInputConcave,
00783                                                                      const bool fOutputConcave )
00784 {
00785     // Let us alos write out the TempestRemap equivalent so that we can do some verification checks
00786     if( fAllParallel )
00787     {
00788         if( is_root && size == 1 )
00789         {
00790             this->m_source->CalculateFaceAreas( fInputConcave );
00791             this->m_target->CalculateFaceAreas( fOutputConcave );
00792             this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
00793         }
00794         else
00795         {
00796             // Perform reduction and write from root processor
00797             // if ( is_root )
00798             //     std::cout << "--- PARALLEL IMPLEMENTATION is NOT AVAILABLE yet ---\n";
00799 
00800             this->m_source->CalculateFaceAreas( fInputConcave );
00801             this->m_covering_source->CalculateFaceAreas( fInputConcave );
00802             this->m_target->CalculateFaceAreas( fOutputConcave );
00803             this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
00804         }
00805     }
00806     else
00807     {
00808         this->m_source->CalculateFaceAreas( fInputConcave );
00809         this->m_target->CalculateFaceAreas( fOutputConcave );
00810         this->m_overlap->Write( strOutputFileName.c_str(), NcFile::Netcdf4 );
00811     }
00812 
00813     return moab::MB_SUCCESS;
00814 }
00815 void TempestRemapper::SetMeshSet( Remapper::IntersectionContext ctx /* Remapper::CoveringMesh*/,
00816                                   moab::EntityHandle mset,
00817                                   moab::Range& entities )
00818 {
00819 
00820     if( ctx == Remapper::SourceMesh )  // should not be used
00821     {
00822         m_source_entities = entities;
00823         m_source_set      = mset;
00824     }
00825     else if( ctx == Remapper::TargetMesh )
00826     {
00827         m_target_entities = entities;
00828         m_target_set      = mset;
00829     }
00830     else if( ctx == Remapper::CoveringMesh )
00831     {
00832         m_covering_source_entities = entities;
00833         m_covering_source_set      = mset;
00834     }
00835     else
00836     {
00837         // some error
00838     }
00839     return;
00840 }
00841 ///////////////////////////////////////////////////////////////////////////////////
00842 
00843 #ifndef MOAB_HAVE_MPI
00844 ErrorCode TempestRemapper::assign_vertex_element_IDs( Tag idtag,
00845                                                       EntityHandle this_set,
00846                                                       const int dimension,
00847                                                       const int start_id )
00848 {
00849     assert( idtag );
00850 
00851     ErrorCode rval;
00852     Range entities;
00853     rval = m_interface->get_entities_by_dimension( this_set, dimension, entities );MB_CHK_SET_ERR( rval, "Failed to get entities" );
00854 
00855     if( entities.size() == 0 ) return moab::MB_SUCCESS;
00856 
00857     int idoffset = start_id;
00858     std::vector< int > gid( entities.size() );
00859     for( unsigned i = 0; i < entities.size(); ++i )
00860         gid[i] = idoffset++;
00861 
00862     rval = m_interface->tag_set_data( idtag, entities, &gid[0] );MB_CHK_ERR( rval );
00863 
00864     return moab::MB_SUCCESS;
00865 }
00866 #endif
00867 
00868 ///////////////////////////////////////////////////////////////////////////////
00869 
00870 // Create a custom comparator for Nodes
00871 bool operator<( Node const& lhs, Node const& rhs )
00872 {
00873     return std::pow( lhs.x - rhs.x, 2.0 ) + std::pow( lhs.y - rhs.y, 2.0 ) + std::pow( lhs.z - rhs.z, 2.0 );
00874 }
00875 
00876 ErrorCode TempestRemapper::GenerateCSMeshMetadata( const int ntot_elements,
00877                                                    moab::Range& ents,
00878                                                    moab::Range* secondary_ents,
00879                                                    const std::string& dofTagName,
00880                                                    int nP )
00881 {
00882     Mesh csMesh;
00883     int err;
00884     moab::ErrorCode rval;
00885 
00886     const int res = std::sqrt( ntot_elements / 6 );
00887 
00888     // create a temporary CS mesh
00889     // NOTE: This will not work for RRM grids. Need to run HOMME for that case anyway
00890     err = GenerateCSMesh( csMesh, res, "", "NetCDF4" );
00891     if( err )
00892     {
00893         MB_CHK_SET_ERR( MB_FAILURE, "Failed to generate CS mesh through TempestRemap" );
00894         ;
00895     }
00896 
00897     rval = this->GenerateMeshMetadata( csMesh, ntot_elements, ents, secondary_ents, dofTagName, nP );MB_CHK_SET_ERR( rval, "Failed in call to GenerateMeshMetadata" );
00898 
00899     return moab::MB_SUCCESS;
00900 }
00901 
00902 ErrorCode TempestRemapper::GenerateMeshMetadata( Mesh& csMesh,
00903                                                  const int ntot_elements,
00904                                                  moab::Range& ents,
00905                                                  moab::Range* secondary_ents,
00906                                                  const std::string dofTagName,
00907                                                  int nP )
00908 {
00909     moab::ErrorCode rval;
00910 
00911     Tag dofTag;
00912     bool created = false;
00913     rval         = m_interface->tag_get_handle( dofTagName.c_str(), nP * nP, MB_TYPE_INTEGER, dofTag,
00914                                         MB_TAG_DENSE | MB_TAG_CREAT, 0, &created );MB_CHK_SET_ERR( rval, "Failed creating DoF tag" );
00915 
00916     // Number of Faces
00917     int nElements = static_cast< int >( csMesh.faces.size() );
00918 
00919     assert( nElements == ntot_elements );
00920 
00921     // Initialize data structures
00922     DataArray3D< int > dataGLLnodes;
00923     dataGLLnodes.Allocate( nP, nP, nElements );
00924 
00925     std::map< Node, int > mapNodes;
00926     std::map< Node, moab::EntityHandle > mapLocalMBNodes;
00927 
00928     // GLL Quadrature nodes
00929     DataArray1D< double > dG;
00930     DataArray1D< double > dW;
00931     GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
00932 
00933     moab::Range entities( ents );
00934     if( secondary_ents ) entities.insert( secondary_ents->begin(), secondary_ents->end() );
00935     double elcoords[3];
00936     for( unsigned iel = 0; iel < entities.size(); ++iel )
00937     {
00938         EntityHandle eh = entities[iel];
00939         rval            = m_interface->get_coords( &eh, 1, elcoords );
00940         Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
00941         mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
00942     }
00943 
00944     // Build a Kd-tree for local mesh (nearest neighbor searches)
00945     // Loop over all elements in CS-Mesh
00946     // Then find if current centroid is in an element
00947     //     If yes - then let us compute the DoF numbering and set to tag data
00948     //     If no - then compute DoF numbering BUT DO NOT SET to tag data
00949     // continue
00950     int* dofIDs = new int[nP * nP];
00951 
00952     // Write metadata
00953     for( int k = 0; k < nElements; k++ )
00954     {
00955         const Face& face        = csMesh.faces[k];
00956         const NodeVector& nodes = csMesh.nodes;
00957 
00958         if( face.edges.size() != 4 )
00959         {
00960             _EXCEPTIONT( "Mesh must only contain quadrilateral elements" );
00961         }
00962 
00963         Node centroid;
00964         centroid.x = centroid.y = centroid.z = 0.0;
00965         for( unsigned l = 0; l < face.edges.size(); ++l )
00966         {
00967             centroid.x += nodes[face[l]].x;
00968             centroid.y += nodes[face[l]].y;
00969             centroid.z += nodes[face[l]].z;
00970         }
00971         const double factor = 1.0 / face.edges.size();
00972         centroid.x *= factor;
00973         centroid.y *= factor;
00974         centroid.z *= factor;
00975 
00976         bool locElem = false;
00977         EntityHandle current_eh;
00978         if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
00979         {
00980             locElem    = true;
00981             current_eh = mapLocalMBNodes[centroid];
00982         }
00983 
00984         for( int j = 0; j < nP; j++ )
00985         {
00986             for( int i = 0; i < nP; i++ )
00987             {
00988 
00989                 // Get local map vectors
00990                 Node nodeGLL;
00991                 Node dDx1G;
00992                 Node dDx2G;
00993 
00994                 // ApplyLocalMap(
00995                 //     face,
00996                 //     nodevec,
00997                 //     dG[i],
00998                 //     dG[j],
00999                 //     nodeGLL,
01000                 //     dDx1G,
01001                 //     dDx2G);
01002                 const double& dAlpha = dG[i];
01003                 const double& dBeta  = dG[j];
01004 
01005                 // Calculate nodal locations on the plane
01006                 double dXc = nodes[face[0]].x * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
01007                              nodes[face[1]].x * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].x * dAlpha * dBeta +
01008                              nodes[face[3]].x * ( 1.0 - dAlpha ) * dBeta;
01009 
01010                 double dYc = nodes[face[0]].y * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
01011                              nodes[face[1]].y * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].y * dAlpha * dBeta +
01012                              nodes[face[3]].y * ( 1.0 - dAlpha ) * dBeta;
01013 
01014                 double dZc = nodes[face[0]].z * ( 1.0 - dAlpha ) * ( 1.0 - dBeta ) +
01015                              nodes[face[1]].z * dAlpha * ( 1.0 - dBeta ) + nodes[face[2]].z * dAlpha * dBeta +
01016                              nodes[face[3]].z * ( 1.0 - dAlpha ) * dBeta;
01017 
01018                 double dR = sqrt( dXc * dXc + dYc * dYc + dZc * dZc );
01019 
01020                 // Mapped node location
01021                 nodeGLL.x = dXc / dR;
01022                 nodeGLL.y = dYc / dR;
01023                 nodeGLL.z = dZc / dR;
01024 
01025                 // Determine if this is a unique Node
01026                 std::map< Node, int >::const_iterator iter = mapNodes.find( nodeGLL );
01027                 if( iter == mapNodes.end() )
01028                 {
01029                     // Insert new unique node into map
01030                     int ixNode = static_cast< int >( mapNodes.size() );
01031                     mapNodes.insert( std::pair< Node, int >( nodeGLL, ixNode ) );
01032                     dataGLLnodes[j][i][k] = ixNode + 1;
01033                 }
01034                 else
01035                 {
01036                     dataGLLnodes[j][i][k] = iter->second + 1;
01037                 }
01038 
01039                 dofIDs[j * nP + i] = dataGLLnodes[j][i][k];
01040             }
01041         }
01042 
01043         if( locElem )
01044         {
01045             rval = m_interface->tag_set_data( dofTag, &current_eh, 1, dofIDs );MB_CHK_SET_ERR( rval, "Failed to tag_set_data for DoFs" );
01046         }
01047     }
01048 
01049     // clear memory
01050     delete[] dofIDs;
01051     mapLocalMBNodes.clear();
01052     mapNodes.clear();
01053 
01054     return moab::MB_SUCCESS;
01055 }
01056 
01057 ///////////////////////////////////////////////////////////////////////////////////
01058 
01059 ErrorCode TempestRemapper::ConstructCoveringSet( double tolerance,
01060                                                  double radius_src,
01061                                                  double radius_tgt,
01062                                                  double boxeps,
01063                                                  bool regional_mesh )
01064 {
01065     ErrorCode rval;
01066 
01067     rrmgrids = regional_mesh;
01068     moab::Range local_verts;
01069 
01070     // Initialize intersection context
01071     mbintx = new moab::Intx2MeshOnSphere( m_interface );
01072 
01073     mbintx->set_error_tolerance( tolerance );
01074     mbintx->set_radius_source_mesh( radius_src );
01075     mbintx->set_radius_destination_mesh( radius_tgt );
01076     mbintx->set_box_error( boxeps );
01077 #ifdef MOAB_HAVE_MPI
01078     mbintx->set_parallel_comm( m_pcomm );
01079 #endif
01080 
01081     // compute the maxiumum edges in elements comprising source and target mesh
01082     rval = mbintx->FindMaxEdges( m_source_set, m_target_set );MB_CHK_ERR( rval );
01083 
01084     this->max_source_edges = mbintx->max_edges_1;
01085     this->max_target_edges = mbintx->max_edges_2;
01086 
01087     // Note: lots of communication possible, if mesh is distributed very differently
01088 #ifdef MOAB_HAVE_MPI
01089     if( is_parallel )
01090     {
01091         rval = mbintx->build_processor_euler_boxes( m_target_set, local_verts );MB_CHK_ERR( rval );
01092 
01093         rval = m_interface->create_meshset( moab::MESHSET_SET, m_covering_source_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
01094 
01095         rval = mbintx->construct_covering_set( m_source_set, m_covering_source_set );MB_CHK_ERR( rval );
01096         // if (rank == 1)
01097         // {
01098         //     moab::Range ents;
01099         //     m_interface->get_entities_by_dimension(m_covering_source_set, 2, ents);
01100         //     m_interface->remove_entities(m_covering_source_set, ents);
01101         // }
01102     }
01103     else
01104     {
01105 #endif
01106         if( rrmgrids )
01107         {
01108             rval = m_interface->create_meshset( moab::MESHSET_SET, m_covering_source_set );MB_CHK_SET_ERR( rval, "Can't create new set" );
01109 
01110             double tolerance = 1e-6, btolerance = 1e-3;
01111             moab::AdaptiveKDTree tree( m_interface );
01112             moab::Range targetVerts;
01113 
01114             rval = m_interface->get_connectivity( m_target_entities, targetVerts, true );MB_CHK_ERR( rval );
01115 
01116             rval = tree.build_tree( m_source_entities, &m_source_set );MB_CHK_ERR( rval );
01117 
01118             for( unsigned ie = 0; ie < targetVerts.size(); ++ie )
01119             {
01120                 EntityHandle el = targetVerts[ie], leaf;
01121                 double point[3];
01122 
01123                 // Get the element centroid to be queried
01124                 rval = m_interface->get_coords( &el, 1, point );MB_CHK_ERR( rval );
01125 
01126                 // Search for the closest source element in the master mesh corresponding
01127                 // to the target element centroid in the slave mesh
01128                 rval = tree.point_search( point, leaf, tolerance, btolerance );MB_CHK_ERR( rval );
01129 
01130                 if( leaf == 0 )
01131                 {
01132                     leaf = m_source_set;  // no hint
01133                 }
01134 
01135                 std::vector< moab::EntityHandle > leaf_elems;
01136                 // We only care about the dimension that the user specified.
01137                 // MOAB partitions are ordered by elements anyway.
01138                 rval = m_interface->get_entities_by_dimension( leaf, 2, leaf_elems );MB_CHK_ERR( rval );
01139 
01140                 if( !leaf_elems.size() )
01141                 {
01142                     // std::cout << ie << ": " << " No leaf elements found." << std::endl;
01143                     continue;
01144                 }
01145 
01146                 // Now get the master element centroids so that we can compute
01147                 // the minimum distance to the target point
01148                 std::vector< double > centroids( leaf_elems.size() * 3 );
01149                 rval = m_interface->get_coords( &leaf_elems[0], leaf_elems.size(), &centroids[0] );MB_CHK_ERR( rval );
01150 
01151                 double dist = 1e5;
01152                 int pinelem = -1;
01153                 for( size_t il = 0; il < leaf_elems.size(); ++il )
01154                 {
01155                     const double* centroid = &centroids[il * 3];
01156                     const double locdist   = std::pow( point[0] - centroid[0], 2 ) +
01157                                            std::pow( point[1] - centroid[1], 2 ) +
01158                                            std::pow( point[2] - centroid[2], 2 );
01159 
01160                     if( locdist < dist )
01161                     {
01162                         dist    = locdist;
01163                         pinelem = il;
01164                         m_covering_source_entities.insert( leaf_elems[il] );
01165                     }
01166                 }
01167 
01168                 if( pinelem < 0 )
01169                 {
01170                     std::cout << ie
01171                               << ": [Error] - Could not find a minimum distance within the leaf "
01172                                  "nodes. Dist = "
01173                               << dist << std::endl;
01174                 }
01175             }
01176             // rval = tree.reset_tree();MB_CHK_ERR(rval);
01177             std::cout << "[INFO] - Total covering source entities = " << m_covering_source_entities.size() << std::endl;
01178             rval = m_interface->add_entities( m_covering_source_set, m_covering_source_entities );MB_CHK_ERR( rval );
01179         }
01180         else
01181         {
01182             m_covering_source_set      = m_source_set;
01183             m_covering_source          = m_source;
01184             m_covering_source_entities = m_source_entities;  // this is a tempest mesh object; careful about
01185                                                              // incrementing the reference?
01186             m_covering_source_vertices = m_source_vertices;  // this is a tempest mesh object; careful about
01187                                                              // incrementing the reference?
01188         }
01189 #ifdef MOAB_HAVE_MPI
01190     }
01191 #endif
01192 
01193     return rval;
01194 }
01195 
01196 ErrorCode TempestRemapper::ComputeOverlapMesh( bool kdtree_search, bool use_tempest )
01197 {
01198     ErrorCode rval;
01199     const bool outputEnabled = ( this->rank == 0 );
01200     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
01201     dbgprint.set_prefix( "[ComputeOverlapMesh]: " );
01202 
01203     // const double radius = 1.0 /*2.0*acos(-1.0)*/;
01204     // const double boxeps = 0.1;
01205     // Create the intersection on the sphere object and set up necessary parameters
01206 
01207     // First, split based on whether to use Tempest or MOAB
01208     // If Tempest
01209     //   1) Check for valid Mesh and pointers to objects for source/target
01210     //   2) Invoke GenerateOverlapWithMeshes routine from Tempest library
01211     // If MOAB
01212     //   1) Check for valid source and target meshsets (and entities)
01213     //   2) Build processor bounding boxes and construct a covering set
01214     //   3) Perform intersection between the source (covering) and target entities
01215     if( use_tempest )
01216     {
01217         // Now let us construct the overlap mesh, by calling TempestRemap interface directly
01218         // For the overlap method, choose between: "fuzzy", "exact" or "mixed"
01219         assert( m_source != NULL );
01220         assert( m_target != NULL );
01221         if( m_overlap != NULL ) delete m_overlap;
01222         m_overlap         = new Mesh();
01223         bool concaveMeshA = false, concaveMeshB = false;
01224         int err = GenerateOverlapWithMeshes( *m_covering_source, *m_target, *m_overlap, "" /*outFilename*/, "Netcdf4",
01225                                              "exact", concaveMeshA, concaveMeshB, false );
01226         if( err )
01227         {
01228             MB_CHK_SET_ERR( MB_FAILURE, "TempestRemap: Can't compute the intersection of meshes on the sphere" );
01229         }
01230     }
01231     else
01232     {
01233         Tag gidtag = m_interface->globalId_tag();
01234         moab::EntityHandle subrange[2];
01235         int gid[2];
01236         if( m_source_entities.size() > 1 )
01237         {  // Let us do some sanity checking to fix ID if they have are setup incorrectly
01238             subrange[0] = m_source_entities[0];
01239             subrange[1] = m_source_entities[1];
01240             rval        = m_interface->tag_get_data( gidtag, subrange, 2, gid );MB_CHK_ERR( rval );
01241 
01242             // Check if we need to impose Global ID numbering for vertices and elements. This may be
01243             // needed if we load the meshes from exodus or some other formats that may not have a
01244             // numbering forced.
01245             if( gid[0] + gid[1] == 0 )  // this implies first two elements have GID = 0
01246             {
01247 #ifdef MOAB_HAVE_MPI
01248                 rval = m_pcomm->assign_global_ids( m_source_set, 2, 1, false, true, false );MB_CHK_ERR( rval );
01249 #else
01250                 rval = this->assign_vertex_element_IDs( gidtag, m_source_set, 2, 1 );MB_CHK_ERR( rval );
01251 #endif
01252             }
01253         }
01254         if( m_target_entities.size() > 1 )
01255         {
01256             subrange[0] = m_target_entities[0];
01257             subrange[1] = m_target_entities[1];
01258             rval        = m_interface->tag_get_data( gidtag, subrange, 2, gid );MB_CHK_ERR( rval );
01259 
01260             // Check if we need to impose Global ID numbering for vertices and elements. This may be
01261             // needed if we load the meshes from exodus or some other formats that may not have a
01262             // numbering forced.
01263             if( gid[0] + gid[1] == 0 )  // this implies first two elements have GID = 0
01264             {
01265 #ifdef MOAB_HAVE_MPI
01266                 rval = m_pcomm->assign_global_ids( m_target_set, 2, 1, false, true, false );MB_CHK_ERR( rval );
01267 #else
01268                 rval = this->assign_vertex_element_IDs( gidtag, m_target_set, 2, 1 );MB_CHK_ERR( rval );
01269 #endif
01270             }
01271         }
01272 
01273         // Now perform the actual parallel intersection between the source and the target meshes
01274         if( kdtree_search )
01275         {
01276             if( outputEnabled ) dbgprint.printf( 0, "Computing intersection mesh with the Kd-tree search algorithm" );
01277             rval = mbintx->intersect_meshes_kdtree( m_covering_source_set, m_target_set, m_overlap_set );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere with brute-force" );
01278         }
01279         else
01280         {
01281             if( outputEnabled )
01282                 dbgprint.printf( 0, "Computing intersection mesh with the advancing-front propagation algorithm" );
01283             rval = mbintx->intersect_meshes( m_covering_source_set, m_target_set, m_overlap_set );MB_CHK_SET_ERR( rval, "Can't compute the intersection of meshes on the sphere" );
01284         }
01285 
01286 #ifdef MOAB_HAVE_MPI
01287         if( is_parallel || rrmgrids )
01288         {
01289 #ifdef VERBOSE
01290 
01291             std::stringstream ffc, fft, ffo;
01292             ffc << "cover_" << size << "_" << rank << ".h5m";
01293             fft << "target_" << size << "_" << rank << ".h5m";
01294             ffo << "intx_" << size << "_" << rank << ".h5m";
01295             rval = m_interface->write_mesh( ffc.str().c_str(), &m_covering_source_set, 1 );MB_CHK_ERR( rval );
01296             rval = m_interface->write_mesh( fft.str().c_str(), &m_target_set, 1 );MB_CHK_ERR( rval );
01297             rval = m_interface->write_mesh( ffo.str().c_str(), &m_overlap_set, 1 );MB_CHK_ERR( rval );
01298 
01299             // write the intx mesh, only in serial
01300             std::ostringstream opts;
01301             opts << "PARALLEL=WRITE_PART";
01302             std::ostringstream file_name;
01303             file_name << "intx_gl_" << size << ".h5m";
01304             if( size == 1 )
01305             {
01306                 rval = m_interface->write_file( file_name.str().c_str(), 0, opts.str().c_str(), &m_overlap_set, 1 );MB_CHK_ERR( rval );
01307             }
01308 
01309 #endif
01310             // because we do not want to work with elements in coverage set that do not participate
01311             // in intersection, remove them from the coverage set we will not delete them yet, just
01312             // remove from the set !
01313             if( !point_cloud_target )
01314             {
01315                 Range covEnts;
01316                 rval = m_interface->get_entities_by_dimension( m_covering_source_set, 2, covEnts );MB_CHK_ERR( rval );
01317                 Tag gidtag = m_interface->globalId_tag();
01318 
01319                 std::map< int, int > loc_gid_to_lid_covsrc;
01320                 std::vector< int > gids( covEnts.size(), -1 );
01321                 rval = m_interface->tag_get_data( gidtag, covEnts, &gids[0] );MB_CHK_ERR( rval );
01322                 for( unsigned ie = 0; ie < gids.size(); ++ie )
01323                 {
01324                     loc_gid_to_lid_covsrc[gids[ie]] = ie;
01325                 }
01326 
01327                 std::set< EntityHandle > intxCov;
01328                 Range intxCells;
01329                 Tag srcParentTag;
01330                 rval = m_interface->tag_get_handle( "SourceParent", srcParentTag );MB_CHK_ERR( rval );
01331                 rval = m_interface->get_entities_by_dimension( m_overlap_set, 2, intxCells );MB_CHK_ERR( rval );
01332                 for( Range::iterator it = intxCells.begin(); it != intxCells.end(); it++ )
01333                 {
01334                     EntityHandle intxCell = *it;
01335                     int blueParent        = -1;
01336                     rval                  = m_interface->tag_get_data( srcParentTag, &intxCell, 1, &blueParent );MB_CHK_ERR( rval );
01337                     // if (is_root) std::cout << "Found intersecting element: " << blueParent << ",
01338                     // " << gid_to_lid_covsrc[blueParent] << "\n";
01339                     assert( blueParent >= 0 );
01340                     intxCov.insert( covEnts[loc_gid_to_lid_covsrc[blueParent]] );
01341                 }
01342 
01343                 Range intxCovRange;
01344                 std::copy( intxCov.rbegin(), intxCov.rend(), range_inserter( intxCovRange ) );
01345                 Range notNeededCovCells = moab::subtract( covEnts, intxCovRange );
01346 
01347                 rval = m_interface->remove_entities( m_covering_source_set, notNeededCovCells );MB_CHK_ERR( rval );
01348                 covEnts = moab::subtract( covEnts, notNeededCovCells );
01349 #ifdef VERBOSE
01350                 std::cout << " total participating elements in the covering set: " << intxCov.size() << "\n";
01351                 std::cout << " remove from coverage set elements that are not intersected: " << notNeededCovCells.size()
01352                           << "\n";
01353 #endif
01354                 if( size > 1 )
01355                 {
01356                     // some source elements cover multiple target partitions; the conservation logic
01357                     // requires to know all overlap elements for a source element; they need to be
01358                     // communicated from the other target partitions
01359                     //
01360                     // so first we have to identify source (coverage) elements that cover multiple
01361                     // target partitions
01362 
01363                     // we will then mark the source, we will need to migrate the overlap elements
01364                     // that cover this to the original source for the source element; then
01365                     // distribute the overlap elements to all processors that have the coverage mesh
01366                     // used
01367 
01368                     rval = augment_overlap_set();MB_CHK_ERR( rval );
01369                 }
01370             }
01371 
01372             // m_covering_source = new Mesh();
01373             // rval = convert_mesh_to_tempest_private ( m_covering_source, m_covering_source_set,
01374             // m_covering_source_entities, &m_covering_source_vertices ); MB_CHK_SET_ERR ( rval,
01375             // "Can't convert source Tempest mesh" );
01376         }
01377 #endif
01378 
01379         // Fix any inconsistencies in the overlap mesh
01380         {
01381             IntxAreaUtils areaAdaptor;
01382             rval = IntxUtils::fix_degenerate_quads( m_interface, m_overlap_set );MB_CHK_ERR( rval );
01383             rval = areaAdaptor.positive_orientation( m_interface, m_overlap_set, 1.0 /*radius*/, rank );MB_CHK_ERR( rval );
01384         }
01385 
01386         // Now let us re-convert the MOAB mesh back to Tempest representation
01387         rval = this->ComputeGlobalLocalMaps();MB_CHK_ERR( rval );
01388 
01389         rval = this->convert_overlap_mesh_sorted_by_source();MB_CHK_ERR( rval );
01390         // free the memory
01391         delete mbintx;
01392     }
01393 
01394     return MB_SUCCESS;
01395 }
01396 
01397 #ifdef MOAB_HAVE_MPI
01398 // this function is called only in parallel
01399 ///////////////////////////////////////////////////////////////////////////////////
01400 ErrorCode TempestRemapper::augment_overlap_set()
01401 {
01402     /*
01403      * overall strategy:
01404      *
01405      * 1) collect all boundary target cells on the current task, affected by the partition boundary;
01406      *    note: not only partition boundary, we need all boundary (all coastal lines) and partition
01407      * boundary targetBoundaryIds is the set of target boundary cell IDs
01408      *
01409      * 2) collect all source cells that are intersecting boundary cells (call them
01410      * affectedSourceCellsIds)
01411      *
01412      * 3) collect overlap, that is accumulate all overlap cells that have source target in
01413      * affectedSourceCellsIds
01414      */
01415     // first, get all edges on the partition boundary, on the target mesh, then all the target
01416     // elements that border the partition boundary
01417     ErrorCode rval;
01418     Skinner skinner( m_interface );
01419     Range targetCells, boundaryEdges;
01420     rval = m_interface->get_entities_by_dimension( m_target_set, 2, targetCells );MB_CHK_ERR( rval );
01421     /// find all boundary edges
01422     rval = skinner.find_skin( 0, targetCells, false, boundaryEdges );MB_CHK_ERR( rval );
01423     // filter boundary edges that are on partition boundary, not on boundary
01424     // find all cells adjacent to these boundary edges, from target set
01425     Range boundaryCells;  // these will be filtered from target_set
01426     rval = m_interface->get_adjacencies( boundaryEdges, 2, false, boundaryCells, Interface::UNION );MB_CHK_ERR( rval );
01427     boundaryCells = intersect( boundaryCells, targetCells );
01428 #ifdef VERBOSE
01429     EntityHandle tmpSet;
01430     rval = m_interface->create_meshset( MESHSET_SET, tmpSet );MB_CHK_SET_ERR( rval, "Can't create temporary set" );
01431     // add the boundary set and edges, and save it to a file
01432     rval = m_interface->add_entities( tmpSet, boundaryCells );MB_CHK_SET_ERR( rval, "Can't add entities" );
01433     rval = m_interface->add_entities( tmpSet, boundaryEdges );MB_CHK_SET_ERR( rval, "Can't add edges" );
01434     std::stringstream ffs;
01435     ffs << "boundaryCells_0" << rank << ".h5m";
01436     rval = m_interface->write_mesh( ffs.str().c_str(), &tmpSet, 1 );MB_CHK_ERR( rval );
01437 #endif
01438 
01439     // now that we have the boundary cells, see which overlap polys have have these as parents;
01440     //   find the ids of the boundary cells;
01441     Tag gid = m_interface->globalId_tag();
01442     std::set< int > targetBoundaryIds;
01443     for( Range::iterator it = boundaryCells.begin(); it != boundaryCells.end(); it++ )
01444     {
01445         int tid;
01446         EntityHandle targetCell = *it;
01447         rval                    = m_interface->tag_get_data( gid, &targetCell, 1, &tid );MB_CHK_SET_ERR( rval, "Can't get global id tag on target cell" );
01448         if( tid < 0 ) std::cout << " incorrect id for a target cell\n";
01449         targetBoundaryIds.insert( tid );
01450     }
01451 
01452     Range overlapCells;
01453     rval = m_interface->get_entities_by_dimension( m_overlap_set, 2, overlapCells );MB_CHK_ERR( rval );
01454 
01455     std::set< int > affectedSourceCellsIds;
01456     Tag targetParentTag, sourceParentTag;  // do not use blue/red, as it is more confusing
01457     rval = m_interface->tag_get_handle( "TargetParent", targetParentTag );MB_CHK_ERR( rval );
01458     rval = m_interface->tag_get_handle( "SourceParent", sourceParentTag );MB_CHK_ERR( rval );
01459     for( Range::iterator it = overlapCells.begin(); it != overlapCells.end(); it++ )
01460     {
01461         EntityHandle intxCell = *it;
01462         int targetParentID, sourceParentID;
01463         rval = m_interface->tag_get_data( targetParentTag, &intxCell, 1, &targetParentID );MB_CHK_ERR( rval );
01464         if( targetBoundaryIds.find( targetParentID ) != targetBoundaryIds.end() )
01465         {
01466             // this means that the source element is affected
01467             rval = m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID );MB_CHK_ERR( rval );
01468             affectedSourceCellsIds.insert( sourceParentID );
01469         }
01470     }
01471 
01472     // now find all source cells affected, based on their id;
01473     //  (we do not have yet the mapping gid_to_lid_covsrc)
01474     std::map< int, EntityHandle > affectedCovCellFromID;  // map from source cell id to the eh; it is needed to find out
01475                                                           // the original processor
01476     // this one came from , so to know where to send the overlap elements
01477 
01478     // use std::set<EntityHandle> instead of moab::Range for collecting cells, either on coverage or
01479     // target or intx cells
01480     std::set< EntityHandle > affectedCovCells;  // their overlap cells will be sent to their
01481                                                 // original task, then distributed to all
01482     // other processes that might need them to compute conservation
01483 
01484     Range covCells;
01485     rval = m_interface->get_entities_by_dimension( m_covering_source_set, 2, covCells );MB_CHK_ERR( rval );
01486     // loop thru all cov cells, to find the ones with global ids in affectedSourceCellsIds
01487     for( Range::iterator it = covCells.begin(); it != covCells.end(); it++ )
01488     {
01489         EntityHandle covCell = *it;  //
01490         int covID;
01491         rval = m_interface->tag_get_data( gid, &covCell, 1, &covID );
01492         if( affectedSourceCellsIds.find( covID ) != affectedSourceCellsIds.end() )
01493         {
01494             // this source cell is affected;
01495             affectedCovCellFromID[covID] = covCell;
01496             affectedCovCells.insert( covCell );
01497         }
01498     }
01499 
01500     // now loop again over all overlap cells, to see if their source parent is "affected"
01501     // store in ranges the overlap cells that need to be sent to original task of the source cell
01502     // from there, they will be redistributed to the tasks that need that coverage cell
01503     Tag sendProcTag;
01504     rval = m_interface->tag_get_handle( "sending_processor", 1, MB_TYPE_INTEGER, sendProcTag );
01505 
01506     // basically a map from original processor task to the set of overlap cells to be sent there
01507     std::map< int, std::set< EntityHandle > > overlapCellsForTask;
01508     // this set will contain all intx cells that will need to be sent ( a union of above sets ,
01509     //   that are organized per task on the above map )
01510     std::set< EntityHandle > overlapCellsToSend;
01511 
01512     for( Range::iterator it = overlapCells.begin(); it != overlapCells.end(); it++ )
01513     {
01514         EntityHandle intxCell = *it;
01515         int sourceParentID;
01516         rval = m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID );MB_CHK_ERR( rval );
01517         if( affectedSourceCellsIds.find( sourceParentID ) != affectedSourceCellsIds.end() )
01518         {
01519             EntityHandle covCell = affectedCovCellFromID[sourceParentID];
01520             int orgTask;
01521             rval = m_interface->tag_get_data( sendProcTag, &covCell, 1, &orgTask );MB_CHK_ERR( rval );
01522             overlapCellsForTask[orgTask].insert(
01523                 intxCell );                         // put the overlap cell in corresponding range (set<EntityHandle>)
01524             overlapCellsToSend.insert( intxCell );  // also put it in this range, for debugging mostly
01525         }
01526     }
01527 
01528     // now prepare to send; will use crystal router, as the buffers in ParComm are prepared only
01529     // for neighbors; coverage mesh was also migrated with crystal router, so here we go again :(
01530 
01531     // find out the maximum number of edges of the polygons needed to be sent
01532     // we could we conservative and use a big number, or the number from intx, if we store it then?
01533     int maxEdges = 0;
01534     for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); it++ )
01535     {
01536         EntityHandle intxCell = *it;
01537         int nnodes;
01538         const EntityHandle* conn;
01539         rval = m_interface->get_connectivity( intxCell, conn, nnodes );MB_CHK_ERR( rval );
01540         if( maxEdges < nnodes ) maxEdges = nnodes;
01541     }
01542 
01543     // find the maximum among processes in intersection
01544     int globalMaxEdges;
01545     if( m_pcomm )
01546         MPI_Allreduce( &maxEdges, &globalMaxEdges, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
01547     else
01548         globalMaxEdges = maxEdges;
01549 
01550     //#ifdef VERBOSE
01551     if( is_root ) std::cout << "maximum number of edges for polygons to send is " << globalMaxEdges << "\n";
01552         //#endif
01553 
01554 #ifdef VERBOSE
01555     EntityHandle tmpSet2;
01556     rval = m_interface->create_meshset( MESHSET_SET, tmpSet2 );MB_CHK_SET_ERR( rval, "Can't create temporary set2" );
01557     // add the affected source and overlap elements
01558     for( std::set< EntityHandle >::iterator it = overlapCellsToSend.begin(); it != overlapCellsToSend.end(); it++ )
01559     {
01560         EntityHandle intxCell = *it;
01561         rval                  = m_interface->add_entities( tmpSet2, &intxCell, 1 );MB_CHK_SET_ERR( rval, "Can't add entities" );
01562     }
01563     for( std::set< EntityHandle >::iterator it = affectedCovCells.begin(); it != affectedCovCells.end(); it++ )
01564     {
01565         EntityHandle covCell = *it;
01566         rval                 = m_interface->add_entities( tmpSet2, &covCell, 1 );MB_CHK_SET_ERR( rval, "Can't add entities" );
01567     }
01568     std::stringstream ffs2;
01569     // these will contain coverage cells and intx cells on the boundary
01570     ffs2 << "affectedCells_" << m_pcomm->rank() << ".h5m";
01571     rval = m_interface->write_mesh( ffs2.str().c_str(), &tmpSet2, 1 );MB_CHK_ERR( rval );
01572 #endif
01573     // form tuple lists to send vertices and cells;
01574     // the problem is that the lists of vertices will need to have other information, like the
01575     // processor it comes from, and its index in that list; we may have to duplicate vertices, but
01576     // we do not care much; we will not duplicate overlap elements, just the vertices, as they may
01577     // come from different cells and different processes each vertex will have a local index and a
01578     // processor task it is coming from
01579 
01580     // look through the std::set's to be sent to other processes, and form the vertex tuples and
01581     // cell tuples
01582     //
01583     std::map< int, std::set< EntityHandle > > verticesOverlapForTask;
01584     // Range allVerticesToSend;
01585     std::set< EntityHandle > allVerticesToSend;
01586     std::map< EntityHandle, int > allVerticesToSendMap;
01587     int numVerts        = 0;
01588     int numOverlapCells = 0;
01589     for( std::map< int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
01590          it != overlapCellsForTask.end(); it++ )
01591     {
01592         int sendToProc                                = it->first;
01593         std::set< EntityHandle >& overlapCellsToSend2 = it->second;  // organize vertices in std::set per processor
01594         // Range vertices;
01595         std::set< EntityHandle > vertices;  // collect all vertices connected to overlapCellsToSend2
01596         for( std::set< EntityHandle >::iterator set_it = overlapCellsToSend2.begin();
01597              set_it != overlapCellsToSend2.end(); ++set_it )
01598         {
01599             int nnodes_local          = 0;
01600             const EntityHandle* conn1 = NULL;
01601             rval                      = m_interface->get_connectivity( *set_it, conn1, nnodes_local );MB_CHK_ERR( rval );
01602             for( int k = 0; k < nnodes_local; k++ )
01603                 vertices.insert( conn1[k] );
01604         }
01605         verticesOverlapForTask[sendToProc] = vertices;
01606         numVerts += (int)vertices.size();
01607         numOverlapCells += (int)overlapCellsToSend2.size();
01608         allVerticesToSend.insert( vertices.begin(), vertices.end() );
01609     }
01610     // build the index map, from entity handle to index in all vert set
01611     int j = 0;
01612     for( std::set< EntityHandle >::iterator vert_it = allVerticesToSend.begin(); vert_it != allVerticesToSend.end();
01613          vert_it++, j++ )
01614     {
01615         EntityHandle vert          = *vert_it;
01616         allVerticesToSendMap[vert] = j;
01617     }
01618 
01619     // first send vertices in a tuple list, then send overlap cells, according to requests
01620     // overlap cells need to send info about the blue and red parent tags, too
01621     TupleList TLv;                           //
01622     TLv.initialize( 2, 0, 0, 3, numVerts );  // to proc, index in all range, DP points
01623     TLv.enableWriteAccess();
01624 
01625     for( std::map< int, std::set< EntityHandle > >::iterator it = verticesOverlapForTask.begin();
01626          it != verticesOverlapForTask.end(); it++ )
01627     {
01628         int sendToProc                     = it->first;
01629         std::set< EntityHandle >& vertices = it->second;
01630         int i                              = 0;
01631         for( std::set< EntityHandle >::iterator it2 = vertices.begin(); it2 != vertices.end(); it2++, i++ )
01632         {
01633             int n                = TLv.get_n();
01634             TLv.vi_wr[2 * n]     = sendToProc;  // send to processor
01635             EntityHandle v       = *it2;
01636             int indexInAllVert   = allVerticesToSendMap[v];
01637             TLv.vi_wr[2 * n + 1] = indexInAllVert;  // will be orgProc, to differentiate indices
01638                                                     // of vertices sent to "sentToProc"
01639             double coords[3];
01640             rval = m_interface->get_coords( &v, 1, coords );MB_CHK_ERR( rval );
01641             TLv.vr_wr[3 * n]     = coords[0];  // departure position, of the node local_verts[i]
01642             TLv.vr_wr[3 * n + 1] = coords[1];
01643             TLv.vr_wr[3 * n + 2] = coords[2];
01644             TLv.inc_n();
01645         }
01646     }
01647 
01648     TupleList TLc;
01649     int sizeTuple = 4 + globalMaxEdges;
01650     // total number of overlap cells to send
01651     TLc.initialize( sizeTuple, 0, 0, 0,
01652                     numOverlapCells );  // to proc, blue parent ID, red parent ID, nvert,
01653                                         // connectivity[globalMaxEdges] (global ID v), local eh)
01654     TLc.enableWriteAccess();
01655 
01656     for( std::map< int, std::set< EntityHandle > >::iterator it = overlapCellsForTask.begin();
01657          it != overlapCellsForTask.end(); it++ )
01658     {
01659         int sendToProc                                = it->first;
01660         std::set< EntityHandle >& overlapCellsToSend2 = it->second;
01661         // send also the target and source parents for these overlap cells
01662         for( std::set< EntityHandle >::iterator it2 = overlapCellsToSend2.begin(); it2 != overlapCellsToSend2.end();
01663              it2++ )
01664         {
01665             EntityHandle intxCell = *it2;
01666             int sourceParentID, targetParentID;
01667             rval = m_interface->tag_get_data( targetParentTag, &intxCell, 1, &targetParentID );MB_CHK_ERR( rval );
01668             rval = m_interface->tag_get_data( sourceParentTag, &intxCell, 1, &sourceParentID );MB_CHK_ERR( rval );
01669             int n                        = TLc.get_n();
01670             TLc.vi_wr[sizeTuple * n]     = sendToProc;
01671             TLc.vi_wr[sizeTuple * n + 1] = sourceParentID;
01672             TLc.vi_wr[sizeTuple * n + 2] = targetParentID;
01673             int nnodes;
01674             const EntityHandle* conn = NULL;
01675             rval                     = m_interface->get_connectivity( intxCell, conn, nnodes );MB_CHK_ERR( rval );
01676             TLc.vi_wr[sizeTuple * n + 3] = nnodes;
01677             for( int i = 0; i < nnodes; i++ )
01678             {
01679                 int indexVertex = allVerticesToSendMap[conn[i]];
01680                 ;  // the vertex index will be now unique per original proc
01681                 if( -1 == indexVertex ) MB_CHK_SET_ERR( MB_FAILURE, "Can't find vertex in range of vertices to send" );
01682                 TLc.vi_wr[sizeTuple * n + 4 + i] = indexVertex;
01683             }
01684             // fill the rest with 0, just because we do not like uninitialized data
01685             for( int i = nnodes; i < globalMaxEdges; i++ )
01686                 TLc.vi_wr[sizeTuple * n + 4 + i] = 0;
01687 
01688             TLc.inc_n();
01689         }
01690     }
01691 
01692     // send first the vertices and overlap cells to original task for coverage cells
01693     // now we are done populating the tuples; route them to the appropriate processors
01694 #ifdef VERBOSE
01695     std::stringstream ff1;
01696     ff1 << "TLc_" << rank << ".txt";
01697     TLc.print_to_file( ff1.str().c_str() );
01698     std::stringstream ffv;
01699     ffv << "TLv_" << rank << ".txt";
01700     TLv.print_to_file( ffv.str().c_str() );
01701 #endif
01702     ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv, 0 );
01703     ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc, 0 );
01704 
01705 #ifdef VERBOSE
01706     TLc.print_to_file( ff1.str().c_str() );  // will append to existing file
01707     TLv.print_to_file( ffv.str().c_str() );
01708 #endif
01709     // first phase of transfer complete
01710     // now look at TLc, and sort by the source parent (index 1)
01711 
01712     TupleList::buffer buffer;
01713     buffer.buffer_init( sizeTuple * TLc.get_n() * 2 );  // allocate memory for sorting !! double
01714     TLc.sort( 1, &buffer );
01715 #ifdef VERBOSE
01716     TLc.print_to_file( ff1.str().c_str() );
01717 #endif
01718 
01719     // will keep a map with vertices per processor that will need to be used in TLv2;
01720     // so, availVertexIndicesPerProcessor[proc] is a map from vertex indices that are available from
01721     // this processor to the index in the local TLv; the used vertices will have to be sent to the
01722     // tasks that need them
01723 
01724     // connectivity of a cell is given by sending proc and index in original list of vertices from
01725     // that proc
01726 
01727     std::map< int, std::map< int, int > > availVertexIndicesPerProcessor;
01728     int nv = TLv.get_n();
01729     for( int i = 0; i < nv; i++ )
01730     {
01731         // int proc=TLv.vi_rd[3*i]; // it is coming from this processor
01732         int orgProc   = TLv.vi_rd[2 * i];  // this is the original processor, for index vertex consideration
01733         int indexVert = TLv.vi_rd[2 * i + 1];
01734         availVertexIndicesPerProcessor[orgProc][indexVert] = i;
01735     }
01736 
01737     // now we have sorted the incoming overlap elements by the source element;
01738     // if we have overlap elements for one source coming from 2 or more processes, we need to send
01739     // back to the processes that do not have that overlap cell;
01740 
01741     // form new TLc2, TLv2, that will be distributed to necessary processes
01742     // first count source elements that are "spread" over multiple processes
01743     // TLc is ordered now by source ID; loop over them
01744     int n = TLc.get_n();  // total number of overlap elements received on current task;
01745     std::map< int, int > currentProcsCount;
01746     // form a map from proc to sets of vertex indices that will be sent using TLv2
01747     // will form a map between a source cell ID and tasks/targets that are partially overlapped by
01748     // these sources
01749     std::map< int, std::set< int > > sourcesForTasks;
01750     int sizeOfTLc2 = 0;  // only increase when we will have to send data
01751     if( n > 0 )
01752     {
01753         int currentSourceID      = TLc.vi_rd[sizeTuple * 0 + 1];  // we  have written sizeTuple*0 for "clarity"
01754         int proc0                = TLc.vi_rd[sizeTuple * 0];
01755         currentProcsCount[proc0] = 1;  //
01756 
01757         for( int i = 1; i < n; i++ )
01758         {
01759             int proc     = TLc.vi_rd[sizeTuple * i];
01760             int sourceID = TLc.vi_rd[sizeTuple * i + 1];
01761             if( sourceID == currentSourceID )
01762             {
01763                 if( currentProcsCount.find( proc ) == currentProcsCount.end() )
01764                 {
01765                     currentProcsCount[proc] = 1;
01766                 }
01767                 else
01768                     currentProcsCount[proc]++;
01769             }
01770             if( sourceID != currentSourceID || ( ( n - 1 ) == i ) )  // we study the current source if we reach the last
01771             {
01772                 // we have found a new source id, need to reset the proc counts, and establish if we
01773                 // need to send data
01774                 if( currentProcsCount.size() > 1 )
01775                 {
01776 #ifdef VERBOSE
01777                     std::cout << " source element " << currentSourceID << " intersects with "
01778                               << currentProcsCount.size() << " target partitions\n";
01779                     for( std::map< int, int >::iterator it = currentProcsCount.begin(); it != currentProcsCount.end();
01780                          it++ )
01781                     {
01782                         int procID       = it->first;
01783                         int numOverCells = it->second;
01784                         std::cout << "   task:" << procID << " " << numOverCells << " cells\n";
01785                     }
01786 
01787 #endif
01788                     // estimate what we need to send
01789                     for( std::map< int, int >::iterator it1 = currentProcsCount.begin(); it1 != currentProcsCount.end();
01790                          it1++ )
01791                     {
01792                         int proc1 = it1->first;
01793                         sourcesForTasks[currentSourceID].insert( proc1 );
01794                         for( std::map< int, int >::iterator it2 = currentProcsCount.begin();
01795                              it2 != currentProcsCount.end(); it2++ )
01796                         {
01797                             int proc2 = it2->first;
01798                             if( proc1 != proc2 ) sizeOfTLc2 += it2->second;
01799                         }
01800                     }
01801                     // mark vertices in TLv tuple that need to be sent
01802                 }
01803                 if( sourceID != currentSourceID )  // maybe we are not at the end, so continue on
01804                 {
01805                     currentSourceID = sourceID;
01806                     currentProcsCount.clear();
01807                     currentProcsCount[proc] = 1;
01808                 }
01809             }
01810         }
01811     }
01812     // begin a loop to send the needed cells to the processes; also mark the vertices that need to
01813     // be sent, put them in a set
01814 
01815 #ifdef VERBOSE
01816     std::cout << " need to initialize TLc2 with " << sizeOfTLc2 << " cells\n ";
01817 #endif
01818 
01819     TupleList TLc2;
01820     int sizeTuple2 = 5 + globalMaxEdges;  // send to, original proc for intx cell, source parent id,
01821                                           // target parent id,
01822     // number of vertices, then connectivity in terms of indices in vertex lists from original proc
01823     TLc2.initialize( sizeTuple2, 0, 0, 0, sizeOfTLc2 );
01824     TLc2.enableWriteAccess();
01825     // loop again through TLc, and select intx cells that have the problem sources;
01826 
01827     std::map< int, std::set< int > > verticesToSendForProc;  // will look at indices in the TLv list
01828     // will form for each processor, the index list from TLv
01829     for( int i = 0; i < n; i++ )
01830     {
01831         int sourceID = TLc.vi_rd[sizeTuple * i + 1];
01832         if( sourcesForTasks.find( sourceID ) != sourcesForTasks.end() )
01833         {
01834             // it means this intx cell needs to be sent to any proc that is not "original" to it
01835             std::set< int > procs = sourcesForTasks[sourceID];  // set of processors involved with this source
01836             if( procs.size() < 2 ) MB_CHK_SET_ERR( MB_FAILURE, " not enough processes involved with a sourceID cell" );
01837 
01838             int orgProc = TLc.vi_rd[sizeTuple * i];  // this intx cell was sent from this orgProc, originally
01839             // will need to be sent to all other procs from above set; also, need to mark the vertex
01840             // indices for that proc, and check that they are available to populate TLv2
01841             std::map< int, int >& availableVerticesFromThisProc = availVertexIndicesPerProcessor[orgProc];
01842             for( std::set< int >::iterator setIt = procs.begin(); setIt != procs.end(); setIt++ )
01843             {
01844                 int procID = *setIt;
01845                 // send this cell to the other processors, not to orgProc this cell is coming from
01846 
01847                 if( procID != orgProc )
01848                 {
01849                     // send the cell to this processor;
01850                     int n2 = TLc2.get_n();
01851                     if( n2 >= sizeOfTLc2 ) MB_CHK_SET_ERR( MB_FAILURE, " memory overflow" );
01852                     //
01853                     std::set< int >& indexVerticesInTLv = verticesToSendForProc[procID];
01854                     TLc2.vi_wr[n2 * sizeTuple2]         = procID;                    // send to
01855                     TLc2.vi_wr[n2 * sizeTuple2 + 1]     = orgProc;                   // this cell is coming from here
01856                     TLc2.vi_wr[n2 * sizeTuple2 + 2]     = sourceID;                  // source parent of the intx cell
01857                     TLc2.vi_wr[n2 * sizeTuple2 + 3] = TLc.vi_rd[sizeTuple * i + 2];  // target parent of the intx cell
01858                         // number of vertices of the intx cell
01859                     int nvert                       = TLc.vi_rd[sizeTuple * i + 3];
01860                     TLc2.vi_wr[n2 * sizeTuple2 + 4] = nvert;
01861                     // now loop through the connectivity, and make sure the vertices are available;
01862                     // mark them, to populate later the TLv2 tuple list
01863 
01864                     // just copy the vertices, including 0 ones
01865                     for( int j = 0; j < nvert; j++ )
01866                     {
01867                         int vertexIndex = TLc.vi_rd[i * sizeTuple + 4 + j];
01868                         // is this vertex available from org proc?
01869                         if( availableVerticesFromThisProc.find( vertexIndex ) == availableVerticesFromThisProc.end() )
01870                         {
01871                             MB_CHK_SET_ERR( MB_FAILURE, " vertex index not available from processor" );
01872                         }
01873                         TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = vertexIndex;
01874                         int indexInTLv                      = availVertexIndicesPerProcessor[orgProc][vertexIndex];
01875                         indexVerticesInTLv.insert( indexInTLv );
01876                     }
01877 
01878                     for( int j = nvert; j < globalMaxEdges; j++ )
01879                     {
01880                         TLc2.vi_wr[n2 * sizeTuple2 + 5 + j] = 0;  // or mark them 0
01881                     }
01882                     TLc2.inc_n();
01883                 }
01884             }
01885         }
01886     }
01887 
01888     // now we have to populate TLv2, with original source proc, index of vertex, and coordinates
01889     // from TLv use the verticesToSendForProc sets from above, and the map from index in proc to the
01890     // index in TLv
01891     TupleList TLv2;
01892     int numVerts2 = 0;
01893     // how many vertices to send?
01894     for( std::map< int, std::set< int > >::iterator it = verticesToSendForProc.begin();
01895          it != verticesToSendForProc.end(); it++ )
01896     {
01897         std::set< int >& indexInTLvSet = it->second;
01898         numVerts2 += (int)indexInTLvSet.size();
01899     }
01900     TLv2.initialize( 3, 0, 0, 3,
01901                      numVerts2 );  // send to, original proc, index in original proc, and 3 coords
01902     TLv2.enableWriteAccess();
01903     for( std::map< int, std::set< int > >::iterator it = verticesToSendForProc.begin();
01904          it != verticesToSendForProc.end(); it++ )
01905     {
01906         int sendToProc                 = it->first;
01907         std::set< int >& indexInTLvSet = it->second;
01908         // now, look at indices in TLv, to find out the original proc, and the index in that list
01909         for( std::set< int >::iterator itSet = indexInTLvSet.begin(); itSet != indexInTLvSet.end(); itSet++ )
01910         {
01911             int indexInTLv           = *itSet;
01912             int orgProc              = TLv.vi_rd[2 * indexInTLv];
01913             int indexVertexInOrgProc = TLv.vi_rd[2 * indexInTLv + 1];
01914             int nv2                  = TLv2.get_n();
01915             TLv2.vi_wr[3 * nv2]      = sendToProc;
01916             TLv2.vi_wr[3 * nv2 + 1]  = orgProc;
01917             TLv2.vi_wr[3 * nv2 + 2]  = indexVertexInOrgProc;
01918             for( int j = 0; j < 3; j++ )
01919                 TLv2.vr_wr[3 * nv2 + j] =
01920                     TLv.vr_rd[3 * indexInTLv + j];  // departure position, of the node local_verts[i]
01921             TLv2.inc_n();
01922         }
01923     }
01924     // now, finally, transfer the vertices and the intx cells;
01925     ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLv2, 0 );
01926     ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, TLc2, 0 );
01927     // now, look at vertices from TLv2, and create them
01928     // we should have in TLv2 only vertices with orgProc different from current task
01929 #ifdef VERBOSE
01930     std::stringstream ff2;
01931     ff2 << "TLc2_" << rank << ".txt";
01932     TLc2.print_to_file( ff2.str().c_str() );
01933     std::stringstream ffv2;
01934     ffv2 << "TLv2_" << rank << ".txt";
01935     TLv2.print_to_file( ffv2.str().c_str() );
01936 #endif
01937     // first create vertices, and make a map from origin processor, and index, to entity handle
01938     // (index in TLv2 )
01939     Tag ghostTag;
01940     int orig_proc = -1;
01941     rval          = m_interface->tag_get_handle( "ORIG_PROC", 1, MB_TYPE_INTEGER, ghostTag, MB_TAG_DENSE | MB_TAG_CREAT,
01942                                         &orig_proc );MB_CHK_ERR( rval );
01943 
01944     int nvNew = TLv2.get_n();
01945     // if number of vertices to be created is 0, it means there is no need of ghost intx cells,
01946     // because everything matched perfectly (it can happen in manufactured cases)
01947     if( 0 == nvNew ) return MB_SUCCESS;
01948     // create a vertex h for each coordinate
01949     Range newVerts;
01950     rval = m_interface->create_vertices( &( TLv2.vr_rd[0] ), nvNew, newVerts );MB_CHK_ERR( rval );
01951     // now create a map from index , org proc, to actual entity handle corresponding to it
01952     std::map< int, std::map< int, EntityHandle > > vertexPerProcAndIndex;
01953     for( int i = 0; i < nvNew; i++ )
01954     {
01955         int orgProc                                 = TLv2.vi_rd[3 * i + 1];
01956         int indexInVert                             = TLv2.vi_rd[3 * i + 2];
01957         vertexPerProcAndIndex[orgProc][indexInVert] = newVerts[i];
01958     }
01959 
01960     // new polygons will receive a dense tag, with default value -1, with the processor task they
01961     // originally belonged to
01962 
01963     // now form the needed cells, in order
01964     Range newPolygons;
01965     int ne = TLc2.get_n();
01966     for( int i = 0; i < ne; i++ )
01967     {
01968         int orgProc  = TLc2.vi_rd[i * sizeTuple2 + 1];  // this cell is coming from here, originally
01969         int sourceID = TLc2.vi_rd[i * sizeTuple2 + 2];  // source parent of the intx cell
01970         int targetID = TLc2.vi_wr[i * sizeTuple2 + 3];  // target parent of intx cell
01971         int nve      = TLc2.vi_wr[i * sizeTuple2 + 4];  // number of vertices for the polygon
01972         std::vector< EntityHandle > conn;
01973         conn.resize( nve );
01974         for( int j = 0; j < nve; j++ )
01975         {
01976             int indexV      = TLc2.vi_wr[i * sizeTuple2 + 5 + j];
01977             EntityHandle vh = vertexPerProcAndIndex[orgProc][indexV];
01978             conn[j]         = vh;
01979         }
01980         EntityHandle polyNew;
01981         rval = m_interface->create_element( MBPOLYGON, &conn[0], nve, polyNew );MB_CHK_ERR( rval );
01982         newPolygons.insert( polyNew );
01983         rval = m_interface->tag_set_data( targetParentTag, &polyNew, 1, &targetID );MB_CHK_ERR( rval );
01984         rval = m_interface->tag_set_data( sourceParentTag, &polyNew, 1, &sourceID );MB_CHK_ERR( rval );
01985         rval = m_interface->tag_set_data( ghostTag, &polyNew, 1, &orgProc );MB_CHK_ERR( rval );
01986     }
01987 
01988 #ifdef VERBOSE
01989     EntityHandle tmpSet3;
01990     rval = m_interface->create_meshset( MESHSET_SET, tmpSet3 );MB_CHK_SET_ERR( rval, "Can't create temporary set3" );
01991     // add the boundary set and edges, and save it to a file
01992     rval = m_interface->add_entities( tmpSet3, newPolygons );MB_CHK_SET_ERR( rval, "Can't add entities" );
01993 
01994     std::stringstream ffs4;
01995     ffs4 << "extraIntxCells" << rank << ".h5m";
01996     rval = m_interface->write_mesh( ffs4.str().c_str(), &tmpSet3, 1 );MB_CHK_ERR( rval );
01997 #endif
01998 
01999     // add the new polygons to the overlap set
02000     // these will be ghosted, so will participate in conservation only
02001     rval = m_interface->add_entities( m_overlap_set, newPolygons );MB_CHK_ERR( rval );
02002 #ifdef VERBOSE
02003     if( !rank )
02004     {
02005         std::cout << "Augmenting: add " << newPolygons.size() << " polygons on root task \n";
02006     }
02007 #endif
02008     return MB_SUCCESS;
02009 }
02010 
02011 #endif
02012 
02013 ErrorCode TempestRemapper::GetIMasks( Remapper::IntersectionContext ctx, std::vector< int >& masks )
02014 {
02015     Tag maskTag;
02016     // it should have been created already, if not, we might have a problem
02017     int def_val = 1;
02018     ErrorCode rval =
02019         m_interface->tag_get_handle( "GRID_IMASK", 1, MB_TYPE_INTEGER, maskTag, MB_TAG_DENSE | MB_TAG_CREAT, &def_val );MB_CHK_SET_ERR( rval, "Trouble creating GRID_IMASK tag" );
02020 
02021     switch( ctx )
02022     {
02023         case Remapper::SourceMesh: {
02024             if( point_cloud_source )
02025             {
02026                 masks.resize( m_source_vertices.size() );
02027                 rval = m_interface->tag_get_data( maskTag, m_source_vertices, &masks[0] );MB_CHK_SET_ERR( rval, "Trouble getting GRID_IMASK tag" );
02028             }
02029             else
02030             {
02031                 masks.resize( m_source_entities.size() );
02032                 rval = m_interface->tag_get_data( maskTag, m_source_entities, &masks[0] );MB_CHK_SET_ERR( rval, "Trouble getting GRID_IMASK tag" );
02033             }
02034             return MB_SUCCESS;
02035         }
02036         case Remapper::TargetMesh: {
02037             if( point_cloud_target )
02038             {
02039                 masks.resize( m_target_vertices.size() );
02040                 rval = m_interface->tag_get_data( maskTag, m_target_vertices, &masks[0] );MB_CHK_SET_ERR( rval, "Trouble getting GRID_IMASK tag" );
02041             }
02042             else
02043             {
02044                 masks.resize( m_target_entities.size() );
02045                 rval = m_interface->tag_get_data( maskTag, m_target_entities, &masks[0] );MB_CHK_SET_ERR( rval, "Trouble getting GRID_IMASK tag" );
02046             }
02047             return MB_SUCCESS;
02048         }
02049         case Remapper::CoveringMesh:
02050         case Remapper::OverlapMesh:
02051         default:
02052             return MB_SUCCESS;
02053     }
02054 }
02055 
02056 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines