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