|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.4.1)
|
An offline map between two Meshes. More...
#include <TempestOnlineMap.hpp>
Collaboration diagram for moab::TempestOnlineMap:Public Types | |
| enum | DiscretizationType { DiscretizationType_FV, DiscretizationType_CGLL, DiscretizationType_DGLL, DiscretizationType_PCLOUD } |
| typedef double(* | sample_function )(double, double) |
Public Member Functions | |
| TempestOnlineMap (moab::TempestRemapper *remapper) | |
| Generate the metadata associated with the offline map. | |
| virtual | ~TempestOnlineMap () |
| Define a virtual destructor. | |
| moab::ErrorCode | GenerateRemappingWeights (std::string strInputType, std::string strOutputType, const GenerateOfflineMapAlgorithmOptions &mapOptions, const std::string &srcDofTagName="GLOBAL_ID", const std::string &tgtDofTagName="GLOBAL_ID") |
| Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix. | |
| moab::ErrorCode | ReadParallelMap (const char *strSource, const std::vector< int > &owned_dof_ids, bool row_major_ownership=true) |
| Generate the metadata associated with the offline map. | |
| moab::ErrorCode | WriteParallelMap (const std::string &strTarget) |
| Write the TempestOnlineMap to a parallel NetCDF file. | |
| virtual int | IsConsistent (double dTolerance) |
| Determine if the map is first-order accurate. | |
| virtual int | IsConservative (double dTolerance) |
| Determine if the map is conservative. | |
| virtual int | IsMonotone (double dTolerance) |
| Determine if the map is monotone. | |
| const DataArray1D< double > & | GetGlobalSourceAreas () const |
| If we computed the reduction, get the vector representing the source areas for all entities in the mesh. | |
| const DataArray1D< double > & | GetGlobalTargetAreas () const |
| If we computed the reduction, get the vector representing the target areas for all entities in the mesh. | |
| moab::ErrorCode | SetDOFmapTags (const std::string srcDofTagName, const std::string tgtDofTagName) |
| Store the tag names associated with global DoF ids for source and target meshes. | |
| moab::ErrorCode | SetDOFmapAssociation (DiscretizationType srcType, bool isSrcContinuous, DataArray3D< int > *srcdataGLLNodes, DataArray3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, bool isDestContinuous, DataArray3D< int > *tgtdataGLLNodes) |
| Compute the association between the solution tag global DoF numbering and the local matrix numbering so that matvec operations can be performed consistently. | |
| int | GetSourceGlobalNDofs () |
| Get the number of total Degrees-Of-Freedom defined on the source mesh. | |
| int | GetDestinationGlobalNDofs () |
| Get the number of total Degrees-Of-Freedom defined on the destination mesh. | |
| int | GetSourceLocalNDofs () |
| Get the number of local Degrees-Of-Freedom defined on the source mesh. | |
| int | GetDestinationLocalNDofs () |
| Get the number of local Degrees-Of-Freedom defined on the destination mesh. | |
| int | GetSourceNDofsPerElement () |
| Get the number of Degrees-Of-Freedom per element on the source mesh. | |
| int | GetDestinationNDofsPerElement () |
| Get the number of Degrees-Of-Freedom per element on the destination mesh. | |
| void | SetSourceNDofsPerElement (int ns) |
| Set the number of Degrees-Of-Freedom per element on the source mesh. | |
| void | SetDestinationNDofsPerElement (int nt) |
| Get the number of Degrees-Of-Freedom per element on the destination mesh. | |
| int | GetRowGlobalDoF (int localID) const |
| Get the global Degrees-Of-Freedom ID on the destination mesh. | |
| int | GetIndexOfRowGlobalDoF (int globalRowDoF) const |
| Get the index of globaRowDoF. | |
| int | GetColGlobalDoF (int localID) const |
| Get the global Degrees-Of-Freedom ID on the source mesh. | |
| int | GetIndexOfColGlobalDoF (int globalColDoF) const |
| Get the index of globaColDoF. | |
| moab::ErrorCode | ApplyWeights (std::vector< double > &srcVals, std::vector< double > &tgtVals, bool transpose=false) |
Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * , or if (transpose) tgtVals = [A(T->S)]^T * . | |
| moab::ErrorCode | ApplyWeights (moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose=false) |
Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * , or if (transpose) tgtVals = [A(T->S)]^T * . | |
| moab::ErrorCode | DefineAnalyticalSolution (moab::Tag &exactSolnTag, const std::string &solnName, Remapper::IntersectionContext ctx, sample_function testFunction, moab::Tag *clonedSolnTag=NULL, std::string cloneSolnName="") |
| Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user. | |
| moab::ErrorCode | ComputeMetrics (Remapper::IntersectionContext ctx, moab::Tag &exactTag, moab::Tag &approxTag, std::map< std::string, double > &metrics, bool verbose=true) |
| Compute the error between a sampled (exact) solution and a projected solution in various error norms. | |
| moab::ErrorCode | fill_row_ids (std::vector< int > &ids_of_interest) |
| moab::ErrorCode | fill_col_ids (std::vector< int > &ids_of_interest) |
| moab::ErrorCode | set_col_dc_dofs (std::vector< int > &values_entities) |
| moab::ErrorCode | set_row_dc_dofs (std::vector< int > &values_entities) |
Private Member Functions | |
| moab::ErrorCode | LinearRemapNN_MOAB (bool use_GID_matching=false, bool strict_check=false) |
| Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh. | |
| void | LinearRemapFVtoFV_Tempest_MOAB (int nOrder) |
| Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh. | |
| void | LinearRemapSE0_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian) |
| Generate the OfflineMap for linear conserative element-average spectral element to element average remapping. | |
| void | LinearRemapSE4_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, int nMonotoneType, bool fContinuousIn, bool fNoConservation) |
| Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping. | |
| void | LinearRemapFVtoGLL_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, const DataArray1D< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite volumes to finite elements. | |
| void | LinearRemapGLLtoGLL2_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite elements to finite elements. | |
| void | LinearRemapGLLtoGLL2_Pointwise_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut) |
| Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation). | |
| moab::ErrorCode | WriteSCRIPMapFile (const std::string &strOutputFile) |
| Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections. | |
| moab::ErrorCode | WriteHDF5MapFile (const std::string &filename) |
| Parallel I/O with NetCDF to write out the SCRIP file from multiple processors. | |
| void | setup_sizes_dimensions () |
Private Attributes | |
| moab::TempestRemapper * | m_remapper |
| The fundamental remapping operator object. | |
| moab::Interface * | m_interface |
| The reference to the moab::Core object that contains source/target and overlap sets. | |
| moab::Tag | m_dofTagSrc |
| The original tag data and local to global DoF mapping to associate matrix values to solution. | |
| moab::Tag | m_dofTagDest |
| std::vector< unsigned > | row_gdofmap |
| std::vector< unsigned > | col_gdofmap |
| std::vector< unsigned > | srccol_gdofmap |
| std::vector< int > | row_dtoc_dofmap |
| std::vector< int > | col_dtoc_dofmap |
| std::vector< int > | srccol_dtoc_dofmap |
| std::map< int, int > | rowMap |
| std::map< int, int > | colMap |
| DataArray3D< int > | dataGLLNodesSrc |
| DataArray3D< int > | dataGLLNodesSrcCov |
| DataArray3D< int > | dataGLLNodesDest |
| DiscretizationType | m_srcDiscType |
| DiscretizationType | m_destDiscType |
| int | m_nTotDofs_Src |
| int | m_nTotDofs_SrcCov |
| int | m_nTotDofs_Dest |
| int | m_nDofsPEl_Src |
| int | m_nDofsPEl_Dest |
| DiscretizationType | m_eInputType |
| DiscretizationType | m_eOutputType |
| bool | m_bConserved |
| int | m_iMonotonicity |
| Mesh * | m_meshInput |
| Mesh * | m_meshInputCov |
| Mesh * | m_meshOutput |
| Mesh * | m_meshOverlap |
| bool | is_parallel |
| bool | is_root |
| int | rank |
| int | size |
| int | root_proc |
An offline map between two Meshes.
Definition at line 48 of file TempestOnlineMap.hpp.
| typedef double( * moab::TempestOnlineMap::sample_function)(double, double) |
Definition at line 347 of file TempestOnlineMap.hpp.
| DiscretizationType_FV | |
| DiscretizationType_CGLL | |
| DiscretizationType_DGLL | |
| DiscretizationType_PCLOUD |
Definition at line 64 of file TempestOnlineMap.hpp.
Generate the metadata associated with the offline map.
Definition at line 55 of file TempestOnlineMap.cpp.
References moab::Remapper::get_interface(), moab::TempestRemapper::GetCoveringMesh(), moab::TempestRemapper::GetMesh(), is_parallel, is_root, m_interface, m_meshInput, m_meshInputCov, m_meshOutput, m_meshOverlap, m_remapper, moab::Remapper::OverlapMesh, rank, root_proc, setup_sizes_dimensions(), size, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
: OfflineMap(), m_remapper( remapper ) { // Get the references for the MOAB core objects m_interface = m_remapper->get_interface(); #ifdef MOAB_HAVE_MPI m_pcomm = m_remapper->get_parallel_communicator(); #endif // Update the references to the meshes m_meshInput = remapper->GetMesh( moab::Remapper::SourceMesh ); m_meshInputCov = remapper->GetCoveringMesh(); m_meshOutput = remapper->GetMesh( moab::Remapper::TargetMesh ); m_meshOverlap = remapper->GetMesh( moab::Remapper::OverlapMesh ); is_parallel = false; is_root = true; rank = 0; root_proc = rank; size = 1; #ifdef MOAB_HAVE_MPI int flagInit; MPI_Initialized( &flagInit ); if( flagInit ) { is_parallel = true; assert( m_pcomm != NULL ); rank = m_pcomm->rank(); size = m_pcomm->size(); is_root = ( rank == 0 ); } #endif // Compute and store the total number of source and target DoFs corresponding // to number of rows and columns in the mapping. // Initialize dimension information from file this->setup_sizes_dimensions(); // Build a matrix of source and target discretization so that we know how to assign // the global DoFs in parallel for the mapping weights // For example, FV->FV: rows X cols = faces_source X faces_target }
| moab::TempestOnlineMap::~TempestOnlineMap | ( | ) | [virtual] |
Define a virtual destructor.
Definition at line 145 of file TempestOnlineMap.cpp.
{
m_interface = NULL;
#ifdef MOAB_HAVE_MPI
m_pcomm = NULL;
#endif
m_meshInput = NULL;
m_meshOutput = NULL;
m_meshOverlap = NULL;
}
| moab::ErrorCode moab::TempestOnlineMap::ApplyWeights | ( | std::vector< double > & | srcVals, |
| std::vector< double > & | tgtVals, | ||
| bool | transpose = false |
||
| ) |
Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * , or if (transpose) tgtVals = [A(T->S)]^T * .
Referenced by main().
| moab::ErrorCode moab::TempestOnlineMap::ApplyWeights | ( | moab::Tag | srcSolutionTag, |
| moab::Tag | tgtSolutionTag, | ||
| bool | transpose = false |
||
| ) |
Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * , or if (transpose) tgtVals = [A(T->S)]^T * .
Definition at line 1700 of file TempestOnlineMap.cpp.
References moab::Remapper::CoveringMesh, ErrorCode, MB_CHK_SET_ERR, MB_SUCCESS, moab::Range::size(), and moab::Remapper::TargetMesh.
{
moab::ErrorCode rval;
std::vector< double > solSTagVals;
std::vector< double > solTTagVals;
moab::Range sents, tents;
if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
{
if( m_remapper->point_cloud_source )
{
moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
solSTagVals.resize( covSrcEnts.size(), -1.0 );
sents = covSrcEnts;
}
else
{
moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
-1.0 );
sents = covSrcEnts;
}
if( m_remapper->point_cloud_target )
{
moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
solTTagVals.resize( tgtEnts.size(), -1.0 );
tents = tgtEnts;
}
else
{
moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
solTTagVals.resize(
tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
tents = tgtEnts;
}
}
else
{
moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
-1.0 );
solTTagVals.resize(
tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
sents = covSrcEnts;
tents = tgtEnts;
}
// The tag data is np*np*n_el_src
rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );
// Compute the application of weights on the suorce solution data and store it in the
// destination solution vector data Optionally, can also perform the transpose application of
// the weight matrix. Set the 3rd argument to true if this is needed
rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );
// The tag data is np*np*n_el_dest
rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics | ( | Remapper::IntersectionContext | ctx, |
| moab::Tag & | exactTag, | ||
| moab::Tag & | approxTag, | ||
| std::map< std::string, double > & | metrics, | ||
| bool | verbose = true |
||
| ) |
Compute the error between a sampled (exact) solution and a projected solution in various error norms.
Definition at line 2128 of file TempestOnlineMap.cpp.
References entities, moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, moab::Range::size(), moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
Referenced by main().
{
moab::ErrorCode rval;
const bool outputEnabled = ( is_root );
int discOrder;
DiscretizationType discMethod;
moab::EntityHandle meshset;
moab::Range entities;
Mesh* trmesh;
switch( ctx )
{
case Remapper::SourceMesh:
meshset = m_remapper->m_covering_source_set;
trmesh = m_remapper->m_covering_source;
entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
: m_remapper->m_covering_source_entities );
discOrder = m_nDofsPEl_Src;
discMethod = m_eInputType;
break;
case Remapper::TargetMesh:
meshset = m_remapper->m_target_set;
trmesh = m_remapper->m_target;
entities =
( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
discOrder = m_nDofsPEl_Dest;
discMethod = m_eOutputType;
break;
default:
if( outputEnabled )
std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
return moab::MB_FAILURE;
}
// Let us create teh solution tag with appropriate information for name, discretization order
// (DoF space)
std::string exactTagName, projTagName;
const int ntotsize = entities.size() * discOrder * discOrder;
int ntotsize_glob = 0;
std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 ); // L1Err, L2Err, LinfErr
for( int i = 0; i < ntotsize; ++i )
{
const double error = fabs( exactSolution[i] - projSolution[i] );
errnorms[0] += error;
errnorms[1] += error * error;
errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] );
}
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
{
MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
}
#else
ntotsize_glob = ntotsize;
globerrnorms = errnorms;
#endif
globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob );
globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob );
metrics.clear();
metrics["L1Error"] = globerrnorms[0];
metrics["L2Error"] = globerrnorms[1];
metrics["LinfError"] = globerrnorms[2];
if( verbose && is_root )
{
std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl;
}
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution | ( | moab::Tag & | exactSolnTag, |
| const std::string & | solnName, | ||
| Remapper::IntersectionContext | ctx, | ||
| sample_function | testFunction, | ||
| moab::Tag * | clonedSolnTag = NULL, |
||
| std::string | cloneSolnName = "" |
||
| ) |
Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user.
Definition at line 1766 of file TempestOnlineMap.cpp.
References entities, ErrorCode, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, moab::Range::size(), moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
Referenced by main().
{
moab::ErrorCode rval;
const bool outputEnabled = ( is_root );
int discOrder;
DiscretizationType discMethod;
moab::EntityHandle meshset;
moab::Range entities;
Mesh* trmesh;
switch( ctx )
{
case Remapper::SourceMesh:
meshset = m_remapper->m_covering_source_set;
trmesh = m_remapper->m_covering_source;
entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
: m_remapper->m_covering_source_entities );
discOrder = m_nDofsPEl_Src;
discMethod = m_eInputType;
break;
case Remapper::TargetMesh:
meshset = m_remapper->m_target_set;
trmesh = m_remapper->m_target;
entities =
( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
discOrder = m_nDofsPEl_Dest;
discMethod = m_eOutputType;
break;
default:
if( outputEnabled )
std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
return moab::MB_FAILURE;
}
// Let us create teh solution tag with appropriate information for name, discretization order
// (DoF space)
rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
if( clonedSolnTag != NULL )
{
if( cloneSolnName.size() == 0 )
{
cloneSolnName = solnName + std::string( "Cloned" );
}
rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
*clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
}
// Triangular quadrature rule
const int TriQuadratureOrder = 10;
if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
TriangularQuadratureRule triquadrule( TriQuadratureOrder );
const int TriQuadraturePoints = triquadrule.GetPoints();
const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
// Output data
DataArray1D< double > dVar;
DataArray1D< double > dVarMB; // re-arranged local MOAB vector
// Nodal geometric area
DataArray1D< double > dNodeArea;
// Calculate element areas
// trmesh->CalculateFaceAreas(fContainsConcaveFaces);
if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
{
/* Get the spectral points and sample the functionals accordingly */
const bool fGLL = true;
const bool fGLLIntegrate = false;
// Generate grid metadata
DataArray3D< int > dataGLLNodes;
DataArray3D< double > dataGLLJacobian;
GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
// Number of elements
int nElements = trmesh->faces.size();
// Verify all elements are quadrilaterals
for( int k = 0; k < nElements; k++ )
{
const Face& face = trmesh->faces[k];
if( face.edges.size() != 4 )
{
_EXCEPTIONT( "Non-quadrilateral face detected; "
"incompatible with --gll" );
}
}
// Number of unique nodes
int iMaxNode = 0;
for( int i = 0; i < discOrder; i++ )
{
for( int j = 0; j < discOrder; j++ )
{
for( int k = 0; k < nElements; k++ )
{
if( dataGLLNodes[i][j][k] > iMaxNode )
{
iMaxNode = dataGLLNodes[i][j][k];
}
}
}
}
// Get Gauss-Lobatto quadrature nodes
DataArray1D< double > dG;
DataArray1D< double > dW;
GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
// Get Gauss quadrature nodes
const int nGaussP = 10;
DataArray1D< double > dGaussG;
DataArray1D< double > dGaussW;
GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
// Allocate data
dVar.Allocate( iMaxNode );
dVarMB.Allocate( discOrder * discOrder * nElements );
dNodeArea.Allocate( iMaxNode );
// Sample data
for( int k = 0; k < nElements; k++ )
{
const Face& face = trmesh->faces[k];
// Sample data at GLL nodes
if( fGLL )
{
for( int i = 0; i < discOrder; i++ )
{
for( int j = 0; j < discOrder; j++ )
{
// Apply local map
Node node;
Node dDx1G;
Node dDx2G;
ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
// Sample data at this point
double dNodeLon = atan2( node.y, node.x );
if( dNodeLon < 0.0 )
{
dNodeLon += 2.0 * M_PI;
}
double dNodeLat = asin( node.z );
double dSample = ( *testFunction )( dNodeLon, dNodeLat );
dVar[dataGLLNodes[j][i][k] - 1] = dSample;
}
}
// High-order Gaussian integration over basis function
}
else
{
DataArray2D< double > dCoeff( discOrder, discOrder );
for( int p = 0; p < nGaussP; p++ )
{
for( int q = 0; q < nGaussP; q++ )
{
// Apply local map
Node node;
Node dDx1G;
Node dDx2G;
ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
// Cross product gives local Jacobian
Node nodeCross = CrossProduct( dDx1G, dDx2G );
double dJacobian =
sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
// Find components of quadrature point in basis
// of the first Face
SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
// Sample data at this point
double dNodeLon = atan2( node.y, node.x );
if( dNodeLon < 0.0 )
{
dNodeLon += 2.0 * M_PI;
}
double dNodeLat = asin( node.z );
double dSample = ( *testFunction )( dNodeLon, dNodeLat );
// Integrate
for( int i = 0; i < discOrder; i++ )
{
for( int j = 0; j < discOrder; j++ )
{
double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
}
}
}
}
}
}
// Divide by area
if( fGLLIntegrate )
{
for( size_t i = 0; i < dVar.GetRows(); i++ )
{
dVar[i] /= dNodeArea[i];
}
}
// Let us rearrange the data based on DoF ID specification
if( ctx == Remapper::SourceMesh )
{
for( unsigned j = 0; j < entities.size(); j++ )
for( int p = 0; p < discOrder; p++ )
for( int q = 0; q < discOrder; q++ )
{
const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
}
}
else
{
for( unsigned j = 0; j < entities.size(); j++ )
for( int p = 0; p < discOrder; p++ )
for( int q = 0; q < discOrder; q++ )
{
const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
}
}
// Set the tag data
rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
}
else
{
// assert( discOrder == 1 );
if( discMethod == DiscretizationType_FV )
{
/* Compute an element-wise integral to store the sampled solution based on Quadrature
* rules */
// Resize the array
dVar.Allocate( trmesh->faces.size() );
std::vector< Node >& nodes = trmesh->nodes;
// Loop through all Faces
for( size_t i = 0; i < trmesh->faces.size(); i++ )
{
const Face& face = trmesh->faces[i];
// Loop through all sub-triangles
for( size_t j = 0; j < face.edges.size() - 2; j++ )
{
const Node& node0 = nodes[face[0]];
const Node& node1 = nodes[face[j + 1]];
const Node& node2 = nodes[face[j + 2]];
// Triangle area
Face faceTri( 3 );
faceTri.SetNode( 0, face[0] );
faceTri.SetNode( 1, face[j + 1] );
faceTri.SetNode( 2, face[j + 2] );
double dTriangleArea = CalculateFaceArea( faceTri, nodes );
// Calculate the element average
double dTotalSample = 0.0;
// Loop through all quadrature points
for( int k = 0; k < TriQuadraturePoints; k++ )
{
Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
TriQuadratureG[k][2] * node2.x,
TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
TriQuadratureG[k][2] * node2.y,
TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
TriQuadratureG[k][2] * node2.z );
double dMagnitude = node.Magnitude();
node.x /= dMagnitude;
node.y /= dMagnitude;
node.z /= dMagnitude;
double dLon = atan2( node.y, node.x );
if( dLon < 0.0 )
{
dLon += 2.0 * M_PI;
}
double dLat = asin( node.z );
double dSample = ( *testFunction )( dLon, dLat );
dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
}
dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
}
}
rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
}
else /* discMethod == DiscretizationType_PCLOUD */
{
/* Get the coordinates of the vertices and sample the functionals accordingly */
std::vector< Node >& nodes = trmesh->nodes;
// Resize the array
dVar.Allocate( nodes.size() );
for( size_t j = 0; j < nodes.size(); j++ )
{
Node& node = nodes[j];
double dMagnitude = node.Magnitude();
node.x /= dMagnitude;
node.y /= dMagnitude;
node.z /= dMagnitude;
double dLon = atan2( node.y, node.x );
if( dLon < 0.0 )
{
dLon += 2.0 * M_PI;
}
double dLat = asin( node.z );
double dSample = ( *testFunction )( dLon, dLat );
dVar[j] = dSample;
}
rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
}
}
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::fill_col_ids | ( | std::vector< int > & | ids_of_interest | ) | [inline] |
Definition at line 382 of file TempestOnlineMap.hpp.
References col_gdofmap, and MB_SUCCESS.
{
ids_of_interest.reserve( col_gdofmap.size() );
// need to add 1
for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
ids_of_interest.push_back( *it + 1 );
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::fill_row_ids | ( | std::vector< int > & | ids_of_interest | ) | [inline] |
Definition at line 372 of file TempestOnlineMap.hpp.
References MB_SUCCESS, and row_gdofmap.
{
ids_of_interest.reserve( row_gdofmap.size() );
// need to add 1
for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
ids_of_interest.push_back( *it + 1 );
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights | ( | std::string | strInputType, |
| std::string | strOutputType, | ||
| const GenerateOfflineMapAlgorithmOptions & | mapOptions, | ||
| const std::string & | srcDofTagName = "GLOBAL_ID", |
||
| const std::string & | tgtDofTagName = "GLOBAL_ID" |
||
| ) |
Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix.
the tag should be created already in the e3sm workflow; if not, create it here
Definition at line 731 of file TempestOnlineMap.cpp.
References dbgprint, moab::error(), ErrorCode, MB_ALREADY_ALLOCATED, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_EXCL, MB_TYPE_DOUBLE, moab::Remapper::OverlapMesh, moab::DebugOutput::printf(), rank, moab::DebugOutput::set_prefix(), and size.
Referenced by main().
{
NcError error( NcError::silent_nonfatal );
moab::DebugOutput dbgprint( std::cout, rank, 0 );
dbgprint.set_prefix( "[TempestOnlineMap]: " );
moab::ErrorCode rval;
const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
try
{
// Check command line parameters (data type arguments)
STLStringHelper::ToLower( strInputType );
STLStringHelper::ToLower( strOutputType );
DiscretizationType eInputType;
DiscretizationType eOutputType;
if( strInputType == "fv" )
{
eInputType = DiscretizationType_FV;
}
else if( strInputType == "cgll" )
{
eInputType = DiscretizationType_CGLL;
}
else if( strInputType == "dgll" )
{
eInputType = DiscretizationType_DGLL;
}
else if( strInputType == "pcloud" )
{
eInputType = DiscretizationType_PCLOUD;
}
else
{
_EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
}
if( strOutputType == "fv" )
{
eOutputType = DiscretizationType_FV;
}
else if( strOutputType == "cgll" )
{
eOutputType = DiscretizationType_CGLL;
}
else if( strOutputType == "dgll" )
{
eOutputType = DiscretizationType_DGLL;
}
else if( strOutputType == "pcloud" )
{
eOutputType = DiscretizationType_PCLOUD;
}
else
{
_EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
}
// set all required input params
m_bConserved = !mapOptions.fNoConservation;
m_eInputType = eInputType;
m_eOutputType = eOutputType;
// Method flags
std::string strMapAlgorithm( "" );
int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
// Make an index of method arguments
std::set< std::string > setMethodStrings;
{
int iLast = 0;
for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
{
if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
{
std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
if( strMethodString.length() > 0 )
{
setMethodStrings.insert( strMethodString );
}
iLast = i + 1;
}
}
}
for( auto it : setMethodStrings )
{
// Piecewise constant monotonicity
if( it == "mono2" )
{
if( nMonotoneType != 0 )
{
_EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
}
if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
}
nMonotoneType = 2;
// Piecewise linear monotonicity
}
else if( it == "mono3" )
{
if( nMonotoneType != 0 )
{
_EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
}
if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
}
nMonotoneType = 3;
// Volumetric remapping from FV to GLL
}
else if( it == "volumetric" )
{
if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
}
strMapAlgorithm = "volumetric";
// Inverse distance mapping
}
else if( it == "invdist" )
{
if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
}
strMapAlgorithm = "invdist";
// Delaunay triangulation mapping
}
else if( it == "delaunay" )
{
if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
}
strMapAlgorithm = "delaunay";
// Bilinear
}
else if( it == "bilin" )
{
if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
}
strMapAlgorithm = "fvbilin";
// Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
}
else if( it == "intbilin" )
{
if( m_eOutputType != DiscretizationType_FV )
{
_EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
}
if( m_eInputType == DiscretizationType_FV )
{
strMapAlgorithm = "fvintbilin";
}
else
{
strMapAlgorithm = "mono3";
}
// Integrated bilinear with generalized Barycentric coordinates
}
else if( it == "intbilingb" )
{
if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
{
_EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
}
strMapAlgorithm = "fvintbilingb";
}
else
{
_EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
}
}
m_nDofsPEl_Src =
( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
: mapOptions.nPin );
m_nDofsPEl_Dest =
( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
: mapOptions.nPout );
rval = SetDOFmapTags( srcDofTagName, tgtDofTagName );MB_CHK_ERR( rval );
/// the tag should be created already in the e3sm workflow; if not, create it here
Tag areaTag;
rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT );
if ( MB_ALREADY_ALLOCATED == rval )
{
if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
}
double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0;
if( !m_bPointCloudSource )
{
// Calculate Input Mesh Face areas
if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
rval = m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea );MB_CHK_ERR( rval );
dTotalAreaInput = dTotalAreaInput_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
// Input mesh areas
m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
// we do not need to set the area on coverage mesh, only on source and target meshes
// rval = m_interface->tag_set_data( areaTag, m_remapper->m_covering_source_entities, m_meshInputCov->vecFaceArea );MB_CHK_ERR( rval );
}
if( !m_bPointCloudTarget )
{
// Calculate Output Mesh Face areas
if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
dTotalAreaOutput = dTotalAreaOutput_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
rval = m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea );MB_CHK_ERR( rval );
}
if( !m_bPointCloud )
{
// Verify that overlap mesh is in the correct order, only if size == 1
if( 1 == size )
{
int ixSourceFaceMax = ( -1 );
int ixTargetFaceMax = ( -1 );
if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() )
{
_EXCEPTIONT( "Invalid overlap mesh:\n"
" Possible mesh file corruption?" );
}
for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ )
{
if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;
if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
}
// Check for forward correspondence in overlap mesh
if( m_meshInputCov->faces.size() - ixSourceFaceMax == 0 )
{
if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" );
}
else if( m_meshOutput->faces.size() - ixSourceFaceMax == 0 )
{ // Check for reverse correspondence in overlap mesh
if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" );
// Reorder overlap mesh
m_meshOverlap->ExchangeFirstAndSecondMesh();
}
else
{ // No correspondence found
_EXCEPTION4( "Invalid overlap mesh:\n No correspondence found with input and output meshes "
"(%i,%i) vs (%i,%i)",
m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax,
ixTargetFaceMax );
}
}
// Calculate Face areas
if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false );
double dTotalAreaOverlap = dTotalAreaOverlap_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
// Correct areas to match the areas calculated in the overlap mesh
// if (fCorrectAreas) // In MOAB-TempestRemap, we will always keep this to be true
{
if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
{
if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
continue; // skip this cell since it is ghosted
// let us recompute the source/target areas based on overlap mesh areas
assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
}
for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
{
if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
{
m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
}
}
for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
{
if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
{
m_meshOutput->vecFaceArea[i] = dTargetArea[i];
}
}
}
// Set source mesh areas in map
if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
{
this->SetSourceAreas( m_meshInputCov->vecFaceArea );
if( m_meshInputCov->vecMask.size() )
{
this->SetSourceMask( m_meshInputCov->vecMask );
}
}
// Set target mesh areas in map
if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
{
this->SetTargetAreas( m_meshOutput->vecFaceArea );
if( m_meshOutput->vecMask.size() )
{
this->SetTargetMask( m_meshOutput->vecMask );
}
}
/*
// Recalculate input mesh area from overlap mesh
if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
dbgprint.printf(0, "Recalculating source mesh areas\n");
dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
}
*/
}
// Finite volume input / Finite volume output
if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
{
// Generate reverse node array and edge map
if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
// Initialize coordinates for map
this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
// Finite volume input / Finite element output
rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, false, NULL );MB_CHK_ERR( rval );
// Construct remap for FV-FV
if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
// Construct OfflineMap
if( strMapAlgorithm == "invdist" )
{
if( is_root ) AnnounceStartBlock( "Calculating map (invdist)" );
LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
}
else if( strMapAlgorithm == "delaunay" )
{
if( is_root ) AnnounceStartBlock( "Calculating map (delaunay)" );
LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
}
else if( strMapAlgorithm == "fvintbilin" )
{
if( is_root ) AnnounceStartBlock( "Calculating map (intbilin)" );
LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
}
else if( strMapAlgorithm == "fvintbilingb" )
{
if( is_root ) AnnounceStartBlock( "Calculating map (intbilingb)" );
LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
}
else if( strMapAlgorithm == "fvbilin" )
{
if( is_root ) AnnounceStartBlock( "Calculating map (bilin)" );
LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
}
else
{
if( is_root ) AnnounceStartBlock( "Calculating map (default)" );
// LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
// ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
}
}
else if( eInputType == DiscretizationType_FV )
{
DataArray3D< double > dataGLLJacobian;
if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
dataGLLNodesDest, dataGLLJacobian );
double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
// Initialize coordinates for map
this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
// Generate the continuous Jacobian
bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
if( eOutputType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
}
// Generate reverse node array and edge map
if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
// Finite volume input / Finite element output
rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType,
( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
// Generate remap weights
if( strMapAlgorithm == "volumetric" )
{
if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
nMonotoneType, fContinuous, mapOptions.fNoConservation );
}
else
{
if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
mapOptions.fNoConservation );
}
}
else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
{
DataArray3D< double > dataGLLJacobian;
if( !m_bPointCloudSource )
{
// Generate reverse node array and edge map
if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
// Initialize coordinates for map
if( eInputType == DiscretizationType_FV )
{
this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
}
else
{
if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
DataArray3D< double > dataGLLJacobianSrc;
GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
dataGLLJacobian );
GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
dataGLLJacobianSrc );
}
}
// else { /* Source is a point cloud dataset */ }
if( !m_bPointCloudTarget )
{
// Generate reverse node array and edge map
if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
// Initialize coordinates for map
if( eOutputType == DiscretizationType_FV )
{
this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
}
else
{
if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
dataGLLJacobian );
}
}
// else { /* Target is a point cloud dataset */ }
// Finite volume input / Finite element output
rval = this->SetDOFmapAssociation(
eInputType, ( eInputType == DiscretizationType_CGLL ),
( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrcCov ),
( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrc ), eOutputType,
( eOutputType == DiscretizationType_CGLL ), ( m_bPointCloudTarget ? NULL : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
// Construct remap
if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
}
else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
{
DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
// double dNumericalAreaCov_loc =
GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
dataGLLJacobian );
double dNumericalArea_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
dataGLLNodesSrc, dataGLLJacobianSrc );
// if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area:
// %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
// assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);
double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root )
{
dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );
if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
{
dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
"numerical area and geometric area\n" );
}
}
if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
{
_EXCEPTIONT( "Number of element does not match between metadata and "
"input mesh" );
}
// Initialize coordinates for map
this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
// Generate the continuous Jacobian for input mesh
bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
if( eInputType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
}
// Finite element input / Finite volume output
rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
&dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, false, NULL );MB_CHK_ERR( rval );
// Generate remap
if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
if( strMapAlgorithm == "volumetric" )
{
_EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
"GLL input mesh" );
}
LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
mapOptions.fNoConservation );
}
else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
{
DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
DataArray3D< double > dataGLLJacobianOut;
// Input metadata
if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
dataGLLJacobianIn );
double dNumericalAreaSrc_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
dataGLLNodesSrc, dataGLLJacobianSrc );
double dNumericalAreaIn = dNumericalAreaSrc_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root )
{
dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn );
if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
{
dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
"numerical area and geometric area\n" );
}
}
// Output metadata
if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
double dNumericalAreaOut_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
dataGLLNodesDest, dataGLLJacobianOut );
double dNumericalAreaOut = dNumericalAreaOut_loc;
#ifdef MOAB_HAVE_MPI
if( m_pcomm )
MPI_Reduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
if( is_root )
{
dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut );
if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
{
if( is_root )
dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh "
"numerical area and geometric area\n" );
}
}
// Initialize coordinates for map
this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
// Generate the continuous Jacobian for input mesh
bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
if( eInputType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
}
// Generate the continuous Jacobian for output mesh
bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
if( eOutputType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
}
// Input Finite Element to Output Finite Element
rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
&dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType,
( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
// Generate remap
if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
}
else
{
_EXCEPTIONT( "Not implemented" );
}
#ifdef MOAB_HAVE_EIGEN3
copy_tempest_sparsemat_to_eigen3();
#endif
#ifdef MOAB_HAVE_MPI
{
// Remove ghosted entities from overlap set
moab::Range ghostedEnts;
rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
}
#endif
// Verify consistency, conservation and monotonicity, globally
if( !mapOptions.fNoCheck )
{
if( is_root ) dbgprint.printf( 0, "Verifying map" );
this->IsConsistent( 1.0e-8 );
if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
if( nMonotoneType != 0 )
{
this->IsMonotone( 1.0e-12 );
}
}
}
catch( Exception& e )
{
dbgprint.printf( 0, "%s", e.ToString().c_str() );
return ( moab::MB_FAILURE );
}
catch( ... )
{
return ( moab::MB_FAILURE );
}
return moab::MB_SUCCESS;
}
| int moab::TempestOnlineMap::GetColGlobalDoF | ( | int | localID | ) | const [inline] |
Get the global Degrees-Of-Freedom ID on the source mesh.
Definition at line 479 of file TempestOnlineMap.hpp.
{
return col_gdofmap[localColID];
}
Get the number of total Degrees-Of-Freedom defined on the destination mesh.
Get the number of local Degrees-Of-Freedom defined on the destination mesh.
| int moab::TempestOnlineMap::GetDestinationNDofsPerElement | ( | ) | [inline] |
Get the number of Degrees-Of-Freedom per element on the destination mesh.
Definition at line 497 of file TempestOnlineMap.hpp.
{
return m_nDofsPEl_Dest;
}
| const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalSourceAreas | ( | ) | const |
If we computed the reduction, get the vector representing the source areas for all entities in the mesh.
| const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalTargetAreas | ( | ) | const |
If we computed the reduction, get the vector representing the target areas for all entities in the mesh.
| int moab::TempestOnlineMap::GetIndexOfColGlobalDoF | ( | int | globalColDoF | ) | const [inline] |
Get the index of globaColDoF.
Definition at line 484 of file TempestOnlineMap.hpp.
{
return globalColDoF + 1; // temporary
}
| int moab::TempestOnlineMap::GetIndexOfRowGlobalDoF | ( | int | globalRowDoF | ) | const [inline] |
Get the index of globaRowDoF.
Definition at line 473 of file TempestOnlineMap.hpp.
{
return globalRowDoF + 1;
}
| int moab::TempestOnlineMap::GetRowGlobalDoF | ( | int | localID | ) | const [inline] |
Get the global Degrees-Of-Freedom ID on the destination mesh.
Definition at line 468 of file TempestOnlineMap.hpp.
{
return row_gdofmap[localRowID];
}
Get the number of total Degrees-Of-Freedom defined on the source mesh.
Get the number of local Degrees-Of-Freedom defined on the source mesh.
| int moab::TempestOnlineMap::GetSourceNDofsPerElement | ( | ) | [inline] |
Get the number of Degrees-Of-Freedom per element on the source mesh.
Definition at line 490 of file TempestOnlineMap.hpp.
{
return m_nDofsPEl_Src;
}
| int moab::TempestOnlineMap::IsConservative | ( | double | dTolerance | ) | [virtual] |
Determine if the map is conservative.
Definition at line 1520 of file TempestOnlineMap.cpp.
References ierr, rank, and size.
{
#ifndef MOAB_HAVE_MPI
return OfflineMap::IsConservative( dTolerance );
#else
// return OfflineMap::IsConservative(dTolerance);
int ierr;
// Get map entries
DataArray1D< int > dataRows;
DataArray1D< int > dataCols;
DataArray1D< double > dataEntries;
const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
// Calculate column sums
std::vector< int > dColumnsUnique;
std::vector< double > dColumnSums;
int nColumns = m_mapRemap.GetColumns();
m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
{
dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
assert( dataCols[i] < m_nTotDofs_SrcCov );
// GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
// int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
dColumnsUnique[dataCols[i]] = colGID;
// std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
// std::endl;
}
int rootProc = 0;
std::vector< int > nElementsInProc;
const int nDATA = 3;
if( !rank ) nElementsInProc.resize( size * nDATA );
int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
std::vector< int > dColumnIndices;
std::vector< double > dColumnSourceAreas;
std::vector< double > dColumnSumsTotal;
std::vector< int > displs, rcount;
if( rank == rootProc )
{
displs.resize( size + 1, 0 );
rcount.resize( size, 0 );
int gsum = 0;
for( int ir = 0; ir < size; ++ir )
{
nTotVals += nElementsInProc[ir * nDATA];
nTotColumns += nElementsInProc[ir * nDATA + 1];
nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
displs[ir] = gsum;
rcount[ir] = nElementsInProc[ir * nDATA + 1];
gsum += rcount[ir];
// printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns,
// displs[ir], rcount[ir], gsum);
}
dColumnIndices.resize( nTotColumns, -1 );
dColumnSumsTotal.resize( nTotColumns, 0.0 );
// dColumnSourceAreas.resize ( nTotColumns, 0.0 );
}
// Gather all ColumnSums to root process and accumulate
// We expect that the sums of all columns equate to 1.0 within user specified tolerance
// Need to do a gatherv here since different processes have different number of elements
// MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
// MPI_SUM, 0, m_pcomm->comm());
ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
// ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
// rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
// MPI_SUCCESS ) return -1;
// Clean out unwanted arrays now
dColumnSums.clear();
dColumnsUnique.clear();
// Verify all column sums equal the input Jacobian
int fConservative = 0;
if( rank == rootProc )
{
displs[size] = ( nTotColumns );
// std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
std::map< int, double > dColumnSumsOnRoot;
// std::map<int, double> dColumnSourceAreasOnRoot;
for( int ir = 0; ir < size; ir++ )
{
for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
{
if( dColumnIndices[ips] < 0 ) continue;
// printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
assert( dColumnIndices[ips] < nTotColumnsUnq );
dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
// dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
// dColumnSourceAreas[ dColumnIndices[ips] ]
}
}
for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
{
// if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
if( fabs( it->second - 1.0 ) > dTolerance )
{
fConservative++;
Announce( "TempestOnlineMap is not conservative in column "
// "%i (%1.15e)", it->first, it->second );
"%i (%1.15e)",
it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
}
}
}
// TODO: Just do a broadcast from root instead of a reduction
ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
return fConservative;
#endif
}
| int moab::TempestOnlineMap::IsConsistent | ( | double | dTolerance | ) | [virtual] |
Determine if the map is first-order accurate.
Definition at line 1474 of file TempestOnlineMap.cpp.
References ierr.
Referenced by read_map_from_disk().
{
#ifndef MOAB_HAVE_MPI
return OfflineMap::IsConsistent( dTolerance );
#else
// Get map entries
DataArray1D< int > dataRows;
DataArray1D< int > dataCols;
DataArray1D< double > dataEntries;
// Calculate row sums
DataArray1D< double > dRowSums;
m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
dRowSums.Allocate( m_mapRemap.GetRows() );
for( unsigned i = 0; i < dataRows.GetRows(); i++ )
{
dRowSums[dataRows[i]] += dataEntries[i];
}
// Verify all row sums are equal to 1
int fConsistent = 0;
for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
{
if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
{
fConsistent++;
int rowGID = row_gdofmap[i];
Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
}
}
int ierr;
int fConsistentGlobal = 0;
ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
return fConsistentGlobal;
#endif
}
| int moab::TempestOnlineMap::IsMonotone | ( | double | dTolerance | ) | [virtual] |
Determine if the map is monotone.
Definition at line 1662 of file TempestOnlineMap.cpp.
References ierr.
Referenced by read_map_from_disk().
{
#ifndef MOAB_HAVE_MPI
return OfflineMap::IsMonotone( dTolerance );
#else
// Get map entries
DataArray1D< int > dataRows;
DataArray1D< int > dataCols;
DataArray1D< double > dataEntries;
m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
// Verify all entries are in the range [0,1]
int fMonotone = 0;
for( unsigned i = 0; i < dataRows.GetRows(); i++ )
{
if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
{
fMonotone++;
Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
}
}
int ierr;
int fMonotoneGlobal = 0;
ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
if( ierr != MPI_SUCCESS ) return -1;
return fMonotoneGlobal;
#endif
}
| void moab::TempestOnlineMap::LinearRemapFVtoFV_Tempest_MOAB | ( | int | nOrder | ) | [private] |
Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh.
Definition at line 115 of file TempestLinearRemap.cpp.
References dbgprint, moab::DebugOutput::printf(), rank, and moab::DebugOutput::set_prefix().
{
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 4;
// Verify ReverseNodeArray has been calculated
if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
{
_EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInput" );
}
// Triangular quadrature rule
TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
// Number of coefficients needed at this order
#ifdef RECTANGULAR_TRUNCATION
int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif
// Number of faces you need
const int nRequiredFaceSetSize = nCoefficients;
// Fit weight exponent
const int nFitWeightsExponent = nOrder + 2;
// Announcemnets
moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
if( is_root )
{
dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
}
// Current overlap face
int ixOverlap = 0;
const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
DataArray2D< double > dIntArray;
DataArray1D< double > dConstraint( nCoefficients );
// Loop through all faces on m_meshInput
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
// Output every 1000 elements
#ifdef VERBOSE
if( ixFirst % outputFrequency == 0 && is_root )
{
dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// Find the set of Faces that overlap faceFirst
int ixOverlapBegin = ixOverlap;
unsigned ixOverlapEnd = ixOverlapBegin;
for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
{
if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 )
break;
}
unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
if( nOverlapFaces == 0 ) continue;
// Build integration array
BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
nOrder, dIntArray );
// Set of Faces to use in building the reconstruction and associated
// distance metric.
AdjacentFaceVector vecAdjFaces;
GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
// Number of adjacent Faces
int nAdjFaces = vecAdjFaces.size();
// Determine the conservative constraint equation
double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
dConstraint.Zero();
for( int p = 0; p < nCoefficients; p++ )
{
for( unsigned j = 0; j < nOverlapFaces; j++ )
{
dConstraint[p] += dIntArray[p][j];
}
dConstraint[p] /= dFirstArea;
}
// Build the fit array from the integration operator
DataArray2D< double > dFitArray;
DataArray1D< double > dFitWeights;
DataArray2D< double > dFitArrayPlus;
BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
dFitArray, dFitWeights );
// Compute the inverse fit array
bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
// Multiply integration array and fit array
DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
if( fSuccess )
{
// Multiply integration array and inverse fit array
for( int i = 0; i < nAdjFaces; i++ )
{
for( size_t j = 0; j < nOverlapFaces; j++ )
{
for( int k = 0; k < nCoefficients; k++ )
{
dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
}
}
}
// Unable to invert fit array, drop to 1st order. In this case
// dFitArrayPlus(0,0) = 1 and all other entries are zero.
}
else
{
dComposedArray.Zero();
for( size_t j = 0; j < nOverlapFaces; j++ )
{
dComposedArray( 0, j ) += dIntArray( 0, j );
}
}
// Put composed array into map
for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
{
for( unsigned j = 0; j < nOverlapFaces; j++ )
{
int& ixFirstFaceLoc = vecAdjFaces[i].first;
int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
// int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
// ixFirstFaceLoc); int ixSecondFaceGlob =
// m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
// signal to not participate, because it is a ghost target
if( ixSecondFaceLoc < 0 ) continue; // do not do anything
m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
return;
}
| void moab::TempestOnlineMap::LinearRemapFVtoGLL_MOAB | ( | const DataArray3D< int > & | dataGLLNodes, |
| const DataArray3D< double > & | dataGLLJacobian, | ||
| const DataArray1D< double > & | dataGLLNodalArea, | ||
| int | nOrder, | ||
| int | nMonotoneType, | ||
| bool | fContinuous, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite volumes to finite elements.
| void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB | ( | const DataArray3D< int > & | dataGLLNodesIn, |
| const DataArray3D< double > & | dataGLLJacobianIn, | ||
| const DataArray3D< int > & | dataGLLNodesOut, | ||
| const DataArray3D< double > & | dataGLLJacobianOut, | ||
| const DataArray1D< double > & | dataNodalAreaOut, | ||
| int | nPin, | ||
| int | nPout, | ||
| int | nMonotoneType, | ||
| bool | fContinuousIn, | ||
| bool | fContinuousOut, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite elements to finite elements.
Definition at line 894 of file TempestLinearRemap.cpp.
References center(), dbgprint, ForceIntArrayConsistencyConservation(), moab::DebugOutput::printf(), moab::DebugOutput::set_prefix(), and t.
{
// Triangular quadrature rule
TriangularQuadratureRule triquadrule( 8 );
const DataArray2D< double >& dG = triquadrule.GetG();
const DataArray1D< double >& dW = triquadrule.GetW();
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix< double >& smatMap = this->GetSparseMatrix();
// Sample coefficients
DataArray2D< double > dSampleCoeffIn( nPin, nPin );
DataArray2D< double > dSampleCoeffOut( nPout, nPout );
// Announcemnets
moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
if( is_root )
{
dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
}
// Build the integration array for each element on m_meshOverlap
DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
// Number of overlap Faces per source Face
DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
int ixOverlap = 0;
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
// Determine how many overlap Faces and triangles are present
int nOverlapFaces = 0;
size_t ixOverlapTemp = ixOverlap;
for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
{
// const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
{
break;
}
nOverlapFaces++;
}
nAllOverlapFaces[ixFirst] = nOverlapFaces;
// Increment the current overlap index
ixOverlap += nAllOverlapFaces[ixFirst];
}
// Geometric area of each output node
DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
// Area of each overlap element in the output basis
DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
// Loop through all faces on m_meshInput
ixOverlap = 0;
const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
// generic triangle used for area computation, for triangles around the center of overlap face;
// used for overlap faces with more than 4 edges;
// nodes array will be set for each triangle;
// these triangles are not part of the mesh structure, they are just temporary during
// aforementioned decomposition.
Face faceTri( 3 );
NodeVector nodes( 3 );
faceTri.SetNode( 0, 0 );
faceTri.SetNode( 1, 1 );
faceTri.SetNode( 2, 2 );
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if( ixFirst % outputFrequency == 0 && is_root )
{
dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// Quantities from the First Mesh
const Face& faceFirst = m_meshInputCov->faces[ixFirst];
const NodeVector& nodesFirst = m_meshInputCov->nodes;
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
if( !nOverlapFaces ) continue;
// // Calculate total element Jacobian
// double dTotalJacobian = 0.0;
// for (int s = 0; s < nPin; s++) {
// for (int t = 0; t < nPin; t++) {
// dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
// }
// }
// Loop through all Overlap Faces
for( int i = 0; i < nOverlapFaces; i++ )
{
// Quantities from the overlap Mesh
const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
const NodeVector& nodesOverlap = m_meshOverlap->nodes;
// Quantities from the Second Mesh
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
// signal to not participate, because it is a ghost target
if( ixSecond < 0 ) continue; // do not do anything
const NodeVector& nodesSecond = m_meshOutput->nodes;
const Face& faceSecond = m_meshOutput->faces[ixSecond];
int nbEdges = faceOverlap.edges.size();
int nOverlapTriangles = 1;
Node center; // not used if nbEdges == 3
if( nbEdges > 3 )
{ // decompose from center in this case
nOverlapTriangles = nbEdges;
for( int k = 0; k < nbEdges; k++ )
{
const Node& node = nodesOverlap[faceOverlap[k]];
center = center + node;
}
center = center / nbEdges;
center = center.Normalized(); // project back on sphere of radius 1
}
Node node0, node1, node2;
double dTriArea;
// Loop over all sub-triangles of this Overlap Face
for( int j = 0; j < nOverlapTriangles; j++ )
{
if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
{
node0 = nodesOverlap[faceOverlap[0]];
node1 = nodesOverlap[faceOverlap[1]];
node2 = nodesOverlap[faceOverlap[2]];
dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
}
else // decompose polygon in triangles around the center
{
node0 = center;
node1 = nodesOverlap[faceOverlap[j]];
int j1 = ( j + 1 ) % nbEdges;
node2 = nodesOverlap[faceOverlap[j1]];
nodes[0] = center;
nodes[1] = node1;
nodes[2] = node2;
dTriArea = CalculateFaceArea( faceTri, nodes );
}
for( int k = 0; k < triquadrule.GetPoints(); k++ )
{
// Get the nodal location of this point
double dX[3];
dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
dX[0] /= dMag;
dX[1] /= dMag;
dX[2] /= dMag;
Node nodeQuadrature( dX[0], dX[1], dX[2] );
// Find the components of this quadrature point in the basis
// of the first Face.
double dAlphaIn;
double dBetaIn;
ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
// Find the components of this quadrature point in the basis
// of the second Face.
double dAlphaOut;
double dBetaOut;
ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
/*
// Check inverse map value
if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
(dBetaIn < 0.0) || (dBetaIn > 1.0)
) {
_EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
dAlphaIn, dBetaIn);
}
// Check inverse map value
if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
(dBetaOut < 0.0) || (dBetaOut > 1.0)
) {
_EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
dAlphaOut, dBetaOut);
}
*/
// Sample the First finite element at this point
SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
// Sample the Second finite element at this point
SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
// Overlap output area
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
}
}
// Compute overlap integral
int ixp = 0;
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
int ixs = 0;
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
// Sample the Second finite element at this point
dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
ixs++;
}
}
ixp++;
}
}
}
}
}
// Coefficients
DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
for( int i = 0; i < nOverlapFaces; i++ )
{
// int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
int ixp = 0;
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
int ixs = 0;
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
dOverlapOutputArea[ixOverlap + i][s * nPout + t];
ixs++;
}
}
ixp++;
}
}
}
// Source areas
DataArray1D< double > vecSourceArea( nPin * nPin );
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
}
}
// Target areas
DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
for( int i = 0; i < nOverlapFaces; i++ )
{
// int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
int ixs = 0;
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
ixs++;
}
}
}
// Force consistency and conservation
if( !fNoConservation )
{
ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
}
// Update global coefficients
for( int i = 0; i < nOverlapFaces; i++ )
{
int ixp = 0;
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
int ixs = 0;
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
dGlobalIntArray[ixp][ixOverlap + i][ixs] =
dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
ixs++;
}
}
ixp++;
}
}
}
#ifdef VVERBOSE
// Check column sums (conservation)
for( int i = 0; i < nPin * nPin; i++ )
{
double dColSum = 0.0;
for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
{
dColSum += dCoeff[j][i] * vecTargetArea[j];
}
printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
}
// Check row sums (consistency)
for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
{
double dRowSum = 0.0;
for( int i = 0; i < nPin * nPin; i++ )
{
dRowSum += dCoeff[j][i];
}
printf( "Row %i: %1.15e\n", j, dRowSum );
}
#endif
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
// Build redistribution map within target element
if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
DataArray1D< double > dRedistSourceArea( nPout * nPout );
DataArray1D< double > dRedistTargetArea( nPout * nPout );
std::vector< DataArray2D< double > > dRedistributionMaps;
dRedistributionMaps.resize( m_meshOutput->faces.size() );
for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
{
dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
for( int i = 0; i < nPout * nPout; i++ )
{
dRedistributionMaps[ixSecond][i][i] = 1.0;
}
for( int s = 0; s < nPout * nPout; s++ )
{
dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
}
for( int s = 0; s < nPout * nPout; s++ )
{
dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
}
if( !fNoConservation )
{
ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
( nMonotoneType != 0 ) );
for( int s = 0; s < nPout * nPout; s++ )
{
for( int t = 0; t < nPout * nPout; t++ )
{
dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
}
}
}
}
// Construct the total geometric area
DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
{
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
dGeometricOutputArea[ixSecond][s * nPout + t];
}
}
}
// Compose the integration operator with the output map
ixOverlap = 0;
if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
// Map from source DOFs to target DOFs with redistribution applied
DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if( ixFirst % outputFrequency == 0 && is_root )
{
dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
if( !nOverlapFaces ) continue;
// Put composed array into map
for( int j = 0; j < nOverlapFaces; j++ )
{
int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
// signal to not participate, because it is a ghost target
if( ixSecondFace < 0 ) continue; // do not do anything
dRedistributedOp.Zero();
for( int p = 0; p < nPin * nPin; p++ )
{
for( int s = 0; s < nPout * nPout; s++ )
{
for( int t = 0; t < nPout * nPout; t++ )
{
dRedistributedOp[p][s] +=
dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
}
}
}
int ixp = 0;
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
int ixFirstNode;
if( fContinuousIn )
{
ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
}
else
{
ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
}
int ixs = 0;
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
int ixSecondNode;
if( fContinuousOut )
{
ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
if( !fNoConservation )
{
smatMap( ixSecondNode, ixFirstNode ) +=
dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
}
else
{
smatMap( ixSecondNode, ixFirstNode ) +=
dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
}
}
else
{
ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
if( !fNoConservation )
{
smatMap( ixSecondNode, ixFirstNode ) +=
dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
}
else
{
smatMap( ixSecondNode, ixFirstNode ) +=
dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
}
}
ixs++;
}
}
ixp++;
}
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
return;
}
| void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB | ( | const DataArray3D< int > & | dataGLLNodesIn, |
| const DataArray3D< double > & | dataGLLJacobianIn, | ||
| const DataArray3D< int > & | dataGLLNodesOut, | ||
| const DataArray3D< double > & | dataGLLJacobianOut, | ||
| const DataArray1D< double > & | dataNodalAreaOut, | ||
| int | nPin, | ||
| int | nPout, | ||
| int | nMonotoneType, | ||
| bool | fContinuousIn, | ||
| bool | fContinuousOut | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation).
Definition at line 1444 of file TempestLinearRemap.cpp.
References dbgprint, moab::DebugOutput::printf(), moab::DebugOutput::set_prefix(), and t.
{
// Gauss-Lobatto quadrature within Faces
DataArray1D< double > dGL;
DataArray1D< double > dWL;
GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix< double >& smatMap = this->GetSparseMatrix();
// Sample coefficients
DataArray2D< double > dSampleCoeffIn( nPin, nPin );
// Announcemnets
moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
if( is_root )
{
dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
}
// Number of overlap Faces per source Face
DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
int ixOverlap = 0;
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
size_t ixOverlapTemp = ixOverlap;
for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
{
// const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
nAllOverlapFaces[ixFirst]++;
}
// Increment the current overlap index
ixOverlap += nAllOverlapFaces[ixFirst];
}
// Number of times this point was found
DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
ixOverlap = 0;
const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
// Loop through all faces on m_meshInputCov
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if( ixFirst % outputFrequency == 0 && is_root )
{
dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// Quantities from the First Mesh
const Face& faceFirst = m_meshInputCov->faces[ixFirst];
const NodeVector& nodesFirst = m_meshInputCov->nodes;
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
// Loop through all Overlap Faces
for( int i = 0; i < nOverlapFaces; i++ )
{
// Quantities from the Second Mesh
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
// signal to not participate, because it is a ghost target
if( ixSecond < 0 ) continue; // do not do anything
const NodeVector& nodesSecond = m_meshOutput->nodes;
const Face& faceSecond = m_meshOutput->faces[ixSecond];
// Loop through all nodes on the second face
for( int s = 0; s < nPout; s++ )
{
for( int t = 0; t < nPout; t++ )
{
size_t ixSecondNode;
if( fContinuousOut )
{
ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
}
else
{
ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
}
if( ixSecondNode >= fSecondNodeFound.GetRows() )
{
_EXCEPTIONT( "Logic error" );
}
// Check if this node has been found already
if( fSecondNodeFound[ixSecondNode] )
{
continue;
}
// Check this node
Node node;
Node dDx1G;
Node dDx2G;
ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
// Find the components of this quadrature point in the basis
// of the first Face.
double dAlphaIn;
double dBetaIn;
ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
// Check if this node is within the first Face
if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
( dBetaIn > 1.0 + 1.0e-10 ) )
{
continue;
}
// Node is within the overlap region, mark as found
fSecondNodeFound[ixSecondNode] = true;
// Sample the First finite element at this point
SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
// Add to map
for( int p = 0; p < nPin; p++ )
{
for( int q = 0; q < nPin; q++ )
{
int ixFirstNode;
if( fContinuousIn )
{
ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
}
else
{
ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
}
smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
}
}
}
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
// Check for missing samples
for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
{
if( !fSecondNodeFound[i] )
{
_EXCEPTION1( "Can't sample point %i", i );
}
}
return;
}
| moab::ErrorCode moab::TempestOnlineMap::LinearRemapNN_MOAB | ( | bool | use_GID_matching = false, |
| bool | strict_check = false |
||
| ) | [private] |
Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh.
Definition at line 52 of file TempestLinearRemap.cpp.
References col_gdofmap, m_nTotDofs_Dest, m_nTotDofs_SrcCov, MB_SUCCESS, and row_gdofmap.
{
/* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov) */
#ifdef VVERBOSE
{
std::ofstream output_file( "rowcolindices.txt", std::ios::out );
output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
<< row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
output_file << "Rows \n";
for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
output_file << "Cols \n";
for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
output_file.flush(); // required here
output_file.close();
}
#endif
if( use_GID_matching )
{
std::map< unsigned, unsigned > src_gl;
for( unsigned it = 0; it < col_gdofmap.size(); ++it )
src_gl[col_gdofmap[it]] = it;
std::map< unsigned, unsigned >::iterator iter;
for( unsigned it = 0; it < row_gdofmap.size(); ++it )
{
unsigned row = row_gdofmap[it];
iter = src_gl.find( row );
if( strict_check && iter == src_gl.end() )
{
std::cout << "Searching for global target DOF " << row
<< " but could not find correspondence in source mesh.\n";
assert( false );
}
else if( iter == src_gl.end() )
{
continue;
}
else
{
unsigned icol = src_gl[row];
unsigned irow = it;
// Set the permutation matrix in local space
m_mapRemap( irow, icol ) = 1.0;
}
}
return moab::MB_SUCCESS;
}
else
{
/* Create a Kd-tree to perform local queries to find nearest neighbors */
return moab::MB_SUCCESS;
}
}
| void moab::TempestOnlineMap::LinearRemapSE0_Tempest_MOAB | ( | const DataArray3D< int > & | dataGLLNodes, |
| const DataArray3D< double > & | dataGLLJacobian | ||
| ) | [private] |
Generate the OfflineMap for linear conserative element-average spectral element to element average remapping.
| void moab::TempestOnlineMap::LinearRemapSE4_Tempest_MOAB | ( | const DataArray3D< int > & | dataGLLNodes, |
| const DataArray3D< double > & | dataGLLJacobian, | ||
| int | nMonotoneType, | ||
| bool | fContinuousIn, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping.
Definition at line 469 of file TempestLinearRemap.cpp.
References center(), dbgprint, ForceConsistencyConservation3(), moab::DebugOutput::printf(), and moab::DebugOutput::set_prefix().
{
// Order of the polynomial interpolant
int nP = dataGLLNodes.GetRows();
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 4;
// Triangular quadrature rule
TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
int TriQuadraturePoints = triquadrule.GetPoints();
const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
// Sample coefficients
DataArray2D< double > dSampleCoeff( nP, nP );
// GLL Quadrature nodes on quadrilateral elements
DataArray1D< double > dG;
DataArray1D< double > dW;
GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
// Announcements
moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
if( is_root )
{
dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
}
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix< double >& smatMap = this->GetSparseMatrix();
// NodeVector from m_meshOverlap
const NodeVector& nodesOverlap = m_meshOverlap->nodes;
const NodeVector& nodesFirst = m_meshInputCov->nodes;
// Vector of source areas
DataArray1D< double > vecSourceArea( nP * nP );
DataArray1D< double > vecTargetArea;
DataArray2D< double > dCoeff;
#ifdef VERBOSE
std::stringstream sstr;
sstr << "remapdata_" << rank << ".txt";
std::ofstream output_file( sstr.str() );
#endif
// Current Overlap Face
int ixOverlap = 0;
const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
// generic triangle used for area computation, for triangles around the center of overlap face;
// used for overlap faces with more than 4 edges;
// nodes array will be set for each triangle;
// these triangles are not part of the mesh structure, they are just temporary during
// aforementioned decomposition.
Face faceTri( 3 );
NodeVector nodes( 3 );
faceTri.SetNode( 0, 0 );
faceTri.SetNode( 1, 1 );
faceTri.SetNode( 2, 2 );
// Loop over all input Faces
for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
const Face& faceFirst = m_meshInputCov->faces[ixFirst];
if( faceFirst.edges.size() != 4 )
{
_EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
}
#ifdef VERBOSE
// Announce computation progress
if( ixFirst % outputFrequency == 0 && is_root )
{
dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
// on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
// However, the relation with MOAB and Tempest will go out of the roof
// Determine how many overlap Faces and triangles are present
int nOverlapFaces = 0;
size_t ixOverlapTemp = ixOverlap;
for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
{
// const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
{
break;
}
nOverlapFaces++;
}
// No overlaps
if( nOverlapFaces == 0 )
{
continue;
}
// Allocate remap coefficients array for meshFirst Face
DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
// Find the local remap coefficients
for( int j = 0; j < nOverlapFaces; j++ )
{
const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
{
Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
m_meshOverlap->vecFaceArea[ixOverlap + j] );
int n = faceOverlap.edges.size();
Announce( "Number nodes: %d", n );
for( int k = 0; k < n; k++ )
{
Node nd = nodesOverlap[faceOverlap[k]];
Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
}
continue;
}
// #ifdef VERBOSE
// if ( is_root )
// Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
// m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
// m_meshOverlap->vecFaceArea[ixOverlap + j] );
// #endif
int nbEdges = faceOverlap.edges.size();
int nOverlapTriangles = 1;
Node center; // not used if nbEdges == 3
if( nbEdges > 3 )
{ // decompose from center in this case
nOverlapTriangles = nbEdges;
for( int k = 0; k < nbEdges; k++ )
{
const Node& node = nodesOverlap[faceOverlap[k]];
center = center + node;
}
center = center / nbEdges;
center = center.Normalized(); // project back on sphere of radius 1
}
Node node0, node1, node2;
double dTriangleArea;
// Loop over all sub-triangles of this Overlap Face
for( int k = 0; k < nOverlapTriangles; k++ )
{
if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
{
node0 = nodesOverlap[faceOverlap[0]];
node1 = nodesOverlap[faceOverlap[1]];
node2 = nodesOverlap[faceOverlap[2]];
dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
}
else // decompose polygon in triangles around the center
{
node0 = center;
node1 = nodesOverlap[faceOverlap[k]];
int k1 = ( k + 1 ) % nbEdges;
node2 = nodesOverlap[faceOverlap[k1]];
nodes[0] = center;
nodes[1] = node1;
nodes[2] = node2;
dTriangleArea = CalculateFaceArea( faceTri, nodes );
}
// Coordinates of quadrature Node
for( int l = 0; l < TriQuadraturePoints; l++ )
{
Node nodeQuadrature;
nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
TriQuadratureG[l][2] * node2.x;
nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
TriQuadratureG[l][2] * node2.y;
nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
TriQuadratureG[l][2] * node2.z;
nodeQuadrature = nodeQuadrature.Normalized();
// Find components of quadrature point in basis
// of the first Face
double dAlpha;
double dBeta;
ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
// Check inverse map value
if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
( dBeta > 1.0 + 1.0e-13 ) )
{
_EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
"(%1.5e %1.5e)",
j, l, dAlpha, dBeta );
}
// Sample the finite element at this point
SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
// Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
m_meshOverlap->vecFaceArea[ixOverlap + j];
}
}
}
}
}
#ifdef VERBOSE
output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
for( int j = 0; j < nOverlapFaces; j++ )
{
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
output_file << dRemapCoeff[p][q][j] << " ";
}
}
}
output_file << std::endl;
#endif
// Force consistency and conservation
if( !fNoConservation )
{
double dTargetArea = 0.0;
for( int j = 0; j < nOverlapFaces; j++ )
{
dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
}
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
}
}
const double areaTolerance = 1e-10;
// Source elements are completely covered by target volumes
if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
{
vecTargetArea.Allocate( nOverlapFaces );
for( int j = 0; j < nOverlapFaces; j++ )
{
vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
}
dCoeff.Allocate( nOverlapFaces, nP * nP );
for( int j = 0; j < nOverlapFaces; j++ )
{
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
}
}
}
// Target volumes only partially cover source elements
}
else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
{
double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
vecTargetArea.Allocate( nOverlapFaces + 1 );
for( int j = 0; j < nOverlapFaces; j++ )
{
vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
}
vecTargetArea[nOverlapFaces] = dExtraneousArea;
#ifdef VERBOSE
Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
m_meshInputCov->vecFaceArea[ixFirst] );
#endif
if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
{
_EXCEPTIONT( "Partial element area exceeds total element area" );
}
dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
for( int j = 0; j < nOverlapFaces; j++ )
{
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
}
}
}
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
}
}
for( int j = 0; j < nOverlapFaces; j++ )
{
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dCoeff[nOverlapFaces][p * nP + q] -=
dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
}
}
}
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
}
}
// Source elements only partially cover target volumes
}
else
{
Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
_EXCEPTIONT( "Target grid must be a subset of source grid" );
}
ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
/*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
for( int j = 0; j < nOverlapFaces; j++ )
{
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
}
}
}
}
#ifdef VERBOSE
// output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
// for ( int j = 0; j < nOverlapFaces; j++ )
// {
// for ( int p = 0; p < nP; p++ )
// {
// for ( int q = 0; q < nP; q++ )
// {
// output_file << dRemapCoeff[p][q][j] << " ";
// }
// }
// }
// output_file << std::endl;
#endif
// Put these remap coefficients into the SparseMatrix map
for( int j = 0; j < nOverlapFaces; j++ )
{
int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
// signal to not participate, because it is a ghost target
if( ixSecondFace < 0 ) continue; // do not do anything
for( int p = 0; p < nP; p++ )
{
for( int q = 0; q < nP; q++ )
{
if( fContinuousIn )
{
int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
m_meshOverlap->vecFaceArea[ixOverlap + j] /
m_meshOutput->vecFaceArea[ixSecondFace];
}
else
{
int ixFirstNode = ixFirst * nP * nP + p * nP + q;
smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
m_meshOverlap->vecFaceArea[ixOverlap + j] /
m_meshOutput->vecFaceArea[ixSecondFace];
}
}
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
#ifdef VERBOSE
output_file.flush(); // required here
output_file.close();
#endif
return;
}
| moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap | ( | const char * | strSource, |
| const std::vector< int > & | owned_dof_ids, | ||
| bool | row_major_ownership = true |
||
| ) |
Generate the metadata associated with the offline map.
Read the OfflineMap from a NetCDF file.
Definition at line 1169 of file TempestOnlineMapIO.cpp.
References CHECK_EXCEPTION, moab::TupleList::enableWriteAccess(), moab::error(), moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_SUCCESS, rank, moab::TupleList::reset(), size, moab::TupleList::sort(), moab::TupleList::vi_rd, moab::TupleList::vi_wr, moab::TupleList::vr_rd, and moab::TupleList::vr_wr.
Referenced by read_map_from_disk().
{
NcError error( NcError::silent_nonfatal );
NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
int nS = 0, nA = 0, nB = 0;
#ifdef MOAB_HAVE_PNETCDF
// some variables will be used just in the case netcdfpar reader fails
int ncfile = -1, ret = 0;
int ndims, nvars, ngatts, unlimited;
#endif
#ifdef MOAB_HAVE_NETCDFPAR
bool is_independent = true;
ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
// ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
NcFile ncMap( strSource, NcFile::ReadOnly );
#endif
#define CHECK_EXCEPTION( obj, type, varstr ) \
{ \
if( obj == NULL ) \
{ \
_EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
} \
}
// Read SparseMatrix entries
if( ncMap.is_valid() )
{
NcDim* dimNS = ncMap.get_dim( "n_s" );
CHECK_EXCEPTION( dimNS, "dimension", "n_s" );
NcDim* dimNA = ncMap.get_dim( "n_a" );
CHECK_EXCEPTION( dimNA, "dimension", "n_a" );
NcDim* dimNB = ncMap.get_dim( "n_b" );
CHECK_EXCEPTION( dimNB, "dimension", "n_b" );
// store total number of nonzeros
nS = dimNS->size();
nA = dimNA->size();
nB = dimNB->size();
varRow = ncMap.get_var( "row" );
CHECK_EXCEPTION( varRow, "variable", "row" );
varCol = ncMap.get_var( "col" );
CHECK_EXCEPTION( varCol, "variable", "col" );
varS = ncMap.get_var( "S" );
CHECK_EXCEPTION( varS, "variable", "S" );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varRow, is_independent );
ncMap.enable_var_par_access( varCol, is_independent );
ncMap.enable_var_par_access( varS, is_independent );
#endif
}
else
{
#ifdef MOAB_HAVE_PNETCDF
// read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
// why build wth pnetcdf without MPI ?
// ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
ret = ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile );
ERR_PARNC( ret ); // bail out completely
ret = ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited );
ERR_PARNC( ret );
// find dimension ids for n_S
int ins;
ret = ncmpi_inq_dimid( ncfile, "n_s", &ins );
ERR_PARNC( ret );
MPI_Offset leng;
ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
ERR_PARNC( ret );
nS = (int)leng;
ret = ncmpi_inq_dimid( ncfile, "n_b", &ins );
ERR_PARNC( ret );
ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
ERR_PARNC( ret );
nB = (int)leng;
#else
_EXCEPTION1( "cannot read the file %s", strSource );
#endif
}
// Let us declare the map object for every process
SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
int localSize = nS / size;
long offsetRead = rank * localSize;
// leftovers on last rank
if( rank == size - 1 )
{
localSize += nS % size;
}
std::vector< int > vecRow, vecCol;
std::vector< double > vecS;
vecRow.resize( localSize );
vecCol.resize( localSize );
vecS.resize( localSize );
if( ncMap.is_valid() )
{
varRow->set_cur( (long)( offsetRead ) );
varRow->get( &( vecRow[0] ), localSize );
varCol->set_cur( (long)( offsetRead ) );
varCol->get( &( vecCol[0] ), localSize );
varS->set_cur( (long)( offsetRead ) );
varS->get( &( vecS[0] ), localSize );
ncMap.close();
}
else
{
#ifdef MOAB_HAVE_PNETCDF
// fill the local vectors with the variables from pnetcdf file; first inquire, then fill
MPI_Offset start = (MPI_Offset)offsetRead;
MPI_Offset count = (MPI_Offset)localSize;
int varid;
ret = ncmpi_inq_varid( ncfile, "S", &varid );
ERR_PARNC( ret );
ret = ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] );
ERR_PARNC( ret );
ret = ncmpi_inq_varid( ncfile, "row", &varid );
ERR_PARNC( ret );
ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] );
ERR_PARNC( ret );
ret = ncmpi_inq_varid( ncfile, "col", &varid );
ERR_PARNC( ret );
ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] );
ERR_PARNC( ret );
ret = ncmpi_close( ncfile );
ERR_PARNC( ret );
#endif
}
// Now let us set the necessary global-to-local ID maps so that A*x operations
// can be performed cleanly as if map was computed online
row_dtoc_dofmap.clear();
// row_dtoc_dofmap.reserve( nB / size );
col_dtoc_dofmap.clear();
rowMap.clear();
colMap.clear();
// col_dtoc_dofmap.reserve( 2 * nA / size );
// row_dtoc_dofmap.resize( m_nTotDofs_Dest, UINT_MAX );
// col_dtoc_dofmap.resize( m_nTotDofs_SrcCov, UINT_MAX );
#ifdef MOAB_HAVE_MPI
// bother with tuple list only if size > 1
// otherwise, just fill the sparse matrix
if( size > 1 )
{
std::vector< int > ownership;
// the default trivial partitioning scheme
int nDofs = nB; // this is for row partitioning
if( !row_partition ) nDofs = nA; // column partitioning
// assert(row_major_ownership == true); // this block is valid only for row-based partitioning
ownership.resize( size );
int nPerPart = nDofs / size;
int nRemainder = nDofs % size; // Keep the remainder in root
ownership[0] = nPerPart + nRemainder;
for( int ip = 1, roffset = ownership[0]; ip < size; ++ip )
{
roffset += nPerPart;
ownership[ip] = roffset;
}
moab::TupleList* tl = new moab::TupleList;
unsigned numr = 1; //
tl->initialize( 3, 0, 0, numr, localSize ); // to proc, row, col, value
tl->enableWriteAccess();
// populate
for( int i = 0; i < localSize; i++ )
{
int rowval = vecRow[i] - 1; // dofs are 1 based in the file
int colval = vecCol[i] - 1;
int to_proc = -1;
int dof_val = colval;
if( row_partition ) dof_val = rowval;
if( ownership[0] > dof_val )
to_proc = 0;
else
{
for( int ip = 1; ip < size; ++ip )
{
if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
{
to_proc = ip;
break;
}
}
}
int n = tl->get_n();
tl->vi_wr[3 * n] = to_proc;
tl->vi_wr[3 * n + 1] = rowval;
tl->vi_wr[3 * n + 2] = colval;
tl->vr_wr[n] = vecS[i];
tl->inc_n();
}
// heavy communication
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl, 0 );
if( owned_dof_ids.size() > 0 )
{
// we need to send desired dof to the rendez_vous point
moab::TupleList tl_re; //
tl_re.initialize( 2, 0, 0, 0, owned_dof_ids.size() ); // to proc, value
tl_re.enableWriteAccess();
// send first to rendez_vous point, decided by trivial partitioning
for( size_t i = 0; i < owned_dof_ids.size(); i++ )
{
int to_proc = -1;
int dof_val = owned_dof_ids[i] - 1; // dofs are 1 based in the file, partition from 0 ?
if( ownership[0] > dof_val )
to_proc = 0;
else
{
for( int ip = 1; ip < size; ++ip )
{
if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
{
to_proc = ip;
break;
}
}
}
int n = tl_re.get_n();
tl_re.vi_wr[2 * n] = to_proc;
tl_re.vi_wr[2 * n + 1] = dof_val;
tl_re.inc_n();
}
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl_re, 0 );
// now we know in tl_re where do we need to send back dof_val
moab::TupleList::buffer sort_buffer;
sort_buffer.buffer_init( tl_re.get_n() );
tl_re.sort( 1, &sort_buffer ); // so now we order by value
sort_buffer.buffer_init( tl->get_n() );
int indexOrder = 2; // colVal
if( row_partition ) indexOrder = 1; // rowVal
//tl->sort( indexOrder, &sort_buffer );
std::map< int, int > startDofIndex, endDofIndex; // indices in tl_re for values we want
int dofVal = -1;
if( tl_re.get_n() > 0 ) dofVal = tl_re.vi_rd[1]; // first dof val on this rank
startDofIndex[dofVal] = 0;
endDofIndex[dofVal] = 0; // start and end
for( unsigned k = 1; k < tl_re.get_n(); k++ )
{
int newDof = tl_re.vi_rd[2 * k + 1];
if( dofVal == newDof )
{
endDofIndex[dofVal] = k; // increment by 1 actually
}
else
{
dofVal = newDof;
startDofIndex[dofVal] = k;
endDofIndex[dofVal] = k;
}
}
// basically, for each value we are interested in, index in tl_re with those values are
// tl_re.vi_rd[2*startDofIndex+1] == valDof == tl_re.vi_rd[2*endDofIndex+1]
// so now we have ordered
// tl_re shows to what proc do we need to send the tuple (row, col, val)
moab::TupleList* tl_back = new moab::TupleList;
unsigned numr = 1; //
// localSize is a good guess, but maybe it should be bigger ?
// this could be bigger for repeated dofs
tl_back->initialize( 3, 0, 0, numr, tl->get_n() ); // to proc, row, col, value
tl_back->enableWriteAccess();
// now loop over tl and tl_re to see where to send
// form the new tuple, which will contain the desired dofs per task, per row or column distribution
for( unsigned k = 0; k < tl->get_n(); k++ )
{
int valDof = tl->vi_rd[3 * k + indexOrder]; // 1 for row, 2 for column // first value, it should be
for( int ire = startDofIndex[valDof]; ire <= endDofIndex[valDof]; ire++ )
{
int to_proc = tl_re.vi_rd[2 * ire];
int n = tl_back->get_n();
tl_back->vi_wr[3 * n] = to_proc;
tl_back->vi_wr[3 * n + 1] = tl->vi_rd[3 * k + 1]; // row
tl_back->vi_wr[3 * n + 2] = tl->vi_rd[3 * k + 2]; // col
tl_back->vr_wr[n] = tl->vr_rd[k];
tl_back->inc_n();
}
}
// now communicate to the desired tasks:
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl_back, 0 );
tl_re.reset(); // clear memory, although this will go out of scope
tl->reset();
tl = tl_back;
}
int rindexMax = 0, cindexMax = 0;
// populate the sparsematrix, using rowMap and colMap
int n = tl->get_n();
for( int i = 0; i < n; i++ )
{
int rindex, cindex;
const int& vecRowValue = tl->vi_wr[3 * i + 1];
const int& vecColValue = tl->vi_wr[3 * i + 2];
std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
if( riter == rowMap.end() )
{
rowMap[vecRowValue] = rindexMax;
rindex = rindexMax;
row_gdofmap.push_back( vecRowValue );
// row_dtoc_dofmap.push_back( vecRowValue );
rindexMax++;
}
else
rindex = riter->second;
std::map< int, int >::iterator citer = colMap.find( vecColValue );
if( citer == colMap.end() )
{
colMap[vecColValue] = cindexMax;
cindex = cindexMax;
col_gdofmap.push_back( vecColValue );
// col_dtoc_dofmap.push_back( vecColValue );
cindexMax++;
}
else
cindex = citer->second;
sparseMatrix( rindex, cindex ) = tl->vr_wr[i];
}
tl->reset();
}
else
#endif
{
int rindexMax = 0, cindexMax = 0;
for( int i = 0; i < nS; i++ )
{
int rindex, cindex;
const int& vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
const int& vecColValue = vecCol[i] - 1;
std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
if( riter == rowMap.end() )
{
rowMap[vecRowValue] = rindexMax;
rindex = rindexMax;
row_gdofmap.push_back( vecRowValue );
// row_dtoc_dofmap.push_back( vecRowValue );
rindexMax++;
}
else
rindex = riter->second;
std::map< int, int >::iterator citer = colMap.find( vecColValue );
if( citer == colMap.end() )
{
colMap[vecColValue] = cindexMax;
cindex = cindexMax;
col_gdofmap.push_back( vecColValue );
// col_dtoc_dofmap.push_back( vecColValue );
cindexMax++;
}
else
cindex = citer->second;
sparseMatrix( rindex, cindex ) = vecS[i];
}
}
m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
m_nTotDofs_Dest = sparseMatrix.GetRows();
#ifdef MOAB_HAVE_EIGEN3
this->copy_tempest_sparsemat_to_eigen3();
#endif
// Reset the source and target data first
m_rowVector.setZero();
m_colVector.setZero();
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs | ( | std::vector< int > & | values_entities | ) |
Definition at line 687 of file TempestOnlineMap.cpp.
References MB_SUCCESS.
{
// col_gdofmap has global dofs , that should be in the list of values, such that
// row_dtoc_dofmap[offsetDOF] = localDOF;
// // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
// we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
// form first inverse
col_dtoc_dofmap.resize( values_entities.size() );
for( int j = 0; j < (int)values_entities.size(); j++ )
{
if( colMap.find( values_entities[j] - 1 ) != colMap.end() )
col_dtoc_dofmap[j] = colMap[values_entities[j] - 1];
else
{
col_dtoc_dofmap[j] = -1; // signal that this value should not be used in
// std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
// found in colMap \n";
}
}
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs | ( | std::vector< int > & | values_entities | ) |
Definition at line 710 of file TempestOnlineMap.cpp.
References MB_SUCCESS.
{
// row_dtoc_dofmap = values_entities; // needs to point to local
// we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
row_dtoc_dofmap.resize( values_entities.size() );
for( int j = 0; j < (int)values_entities.size(); j++ )
{
if( rowMap.find( values_entities[j] - 1 ) != rowMap.end() )
row_dtoc_dofmap[j] = rowMap[values_entities[j] - 1]; // values are 1 based, but rowMap, colMap are not
else
{
row_dtoc_dofmap[j] = -1; // not all values are used
// std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
// found in rowMap \n";
}
}
return moab::MB_SUCCESS;
}
| void moab::TempestOnlineMap::SetDestinationNDofsPerElement | ( | int | nt | ) | [inline] |
Get the number of Degrees-Of-Freedom per element on the destination mesh.
Definition at line 508 of file TempestOnlineMap.hpp.
{
m_nDofsPEl_Dest = nt;
}
| moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation | ( | DiscretizationType | srcType, |
| bool | isSrcContinuous, | ||
| DataArray3D< int > * | srcdataGLLNodes, | ||
| DataArray3D< int > * | srcdataGLLNodesSrc, | ||
| DiscretizationType | destType, | ||
| bool | isDestContinuous, | ||
| DataArray3D< int > * | tgtdataGLLNodes | ||
| ) |
Compute the association between the solution tag global DoF numbering and the local matrix numbering so that matvec operations can be performed consistently.
Definition at line 214 of file TempestOnlineMap.cpp.
References ErrorCode, MB_CHK_ERR, MB_SUCCESS, and rank.
{
moab::ErrorCode rval;
std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
// We are assuming that these are element based tags that are sized: np * np
m_srcDiscType = srcType;
m_destDiscType = destType;
bool vprint = is_root && false;
#ifdef VVERBOSE
{
src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
if( is_root )
{
{
std::ofstream output_file( "sourcecov-gids-0.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
output_file << i << ", " << src_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
m_nTotDofs_SrcCov = 0;
if( isSrcContinuous )
dgll_cgll_covcol_ldofmap.resize(
m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
{
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[localDOF] = true;
}
output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
<< ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
}
}
}
output_file.flush(); // required here
output_file.close();
dgll_cgll_covcol_ldofmap.clear();
}
{
std::ofstream output_file( "source-gids-0.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
m_nTotDofs_Src = 0;
if( isSrcContinuous )
dgll_cgll_col_ldofmap.resize(
m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
{
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[localDOF] = true;
}
output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
<< ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
}
}
}
output_file.flush(); // required here
output_file.close();
dgll_cgll_col_ldofmap.clear();
}
{
std::ofstream output_file( "target-gids-0.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, GDOF, NDOF\n";
m_nTotDofs_Dest = 0;
for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
{
output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
<< m_nTotDofs_Dest << "\n";
m_nTotDofs_Dest++;
}
output_file.flush(); // required here
output_file.close();
}
}
else
{
{
std::ofstream output_file( "sourcecov-gids-1.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
output_file << i << ", " << src_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
m_nTotDofs_SrcCov = 0;
if( isSrcContinuous )
dgll_cgll_covcol_ldofmap.resize(
m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
{
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[localDOF] = true;
}
output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
<< ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
}
}
}
output_file.flush(); // required here
output_file.close();
dgll_cgll_covcol_ldofmap.clear();
}
{
std::ofstream output_file( "source-gids-1.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
m_nTotDofs_Src = 0;
if( isSrcContinuous )
dgll_cgll_col_ldofmap.resize(
m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
{
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[localDOF] = true;
}
output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
<< ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
}
}
}
output_file.flush(); // required here
output_file.close();
dgll_cgll_col_ldofmap.clear();
}
{
std::ofstream output_file( "target-gids-1.txt" );
output_file << "I, GDOF\n";
for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
output_file << "ELEMID, IDOF, GDOF, NDOF\n";
m_nTotDofs_Dest = 0;
for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
{
output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
<< m_nTotDofs_Dest << "\n";
m_nTotDofs_Dest++;
}
output_file.flush(); // required here
output_file.close();
}
}
}
#endif
// Now compute the mapping and store it for the covering mesh
int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
if( m_remapper->point_cloud_source )
{
assert( m_nDofsPEl_Src == 1 );
col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
srcTagSize = 1;
}
else
{
col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
}
// std::cout << "TOnlineMap: Process: " << rank << " and covering entities = [" <<
// col_dofmap.size() << ", " << src_soln_gdofs.size() << "]\n"; MPI_Barrier(MPI_COMM_WORLD);
#ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
unsigned maxSrcIndx = 0;
// for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
std::vector< int > locdofs( srcTagSize );
std::map< Node, moab::EntityHandle > mapLocalMBNodes;
double elcoords[3];
for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
{
EntityHandle eh = m_remapper->m_covering_source_entities[iel];
rval = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
}
const NodeVector& nodes = m_remapper->m_covering_source->nodes;
for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
{
const Face& face = m_remapper->m_covering_source->faces[j];
Node centroid;
centroid.x = centroid.y = centroid.z = 0.0;
for( unsigned l = 0; l < face.edges.size(); ++l )
{
centroid.x += nodes[face[l]].x;
centroid.y += nodes[face[l]].y;
centroid.z += nodes[face[l]].z;
}
const double factor = 1.0 / face.edges.size();
centroid.x *= factor;
centroid.y *= factor;
centroid.z *= factor;
EntityHandle current_eh;
if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
{
current_eh = mapLocalMBNodes[centroid];
}
rval = m_interface->tag_get_data( m_dofTagSrc, ¤t_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
const int offsetDOF = p * m_nDofsPEl_Src + q;
maxSrcIndx = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
<< ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
}
}
}
#endif
m_nTotDofs_SrcCov = 0;
if( srcdataGLLNodes == NULL )
{ /* we only have a mapping for elements as DoFs */
for( unsigned i = 0; i < col_gdofmap.size(); ++i )
{
assert( src_soln_gdofs[i] > 0 );
col_gdofmap[i] = src_soln_gdofs[i] - 1;
col_dtoc_dofmap[i] = i;
if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
m_nTotDofs_SrcCov++;
}
}
else
{
if( isSrcContinuous )
dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
// Put these remap coefficients into the SparseMatrix map
for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
{
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[localDOF] = true;
}
if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
assert( src_soln_gdofs[offsetDOF] > 0 );
col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
col_dtoc_dofmap[offsetDOF] = localDOF;
if( vprint )
std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", "
<< localDOF << ", " << col_gdofmap[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
}
}
}
}
if( m_remapper->point_cloud_source )
{
assert( m_nDofsPEl_Src == 1 );
srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
}
else
{
srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
}
// Now compute the mapping and store it for the original source mesh
m_nTotDofs_Src = 0;
if( srcdataGLLNodesSrc == NULL )
{ /* we only have a mapping for elements as DoFs */
for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
{
assert( locsrc_soln_gdofs[i] > 0 );
srccol_gdofmap[i] = locsrc_soln_gdofs[i] - 1;
srccol_dtoc_dofmap[i] = i;
m_nTotDofs_Src++;
}
}
else
{
if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
// Put these remap coefficients into the SparseMatrix map
for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Src; p++ )
{
for( int q = 0; q < m_nDofsPEl_Src; q++ )
{
const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
{
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[localDOF] = true;
}
if( !isSrcContinuous ) m_nTotDofs_Src++;
assert( locsrc_soln_gdofs[offsetDOF] > 0 );
srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
srccol_dtoc_dofmap[offsetDOF] = localDOF;
}
}
}
}
int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
if( m_remapper->point_cloud_target )
{
assert( m_nDofsPEl_Dest == 1 );
row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
tgtTagSize = 1;
}
else
{
row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
}
// Now compute the mapping and store it for the target mesh
// To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
m_nTotDofs_Dest = 0;
if( tgtdataGLLNodes == NULL )
{ /* we only have a mapping for elements as DoFs */
for( unsigned i = 0; i < row_gdofmap.size(); ++i )
{
assert( tgt_soln_gdofs[i] > 0 );
row_gdofmap[i] = tgt_soln_gdofs[i] - 1;
row_dtoc_dofmap[i] = i;
if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
m_nTotDofs_Dest++;
}
}
else
{
if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
// Put these remap coefficients into the SparseMatrix map
for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
{
for( int p = 0; p < m_nDofsPEl_Dest; p++ )
{
for( int q = 0; q < m_nDofsPEl_Dest; q++ )
{
const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
{
m_nTotDofs_Dest++;
dgll_cgll_row_ldofmap[localDOF] = true;
}
if( !isTgtContinuous ) m_nTotDofs_Dest++;
assert( tgt_soln_gdofs[offsetDOF] > 0 );
row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
row_dtoc_dofmap[offsetDOF] = localDOF;
if( vprint )
std::cout << "Row: " << m_remapper->lid_to_gid_tgt[j] << ", " << offsetDOF << ", " << localDOF
<< ", " << row_gdofmap[offsetDOF] << ", " << m_nTotDofs_Dest << "\n";
}
}
}
}
// Let us also allocate the local representation of the sparse matrix
#if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
if( vprint )
{
std::cout << "[" << rank << "]"
<< "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size() << ", col = " << m_nTotDofs_Src
<< ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
// std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
}
#endif
// check monotonicity of row_gdofmap and col_gdofmap
#ifdef CHECK_INCREASING_DOF
for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
{
if( row_gdofmap[i] > row_gdofmap[i + 1] )
std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
<< i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
}
for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
{
if( col_gdofmap[i] > col_gdofmap[i + 1] )
std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
<< i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
}
#endif
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags | ( | const std::string | srcDofTagName, |
| const std::string | tgtDofTagName | ||
| ) |
Store the tag names associated with global DoF ids for source and target meshes.
Definition at line 158 of file TempestOnlineMap.cpp.
References ErrorCode, ierr, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_ANY, MB_TAG_NOT_FOUND, and MB_TYPE_INTEGER.
{
moab::ErrorCode rval;
int tagSize = 0;
tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
rval =
m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
{
int ntot_elements = 0, nelements = m_remapper->m_source_entities.size();
#ifdef MOAB_HAVE_MPI
int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
#else
ntot_elements = nelements;
#endif
rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_covering_source_entities,
&m_remapper->m_source_entities, srcDofTagName, m_nDofsPEl_Src );MB_CHK_ERR( rval );
rval = m_interface->tag_get_handle( srcDofTagName.c_str(), m_nDofsPEl_Src * m_nDofsPEl_Src, MB_TYPE_INTEGER,
this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR( rval );
}
else
MB_CHK_ERR( rval );
tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
rval =
m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
{
int ntot_elements = 0, nelements = m_remapper->m_target_entities.size();
#ifdef MOAB_HAVE_MPI
int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
#else
ntot_elements = nelements;
#endif
rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_target_entities, NULL, tgtDofTagName,
m_nDofsPEl_Dest );MB_CHK_ERR( rval );
rval = m_interface->tag_get_handle( tgtDofTagName.c_str(), m_nDofsPEl_Dest * m_nDofsPEl_Dest, MB_TYPE_INTEGER,
this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR( rval );
}
else
MB_CHK_ERR( rval );
return moab::MB_SUCCESS;
}
| void moab::TempestOnlineMap::SetSourceNDofsPerElement | ( | int | ns | ) | [inline] |
Set the number of Degrees-Of-Freedom per element on the source mesh.
Definition at line 503 of file TempestOnlineMap.hpp.
{
m_nDofsPEl_Src = ns;
}
| void moab::TempestOnlineMap::setup_sizes_dimensions | ( | ) | [private] |
Definition at line 98 of file TempestOnlineMap.cpp.
References moab::TempestRemapper::RLL.
Referenced by TempestOnlineMap().
{
if( m_meshInputCov )
{
std::vector< std::string > dimNames;
std::vector< int > dimSizes;
if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
{
dimNames.push_back( "lat" );
dimNames.push_back( "lon" );
dimSizes.resize( 2, 0 );
dimSizes[0] = m_remapper->m_source_metadata[1];
dimSizes[1] = m_remapper->m_source_metadata[2];
}
else
{
dimNames.push_back( "num_elem" );
dimSizes.push_back( m_meshInputCov->faces.size() );
}
this->InitializeSourceDimensions( dimNames, dimSizes );
}
if( m_meshOutput )
{
std::vector< std::string > dimNames;
std::vector< int > dimSizes;
if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
{
dimNames.push_back( "lat" );
dimNames.push_back( "lon" );
dimSizes.resize( 2, 0 );
dimSizes[0] = m_remapper->m_target_metadata[1];
dimSizes[1] = m_remapper->m_target_metadata[2];
}
else
{
dimNames.push_back( "num_elem" );
dimSizes.push_back( m_meshOutput->faces.size() );
}
this->InitializeTargetDimensions( dimNames, dimSizes );
}
}
| moab::ErrorCode moab::TempestOnlineMap::WriteHDF5MapFile | ( | const std::string & | filename | ) | [private] |
Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.
Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.
Definition at line 764 of file TempestOnlineMapIO.cpp.
References ErrorCode, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TAG_VARLEN, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, rank, moab::TempestRemapper::RLL, root_set, size, and t.
{
moab::ErrorCode rval;
/**
* Need to get the global maximum of number of vertices per element
* Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
*dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
*when writing this out, other processes may have a different size for the same array. This is
*hence a mess to consolidate in h5mtoscrip eventually.
**/
/* Let us compute all relevant data for the current original source mesh on the process */
DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
{
this->InitializeCoordinatesFromMeshFV(
*m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
m_remapper->max_source_edges );
vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
}
else
{
DataArray3D< double > dataGLLJacobianSrc;
this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
// Generate the continuous Jacobian for input mesh
GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
if( m_srcDiscType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, m_meshInput->vecFaceArea );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianSrc, m_meshInput->vecFaceArea );
}
vecSourceFaceArea.Allocate( m_meshInput->faces.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
int offset = 0;
for( size_t e = 0; e < m_meshInput->faces.size(); e++ )
{
for( int s = 0; s < m_nDofsPEl_Src; s++ )
{
for( int t = 0; t < m_nDofsPEl_Src; t++ )
{
vecSourceFaceArea[srccol_dtoc_dofmap[offset + s * m_nDofsPEl_Src + t]] =
dataGLLJacobianSrc[s][t][e];
}
}
offset += m_nDofsPEl_Src * m_nDofsPEl_Src;
}
}
if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
{
this->InitializeCoordinatesFromMeshFV(
*m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
m_remapper->max_target_edges );
vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
}
else
{
DataArray3D< double > dataGLLJacobianDest;
this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
// Generate the continuous Jacobian for input mesh
GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
if( m_destDiscType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, m_meshOutput->vecFaceArea );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianDest, m_meshOutput->vecFaceArea );
}
vecTargetFaceArea.Allocate( m_meshOutput->faces.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
int offset = 0;
for( size_t e = 0; e < m_meshOutput->faces.size(); e++ )
{
for( int s = 0; s < m_nDofsPEl_Dest; s++ )
{
for( int t = 0; t < m_nDofsPEl_Dest; t++ )
{
vecTargetFaceArea[row_dtoc_dofmap[offset + s * m_nDofsPEl_Dest + t]] = dataGLLJacobianDest[s][t][e];
}
}
offset += m_nDofsPEl_Dest * m_nDofsPEl_Dest;
}
}
moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
int tot_src_ents = m_remapper->m_source_entities.size();
int tot_tgt_ents = m_remapper->m_target_entities.size();
int tot_src_size = dSourceCenterLon.GetRows();
int tot_tgt_size = m_dTargetCenterLon.GetRows();
int tot_vsrc_size = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
int tot_vtgt_size = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
const int weightMatNNZ = m_weightMatrix.nonZeros();
moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
moab::Tag srcAreaValues, tgtAreaValues;
rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
moab::Tag srcMaskValues, tgtMaskValues;
if( m_iSourceMask.IsAttached() )
{
rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
}
if( m_iTargetMask.IsAttached() )
{
rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
}
std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
std::vector< double > smatvals( weightMatNNZ );
// const double* smatvals = m_weightMatrix.valuePtr();
// Loop over the matrix entries and find the max global ID for rows and columns
for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
{
for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
{
smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
smatvals[offset] = it.value();
}
}
/* Set the global IDs for the DoFs */
////
// col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
// row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
////
int maxrow = 0, maxcol = 0;
std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
for( int i = 0; i < tot_src_size; ++i )
{
src_global_dofs[i] = srccol_gdofmap[i];
maxcol = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
}
for( int i = 0; i < tot_tgt_size; ++i )
{
tgt_global_dofs[i] = row_gdofmap[i];
maxrow = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
}
///////////////////////////////////////////////////////////////////////////
// The metadata in H5M file contains the following data:
//
// 1. n_a: Total source entities: (number of elements in source mesh)
// 2. n_b: Total target entities: (number of elements in target mesh)
// 3. nv_a: Max edge size of elements in source mesh
// 4. nv_b: Max edge size of elements in target mesh
// 5. maxrows: Number of rows in remap weight matrix
// 6. maxcols: Number of cols in remap weight matrix
// 7. nnz: Number of total nnz in sparse remap weight matrix
// 8. np_a: The order of the field description on the source mesh: >= 1
// 9. np_b: The order of the field description on the target mesh: >= 1
// 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
// dGLL]
// 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
// dGLL]
// 12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
// 1]
// 13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
// [0, 1, 2]
//
///////////////////////////////////////////////////////////////////////////
int map_disc_details[6];
map_disc_details[0] = m_nDofsPEl_Src;
map_disc_details[1] = m_nDofsPEl_Dest;
map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
? 0
: ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
? 0
: ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
map_disc_details[4] = ( m_bConserved ? 1 : 0 );
map_disc_details[5] = m_iMonotonicity;
#ifdef MOAB_HAVE_MPI
int loc_smatmetadata[13] = { tot_src_ents,
tot_tgt_ents,
m_remapper->max_source_edges,
m_remapper->max_target_edges,
maxrow + 1,
maxcol + 1,
weightMatNNZ,
map_disc_details[0],
map_disc_details[1],
map_disc_details[2],
map_disc_details[3],
map_disc_details[4],
map_disc_details[5] };
rval = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
int glb_smatmetadata[13] = { 0,
0,
0,
0,
0,
0,
0,
map_disc_details[0],
map_disc_details[1],
map_disc_details[2],
map_disc_details[3],
map_disc_details[4],
map_disc_details[5] };
int loc_buf[7] = {
tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
maxrow, maxcol };
int glb_buf[4] = { 0, 0, 0, 0 };
MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
glb_smatmetadata[0] = glb_buf[0];
glb_smatmetadata[1] = glb_buf[1];
glb_smatmetadata[6] = glb_buf[2];
MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
glb_smatmetadata[2] = glb_buf[0];
glb_smatmetadata[3] = glb_buf[1];
glb_smatmetadata[4] = glb_buf[2];
glb_smatmetadata[5] = glb_buf[3];
#else
int glb_smatmetadata[13] = { tot_src_ents,
tot_tgt_ents,
m_remapper->max_source_edges,
m_remapper->max_target_edges,
maxrow,
maxcol,
weightMatNNZ,
map_disc_details[0],
map_disc_details[1],
map_disc_details[2],
map_disc_details[3],
map_disc_details[4],
map_disc_details[5] };
#endif
// These values represent number of rows and columns. So should be 1-based.
glb_smatmetadata[4]++;
glb_smatmetadata[5]++;
if( this->is_root )
{
std::cout << " " << this->rank << " Writing remap weights with size [" << glb_smatmetadata[4] << " X "
<< glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
EntityHandle root_set = 0;
rval = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
}
int dsize;
const int numval = weightMatNNZ;
const void* smatrowvals_d = smatrowvals.data();
const void* smatcolvals_d = smatcolvals.data();
const void* smatvals_d = smatvals.data();
rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
/* Set the global IDs for the DoFs */
const void* srceleidvals_d = src_global_dofs.data();
const void* tgteleidvals_d = tgt_global_dofs.data();
dsize = src_global_dofs.size();
rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
dsize = tgt_global_dofs.size();
rval = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
/* Set the source and target areas */
const void* srcareavals_d = vecSourceFaceArea;
const void* tgtareavals_d = vecTargetFaceArea;
dsize = tot_src_size;
rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
dsize = tot_tgt_size;
rval = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
/* Set the coordinates for source and target center vertices */
const void* srccoordsclonvals_d = &dSourceCenterLon[0];
const void* srccoordsclatvals_d = &dSourceCenterLat[0];
dsize = dSourceCenterLon.GetRows();
rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
dsize = vecTargetFaceArea.GetRows();
rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
/* Set the coordinates for source and target element vertices */
const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
dsize = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
dsize = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
/* Set the masks for source and target meshes if available */
if( m_iSourceMask.IsAttached() )
{
const void* srcmaskvals_d = m_iSourceMask;
dsize = m_iSourceMask.GetRows();
rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
}
if( m_iTargetMask.IsAttached() )
{
const void* tgtmaskvals_d = m_iTargetMask;
dsize = m_iTargetMask.GetRows();
rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
}
#ifdef MOAB_HAVE_MPI
const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
const char* writeOptions = "";
#endif
// EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
EntityHandle sets[1] = { m_remapper->m_overlap_set };
rval = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
#ifdef WRITE_SCRIP_FILE
sstr.str( "" );
sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
std::map< std::string, std::string > mapAttributes;
mapAttributes["Creator"] = "MOAB mbtempest workflow";
if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
sstr.str( "" );
#endif
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap | ( | const std::string & | strTarget | ) |
Write the TempestOnlineMap to a parallel NetCDF file.
Definition at line 156 of file TempestOnlineMapIO.cpp.
References ErrorCode, and MB_CHK_ERR.
Referenced by main().
{
moab::ErrorCode rval;
size_t lastindex = strFilename.find_last_of( "." );
std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
// Write the map file to disk in parallel
if( extension == "nc" )
{
/* Invoke the actual call to write the parallel map to disk in SCRIP format */
rval = this->WriteSCRIPMapFile( strFilename.c_str() );MB_CHK_ERR( rval );
}
else
{
/* Write to the parallel H5M format */
rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
}
return rval;
}
| moab::ErrorCode moab::TempestOnlineMap::WriteSCRIPMapFile | ( | const std::string & | strOutputFile | ) | [private] |
Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections.
Parallel I/O with HDF5 to write out the remapping weights from multiple processors.
Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.
Definition at line 180 of file TempestOnlineMapIO.cpp.
References moab::TupleList::enableWriteAccess(), moab::error(), ErrorCode, moab::TupleList::get_n(), ierr, moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_CHK_SET_ERR, MB_SUCCESS, nr, rank, moab::TempestRemapper::RLL, size, moab::Remapper::SourceMesh, moab::Remapper::TargetMesh, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.
{
NcError error( NcError::silent_nonfatal );
#ifdef MOAB_HAVE_NETCDFPAR
bool is_independent = true;
ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
// ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
NcFile ncMap( strFilename.c_str(), NcFile::Replace );
#endif
if( !ncMap.is_valid() )
{
_EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
}
// Attributes
ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
/**
* Need to get the global maximum of number of vertices per element
* Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
*dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
*when writing this out, other processes may have a different size for the same array. This is
*hence a mess to consolidate in h5mtoscrip eventually.
**/
/* Let us compute all relevant data for the current original source mesh on the process */
DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
{
this->InitializeCoordinatesFromMeshFV(
*m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
m_remapper->max_source_edges );
vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
}
else
{
DataArray3D< double > dataGLLJacobianSrc;
this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
// Generate the continuous Jacobian for input mesh
GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
if( m_srcDiscType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
}
}
if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
{
this->InitializeCoordinatesFromMeshFV(
*m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
m_remapper->max_target_edges );
vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
{
vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
}
}
else
{
DataArray3D< double > dataGLLJacobianDest;
this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
// Generate the continuous Jacobian for input mesh
GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
if( m_destDiscType == DiscretizationType_CGLL )
{
GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
}
else
{
GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
}
}
// Map dimensions
unsigned nA = ( vecSourceFaceArea.GetRows() );
unsigned nB = ( vecTargetFaceArea.GetRows() );
std::vector< int > masksA, masksB;
ErrorCode rval = m_remapper->GetIMasks( moab::Remapper::SourceMesh, masksA );MB_CHK_SET_ERR( rval, "Trouble getting masks for source" );
rval = m_remapper->GetIMasks( moab::Remapper::TargetMesh, masksB );MB_CHK_SET_ERR( rval, "Trouble getting masks for target" );
// Number of nodes per Face
int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
// if source or target cells have triangles at poles, center of those triangles need to come from
// the original quad, not from center in 3d, converted to 2d again
// start copy OnlineMap.cpp tempestremap
// right now, do this only for source mesh; copy the logic for target mesh
for( unsigned i = 0; i < nA; i++ )
{
const Face& face = m_meshInput->faces[i];
int nNodes = face.edges.size();
int indexNodeAtPole = -1;
if( 3 == nNodes ) // check if one node at the poles
{
for( int j = 0; j < nNodes; j++ )
if( fabs( fabs( dSourceVertexLat[i][j] ) - 90.0 ) < 1.0e-12 )
{
indexNodeAtPole = j;
break;
}
}
if( indexNodeAtPole < 0 ) continue; // continue i loop, do nothing
// recompute center of cell, from 3d data; add one 2 nodes at pole, and average
int nodeAtPole = face[indexNodeAtPole]; // use the overloaded operator
Node nodePole = m_meshInput->nodes[nodeAtPole];
Node newCenter = nodePole * 2;
for( int j = 1; j < nNodes; j++ )
{
int indexi = ( indexNodeAtPole + j ) % nNodes; // nNodes is 3 !
const Node& node = m_meshInput->nodes[face[indexi]];
newCenter = newCenter + node;
}
newCenter = newCenter * 0.25;
newCenter = newCenter.Normalized();
#ifdef VERBOSE
double iniLon = dSourceCenterLon[i], iniLat = dSourceCenterLat[i];
#endif
// dSourceCenterLon, dSourceCenterLat
XYZtoRLL_Deg( newCenter.x, newCenter.y, newCenter.z, dSourceCenterLon[i], dSourceCenterLat[i] );
#ifdef VERBOSE
std::cout << " modify center of triangle from " << iniLon << " " << iniLat << " to " << dSourceCenterLon[i]
<< " " << dSourceCenterLat[i] << "\n";
#endif
}
// first move data if in parallel
#if defined( MOAB_HAVE_MPI )
int max_row_dof, max_col_dof; // output; arrays will be re-distributed in chunks [maxdof/size]
// if (size > 1)
{
int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
dSourceVertexLon, dSourceVertexLat, masksA, nA, nSourceNodesPerFace,
max_col_dof ); // now nA will be close to maxdof/size
if( ierr != 0 )
{
_EXCEPTION1( "Unable to arrange source data %d ", nA );
}
// rearrange target data: (nB)
//
ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
dTargetVertexLon, dTargetVertexLat, masksB, nB, nTargetNodesPerFace,
max_row_dof ); // now nA will be close to maxdof/size
if( ierr != 0 )
{
_EXCEPTION1( "Unable to arrange target data %d ", nB );
}
}
#endif
// Number of non-zeros in the remap matrix operator
int nS = m_weightMatrix.nonZeros();
#if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
int offbuf[3] = { 0, 0, 0 };
int globuf[5] = { 0, 0, 0, 0, 0 };
MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
// MPI_Scan is inclusive of data in current rank; modify accordingly.
offbuf[0] -= nA;
offbuf[1] -= nB;
offbuf[2] -= nS;
#else
int offbuf[3] = { 0, 0, 0 };
int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
#endif
// Write output dimensions entries
unsigned nSrcGridDims = ( m_vecSourceDimSizes.size() );
unsigned nDstGridDims = ( m_vecTargetDimSizes.size() );
NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varSrcGridDims, is_independent );
ncMap.enable_var_par_access( varDstGridDims, is_independent );
#endif
char szDim[64];
if( ( nSrcGridDims == 1 ) && ( m_vecSourceDimSizes[0] != (int)nA ) )
{
varSrcGridDims->put( &globuf[0], 1 );
varSrcGridDims->add_att( "name0", "num_dof" );
}
else
{
for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
{
int tmp = ( i == 0 ? globuf[0] : m_vecSourceDimSizes[i] );
varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
varSrcGridDims->put( &( tmp ), 1 );
}
for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
{
sprintf( szDim, "name%i", i );
varSrcGridDims->add_att( szDim, m_vecSourceDimNames[nSrcGridDims - i - 1].c_str() );
}
}
if( ( nDstGridDims == 1 ) && ( m_vecTargetDimSizes[0] != (int)nB ) )
{
varDstGridDims->put( &globuf[1], 1 );
varDstGridDims->add_att( "name0", "num_dof" );
}
else
{
for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
{
int tmp = ( i == 0 ? globuf[1] : m_vecTargetDimSizes[i] );
varDstGridDims->set_cur( nDstGridDims - i - 1 );
varDstGridDims->put( &( tmp ), 1 );
}
for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
{
sprintf( szDim, "name%i", i );
varDstGridDims->add_att( szDim, m_vecTargetDimNames[nDstGridDims - i - 1].c_str() );
}
}
// Source and Target mesh resolutions
NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
// Number of nodes per Face
NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
// Write coordinates
NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
// Write masks
NcVar* varMaskA = ncMap.add_var( "mask_a", ncInt, dimNA );
NcVar* varMaskB = ncMap.add_var( "mask_b", ncInt, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varYCA, is_independent );
ncMap.enable_var_par_access( varYCB, is_independent );
ncMap.enable_var_par_access( varXCA, is_independent );
ncMap.enable_var_par_access( varXCB, is_independent );
ncMap.enable_var_par_access( varYVA, is_independent );
ncMap.enable_var_par_access( varYVB, is_independent );
ncMap.enable_var_par_access( varXVA, is_independent );
ncMap.enable_var_par_access( varXVB, is_independent );
ncMap.enable_var_par_access( varMaskA, is_independent );
ncMap.enable_var_par_access( varMaskB, is_independent );
#endif
varYCA->add_att( "units", "degrees" );
varYCB->add_att( "units", "degrees" );
varXCA->add_att( "units", "degrees" );
varXCB->add_att( "units", "degrees" );
varYVA->add_att( "units", "degrees" );
varYVB->add_att( "units", "degrees" );
varXVA->add_att( "units", "degrees" );
varXVB->add_att( "units", "degrees" );
// Verify dimensionality
if( dSourceCenterLon.GetRows() != nA )
{
_EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" );
}
if( dSourceCenterLat.GetRows() != nA )
{
_EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" );
}
if( dTargetCenterLon.GetRows() != nB )
{
_EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" );
}
if( dTargetCenterLat.GetRows() != nB )
{
_EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" );
}
if( dSourceVertexLon.GetRows() != nA )
{
_EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" );
}
if( dSourceVertexLat.GetRows() != nA )
{
_EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" );
}
if( dTargetVertexLon.GetRows() != nB )
{
_EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" );
}
if( dTargetVertexLat.GetRows() != nB )
{
_EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" );
}
varYCA->set_cur( (long)offbuf[0] );
varYCA->put( &( dSourceCenterLat[0] ), nA );
varYCB->set_cur( (long)offbuf[1] );
varYCB->put( &( dTargetCenterLat[0] ), nB );
varXCA->set_cur( (long)offbuf[0] );
varXCA->put( &( dSourceCenterLon[0] ), nA );
varXCB->set_cur( (long)offbuf[1] );
varXCB->put( &( dTargetCenterLon[0] ), nB );
varYVA->set_cur( (long)offbuf[0] );
varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
varYVB->set_cur( (long)offbuf[1] );
varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
varXVA->set_cur( (long)offbuf[0] );
varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
varXVB->set_cur( (long)offbuf[1] );
varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
varMaskA->set_cur( (long)offbuf[0] );
varMaskA->put( &( masksA[0] ), nA );
varMaskB->set_cur( (long)offbuf[1] );
varMaskB->put( &( masksB[0] ), nB );
// Write areas
NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varAreaA, is_independent );
#endif
varAreaA->set_cur( (long)offbuf[0] );
varAreaA->put( &( vecSourceFaceArea[0] ), nA );
NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varAreaB, is_independent );
#endif
varAreaB->set_cur( (long)offbuf[1] );
varAreaB->put( &( vecTargetFaceArea[0] ), nB );
// Write SparseMatrix entries
DataArray1D< int > vecRow( nS );
DataArray1D< int > vecCol( nS );
DataArray1D< double > vecS( nS );
DataArray1D< double > dFracA( nA );
DataArray1D< double > dFracB( nB );
moab::TupleList tlValRow, tlValCol;
unsigned numr = 1; //
// value has to be sent to processor row/nB for for fracA and col/nA for fracB
// vecTargetArea (indexRow ) has to be sent for fracA (index col?)
// vecTargetFaceArea will have to be sent to col index, with its index !
tlValRow.initialize( 2, 0, 0, numr, nS ); // to proc(row), global row , value
tlValCol.initialize( 3, 0, 0, numr, nS ); // to proc(col), global row / col, value
tlValRow.enableWriteAccess();
tlValCol.enableWriteAccess();
/*
*
dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
dFracB[ row ] += val ;
*/
int offset = 0;
#if defined( MOAB_HAVE_MPI )
int nAbase = ( max_col_dof + 1 ) / size; // it is nA, except last rank ( == size - 1 )
int nBbase = ( max_row_dof + 1 ) / size; // it is nB, except last rank ( == size - 1 )
#endif
for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
{
for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
{
vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() ); // row index
vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() ); // col index
vecS[offset] = it.value(); // value
#if defined( MOAB_HAVE_MPI )
{
// value M(row, col) will contribute to procRow and procCol values for fracA and fracB
int procRow = ( vecRow[offset] - 1 ) / nBbase;
if( procRow >= size ) procRow = size - 1;
int procCol = ( vecCol[offset] - 1 ) / nAbase;
if( procCol >= size ) procCol = size - 1;
int nrInd = tlValRow.get_n();
tlValRow.vi_wr[2 * nrInd] = procRow;
tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
tlValRow.vr_wr[nrInd] = vecS[offset];
tlValRow.inc_n();
int ncInd = tlValCol.get_n();
tlValCol.vi_wr[3 * ncInd] = procCol;
tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1; // this is column
tlValCol.vr_wr[ncInd] = vecS[offset];
tlValCol.inc_n();
}
#endif
offset++;
}
}
#if defined( MOAB_HAVE_MPI )
// need to send values for their row and col processors, to compute fractions there
// now do the heavy communication
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );
// we have now, for example, dFracB[ row ] += val ;
// so we know that on current task, we received tlValRow
// reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
// dFracB[ row ] += val ;
for( unsigned i = 0; i < tlValRow.get_n(); i++ )
{
// int fromProc = tlValRow.vi_wr[2 * i];
int gRowInd = tlValRow.vi_wr[2 * i + 1];
int localIndexRow = gRowInd - nBbase * rank; // modulo nBbase rank is from 0 to size - 1;
double wgt = tlValRow.vr_wr[i];
assert( localIndexRow >= 0 );
assert( nB - localIndexRow > 0 );
dFracB[localIndexRow] += wgt;
}
// to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from
std::set< int > neededRows;
for( unsigned i = 0; i < tlValCol.get_n(); i++ )
{
int rRowInd = tlValCol.vi_wr[3 * i + 1];
neededRows.insert( rRowInd );
// we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
}
moab::TupleList tgtAreaReq;
tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
tgtAreaReq.enableWriteAccess();
for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
{
int neededRow = *sit;
int procRow = neededRow / nBbase;
if( procRow >= size ) procRow = size - 1;
int nr = tgtAreaReq.get_n();
tgtAreaReq.vi_wr[2 * nr] = procRow;
tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
tgtAreaReq.inc_n();
}
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
// we need to send back the tgtArea corresponding to row
moab::TupleList tgtAreaInfo; // load it with tgtArea at row
tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
tgtAreaInfo.enableWriteAccess();
for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
{
int from_proc = tgtAreaReq.vi_wr[2 * i];
int row = tgtAreaReq.vi_wr[2 * i + 1];
int locaIndexRow = row - rank * nBbase;
double areaToSend = vecTargetFaceArea[locaIndexRow];
// int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;
tgtAreaInfo.vi_wr[2 * i] = from_proc; // send back requested info
tgtAreaInfo.vi_wr[2 * i + 1] = row;
tgtAreaInfo.vr_wr[i] = areaToSend; // this will be tgt area at row
tgtAreaInfo.inc_n();
}
( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );
std::map< int, double > areaAtRow;
for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
{
// we have received from proc, value for row !
int row = tgtAreaInfo.vi_wr[2 * i + 1];
areaAtRow[row] = tgtAreaInfo.vr_wr[i];
}
// we have now for rows the
// it is ordered by index, so:
// now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
// tgtAreaInfo will have at index i the area we need (from row)
// there should be an easier way :(
for( unsigned i = 0; i < tlValCol.get_n(); i++ )
{
int rRowInd = tlValCol.vi_wr[3 * i + 1];
int colInd = tlValCol.vi_wr[3 * i + 2];
double val = tlValCol.vr_wr[i];
int localColInd = colInd - rank * nAbase; // < local nA
// we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
auto itMap = areaAtRow.find( rRowInd ); // it should be different from end
if( itMap != areaAtRow.end() )
{
double areaRow = itMap->second; // we fished a lot for this !
dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
}
}
#endif
// Load in data
NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varRow, is_independent );
ncMap.enable_var_par_access( varCol, is_independent );
ncMap.enable_var_par_access( varS, is_independent );
#endif
varRow->set_cur( (long)offbuf[2] );
varRow->put( vecRow, nS );
varCol->set_cur( (long)offbuf[2] );
varCol->put( vecCol, nS );
varS->set_cur( (long)offbuf[2] );
varS->put( &( vecS[0] ), nS );
// Calculate and write fractional coverage arrays
NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varFracA, is_independent );
#endif
varFracA->add_att( "name", "fraction of target coverage of source dof" );
varFracA->add_att( "units", "unitless" );
varFracA->set_cur( (long)offbuf[0] );
varFracA->put( &( dFracA[0] ), nA );
NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
ncMap.enable_var_par_access( varFracB, is_independent );
#endif
varFracB->add_att( "name", "fraction of source coverage of target dof" );
varFracB->add_att( "units", "unitless" );
varFracB->set_cur( (long)offbuf[1] );
varFracB->put( &( dFracB[0] ), nB );
// Add global attributes
// std::map<std::string, std::string>::const_iterator iterAttributes =
// mapAttributes.begin();
// for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
// ncMap.add_att(
// iterAttributes->first.c_str(),
// iterAttributes->second.c_str());
// }
ncMap.close();
return moab::MB_SUCCESS;
}
std::vector< int > moab::TempestOnlineMap::col_dtoc_dofmap [private] |
Definition at line 443 of file TempestOnlineMap.hpp.
std::vector< unsigned > moab::TempestOnlineMap::col_gdofmap [private] |
Definition at line 440 of file TempestOnlineMap.hpp.
Referenced by fill_col_ids(), and LinearRemapNN_MOAB().
std::map< int, int > moab::TempestOnlineMap::colMap [private] |
Definition at line 445 of file TempestOnlineMap.hpp.
DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesDest [private] |
Definition at line 447 of file TempestOnlineMap.hpp.
DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrc [private] |
Definition at line 447 of file TempestOnlineMap.hpp.
DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrcCov [private] |
Definition at line 447 of file TempestOnlineMap.hpp.
bool moab::TempestOnlineMap::is_parallel [private] |
Definition at line 462 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
bool moab::TempestOnlineMap::is_root [private] |
Definition at line 462 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
bool moab::TempestOnlineMap::m_bConserved [private] |
Definition at line 454 of file TempestOnlineMap.hpp.
Definition at line 448 of file TempestOnlineMap.hpp.
Definition at line 439 of file TempestOnlineMap.hpp.
moab::Tag moab::TempestOnlineMap::m_dofTagSrc [private] |
The original tag data and local to global DoF mapping to associate matrix values to solution.
Definition at line 439 of file TempestOnlineMap.hpp.
Definition at line 453 of file TempestOnlineMap.hpp.
Definition at line 453 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::m_iMonotonicity [private] |
Definition at line 455 of file TempestOnlineMap.hpp.
The reference to the moab::Core object that contains source/target and overlap sets.
Definition at line 427 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
Mesh* moab::TempestOnlineMap::m_meshInput [private] |
Definition at line 457 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
Mesh* moab::TempestOnlineMap::m_meshInputCov [private] |
Definition at line 458 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
Mesh* moab::TempestOnlineMap::m_meshOutput [private] |
Definition at line 459 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
Mesh* moab::TempestOnlineMap::m_meshOverlap [private] |
Definition at line 460 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
int moab::TempestOnlineMap::m_nDofsPEl_Dest [private] |
Definition at line 452 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::m_nDofsPEl_Src [private] |
Definition at line 452 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::m_nTotDofs_Dest [private] |
Definition at line 449 of file TempestOnlineMap.hpp.
Referenced by LinearRemapNN_MOAB().
int moab::TempestOnlineMap::m_nTotDofs_Src [private] |
Definition at line 449 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::m_nTotDofs_SrcCov [private] |
Definition at line 449 of file TempestOnlineMap.hpp.
Referenced by LinearRemapNN_MOAB().
The fundamental remapping operator object.
Definition at line 414 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
Definition at line 448 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::rank [private] |
Definition at line 463 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
int moab::TempestOnlineMap::root_proc [private] |
Definition at line 463 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
std::vector< int > moab::TempestOnlineMap::row_dtoc_dofmap [private] |
Definition at line 443 of file TempestOnlineMap.hpp.
std::vector< unsigned > moab::TempestOnlineMap::row_gdofmap [private] |
Definition at line 440 of file TempestOnlineMap.hpp.
Referenced by fill_row_ids(), and LinearRemapNN_MOAB().
std::map< int, int > moab::TempestOnlineMap::rowMap [private] |
Definition at line 445 of file TempestOnlineMap.hpp.
int moab::TempestOnlineMap::size [private] |
Definition at line 463 of file TempestOnlineMap.hpp.
Referenced by TempestOnlineMap().
std::vector< int > moab::TempestOnlineMap::srccol_dtoc_dofmap [private] |
Definition at line 443 of file TempestOnlineMap.hpp.
std::vector< unsigned > moab::TempestOnlineMap::srccol_gdofmap [private] |
Definition at line 440 of file TempestOnlineMap.hpp.
1.7.6.1