MOAB: Mesh Oriented datABase  (version 5.4.1)
TempestOnlineMap.cpp
Go to the documentation of this file.
00001 /*
00002  * =====================================================================================
00003  *
00004  *       Filename:  TempestOnlineMap.hpp
00005  *
00006  *    Description:  Interface to the TempestRemap library to compute the consistent,
00007  *                  and accurate high-order conservative remapping weights for overlap
00008  *                  grids on the sphere in climate simulations.
00009  *
00010  *         Author:  Vijay S. Mahadevan (vijaysm), [email protected]
00011  *
00012  * =====================================================================================
00013  */
00014 
00015 #include "Announce.h"
00016 #include "DataArray3D.h"
00017 #include "FiniteElementTools.h"
00018 #include "TriangularQuadrature.h"
00019 #include "GaussQuadrature.h"
00020 #include "GaussLobattoQuadrature.h"
00021 #include "SparseMatrix.h"
00022 #include "STLStringHelper.h"
00023 #include "LinearRemapFV.h"
00024 
00025 #include "moab/Remapping/TempestOnlineMap.hpp"
00026 #include "DebugOutput.hpp"
00027 #include "moab/TupleList.hpp"
00028 
00029 #include <fstream>
00030 #include <cmath>
00031 #include <cstdlib>
00032 
00033 #ifdef MOAB_HAVE_NETCDFPAR
00034 #include "netcdfcpp_par.hpp"
00035 #else
00036 #include "netcdfcpp.h"
00037 #endif
00038 
00039 ///////////////////////////////////////////////////////////////////////////////
00040 
00041 // #define VERBOSE
00042 // #define VVERBOSE
00043 // #define CHECK_INCREASING_DOF
00044 
00045 ///////////////////////////////////////////////////////////////////////////////
00046 
00047 #define MPI_CHK_ERR( err )                                          \
00048     if( err )                                                       \
00049     {                                                               \
00050         std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
00051         std::cout << "\nMPI Aborting... \n";                        \
00052         return moab::MB_FAILURE;                                    \
00053     }
00054 
00055 moab::TempestOnlineMap::TempestOnlineMap( moab::TempestRemapper* remapper ) : OfflineMap(), m_remapper( remapper )
00056 {
00057     // Get the references for the MOAB core objects
00058     m_interface = m_remapper->get_interface();
00059 #ifdef MOAB_HAVE_MPI
00060     m_pcomm = m_remapper->get_parallel_communicator();
00061 #endif
00062 
00063     // Update the references to the meshes
00064     m_meshInput    = remapper->GetMesh( moab::Remapper::SourceMesh );
00065     m_meshInputCov = remapper->GetCoveringMesh();
00066     m_meshOutput   = remapper->GetMesh( moab::Remapper::TargetMesh );
00067     m_meshOverlap  = remapper->GetMesh( moab::Remapper::OverlapMesh );
00068 
00069     is_parallel = false;
00070     is_root     = true;
00071     rank        = 0;
00072     root_proc   = rank;
00073     size        = 1;
00074 #ifdef MOAB_HAVE_MPI
00075     int flagInit;
00076     MPI_Initialized( &flagInit );
00077     if( flagInit )
00078     {
00079         is_parallel = true;
00080         assert( m_pcomm != NULL );
00081         rank    = m_pcomm->rank();
00082         size    = m_pcomm->size();
00083         is_root = ( rank == 0 );
00084     }
00085 #endif
00086 
00087     // Compute and store the total number of source and target DoFs corresponding
00088     // to number of rows and columns in the mapping.
00089 
00090     // Initialize dimension information from file
00091     this->setup_sizes_dimensions();
00092 
00093     // Build a matrix of source and target discretization so that we know how to assign
00094     // the global DoFs in parallel for the mapping weights
00095     // For example, FV->FV: rows X cols = faces_source X faces_target
00096 }
00097 
00098 void moab::TempestOnlineMap::setup_sizes_dimensions()
00099 {
00100     if( m_meshInputCov )
00101     {
00102         std::vector< std::string > dimNames;
00103         std::vector< int > dimSizes;
00104         if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
00105         {
00106             dimNames.push_back( "lat" );
00107             dimNames.push_back( "lon" );
00108             dimSizes.resize( 2, 0 );
00109             dimSizes[0] = m_remapper->m_source_metadata[1];
00110             dimSizes[1] = m_remapper->m_source_metadata[2];
00111         }
00112         else
00113         {
00114             dimNames.push_back( "num_elem" );
00115             dimSizes.push_back( m_meshInputCov->faces.size() );
00116         }
00117 
00118         this->InitializeSourceDimensions( dimNames, dimSizes );
00119     }
00120 
00121     if( m_meshOutput )
00122     {
00123         std::vector< std::string > dimNames;
00124         std::vector< int > dimSizes;
00125         if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
00126         {
00127             dimNames.push_back( "lat" );
00128             dimNames.push_back( "lon" );
00129             dimSizes.resize( 2, 0 );
00130             dimSizes[0] = m_remapper->m_target_metadata[1];
00131             dimSizes[1] = m_remapper->m_target_metadata[2];
00132         }
00133         else
00134         {
00135             dimNames.push_back( "num_elem" );
00136             dimSizes.push_back( m_meshOutput->faces.size() );
00137         }
00138 
00139         this->InitializeTargetDimensions( dimNames, dimSizes );
00140     }
00141 }
00142 
00143 ///////////////////////////////////////////////////////////////////////////////
00144 
00145 moab::TempestOnlineMap::~TempestOnlineMap()
00146 {
00147     m_interface = NULL;
00148 #ifdef MOAB_HAVE_MPI
00149     m_pcomm = NULL;
00150 #endif
00151     m_meshInput   = NULL;
00152     m_meshOutput  = NULL;
00153     m_meshOverlap = NULL;
00154 }
00155 
00156 ///////////////////////////////////////////////////////////////////////////////
00157 
00158 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags( const std::string srcDofTagName,
00159                                                        const std::string tgtDofTagName )
00160 {
00161     moab::ErrorCode rval;
00162 
00163     int tagSize = 0;
00164     tagSize     = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
00165     rval =
00166         m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
00167 
00168     if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
00169     {
00170         int ntot_elements = 0, nelements = m_remapper->m_source_entities.size();
00171 #ifdef MOAB_HAVE_MPI
00172         int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
00173         if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
00174 #else
00175         ntot_elements = nelements;
00176 #endif
00177 
00178         rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_covering_source_entities,
00179                                                    &m_remapper->m_source_entities, srcDofTagName, m_nDofsPEl_Src );MB_CHK_ERR( rval );
00180 
00181         rval = m_interface->tag_get_handle( srcDofTagName.c_str(), m_nDofsPEl_Src * m_nDofsPEl_Src, MB_TYPE_INTEGER,
00182                                             this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR( rval );
00183     }
00184     else
00185         MB_CHK_ERR( rval );
00186 
00187     tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
00188     rval =
00189         m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
00190     if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
00191     {
00192         int ntot_elements = 0, nelements = m_remapper->m_target_entities.size();
00193 #ifdef MOAB_HAVE_MPI
00194         int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
00195         if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
00196 #else
00197         ntot_elements = nelements;
00198 #endif
00199 
00200         rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_target_entities, NULL, tgtDofTagName,
00201                                                    m_nDofsPEl_Dest );MB_CHK_ERR( rval );
00202 
00203         rval = m_interface->tag_get_handle( tgtDofTagName.c_str(), m_nDofsPEl_Dest * m_nDofsPEl_Dest, MB_TYPE_INTEGER,
00204                                             this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR( rval );
00205     }
00206     else
00207         MB_CHK_ERR( rval );
00208 
00209     return moab::MB_SUCCESS;
00210 }
00211 
00212 ///////////////////////////////////////////////////////////////////////////////
00213 
00214 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation( DiscretizationType srcType,
00215                                                               bool isSrcContinuous,
00216                                                               DataArray3D< int >* srcdataGLLNodes,
00217                                                               DataArray3D< int >* srcdataGLLNodesSrc,
00218                                                               DiscretizationType destType,
00219                                                               bool isTgtContinuous,
00220                                                               DataArray3D< int >* tgtdataGLLNodes )
00221 {
00222     moab::ErrorCode rval;
00223     std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
00224     std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
00225 
00226     // We are assuming that these are element based tags that are sized: np * np
00227     m_srcDiscType  = srcType;
00228     m_destDiscType = destType;
00229 
00230     bool vprint = is_root && false;
00231 
00232 #ifdef VVERBOSE
00233     {
00234         src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
00235         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
00236         locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
00237         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
00238         tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
00239         rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
00240 
00241         if( is_root )
00242         {
00243             {
00244                 std::ofstream output_file( "sourcecov-gids-0.txt" );
00245                 output_file << "I, GDOF\n";
00246                 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
00247                     output_file << i << ", " << src_soln_gdofs[i] << "\n";
00248 
00249                 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
00250                 m_nTotDofs_SrcCov = 0;
00251                 if( isSrcContinuous )
00252                     dgll_cgll_covcol_ldofmap.resize(
00253                         m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
00254                 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
00255                 {
00256                     for( int p = 0; p < m_nDofsPEl_Src; p++ )
00257                     {
00258                         for( int q = 0; q < m_nDofsPEl_Src; q++ )
00259                         {
00260                             const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
00261                             const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
00262                             if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
00263                             {
00264                                 m_nTotDofs_SrcCov++;
00265                                 dgll_cgll_covcol_ldofmap[localDOF] = true;
00266                             }
00267                             output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
00268                                         << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
00269                         }
00270                     }
00271                 }
00272                 output_file.flush();  // required here
00273                 output_file.close();
00274                 dgll_cgll_covcol_ldofmap.clear();
00275             }
00276 
00277             {
00278                 std::ofstream output_file( "source-gids-0.txt" );
00279                 output_file << "I, GDOF\n";
00280                 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
00281                     output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
00282 
00283                 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
00284                 m_nTotDofs_Src = 0;
00285                 if( isSrcContinuous )
00286                     dgll_cgll_col_ldofmap.resize(
00287                         m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
00288                 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
00289                 {
00290                     for( int p = 0; p < m_nDofsPEl_Src; p++ )
00291                     {
00292                         for( int q = 0; q < m_nDofsPEl_Src; q++ )
00293                         {
00294                             const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
00295                             const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
00296                             if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
00297                             {
00298                                 m_nTotDofs_Src++;
00299                                 dgll_cgll_col_ldofmap[localDOF] = true;
00300                             }
00301                             output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
00302                                         << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
00303                         }
00304                     }
00305                 }
00306                 output_file.flush();  // required here
00307                 output_file.close();
00308                 dgll_cgll_col_ldofmap.clear();
00309             }
00310 
00311             {
00312                 std::ofstream output_file( "target-gids-0.txt" );
00313                 output_file << "I, GDOF\n";
00314                 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
00315                     output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
00316 
00317                 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
00318                 m_nTotDofs_Dest = 0;
00319 
00320                 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
00321                 {
00322                     output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
00323                                 << m_nTotDofs_Dest << "\n";
00324                     m_nTotDofs_Dest++;
00325                 }
00326 
00327                 output_file.flush();  // required here
00328                 output_file.close();
00329             }
00330         }
00331         else
00332         {
00333             {
00334                 std::ofstream output_file( "sourcecov-gids-1.txt" );
00335                 output_file << "I, GDOF\n";
00336                 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
00337                     output_file << i << ", " << src_soln_gdofs[i] << "\n";
00338 
00339                 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
00340                 m_nTotDofs_SrcCov = 0;
00341                 if( isSrcContinuous )
00342                     dgll_cgll_covcol_ldofmap.resize(
00343                         m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
00344                 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
00345                 {
00346                     for( int p = 0; p < m_nDofsPEl_Src; p++ )
00347                     {
00348                         for( int q = 0; q < m_nDofsPEl_Src; q++ )
00349                         {
00350                             const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
00351                             const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
00352                             if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
00353                             {
00354                                 m_nTotDofs_SrcCov++;
00355                                 dgll_cgll_covcol_ldofmap[localDOF] = true;
00356                             }
00357                             output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
00358                                         << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
00359                         }
00360                     }
00361                 }
00362                 output_file.flush();  // required here
00363                 output_file.close();
00364                 dgll_cgll_covcol_ldofmap.clear();
00365             }
00366 
00367             {
00368                 std::ofstream output_file( "source-gids-1.txt" );
00369                 output_file << "I, GDOF\n";
00370                 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
00371                     output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
00372 
00373                 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
00374                 m_nTotDofs_Src = 0;
00375                 if( isSrcContinuous )
00376                     dgll_cgll_col_ldofmap.resize(
00377                         m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
00378                 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
00379                 {
00380                     for( int p = 0; p < m_nDofsPEl_Src; p++ )
00381                     {
00382                         for( int q = 0; q < m_nDofsPEl_Src; q++ )
00383                         {
00384                             const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
00385                             const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
00386                             if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
00387                             {
00388                                 m_nTotDofs_Src++;
00389                                 dgll_cgll_col_ldofmap[localDOF] = true;
00390                             }
00391                             output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
00392                                         << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
00393                         }
00394                     }
00395                 }
00396                 output_file.flush();  // required here
00397                 output_file.close();
00398                 dgll_cgll_col_ldofmap.clear();
00399             }
00400 
00401             {
00402                 std::ofstream output_file( "target-gids-1.txt" );
00403                 output_file << "I, GDOF\n";
00404                 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
00405                     output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
00406 
00407                 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
00408                 m_nTotDofs_Dest = 0;
00409 
00410                 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
00411                 {
00412                     output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
00413                                 << m_nTotDofs_Dest << "\n";
00414                     m_nTotDofs_Dest++;
00415                 }
00416 
00417                 output_file.flush();  // required here
00418                 output_file.close();
00419             }
00420         }
00421     }
00422 #endif
00423 
00424     // Now compute the mapping and store it for the covering mesh
00425     int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
00426     if( m_remapper->point_cloud_source )
00427     {
00428         assert( m_nDofsPEl_Src == 1 );
00429         col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
00430         col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
00431         src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
00432         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
00433         srcTagSize = 1;
00434     }
00435     else
00436     {
00437         col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
00438         col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
00439         src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
00440         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
00441     }
00442 
00443     // std::cout << "TOnlineMap: Process: " << rank << " and covering entities = [" <<
00444     // col_dofmap.size() << ", " << src_soln_gdofs.size() << "]\n"; MPI_Barrier(MPI_COMM_WORLD);
00445 
00446 #ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
00447     unsigned maxSrcIndx = 0;
00448 
00449     // for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
00450     std::vector< int > locdofs( srcTagSize );
00451     std::map< Node, moab::EntityHandle > mapLocalMBNodes;
00452     double elcoords[3];
00453     for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
00454     {
00455         EntityHandle eh = m_remapper->m_covering_source_entities[iel];
00456         rval            = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
00457         Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
00458         mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
00459     }
00460 
00461     const NodeVector& nodes = m_remapper->m_covering_source->nodes;
00462     for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
00463     {
00464         const Face& face = m_remapper->m_covering_source->faces[j];
00465 
00466         Node centroid;
00467         centroid.x = centroid.y = centroid.z = 0.0;
00468         for( unsigned l = 0; l < face.edges.size(); ++l )
00469         {
00470             centroid.x += nodes[face[l]].x;
00471             centroid.y += nodes[face[l]].y;
00472             centroid.z += nodes[face[l]].z;
00473         }
00474         const double factor = 1.0 / face.edges.size();
00475         centroid.x *= factor;
00476         centroid.y *= factor;
00477         centroid.z *= factor;
00478 
00479         EntityHandle current_eh;
00480         if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
00481         {
00482             current_eh = mapLocalMBNodes[centroid];
00483         }
00484 
00485         rval = m_interface->tag_get_data( m_dofTagSrc, &current_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
00486         for( int p = 0; p < m_nDofsPEl_Src; p++ )
00487         {
00488             for( int q = 0; q < m_nDofsPEl_Src; q++ )
00489             {
00490                 const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
00491                 const int offsetDOF = p * m_nDofsPEl_Src + q;
00492                 maxSrcIndx          = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
00493                 std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
00494                           << ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
00495             }
00496         }
00497     }
00498 #endif
00499 
00500     m_nTotDofs_SrcCov = 0;
00501     if( srcdataGLLNodes == NULL )
00502     { /* we only have a mapping for elements as DoFs */
00503         for( unsigned i = 0; i < col_gdofmap.size(); ++i )
00504         {
00505             assert( src_soln_gdofs[i] > 0 );
00506             col_gdofmap[i]     = src_soln_gdofs[i] - 1;
00507             col_dtoc_dofmap[i] = i;
00508             if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
00509             m_nTotDofs_SrcCov++;
00510         }
00511     }
00512     else
00513     {
00514         if( isSrcContinuous )
00515             dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
00516         // Put these remap coefficients into the SparseMatrix map
00517         for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
00518         {
00519             for( int p = 0; p < m_nDofsPEl_Src; p++ )
00520             {
00521                 for( int q = 0; q < m_nDofsPEl_Src; q++ )
00522                 {
00523                     const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
00524                     const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
00525                     if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
00526                     {
00527                         m_nTotDofs_SrcCov++;
00528                         dgll_cgll_covcol_ldofmap[localDOF] = true;
00529                     }
00530                     if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
00531                     assert( src_soln_gdofs[offsetDOF] > 0 );
00532                     col_gdofmap[localDOF]      = src_soln_gdofs[offsetDOF] - 1;
00533                     col_dtoc_dofmap[offsetDOF] = localDOF;
00534                     if( vprint )
00535                         std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", "
00536                                   << localDOF << ", " << col_gdofmap[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
00537                 }
00538             }
00539         }
00540     }
00541 
00542     if( m_remapper->point_cloud_source )
00543     {
00544         assert( m_nDofsPEl_Src == 1 );
00545         srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
00546         srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
00547         locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
00548         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
00549     }
00550     else
00551     {
00552         srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
00553         srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
00554         locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
00555         rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
00556     }
00557 
00558     // Now compute the mapping and store it for the original source mesh
00559     m_nTotDofs_Src = 0;
00560     if( srcdataGLLNodesSrc == NULL )
00561     { /* we only have a mapping for elements as DoFs */
00562         for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
00563         {
00564             assert( locsrc_soln_gdofs[i] > 0 );
00565             srccol_gdofmap[i]     = locsrc_soln_gdofs[i] - 1;
00566             srccol_dtoc_dofmap[i] = i;
00567             m_nTotDofs_Src++;
00568         }
00569     }
00570     else
00571     {
00572         if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
00573         // Put these remap coefficients into the SparseMatrix map
00574         for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
00575         {
00576             for( int p = 0; p < m_nDofsPEl_Src; p++ )
00577             {
00578                 for( int q = 0; q < m_nDofsPEl_Src; q++ )
00579                 {
00580                     const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
00581                     const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
00582                     if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
00583                     {
00584                         m_nTotDofs_Src++;
00585                         dgll_cgll_col_ldofmap[localDOF] = true;
00586                     }
00587                     if( !isSrcContinuous ) m_nTotDofs_Src++;
00588                     assert( locsrc_soln_gdofs[offsetDOF] > 0 );
00589                     srccol_gdofmap[localDOF]      = locsrc_soln_gdofs[offsetDOF] - 1;
00590                     srccol_dtoc_dofmap[offsetDOF] = localDOF;
00591                 }
00592             }
00593         }
00594     }
00595 
00596     int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
00597     if( m_remapper->point_cloud_target )
00598     {
00599         assert( m_nDofsPEl_Dest == 1 );
00600         row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
00601         row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
00602         tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
00603         rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
00604         tgtTagSize = 1;
00605     }
00606     else
00607     {
00608         row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
00609         row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
00610         tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
00611         rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
00612     }
00613 
00614     // Now compute the mapping and store it for the target mesh
00615     // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
00616     m_nTotDofs_Dest = 0;
00617     if( tgtdataGLLNodes == NULL )
00618     { /* we only have a mapping for elements as DoFs */
00619         for( unsigned i = 0; i < row_gdofmap.size(); ++i )
00620         {
00621             assert( tgt_soln_gdofs[i] > 0 );
00622             row_gdofmap[i]     = tgt_soln_gdofs[i] - 1;
00623             row_dtoc_dofmap[i] = i;
00624             if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
00625             m_nTotDofs_Dest++;
00626         }
00627     }
00628     else
00629     {
00630         if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
00631         // Put these remap coefficients into the SparseMatrix map
00632         for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
00633         {
00634             for( int p = 0; p < m_nDofsPEl_Dest; p++ )
00635             {
00636                 for( int q = 0; q < m_nDofsPEl_Dest; q++ )
00637                 {
00638                     const int localDOF  = ( *tgtdataGLLNodes )[p][q][j] - 1;
00639                     const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
00640                     if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
00641                     {
00642                         m_nTotDofs_Dest++;
00643                         dgll_cgll_row_ldofmap[localDOF] = true;
00644                     }
00645                     if( !isTgtContinuous ) m_nTotDofs_Dest++;
00646                     assert( tgt_soln_gdofs[offsetDOF] > 0 );
00647                     row_gdofmap[localDOF]      = tgt_soln_gdofs[offsetDOF] - 1;
00648                     row_dtoc_dofmap[offsetDOF] = localDOF;
00649                     if( vprint )
00650                         std::cout << "Row: " << m_remapper->lid_to_gid_tgt[j] << ", " << offsetDOF << ", " << localDOF
00651                                   << ", " << row_gdofmap[offsetDOF] << ", " << m_nTotDofs_Dest << "\n";
00652                 }
00653             }
00654         }
00655     }
00656 
00657     // Let us also allocate the local representation of the sparse matrix
00658 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
00659     if( vprint )
00660     {
00661         std::cout << "[" << rank << "]"
00662                   << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size() << ", col = " << m_nTotDofs_Src
00663                   << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
00664         // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
00665     }
00666 #endif
00667 
00668     // check monotonicity of row_gdofmap and col_gdofmap
00669 #ifdef CHECK_INCREASING_DOF
00670     for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
00671     {
00672         if( row_gdofmap[i] > row_gdofmap[i + 1] )
00673             std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
00674                       << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
00675     }
00676     for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
00677     {
00678         if( col_gdofmap[i] > col_gdofmap[i + 1] )
00679             std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
00680                       << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
00681     }
00682 #endif
00683 
00684     return moab::MB_SUCCESS;
00685 }
00686 
00687 moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs( std::vector< int >& values_entities )
00688 {
00689     // col_gdofmap has global dofs , that should be in the list of values, such that
00690     // row_dtoc_dofmap[offsetDOF] = localDOF;
00691     // //  we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
00692     // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
00693     // form first inverse
00694 
00695     col_dtoc_dofmap.resize( values_entities.size() );
00696     for( int j = 0; j < (int)values_entities.size(); j++ )
00697     {
00698         if( colMap.find( values_entities[j] - 1 ) != colMap.end() )
00699             col_dtoc_dofmap[j] = colMap[values_entities[j] - 1];
00700         else
00701         {
00702             col_dtoc_dofmap[j] = -1;  // signal that this value should not be used in
00703             // std::cout <<"values_entities[j] -  1: " << values_entities[j] -  1 <<" at index j = " << j <<  " not
00704             // found in colMap \n";
00705         }
00706     }
00707     return moab::MB_SUCCESS;
00708 }
00709 
00710 moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs( std::vector< int >& values_entities )
00711 {
00712     // row_dtoc_dofmap = values_entities; // needs to point to local
00713     //  we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
00714 
00715     row_dtoc_dofmap.resize( values_entities.size() );
00716     for( int j = 0; j < (int)values_entities.size(); j++ )
00717     {
00718         if( rowMap.find( values_entities[j] - 1 ) != rowMap.end() )
00719             row_dtoc_dofmap[j] = rowMap[values_entities[j] - 1];  // values are 1 based, but rowMap, colMap are not
00720         else
00721         {
00722             row_dtoc_dofmap[j] = -1;  // not all values are used
00723             // std::cout <<"values_entities[j] -  1: " << values_entities[j] -  1 <<" at index j = " << j <<  " not
00724             // found in rowMap \n";
00725         }
00726     }
00727     return moab::MB_SUCCESS;
00728 }
00729 ///////////////////////////////////////////////////////////////////////////////
00730 
00731 moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights( std::string strInputType,
00732                                                                   std::string strOutputType,
00733                                                                   const GenerateOfflineMapAlgorithmOptions& mapOptions,
00734                                                                   const std::string& srcDofTagName,
00735                                                                   const std::string& tgtDofTagName )
00736 {
00737     NcError error( NcError::silent_nonfatal );
00738 
00739     moab::DebugOutput dbgprint( std::cout, rank, 0 );
00740     dbgprint.set_prefix( "[TempestOnlineMap]: " );
00741     moab::ErrorCode rval;
00742 
00743     const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
00744     const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
00745     const bool m_bPointCloud       = m_bPointCloudSource || m_bPointCloudTarget;
00746 
00747     try
00748     {
00749         // Check command line parameters (data type arguments)
00750         STLStringHelper::ToLower( strInputType );
00751         STLStringHelper::ToLower( strOutputType );
00752 
00753         DiscretizationType eInputType;
00754         DiscretizationType eOutputType;
00755 
00756         if( strInputType == "fv" )
00757         {
00758             eInputType = DiscretizationType_FV;
00759         }
00760         else if( strInputType == "cgll" )
00761         {
00762             eInputType = DiscretizationType_CGLL;
00763         }
00764         else if( strInputType == "dgll" )
00765         {
00766             eInputType = DiscretizationType_DGLL;
00767         }
00768         else if( strInputType == "pcloud" )
00769         {
00770             eInputType = DiscretizationType_PCLOUD;
00771         }
00772         else
00773         {
00774             _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
00775         }
00776 
00777         if( strOutputType == "fv" )
00778         {
00779             eOutputType = DiscretizationType_FV;
00780         }
00781         else if( strOutputType == "cgll" )
00782         {
00783             eOutputType = DiscretizationType_CGLL;
00784         }
00785         else if( strOutputType == "dgll" )
00786         {
00787             eOutputType = DiscretizationType_DGLL;
00788         }
00789         else if( strOutputType == "pcloud" )
00790         {
00791             eOutputType = DiscretizationType_PCLOUD;
00792         }
00793         else
00794         {
00795             _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
00796         }
00797 
00798         // set all required input params
00799         m_bConserved  = !mapOptions.fNoConservation;
00800         m_eInputType  = eInputType;
00801         m_eOutputType = eOutputType;
00802 
00803         // Method flags
00804         std::string strMapAlgorithm( "" );
00805         int nMonotoneType    = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
00806 
00807         // Make an index of method arguments
00808         std::set< std::string > setMethodStrings;
00809         {
00810             int iLast = 0;
00811             for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
00812             {
00813                 if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
00814                 {
00815                     std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
00816                     STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
00817                     if( strMethodString.length() > 0 )
00818                     {
00819                         setMethodStrings.insert( strMethodString );
00820                     }
00821                     iLast = i + 1;
00822                 }
00823             }
00824         }
00825 
00826         for( auto it : setMethodStrings )
00827         {
00828             // Piecewise constant monotonicity
00829             if( it == "mono2" )
00830             {
00831                 if( nMonotoneType != 0 )
00832                 {
00833                     _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
00834                 }
00835                 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
00836                 {
00837                     _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
00838                 }
00839                 nMonotoneType = 2;
00840 
00841                 // Piecewise linear monotonicity
00842             }
00843             else if( it == "mono3" )
00844             {
00845                 if( nMonotoneType != 0 )
00846                 {
00847                     _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
00848                 }
00849                 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
00850                 {
00851                     _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
00852                 }
00853                 nMonotoneType = 3;
00854 
00855                 // Volumetric remapping from FV to GLL
00856             }
00857             else if( it == "volumetric" )
00858             {
00859                 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
00860                 {
00861                     _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
00862                 }
00863                 strMapAlgorithm = "volumetric";
00864 
00865                 // Inverse distance mapping
00866             }
00867             else if( it == "invdist" )
00868             {
00869                 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
00870                 {
00871                     _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
00872                 }
00873                 strMapAlgorithm = "invdist";
00874 
00875                 // Delaunay triangulation mapping
00876             }
00877             else if( it == "delaunay" )
00878             {
00879                 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
00880                 {
00881                     _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
00882                 }
00883                 strMapAlgorithm = "delaunay";
00884 
00885                 // Bilinear
00886             }
00887             else if( it == "bilin" )
00888             {
00889                 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
00890                 {
00891                     _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
00892                 }
00893                 strMapAlgorithm = "fvbilin";
00894 
00895                 // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
00896             }
00897             else if( it == "intbilin" )
00898             {
00899                 if( m_eOutputType != DiscretizationType_FV )
00900                 {
00901                     _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
00902                 }
00903                 if( m_eInputType == DiscretizationType_FV )
00904                 {
00905                     strMapAlgorithm = "fvintbilin";
00906                 }
00907                 else
00908                 {
00909                     strMapAlgorithm = "mono3";
00910                 }
00911 
00912                 // Integrated bilinear with generalized Barycentric coordinates
00913             }
00914             else if( it == "intbilingb" )
00915             {
00916                 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
00917                 {
00918                     _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
00919                 }
00920                 strMapAlgorithm = "fvintbilingb";
00921             }
00922             else
00923             {
00924                 _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
00925             }
00926         }
00927 
00928         m_nDofsPEl_Src =
00929             ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
00930                                                                                                  : mapOptions.nPin );
00931         m_nDofsPEl_Dest =
00932             ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
00933                                                                                                    : mapOptions.nPout );
00934 
00935         rval = SetDOFmapTags( srcDofTagName, tgtDofTagName );MB_CHK_ERR( rval );
00936 
00937         ///   the tag should be created already in the e3sm workflow; if not, create it here
00938         Tag areaTag;
00939         rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
00940                                              MB_TAG_DENSE |  MB_TAG_EXCL | MB_TAG_CREAT );
00941         if ( MB_ALREADY_ALLOCATED == rval )
00942         {
00943             if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
00944         }
00945 
00946         double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0;
00947         if( !m_bPointCloudSource )
00948         {
00949             // Calculate Input Mesh Face areas
00950             if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
00951             double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
00952             rval = m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea );MB_CHK_ERR( rval ); 
00953             dTotalAreaInput            = dTotalAreaInput_loc;
00954 #ifdef MOAB_HAVE_MPI
00955             if( m_pcomm )
00956                 MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
00957 #endif
00958             if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
00959 
00960             // Input mesh areas
00961             m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
00962             // we do not need to set the area on coverage mesh, only on source and target meshes
00963             // rval = m_interface->tag_set_data( areaTag, m_remapper->m_covering_source_entities, m_meshInputCov->vecFaceArea );MB_CHK_ERR( rval );
00964         }
00965 
00966         if( !m_bPointCloudTarget )
00967         {
00968             // Calculate Output Mesh Face areas
00969             if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
00970             double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
00971             dTotalAreaOutput            = dTotalAreaOutput_loc;
00972 #ifdef MOAB_HAVE_MPI
00973             if( m_pcomm )
00974                 MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
00975 #endif
00976             if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
00977             rval = m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea );MB_CHK_ERR( rval ); 
00978         }
00979 
00980         if( !m_bPointCloud )
00981         {
00982             // Verify that overlap mesh is in the correct order, only if size == 1
00983             if( 1 == size )
00984             {
00985                 int ixSourceFaceMax = ( -1 );
00986                 int ixTargetFaceMax = ( -1 );
00987 
00988                 if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() )
00989                 {
00990                     _EXCEPTIONT( "Invalid overlap mesh:\n"
00991                                  "    Possible mesh file corruption?" );
00992                 }
00993 
00994                 for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ )
00995                 {
00996                     if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
00997                         ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;
00998 
00999                     if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
01000                         ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
01001                 }
01002 
01003                 // Check for forward correspondence in overlap mesh
01004                 if( m_meshInputCov->faces.size() - ixSourceFaceMax == 0 )
01005                 {
01006                     if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" );
01007                 }
01008                 else if( m_meshOutput->faces.size() - ixSourceFaceMax == 0 )
01009                 {  // Check for reverse correspondence in overlap mesh
01010                     if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" );
01011 
01012                     // Reorder overlap mesh
01013                     m_meshOverlap->ExchangeFirstAndSecondMesh();
01014                 }
01015                 else
01016                 {  // No correspondence found
01017                     _EXCEPTION4( "Invalid overlap mesh:\n    No correspondence found with input and output meshes "
01018                                  "(%i,%i) vs (%i,%i)",
01019                                  m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax,
01020                                  ixTargetFaceMax );
01021                 }
01022             }
01023 
01024             // Calculate Face areas
01025             if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
01026             double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false );
01027             double dTotalAreaOverlap     = dTotalAreaOverlap_loc;
01028 #ifdef MOAB_HAVE_MPI
01029             if( m_pcomm )
01030                 MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
01031 #endif
01032             if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
01033 
01034             // Correct areas to match the areas calculated in the overlap mesh
01035             // if (fCorrectAreas) // In MOAB-TempestRemap, we will always keep this to be true
01036             {
01037                 if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
01038                 DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
01039                 DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
01040 
01041                 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
01042                 assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
01043                 assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
01044 
01045                 assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
01046                 assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
01047 
01048                 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
01049                 {
01050                     if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
01051                         continue;  // skip this cell since it is ghosted
01052 
01053                     // let us recompute the source/target areas based on overlap mesh areas
01054                     assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
01055                     dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
01056                     assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
01057                     dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
01058                 }
01059 
01060                 for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
01061                 {
01062                     if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
01063                     {
01064                         m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
01065                     }
01066                 }
01067                 for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
01068                 {
01069                     if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
01070                     {
01071                         m_meshOutput->vecFaceArea[i] = dTargetArea[i];
01072                     }
01073                 }
01074             }
01075 
01076             // Set source mesh areas in map
01077             if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
01078             {
01079                 this->SetSourceAreas( m_meshInputCov->vecFaceArea );
01080                 if( m_meshInputCov->vecMask.size() )
01081                 {
01082                     this->SetSourceMask( m_meshInputCov->vecMask );
01083                 }
01084             }
01085 
01086             // Set target mesh areas in map
01087             if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
01088             {
01089                 this->SetTargetAreas( m_meshOutput->vecFaceArea );
01090                 if( m_meshOutput->vecMask.size() )
01091                 {
01092                     this->SetTargetMask( m_meshOutput->vecMask );
01093                 }
01094             }
01095 
01096             /*
01097                 // Recalculate input mesh area from overlap mesh
01098                 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
01099                     dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
01100                     dbgprint.printf(0, "Recalculating source mesh areas\n");
01101                     dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
01102                     dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
01103                 }
01104             */
01105         }
01106 
01107         // Finite volume input / Finite volume output
01108         if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
01109         {
01110             // Generate reverse node array and edge map
01111             if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
01112             if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
01113 
01114             // Initialize coordinates for map
01115             this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
01116             this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
01117 
01118             // Finite volume input / Finite element output
01119             rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, false, NULL );MB_CHK_ERR( rval );
01120 
01121             // Construct remap for FV-FV
01122             if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
01123 
01124             // Construct OfflineMap
01125             if( strMapAlgorithm == "invdist" )
01126             {
01127                 if( is_root ) AnnounceStartBlock( "Calculating map (invdist)" );
01128                 LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
01129             }
01130             else if( strMapAlgorithm == "delaunay" )
01131             {
01132                 if( is_root ) AnnounceStartBlock( "Calculating map (delaunay)" );
01133                 LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
01134             }
01135             else if( strMapAlgorithm == "fvintbilin" )
01136             {
01137                 if( is_root ) AnnounceStartBlock( "Calculating map (intbilin)" );
01138                 LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
01139             }
01140             else if( strMapAlgorithm == "fvintbilingb" )
01141             {
01142                 if( is_root ) AnnounceStartBlock( "Calculating map (intbilingb)" );
01143                 LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
01144             }
01145             else if( strMapAlgorithm == "fvbilin" )
01146             {
01147                 if( is_root ) AnnounceStartBlock( "Calculating map (bilin)" );
01148                 LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
01149             }
01150             else
01151             {
01152                 if( is_root ) AnnounceStartBlock( "Calculating map (default)" );
01153                 // LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
01154                 //                   ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
01155                 LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
01156             }
01157         }
01158         else if( eInputType == DiscretizationType_FV )
01159         {
01160             DataArray3D< double > dataGLLJacobian;
01161 
01162             if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
01163             double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
01164                                                           dataGLLNodesDest, dataGLLJacobian );
01165 
01166             double dNumericalArea = dNumericalArea_loc;
01167 #ifdef MOAB_HAVE_MPI
01168             if( m_pcomm )
01169                 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
01170 #endif
01171             if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
01172 
01173             // Initialize coordinates for map
01174             this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
01175             this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
01176 
01177             // Generate the continuous Jacobian
01178             bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
01179 
01180             if( eOutputType == DiscretizationType_CGLL )
01181             {
01182                 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
01183             }
01184             else
01185             {
01186                 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
01187             }
01188 
01189             // Generate reverse node array and edge map
01190             if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
01191             if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
01192 
01193             // Finite volume input / Finite element output
01194             rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType,
01195                                                ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
01196 
01197             // Generate remap weights
01198             if( strMapAlgorithm == "volumetric" )
01199             {
01200                 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
01201                 LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
01202                                                dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
01203                                                nMonotoneType, fContinuous, mapOptions.fNoConservation );
01204             }
01205             else
01206             {
01207                 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
01208                 LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
01209                                     this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
01210                                     mapOptions.fNoConservation );
01211             }
01212         }
01213         else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
01214         {
01215             DataArray3D< double > dataGLLJacobian;
01216             if( !m_bPointCloudSource )
01217             {
01218                 // Generate reverse node array and edge map
01219                 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
01220                 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
01221 
01222                 // Initialize coordinates for map
01223                 if( eInputType == DiscretizationType_FV )
01224                 {
01225                     this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
01226                 }
01227                 else
01228                 {
01229                     if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
01230                     DataArray3D< double > dataGLLJacobianSrc;
01231                     GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
01232                                       dataGLLJacobian );
01233                     GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
01234                                       dataGLLJacobianSrc );
01235                 }
01236             }
01237             // else { /* Source is a point cloud dataset */ }
01238 
01239             if( !m_bPointCloudTarget )
01240             {
01241                 // Generate reverse node array and edge map
01242                 if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
01243                 if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
01244 
01245                 // Initialize coordinates for map
01246                 if( eOutputType == DiscretizationType_FV )
01247                 {
01248                     this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
01249                 }
01250                 else
01251                 {
01252                     if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
01253                     GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
01254                                       dataGLLJacobian );
01255                 }
01256             }
01257             // else { /* Target is a point cloud dataset */ }
01258 
01259             // Finite volume input / Finite element output
01260             rval = this->SetDOFmapAssociation(
01261                 eInputType, ( eInputType == DiscretizationType_CGLL ),
01262                 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrcCov ),
01263                 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrc ), eOutputType,
01264                 ( eOutputType == DiscretizationType_CGLL ), ( m_bPointCloudTarget ? NULL : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
01265 
01266             // Construct remap
01267             if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
01268             rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
01269         }
01270         else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
01271         {
01272             DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
01273 
01274             if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
01275             // double dNumericalAreaCov_loc =
01276             GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
01277                               dataGLLJacobian );
01278 
01279             double dNumericalArea_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
01280                                                           dataGLLNodesSrc, dataGLLJacobianSrc );
01281 
01282             // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area:
01283             // %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
01284             // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);
01285 
01286             double dNumericalArea = dNumericalArea_loc;
01287 #ifdef MOAB_HAVE_MPI
01288             if( m_pcomm )
01289                 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
01290 #endif
01291             if( is_root )
01292             {
01293                 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );
01294 
01295                 if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
01296                 {
01297                     dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
01298                                         "numerical area and geometric area\n" );
01299                 }
01300             }
01301 
01302             if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
01303             {
01304                 _EXCEPTIONT( "Number of element does not match between metadata and "
01305                              "input mesh" );
01306             }
01307 
01308             // Initialize coordinates for map
01309             this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
01310             this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
01311 
01312             // Generate the continuous Jacobian for input mesh
01313             bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
01314 
01315             if( eInputType == DiscretizationType_CGLL )
01316             {
01317                 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
01318             }
01319             else
01320             {
01321                 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
01322             }
01323 
01324             // Finite element input / Finite volume output
01325             rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
01326                                                &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, false, NULL );MB_CHK_ERR( rval );
01327 
01328             // Generate remap
01329             if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
01330 
01331             if( strMapAlgorithm == "volumetric" )
01332             {
01333                 _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
01334                              "GLL input mesh" );
01335             }
01336 
01337             LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
01338                                          mapOptions.fNoConservation );
01339         }
01340         else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
01341         {
01342             DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
01343             DataArray3D< double > dataGLLJacobianOut;
01344 
01345             // Input metadata
01346             if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
01347             GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
01348                               dataGLLJacobianIn );
01349 
01350             double dNumericalAreaSrc_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
01351                                                              dataGLLNodesSrc, dataGLLJacobianSrc );
01352             double dNumericalAreaIn      = dNumericalAreaSrc_loc;
01353 #ifdef MOAB_HAVE_MPI
01354             if( m_pcomm )
01355                 MPI_Reduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
01356 #endif
01357             if( is_root )
01358             {
01359                 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn );
01360 
01361                 if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
01362                 {
01363                     dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
01364                                         "numerical area and geometric area\n" );
01365                 }
01366             }
01367 
01368             // Output metadata
01369             if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
01370             double dNumericalAreaOut_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
01371                                                              dataGLLNodesDest, dataGLLJacobianOut );
01372             double dNumericalAreaOut     = dNumericalAreaOut_loc;
01373 #ifdef MOAB_HAVE_MPI
01374             if( m_pcomm )
01375                 MPI_Reduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
01376 #endif
01377             if( is_root )
01378             {
01379                 dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut );
01380 
01381                 if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
01382                 {
01383                     if( is_root )
01384                         dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh "
01385                                             "numerical area and geometric area\n" );
01386                 }
01387             }
01388 
01389             // Initialize coordinates for map
01390             this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
01391             this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
01392 
01393             // Generate the continuous Jacobian for input mesh
01394             bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
01395 
01396             if( eInputType == DiscretizationType_CGLL )
01397             {
01398                 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
01399             }
01400             else
01401             {
01402                 GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
01403             }
01404 
01405             // Generate the continuous Jacobian for output mesh
01406             bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
01407 
01408             if( eOutputType == DiscretizationType_CGLL )
01409             {
01410                 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
01411             }
01412             else
01413             {
01414                 GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
01415             }
01416 
01417             // Input Finite Element to Output Finite Element
01418             rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
01419                                                &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType,
01420                                                ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
01421 
01422             // Generate remap
01423             if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
01424 
01425             LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
01426                                        this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
01427                                        fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
01428         }
01429         else
01430         {
01431             _EXCEPTIONT( "Not implemented" );
01432         }
01433 
01434 #ifdef MOAB_HAVE_EIGEN3
01435         copy_tempest_sparsemat_to_eigen3();
01436 #endif
01437 
01438 #ifdef MOAB_HAVE_MPI
01439         {
01440             // Remove ghosted entities from overlap set
01441             moab::Range ghostedEnts;
01442             rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
01443             moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
01444             rval                                = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
01445         }
01446 #endif
01447         // Verify consistency, conservation and monotonicity, globally
01448         if( !mapOptions.fNoCheck )
01449         {
01450             if( is_root ) dbgprint.printf( 0, "Verifying map" );
01451             this->IsConsistent( 1.0e-8 );
01452             if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
01453 
01454             if( nMonotoneType != 0 )
01455             {
01456                 this->IsMonotone( 1.0e-12 );
01457             }
01458         }
01459     }
01460     catch( Exception& e )
01461     {
01462         dbgprint.printf( 0, "%s", e.ToString().c_str() );
01463         return ( moab::MB_FAILURE );
01464     }
01465     catch( ... )
01466     {
01467         return ( moab::MB_FAILURE );
01468     }
01469     return moab::MB_SUCCESS;
01470 }
01471 
01472 ///////////////////////////////////////////////////////////////////////////////
01473 
01474 int moab::TempestOnlineMap::IsConsistent( double dTolerance )
01475 {
01476 #ifndef MOAB_HAVE_MPI
01477 
01478     return OfflineMap::IsConsistent( dTolerance );
01479 
01480 #else
01481 
01482     // Get map entries
01483     DataArray1D< int > dataRows;
01484     DataArray1D< int > dataCols;
01485     DataArray1D< double > dataEntries;
01486 
01487     // Calculate row sums
01488     DataArray1D< double > dRowSums;
01489     m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
01490     dRowSums.Allocate( m_mapRemap.GetRows() );
01491 
01492     for( unsigned i = 0; i < dataRows.GetRows(); i++ )
01493     {
01494         dRowSums[dataRows[i]] += dataEntries[i];
01495     }
01496 
01497     // Verify all row sums are equal to 1
01498     int fConsistent = 0;
01499     for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
01500     {
01501         if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
01502         {
01503             fConsistent++;
01504             int rowGID = row_gdofmap[i];
01505             Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
01506         }
01507     }
01508 
01509     int ierr;
01510     int fConsistentGlobal = 0;
01511     ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
01512     if( ierr != MPI_SUCCESS ) return -1;
01513 
01514     return fConsistentGlobal;
01515 #endif
01516 }
01517 
01518 ///////////////////////////////////////////////////////////////////////////////
01519 
01520 int moab::TempestOnlineMap::IsConservative( double dTolerance )
01521 {
01522 #ifndef MOAB_HAVE_MPI
01523 
01524     return OfflineMap::IsConservative( dTolerance );
01525 
01526 #else
01527     // return OfflineMap::IsConservative(dTolerance);
01528 
01529     int ierr;
01530     // Get map entries
01531     DataArray1D< int > dataRows;
01532     DataArray1D< int > dataCols;
01533     DataArray1D< double > dataEntries;
01534     const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
01535     const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
01536 
01537     // Calculate column sums
01538     std::vector< int > dColumnsUnique;
01539     std::vector< double > dColumnSums;
01540 
01541     int nColumns = m_mapRemap.GetColumns();
01542     m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
01543     dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
01544     dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
01545 
01546     for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
01547     {
01548         dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
01549 
01550         assert( dataCols[i] < m_nTotDofs_SrcCov );
01551 
01552         // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
01553         int colGID = this->GetColGlobalDoF( dataCols[i] );  // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
01554         // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
01555         dColumnsUnique[dataCols[i]] = colGID;
01556 
01557         // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
01558         // std::endl;
01559     }
01560 
01561     int rootProc = 0;
01562     std::vector< int > nElementsInProc;
01563     const int nDATA = 3;
01564     if( !rank ) nElementsInProc.resize( size * nDATA );
01565     int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
01566     ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
01567     if( ierr != MPI_SUCCESS ) return -1;
01568 
01569     int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
01570     std::vector< int > dColumnIndices;
01571     std::vector< double > dColumnSourceAreas;
01572     std::vector< double > dColumnSumsTotal;
01573     std::vector< int > displs, rcount;
01574     if( rank == rootProc )
01575     {
01576         displs.resize( size + 1, 0 );
01577         rcount.resize( size, 0 );
01578         int gsum = 0;
01579         for( int ir = 0; ir < size; ++ir )
01580         {
01581             nTotVals += nElementsInProc[ir * nDATA];
01582             nTotColumns += nElementsInProc[ir * nDATA + 1];
01583             nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
01584 
01585             displs[ir] = gsum;
01586             rcount[ir] = nElementsInProc[ir * nDATA + 1];
01587             gsum += rcount[ir];
01588 
01589             // printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns,
01590             // displs[ir], rcount[ir], gsum);
01591         }
01592 
01593         dColumnIndices.resize( nTotColumns, -1 );
01594         dColumnSumsTotal.resize( nTotColumns, 0.0 );
01595         // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
01596     }
01597 
01598     // Gather all ColumnSums to root process and accumulate
01599     // We expect that the sums of all columns equate to 1.0 within user specified tolerance
01600     // Need to do a gatherv here since different processes have different number of elements
01601     // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
01602     // MPI_SUM, 0, m_pcomm->comm());
01603     ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
01604                         displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
01605     if( ierr != MPI_SUCCESS ) return -1;
01606     ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
01607                         displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
01608     if( ierr != MPI_SUCCESS ) return -1;
01609     // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
01610     // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
01611     // MPI_SUCCESS ) return -1;
01612 
01613     // Clean out unwanted arrays now
01614     dColumnSums.clear();
01615     dColumnsUnique.clear();
01616 
01617     // Verify all column sums equal the input Jacobian
01618     int fConservative = 0;
01619     if( rank == rootProc )
01620     {
01621         displs[size] = ( nTotColumns );
01622         // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
01623         std::map< int, double > dColumnSumsOnRoot;
01624         // std::map<int, double> dColumnSourceAreasOnRoot;
01625         for( int ir = 0; ir < size; ir++ )
01626         {
01627             for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
01628             {
01629                 if( dColumnIndices[ips] < 0 ) continue;
01630                 // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
01631                 assert( dColumnIndices[ips] < nTotColumnsUnq );
01632                 dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips];  // / dColumnSourceAreas[ips];
01633                 // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
01634                 // dColumnSourceAreas[ dColumnIndices[ips] ]
01635             }
01636         }
01637 
01638         for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
01639         {
01640             // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
01641             if( fabs( it->second - 1.0 ) > dTolerance )
01642             {
01643                 fConservative++;
01644                 Announce( "TempestOnlineMap is not conservative in column "
01645                           // "%i (%1.15e)", it->first, it->second );
01646                           "%i (%1.15e)",
01647                           it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
01648             }
01649         }
01650     }
01651 
01652     // TODO: Just do a broadcast from root instead of a reduction
01653     ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
01654     if( ierr != MPI_SUCCESS ) return -1;
01655 
01656     return fConservative;
01657 #endif
01658 }
01659 
01660 ///////////////////////////////////////////////////////////////////////////////
01661 
01662 int moab::TempestOnlineMap::IsMonotone( double dTolerance )
01663 {
01664 #ifndef MOAB_HAVE_MPI
01665 
01666     return OfflineMap::IsMonotone( dTolerance );
01667 
01668 #else
01669 
01670     // Get map entries
01671     DataArray1D< int > dataRows;
01672     DataArray1D< int > dataCols;
01673     DataArray1D< double > dataEntries;
01674 
01675     m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
01676 
01677     // Verify all entries are in the range [0,1]
01678     int fMonotone = 0;
01679     for( unsigned i = 0; i < dataRows.GetRows(); i++ )
01680     {
01681         if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
01682         {
01683             fMonotone++;
01684 
01685             Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
01686         }
01687     }
01688 
01689     int ierr;
01690     int fMonotoneGlobal = 0;
01691     ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
01692     if( ierr != MPI_SUCCESS ) return -1;
01693 
01694     return fMonotoneGlobal;
01695 #endif
01696 }
01697 
01698 ///////////////////////////////////////////////////////////////////////////////
01699 
01700 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( moab::Tag srcSolutionTag,
01701                                                       moab::Tag tgtSolutionTag,
01702                                                       bool transpose )
01703 {
01704     moab::ErrorCode rval;
01705 
01706     std::vector< double > solSTagVals;
01707     std::vector< double > solTTagVals;
01708 
01709     moab::Range sents, tents;
01710     if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
01711     {
01712         if( m_remapper->point_cloud_source )
01713         {
01714             moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
01715             solSTagVals.resize( covSrcEnts.size(), -1.0 );
01716             sents = covSrcEnts;
01717         }
01718         else
01719         {
01720             moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
01721             solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
01722                                 -1.0 );
01723             sents = covSrcEnts;
01724         }
01725         if( m_remapper->point_cloud_target )
01726         {
01727             moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
01728             solTTagVals.resize( tgtEnts.size(), -1.0 );
01729             tents = tgtEnts;
01730         }
01731         else
01732         {
01733             moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
01734             solTTagVals.resize(
01735                 tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
01736             tents = tgtEnts;
01737         }
01738     }
01739     else
01740     {
01741         moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
01742         moab::Range& tgtEnts    = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
01743         solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
01744                             -1.0 );
01745         solTTagVals.resize(
01746             tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
01747 
01748         sents = covSrcEnts;
01749         tents = tgtEnts;
01750     }
01751 
01752     // The tag data is np*np*n_el_src
01753     rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );
01754 
01755     // Compute the application of weights on the suorce solution data and store it in the
01756     // destination solution vector data Optionally, can also perform the transpose application of
01757     // the weight matrix. Set the 3rd argument to true if this is needed
01758     rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );
01759 
01760     // The tag data is np*np*n_el_dest
01761     rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
01762 
01763     return moab::MB_SUCCESS;
01764 }
01765 
01766 moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution( moab::Tag& solnTag,
01767                                                                   const std::string& solnName,
01768                                                                   moab::Remapper::IntersectionContext ctx,
01769                                                                   sample_function testFunction,
01770                                                                   moab::Tag* clonedSolnTag,
01771                                                                   std::string cloneSolnName )
01772 {
01773     moab::ErrorCode rval;
01774     const bool outputEnabled = ( is_root );
01775     int discOrder;
01776     DiscretizationType discMethod;
01777     moab::EntityHandle meshset;
01778     moab::Range entities;
01779     Mesh* trmesh;
01780     switch( ctx )
01781     {
01782         case Remapper::SourceMesh:
01783             meshset    = m_remapper->m_covering_source_set;
01784             trmesh     = m_remapper->m_covering_source;
01785             entities   = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
01786                                                           : m_remapper->m_covering_source_entities );
01787             discOrder  = m_nDofsPEl_Src;
01788             discMethod = m_eInputType;
01789             break;
01790 
01791         case Remapper::TargetMesh:
01792             meshset = m_remapper->m_target_set;
01793             trmesh  = m_remapper->m_target;
01794             entities =
01795                 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
01796             discOrder  = m_nDofsPEl_Dest;
01797             discMethod = m_eOutputType;
01798             break;
01799 
01800         default:
01801             if( outputEnabled )
01802                 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
01803             return moab::MB_FAILURE;
01804     }
01805 
01806     // Let us create teh solution tag with appropriate information for name, discretization order
01807     // (DoF space)
01808     rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
01809                                         MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
01810     if( clonedSolnTag != NULL )
01811     {
01812         if( cloneSolnName.size() == 0 )
01813         {
01814             cloneSolnName = solnName + std::string( "Cloned" );
01815         }
01816         rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
01817                                             *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
01818     }
01819 
01820     // Triangular quadrature rule
01821     const int TriQuadratureOrder = 10;
01822 
01823     if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
01824 
01825     TriangularQuadratureRule triquadrule( TriQuadratureOrder );
01826 
01827     const int TriQuadraturePoints = triquadrule.GetPoints();
01828 
01829     const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
01830     const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
01831 
01832     // Output data
01833     DataArray1D< double > dVar;
01834     DataArray1D< double > dVarMB;  // re-arranged local MOAB vector
01835 
01836     // Nodal geometric area
01837     DataArray1D< double > dNodeArea;
01838 
01839     // Calculate element areas
01840     // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
01841 
01842     if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
01843     {
01844         /* Get the spectral points and sample the functionals accordingly */
01845         const bool fGLL          = true;
01846         const bool fGLLIntegrate = false;
01847 
01848         // Generate grid metadata
01849         DataArray3D< int > dataGLLNodes;
01850         DataArray3D< double > dataGLLJacobian;
01851 
01852         GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
01853 
01854         // Number of elements
01855         int nElements = trmesh->faces.size();
01856 
01857         // Verify all elements are quadrilaterals
01858         for( int k = 0; k < nElements; k++ )
01859         {
01860             const Face& face = trmesh->faces[k];
01861 
01862             if( face.edges.size() != 4 )
01863             {
01864                 _EXCEPTIONT( "Non-quadrilateral face detected; "
01865                              "incompatible with --gll" );
01866             }
01867         }
01868 
01869         // Number of unique nodes
01870         int iMaxNode = 0;
01871         for( int i = 0; i < discOrder; i++ )
01872         {
01873             for( int j = 0; j < discOrder; j++ )
01874             {
01875                 for( int k = 0; k < nElements; k++ )
01876                 {
01877                     if( dataGLLNodes[i][j][k] > iMaxNode )
01878                     {
01879                         iMaxNode = dataGLLNodes[i][j][k];
01880                     }
01881                 }
01882             }
01883         }
01884 
01885         // Get Gauss-Lobatto quadrature nodes
01886         DataArray1D< double > dG;
01887         DataArray1D< double > dW;
01888 
01889         GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
01890 
01891         // Get Gauss quadrature nodes
01892         const int nGaussP = 10;
01893 
01894         DataArray1D< double > dGaussG;
01895         DataArray1D< double > dGaussW;
01896 
01897         GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
01898 
01899         // Allocate data
01900         dVar.Allocate( iMaxNode );
01901         dVarMB.Allocate( discOrder * discOrder * nElements );
01902         dNodeArea.Allocate( iMaxNode );
01903 
01904         // Sample data
01905         for( int k = 0; k < nElements; k++ )
01906         {
01907             const Face& face = trmesh->faces[k];
01908 
01909             // Sample data at GLL nodes
01910             if( fGLL )
01911             {
01912                 for( int i = 0; i < discOrder; i++ )
01913                 {
01914                     for( int j = 0; j < discOrder; j++ )
01915                     {
01916 
01917                         // Apply local map
01918                         Node node;
01919                         Node dDx1G;
01920                         Node dDx2G;
01921 
01922                         ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
01923 
01924                         // Sample data at this point
01925                         double dNodeLon = atan2( node.y, node.x );
01926                         if( dNodeLon < 0.0 )
01927                         {
01928                             dNodeLon += 2.0 * M_PI;
01929                         }
01930                         double dNodeLat = asin( node.z );
01931 
01932                         double dSample = ( *testFunction )( dNodeLon, dNodeLat );
01933 
01934                         dVar[dataGLLNodes[j][i][k] - 1] = dSample;
01935                     }
01936                 }
01937                 // High-order Gaussian integration over basis function
01938             }
01939             else
01940             {
01941                 DataArray2D< double > dCoeff( discOrder, discOrder );
01942 
01943                 for( int p = 0; p < nGaussP; p++ )
01944                 {
01945                     for( int q = 0; q < nGaussP; q++ )
01946                     {
01947 
01948                         // Apply local map
01949                         Node node;
01950                         Node dDx1G;
01951                         Node dDx2G;
01952 
01953                         ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
01954 
01955                         // Cross product gives local Jacobian
01956                         Node nodeCross = CrossProduct( dDx1G, dDx2G );
01957 
01958                         double dJacobian =
01959                             sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
01960 
01961                         // Find components of quadrature point in basis
01962                         // of the first Face
01963                         SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
01964 
01965                         // Sample data at this point
01966                         double dNodeLon = atan2( node.y, node.x );
01967                         if( dNodeLon < 0.0 )
01968                         {
01969                             dNodeLon += 2.0 * M_PI;
01970                         }
01971                         double dNodeLat = asin( node.z );
01972 
01973                         double dSample = ( *testFunction )( dNodeLon, dNodeLat );
01974 
01975                         // Integrate
01976                         for( int i = 0; i < discOrder; i++ )
01977                         {
01978                             for( int j = 0; j < discOrder; j++ )
01979                             {
01980 
01981                                 double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
01982 
01983                                 dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
01984 
01985                                 dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
01986                             }
01987                         }
01988                     }
01989                 }
01990             }
01991         }
01992 
01993         // Divide by area
01994         if( fGLLIntegrate )
01995         {
01996             for( size_t i = 0; i < dVar.GetRows(); i++ )
01997             {
01998                 dVar[i] /= dNodeArea[i];
01999             }
02000         }
02001 
02002         // Let us rearrange the data based on DoF ID specification
02003         if( ctx == Remapper::SourceMesh )
02004         {
02005             for( unsigned j = 0; j < entities.size(); j++ )
02006                 for( int p = 0; p < discOrder; p++ )
02007                     for( int q = 0; q < discOrder; q++ )
02008                     {
02009                         const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
02010                         dVarMB[offsetDOF]   = dVar[col_dtoc_dofmap[offsetDOF]];
02011                     }
02012         }
02013         else
02014         {
02015             for( unsigned j = 0; j < entities.size(); j++ )
02016                 for( int p = 0; p < discOrder; p++ )
02017                     for( int q = 0; q < discOrder; q++ )
02018                     {
02019                         const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
02020                         dVarMB[offsetDOF]   = dVar[row_dtoc_dofmap[offsetDOF]];
02021                     }
02022         }
02023 
02024         // Set the tag data
02025         rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
02026     }
02027     else
02028     {
02029         // assert( discOrder == 1 );
02030         if( discMethod == DiscretizationType_FV )
02031         {
02032             /* Compute an element-wise integral to store the sampled solution based on Quadrature
02033              * rules */
02034             // Resize the array
02035             dVar.Allocate( trmesh->faces.size() );
02036 
02037             std::vector< Node >& nodes = trmesh->nodes;
02038 
02039             // Loop through all Faces
02040             for( size_t i = 0; i < trmesh->faces.size(); i++ )
02041             {
02042                 const Face& face = trmesh->faces[i];
02043 
02044                 // Loop through all sub-triangles
02045                 for( size_t j = 0; j < face.edges.size() - 2; j++ )
02046                 {
02047 
02048                     const Node& node0 = nodes[face[0]];
02049                     const Node& node1 = nodes[face[j + 1]];
02050                     const Node& node2 = nodes[face[j + 2]];
02051 
02052                     // Triangle area
02053                     Face faceTri( 3 );
02054                     faceTri.SetNode( 0, face[0] );
02055                     faceTri.SetNode( 1, face[j + 1] );
02056                     faceTri.SetNode( 2, face[j + 2] );
02057 
02058                     double dTriangleArea = CalculateFaceArea( faceTri, nodes );
02059 
02060                     // Calculate the element average
02061                     double dTotalSample = 0.0;
02062 
02063                     // Loop through all quadrature points
02064                     for( int k = 0; k < TriQuadraturePoints; k++ )
02065                     {
02066                         Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
02067                                        TriQuadratureG[k][2] * node2.x,
02068                                    TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
02069                                        TriQuadratureG[k][2] * node2.y,
02070                                    TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
02071                                        TriQuadratureG[k][2] * node2.z );
02072 
02073                         double dMagnitude = node.Magnitude();
02074                         node.x /= dMagnitude;
02075                         node.y /= dMagnitude;
02076                         node.z /= dMagnitude;
02077 
02078                         double dLon = atan2( node.y, node.x );
02079                         if( dLon < 0.0 )
02080                         {
02081                             dLon += 2.0 * M_PI;
02082                         }
02083                         double dLat = asin( node.z );
02084 
02085                         double dSample = ( *testFunction )( dLon, dLat );
02086 
02087                         dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
02088                     }
02089 
02090                     dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
02091                 }
02092             }
02093             rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
02094         }
02095         else /* discMethod == DiscretizationType_PCLOUD */
02096         {
02097             /* Get the coordinates of the vertices and sample the functionals accordingly */
02098             std::vector< Node >& nodes = trmesh->nodes;
02099 
02100             // Resize the array
02101             dVar.Allocate( nodes.size() );
02102 
02103             for( size_t j = 0; j < nodes.size(); j++ )
02104             {
02105                 Node& node        = nodes[j];
02106                 double dMagnitude = node.Magnitude();
02107                 node.x /= dMagnitude;
02108                 node.y /= dMagnitude;
02109                 node.z /= dMagnitude;
02110                 double dLon = atan2( node.y, node.x );
02111                 if( dLon < 0.0 )
02112                 {
02113                     dLon += 2.0 * M_PI;
02114                 }
02115                 double dLat = asin( node.z );
02116 
02117                 double dSample = ( *testFunction )( dLon, dLat );
02118                 dVar[j]        = dSample;
02119             }
02120 
02121             rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
02122         }
02123     }
02124 
02125     return moab::MB_SUCCESS;
02126 }
02127 
02128 moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics( moab::Remapper::IntersectionContext ctx,
02129                                                         moab::Tag& exactTag,
02130                                                         moab::Tag& approxTag,
02131                                                         std::map< std::string, double >& metrics,
02132                                                         bool verbose )
02133 {
02134     moab::ErrorCode rval;
02135     const bool outputEnabled = ( is_root );
02136     int discOrder;
02137     DiscretizationType discMethod;
02138     moab::EntityHandle meshset;
02139     moab::Range entities;
02140     Mesh* trmesh;
02141     switch( ctx )
02142     {
02143         case Remapper::SourceMesh:
02144             meshset    = m_remapper->m_covering_source_set;
02145             trmesh     = m_remapper->m_covering_source;
02146             entities   = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
02147                                                           : m_remapper->m_covering_source_entities );
02148             discOrder  = m_nDofsPEl_Src;
02149             discMethod = m_eInputType;
02150             break;
02151 
02152         case Remapper::TargetMesh:
02153             meshset = m_remapper->m_target_set;
02154             trmesh  = m_remapper->m_target;
02155             entities =
02156                 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
02157             discOrder  = m_nDofsPEl_Dest;
02158             discMethod = m_eOutputType;
02159             break;
02160 
02161         default:
02162             if( outputEnabled )
02163                 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
02164             return moab::MB_FAILURE;
02165     }
02166 
02167     // Let us create teh solution tag with appropriate information for name, discretization order
02168     // (DoF space)
02169     std::string exactTagName, projTagName;
02170     const int ntotsize = entities.size() * discOrder * discOrder;
02171     int ntotsize_glob  = 0;
02172     std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
02173     rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
02174     rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
02175     rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
02176     rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
02177 
02178     std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 );  //  L1Err, L2Err, LinfErr
02179     for( int i = 0; i < ntotsize; ++i )
02180     {
02181         const double error = fabs( exactSolution[i] - projSolution[i] );
02182         errnorms[0] += error;
02183         errnorms[1] += error * error;
02184         errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] );
02185     }
02186 #ifdef MOAB_HAVE_MPI
02187     if( m_pcomm )
02188     {
02189         MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
02190         MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
02191         MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
02192     }
02193 #else
02194     ntotsize_glob = ntotsize;
02195     globerrnorms = errnorms;
02196 #endif
02197     globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob );
02198     globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob );
02199 
02200     metrics.clear();
02201     metrics["L1Error"]   = globerrnorms[0];
02202     metrics["L2Error"]   = globerrnorms[1];
02203     metrics["LinfError"] = globerrnorms[2];
02204 
02205     if( verbose && is_root )
02206     {
02207         std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
02208         std::cout << "\t L_1 error   = " << globerrnorms[0] << std::endl;
02209         std::cout << "\t L_2 error   = " << globerrnorms[1] << std::endl;
02210         std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl;
02211     }
02212 
02213     return moab::MB_SUCCESS;
02214 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines