|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
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 /// the tag should be created already in the e3sm workflow; if not, create it here 00938 Tag areaTag; 00939 rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag, 00940 MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT ); 00941 if ( MB_ALREADY_ALLOCATED == rval ) 00942 { 00943 if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" ); 00944 } 00945 00946 double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0; 00947 if( !m_bPointCloudSource ) 00948 { 00949 // Calculate Input Mesh Face areas 00950 if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" ); 00951 double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave ); 00952 rval = m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea );MB_CHK_ERR( rval ); 00953 dTotalAreaInput = dTotalAreaInput_loc; 00954 #ifdef MOAB_HAVE_MPI 00955 if( m_pcomm ) 00956 MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 00957 #endif 00958 if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput ); 00959 00960 // Input mesh areas 00961 m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave ); 00962 // we do not need to set the area on coverage mesh, only on source and target meshes 00963 // rval = m_interface->tag_set_data( areaTag, m_remapper->m_covering_source_entities, m_meshInputCov->vecFaceArea );MB_CHK_ERR( rval ); 00964 } 00965 00966 if( !m_bPointCloudTarget ) 00967 { 00968 // Calculate Output Mesh Face areas 00969 if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" ); 00970 double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave ); 00971 dTotalAreaOutput = dTotalAreaOutput_loc; 00972 #ifdef MOAB_HAVE_MPI 00973 if( m_pcomm ) 00974 MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 00975 #endif 00976 if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput ); 00977 rval = m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea );MB_CHK_ERR( rval ); 00978 } 00979 00980 if( !m_bPointCloud ) 00981 { 00982 // Verify that overlap mesh is in the correct order, only if size == 1 00983 if( 1 == size ) 00984 { 00985 int ixSourceFaceMax = ( -1 ); 00986 int ixTargetFaceMax = ( -1 ); 00987 00988 if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() ) 00989 { 00990 _EXCEPTIONT( "Invalid overlap mesh:\n" 00991 " Possible mesh file corruption?" ); 00992 } 00993 00994 for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ ) 00995 { 00996 if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax ) 00997 ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1; 00998 00999 if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax ) 01000 ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1; 01001 } 01002 01003 // Check for forward correspondence in overlap mesh 01004 if( m_meshInputCov->faces.size() - ixSourceFaceMax == 0 ) 01005 { 01006 if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" ); 01007 } 01008 else if( m_meshOutput->faces.size() - ixSourceFaceMax == 0 ) 01009 { // Check for reverse correspondence in overlap mesh 01010 if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" ); 01011 01012 // Reorder overlap mesh 01013 m_meshOverlap->ExchangeFirstAndSecondMesh(); 01014 } 01015 else 01016 { // No correspondence found 01017 _EXCEPTION4( "Invalid overlap mesh:\n No correspondence found with input and output meshes " 01018 "(%i,%i) vs (%i,%i)", 01019 m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax, 01020 ixTargetFaceMax ); 01021 } 01022 } 01023 01024 // Calculate Face areas 01025 if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" ); 01026 double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false ); 01027 double dTotalAreaOverlap = dTotalAreaOverlap_loc; 01028 #ifdef MOAB_HAVE_MPI 01029 if( m_pcomm ) 01030 MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 01031 #endif 01032 if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap ); 01033 01034 // Correct areas to match the areas calculated in the overlap mesh 01035 // if (fCorrectAreas) // In MOAB-TempestRemap, we will always keep this to be true 01036 { 01037 if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" ); 01038 DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() ); 01039 DataArray1D< double > dTargetArea( m_meshOutput->faces.size() ); 01040 01041 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() ); 01042 assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() ); 01043 assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() ); 01044 01045 assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() ); 01046 assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() ); 01047 01048 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ ) 01049 { 01050 if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 ) 01051 continue; // skip this cell since it is ghosted 01052 01053 // let us recompute the source/target areas based on overlap mesh areas 01054 assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() ); 01055 dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i]; 01056 assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() ); 01057 dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i]; 01058 } 01059 01060 for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ ) 01061 { 01062 if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 ) 01063 { 01064 m_meshInputCov->vecFaceArea[i] = dSourceArea[i]; 01065 } 01066 } 01067 for( size_t i = 0; i < m_meshOutput->faces.size(); i++ ) 01068 { 01069 if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 ) 01070 { 01071 m_meshOutput->vecFaceArea[i] = dTargetArea[i]; 01072 } 01073 } 01074 } 01075 01076 // Set source mesh areas in map 01077 if( !m_bPointCloudSource && eInputType == DiscretizationType_FV ) 01078 { 01079 this->SetSourceAreas( m_meshInputCov->vecFaceArea ); 01080 if( m_meshInputCov->vecMask.size() ) 01081 { 01082 this->SetSourceMask( m_meshInputCov->vecMask ); 01083 } 01084 } 01085 01086 // Set target mesh areas in map 01087 if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV ) 01088 { 01089 this->SetTargetAreas( m_meshOutput->vecFaceArea ); 01090 if( m_meshOutput->vecMask.size() ) 01091 { 01092 this->SetTargetMask( m_meshOutput->vecMask ); 01093 } 01094 } 01095 01096 /* 01097 // Recalculate input mesh area from overlap mesh 01098 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) { 01099 dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n"); 01100 dbgprint.printf(0, "Recalculating source mesh areas\n"); 01101 dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap); 01102 dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput); 01103 } 01104 */ 01105 } 01106 01107 // Finite volume input / Finite volume output 01108 if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) ) 01109 { 01110 // Generate reverse node array and edge map 01111 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray(); 01112 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false ); 01113 01114 // Initialize coordinates for map 01115 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov ); 01116 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput ); 01117 01118 // Finite volume input / Finite element output 01119 rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, false, NULL );MB_CHK_ERR( rval ); 01120 01121 // Construct remap for FV-FV 01122 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" ); 01123 01124 // Construct OfflineMap 01125 if( strMapAlgorithm == "invdist" ) 01126 { 01127 if( is_root ) AnnounceStartBlock( "Calculating map (invdist)" ); 01128 LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this ); 01129 } 01130 else if( strMapAlgorithm == "delaunay" ) 01131 { 01132 if( is_root ) AnnounceStartBlock( "Calculating map (delaunay)" ); 01133 LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this ); 01134 } 01135 else if( strMapAlgorithm == "fvintbilin" ) 01136 { 01137 if( is_root ) AnnounceStartBlock( "Calculating map (intbilin)" ); 01138 LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this ); 01139 } 01140 else if( strMapAlgorithm == "fvintbilingb" ) 01141 { 01142 if( is_root ) AnnounceStartBlock( "Calculating map (intbilingb)" ); 01143 LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this ); 01144 } 01145 else if( strMapAlgorithm == "fvbilin" ) 01146 { 01147 if( is_root ) AnnounceStartBlock( "Calculating map (bilin)" ); 01148 LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this ); 01149 } 01150 else 01151 { 01152 if( is_root ) AnnounceStartBlock( "Calculating map (default)" ); 01153 // LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, 01154 // ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this ); 01155 LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) ); 01156 } 01157 } 01158 else if( eInputType == DiscretizationType_FV ) 01159 { 01160 DataArray3D< double > dataGLLJacobian; 01161 01162 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" ); 01163 double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, 01164 dataGLLNodesDest, dataGLLJacobian ); 01165 01166 double dNumericalArea = dNumericalArea_loc; 01167 #ifdef MOAB_HAVE_MPI 01168 if( m_pcomm ) 01169 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 01170 #endif 01171 if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea ); 01172 01173 // Initialize coordinates for map 01174 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov ); 01175 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest ); 01176 01177 // Generate the continuous Jacobian 01178 bool fContinuous = ( eOutputType == DiscretizationType_CGLL ); 01179 01180 if( eOutputType == DiscretizationType_CGLL ) 01181 { 01182 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() ); 01183 } 01184 else 01185 { 01186 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() ); 01187 } 01188 01189 // Generate reverse node array and edge map 01190 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray(); 01191 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false ); 01192 01193 // Finite volume input / Finite element output 01194 rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, 01195 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval ); 01196 01197 // Generate remap weights 01198 if( strMapAlgorithm == "volumetric" ) 01199 { 01200 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" ); 01201 LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, 01202 dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this, 01203 nMonotoneType, fContinuous, mapOptions.fNoConservation ); 01204 } 01205 else 01206 { 01207 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" ); 01208 LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian, 01209 this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous, 01210 mapOptions.fNoConservation ); 01211 } 01212 } 01213 else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) ) 01214 { 01215 DataArray3D< double > dataGLLJacobian; 01216 if( !m_bPointCloudSource ) 01217 { 01218 // Generate reverse node array and edge map 01219 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray(); 01220 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false ); 01221 01222 // Initialize coordinates for map 01223 if( eInputType == DiscretizationType_FV ) 01224 { 01225 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov ); 01226 } 01227 else 01228 { 01229 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" ); 01230 DataArray3D< double > dataGLLJacobianSrc; 01231 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov, 01232 dataGLLJacobian ); 01233 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc, 01234 dataGLLJacobianSrc ); 01235 } 01236 } 01237 // else { /* Source is a point cloud dataset */ } 01238 01239 if( !m_bPointCloudTarget ) 01240 { 01241 // Generate reverse node array and edge map 01242 if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray(); 01243 if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false ); 01244 01245 // Initialize coordinates for map 01246 if( eOutputType == DiscretizationType_FV ) 01247 { 01248 this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput ); 01249 } 01250 else 01251 { 01252 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" ); 01253 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest, 01254 dataGLLJacobian ); 01255 } 01256 } 01257 // else { /* Target is a point cloud dataset */ } 01258 01259 // Finite volume input / Finite element output 01260 rval = this->SetDOFmapAssociation( 01261 eInputType, ( eInputType == DiscretizationType_CGLL ), 01262 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrcCov ), 01263 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrc ), eOutputType, 01264 ( eOutputType == DiscretizationType_CGLL ), ( m_bPointCloudTarget ? NULL : &dataGLLNodesDest ) );MB_CHK_ERR( rval ); 01265 01266 // Construct remap 01267 if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" ); 01268 rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval ); 01269 } 01270 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) ) 01271 { 01272 DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian; 01273 01274 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" ); 01275 // double dNumericalAreaCov_loc = 01276 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov, 01277 dataGLLJacobian ); 01278 01279 double dNumericalArea_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, 01280 dataGLLNodesSrc, dataGLLJacobianSrc ); 01281 01282 // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area: 01283 // %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc ); 01284 // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc); 01285 01286 double dNumericalArea = dNumericalArea_loc; 01287 #ifdef MOAB_HAVE_MPI 01288 if( m_pcomm ) 01289 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 01290 #endif 01291 if( is_root ) 01292 { 01293 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea ); 01294 01295 if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 ) 01296 { 01297 dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh " 01298 "numerical area and geometric area\n" ); 01299 } 01300 } 01301 01302 if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() ) 01303 { 01304 _EXCEPTIONT( "Number of element does not match between metadata and " 01305 "input mesh" ); 01306 } 01307 01308 // Initialize coordinates for map 01309 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov ); 01310 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput ); 01311 01312 // Generate the continuous Jacobian for input mesh 01313 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL ); 01314 01315 if( eInputType == DiscretizationType_CGLL ) 01316 { 01317 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() ); 01318 } 01319 else 01320 { 01321 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() ); 01322 } 01323 01324 // Finite element input / Finite volume output 01325 rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ), 01326 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, false, NULL );MB_CHK_ERR( rval ); 01327 01328 // Generate remap 01329 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" ); 01330 01331 if( strMapAlgorithm == "volumetric" ) 01332 { 01333 _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for" 01334 "GLL input mesh" ); 01335 } 01336 01337 LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn, 01338 mapOptions.fNoConservation ); 01339 } 01340 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) ) 01341 { 01342 DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc; 01343 DataArray3D< double > dataGLLJacobianOut; 01344 01345 // Input metadata 01346 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" ); 01347 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov, 01348 dataGLLJacobianIn ); 01349 01350 double dNumericalAreaSrc_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, 01351 dataGLLNodesSrc, dataGLLJacobianSrc ); 01352 double dNumericalAreaIn = dNumericalAreaSrc_loc; 01353 #ifdef MOAB_HAVE_MPI 01354 if( m_pcomm ) 01355 MPI_Reduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 01356 #endif 01357 if( is_root ) 01358 { 01359 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn ); 01360 01361 if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 ) 01362 { 01363 dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh " 01364 "numerical area and geometric area\n" ); 01365 } 01366 } 01367 01368 // Output metadata 01369 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" ); 01370 double dNumericalAreaOut_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, 01371 dataGLLNodesDest, dataGLLJacobianOut ); 01372 double dNumericalAreaOut = dNumericalAreaOut_loc; 01373 #ifdef MOAB_HAVE_MPI 01374 if( m_pcomm ) 01375 MPI_Reduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 01376 #endif 01377 if( is_root ) 01378 { 01379 dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut ); 01380 01381 if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 ) 01382 { 01383 if( is_root ) 01384 dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh " 01385 "numerical area and geometric area\n" ); 01386 } 01387 } 01388 01389 // Initialize coordinates for map 01390 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov ); 01391 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest ); 01392 01393 // Generate the continuous Jacobian for input mesh 01394 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL ); 01395 01396 if( eInputType == DiscretizationType_CGLL ) 01397 { 01398 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() ); 01399 } 01400 else 01401 { 01402 GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() ); 01403 } 01404 01405 // Generate the continuous Jacobian for output mesh 01406 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL ); 01407 01408 if( eOutputType == DiscretizationType_CGLL ) 01409 { 01410 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() ); 01411 } 01412 else 01413 { 01414 GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() ); 01415 } 01416 01417 // Input Finite Element to Output Finite Element 01418 rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ), 01419 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, 01420 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval ); 01421 01422 // Generate remap 01423 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" ); 01424 01425 LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut, 01426 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType, 01427 fContinuousIn, fContinuousOut, mapOptions.fNoConservation ); 01428 } 01429 else 01430 { 01431 _EXCEPTIONT( "Not implemented" ); 01432 } 01433 01434 #ifdef MOAB_HAVE_EIGEN3 01435 copy_tempest_sparsemat_to_eigen3(); 01436 #endif 01437 01438 #ifdef MOAB_HAVE_MPI 01439 { 01440 // Remove ghosted entities from overlap set 01441 moab::Range ghostedEnts; 01442 rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval ); 01443 moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh ); 01444 rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" ); 01445 } 01446 #endif 01447 // Verify consistency, conservation and monotonicity, globally 01448 if( !mapOptions.fNoCheck ) 01449 { 01450 if( is_root ) dbgprint.printf( 0, "Verifying map" ); 01451 this->IsConsistent( 1.0e-8 ); 01452 if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 ); 01453 01454 if( nMonotoneType != 0 ) 01455 { 01456 this->IsMonotone( 1.0e-12 ); 01457 } 01458 } 01459 } 01460 catch( Exception& e ) 01461 { 01462 dbgprint.printf( 0, "%s", e.ToString().c_str() ); 01463 return ( moab::MB_FAILURE ); 01464 } 01465 catch( ... ) 01466 { 01467 return ( moab::MB_FAILURE ); 01468 } 01469 return moab::MB_SUCCESS; 01470 } 01471 01472 /////////////////////////////////////////////////////////////////////////////// 01473 01474 int moab::TempestOnlineMap::IsConsistent( double dTolerance ) 01475 { 01476 #ifndef MOAB_HAVE_MPI 01477 01478 return OfflineMap::IsConsistent( dTolerance ); 01479 01480 #else 01481 01482 // Get map entries 01483 DataArray1D< int > dataRows; 01484 DataArray1D< int > dataCols; 01485 DataArray1D< double > dataEntries; 01486 01487 // Calculate row sums 01488 DataArray1D< double > dRowSums; 01489 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries ); 01490 dRowSums.Allocate( m_mapRemap.GetRows() ); 01491 01492 for( unsigned i = 0; i < dataRows.GetRows(); i++ ) 01493 { 01494 dRowSums[dataRows[i]] += dataEntries[i]; 01495 } 01496 01497 // Verify all row sums are equal to 1 01498 int fConsistent = 0; 01499 for( unsigned i = 0; i < dRowSums.GetRows(); i++ ) 01500 { 01501 if( fabs( dRowSums[i] - 1.0 ) > dTolerance ) 01502 { 01503 fConsistent++; 01504 int rowGID = row_gdofmap[i]; 01505 Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] ); 01506 } 01507 } 01508 01509 int ierr; 01510 int fConsistentGlobal = 0; 01511 ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() ); 01512 if( ierr != MPI_SUCCESS ) return -1; 01513 01514 return fConsistentGlobal; 01515 #endif 01516 } 01517 01518 /////////////////////////////////////////////////////////////////////////////// 01519 01520 int moab::TempestOnlineMap::IsConservative( double dTolerance ) 01521 { 01522 #ifndef MOAB_HAVE_MPI 01523 01524 return OfflineMap::IsConservative( dTolerance ); 01525 01526 #else 01527 // return OfflineMap::IsConservative(dTolerance); 01528 01529 int ierr; 01530 // Get map entries 01531 DataArray1D< int > dataRows; 01532 DataArray1D< int > dataCols; 01533 DataArray1D< double > dataEntries; 01534 const DataArray1D< double >& dTargetAreas = this->GetTargetAreas(); 01535 const DataArray1D< double >& dSourceAreas = this->GetSourceAreas(); 01536 01537 // Calculate column sums 01538 std::vector< int > dColumnsUnique; 01539 std::vector< double > dColumnSums; 01540 01541 int nColumns = m_mapRemap.GetColumns(); 01542 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries ); 01543 dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 ); 01544 dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 ); 01545 01546 for( unsigned i = 0; i < dataEntries.GetRows(); i++ ) 01547 { 01548 dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]]; 01549 01550 assert( dataCols[i] < m_nTotDofs_SrcCov ); 01551 01552 // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ] 01553 int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ]; 01554 // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ]; 01555 dColumnsUnique[dataCols[i]] = colGID; 01556 01557 // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID << 01558 // std::endl; 01559 } 01560 01561 int rootProc = 0; 01562 std::vector< int > nElementsInProc; 01563 const int nDATA = 3; 01564 if( !rank ) nElementsInProc.resize( size * nDATA ); 01565 int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src }; 01566 ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() ); 01567 if( ierr != MPI_SUCCESS ) return -1; 01568 01569 int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0; 01570 std::vector< int > dColumnIndices; 01571 std::vector< double > dColumnSourceAreas; 01572 std::vector< double > dColumnSumsTotal; 01573 std::vector< int > displs, rcount; 01574 if( rank == rootProc ) 01575 { 01576 displs.resize( size + 1, 0 ); 01577 rcount.resize( size, 0 ); 01578 int gsum = 0; 01579 for( int ir = 0; ir < size; ++ir ) 01580 { 01581 nTotVals += nElementsInProc[ir * nDATA]; 01582 nTotColumns += nElementsInProc[ir * nDATA + 1]; 01583 nTotColumnsUnq += nElementsInProc[ir * nDATA + 2]; 01584 01585 displs[ir] = gsum; 01586 rcount[ir] = nElementsInProc[ir * nDATA + 1]; 01587 gsum += rcount[ir]; 01588 01589 // printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns, 01590 // displs[ir], rcount[ir], gsum); 01591 } 01592 01593 dColumnIndices.resize( nTotColumns, -1 ); 01594 dColumnSumsTotal.resize( nTotColumns, 0.0 ); 01595 // dColumnSourceAreas.resize ( nTotColumns, 0.0 ); 01596 } 01597 01598 // Gather all ColumnSums to root process and accumulate 01599 // We expect that the sums of all columns equate to 1.0 within user specified tolerance 01600 // Need to do a gatherv here since different processes have different number of elements 01601 // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE, 01602 // MPI_SUM, 0, m_pcomm->comm()); 01603 ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(), 01604 displs.data(), MPI_INT, rootProc, m_pcomm->comm() ); 01605 if( ierr != MPI_SUCCESS ) return -1; 01606 ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(), 01607 displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); 01608 if( ierr != MPI_SUCCESS ) return -1; 01609 // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0], 01610 // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr != 01611 // MPI_SUCCESS ) return -1; 01612 01613 // Clean out unwanted arrays now 01614 dColumnSums.clear(); 01615 dColumnsUnique.clear(); 01616 01617 // Verify all column sums equal the input Jacobian 01618 int fConservative = 0; 01619 if( rank == rootProc ) 01620 { 01621 displs[size] = ( nTotColumns ); 01622 // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0); 01623 std::map< int, double > dColumnSumsOnRoot; 01624 // std::map<int, double> dColumnSourceAreasOnRoot; 01625 for( int ir = 0; ir < size; ir++ ) 01626 { 01627 for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ ) 01628 { 01629 if( dColumnIndices[ips] < 0 ) continue; 01630 // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]); 01631 assert( dColumnIndices[ips] < nTotColumnsUnq ); 01632 dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips]; 01633 // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips]; 01634 // dColumnSourceAreas[ dColumnIndices[ips] ] 01635 } 01636 } 01637 01638 for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it ) 01639 { 01640 // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance ) 01641 if( fabs( it->second - 1.0 ) > dTolerance ) 01642 { 01643 fConservative++; 01644 Announce( "TempestOnlineMap is not conservative in column " 01645 // "%i (%1.15e)", it->first, it->second ); 01646 "%i (%1.15e)", 01647 it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ ); 01648 } 01649 } 01650 } 01651 01652 // TODO: Just do a broadcast from root instead of a reduction 01653 ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() ); 01654 if( ierr != MPI_SUCCESS ) return -1; 01655 01656 return fConservative; 01657 #endif 01658 } 01659 01660 /////////////////////////////////////////////////////////////////////////////// 01661 01662 int moab::TempestOnlineMap::IsMonotone( double dTolerance ) 01663 { 01664 #ifndef MOAB_HAVE_MPI 01665 01666 return OfflineMap::IsMonotone( dTolerance ); 01667 01668 #else 01669 01670 // Get map entries 01671 DataArray1D< int > dataRows; 01672 DataArray1D< int > dataCols; 01673 DataArray1D< double > dataEntries; 01674 01675 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries ); 01676 01677 // Verify all entries are in the range [0,1] 01678 int fMonotone = 0; 01679 for( unsigned i = 0; i < dataRows.GetRows(); i++ ) 01680 { 01681 if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) ) 01682 { 01683 fMonotone++; 01684 01685 Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] ); 01686 } 01687 } 01688 01689 int ierr; 01690 int fMonotoneGlobal = 0; 01691 ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() ); 01692 if( ierr != MPI_SUCCESS ) return -1; 01693 01694 return fMonotoneGlobal; 01695 #endif 01696 } 01697 01698 /////////////////////////////////////////////////////////////////////////////// 01699 01700 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( moab::Tag srcSolutionTag, 01701 moab::Tag tgtSolutionTag, 01702 bool transpose ) 01703 { 01704 moab::ErrorCode rval; 01705 01706 std::vector< double > solSTagVals; 01707 std::vector< double > solTTagVals; 01708 01709 moab::Range sents, tents; 01710 if( m_remapper->point_cloud_source || m_remapper->point_cloud_target ) 01711 { 01712 if( m_remapper->point_cloud_source ) 01713 { 01714 moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh ); 01715 solSTagVals.resize( covSrcEnts.size(), -1.0 ); 01716 sents = covSrcEnts; 01717 } 01718 else 01719 { 01720 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh ); 01721 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(), 01722 -1.0 ); 01723 sents = covSrcEnts; 01724 } 01725 if( m_remapper->point_cloud_target ) 01726 { 01727 moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh ); 01728 solTTagVals.resize( tgtEnts.size(), -1.0 ); 01729 tents = tgtEnts; 01730 } 01731 else 01732 { 01733 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh ); 01734 solTTagVals.resize( 01735 tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 ); 01736 tents = tgtEnts; 01737 } 01738 } 01739 else 01740 { 01741 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh ); 01742 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh ); 01743 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(), 01744 -1.0 ); 01745 solTTagVals.resize( 01746 tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 ); 01747 01748 sents = covSrcEnts; 01749 tents = tgtEnts; 01750 } 01751 01752 // The tag data is np*np*n_el_src 01753 rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" ); 01754 01755 // Compute the application of weights on the suorce solution data and store it in the 01756 // destination solution vector data Optionally, can also perform the transpose application of 01757 // the weight matrix. Set the 3rd argument to true if this is needed 01758 rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" ); 01759 01760 // The tag data is np*np*n_el_dest 01761 rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" ); 01762 01763 return moab::MB_SUCCESS; 01764 } 01765 01766 moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution( moab::Tag& solnTag, 01767 const std::string& solnName, 01768 moab::Remapper::IntersectionContext ctx, 01769 sample_function testFunction, 01770 moab::Tag* clonedSolnTag, 01771 std::string cloneSolnName ) 01772 { 01773 moab::ErrorCode rval; 01774 const bool outputEnabled = ( is_root ); 01775 int discOrder; 01776 DiscretizationType discMethod; 01777 moab::EntityHandle meshset; 01778 moab::Range entities; 01779 Mesh* trmesh; 01780 switch( ctx ) 01781 { 01782 case Remapper::SourceMesh: 01783 meshset = m_remapper->m_covering_source_set; 01784 trmesh = m_remapper->m_covering_source; 01785 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices 01786 : m_remapper->m_covering_source_entities ); 01787 discOrder = m_nDofsPEl_Src; 01788 discMethod = m_eInputType; 01789 break; 01790 01791 case Remapper::TargetMesh: 01792 meshset = m_remapper->m_target_set; 01793 trmesh = m_remapper->m_target; 01794 entities = 01795 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities ); 01796 discOrder = m_nDofsPEl_Dest; 01797 discMethod = m_eOutputType; 01798 break; 01799 01800 default: 01801 if( outputEnabled ) 01802 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl; 01803 return moab::MB_FAILURE; 01804 } 01805 01806 // Let us create teh solution tag with appropriate information for name, discretization order 01807 // (DoF space) 01808 rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag, 01809 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval ); 01810 if( clonedSolnTag != NULL ) 01811 { 01812 if( cloneSolnName.size() == 0 ) 01813 { 01814 cloneSolnName = solnName + std::string( "Cloned" ); 01815 } 01816 rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, 01817 *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval ); 01818 } 01819 01820 // Triangular quadrature rule 01821 const int TriQuadratureOrder = 10; 01822 01823 if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl; 01824 01825 TriangularQuadratureRule triquadrule( TriQuadratureOrder ); 01826 01827 const int TriQuadraturePoints = triquadrule.GetPoints(); 01828 01829 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG(); 01830 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW(); 01831 01832 // Output data 01833 DataArray1D< double > dVar; 01834 DataArray1D< double > dVarMB; // re-arranged local MOAB vector 01835 01836 // Nodal geometric area 01837 DataArray1D< double > dNodeArea; 01838 01839 // Calculate element areas 01840 // trmesh->CalculateFaceAreas(fContainsConcaveFaces); 01841 01842 if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL ) 01843 { 01844 /* Get the spectral points and sample the functionals accordingly */ 01845 const bool fGLL = true; 01846 const bool fGLLIntegrate = false; 01847 01848 // Generate grid metadata 01849 DataArray3D< int > dataGLLNodes; 01850 DataArray3D< double > dataGLLJacobian; 01851 01852 GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian ); 01853 01854 // Number of elements 01855 int nElements = trmesh->faces.size(); 01856 01857 // Verify all elements are quadrilaterals 01858 for( int k = 0; k < nElements; k++ ) 01859 { 01860 const Face& face = trmesh->faces[k]; 01861 01862 if( face.edges.size() != 4 ) 01863 { 01864 _EXCEPTIONT( "Non-quadrilateral face detected; " 01865 "incompatible with --gll" ); 01866 } 01867 } 01868 01869 // Number of unique nodes 01870 int iMaxNode = 0; 01871 for( int i = 0; i < discOrder; i++ ) 01872 { 01873 for( int j = 0; j < discOrder; j++ ) 01874 { 01875 for( int k = 0; k < nElements; k++ ) 01876 { 01877 if( dataGLLNodes[i][j][k] > iMaxNode ) 01878 { 01879 iMaxNode = dataGLLNodes[i][j][k]; 01880 } 01881 } 01882 } 01883 } 01884 01885 // Get Gauss-Lobatto quadrature nodes 01886 DataArray1D< double > dG; 01887 DataArray1D< double > dW; 01888 01889 GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW ); 01890 01891 // Get Gauss quadrature nodes 01892 const int nGaussP = 10; 01893 01894 DataArray1D< double > dGaussG; 01895 DataArray1D< double > dGaussW; 01896 01897 GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW ); 01898 01899 // Allocate data 01900 dVar.Allocate( iMaxNode ); 01901 dVarMB.Allocate( discOrder * discOrder * nElements ); 01902 dNodeArea.Allocate( iMaxNode ); 01903 01904 // Sample data 01905 for( int k = 0; k < nElements; k++ ) 01906 { 01907 const Face& face = trmesh->faces[k]; 01908 01909 // Sample data at GLL nodes 01910 if( fGLL ) 01911 { 01912 for( int i = 0; i < discOrder; i++ ) 01913 { 01914 for( int j = 0; j < discOrder; j++ ) 01915 { 01916 01917 // Apply local map 01918 Node node; 01919 Node dDx1G; 01920 Node dDx2G; 01921 01922 ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G ); 01923 01924 // Sample data at this point 01925 double dNodeLon = atan2( node.y, node.x ); 01926 if( dNodeLon < 0.0 ) 01927 { 01928 dNodeLon += 2.0 * M_PI; 01929 } 01930 double dNodeLat = asin( node.z ); 01931 01932 double dSample = ( *testFunction )( dNodeLon, dNodeLat ); 01933 01934 dVar[dataGLLNodes[j][i][k] - 1] = dSample; 01935 } 01936 } 01937 // High-order Gaussian integration over basis function 01938 } 01939 else 01940 { 01941 DataArray2D< double > dCoeff( discOrder, discOrder ); 01942 01943 for( int p = 0; p < nGaussP; p++ ) 01944 { 01945 for( int q = 0; q < nGaussP; q++ ) 01946 { 01947 01948 // Apply local map 01949 Node node; 01950 Node dDx1G; 01951 Node dDx2G; 01952 01953 ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G ); 01954 01955 // Cross product gives local Jacobian 01956 Node nodeCross = CrossProduct( dDx1G, dDx2G ); 01957 01958 double dJacobian = 01959 sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z ); 01960 01961 // Find components of quadrature point in basis 01962 // of the first Face 01963 SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff ); 01964 01965 // Sample data at this point 01966 double dNodeLon = atan2( node.y, node.x ); 01967 if( dNodeLon < 0.0 ) 01968 { 01969 dNodeLon += 2.0 * M_PI; 01970 } 01971 double dNodeLat = asin( node.z ); 01972 01973 double dSample = ( *testFunction )( dNodeLon, dNodeLat ); 01974 01975 // Integrate 01976 for( int i = 0; i < discOrder; i++ ) 01977 { 01978 for( int j = 0; j < discOrder; j++ ) 01979 { 01980 01981 double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian; 01982 01983 dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea; 01984 01985 dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea; 01986 } 01987 } 01988 } 01989 } 01990 } 01991 } 01992 01993 // Divide by area 01994 if( fGLLIntegrate ) 01995 { 01996 for( size_t i = 0; i < dVar.GetRows(); i++ ) 01997 { 01998 dVar[i] /= dNodeArea[i]; 01999 } 02000 } 02001 02002 // Let us rearrange the data based on DoF ID specification 02003 if( ctx == Remapper::SourceMesh ) 02004 { 02005 for( unsigned j = 0; j < entities.size(); j++ ) 02006 for( int p = 0; p < discOrder; p++ ) 02007 for( int q = 0; q < discOrder; q++ ) 02008 { 02009 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q; 02010 dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]]; 02011 } 02012 } 02013 else 02014 { 02015 for( unsigned j = 0; j < entities.size(); j++ ) 02016 for( int p = 0; p < discOrder; p++ ) 02017 for( int q = 0; q < discOrder; q++ ) 02018 { 02019 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q; 02020 dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]]; 02021 } 02022 } 02023 02024 // Set the tag data 02025 rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval ); 02026 } 02027 else 02028 { 02029 // assert( discOrder == 1 ); 02030 if( discMethod == DiscretizationType_FV ) 02031 { 02032 /* Compute an element-wise integral to store the sampled solution based on Quadrature 02033 * rules */ 02034 // Resize the array 02035 dVar.Allocate( trmesh->faces.size() ); 02036 02037 std::vector< Node >& nodes = trmesh->nodes; 02038 02039 // Loop through all Faces 02040 for( size_t i = 0; i < trmesh->faces.size(); i++ ) 02041 { 02042 const Face& face = trmesh->faces[i]; 02043 02044 // Loop through all sub-triangles 02045 for( size_t j = 0; j < face.edges.size() - 2; j++ ) 02046 { 02047 02048 const Node& node0 = nodes[face[0]]; 02049 const Node& node1 = nodes[face[j + 1]]; 02050 const Node& node2 = nodes[face[j + 2]]; 02051 02052 // Triangle area 02053 Face faceTri( 3 ); 02054 faceTri.SetNode( 0, face[0] ); 02055 faceTri.SetNode( 1, face[j + 1] ); 02056 faceTri.SetNode( 2, face[j + 2] ); 02057 02058 double dTriangleArea = CalculateFaceArea( faceTri, nodes ); 02059 02060 // Calculate the element average 02061 double dTotalSample = 0.0; 02062 02063 // Loop through all quadrature points 02064 for( int k = 0; k < TriQuadraturePoints; k++ ) 02065 { 02066 Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x + 02067 TriQuadratureG[k][2] * node2.x, 02068 TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y + 02069 TriQuadratureG[k][2] * node2.y, 02070 TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z + 02071 TriQuadratureG[k][2] * node2.z ); 02072 02073 double dMagnitude = node.Magnitude(); 02074 node.x /= dMagnitude; 02075 node.y /= dMagnitude; 02076 node.z /= dMagnitude; 02077 02078 double dLon = atan2( node.y, node.x ); 02079 if( dLon < 0.0 ) 02080 { 02081 dLon += 2.0 * M_PI; 02082 } 02083 double dLat = asin( node.z ); 02084 02085 double dSample = ( *testFunction )( dLon, dLat ); 02086 02087 dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea; 02088 } 02089 02090 dVar[i] += dTotalSample / trmesh->vecFaceArea[i]; 02091 } 02092 } 02093 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval ); 02094 } 02095 else /* discMethod == DiscretizationType_PCLOUD */ 02096 { 02097 /* Get the coordinates of the vertices and sample the functionals accordingly */ 02098 std::vector< Node >& nodes = trmesh->nodes; 02099 02100 // Resize the array 02101 dVar.Allocate( nodes.size() ); 02102 02103 for( size_t j = 0; j < nodes.size(); j++ ) 02104 { 02105 Node& node = nodes[j]; 02106 double dMagnitude = node.Magnitude(); 02107 node.x /= dMagnitude; 02108 node.y /= dMagnitude; 02109 node.z /= dMagnitude; 02110 double dLon = atan2( node.y, node.x ); 02111 if( dLon < 0.0 ) 02112 { 02113 dLon += 2.0 * M_PI; 02114 } 02115 double dLat = asin( node.z ); 02116 02117 double dSample = ( *testFunction )( dLon, dLat ); 02118 dVar[j] = dSample; 02119 } 02120 02121 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval ); 02122 } 02123 } 02124 02125 return moab::MB_SUCCESS; 02126 } 02127 02128 moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics( moab::Remapper::IntersectionContext ctx, 02129 moab::Tag& exactTag, 02130 moab::Tag& approxTag, 02131 std::map< std::string, double >& metrics, 02132 bool verbose ) 02133 { 02134 moab::ErrorCode rval; 02135 const bool outputEnabled = ( is_root ); 02136 int discOrder; 02137 DiscretizationType discMethod; 02138 moab::EntityHandle meshset; 02139 moab::Range entities; 02140 Mesh* trmesh; 02141 switch( ctx ) 02142 { 02143 case Remapper::SourceMesh: 02144 meshset = m_remapper->m_covering_source_set; 02145 trmesh = m_remapper->m_covering_source; 02146 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices 02147 : m_remapper->m_covering_source_entities ); 02148 discOrder = m_nDofsPEl_Src; 02149 discMethod = m_eInputType; 02150 break; 02151 02152 case Remapper::TargetMesh: 02153 meshset = m_remapper->m_target_set; 02154 trmesh = m_remapper->m_target; 02155 entities = 02156 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities ); 02157 discOrder = m_nDofsPEl_Dest; 02158 discMethod = m_eOutputType; 02159 break; 02160 02161 default: 02162 if( outputEnabled ) 02163 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl; 02164 return moab::MB_FAILURE; 02165 } 02166 02167 // Let us create teh solution tag with appropriate information for name, discretization order 02168 // (DoF space) 02169 std::string exactTagName, projTagName; 02170 const int ntotsize = entities.size() * discOrder * discOrder; 02171 int ntotsize_glob = 0; 02172 std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 ); 02173 rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval ); 02174 rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval ); 02175 rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval ); 02176 rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval ); 02177 02178 std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 ); // L1Err, L2Err, LinfErr 02179 for( int i = 0; i < ntotsize; ++i ) 02180 { 02181 const double error = fabs( exactSolution[i] - projSolution[i] ); 02182 errnorms[0] += error; 02183 errnorms[1] += error * error; 02184 errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] ); 02185 } 02186 #ifdef MOAB_HAVE_MPI 02187 if( m_pcomm ) 02188 { 02189 MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() ); 02190 MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() ); 02191 MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() ); 02192 } 02193 #else 02194 ntotsize_glob = ntotsize; 02195 globerrnorms = errnorms; 02196 #endif 02197 globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob ); 02198 globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob ); 02199 02200 metrics.clear(); 02201 metrics["L1Error"] = globerrnorms[0]; 02202 metrics["L2Error"] = globerrnorms[1]; 02203 metrics["LinfError"] = globerrnorms[2]; 02204 02205 if( verbose && is_root ) 02206 { 02207 std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl; 02208 std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl; 02209 std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl; 02210 std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl; 02211 } 02212 02213 return moab::MB_SUCCESS; 02214 }
1.7.6.1 
