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