MOAB: Mesh Oriented datABase
(version 5.1.1)
|
00001 /* 00002 * ===================================================================================== 00003 * 00004 * Filename: TempestOfflineMap.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), [email protected] 00011 * 00012 * ===================================================================================== 00013 */ 00014 00015 #include "Announce.h" 00016 #include "DataMatrix3D.h" 00017 #include "FiniteElementTools.h" 00018 #include "SparseMatrix.h" 00019 #include "STLStringHelper.h" 00020 00021 #include "moab/Remapping/TempestOfflineMap.hpp" 00022 #include "DebugOutput.hpp" 00023 00024 #include <fstream> 00025 #include <cmath> 00026 #include <cstdlib> 00027 00028 /////////////////////////////////////////////////////////////////////////////// 00029 00030 // #define VERBOSE 00031 // #define VVERBOSE 00032 00033 /////////////////////////////////////////////////////////////////////////////// 00034 00035 moab::TempestOfflineMap::TempestOfflineMap ( moab::TempestRemapper* remapper ) : OfflineMap(), m_remapper ( remapper ) 00036 { 00037 // Get the references for the MOAB core objects 00038 mbCore = m_remapper->get_interface(); 00039 #ifdef MOAB_HAVE_MPI 00040 pcomm = m_remapper->get_parallel_communicator(); 00041 #endif 00042 00043 // Update the references to the meshes 00044 m_meshInput = remapper->GetMesh ( moab::Remapper::SourceMesh ); 00045 m_meshInputCov = remapper->GetCoveringMesh(); 00046 m_meshOutput = remapper->GetMesh ( moab::Remapper::TargetMesh ); 00047 m_meshOverlap = remapper->GetMesh ( moab::Remapper::IntersectedMesh ); 00048 m_globalMapAvailable = false; 00049 00050 is_parallel = false; 00051 is_root = true; 00052 rank = 0; 00053 size = 1; 00054 #ifdef MOAB_HAVE_MPI 00055 int flagInit; 00056 MPI_Initialized( &flagInit ); 00057 if (flagInit) { 00058 is_parallel = true; 00059 assert(pcomm != NULL); 00060 rank = pcomm->rank(); 00061 size = pcomm->size(); 00062 is_root = (rank == 0); 00063 } 00064 #endif 00065 00066 moab::DebugOutput dbgprint ( std::cout, rank ); 00067 00068 // Compute and store the total number of source and target DoFs corresponding 00069 // to number of rows and columns in the mapping. 00070 00071 // Initialize dimension information from file 00072 std::vector<std::string> dimNames(1); 00073 std::vector<int> dimSizes(1); 00074 dimNames[0] = "num_elem"; 00075 00076 dbgprint.printf ( 0, "Initializing dimensions of map\n" ); 00077 dbgprint.printf ( 0, "Input mesh\n" ); 00078 dimSizes[0] = m_meshInputCov->faces.size(); 00079 this->InitializeSourceDimensions(dimNames, dimSizes); 00080 dbgprint.printf ( 0, "Output mesh\n" ); 00081 dimSizes[0] = m_meshOutput->faces.size(); 00082 this->InitializeTargetDimensions(dimNames, dimSizes); 00083 00084 m_weightMapGlobal = NULL; 00085 00086 // Build a matrix of source and target discretization so that we know how to assign 00087 // the global DoFs in parallel for the mapping weights 00088 // For example, FV->FV: rows X cols = faces_source X faces_target 00089 } 00090 00091 /////////////////////////////////////////////////////////////////////////////// 00092 00093 moab::TempestOfflineMap::~TempestOfflineMap() 00094 { 00095 delete m_weightMapGlobal; 00096 mbCore = NULL; 00097 #ifdef MOAB_HAVE_MPI 00098 pcomm = NULL; 00099 #endif 00100 m_meshInput = NULL; 00101 m_meshOutput = NULL; 00102 m_meshOverlap = NULL; 00103 } 00104 00105 /////////////////////////////////////////////////////////////////////////////// 00106 00107 static void ParseVariableList ( 00108 const std::string & strVariables, 00109 std::vector< std::string > & vecVariableStrings 00110 ) 00111 { 00112 unsigned iVarBegin = 0; 00113 unsigned iVarCurrent = 0; 00114 00115 // Parse variable name 00116 for ( ;; ) 00117 { 00118 if ( ( iVarCurrent >= strVariables.length() ) || 00119 ( strVariables[iVarCurrent] == ',' ) || 00120 ( strVariables[iVarCurrent] == ' ' ) 00121 ) 00122 { 00123 if ( iVarCurrent == iVarBegin ) 00124 { 00125 if ( iVarCurrent >= strVariables.length() ) 00126 { 00127 break; 00128 } 00129 continue; 00130 } 00131 00132 vecVariableStrings.push_back ( 00133 strVariables.substr ( iVarBegin, iVarCurrent - iVarBegin ) ); 00134 00135 iVarBegin = iVarCurrent + 1; 00136 } 00137 00138 iVarCurrent++; 00139 } 00140 } 00141 00142 /////////////////////////////////////////////////////////////////////////////// 00143 00144 moab::ErrorCode moab::TempestOfflineMap::SetDofMapTags(const std::string srcDofTagName, const std::string tgtDofTagName) 00145 { 00146 moab::ErrorCode rval; 00147 00148 rval = mbCore->tag_get_handle ( srcDofTagName.c_str(), m_nDofsPEl_Src*m_nDofsPEl_Src, MB_TYPE_INTEGER, 00149 this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR(rval); 00150 rval = mbCore->tag_get_handle ( tgtDofTagName.c_str(), m_nDofsPEl_Dest*m_nDofsPEl_Dest, MB_TYPE_INTEGER, 00151 this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR(rval); 00152 00153 return moab::MB_SUCCESS; 00154 } 00155 00156 /////////////////////////////////////////////////////////////////////////////// 00157 00158 moab::ErrorCode moab::TempestOfflineMap::SetDofMapAssociation(DiscretizationType srcType, bool isSrcContinuous, DataMatrix3D<int>* srcdataGLLNodes, DataMatrix3D<int>* srcdataGLLNodesSrc, 00159 DiscretizationType destType, bool isTgtContinuous, DataMatrix3D<int>* tgtdataGLLNodes) 00160 { 00161 moab::ErrorCode rval; 00162 std::vector<bool> dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap; 00163 std::vector<int> src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs; 00164 00165 // We are assuming that these are element based tags that are sized: np * np 00166 m_srcDiscType = srcType; 00167 m_destDiscType = destType; 00168 00169 bool vprint = is_root && false; 00170 00171 #ifdef VVERBOSE 00172 { 00173 src_soln_gdofs.resize(m_remapper->m_covering_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, -1); 00174 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR(rval); 00175 locsrc_soln_gdofs.resize(m_remapper->m_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src); 00176 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR(rval); 00177 tgt_soln_gdofs.resize(m_remapper->m_target_entities.size()*m_nDofsPEl_Dest*m_nDofsPEl_Dest); 00178 rval = mbCore->tag_get_data ( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR(rval); 00179 00180 if (is_root) 00181 { 00182 { 00183 std::ofstream output_file ( "sourcecov-gids-0.txt" ); 00184 output_file << "I, GDOF\n"; 00185 for (unsigned i=0; i < src_soln_gdofs.size(); ++i) 00186 output_file << i << ", " << src_soln_gdofs[i] << "\n"; 00187 00188 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n"; 00189 m_nTotDofs_SrcCov=0; 00190 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00191 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ ) 00192 { 00193 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00194 { 00195 for ( int q = 0; q < m_nDofsPEl_Src; q++) 00196 { 00197 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1; 00198 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00199 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) { 00200 m_nTotDofs_SrcCov++; 00201 dgll_cgll_covcol_ldofmap[ldof] = true; 00202 } 00203 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << m_nTotDofs_SrcCov << "\n"; 00204 } 00205 } 00206 } 00207 output_file.flush(); // required here 00208 output_file.close(); 00209 dgll_cgll_covcol_ldofmap.clear(); 00210 } 00211 00212 { 00213 std::ofstream output_file ( "source-gids-0.txt" ); 00214 output_file << "I, GDOF\n"; 00215 for (unsigned i=0; i < locsrc_soln_gdofs.size(); ++i) 00216 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n"; 00217 00218 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n"; 00219 m_nTotDofs_Src=0; 00220 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00221 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ ) 00222 { 00223 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00224 { 00225 for ( int q = 0; q < m_nDofsPEl_Src; q++) 00226 { 00227 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1; 00228 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00229 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) { 00230 m_nTotDofs_Src++; 00231 dgll_cgll_col_ldofmap[ldof] = true; 00232 } 00233 output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << m_nTotDofs_Src << "\n"; 00234 } 00235 } 00236 } 00237 output_file.flush(); // required here 00238 output_file.close(); 00239 dgll_cgll_col_ldofmap.clear(); 00240 } 00241 00242 { 00243 std::ofstream output_file ( "target-gids-0.txt" ); 00244 output_file << "I, GDOF\n"; 00245 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) 00246 output_file << i << ", " << tgt_soln_gdofs[i] << "\n"; 00247 00248 output_file << "ELEMID, IDOF, GDOF, NDOF\n"; 00249 m_nTotDofs_Dest=0; 00250 00251 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) { 00252 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", " << m_nTotDofs_Dest << "\n"; 00253 m_nTotDofs_Dest++; 00254 } 00255 00256 output_file.flush(); // required here 00257 output_file.close(); 00258 } 00259 } 00260 else 00261 { 00262 { 00263 std::ofstream output_file ( "sourcecov-gids-1.txt" ); 00264 output_file << "I, GDOF\n"; 00265 for (unsigned i=0; i < src_soln_gdofs.size(); ++i) 00266 output_file << i << ", " << src_soln_gdofs[i] << "\n"; 00267 00268 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n"; 00269 m_nTotDofs_SrcCov=0; 00270 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00271 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ ) 00272 { 00273 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00274 { 00275 for ( int q = 0; q < m_nDofsPEl_Src; q++) 00276 { 00277 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1; 00278 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00279 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) { 00280 m_nTotDofs_SrcCov++; 00281 dgll_cgll_covcol_ldofmap[ldof] = true; 00282 } 00283 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << m_nTotDofs_SrcCov << "\n"; 00284 } 00285 } 00286 } 00287 output_file.flush(); // required here 00288 output_file.close(); 00289 dgll_cgll_covcol_ldofmap.clear(); 00290 } 00291 00292 { 00293 std::ofstream output_file ( "source-gids-1.txt" ); 00294 output_file << "I, GDOF\n"; 00295 for (unsigned i=0; i < locsrc_soln_gdofs.size(); ++i) 00296 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n"; 00297 00298 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n"; 00299 m_nTotDofs_Src=0; 00300 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00301 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ ) 00302 { 00303 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00304 { 00305 for ( int q = 0; q < m_nDofsPEl_Src; q++) 00306 { 00307 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1; 00308 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00309 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) { 00310 m_nTotDofs_Src++; 00311 dgll_cgll_col_ldofmap[ldof] = true; 00312 } 00313 output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << m_nTotDofs_Src << "\n"; 00314 } 00315 } 00316 } 00317 output_file.flush(); // required here 00318 output_file.close(); 00319 dgll_cgll_col_ldofmap.clear(); 00320 } 00321 00322 { 00323 std::ofstream output_file ( "target-gids-1.txt" ); 00324 output_file << "I, GDOF\n"; 00325 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) 00326 output_file << i << ", " << tgt_soln_gdofs[i] << "\n"; 00327 00328 output_file << "ELEMID, IDOF, GDOF, NDOF\n"; 00329 m_nTotDofs_Dest=0; 00330 00331 for (unsigned i=0; i < tgt_soln_gdofs.size(); ++i) { 00332 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", " << m_nTotDofs_Dest << "\n"; 00333 m_nTotDofs_Dest++; 00334 } 00335 00336 output_file.flush(); // required here 00337 output_file.close(); 00338 } 00339 } 00340 } 00341 #endif 00342 00343 // Now compute the mapping and store it for the covering mesh 00344 col_dofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00345 col_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00346 col_gdofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00347 src_soln_gdofs.resize(m_remapper->m_covering_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX); 00348 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR(rval); 00349 m_nTotDofs_SrcCov = 0; 00350 if (srcdataGLLNodes == NULL || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */ 00351 for (unsigned i=0; i < col_dofmap.size(); ++i) { 00352 col_dofmap[i] = i; 00353 col_ldofmap[i] = i; 00354 col_gdofmap[i] = src_soln_gdofs[i]; 00355 if (vprint) std::cout << "Col: " << i << ", " << src_soln_gdofs[i] << "\n"; 00356 m_nTotDofs_SrcCov++; 00357 } 00358 } 00359 else { 00360 if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00361 // Put these remap coefficients into the SparseMatrix map 00362 for ( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ ) 00363 { 00364 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00365 { 00366 for ( int q = 0; q < m_nDofsPEl_Src; q++) 00367 { 00368 const int ldof = (*srcdataGLLNodes)[p][q][j] - 1; 00369 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00370 if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) { 00371 m_nTotDofs_SrcCov++; 00372 dgll_cgll_covcol_ldofmap[ldof] = true; 00373 } 00374 if ( !isSrcContinuous ) m_nTotDofs_SrcCov++; 00375 col_dofmap[ idof ] = ldof; 00376 col_ldofmap[ ldof ] = idof; 00377 assert(src_soln_gdofs[idof] > 0); 00378 col_gdofmap[ idof ] = src_soln_gdofs[idof] - 1; 00379 if (vprint) std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << col_gdofmap[idof] << ", " << m_nTotDofs_SrcCov << "\n"; 00380 } 00381 } 00382 } 00383 } 00384 00385 srccol_dofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00386 srccol_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00387 srccol_gdofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX); 00388 locsrc_soln_gdofs.resize(m_remapper->m_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX); 00389 rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR(rval); 00390 00391 // Now compute the mapping and store it for the original source mesh 00392 m_nTotDofs_Src = 0; 00393 if (srcdataGLLNodesSrc == NULL || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */ 00394 for (unsigned i=0; i < srccol_dofmap.size(); ++i) { 00395 srccol_dofmap[i] = i; 00396 srccol_ldofmap[i] = i; 00397 srccol_gdofmap[i] = locsrc_soln_gdofs[i]; 00398 m_nTotDofs_Src++; 00399 } 00400 } 00401 else { 00402 if (isSrcContinuous) dgll_cgll_col_ldofmap.resize(m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false); 00403 // Put these remap coefficients into the SparseMatrix map 00404 for ( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ ) 00405 { 00406 for ( int p = 0; p < m_nDofsPEl_Src; p++ ) 00407 { 00408 for ( int q = 0; q < m_nDofsPEl_Src; q++ ) 00409 { 00410 const int ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1; 00411 const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q; 00412 if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) { 00413 m_nTotDofs_Src++; 00414 dgll_cgll_col_ldofmap[ldof] = true; 00415 } 00416 if ( !isSrcContinuous ) m_nTotDofs_Src++; 00417 srccol_dofmap[ idof ] = ldof; 00418 srccol_ldofmap[ ldof ] = idof; 00419 srccol_gdofmap[ idof ] = locsrc_soln_gdofs[idof] - 1; 00420 } 00421 } 00422 } 00423 } 00424 00425 row_dofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX); 00426 row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX); 00427 row_gdofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX); 00428 tgt_soln_gdofs.resize(m_remapper->m_target_entities.size()*m_nDofsPEl_Dest*m_nDofsPEl_Dest, INT_MAX); 00429 rval = mbCore->tag_get_data ( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR(rval); 00430 00431 // Now compute the mapping and store it for the target mesh 00432 m_nTotDofs_Dest = 0; 00433 if (tgtdataGLLNodes == NULL || destType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */ 00434 for (unsigned i=0; i < row_dofmap.size(); ++i) { 00435 row_dofmap[i] = i; 00436 row_ldofmap[i] = i; 00437 row_gdofmap[i] = tgt_soln_gdofs[i]; 00438 // if (vprint) std::cout << "Row: " << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << row_dofmap[i] << ", " << row_ldofmap[i] << ", " << tgt_soln_gdofs[i] << ", " << row_gdofmap[i] << "\n"; 00439 m_nTotDofs_Dest++; 00440 } 00441 } 00442 else { 00443 if (isTgtContinuous) dgll_cgll_row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, false); 00444 // Put these remap coefficients into the SparseMatrix map 00445 for ( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ ) 00446 { 00447 for ( int p = 0; p < m_nDofsPEl_Dest; p++ ) 00448 { 00449 for ( int q = 0; q < m_nDofsPEl_Dest; q++ ) 00450 { 00451 const int ldof = (*tgtdataGLLNodes)[p][q][j] - 1; 00452 const int idof = j * m_nDofsPEl_Dest * m_nDofsPEl_Dest + p * m_nDofsPEl_Dest + q; 00453 if ( isTgtContinuous && !dgll_cgll_row_ldofmap[ldof] ) { 00454 m_nTotDofs_Dest++; 00455 dgll_cgll_row_ldofmap[ldof] = true; 00456 } 00457 if ( !isTgtContinuous ) m_nTotDofs_Dest++; 00458 row_dofmap[ idof ] = ldof; 00459 row_ldofmap[ ldof ] = idof; 00460 row_gdofmap[ idof ] = tgt_soln_gdofs[idof] - 1; 00461 if (vprint) std::cout << "Row: " << idof << ", " << ldof << ", " << tgt_soln_gdofs[idof] - 1 << "\n"; 00462 } 00463 } 00464 } 00465 } 00466 00467 // Let us also allocate the local representation of the sparse matrix 00468 #ifdef MOAB_HAVE_EIGEN 00469 // if (vprint) 00470 { 00471 std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_dofmap.size() << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_dofmap.size() << "\n"; 00472 // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n"; 00473 } 00474 #endif 00475 00476 return moab::MB_SUCCESS; 00477 } 00478 00479 /////////////////////////////////////////////////////////////////////////////// 00480 00481 moab::ErrorCode moab::TempestOfflineMap::GenerateOfflineMap ( std::string strInputType, std::string strOutputType, 00482 const int nPin, const int nPout, 00483 bool fBubble, int fMonotoneTypeID, 00484 bool fVolumetric, bool fNoConservation, bool fNoCheck, 00485 const std::string srcDofTagName, const std::string tgtDofTagName, 00486 const std::string strVariables, 00487 const std::string strInputData, const std::string strOutputData, 00488 const std::string strNColName, const bool fOutputDouble, 00489 const std::string strPreserveVariables, const bool fPreserveAll, const double dFillValueOverride, 00490 const bool fInputConcave, const bool fOutputConcave ) 00491 { 00492 NcError error ( NcError::silent_nonfatal ); 00493 00494 moab::DebugOutput dbgprint ( std::cout, ( rank ) ); 00495 moab::ErrorCode rval; 00496 00497 try 00498 { 00499 // Check command line parameters (data arguments) 00500 if ( ( strInputData != "" ) && ( strOutputData == "" ) ) 00501 { 00502 _EXCEPTIONT ( "--in_data specified without --out_data" ); 00503 } 00504 if ( ( strInputData == "" ) && ( strOutputData != "" ) ) 00505 { 00506 _EXCEPTIONT ( "--out_data specified without --in_data" ); 00507 } 00508 00509 // Check command line parameters (data type arguments) 00510 STLStringHelper::ToLower ( strInputType ); 00511 STLStringHelper::ToLower ( strOutputType ); 00512 00513 DiscretizationType eInputType; 00514 DiscretizationType eOutputType; 00515 00516 if ( strInputType == "fv" ) 00517 { 00518 eInputType = DiscretizationType_FV; 00519 } 00520 else if ( strInputType == "cgll" ) 00521 { 00522 eInputType = DiscretizationType_CGLL; 00523 } 00524 else if ( strInputType == "dgll" ) 00525 { 00526 eInputType = DiscretizationType_DGLL; 00527 } 00528 else 00529 { 00530 _EXCEPTION1 ( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", 00531 strInputType.c_str() ); 00532 } 00533 00534 if ( strOutputType == "fv" ) 00535 { 00536 eOutputType = DiscretizationType_FV; 00537 } 00538 else if ( strOutputType == "cgll" ) 00539 { 00540 eOutputType = DiscretizationType_CGLL; 00541 } 00542 else if ( strOutputType == "dgll" ) 00543 { 00544 eOutputType = DiscretizationType_DGLL; 00545 } 00546 else 00547 { 00548 _EXCEPTION1 ( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", 00549 strOutputType.c_str() ); 00550 } 00551 00552 // Monotonicity flags 00553 int nMonotoneType = fMonotoneTypeID; 00554 00555 // Parse variable list 00556 std::vector< std::string > vecVariableStrings; 00557 ParseVariableList ( strVariables, vecVariableStrings ); 00558 00559 // Parse preserve variable list 00560 std::vector< std::string > vecPreserveVariableStrings; 00561 ParseVariableList ( strPreserveVariables, vecPreserveVariableStrings ); 00562 00563 if ( fPreserveAll && ( vecPreserveVariableStrings.size() != 0 ) ) 00564 { 00565 _EXCEPTIONT ( "--preserveall and --preserve cannot both be specified" ); 00566 } 00567 00568 m_nDofsPEl_Src = nPin; 00569 m_nDofsPEl_Dest = nPout; 00570 00571 rval = SetDofMapTags(srcDofTagName, tgtDofTagName); 00572 00573 // Calculate Face areas 00574 if ( is_root ) dbgprint.printf ( 0, "Calculating input mesh Face areas\n" ); 00575 double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas(fInputConcave); 00576 double dTotalAreaInput = dTotalAreaInput_loc; 00577 #ifdef MOAB_HAVE_MPI 00578 if (pcomm) MPI_Allreduce ( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00579 #endif 00580 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput ); 00581 00582 // Input mesh areas 00583 m_meshInputCov->CalculateFaceAreas(fInputConcave); 00584 if ( eInputType == DiscretizationType_FV ) 00585 { 00586 this->SetSourceAreas ( m_meshInputCov->vecFaceArea ); 00587 } 00588 00589 // Calculate Face areas 00590 if ( is_root ) dbgprint.printf ( 0, "Calculating output mesh Face areas\n" ); 00591 double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas(fOutputConcave); 00592 double dTotalAreaOutput = dTotalAreaOutput_loc; 00593 #ifdef MOAB_HAVE_MPI 00594 if (pcomm) MPI_Allreduce ( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00595 #endif 00596 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput ); 00597 00598 // Output mesh areas 00599 if ( eOutputType == DiscretizationType_FV ) 00600 { 00601 this->SetTargetAreas ( m_meshOutput->vecFaceArea ); 00602 } 00603 00604 // Verify that overlap mesh is in the correct order 00605 int ixSourceFaceMax = ( -1 ); 00606 int ixTargetFaceMax = ( -1 ); 00607 00608 if ( m_meshOverlap->vecSourceFaceIx.size() != 00609 m_meshOverlap->vecTargetFaceIx.size() 00610 ) 00611 { 00612 _EXCEPTIONT ( "Invalid overlap mesh:\n" 00613 " Possible mesh file corruption?" ); 00614 } 00615 00616 for ( unsigned i = 0; i < m_meshOverlap->vecSourceFaceIx.size(); i++ ) 00617 { 00618 if ( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax ) 00619 { 00620 ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1; 00621 } 00622 if ( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax ) 00623 { 00624 ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1; 00625 } 00626 } 00627 00628 /* 00629 // Check for forward correspondence in overlap mesh 00630 if ( // m_meshInputCov->faces.size() - ixSourceFaceMax == 0 //&& 00631 ( m_meshOutput->faces.size() - ixTargetFaceMax == 0 ) 00632 ) 00633 { 00634 if ( is_root ) dbgprint.printf ( 0, "Overlap mesh forward correspondence found\n" ); 00635 } 00636 else if ( 00637 // m_meshOutput->faces.size() - ixSourceFaceMax == 0 //&& 00638 ( m_meshInputCov->faces.size() - ixTargetFaceMax == 0 ) 00639 ) 00640 { // Check for reverse correspondence in overlap mesh 00641 if ( is_root ) dbgprint.printf ( 0, "Overlap mesh reverse correspondence found (reversing)\n" ); 00642 00643 // Reorder overlap mesh 00644 m_meshOverlap->ExchangeFirstAndSecondMesh(); 00645 } 00646 else 00647 { // No correspondence found 00648 _EXCEPTION4 ( "Invalid overlap mesh:\n" 00649 " No correspondence found with input and output meshes (%i,%i) vs (%i,%i)", 00650 m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax, ixTargetFaceMax ); 00651 } 00652 */ 00653 00654 // Calculate Face areas 00655 if ( is_root ) dbgprint.printf ( 0, "Calculating overlap mesh Face areas\n" ); 00656 double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas(false); 00657 double dTotalAreaOverlap = dTotalAreaOverlap_loc; 00658 #ifdef MOAB_HAVE_MPI 00659 if (pcomm) MPI_Allreduce ( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00660 #endif 00661 if ( is_root ) dbgprint.printf ( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap ); 00662 00663 // Partial cover 00664 if ( fabs ( dTotalAreaOverlap - dTotalAreaInput ) > 1.0e-10 ) 00665 { 00666 if ( !fNoCheck ) 00667 { 00668 if ( is_root ) dbgprint.printf ( 0, "WARNING: Significant mismatch between overlap mesh area " 00669 "and input mesh area.\n Automatically enabling --nocheck\n" ); 00670 fNoCheck = true; 00671 } 00672 } 00673 00674 /* 00675 // Recalculate input mesh area from overlap mesh 00676 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) { 00677 dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n"); 00678 dbgprint.printf(0, "Recalculating source mesh areas\n"); 00679 dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap); 00680 dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput); 00681 } 00682 */ 00683 // Finite volume input / Finite volume output 00684 if ( ( eInputType == DiscretizationType_FV ) && 00685 ( eOutputType == DiscretizationType_FV ) 00686 ) 00687 { 00688 // Generate reverse node array and edge map 00689 m_meshInputCov->ConstructReverseNodeArray(); 00690 m_meshInputCov->ConstructEdgeMap(); 00691 00692 // Initialize coordinates for map 00693 this->InitializeSourceCoordinatesFromMeshFV ( *m_meshInputCov ); 00694 this->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput ); 00695 00696 // Finite volume input / Finite element output 00697 rval = this->SetDofMapAssociation(eInputType, false, NULL, NULL, eOutputType, false, NULL);MB_CHK_ERR(rval); 00698 00699 // Construct OfflineMap 00700 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" ); 00701 LinearRemapFVtoFV_Tempest_MOAB ( nPin ); 00702 } 00703 else if ( eInputType == DiscretizationType_FV ) 00704 { 00705 DataMatrix3D<double> dataGLLJacobian; 00706 00707 if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data\n" ); 00708 double dNumericalArea_loc = 00709 GenerateMetaData ( 00710 *m_meshOutput, 00711 nPout, 00712 fBubble, 00713 dataGLLNodesDest, 00714 dataGLLJacobian ); 00715 00716 double dNumericalArea = dNumericalArea_loc; 00717 #ifdef MOAB_HAVE_MPI 00718 if (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00719 #endif 00720 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea ); 00721 00722 // Initialize coordinates for map 00723 this->InitializeSourceCoordinatesFromMeshFV ( *m_meshInputCov ); 00724 this->InitializeTargetCoordinatesFromMeshFE ( 00725 *m_meshOutput, nPout, dataGLLNodesDest ); 00726 00727 // Generate the continuous Jacobian 00728 bool fContinuous = ( eOutputType == DiscretizationType_CGLL ); 00729 00730 if ( eOutputType == DiscretizationType_CGLL ) 00731 { 00732 GenerateUniqueJacobian ( 00733 dataGLLNodesDest, 00734 dataGLLJacobian, 00735 this->GetTargetAreas() ); 00736 } 00737 else 00738 { 00739 GenerateDiscontinuousJacobian ( 00740 dataGLLJacobian, 00741 this->GetTargetAreas() ); 00742 } 00743 00744 // Generate reverse node array and edge map 00745 m_meshInputCov->ConstructReverseNodeArray(); 00746 m_meshInputCov->ConstructEdgeMap(); 00747 00748 // Finite volume input / Finite element output 00749 rval = this->SetDofMapAssociation(eInputType, false, NULL, NULL, 00750 eOutputType, (eOutputType == DiscretizationType_CGLL), &dataGLLNodesDest);MB_CHK_ERR(rval); 00751 00752 // Generate remap weights 00753 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" ); 00754 00755 if ( fVolumetric ) 00756 { 00757 LinearRemapFVtoGLL_Volumetric_MOAB ( 00758 dataGLLNodesDest, 00759 dataGLLJacobian, 00760 this->GetTargetAreas(), 00761 nPin, 00762 nMonotoneType, 00763 fContinuous, 00764 fNoConservation ); 00765 } 00766 else 00767 { 00768 LinearRemapFVtoGLL_MOAB ( 00769 dataGLLNodesDest, 00770 dataGLLJacobian, 00771 this->GetTargetAreas(), 00772 nPin, 00773 nMonotoneType, 00774 fContinuous, 00775 fNoConservation ); 00776 } 00777 } 00778 else if ( 00779 ( eInputType != DiscretizationType_FV ) && 00780 ( eOutputType == DiscretizationType_FV ) 00781 ) 00782 { 00783 DataMatrix3D<double> dataGLLJacobianSrc, dataGLLJacobian; 00784 00785 if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data\n" ); 00786 // double dNumericalAreaCov_loc = 00787 GenerateMetaData ( 00788 *m_meshInputCov, 00789 nPin, 00790 fBubble, 00791 dataGLLNodesSrcCov, 00792 dataGLLJacobian ); 00793 00794 double dNumericalArea_loc = 00795 GenerateMetaData ( 00796 *m_meshInput, 00797 nPin, 00798 fBubble, 00799 dataGLLNodesSrc, 00800 dataGLLJacobianSrc ); 00801 00802 // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area: %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc ); 00803 // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc); 00804 00805 double dNumericalArea = dNumericalArea_loc; 00806 #ifdef MOAB_HAVE_MPI 00807 if (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00808 #endif 00809 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea ); 00810 00811 if ( fabs ( dNumericalArea - dTotalAreaInput ) > 1.0e-12 ) 00812 { 00813 dbgprint.printf ( 0, "WARNING: Significant mismatch between input mesh " 00814 "numerical area and geometric area\n" ); 00815 } 00816 00817 if ( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() ) 00818 { 00819 _EXCEPTIONT ( "Number of element does not match between metadata and " 00820 "input mesh" ); 00821 } 00822 00823 // Initialize coordinates for map 00824 this->InitializeSourceCoordinatesFromMeshFE ( 00825 *m_meshInputCov, nPin, dataGLLNodesSrcCov ); 00826 this->InitializeSourceCoordinatesFromMeshFE ( 00827 *m_meshInput, nPin, dataGLLNodesSrc ); 00828 this->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput ); 00829 00830 // Generate the continuous Jacobian for input mesh 00831 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL ); 00832 00833 if ( eInputType == DiscretizationType_CGLL ) 00834 { 00835 GenerateUniqueJacobian ( 00836 dataGLLNodesSrcCov, 00837 dataGLLJacobian, 00838 this->GetSourceAreas() ); 00839 } 00840 else 00841 { 00842 GenerateDiscontinuousJacobian ( 00843 dataGLLJacobian, 00844 this->GetSourceAreas() ); 00845 } 00846 00847 // Finite element input / Finite volume output 00848 rval = this->SetDofMapAssociation(eInputType, (eInputType == DiscretizationType_CGLL), &dataGLLNodesSrcCov, &dataGLLNodesSrc, 00849 eOutputType, false, NULL);MB_CHK_ERR(rval); 00850 00851 // Generate offline map 00852 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" ); 00853 00854 if ( fVolumetric ) 00855 { 00856 _EXCEPTIONT ( "Unimplemented: Volumetric currently unavailable for" 00857 "GLL input mesh" ); 00858 } 00859 00860 LinearRemapSE4_Tempest_MOAB ( 00861 dataGLLNodesSrcCov, 00862 dataGLLJacobian, 00863 nMonotoneType, 00864 fContinuousIn, 00865 fNoConservation 00866 ); 00867 00868 } 00869 else if ( 00870 ( eInputType != DiscretizationType_FV ) && 00871 ( eOutputType != DiscretizationType_FV ) 00872 ) 00873 { 00874 DataMatrix3D<double> dataGLLJacobianIn, dataGLLJacobianSrc; 00875 DataMatrix3D<double> dataGLLJacobianOut; 00876 00877 // Input metadata 00878 if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data" ); 00879 double dNumericalAreaIn_loc = 00880 GenerateMetaData ( 00881 *m_meshInputCov, 00882 nPin, 00883 fBubble, 00884 dataGLLNodesSrcCov, 00885 dataGLLJacobianIn ); 00886 00887 double dNumericalAreaSrc_loc = 00888 GenerateMetaData ( 00889 *m_meshInput, 00890 nPin, 00891 fBubble, 00892 dataGLLNodesSrc, 00893 dataGLLJacobianSrc ); 00894 00895 assert(dNumericalAreaIn_loc >= dNumericalAreaSrc_loc); 00896 00897 double dNumericalAreaIn = dNumericalAreaSrc_loc; 00898 #ifdef MOAB_HAVE_MPI 00899 if (pcomm) MPI_Allreduce ( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00900 #endif 00901 if ( is_root ) dbgprint.printf ( 0, "Input Mesh Numerical Area: %1.15e", dNumericalAreaIn ); 00902 00903 if ( fabs ( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 ) 00904 { 00905 dbgprint.printf ( 0, "WARNING: Significant mismatch between input mesh " 00906 "numerical area and geometric area" ); 00907 } 00908 00909 // Output metadata 00910 if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data" ); 00911 double dNumericalAreaOut_loc = 00912 GenerateMetaData ( 00913 *m_meshOutput, 00914 nPout, 00915 fBubble, 00916 dataGLLNodesDest, 00917 dataGLLJacobianOut ); 00918 00919 double dNumericalAreaOut = dNumericalAreaOut_loc; 00920 #ifdef MOAB_HAVE_MPI 00921 if (pcomm) MPI_Allreduce ( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() ); 00922 #endif 00923 if ( is_root ) dbgprint.printf ( 0, "Output Mesh Numerical Area: %1.15e", dNumericalAreaOut ); 00924 00925 if ( fabs ( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 ) 00926 { 00927 if ( is_root ) dbgprint.printf ( 0, "WARNING: Significant mismatch between output mesh " 00928 "numerical area and geometric area" ); 00929 } 00930 00931 // Initialize coordinates for map 00932 this->InitializeSourceCoordinatesFromMeshFE ( 00933 *m_meshInputCov, nPin, dataGLLNodesSrcCov ); 00934 this->InitializeSourceCoordinatesFromMeshFE ( 00935 *m_meshInput, nPin, dataGLLNodesSrc ); 00936 this->InitializeTargetCoordinatesFromMeshFE ( 00937 *m_meshOutput, nPout, dataGLLNodesDest ); 00938 00939 // Generate the continuous Jacobian for input mesh 00940 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL ); 00941 00942 if ( eInputType == DiscretizationType_CGLL ) 00943 { 00944 GenerateUniqueJacobian ( 00945 dataGLLNodesSrcCov, 00946 dataGLLJacobianIn, 00947 this->GetSourceAreas() ); 00948 } 00949 else 00950 { 00951 GenerateDiscontinuousJacobian ( 00952 dataGLLJacobianIn, 00953 this->GetSourceAreas() ); 00954 } 00955 00956 // Generate the continuous Jacobian for output mesh 00957 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL ); 00958 00959 if ( eOutputType == DiscretizationType_CGLL ) 00960 { 00961 GenerateUniqueJacobian ( 00962 dataGLLNodesDest, 00963 dataGLLJacobianOut, 00964 this->GetTargetAreas() ); 00965 } 00966 else 00967 { 00968 GenerateDiscontinuousJacobian ( 00969 dataGLLJacobianOut, 00970 this->GetTargetAreas() ); 00971 } 00972 00973 // Input Finite Element to Output Finite Element 00974 rval = this->SetDofMapAssociation(eInputType, (eInputType == DiscretizationType_CGLL), &dataGLLNodesSrcCov, &dataGLLNodesSrc, 00975 eOutputType, (eOutputType == DiscretizationType_CGLL), &dataGLLNodesDest);MB_CHK_ERR(rval); 00976 00977 // Generate offline map 00978 if ( is_root ) dbgprint.printf ( 0, "Calculating offline map" ); 00979 00980 LinearRemapGLLtoGLL2_MOAB ( 00981 dataGLLNodesSrcCov, 00982 dataGLLJacobianIn, 00983 dataGLLNodesDest, 00984 dataGLLJacobianOut, 00985 this->GetTargetAreas(), 00986 nPin, 00987 nPout, 00988 nMonotoneType, 00989 fContinuousIn, 00990 fContinuousOut, 00991 fNoConservation 00992 ); 00993 00994 } 00995 else 00996 { 00997 _EXCEPTIONT ( "Not implemented" ); 00998 } 00999 01000 #ifdef MOAB_HAVE_EIGEN 01001 CopyTempestSparseMat_Eigen(); 01002 #endif 01003 01004 // Verify consistency, conservation and monotonicity 01005 // gather weights to root process to perform consistency/conservation checks 01006 if ( !fNoCheck && false) 01007 { 01008 if ( !m_globalMapAvailable && size > 1 ) { 01009 // gather weights to root process to perform consistency/conservation checks 01010 rval = this->GatherAllToRoot();MB_CHK_ERR(rval); 01011 } 01012 01013 if ( is_root ) dbgprint.printf ( 0, "Verifying map" ); 01014 this->IsConsistent ( 1.0e-8 ); 01015 if ( !fNoConservation ) this->IsConservative ( 1.0e-8 ); 01016 01017 if ( nMonotoneType != 0 ) 01018 { 01019 this->IsMonotone ( 1.0e-12 ); 01020 } 01021 } 01022 01023 // Apply Offline Map to data 01024 if ( strInputData != "" ) 01025 { 01026 if ( is_root ) dbgprint.printf ( 0, "Applying offline map to data\n" ); 01027 01028 this->SetFillValueOverride ( static_cast<float> ( dFillValueOverride ) ); 01029 this->Apply ( 01030 strInputData, 01031 strOutputData, 01032 vecVariableStrings, 01033 strNColName, 01034 fOutputDouble, 01035 false ); 01036 } 01037 01038 // Copy variables from input file to output file 01039 if ( ( strInputData != "" ) && ( strOutputData != "" ) ) 01040 { 01041 if ( fPreserveAll ) 01042 { 01043 if ( is_root ) dbgprint.printf ( 0, "Preserving variables" ); 01044 this->PreserveAllVariables ( strInputData, strOutputData ); 01045 01046 } 01047 else if ( vecPreserveVariableStrings.size() != 0 ) 01048 { 01049 if ( is_root ) dbgprint.printf ( 0, "Preserving variables" ); 01050 this->PreserveVariables ( 01051 strInputData, 01052 strOutputData, 01053 vecPreserveVariableStrings ); 01054 } 01055 } 01056 01057 } 01058 catch ( Exception & e ) 01059 { 01060 dbgprint.printf ( 0, "%s", e.ToString().c_str() ); 01061 return ( moab::MB_FAILURE ); 01062 01063 } 01064 catch ( ... ) 01065 { 01066 return ( moab::MB_FAILURE ); 01067 } 01068 return moab::MB_SUCCESS; 01069 } 01070 01071 /////////////////////////////////////////////////////////////////////////////// 01072 #ifdef MOAB_HAVE_EIGEN 01073 01074 int moab::TempestOfflineMap::GetSourceGlobalNDofs() 01075 { 01076 return m_weightMatrix.cols(); // return the global number of rows from the weight matrix 01077 } 01078 01079 // /////////////////////////////////////////////////////////////////////////////// 01080 01081 int moab::TempestOfflineMap::GetDestinationGlobalNDofs() 01082 { 01083 return m_weightMatrix.rows(); // return the global number of columns from the weight matrix 01084 } 01085 01086 /////////////////////////////////////////////////////////////////////////////// 01087 01088 int moab::TempestOfflineMap::GetSourceLocalNDofs() 01089 { 01090 return m_weightMatrix.cols(); // return the local number of rows from the weight matrix 01091 } 01092 01093 /////////////////////////////////////////////////////////////////////////////// 01094 01095 int moab::TempestOfflineMap::GetDestinationLocalNDofs() 01096 { 01097 return m_weightMatrix.rows(); // return the local number of columns from the weight matrix 01098 } 01099 01100 /////////////////////////////////////////////////////////////////////////////// 01101 01102 int moab::TempestOfflineMap::GetSourceNDofsPerElement() 01103 { 01104 return m_nDofsPEl_Src; 01105 } 01106 01107 /////////////////////////////////////////////////////////////////////////////// 01108 01109 int moab::TempestOfflineMap::GetDestinationNDofsPerElement() 01110 { 01111 return m_nDofsPEl_Dest; 01112 } 01113 01114 /////////////////////////////////////////////////////////////////////////////// 01115 01116 void moab::TempestOfflineMap::InitVectors() 01117 { 01118 assert(m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0); 01119 m_rowVector.resize( m_weightMatrix.rows() ); 01120 m_colVector.resize( m_weightMatrix.cols() ); 01121 } 01122 01123 01124 moab::TempestOfflineMap::WeightMatrix& moab::TempestOfflineMap::GetWeightMatrix() 01125 { 01126 assert(m_weightMatrix.rows() != 0 && m_weightMatrix.cols() != 0); 01127 return m_weightMatrix; 01128 } 01129 01130 moab::TempestOfflineMap::WeightRowVector& moab::TempestOfflineMap::GetRowVector() 01131 { 01132 assert(m_rowVector.size() != 0); 01133 return m_rowVector; 01134 } 01135 01136 moab::TempestOfflineMap::WeightColVector& moab::TempestOfflineMap::GetColVector() 01137 { 01138 assert(m_colVector.size() != 0); 01139 return m_colVector; 01140 } 01141 01142 #endif 01143 01144 /////////////////////////////////////////////////////////////////////////////// 01145 01146 bool moab::TempestOfflineMap::IsConsistent ( 01147 double dTolerance 01148 ) 01149 { 01150 // Get map entries 01151 DataVector<int> dataRows; 01152 DataVector<int> dataCols; 01153 DataVector<double> dataEntries; 01154 01155 // Calculate row sums 01156 DataVector<double> dRowSums; 01157 if ( size > 1 ) 01158 { 01159 if ( rank ) return true; 01160 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix(); 01161 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries ); 01162 dRowSums.Initialize ( m_mapRemapGlobal.GetRows() ); 01163 } 01164 else 01165 { 01166 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries ); 01167 dRowSums.Initialize ( m_mapRemap.GetRows() ); 01168 } 01169 01170 for ( unsigned i = 0; i < dataRows.GetRows(); i++ ) 01171 { 01172 dRowSums[dataRows[i]] += dataEntries[i]; 01173 } 01174 01175 // Verify all row sums are equal to 1 01176 bool fConsistent = true; 01177 for ( unsigned i = 0; i < dRowSums.GetRows(); i++ ) 01178 { 01179 if ( fabs ( dRowSums[i] - 1.0 ) > dTolerance ) 01180 { 01181 fConsistent = false; 01182 Announce ( "TempestOfflineMap is not consistent in row %i (%1.15e)", 01183 i, dRowSums[i] ); 01184 } 01185 } 01186 01187 return fConsistent; 01188 } 01189 01190 /////////////////////////////////////////////////////////////////////////////// 01191 01192 bool moab::TempestOfflineMap::IsConservative ( 01193 double dTolerance 01194 ) 01195 { 01196 // Get map entries 01197 DataVector<int> dataRows; 01198 DataVector<int> dataCols; 01199 DataVector<double> dataEntries; 01200 const DataVector<double>& dTargetAreas = this->GetGlobalTargetAreas(); 01201 const DataVector<double>& dSourceAreas = this->GetGlobalSourceAreas(); 01202 01203 // Calculate column sums 01204 DataVector<double> dColumnSums; 01205 if ( size > 1 ) 01206 { 01207 if ( rank ) return true; 01208 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix(); 01209 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries ); 01210 dColumnSums.Initialize ( m_mapRemapGlobal.GetColumns() ); 01211 } 01212 else 01213 { 01214 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries ); 01215 dColumnSums.Initialize ( m_mapRemap.GetColumns() ); 01216 } 01217 01218 for ( unsigned i = 0; i < dataRows.GetRows(); i++ ) 01219 { 01220 dColumnSums[dataCols[i]] += 01221 dataEntries[i] * dTargetAreas[dataRows[i]]; 01222 } 01223 01224 // Verify all column sums equal the input Jacobian 01225 bool fConservative = true; 01226 for ( unsigned i = 0; i < dColumnSums.GetRows(); i++ ) 01227 { 01228 if ( fabs ( dColumnSums[i] - dSourceAreas[i] ) > dTolerance ) 01229 { 01230 fConservative = false; 01231 Announce ( "TempestOfflineMap is not conservative in column " 01232 "%i (%1.15e / %1.15e)", 01233 i, dColumnSums[i], dSourceAreas[i] ); 01234 } 01235 } 01236 01237 return fConservative; 01238 } 01239 01240 /////////////////////////////////////////////////////////////////////////////// 01241 01242 bool moab::TempestOfflineMap::IsMonotone ( 01243 double dTolerance 01244 ) 01245 { 01246 // Get map entries 01247 DataVector<int> dataRows; 01248 DataVector<int> dataCols; 01249 DataVector<double> dataEntries; 01250 01251 if ( size > 1 ) 01252 { 01253 if ( rank ) return true; 01254 SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix(); 01255 m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries ); 01256 } 01257 else 01258 m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries ); 01259 01260 // Verify all entries are in the range [0,1] 01261 bool fMonotone = true; 01262 for ( unsigned i = 0; i < dataRows.GetRows(); i++ ) 01263 { 01264 if ( ( dataEntries[i] < -dTolerance ) || 01265 ( dataEntries[i] > 1.0 + dTolerance ) 01266 ) 01267 { 01268 fMonotone = false; 01269 01270 Announce ( "TempestOfflineMap is not monotone in entry (%i): %1.15e", 01271 i, dataEntries[i] ); 01272 } 01273 } 01274 01275 return fMonotone; 01276 } 01277 01278 01279 /////////////////////////////////////////////////////////////////////////////// 01280 #ifdef MOAB_HAVE_MPI 01281 moab::ErrorCode moab::TempestOfflineMap::GatherAllToRoot() // Collective 01282 { 01283 Mesh globalMesh; 01284 int ierr, rootProc = 0, nprocs = size; 01285 moab::ErrorCode rval; 01286 01287 // Write SparseMatrix entries 01288 DataVector<int> vecRow; 01289 DataVector<int> vecCol; 01290 DataVector<double> vecS; 01291 01292 moab::DebugOutput dbgprint ( std::cout, ( rank ) ); 01293 01294 m_mapRemap.GetEntries ( vecRow, vecCol, vecS ); 01295 const DataVector<double>& dOrigSourceAreas = m_meshInput->vecFaceArea; 01296 const DataVector<double>& dSourceAreas = m_meshInputCov->vecFaceArea; 01297 const DataVector<double>& dTargetAreas = m_meshOutput->vecFaceArea; 01298 01299 // Translate the index in Row and Col to global_id and dump it out 01300 01301 // Communicate the necessary data 01302 const int NDATA = 7; 01303 int gnnz = 0, gsrc = 0, gsrccov = 0, gtar = 0, gsrcdofs = 0, gsrccovdofs = 0, gtgtdofs = 0; 01304 std::vector<int> rootSizesData, rowcolsv; 01305 DataVector<int> rows, cols, srcelmindx, tgtelmindx; 01306 { 01307 // First, accumulate the sizes of rows and columns of the matrix 01308 if ( !rank ) rootSizesData.resize ( size * NDATA ); 01309 01310 int sendarray[NDATA]; 01311 sendarray[0] = vecS.GetRows(); 01312 sendarray[1] = dOrigSourceAreas.GetRows(); 01313 sendarray[2] = dTargetAreas.GetRows(); 01314 sendarray[3] = dSourceAreas.GetRows(); 01315 sendarray[4] = m_nTotDofs_Src; 01316 sendarray[5] = m_nTotDofs_Dest; 01317 sendarray[6] = m_nTotDofs_SrcCov; 01318 01319 ierr = MPI_Gather ( sendarray, NDATA, MPI_INTEGER, rootSizesData.data(), NDATA, MPI_INTEGER, rootProc, pcomm->comm() ); 01320 if ( ierr != MPI_SUCCESS ) return moab::MB_FAILURE; 01321 01322 if ( !rank ) 01323 { 01324 for ( int i = 0; i < nprocs; ++i ) 01325 { 01326 unsigned offset = NDATA * i; 01327 gnnz += rootSizesData[offset]; 01328 gsrc += rootSizesData[offset + 1]; 01329 gtar += rootSizesData[offset + 2]; 01330 gsrccov += rootSizesData[offset + 3]; 01331 gsrcdofs += rootSizesData[offset + 4]; 01332 gtgtdofs += rootSizesData[offset + 5]; 01333 gsrccovdofs += rootSizesData[offset + 6]; 01334 } 01335 rowcolsv.resize ( 2 * gnnz + gsrc + gtar ); 01336 rows.Initialize ( gnnz ); // we are assuming rows = cols 01337 cols.Initialize ( gnnz ); // we are assuming rows = cols 01338 01339 // Let us allocate our global offline map object 01340 m_weightMapGlobal = new OfflineMap(); 01341 std::vector<std::string> dimNames(1); 01342 std::vector<int> dimSizes(1); 01343 dimNames[0] = "num_elem"; 01344 dimSizes[0] = gsrc; 01345 m_weightMapGlobal->InitializeSourceDimensions(dimNames, dimSizes); 01346 dimSizes[0] = gtar; 01347 m_weightMapGlobal->InitializeTargetDimensions(dimNames, dimSizes); 01348 01349 /* 01350 VSM: do we need to gather the entire mesh ?!?! 01351 The following initialization can't work correctly without the mesh on the root process 01352 */ 01353 if (m_srcDiscType == DiscretizationType_FV) /* unsure if we actually care about the type for this */ 01354 m_weightMapGlobal->InitializeSourceCoordinatesFromMeshFV ( *m_meshInput ); 01355 else 01356 m_weightMapGlobal->InitializeSourceCoordinatesFromMeshFE ( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc ); 01357 01358 if (m_destDiscType == DiscretizationType_FV) /* unsure if we actually care about the type for this */ 01359 m_weightMapGlobal->InitializeTargetCoordinatesFromMeshFV ( *m_meshOutput ); 01360 else 01361 m_weightMapGlobal->InitializeTargetCoordinatesFromMeshFE ( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest ); 01362 01363 DataVector<double>& m_areasSrcGlobal = m_weightMapGlobal->GetSourceAreas(); 01364 m_areasSrcGlobal.Initialize ( gsrcdofs ); srcelmindx.Initialize ( gsrc ); 01365 DataVector<double>& m_areasTgtGlobal = m_weightMapGlobal->GetTargetAreas(); 01366 m_areasTgtGlobal.Initialize ( gtgtdofs ); tgtelmindx.Initialize ( gtar ); 01367 01368 #ifdef VERBOSE 01369 dbgprint.printf ( 0, "Received global dimensions: %d, %d\n", vecRow.GetRows(), rows.GetRows() ); 01370 dbgprint.printf ( 0, "Global: n(source) = %d, n(srccov) = %d, and n(target) = %d\n", gsrc, gsrccov, gtar ); 01371 dbgprint.printf ( 0, "Operator size = %d X %d and NNZ = %d\n", gsrcdofs, gtgtdofs, gnnz ); 01372 #endif 01373 } 01374 } 01375 01376 #ifdef VERBOSE 01377 { 01378 std::stringstream sstr; 01379 sstr << "rowscols_" << rank << ".txt"; 01380 std::ofstream output_file ( sstr.str().c_str() ); 01381 output_file << "VALUES\n"; 01382 for ( unsigned ip = 0; ip < vecS.GetRows(); ++ip ) 01383 { 01384 output_file << ip << " (" << row_gdofmap[row_ldofmap[vecRow[ip]]] << ", " << col_gdofmap[col_ldofmap[vecCol[ip]]] << ") = " << vecS[ip] << "\n"; 01385 01386 } 01387 output_file.flush(); // required here 01388 output_file.close(); 01389 } 01390 #endif 01391 01392 { 01393 // Next, accumulate the row and column values for the matrices (in local indexing) 01394 // Let's gather the data in the following order: 01395 // 1) row indices, 01396 // 2) column indices, 01397 // 3) GID for source elements, 01398 // 4) GID for target elements 01399 // 01400 const int nR = 2 * vecS.GetRows() + dOrigSourceAreas.GetRows() + dTargetAreas.GetRows(); 01401 std::vector<int> sendarray ( nR ); 01402 for ( unsigned ix = 0; ix < vecS.GetRows(); ++ix ) 01403 { 01404 sendarray[ix] = row_gdofmap[row_ldofmap[vecRow[ix]]]; 01405 } 01406 for ( unsigned ix = 0, offset = vecS.GetRows(); ix < vecS.GetRows(); ++ix ) 01407 { 01408 sendarray[offset + ix] = col_gdofmap[col_ldofmap[vecCol[ix]]]; 01409 } 01410 01411 { 01412 moab::Tag gidtag; 01413 rval = mbCore->tag_get_handle ( "GLOBAL_ID", gidtag ); MB_CHK_ERR ( rval ); 01414 rval = mbCore->tag_get_data ( gidtag, m_remapper->m_source_entities, &sendarray[2 * vecS.GetRows()] ); MB_CHK_ERR ( rval ); 01415 rval = mbCore->tag_get_data ( gidtag, m_remapper->m_target_entities, &sendarray[2 * vecS.GetRows() + dOrigSourceAreas.GetRows()] ); MB_CHK_ERR ( rval ); 01416 } 01417 01418 std::vector<int> displs, rcount; 01419 if ( !rank ) 01420 { 01421 displs.resize ( size, 0 ); 01422 rcount.resize ( size, 0 ); 01423 int gsum = 0; 01424 for ( int i = 0; i < nprocs; ++i ) 01425 { 01426 displs[i] = gsum; 01427 rcount[i] = 2 * rootSizesData[NDATA * i] + rootSizesData[NDATA * i + 1] + rootSizesData[NDATA * i + 2]; 01428 /* + rootSizesData[NDATA * i + 3] + rootSizesData[NDATA * i + 4] */ 01429 gsum += rcount[i]; 01430 } 01431 assert ( rowcolsv.size() - gsum == 0 ); 01432 } 01433 01434 // Both rows and columns have a size of "rowsize" 01435 ierr = MPI_Gatherv ( &sendarray[0], nR, MPI_INTEGER, &rowcolsv[0], &rcount[0], &displs[0], MPI_INTEGER, rootProc, pcomm->comm() ); 01436 01437 if ( !rank ) 01438 { 01439 #ifdef VERBOSE 01440 std::ofstream output_file ( "rows-cols.txt", std::ios::out ); 01441 output_file << "ROWS\n"; 01442 #endif 01443 for ( int ip = 0, offset = 0; ip < nprocs; ++ip ) 01444 { 01445 int istart = displs[ip], iend = istart + rootSizesData[NDATA * ip]; 01446 for ( int i = istart; i < iend; ++i, ++offset ) 01447 { 01448 rows[offset] = rowcolsv[i]; 01449 #ifdef VERBOSE 01450 output_file << offset << " " << rows[offset] << "\n"; 01451 #endif 01452 } 01453 } 01454 #ifdef VERBOSE 01455 output_file << "COLS\n"; 01456 #endif 01457 for ( int ip = 0, offset = 0; ip < nprocs; ++ip ) 01458 { 01459 int istart = displs[ip] + rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip]; 01460 for ( int i = istart; i < iend; ++i, ++offset ) 01461 { 01462 cols[offset] = rowcolsv[i]; 01463 #ifdef VERBOSE 01464 output_file << offset << " " << cols[offset] << "\n"; 01465 #endif 01466 } 01467 } 01468 #ifdef VERBOSE 01469 output_file.flush(); // required here 01470 output_file.close(); 01471 #endif 01472 for ( int ip = 0, offset = 0; ip < nprocs; ++ip ) 01473 { 01474 int istart = displs[ip] + 2 * rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip + 1]; 01475 for ( int i = istart; i < iend; ++i, ++offset ) 01476 { 01477 srcelmindx[offset] = rowcolsv[i]; 01478 } 01479 } 01480 for ( int ip = 0, offset = 0; ip < nprocs; ++ip ) 01481 { 01482 int istart = displs[ip] + 2 * rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1], iend = istart + rootSizesData[NDATA * ip + 2]; 01483 for ( int i = istart; i < iend; ++i, ++offset ) 01484 { 01485 tgtelmindx[offset] = rowcolsv[i]; 01486 } 01487 } 01488 } /* (!rank) */ 01489 rowcolsv.clear(); 01490 } 01491 01492 int nSc = m_dSourceVertexLon.GetColumns(); 01493 int nTc = m_dTargetVertexLon.GetColumns(); 01494 { 01495 // Next, accumulate the row and column values for the matrices (in local indexing) 01496 // Let's gather the data in the following order: 01497 // 1) matrix weights, 01498 // 2) source element areas, 01499 // 3) target element areas, 01500 // 4) source element: m_dSourceCenterLon 01501 // 5) source element: m_dSourceCenterLat 01502 // 6) target element: m_dTargetCenterLon 01503 // 7) target element: m_dTargetCenterLat 01504 // 01505 const int nR = vecS.GetRows() + dOrigSourceAreas.GetRows() + dTargetAreas.GetRows() + 2*m_nTotDofs_Src*(1+nSc) + 2*m_nTotDofs_Dest*(1+nTc); 01506 01507 std::vector<double> sendarray ( nR ); 01508 int locoffset = 0; 01509 std::copy ( ( double* ) vecS, ( double* ) vecS + vecS.GetRows(), sendarray.begin() + locoffset); locoffset += vecS.GetRows(); 01510 std::copy ( ( const double* ) dOrigSourceAreas, ( const double* ) dOrigSourceAreas + dOrigSourceAreas.GetRows(), sendarray.begin() + locoffset ); locoffset += dOrigSourceAreas.GetRows(); 01511 std::copy ( ( const double* ) dTargetAreas, ( const double* ) dTargetAreas + dTargetAreas.GetRows(), sendarray.begin() + locoffset ); locoffset += dTargetAreas.GetRows(); 01512 01513 std::cout << "[" << rank << "] " << m_nTotDofs_Src << ", m_dSourceCenterLon.size() = " << m_dSourceCenterLon.GetRows() << " and " << m_nTotDofs_Dest << ", m_dTargetCenterLon.size() = " << m_dTargetCenterLon.GetRows() << "\n"; 01514 01515 std::copy ( ( const double* ) m_dSourceCenterLon, ( const double* ) m_dSourceCenterLon + m_dSourceCenterLon.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dSourceCenterLon.GetRows(); 01516 std::copy ( ( const double* ) m_dSourceCenterLat, ( const double* ) m_dSourceCenterLat + m_dSourceCenterLat.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dSourceCenterLat.GetRows(); 01517 01518 double** dSCLon = m_dSourceVertexLon; 01519 std::copy ( &dSCLon[0][0], &dSCLon[0][0] + m_dSourceVertexLon.GetRows()*nSc, sendarray.begin() + locoffset ); locoffset += m_dSourceVertexLon.GetRows()*nSc; 01520 double** dSCLat = m_dSourceVertexLat; 01521 std::copy ( &dSCLat[0][0], &dSCLat[0][0] + m_dSourceVertexLat.GetRows()*nSc, sendarray.begin() + locoffset ); locoffset += m_dSourceVertexLat.GetRows()*nSc; 01522 01523 std::copy ( ( const double* ) m_dTargetCenterLon, ( const double* ) m_dTargetCenterLon + m_dTargetCenterLon.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dTargetCenterLon.GetRows(); 01524 std::copy ( ( const double* ) m_dTargetCenterLat, ( const double* ) m_dTargetCenterLat + m_dTargetCenterLat.GetRows(), sendarray.begin() + locoffset ); locoffset += m_dTargetCenterLat.GetRows(); 01525 01526 double** dTCLon = m_dTargetVertexLon; 01527 std::copy ( &dTCLon[0][0], &dTCLon[0][0] + m_dTargetVertexLon.GetRows()*nTc, sendarray.begin() + locoffset ); locoffset += m_dTargetVertexLon.GetRows()*nTc; 01528 double** dTCLat = m_dTargetVertexLat; 01529 std::copy ( &dTCLat[0][0], &dTCLat[0][0] + m_dTargetVertexLat.GetRows()*nTc, sendarray.begin() + locoffset ); locoffset += m_dTargetVertexLat.GetRows()*nTc; 01530 01531 std::vector<int> displs, rcount; 01532 int gsum = 0; 01533 if ( !rank ) 01534 { 01535 displs.resize ( size, 0 ); 01536 rcount.resize ( size, 0 ); 01537 for ( int i = 0; i < nprocs; ++i ) 01538 { 01539 displs[i] = gsum; 01540 rcount[i] = rootSizesData[NDATA * i] + rootSizesData[NDATA * i + 1] + rootSizesData[NDATA * i + 2] + 01541 2 * rootSizesData[NDATA * i + 4] * (1 + nSc) + 01542 2 * rootSizesData[NDATA * i + 5] * (1 + nTc); 01543 gsum += rcount[i]; 01544 } 01545 } 01546 01547 std::vector<double> rowcolsvals ( (!rank ? gnnz + gsrc + gtar + 2*gsrcdofs*(1+nSc) + 2*gtgtdofs*(1+nTc) : 0 ) ); 01548 // Both rows and columns have a size of "rowsize" 01549 ierr = MPI_Gatherv ( sendarray.data(), sendarray.size(), MPI_DOUBLE, rowcolsvals.data(), &rcount[0], &displs[0], MPI_DOUBLE, rootProc, pcomm->comm() ); 01550 01551 if ( !rank ) 01552 { 01553 SparseMatrix<double>& spmat = m_weightMapGlobal->GetSparseMatrix(); 01554 { 01555 #ifdef VERBOSE 01556 std::ofstream output_file ( "rows-cols.txt", std::ios::app ); 01557 output_file << "VALUES\n"; 01558 #endif 01559 for ( int ip = 0, offset = 0; ip < nprocs; ++ip ) 01560 { 01561 for ( int i = displs[ip]; i < displs[ip] + rootSizesData[NDATA * ip]; ++i, ++offset ) 01562 { 01563 // globvalues[offset] = rowcolsvals[i]; 01564 spmat(rows[offset], cols[offset]) += rowcolsvals[i]; 01565 #ifdef VERBOSE 01566 output_file << offset << " (" << rows[offset] << ", " << cols[offset] << ") = " << rowcolsvals[i] << "\n"; 01567 #endif 01568 } 01569 } 01570 #ifdef VERBOSE 01571 output_file.flush(); // required here 01572 output_file.close(); 01573 #endif 01574 01575 } 01576 // m_mapRemapGlobal.SetEntries(rows, cols, globvalues); 01577 // m_weightMapGlobal->GetSparseMatrix().AddEntries ( rows, cols, globvalues ); 01578 01579 { 01580 DataVector<double>& m_areasSrcGlobal = m_weightMapGlobal->GetSourceAreas(); 01581 DataVector<double>& m_areasTgtGlobal = m_weightMapGlobal->GetTargetAreas(); 01582 // Store the global source and target elements areas 01583 #ifdef VERBOSE 01584 std::ofstream output_file ( "source-target-areas.txt" ); 01585 output_file << "Source areas (nelems = " << gsrc << ")\n"; 01586 #endif 01587 int offset = 0; 01588 for ( int ip = 0; ip < nprocs; ++ip ) 01589 { 01590 int istart = displs[ip] + rootSizesData[NDATA * ip], iend = istart + rootSizesData[NDATA * ip + 1]; 01591 for ( int i = istart; i < iend; ++i, ++offset ) 01592 { 01593 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs); 01594 m_areasSrcGlobal[srcelmindx[offset]] = rowcolsvals[i]; 01595 } 01596 } 01597 01598 #ifdef VERBOSE 01599 for ( unsigned i = 0; i < m_areasSrcGlobal.GetRows(); ++i ) 01600 output_file << srcelmindx[i] << " " << m_areasSrcGlobal[i] << "\n"; 01601 #endif 01602 01603 offset = 0; 01604 for ( int ip = 0; ip < nprocs; ++ip ) 01605 { 01606 int istart = displs[ip] + rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1], iend = istart + rootSizesData[NDATA * ip + 2]; 01607 for ( int i = istart; i < iend; ++i, ++offset ) 01608 { 01609 // assert(offset < gtgtdofs && tgtelmindx[offset] < gtgtdofs); 01610 m_areasTgtGlobal[tgtelmindx[offset]] = rowcolsvals[i]; 01611 } 01612 } 01613 #ifdef VERBOSE 01614 output_file << "Target areas (nelems = " << gtar << ")\n"; 01615 for ( unsigned i = 0; i < m_areasTgtGlobal.GetRows(); ++i ) 01616 output_file << tgtelmindx[i] << " " << m_areasTgtGlobal[i] << "\n"; 01617 #endif 01618 #ifdef VERBOSE 01619 output_file.flush(); // required here 01620 output_file.close(); 01621 #endif 01622 } 01623 01624 DataVector<double>& dSourceCenterLon = m_weightMapGlobal->GetSourceCenterLon(); 01625 DataVector<double>& dSourceCenterLat = m_weightMapGlobal->GetSourceCenterLat(); 01626 dSourceCenterLon.Initialize(gsrcdofs); 01627 dSourceCenterLat.Initialize(gsrcdofs); 01628 { 01629 int offset = 0; 01630 for ( int ip = 0; ip < nprocs; ++ip ) 01631 { 01632 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2]; 01633 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4]; 01634 for ( int i = istart; i < iend; ++i, ++offset ) 01635 { 01636 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs); 01637 // dSourceCenterLon[srcelmindx[offset]] = rowcolsvals[i]; 01638 dSourceCenterLon[offset] = rowcolsvals[i]; 01639 } 01640 } 01641 01642 offset = 0; 01643 for ( int ip = 0; ip < nprocs; ++ip ) 01644 { 01645 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + rootSizesData[NDATA * ip + 4]; 01646 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4]; 01647 for ( int i = istart; i < iend; ++i, ++offset ) 01648 { 01649 // dSourceCenterLat[srcelmindx[offset]] = rowcolsvals[i]; 01650 dSourceCenterLat[offset] = rowcolsvals[i]; 01651 } 01652 } 01653 } 01654 01655 01656 DataMatrix<double>& dSourceVertexLon = m_weightMapGlobal->GetSourceVertexLon(); 01657 DataMatrix<double>& dSourceVertexLat = m_weightMapGlobal->GetSourceVertexLat(); 01658 dSourceVertexLon.Initialize(gsrcdofs, nSc); 01659 dSourceVertexLat.Initialize(gsrcdofs, nSc); 01660 { 01661 int ioffset = 0; 01662 for ( int ip = 0; ip < nprocs; ++ip ) 01663 { 01664 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4]; 01665 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4] * nSc; 01666 for ( int i = istart; i < iend; ++i, ++ioffset ) 01667 { 01668 for ( int joffset = 0; joffset < nSc; ++joffset) 01669 dSourceVertexLon[ioffset][joffset] = rowcolsvals[i+joffset]; 01670 i += nSc; 01671 } 01672 } 01673 01674 ioffset = 0; 01675 for ( int ip = 0; ip < nprocs; ++ip ) 01676 { 01677 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] + rootSizesData[NDATA * ip + 4] * nSc; 01678 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 4] * nSc; 01679 for ( int i = istart; i < iend; ++i, ++ioffset ) 01680 { 01681 for ( int joffset = 0; joffset < nSc; ++joffset) 01682 dSourceVertexLat[ioffset][joffset] = rowcolsvals[i+joffset]; 01683 i += nSc; 01684 } 01685 } 01686 } 01687 01688 DataVector<double>& dTargetCenterLon = m_weightMapGlobal->GetTargetCenterLon(); 01689 DataVector<double>& dTargetCenterLat = m_weightMapGlobal->GetTargetCenterLat(); 01690 dTargetCenterLon.Initialize(gtgtdofs); 01691 dTargetCenterLat.Initialize(gtgtdofs); 01692 { 01693 int offset = 0; 01694 for ( int ip = 0; ip < nprocs; ++ip ) 01695 { 01696 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc); 01697 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]; 01698 for ( int i = istart; i < iend; ++i, ++offset ) 01699 { 01700 // assert(offset < gsrcdofs && srcelmindx[offset] <= gsrcdofs); 01701 // dTargetCenterLon[tgtelmindx[offset]] = rowcolsvals[i]; 01702 dTargetCenterLon[offset] = rowcolsvals[i]; 01703 } 01704 } 01705 01706 offset = 0; 01707 for ( int ip = 0; ip < nprocs; ++ip ) 01708 { 01709 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + rootSizesData[NDATA * ip + 5]; 01710 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]; 01711 for ( int i = istart; i < iend; ++i, ++offset ) 01712 { 01713 // dTargetCenterLat[tgtelmindx[offset]] = rowcolsvals[i]; 01714 dTargetCenterLat[offset] = rowcolsvals[i]; 01715 } 01716 } 01717 } 01718 01719 DataMatrix<double>& dTargetVertexLon = m_weightMapGlobal->GetTargetVertexLon(); 01720 DataMatrix<double>& dTargetVertexLat = m_weightMapGlobal->GetTargetVertexLat(); 01721 dTargetVertexLon.Initialize(gtgtdofs, nTc); 01722 dTargetVertexLat.Initialize(gtgtdofs, nTc); 01723 { 01724 int ioffset = 0; 01725 int nc = dTargetVertexLon.GetColumns(); 01726 for ( int ip = 0; ip < nprocs; ++ip ) 01727 { 01728 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + 2 * rootSizesData[NDATA * ip + 5]; 01729 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]*nTc; 01730 for ( int i = istart; i < iend; ++i, ++ioffset ) 01731 { 01732 for ( int joffset = 0; joffset < nc; ++joffset) 01733 dTargetVertexLon[ioffset][joffset] = rowcolsvals[i+joffset]; 01734 i += nc; 01735 } 01736 } 01737 01738 ioffset = 0; 01739 for ( int ip = 0; ip < nprocs; ++ip ) 01740 { 01741 int ibase = rootSizesData[NDATA * ip] + rootSizesData[NDATA * ip + 1] + rootSizesData[NDATA * ip + 2] + 2 * rootSizesData[NDATA * ip + 4] * (1+nSc) + 2 * rootSizesData[NDATA * ip + 5] + rootSizesData[NDATA * ip + 5]*nTc; 01742 int istart = displs[ip] + ibase, iend = istart + rootSizesData[NDATA * ip + 5]*nTc; 01743 for ( int i = istart; i < iend; ++i, ++ioffset ) 01744 { 01745 for ( int joffset = 0; joffset < nc; ++joffset) 01746 dTargetVertexLat[ioffset][joffset] = rowcolsvals[i+joffset]; 01747 i += nc; 01748 } 01749 } 01750 } 01751 } 01752 01753 } 01754 rootSizesData.clear(); 01755 01756 // Update all processes that we have a global view of the map available 01757 // on the root process 01758 m_globalMapAvailable = true; 01759 01760 #ifdef VERBOSE 01761 if ( !rank ) 01762 { 01763 dbgprint.printf ( 0, "Writing out file outGlobalView.nc\n" ); 01764 m_weightMapGlobal->Write ( "outGlobalView.nc" ); 01765 } 01766 #endif 01767 return moab::MB_SUCCESS; 01768 } 01769 01770 /////////////////////////////////////////////////////////////////////////////// 01771 #endif 01772