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