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