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