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