MOAB: Mesh Oriented datABase  (version 5.1.1)
TempestOfflineMap.cpp
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines