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