|
MOAB: Mesh Oriented datABase
(version 5.1.1)
|
|
MOAB: Mesh Oriented datABase
(version 5.1.1)
|
An offline map between two Meshes. More...
#include <TempestOfflineMap.hpp>
Collaboration diagram for moab::TempestOfflineMap:Public Types | |
| enum | DiscretizationType { DiscretizationType_FV, DiscretizationType_CGLL, DiscretizationType_DGLL } |
Public Member Functions | |
| TempestOfflineMap (moab::TempestRemapper *remapper) | |
| Generate the metadata associated with the offline map. | |
| virtual | ~TempestOfflineMap () |
| Define a virtual destructor. | |
| moab::ErrorCode | GenerateOfflineMap (std::string strInputType="fv", std::string strOutputType="fv", const int nPin=1, const int nPout=1, bool fBubble=false, int fMonotoneTypeID=0, bool fVolumetric=false, bool fNoConservation=false, bool fNoCheck=false, const std::string srcDofTagName="GLOBAL_ID", const std::string tgtDofTagName="GLOBAL_ID", const std::string strVariables="", const std::string strInputData="", const std::string strOutputData="", const std::string strNColName="", const bool fOutputDouble=false, const std::string strPreserveVariables="", const bool fPreserveAll=false, const double dFillValueOverride=0.0, const bool fInputConcave=false, const bool fOutputConcave=false) |
| 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. | |
| virtual bool | IsConsistent (double dTolerance) |
| Generate the metadata associated with the offline map. | |
| virtual bool | IsConservative (double dTolerance) |
| Determine if the map is conservative. | |
| virtual bool | IsMonotone (double dTolerance) |
| Determine if the map is monotone. | |
| const DataVector< double > & | GetGlobalSourceAreas () const |
| If we computed the reduction, get the vector representing the source areas for all entities in the mesh. | |
| const DataVector< double > & | GetGlobalTargetAreas () const |
| If we computed the reduction, get the vector representing the target areas for all entities in the mesh. | |
Public Attributes | |
| moab::TempestRemapper * | m_remapper |
| Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections. | |
| OfflineMap * | m_weightMapGlobal |
| The SparseMatrix representing this operator. | |
| bool | m_globalMapAvailable |
| The boolean flag representing whether the root process has the updated global view. | |
| moab::Interface * | mbCore |
| The DataVector that stores the global (GID-based) areas of the source mesh. | |
| 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 long > | row_dofmap |
| std::vector< unsigned long > | col_dofmap |
| std::vector< unsigned long > | srccol_dofmap |
| std::vector< unsigned long > | row_gdofmap |
| std::vector< unsigned long > | col_gdofmap |
| std::vector< unsigned long > | srccol_gdofmap |
| std::vector< unsigned long > | row_ldofmap |
| std::vector< unsigned long > | col_ldofmap |
| std::vector< unsigned long > | srccol_ldofmap |
| DataMatrix3D< int > | dataGLLNodesSrc |
| DataMatrix3D< int > | dataGLLNodesSrcCov |
| DataMatrix3D< 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 |
| Mesh * | m_meshInput |
| Mesh * | m_meshInputCov |
| Mesh * | m_meshOutput |
| Mesh * | m_meshOverlap |
| bool | is_parallel |
| bool | is_root |
| int | rank |
| int | size |
Private Member Functions | |
| moab::ErrorCode | GatherAllToRoot () |
| Gather the mapping matrix that was computed in different processors and accumulate the data on the root so that OfflineMap can be generated in parallel. | |
| 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 DataMatrix3D< int > &dataGLLNodes, const DataMatrix3D< double > &dataGLLJacobian) |
| Generate the OfflineMap for linear conserative element-average spectral element to element average remapping. | |
| void | LinearRemapSE4_Tempest_MOAB (const DataMatrix3D< int > &dataGLLNodes, const DataMatrix3D< double > &dataGLLJacobian, int nMonotoneType, bool fContinuousIn, bool fNoConservation) |
| Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping. | |
| void | LinearRemapFVtoGLL_Simple_MOAB (const DataMatrix3D< int > &dataGLLNodes, const DataMatrix3D< double > &dataGLLJacobian, const DataVector< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite volumes to finite elements using simple sampling of the FV reconstruction. | |
| void | LinearRemapFVtoGLL_Volumetric_MOAB (const DataMatrix3D< int > &dataGLLNodes, const DataMatrix3D< double > &dataGLLJacobian, const DataVector< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation) |
| Generate the OfflineMap for remapping from finite volumes to finite elements using a new experimental method. | |
| void | LinearRemapFVtoGLL_MOAB (const DataMatrix3D< int > &dataGLLNodes, const DataMatrix3D< double > &dataGLLJacobian, const DataVector< 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 DataMatrix3D< int > &dataGLLNodesIn, const DataMatrix3D< double > &dataGLLJacobianIn, const DataMatrix3D< int > &dataGLLNodesOut, const DataMatrix3D< double > &dataGLLJacobianOut, const DataVector< 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 DataMatrix3D< int > &dataGLLNodesIn, const DataMatrix3D< double > &dataGLLJacobianIn, const DataMatrix3D< int > &dataGLLNodesOut, const DataMatrix3D< double > &dataGLLJacobianOut, const DataVector< 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 | 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, DataMatrix3D< int > *srcdataGLLNodes, DataMatrix3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, bool isDestContinuous, DataMatrix3D< 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. | |
An offline map between two Meshes.
Definition at line 54 of file TempestOfflineMap.hpp.
Definition at line 71 of file TempestOfflineMap.hpp.
Generate the metadata associated with the offline map.
Definition at line 35 of file TempestOfflineMap.cpp.
References dbgprint, moab::Remapper::get_interface(), moab::TempestRemapper::GetCoveringMesh(), moab::TempestRemapper::GetMesh(), moab::Remapper::IntersectedMesh, is_parallel, is_root, m_globalMapAvailable, m_meshInput, m_meshInputCov, m_meshOutput, m_meshOverlap, m_remapper, m_weightMapGlobal, mbCore, moab::DebugOutput::printf(), rank, size, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.
: OfflineMap(), m_remapper ( remapper ) { // Get the references for the MOAB core objects mbCore = m_remapper->get_interface(); #ifdef MOAB_HAVE_MPI 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::IntersectedMesh ); m_globalMapAvailable = false; is_parallel = false; is_root = true; rank = 0; size = 1; #ifdef MOAB_HAVE_MPI int flagInit; MPI_Initialized( &flagInit ); if (flagInit) { is_parallel = true; assert(pcomm != NULL); rank = pcomm->rank(); size = pcomm->size(); is_root = (rank == 0); } #endif moab::DebugOutput dbgprint ( std::cout, rank ); // 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 std::vector<std::string> dimNames(1); std::vector<int> dimSizes(1); dimNames[0] = "num_elem"; dbgprint.printf ( 0, "Initializing dimensions of map\n" ); dbgprint.printf ( 0, "Input mesh\n" ); dimSizes[0] = m_meshInputCov->faces.size(); this->InitializeSourceDimensions(dimNames, dimSizes); dbgprint.printf ( 0, "Output mesh\n" ); dimSizes[0] = m_meshOutput->faces.size(); this->InitializeTargetDimensions(dimNames, dimSizes); m_weightMapGlobal = NULL; // 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::TempestOfflineMap::~TempestOfflineMap | ( | ) | [virtual] |
Define a virtual destructor.
Definition at line 93 of file TempestOfflineMap.cpp.
{
delete m_weightMapGlobal;
mbCore = NULL;
#ifdef MOAB_HAVE_MPI
pcomm = NULL;
#endif
m_meshInput = NULL;
m_meshOutput = NULL;
m_meshOverlap = NULL;
}
| moab::ErrorCode moab::TempestOfflineMap::GatherAllToRoot | ( | ) | [private] |
Gather the mapping matrix that was computed in different processors and accumulate the data on the root so that OfflineMap can be generated in parallel.
| moab::ErrorCode moab::TempestOfflineMap::GenerateOfflineMap | ( | std::string | strInputType = "fv", |
| std::string | strOutputType = "fv", |
||
| const int | nPin = 1, |
||
| const int | nPout = 1, |
||
| bool | fBubble = false, |
||
| int | fMonotoneTypeID = 0, |
||
| bool | fVolumetric = false, |
||
| bool | fNoConservation = false, |
||
| bool | fNoCheck = false, |
||
| const std::string | srcDofTagName = "GLOBAL_ID", |
||
| const std::string | tgtDofTagName = "GLOBAL_ID", |
||
| const std::string | strVariables = "", |
||
| const std::string | strInputData = "", |
||
| const std::string | strOutputData = "", |
||
| const std::string | strNColName = "", |
||
| const bool | fOutputDouble = false, |
||
| const std::string | strPreserveVariables = "", |
||
| const bool | fPreserveAll = false, |
||
| const double | dFillValueOverride = 0.0, |
||
| const bool | fInputConcave = false, |
||
| const bool | fOutputConcave = false |
||
| ) |
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.
Definition at line 481 of file TempestOfflineMap.cpp.
References dbgprint, moab::error(), MB_CHK_ERR, MB_FAILURE, MB_SUCCESS, ParseVariableList(), moab::DebugOutput::printf(), rank, and size.
Referenced by main().
{
NcError error ( NcError::silent_nonfatal );
moab::DebugOutput dbgprint ( std::cout, ( rank ) );
moab::ErrorCode rval;
try
{
// Check command line parameters (data arguments)
if ( ( strInputData != "" ) && ( strOutputData == "" ) )
{
_EXCEPTIONT ( "--in_data specified without --out_data" );
}
if ( ( strInputData == "" ) && ( strOutputData != "" ) )
{
_EXCEPTIONT ( "--out_data specified without --in_data" );
}
// 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
{
_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
{
_EXCEPTION1 ( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]",
strOutputType.c_str() );
}
// Monotonicity flags
int nMonotoneType = fMonotoneTypeID;
// Parse variable list
std::vector< std::string > vecVariableStrings;
ParseVariableList ( strVariables, vecVariableStrings );
// Parse preserve variable list
std::vector< std::string > vecPreserveVariableStrings;
ParseVariableList ( strPreserveVariables, vecPreserveVariableStrings );
if ( fPreserveAll && ( vecPreserveVariableStrings.size() != 0 ) )
{
_EXCEPTIONT ( "--preserveall and --preserve cannot both be specified" );
}
m_nDofsPEl_Src = nPin;
m_nDofsPEl_Dest = nPout;
rval = SetDofMapTags(srcDofTagName, tgtDofTagName);
// Calculate Face areas
if ( is_root ) dbgprint.printf ( 0, "Calculating input mesh Face areas\n" );
double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas(fInputConcave);
double dTotalAreaInput = dTotalAreaInput_loc;
#ifdef MOAB_HAVE_MPI
if (pcomm) MPI_Allreduce ( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
#endif
if ( is_root ) dbgprint.printf ( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
// Input mesh areas
m_meshInputCov->CalculateFaceAreas(fInputConcave);
if ( eInputType == DiscretizationType_FV )
{
this->SetSourceAreas ( m_meshInputCov->vecFaceArea );
}
// Calculate Face areas
if ( is_root ) dbgprint.printf ( 0, "Calculating output mesh Face areas\n" );
double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas(fOutputConcave);
double dTotalAreaOutput = dTotalAreaOutput_loc;
#ifdef MOAB_HAVE_MPI
if (pcomm) MPI_Allreduce ( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
#endif
if ( is_root ) dbgprint.printf ( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
// Output mesh areas
if ( eOutputType == DiscretizationType_FV )
{
this->SetTargetAreas ( m_meshOutput->vecFaceArea );
}
// Verify that overlap mesh is in the correct order
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->vecSourceFaceIx.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 //&&
( m_meshOutput->faces.size() - ixTargetFaceMax == 0 )
)
{
if ( is_root ) dbgprint.printf ( 0, "Overlap mesh forward correspondence found\n" );
}
else if (
// m_meshOutput->faces.size() - ixSourceFaceMax == 0 //&&
( m_meshInputCov->faces.size() - ixTargetFaceMax == 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 (pcomm) MPI_Allreduce ( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
#endif
if ( is_root ) dbgprint.printf ( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
// Partial cover
if ( fabs ( dTotalAreaOverlap - dTotalAreaInput ) > 1.0e-10 )
{
if ( !fNoCheck )
{
if ( is_root ) dbgprint.printf ( 0, "WARNING: Significant mismatch between overlap mesh area "
"and input mesh area.\n Automatically enabling --nocheck\n" );
fNoCheck = true;
}
}
/*
// 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
m_meshInputCov->ConstructReverseNodeArray();
m_meshInputCov->ConstructEdgeMap();
// 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 OfflineMap
if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
LinearRemapFVtoFV_Tempest_MOAB ( nPin );
}
else if ( eInputType == DiscretizationType_FV )
{
DataMatrix3D<double> dataGLLJacobian;
if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data\n" );
double dNumericalArea_loc =
GenerateMetaData (
*m_meshOutput,
nPout,
fBubble,
dataGLLNodesDest,
dataGLLJacobian );
double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
if (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 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, 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
m_meshInputCov->ConstructReverseNodeArray();
m_meshInputCov->ConstructEdgeMap();
// 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 ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
if ( fVolumetric )
{
LinearRemapFVtoGLL_Volumetric_MOAB (
dataGLLNodesDest,
dataGLLJacobian,
this->GetTargetAreas(),
nPin,
nMonotoneType,
fContinuous,
fNoConservation );
}
else
{
LinearRemapFVtoGLL_MOAB (
dataGLLNodesDest,
dataGLLJacobian,
this->GetTargetAreas(),
nPin,
nMonotoneType,
fContinuous,
fNoConservation );
}
}
else if (
( eInputType != DiscretizationType_FV ) &&
( eOutputType == DiscretizationType_FV )
)
{
DataMatrix3D<double> dataGLLJacobianSrc, dataGLLJacobian;
if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data\n" );
// double dNumericalAreaCov_loc =
GenerateMetaData (
*m_meshInputCov,
nPin,
fBubble,
dataGLLNodesSrcCov,
dataGLLJacobian );
double dNumericalArea_loc =
GenerateMetaData (
*m_meshInput,
nPin,
fBubble,
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 (pcomm) MPI_Allreduce ( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 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, nPin, dataGLLNodesSrcCov );
this->InitializeSourceCoordinatesFromMeshFE (
*m_meshInput, nPin, dataGLLNodesSrc );
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 offline map
if ( is_root ) dbgprint.printf ( 0, "Calculating offline map\n" );
if ( fVolumetric )
{
_EXCEPTIONT ( "Unimplemented: Volumetric currently unavailable for"
"GLL input mesh" );
}
LinearRemapSE4_Tempest_MOAB (
dataGLLNodesSrcCov,
dataGLLJacobian,
nMonotoneType,
fContinuousIn,
fNoConservation
);
}
else if (
( eInputType != DiscretizationType_FV ) &&
( eOutputType != DiscretizationType_FV )
)
{
DataMatrix3D<double> dataGLLJacobianIn, dataGLLJacobianSrc;
DataMatrix3D<double> dataGLLJacobianOut;
// Input metadata
if ( is_root ) dbgprint.printf ( 0, "Generating input mesh meta data" );
double dNumericalAreaIn_loc =
GenerateMetaData (
*m_meshInputCov,
nPin,
fBubble,
dataGLLNodesSrcCov,
dataGLLJacobianIn );
double dNumericalAreaSrc_loc =
GenerateMetaData (
*m_meshInput,
nPin,
fBubble,
dataGLLNodesSrc,
dataGLLJacobianSrc );
assert(dNumericalAreaIn_loc >= dNumericalAreaSrc_loc);
double dNumericalAreaIn = dNumericalAreaSrc_loc;
#ifdef MOAB_HAVE_MPI
if (pcomm) MPI_Allreduce ( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
#endif
if ( is_root ) dbgprint.printf ( 0, "Input Mesh Numerical Area: %1.15e", dNumericalAreaIn );
if ( fabs ( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
{
dbgprint.printf ( 0, "WARNING: Significant mismatch between input mesh "
"numerical area and geometric area" );
}
// Output metadata
if ( is_root ) dbgprint.printf ( 0, "Generating output mesh meta data" );
double dNumericalAreaOut_loc =
GenerateMetaData (
*m_meshOutput,
nPout,
fBubble,
dataGLLNodesDest,
dataGLLJacobianOut );
double dNumericalAreaOut = dNumericalAreaOut_loc;
#ifdef MOAB_HAVE_MPI
if (pcomm) MPI_Allreduce ( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, pcomm->comm() );
#endif
if ( is_root ) dbgprint.printf ( 0, "Output Mesh Numerical Area: %1.15e", 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" );
}
// Initialize coordinates for map
this->InitializeSourceCoordinatesFromMeshFE (
*m_meshInputCov, nPin, dataGLLNodesSrcCov );
this->InitializeSourceCoordinatesFromMeshFE (
*m_meshInput, nPin, dataGLLNodesSrc );
this->InitializeTargetCoordinatesFromMeshFE (
*m_meshOutput, 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 offline map
if ( is_root ) dbgprint.printf ( 0, "Calculating offline map" );
LinearRemapGLLtoGLL2_MOAB (
dataGLLNodesSrcCov,
dataGLLJacobianIn,
dataGLLNodesDest,
dataGLLJacobianOut,
this->GetTargetAreas(),
nPin,
nPout,
nMonotoneType,
fContinuousIn,
fContinuousOut,
fNoConservation
);
}
else
{
_EXCEPTIONT ( "Not implemented" );
}
#ifdef MOAB_HAVE_EIGEN
CopyTempestSparseMat_Eigen();
#endif
// Verify consistency, conservation and monotonicity
// gather weights to root process to perform consistency/conservation checks
if ( !fNoCheck && false)
{
if ( !m_globalMapAvailable && size > 1 ) {
// gather weights to root process to perform consistency/conservation checks
rval = this->GatherAllToRoot();MB_CHK_ERR(rval);
}
if ( is_root ) dbgprint.printf ( 0, "Verifying map" );
this->IsConsistent ( 1.0e-8 );
if ( !fNoConservation ) this->IsConservative ( 1.0e-8 );
if ( nMonotoneType != 0 )
{
this->IsMonotone ( 1.0e-12 );
}
}
// Apply Offline Map to data
if ( strInputData != "" )
{
if ( is_root ) dbgprint.printf ( 0, "Applying offline map to data\n" );
this->SetFillValueOverride ( static_cast<float> ( dFillValueOverride ) );
this->Apply (
strInputData,
strOutputData,
vecVariableStrings,
strNColName,
fOutputDouble,
false );
}
// Copy variables from input file to output file
if ( ( strInputData != "" ) && ( strOutputData != "" ) )
{
if ( fPreserveAll )
{
if ( is_root ) dbgprint.printf ( 0, "Preserving variables" );
this->PreserveAllVariables ( strInputData, strOutputData );
}
else if ( vecPreserveVariableStrings.size() != 0 )
{
if ( is_root ) dbgprint.printf ( 0, "Preserving variables" );
this->PreserveVariables (
strInputData,
strOutputData,
vecPreserveVariableStrings );
}
}
}
catch ( Exception & e )
{
dbgprint.printf ( 0, "%s", e.ToString().c_str() );
return ( moab::MB_FAILURE );
}
catch ( ... )
{
return ( moab::MB_FAILURE );
}
return moab::MB_SUCCESS;
}
| const DataVector< double > & moab::TempestOfflineMap::GetGlobalSourceAreas | ( | ) | const [inline] |
If we computed the reduction, get the vector representing the source areas for all entities in the mesh.
Definition at line 432 of file TempestOfflineMap.hpp.
References m_weightMapGlobal.
{
#ifdef MOAB_HAVE_MPI
if (pcomm->size() > 1) {
return m_weightMapGlobal->GetSourceAreas();
}
else {
return this->GetSourceAreas();
}
#else
return this->GetSourceAreas();
#endif
}
| const DataVector< double > & moab::TempestOfflineMap::GetGlobalTargetAreas | ( | ) | const [inline] |
If we computed the reduction, get the vector representing the target areas for all entities in the mesh.
Definition at line 446 of file TempestOfflineMap.hpp.
References m_weightMapGlobal.
{
#ifdef MOAB_HAVE_MPI
if (pcomm->size() > 1) {
return m_weightMapGlobal->GetTargetAreas();
}
else {
return this->GetTargetAreas();
}
#else
return this->GetTargetAreas();
#endif
}
| bool moab::TempestOfflineMap::IsConservative | ( | double | dTolerance | ) | [virtual] |
Determine if the map is conservative.
Definition at line 1192 of file TempestOfflineMap.cpp.
{
// Get map entries
DataVector<int> dataRows;
DataVector<int> dataCols;
DataVector<double> dataEntries;
const DataVector<double>& dTargetAreas = this->GetGlobalTargetAreas();
const DataVector<double>& dSourceAreas = this->GetGlobalSourceAreas();
// Calculate column sums
DataVector<double> dColumnSums;
if ( size > 1 )
{
if ( rank ) return true;
SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
dColumnSums.Initialize ( m_mapRemapGlobal.GetColumns() );
}
else
{
m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
dColumnSums.Initialize ( m_mapRemap.GetColumns() );
}
for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
{
dColumnSums[dataCols[i]] +=
dataEntries[i] * dTargetAreas[dataRows[i]];
}
// Verify all column sums equal the input Jacobian
bool fConservative = true;
for ( unsigned i = 0; i < dColumnSums.GetRows(); i++ )
{
if ( fabs ( dColumnSums[i] - dSourceAreas[i] ) > dTolerance )
{
fConservative = false;
Announce ( "TempestOfflineMap is not conservative in column "
"%i (%1.15e / %1.15e)",
i, dColumnSums[i], dSourceAreas[i] );
}
}
return fConservative;
}
| bool moab::TempestOfflineMap::IsConsistent | ( | double | dTolerance | ) | [virtual] |
Generate the metadata associated with the offline map.
Read the OfflineMap from a NetCDF file. Write the TempestOfflineMap to a parallel NetCDF file. Determine if the map is first-order accurate.
Definition at line 1146 of file TempestOfflineMap.cpp.
{
// Get map entries
DataVector<int> dataRows;
DataVector<int> dataCols;
DataVector<double> dataEntries;
// Calculate row sums
DataVector<double> dRowSums;
if ( size > 1 )
{
if ( rank ) return true;
SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
dRowSums.Initialize ( m_mapRemapGlobal.GetRows() );
}
else
{
m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
dRowSums.Initialize ( m_mapRemap.GetRows() );
}
for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
{
dRowSums[dataRows[i]] += dataEntries[i];
}
// Verify all row sums are equal to 1
bool fConsistent = true;
for ( unsigned i = 0; i < dRowSums.GetRows(); i++ )
{
if ( fabs ( dRowSums[i] - 1.0 ) > dTolerance )
{
fConsistent = false;
Announce ( "TempestOfflineMap is not consistent in row %i (%1.15e)",
i, dRowSums[i] );
}
}
return fConsistent;
}
| bool moab::TempestOfflineMap::IsMonotone | ( | double | dTolerance | ) | [virtual] |
Determine if the map is monotone.
Definition at line 1242 of file TempestOfflineMap.cpp.
{
// Get map entries
DataVector<int> dataRows;
DataVector<int> dataCols;
DataVector<double> dataEntries;
if ( size > 1 )
{
if ( rank ) return true;
SparseMatrix<double>& m_mapRemapGlobal = m_weightMapGlobal->GetSparseMatrix();
m_mapRemapGlobal.GetEntries ( dataRows, dataCols, dataEntries );
}
else
m_mapRemap.GetEntries ( dataRows, dataCols, dataEntries );
// Verify all entries are in the range [0,1]
bool fMonotone = true;
for ( unsigned i = 0; i < dataRows.GetRows(); i++ )
{
if ( ( dataEntries[i] < -dTolerance ) ||
( dataEntries[i] > 1.0 + dTolerance )
)
{
fMonotone = false;
Announce ( "TempestOfflineMap is not monotone in entry (%i): %1.15e",
i, dataEntries[i] );
}
}
return fMonotone;
}
| void moab::TempestOfflineMap::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 251 of file TempestLinearRemap.cpp.
References BuildFitArray(), BuildIntegrationArray(), GetAdjacentFaceVectorByEdge(), InvertFitArray_Corrected(), is_root, ixOverlap, m_meshInputCov, m_meshOutput, and m_meshOverlap.
{
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 4;
// Verify ReverseNodeArray has been calculated
if ( 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
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapFVtoFV_Tempest_MOAB] Finite Volume to Finite Volume Projection" );
Announce ( "Triangular quadrature rule order %i", TriQuadRuleOrder );
Announce ( "Number of coefficients: %i", nCoefficients );
Announce ( "Required adjacency set size: %i", nRequiredFaceSetSize );
Announce ( "Fit weights exponent: %i", nFitWeightsExponent );
}
// Current overlap face
int ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
// Loop through all faces on m_meshInput
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Output every 1000 elements
if ( ixFirst % outputFrequency == 0 )
{
Announce ( "Element %i/%i", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// This Face
// 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 ( is_root ) Announce ( "Element %i / %i :: [%i, %i]", ixFirst, m_meshInputCov->faces.size(), ixOverlapBegin, ixOverlapEnd );
if ( nOverlapFaces == 0 ) continue;
// Build integration array
DataMatrix<double> dIntArray;
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
DataVector<double> dConstraint;
dConstraint.Initialize ( nCoefficients );
double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
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
DataMatrix<double> dFitArray;
DataVector<double> dFitWeights;
DataMatrix<double> dFitArrayPlus;
BuildFitArray (
*m_meshInputCov,
triquadrule,
ixFirst,
vecAdjFaces,
nOrder,
nFitWeightsExponent,
dConstraint,
dFitArray,
dFitWeights
);
// Compute the inverse fit array
InvertFitArray_Corrected (
dConstraint,
dFitArray,
dFitWeights,
dFitArrayPlus
);
// Multiply integration array and fit array
DataMatrix<double> dComposedArray;
dComposedArray.Initialize ( nAdjFaces, nOverlapFaces );
for ( int i = 0; i < nAdjFaces; i++ )
{
for ( unsigned j = 0; j < nOverlapFaces; j++ )
{
for ( int k = 0; k < nCoefficients; k++ )
{
dComposedArray[i][j] += dIntArray[k][j] * dFitArrayPlus[i][k];
}
}
}
// 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);
m_mapRemap ( ixSecondFaceLoc, ixFirstFaceLoc ) +=
dComposedArray[i][j]
/ m_meshOutput->vecFaceArea[ixSecondFaceLoc];
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
return;
}
| void moab::TempestOfflineMap::LinearRemapFVtoGLL_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodes, |
| const DataMatrix3D< double > & | dataGLLJacobian, | ||
| const DataVector< double > & | dataGLLNodalArea, | ||
| int | nOrder, | ||
| int | nMonotoneType, | ||
| bool | fContinuous, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite volumes to finite elements.
Definition at line 2026 of file TempestLinearRemap.cpp.
References BuildFitArray(), dgesv_(), moab::GeomUtil::first(), ForceIntArrayConsistencyConservation(), GetAdjacentFaceVectorByEdge(), GetReferenceNode_MOAB(), InvertFitArray_Corrected(), ixOverlap, n, size, and t.
{
// NOTE: Reducing this quadrature rule order greatly affects error norms
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 8;
// Verify ReverseNodeArray has been calculated
if ( m_meshInputCov->revnodearray.size() == 0 )
{
_EXCEPTIONT ( "ReverseNodeArray has not been calculated for m_meshInput" );
}
if ( m_meshInputCov->edgemap.size() == 0 )
{
_EXCEPTIONT ( "EdgeMap has not been calculated for m_meshInput" );
}
// Triangular quadrature rule
TriangularQuadratureRule triquadrule ( TriQuadRuleOrder );
const DataMatrix<double> & dG = triquadrule.GetG();
const DataVector<double> & dW = triquadrule.GetW();
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix<double> & smatMap = this->GetSparseMatrix();
// Fit weight exponent
int nFitWeightsExponent = nOrder + 2;
// Order of the finite element method
int nP = dataGLLNodes.GetRows();
// Sample coefficients
DataMatrix<double> dSampleCoeff;
dSampleCoeff.Initialize ( nP, nP );
// Number of elements needed
#ifdef RECTANGULAR_TRUNCATION
int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif
int nRequiredFaceSetSize = nCoefficients;
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapFVtoGLL_MOAB] Finite Volume to Finite Element Projection" );
Announce ( "Triangular quadrature rule order %i", TriQuadRuleOrder );
Announce ( "Order of the FE polynomial interpolant: %i", nP );
Announce ( "Number of coefficients: %i", nCoefficients );
Announce ( "Required adjacency set size: %i", nRequiredFaceSetSize );
Announce ( "Fit weights exponent: %i", nFitWeightsExponent );
}
// Current overlap face
int ixOverlap = 0;
// Build the integration array for each element on m_meshOverlap
DataMatrix3D<double> dGlobalIntArray;
dGlobalIntArray.Initialize (
nCoefficients,
m_meshOverlap->faces.size(),
nP * nP );
// Number of overlap Faces per source Face
DataVector<int> nAllOverlapFaces;
nAllOverlapFaces.Initialize ( m_meshInputCov->faces.size() );
// DataVector<int> nAllTotalOverlapTriangles;
// nAllTotalOverlapTriangles.Initialize ( m_meshInputCov->faces.size() );
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]++;
// nAllTotalOverlapTriangles[ixFirst] += faceOverlap.edges.size() - 2;
}
// Increment the current overlap index
ixOverlap += nAllOverlapFaces[ixFirst];
}
// Loop through all faces on m_meshInput
ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if ( ixFirst % outputFrequency == 0 && is_root )
{
Announce ( "Element %i/%i", ixFirst, m_meshInputCov->faces.size() );
}
#endif
// This Face
const Face & faceFirst = m_meshInputCov->faces[ixFirst];
// Area of the First Face
// double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
// Coordinate axes
Node nodeRef = GetReferenceNode_MOAB ( faceFirst, m_meshInputCov->nodes );
Node nodeA1 = m_meshInputCov->nodes[faceFirst[1]] - nodeRef;
Node nodeA2 = m_meshInputCov->nodes[faceFirst[2]] - nodeRef;
Node nodeC = CrossProduct ( nodeA1, nodeA2 );
// Fit matrix
DataMatrix<double> dFit;
dFit.Initialize ( 3, 3 );
dFit[0][0] = nodeA1.x; dFit[0][1] = nodeA1.y; dFit[0][2] = nodeA1.z;
dFit[1][0] = nodeA2.x; dFit[1][1] = nodeA2.y; dFit[1][2] = nodeA2.z;
dFit[2][0] = nodeC.x; dFit[2][1] = nodeC.y; dFit[2][2] = nodeC.z;
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
// int nTotalOverlapTriangles = nAllTotalOverlapTriangles[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;
int nOverlapTriangles = faceOverlap.edges.size() - 2;
// Quantities from the Second Mesh
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
const NodeVector & nodesSecond = m_meshOutput->nodes;
const Face & faceSecond = m_meshOutput->faces[ixSecond];
// Loop over all sub-triangles of this Overlap Face
for ( int j = 0; j < nOverlapTriangles; j++ )
{
// Cornerpoints of triangle
const Node & node0 = nodesOverlap[faceOverlap[0]];
const Node & node1 = nodesOverlap[faceOverlap[j + 1]];
const Node & node2 = nodesOverlap[faceOverlap[j + 2]];
// Calculate the area of the modified Face
Face faceTri ( 3 );
faceTri.SetNode ( 0, faceOverlap[0] );
faceTri.SetNode ( 1, faceOverlap[j + 1] );
faceTri.SetNode ( 2, faceOverlap[j + 2] );
double dTriArea =
CalculateFaceArea ( faceTri, nodesOverlap );
for ( int k = 0; k < triquadrule.GetPoints(); k++ )
{
double * dGL = dG[k];
// Get the nodal location of this point
double dX[3];
dX[0] = dGL[0] * node0.x + dGL[1] * node1.x + dGL[2] * node2.x;
dX[1] = dGL[0] * node0.y + dGL[1] * node1.y + dGL[2] * node2.y;
dX[2] = dGL[0] * node0.z + dGL[1] * node1.z + dGL[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] );
dX[0] -= nodeRef.x;
dX[1] -= nodeRef.y;
dX[2] -= nodeRef.z;
// Find the coefficients for this point of the polynomial
int n = 3;
int nrhs = 1;
int lda = 3;
int ipiv[3];
int ldb = 3;
int info;
DataMatrix<double> dFitTemp;
dFitTemp = dFit;
dgesv_ (
&n, &nrhs, & ( dFitTemp[0][0] ), &lda, ipiv, dX, &ldb, &info );
// Find the components of this quadrature point in the basis
// of the finite element.
double dAlpha;
double dBeta;
ApplyInverseMap (
faceSecond,
nodesSecond,
nodeQuadrature,
dAlpha,
dBeta );
/*
// Check inverse map value
if ((dAlpha < -1.0e-12) || (dAlpha > 1.0 + 1.0e-12) ||
(dBeta < -1.0e-12) || (dBeta > 1.0 + 1.0e-12)
) {
_EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
dAlpha, dBeta);
}
*/
// Sample the finite element at this point
SampleGLLFiniteElement (
nMonotoneType,
nP,
dAlpha,
dBeta,
dSampleCoeff );
// Sample this point in the GLL element
int ixs = 0;
for ( int s = 0; s < nP; s++ )
{
for ( int t = 0; t < nP; t++ )
{
int ixp = 0;
for ( int p = 0; p < nOrder; p++ )
{
#ifdef TRIANGULAR_TRUNCATION
const int redOrder = p;
#else
const int redOrder = 0;
#endif
for ( int q = 0; q < nOrder - redOrder; q++ )
{
dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
dSampleCoeff[s][t]
* IPow ( dX[0], p )
* IPow ( dX[1], q )
* dW[k]
* dTriArea
/ dataGLLJacobian[s][t][ixSecond];
ixp++;
}
}
ixs++;
}
}
}
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
// Reverse map
std::vector< std::vector<int> > vecReverseFaceIx;
vecReverseFaceIx.resize ( m_meshOutput->faces.size() );
for ( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
{
int ixSecond = m_meshOverlap->vecTargetFaceIx[i];
vecReverseFaceIx[ixSecond].push_back ( i );
}
// Force consistency and conservation
for ( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
{
if ( vecReverseFaceIx[ixSecond].size() == 0 )
continue;
DataMatrix<double> dCoeff;
dCoeff.Initialize (
nP * nP,
vecReverseFaceIx[ixSecond].size() );
for ( size_t i = 0; i < vecReverseFaceIx[ixSecond].size(); i++ )
{
int ijxOverlap = vecReverseFaceIx[ixSecond][i];
for ( int s = 0; s < nP * nP; s++ )
{
dCoeff[s][i] = dGlobalIntArray[0][ijxOverlap][s];
}
}
// Target areas
DataVector<double> vecTargetArea;
vecTargetArea.Initialize ( nP * nP );
for ( int s = 0; s < nP * nP; s++ )
{
vecTargetArea[s] =
dataGLLJacobian[s / nP][s % nP][ixSecond];
}
// Source areas
DataVector<double> vecSourceArea;
vecSourceArea.Initialize ( vecReverseFaceIx[ixSecond].size() );
for ( size_t i = 0; i < vecReverseFaceIx[ixSecond].size(); i++ )
{
int ijxOverlap = vecReverseFaceIx[ixSecond][i];
vecSourceArea[i] = m_meshOverlap->vecFaceArea[ijxOverlap];
}
if ( !fNoConservation )
{
ForceIntArrayConsistencyConservation (
vecSourceArea,
vecTargetArea,
dCoeff,
( nMonotoneType != 0 ) );
}
for ( size_t i = 0; i < vecReverseFaceIx[ixSecond].size(); i++ )
{
int ijxOverlap = vecReverseFaceIx[ixSecond][i];
for ( int s = 0; s < nP * nP; s++ )
{
//printf("%1.15e %1.15e\n", dGlobalIntArray[0][ixOverlap][s], dCoeff[s][i]);
dGlobalIntArray[0][ijxOverlap][s] = dCoeff[s][i];
}
}
}
// Construct finite-volume fit matrix and compose with integration operator
ixOverlap = 0;
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if ( ixFirst % outputFrequency == 0 )
{
Announce ( "Element %i", ixFirst );
}
#endif
// Area of the First Face
double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
// Determine the conservative constraint equation
DataVector<double> dConstraint;
dConstraint.Initialize ( nCoefficients );
for ( int p = 0; p < nCoefficients; p++ )
{
for ( int i = 0; i < nOverlapFaces; i++ )
{
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
for ( int s = 0; s < nP * nP; s++ )
{
dConstraint[p] += dGlobalIntArray[p][ixOverlap + i][s]
* dataGLLJacobian[s / nP][s % nP][ixSecond]
/ dFirstArea;
}
}
}
// 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();
for ( int x = 0; x < nAdjFaces; x++ )
{
if ( vecAdjFaces[x].first == ( -1 ) )
{
_EXCEPTION();
}
}
// Build the fit operator
DataMatrix<double> dFitArray;
DataVector<double> dFitWeights;
DataMatrix<double> dFitArrayPlus;
BuildFitArray (
*m_meshInputCov,
triquadrule,
ixFirst,
vecAdjFaces,
nOrder,
nFitWeightsExponent,
dConstraint,
dFitArray,
dFitWeights
);
// Compute the inverse fit array
InvertFitArray_Corrected (
dConstraint,
dFitArray,
dFitWeights,
dFitArrayPlus
);
// Multiply integration array and fit array
DataMatrix<double> dComposedArray;
dComposedArray.Initialize ( nAdjFaces, nOverlapFaces * nP * nP );
for ( int j = 0; j < nOverlapFaces; j++ )
{
// int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
for ( int i = 0; i < nAdjFaces; i++ )
{
for ( int s = 0; s < nP * nP; s++ )
{
for ( int k = 0; k < nCoefficients; k++ )
{
dComposedArray[i][j * nP * nP + s] +=
dGlobalIntArray[k][ixOverlap + j][s]
* dFitArrayPlus[i][k];
}
}
}
}
// Put composed array into map
for ( size_t i = 0; i < vecAdjFaces.size(); i++ )
{
for ( int j = 0; j < nOverlapFaces; j++ )
{
int ixFirstFace = vecAdjFaces[i].first;
int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
for ( int s = 0; s < nP; s++ )
{
for ( int t = 0; t < nP; t++ )
{
int jx = j * nP * nP + s * nP + t;
if ( fContinuous )
{
int ixSecondNode = dataGLLNodes[s][t][ixSecondFace] - 1;
smatMap ( ixSecondNode, ixFirstFace ) +=
dComposedArray[i][jx]
* dataGLLJacobian[s][t][ixSecondFace]
/ dataGLLNodalArea[ixSecondNode];
}
else
{
int ixSecondNode = ixSecondFace * nP * nP + s * nP + t;
smatMap ( ixSecondNode, ixFirstFace ) +=
dComposedArray[i][jx];
}
}
}
}
}
ixOverlap += nOverlapFaces;
}
return;
}
| void moab::TempestOfflineMap::LinearRemapFVtoGLL_Simple_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodes, |
| const DataMatrix3D< double > & | dataGLLJacobian, | ||
| const DataVector< double > & | dataGLLNodalArea, | ||
| int | nOrder, | ||
| int | nMonotoneType, | ||
| bool | fContinuous, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite volumes to finite elements using simple sampling of the FV reconstruction.
Definition at line 1307 of file TempestLinearRemap.cpp.
References BuildFitArray(), dgesv_(), GetAdjacentFaceVectorByEdge(), GetReferenceNode_MOAB(), InvertFitArray_Corrected(), ixOverlap, n, and t.
{
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 8;
// Verify ReverseNodeArray has been calculated
if ( m_meshInputCov->revnodearray.size() == 0 )
{
_EXCEPTIONT ( "ReverseNodeArray has not been calculated for m_meshInput" );
}
if ( m_meshInputCov->edgemap.size() == 0 )
{
_EXCEPTIONT ( "EdgeMap has not been calculated for m_meshInput" );
}
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix<double> & smatMap = this->GetSparseMatrix();
// Fit weight exponent
int nFitWeightsExponent = nOrder + 2;
// Order of the finite element method
int nP = dataGLLNodes.GetRows();
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapFVtoGLL_Simple_MOAB] Finite Volume to Finite Element (Simple) Projection" );
Announce ( "Triangular quadrature rule order %i", TriQuadRuleOrder );
Announce ( "Order of the FE polynomial interpolant: %i", nP );
}
// Mesh utilities
MeshUtilitiesFuzzy meshutil;
// GLL nodes
DataVector<double> dG;
DataVector<double> dW;
GaussLobattoQuadrature::GetPoints ( nP, 0.0, 1.0, dG, dW );
// Triangular quadrature rule
TriangularQuadratureRule triquadrule ( TriQuadRuleOrder );
// Number of elements needed
#ifdef RECTANGULAR_TRUNCATION
int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif
int nRequiredFaceSetSize = nCoefficients;
// Set of found nodes
std::set<int> setFoundNodes;
// Loop through all faces on m_meshInput
int ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if ( ixFirst % outputFrequency == 0 )
{
Announce ( "Element %i", ixFirst );
}
#endif
// This Face
const Face & faceFirst = m_meshInputCov->faces[ixFirst];
// Coordinate axes
Node nodeRef = GetReferenceNode_MOAB ( faceFirst, m_meshInputCov->nodes );
Node nodeA1 = m_meshInputCov->nodes[faceFirst[1]] - nodeRef;
Node nodeA2 = m_meshInputCov->nodes[faceFirst[2]] - nodeRef;
Node nodeC = CrossProduct ( nodeA1, nodeA2 );
// Fit matrix
DataMatrix<double> dFit;
dFit.Initialize ( 3, 3 );
dFit[0][0] = nodeA1.x; dFit[0][1] = nodeA1.y; dFit[0][2] = nodeA1.z;
dFit[1][0] = nodeA2.x; dFit[1][1] = nodeA2.y; dFit[1][2] = nodeA2.z;
dFit[2][0] = nodeC.x; dFit[2][1] = nodeC.y; dFit[2][2] = nodeC.z;
// Set of Faces to use in building the reconstruction and associated
// distance metric.
AdjacentFaceVector vecAdjFaces;
GetAdjacentFaceVectorByEdge (
*m_meshInputCov,
ixFirst,
nRequiredFaceSetSize,
vecAdjFaces );
// Blank constraint
DataVector<double> dConstraint;
// Least squares arrays
DataMatrix<double> dFitArray;
DataVector<double> dFitWeights;
DataMatrix<double> dFitArrayPlus;
BuildFitArray (
*m_meshInputCov,
triquadrule,
ixFirst,
vecAdjFaces,
nOrder,
nFitWeightsExponent,
dConstraint,
dFitArray,
dFitWeights
);
// Compute the inverse fit array
InvertFitArray_Corrected (
dConstraint,
dFitArray,
dFitWeights,
dFitArrayPlus
);
// Number of overlapping Faces
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++;
}
// 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];
// Quantities from the Second Mesh
int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
// const NodeVector & nodesSecond = m_meshOutput->nodes;
const Face & faceSecond = m_meshOutput->faces[ixSecondFace];
for ( int s = 0; s < nP; s++ )
{
for ( int t = 0; t < nP; t++ )
{
// Determine if this Node is in faceFirst
Node node;
Node dDx1G;
Node dDx2G;
ApplyLocalMap (
faceSecond,
m_meshOutput->nodes,
dG[s],
dG[t],
node,
dDx1G,
dDx2G );
Face::NodeLocation loc;
int ixLocation;
meshutil.ContainsNode (
faceFirst,
m_meshInputCov->nodes,
node,
loc,
ixLocation );
if ( loc == Face::NodeLocation_Exterior )
{
continue;
}
// Second node index
int ixSecondNode;
if ( fContinuous )
{
ixSecondNode = dataGLLNodes[t][s][ixSecondFace] - 1;
}
else
{
ixSecondNode = ixSecondFace * nP * nP + s * nP + t;
}
// Avoid duplicates
if ( setFoundNodes.find ( ixSecondNode ) != setFoundNodes.end() )
{
continue;
}
setFoundNodes.insert ( ixSecondNode );
// Find the coefficients for this point
double dX[3];
dX[0] = node.x - nodeRef.x;
dX[1] = node.y - nodeRef.y;
dX[2] = node.z - nodeRef.z;
// if (ixSecondNode < 10) {
// double dLon = atan2(node.y, node.x);
// if (dLon < 0.0) {
// dLon += 2.0 * M_PI;
// }
// double dLat = asin(node.z);
// dLon *= 180.0 / M_PI;
// dLat *= 180.0 / M_PI;
// printf("%i %1.15e %1.15e\n", ixSecondNode, dLon, dLat);
// }
int n = 3;
int nrhs = 1;
int lda = 3;
int ipiv[3];
int ldb = 3;
int info;
DataMatrix<double> dFitTemp;
dFitTemp = dFit;
dgesv_ (
&n, &nrhs, & ( dFitTemp[0][0] ), &lda, ipiv, dX, &ldb, &info );
// Sample the reconstruction at this point
int ixp = 0;
#ifdef RECTANGULAR_TRUNCATION
for ( int p = 0; p < nOrder; p++ )
{
for ( int q = 0; q < nOrder; q++ )
{
#endif
#ifdef TRIANGULAR_TRUNCATION
for ( int p = 0; p < nOrder; p++ )
{
for ( int q = 0; q < nOrder - p; q++ )
{
#endif
for ( size_t nx = 0; nx < vecAdjFaces.size(); nx++ )
{
int ixAdjFace = vecAdjFaces[nx].first;
smatMap ( ixSecondNode, ixAdjFace ) +=
IPow ( dX[0], p )
* IPow ( dX[1], q )
* dFitArrayPlus[nx][ixp];
}
ixp++;
}
}
}
}
}
| void moab::TempestOfflineMap::LinearRemapFVtoGLL_Volumetric_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodes, |
| const DataMatrix3D< double > & | dataGLLJacobian, | ||
| const DataVector< double > & | dataGLLNodalArea, | ||
| int | nOrder, | ||
| int | nMonotoneType, | ||
| bool | fContinuous, | ||
| bool | fNoConservation | ||
| ) | [private] |
Generate the OfflineMap for remapping from finite volumes to finite elements using a new experimental method.
Definition at line 1601 of file TempestLinearRemap.cpp.
References BuildFitArray(), BuildIntegrationArray(), ForceIntArrayConsistencyConservation(), GetAdjacentFaceVectorByEdge(), InvertFitArray_Corrected(), ixOverlap, and n.
{
// Order of triangular quadrature rule
const int TriQuadRuleOrder = 4;
// Verify ReverseNodeArray has been calculated
if ( m_meshInputCov->revnodearray.size() == 0 )
{
_EXCEPTIONT ( "ReverseNodeArray has not been calculated for m_meshInput" );
}
// Triangular quadrature rule
TriangularQuadratureRule triquadrule ( TriQuadRuleOrder );
// Fit weight exponent
int nFitWeightsExponent = nOrder + 2;
// Order of the finite element method
int nP = dataGLLNodes.GetRows();
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapFVtoGLL_Volumetric_MOAB] Finite Volume to Finite Element (Volumetric) Projection" );
Announce ( "Triangular quadrature rule order %i", TriQuadRuleOrder );
Announce ( "Order of the FE polynomial interpolant: %i", nP );
}
// Gauss-Lobatto quadrature nodes and weights
DataVector<double> dG;
DataVector<double> dW;
GaussLobattoQuadrature::GetPoints ( nP, 0.0, 1.0, dG, dW );
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix<double> & smatMap = this->GetSparseMatrix();
// Number of elements needed
#ifdef RECTANGULAR_TRUNCATION
int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif
// #pragma message "This should be a command-line parameter"
int nRequiredFaceSetSize = nCoefficients;
// Accumulated weight vector
DataVector<double> dAccumW ( nP + 1 );
dAccumW[0] = 0.0;
for ( int i = 1; i < nP + 1; i++ )
{
dAccumW[i] = dAccumW[i - 1] + dW[i - 1];
}
if ( fabs ( dAccumW[dAccumW.GetRows() - 1] - 1.0 ) > 1.0e-14 )
{
_EXCEPTIONT ( "Logic error in accumulated weight" );
}
// Create sub-element mesh and redistribution map
Announce ( "Generating sub-element mesh" );
Mesh meshTargetSubElement;
DataVector<double> dFiniteVolumeArea ( nP * nP );
DataVector<double> dQuadratureArea ( nP * nP );
std::vector< DataMatrix<double> > dRedistributionMaps;
dRedistributionMaps.resize ( m_meshOutput->faces.size() );
for ( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
{
const Face & faceSecond = m_meshOutput->faces[ixSecond];
const Node & nodeOutput0 = m_meshOutput->nodes[faceSecond[0]];
const Node & nodeOutput1 = m_meshOutput->nodes[faceSecond[1]];
const Node & nodeOutput2 = m_meshOutput->nodes[faceSecond[2]];
const Node & nodeOutput3 = m_meshOutput->nodes[faceSecond[3]];
for ( int q = 0; q < nP; q++ )
{
for ( int p = 0; p < nP; p++ )
{
Node node0 =
InterpolateQuadrilateralNode (
nodeOutput0, nodeOutput1, nodeOutput2, nodeOutput3,
dAccumW[p], dAccumW[q] );
Node node1 =
InterpolateQuadrilateralNode (
nodeOutput0, nodeOutput1, nodeOutput2, nodeOutput3,
dAccumW[p + 1], dAccumW[q] );
Node node2 =
InterpolateQuadrilateralNode (
nodeOutput0, nodeOutput1, nodeOutput2, nodeOutput3,
dAccumW[p + 1], dAccumW[q + 1] );
Node node3 =
InterpolateQuadrilateralNode (
nodeOutput0, nodeOutput1, nodeOutput2, nodeOutput3,
dAccumW[p], dAccumW[q + 1] );
int nNodeStart = meshTargetSubElement.nodes.size();
meshTargetSubElement.nodes.push_back ( node0 );
meshTargetSubElement.nodes.push_back ( node1 );
meshTargetSubElement.nodes.push_back ( node2 );
meshTargetSubElement.nodes.push_back ( node3 );
Face faceNew ( 4 );
faceNew.SetNode ( 0, nNodeStart );
faceNew.SetNode ( 1, nNodeStart + 1 );
faceNew.SetNode ( 2, nNodeStart + 2 );
faceNew.SetNode ( 3, nNodeStart + 3 );
meshTargetSubElement.faces.push_back ( faceNew );
dFiniteVolumeArea[q * nP + p] =
CalculateFaceArea (
faceNew, meshTargetSubElement.nodes );
dQuadratureArea[q * nP + p] =
dataGLLJacobian[q][p][ixSecond];
}
}
dRedistributionMaps[ixSecond].Initialize ( nP * nP, nP * nP );
for ( int i = 0; i < nP * nP; i++ )
{
dRedistributionMaps[ixSecond][i][i] = 1.0;
}
if ( !fNoConservation )
{
ForceIntArrayConsistencyConservation (
dFiniteVolumeArea,
dQuadratureArea,
dRedistributionMaps[ixSecond],
( nMonotoneType != 0 ) );
}
/*
double dSumQuadArea = 0.0;
double dSumFVArea = 0.0;
for (int i = 0; i < nP * nP; i++) {
dSumQuadArea += dQuadratureArea[i];
dSumFVArea += dFiniteVolumeArea[i];
}
if (fabs(dSumQuadArea - dSumFVArea) > 1.0e-14) {
printf("%1.15e\n", dSumQuadArea - dSumFVArea);
_EXCEPTION();
}
*/
for ( int i = 0; i < nP * nP; i++ )
{
for ( int j = 0; j < nP * nP; j++ )
{
dRedistributionMaps[ixSecond][i][j] *=
dQuadratureArea[i] / dFiniteVolumeArea[j];
}
}
/*
for (int i = 0; i < nP * nP; i++) {
double dSum = 0.0;
for (int j = 0; j < nP * nP; j++) {
dSum += dRedistributionMaps[ixSecond][i][j] * dFiniteVolumeArea[j];
}
printf("%i %1.15e %1.15e\n", i, dSum, dQuadratureArea[i]);
}
_EXCEPTION();
*/
/*
for (int i = 0; i < nP * nP; i++) {
for (int j = 0; j < nP * nP; j++) {
printf("%i %i %1.15e\n", i, j, dRedistributionMaps[ixSecond][i][j]);
}
}
_EXCEPTION();
*/
}
// Current overlap face
int ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
// Loop through all faces on m_meshInput
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if ( ixFirst % outputFrequency == 0 )
{
Announce ( "Element %i", ixFirst );
}
#endif
// This Face
// const Face & faceFirst = m_meshInputCov->faces[ixFirst];
// Find the set of Faces that overlap faceFirst
size_t ixOverlapBegin = ixOverlap;
size_t ixOverlapEnd = ixOverlapBegin;
for ( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
{
if ( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 )
break;
}
size_t nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
// Create a new Mesh representing the division of target finite
// elements associated with this finite volume.
Mesh meshThisElement;
meshThisElement.faces.reserve ( nOverlapFaces * nP * nP );
meshThisElement.vecTargetFaceIx.reserve ( nOverlapFaces * nP * nP );
for ( size_t i = ixOverlapBegin; i < ixOverlapEnd; i++ )
{
int iTargetFace = m_meshOverlap->vecTargetFaceIx[i];
int iSubElementBegin = iTargetFace * nP * nP;
// int iSubElementEnd = ( iTargetFace + 1 ) * nP * nP;
int iSubElement = iSubElementBegin;
for ( int p = 0; p < nP; p++ )
{
for ( int q = 0; q < nP; q++ )
{
// Calculate overlap polygon between sub-element
// and finite volume
NodeVector nodevecOutput;
GenerateOverlapFace<MeshUtilitiesFuzzy, Node> (
*m_meshInputCov,
meshTargetSubElement,
ixFirst,
iSubElementBegin + p * nP + q,
nodevecOutput );
/*
if (nodevecOutput.size() < 3) {
continue;
}
Face faceNew(nodevecOutput.size());
for (int n = 0; n < nodevecOutput.size(); n++) {
faceNew.SetNode(n, n);
}
Real dArea = CalculateFaceArea(faceNew, nodevecOutput);
if (dArea < 1.0e-13) {
continue;
}
*/
/*
if (dataGLLNodes[p][q][iTargetFace] - 1 == 3) {
printf("%1.15e %1.15e %1.15e\n", dArea, dataGLLJacobian[p][q][iTargetFace], dataGLLNodalArea[dataGLLNodes[p][q][iTargetFace] - 1]);
}
*/
Face faceNew ( nodevecOutput.size() );
for ( size_t n = 0; n < nodevecOutput.size(); n++ )
{
meshThisElement.nodes.push_back ( nodevecOutput[n] );
faceNew.SetNode ( n, meshThisElement.nodes.size() - 1 );
}
meshThisElement.faces.push_back ( faceNew );
meshThisElement.vecTargetFaceIx.push_back (
dataGLLNodes[p][q][iTargetFace] - 1 );
iSubElement++;
}
}
}
if ( meshThisElement.faces.size() != nOverlapFaces * nP * nP )
{
_EXCEPTIONT ( "Logic error" );
}
// Build integration array
DataMatrix<double> dIntArray;
BuildIntegrationArray (
*m_meshInputCov,
meshThisElement,
triquadrule,
ixFirst,
0,
meshThisElement.faces.size(),
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
DataVector<double> dConstraint;
dConstraint.Initialize ( nCoefficients );
double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
for ( int p = 0; p < nCoefficients; p++ )
{
for ( size_t j = 0; j < meshThisElement.faces.size(); j++ )
{
dConstraint[p] += dIntArray[p][j];
}
dConstraint[p] /= dFirstArea;
}
// Least squares arrays
DataMatrix<double> dFitArray;
DataVector<double> dFitWeights;
DataMatrix<double> dFitArrayPlus;
BuildFitArray (
*m_meshInputCov,
triquadrule,
ixFirst,
vecAdjFaces,
nOrder,
nFitWeightsExponent,
dConstraint,
dFitArray,
dFitWeights
);
// Compute the inverse fit array
InvertFitArray_Corrected (
dConstraint,
dFitArray,
dFitWeights,
dFitArrayPlus
);
// Multiply integration array and fit array
DataMatrix<double> dComposedArray;
dComposedArray.Initialize ( nAdjFaces, meshThisElement.faces.size() );
for ( int i = 0; i < nAdjFaces; i++ )
{
for ( size_t j = 0; j < meshThisElement.faces.size(); j++ )
{
for ( int k = 0; k < nCoefficients; k++ )
{
dComposedArray[i][j] += dIntArray[k][j] * dFitArrayPlus[i][k];
}
}
}
// Apply redistribution operator
DataMatrix<double> dRedistributedArray;
dRedistributedArray.Initialize ( nAdjFaces, meshThisElement.faces.size() );
for ( int i = 0; i < nAdjFaces; i++ )
{
for ( size_t j = 0; j < meshThisElement.faces.size(); j++ )
{
int ixSubElement = j % ( nP * nP );
int ixElement = j / ( nP * nP );
int ixSecondFace =
m_meshOverlap->vecTargetFaceIx[ixOverlap + ixElement];
for ( int k = 0; k < nP * nP; k++ )
{
dRedistributedArray[i][j] +=
dComposedArray[i][ixElement * nP * nP + k]
* dRedistributionMaps[ixSecondFace][ixSubElement][k];
}
}
}
// Put composed array into map
for ( size_t i = 0; i < vecAdjFaces.size(); i++ )
{
for ( size_t j = 0; j < meshThisElement.faces.size(); j++ )
{
int ixFirstFace = vecAdjFaces[i].first;
int ixSecondNode = meshThisElement.vecTargetFaceIx[j];
smatMap ( ixSecondNode, ixFirstFace ) +=
dRedistributedArray[i][j]
/ dataGLLNodalArea[ixSecondNode];
// meshThisElement.vecFaceArea[j];
// dataGLLJacobian[ixS][ixT][ixSecondElement];
}
}
// Increment the current overlap index
ixOverlap += nOverlapFaces;
//_EXCEPTION();
}
return;
}
| void moab::TempestOfflineMap::LinearRemapGLLtoGLL2_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodesIn, |
| const DataMatrix3D< double > & | dataGLLJacobianIn, | ||
| const DataMatrix3D< int > & | dataGLLNodesOut, | ||
| const DataMatrix3D< double > & | dataGLLJacobianOut, | ||
| const DataVector< 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 2531 of file TempestLinearRemap.cpp.
References ForceIntArrayConsistencyConservation(), ixOverlap, and t.
{
// Triangular quadrature rule
TriangularQuadratureRule triquadrule ( 8 );
const DataMatrix<double> & dG = triquadrule.GetG();
const DataVector<double> & dW = triquadrule.GetW();
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix<double> & smatMap = this->GetSparseMatrix();
// Sample coefficients
DataMatrix<double> dSampleCoeffIn;
dSampleCoeffIn.Initialize ( nPin, nPin );
DataMatrix<double> dSampleCoeffOut;
dSampleCoeffOut.Initialize ( nPout, nPout );
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapGLLtoGLL2_MOAB] Finite Element to Finite Element Projection" );
Announce ( "Order of the input FE polynomial interpolant: %i", nPin );
Announce ( "Order of the output FE polynomial interpolant: %i", nPout );
}
// Build the integration array for each element on m_meshOverlap
DataMatrix3D<double> dGlobalIntArray;
dGlobalIntArray.Initialize (
nPin * nPin,
m_meshOverlap->faces.size(),
nPout * nPout );
// Number of overlap Faces per source Face
DataVector<int> nAllOverlapFaces;
nAllOverlapFaces.Initialize ( 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];
}
// Geometric area of each output node
DataMatrix<double> dGeometricOutputArea;
dGeometricOutputArea.Initialize (
m_meshOutput->faces.size(), nPout * nPout );
// Area of each overlap element in the output basis
DataMatrix<double> dOverlapOutputArea;
dOverlapOutputArea.Initialize (
m_meshOverlap->faces.size(), nPout * nPout );
// Loop through all faces on m_meshInput
ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
if ( is_root )
Announce ( "Building conservative distribution maps" );
for ( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
{
#ifdef VERBOSE
// Announce computation progress
if ( ixFirst % outputFrequency == 0 && is_root )
{
Announce ( "Element %i", ixFirst );
}
#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;
int nOverlapTriangles = faceOverlap.edges.size() - 2;
// Quantities from the Second Mesh
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
const NodeVector & nodesSecond = m_meshOutput->nodes;
const Face & faceSecond = m_meshOutput->faces[ixSecond];
// Loop over all sub-triangles of this Overlap Face
for ( int j = 0; j < nOverlapTriangles; j++ )
{
// Cornerpoints of triangle
const Node & node0 = nodesOverlap[faceOverlap[0]];
const Node & node1 = nodesOverlap[faceOverlap[j + 1]];
const Node & node2 = nodesOverlap[faceOverlap[j + 2]];
// Calculate the area of the modified Face
Face faceTri ( 3 );
faceTri.SetNode ( 0, faceOverlap[0] );
faceTri.SetNode ( 1, faceOverlap[j + 1] );
faceTri.SetNode ( 2, faceOverlap[j + 2] );
double dTriArea =
CalculateFaceArea ( faceTri, nodesOverlap );
for ( int k = 0; k < triquadrule.GetPoints(); k++ )
{
double * dGL = dG[k];
// Get the nodal location of this point
double dX[3];
dX[0] = dGL[0] * node0.x + dGL[1] * node1.x + dGL[2] * node2.x;
dX[1] = dGL[0] * node0.y + dGL[1] * node1.y + dGL[2] * node2.y;
dX[2] = dGL[0] * node0.z + dGL[1] * node1.z + dGL[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
DataMatrix<double> dCoeff;
dCoeff.Initialize ( 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
DataVector<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
DataVector<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++;
}
}
}
/*
// 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);
}
_EXCEPTION();
*/
// Increment the current overlap index
ixOverlap += nOverlapFaces;
}
// Build redistribution map within target element
Announce ( "Building redistribution maps on target mesh" );
DataVector<double> dRedistSourceArea ( nPout * nPout );
DataVector<double> dRedistTargetArea ( nPout * nPout );
std::vector< DataMatrix<double> > dRedistributionMaps;
dRedistributionMaps.resize ( m_meshOutput->faces.size() );
for ( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
{
dRedistributionMaps[ixSecond].Initialize (
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
DataVector<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;
Announce ( "Assembling map" );
// Map from source DOFs to target DOFs with redistribution applied
DataMatrix<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 )
{
Announce ( "Element %i", ixFirst );
}
#endif
// Number of overlapping Faces and triangles
int nOverlapFaces = nAllOverlapFaces[ixFirst];
// Put composed array into map
for ( int j = 0; j < nOverlapFaces; j++ )
{
int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
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::TempestOfflineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodesIn, |
| const DataMatrix3D< double > & | dataGLLJacobianIn, | ||
| const DataMatrix3D< int > & | dataGLLNodesOut, | ||
| const DataMatrix3D< double > & | dataGLLJacobianOut, | ||
| const DataVector< 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 3117 of file TempestLinearRemap.cpp.
{
// Gauss-Lobatto quadrature within Faces
DataVector<double> dGL;
DataVector<double> dWL;
GaussLobattoQuadrature::GetPoints ( nPout, 0.0, 1.0, dGL, dWL );
// Utilities
MeshUtilitiesFuzzy utils;
// Get SparseMatrix represntation of the OfflineMap
SparseMatrix<double> & smatMap = this->GetSparseMatrix();
// Sample coefficients
DataMatrix<double> dSampleCoeffIn;
dSampleCoeffIn.Initialize ( nPin, nPin );
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB] Finite Element to Finite Element (Pointwise) Projection" );
Announce ( "Order of the input FE polynomial interpolant: %i", nPin );
Announce ( "Order of the output FE polynomial interpolant: %i", nPout );
}
// Number of overlap Faces per source Face
DataVector<int> nAllOverlapFaces;
nAllOverlapFaces.Initialize ( 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
DataVector<bool> fSecondNodeFound ( dataNodalAreaOut.GetRows() );
ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
// 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 )
{
Announce ( "Element %i", ixFirst );
}
#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 overlap Mesh
// const Face & faceOverlap = m_meshOverlap->faces[ixOverlap + i];
// const NodeVector & nodesOverlap = m_meshOverlap->nodes;
// size_t nOverlapTriangles = faceOverlap.edges.size() - 2;
// Quantities from the Second Mesh
int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
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;
}
| void moab::TempestOfflineMap::LinearRemapSE0_Tempest_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodes, |
| const DataMatrix3D< double > & | dataGLLJacobian | ||
| ) | [private] |
Generate the OfflineMap for linear conserative element-average spectral element to element average remapping.
| void moab::TempestOfflineMap::LinearRemapSE4_Tempest_MOAB | ( | const DataMatrix3D< int > & | dataGLLNodes, |
| const DataMatrix3D< 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 862 of file TempestLinearRemap.cpp.
References ForceConsistencyConservation3(), ixOverlap, and l.
{
// 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 DataMatrix<double> & TriQuadratureG = triquadrule.GetG();
const DataVector<double> & TriQuadratureW = triquadrule.GetW();
// Sample coefficients
DataMatrix<double> dSampleCoeff;
dSampleCoeff.Initialize ( nP, nP );
// GLL Quadrature nodes on quadrilateral elements
DataVector<double> dG;
DataVector<double> dW;
GaussLobattoQuadrature::GetPoints ( nP, 0.0, 1.0, dG, dW );
// Announcemnets
if ( is_root )
{
Announce ( "[moab::TempestOfflineMap::LinearRemapSE4_Tempest_MOAB] Finite Element to Finite Volume Projection" );
Announce ( "Triangular quadrature rule order %i", TriQuadRuleOrder );
Announce ( "Order of the FE polynomial interpolant: %i", 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
DataVector<double> vecSourceArea;
vecSourceArea.Initialize ( nP * nP );
DataVector<double> vecTargetArea;
DataMatrix<double> dCoeff;
#ifdef VERBOSE
std::stringstream sstr;
sstr << "remapdata_" << rank << ".txt";
std::ofstream output_file ( sstr.str() );
#endif
// Current Overlap Face
int ixOverlap = 0;
#ifdef VERBOSE
const unsigned outputFrequency = (m_meshInputCov->faces.size()/10);
#endif
// 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 )
{
Announce ( "Element %i/%i", 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
DataMatrix3D<double> dRemapCoeff;
dRemapCoeff.Initialize ( nP, nP, nOverlapFaces );
// Find the local remap coefficients
for ( int j = 0; j < nOverlapFaces; j++ )
{
const Face & faceOverlap = m_meshOverlap->faces[ixOverlap + j];
// #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 nOverlapTriangles = faceOverlap.edges.size() - 2;
// #define USE_MININDEX
#ifdef USE_MININDEX
// first find out the minimum node, start there the triangle decomposition
int minIndex = 0;
int nnodes = faceOverlap.edges.size();
for (int j1=1; j1<nnodes; j1++)
{
if ( nodesOverlap[faceOverlap[j1]] < nodesOverlap[faceOverlap[minIndex]] )
{
minIndex = j1;
}
}
#endif
// Loop over all sub-triangles of this Overlap Face
for ( int k = 0; k < nOverlapTriangles; k++ )
{
#ifdef USE_MININDEX
// Cornerpoints of triangle, they start at the minimal Node, for consistency
const Node & node0 = nodesOverlap[faceOverlap[minIndex]];
const Node & node1 = nodesOverlap[faceOverlap[(minIndex + k + 1)%nnodes]];
const Node & node2 = nodesOverlap[faceOverlap[(minIndex + k + 2)%nnodes]];
// Calculate the area of the modified Face
Face faceTri ( 3 );
faceTri.SetNode ( 0, faceOverlap[minIndex] );
faceTri.SetNode ( 1, faceOverlap[(minIndex + k + 1)%nnodes] );
faceTri.SetNode ( 2, faceOverlap[(minIndex + k + 2)%nnodes] );
#else
// Cornerpoints of triangle
const Node & node0 = nodesOverlap[faceOverlap[0]];
const Node & node1 = nodesOverlap[faceOverlap[k+1]];
const Node & node2 = nodesOverlap[faceOverlap[k+2]];
// Calculate the area of the modified Face
Face faceTri(3);
faceTri.SetNode(0, faceOverlap[0]);
faceTri.SetNode(1, faceOverlap[k+1]);
faceTri.SetNode(2, faceOverlap[k+2]);
#endif
double dTriangleArea =
CalculateFaceArea ( faceTri, nodesOverlap );
// 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;
double dMag = sqrt (
nodeQuadrature.x * nodeQuadrature.x
+ nodeQuadrature.y * nodeQuadrature.y
+ nodeQuadrature.z * nodeQuadrature.z );
nodeQuadrature.x /= dMag;
nodeQuadrature.y /= dMag;
nodeQuadrature.z /= dMag;
// 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 )
)
{
_EXCEPTION2 ( "Inverse Map out of range (%1.5e %1.5e)",
dAlpha, dBeta );
}
// Sample the finite element at this point
SampleGLLFiniteElement (
nMonotoneType,
nP,
dAlpha,
dBeta,
dSampleCoeff );
// Add sample coefficients to the map
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.Initialize ( nOverlapFaces );
for ( int j = 0; j < nOverlapFaces; j++ )
{
vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
}
dCoeff.Initialize ( 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.Initialize ( 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.Initialize ( 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];
if (ixSecondFace < 0) // signal to not participate, because it is a ghost target
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::TempestOfflineMap::SetDofMapAssociation | ( | DiscretizationType | srcType, |
| bool | isSrcContinuous, | ||
| DataMatrix3D< int > * | srcdataGLLNodes, | ||
| DataMatrix3D< int > * | srcdataGLLNodesSrc, | ||
| DiscretizationType | destType, | ||
| bool | isDestContinuous, | ||
| DataMatrix3D< int > * | tgtdataGLLNodes | ||
| ) | [private] |
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 158 of file TempestOfflineMap.cpp.
References 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 = mbCore->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 = mbCore->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 = mbCore->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 ldof = (*srcdataGLLNodes)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[ldof] = true;
}
output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << 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 ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[ldof] = true;
}
output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << 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 ldof = (*srcdataGLLNodes)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[ldof] = true;
}
output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << src_soln_gdofs[idof] << ", " << 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 ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[ldof] = true;
}
output_file << m_remapper->lid_to_gid_src[j] << ", " << idof << ", " << ldof << ", " << locsrc_soln_gdofs[idof] << ", " << 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
col_dofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
col_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
col_gdofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
src_soln_gdofs.resize(m_remapper->m_covering_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX);
rval = mbCore->tag_get_data ( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR(rval);
m_nTotDofs_SrcCov = 0;
if (srcdataGLLNodes == NULL || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
for (unsigned i=0; i < col_dofmap.size(); ++i) {
col_dofmap[i] = i;
col_ldofmap[i] = i;
col_gdofmap[i] = src_soln_gdofs[i];
if (vprint) std::cout << "Col: " << i << ", " << src_soln_gdofs[i] << "\n";
m_nTotDofs_SrcCov++;
}
}
else {
if (isSrcContinuous) dgll_cgll_covcol_ldofmap.resize (m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, 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 ldof = (*srcdataGLLNodes)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_covcol_ldofmap[ldof] ) {
m_nTotDofs_SrcCov++;
dgll_cgll_covcol_ldofmap[ldof] = true;
}
if ( !isSrcContinuous ) m_nTotDofs_SrcCov++;
col_dofmap[ idof ] = ldof;
col_ldofmap[ ldof ] = idof;
assert(src_soln_gdofs[idof] > 0);
col_gdofmap[ idof ] = src_soln_gdofs[idof] - 1;
if (vprint) std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << idof << ", " << ldof << ", " << col_gdofmap[idof] << ", " << m_nTotDofs_SrcCov << "\n";
}
}
}
}
srccol_dofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
srccol_ldofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
srccol_gdofmap.resize (m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, ULONG_MAX);
locsrc_soln_gdofs.resize(m_remapper->m_source_entities.size()*m_nDofsPEl_Src*m_nDofsPEl_Src, INT_MAX);
rval = mbCore->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 || srcType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
for (unsigned i=0; i < srccol_dofmap.size(); ++i) {
srccol_dofmap[i] = i;
srccol_ldofmap[i] = i;
srccol_gdofmap[i] = locsrc_soln_gdofs[i];
m_nTotDofs_Src++;
}
}
else {
if (isSrcContinuous) dgll_cgll_col_ldofmap.resize(m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, 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 ldof = (*srcdataGLLNodesSrc)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
if ( isSrcContinuous && !dgll_cgll_col_ldofmap[ldof] ) {
m_nTotDofs_Src++;
dgll_cgll_col_ldofmap[ldof] = true;
}
if ( !isSrcContinuous ) m_nTotDofs_Src++;
srccol_dofmap[ idof ] = ldof;
srccol_ldofmap[ ldof ] = idof;
srccol_gdofmap[ idof ] = locsrc_soln_gdofs[idof] - 1;
}
}
}
}
row_dofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
row_gdofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, ULONG_MAX);
tgt_soln_gdofs.resize(m_remapper->m_target_entities.size()*m_nDofsPEl_Dest*m_nDofsPEl_Dest, INT_MAX);
rval = mbCore->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
m_nTotDofs_Dest = 0;
if (tgtdataGLLNodes == NULL || destType == DiscretizationType_FV) { /* we only have a mapping for elements as DoFs */
for (unsigned i=0; i < row_dofmap.size(); ++i) {
row_dofmap[i] = i;
row_ldofmap[i] = i;
row_gdofmap[i] = tgt_soln_gdofs[i];
// if (vprint) std::cout << "Row: " << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << row_dofmap[i] << ", " << row_ldofmap[i] << ", " << tgt_soln_gdofs[i] << ", " << row_gdofmap[i] << "\n";
m_nTotDofs_Dest++;
}
}
else {
if (isTgtContinuous) dgll_cgll_row_ldofmap.resize (m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest, 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 ldof = (*tgtdataGLLNodes)[p][q][j] - 1;
const int idof = j * m_nDofsPEl_Dest * m_nDofsPEl_Dest + p * m_nDofsPEl_Dest + q;
if ( isTgtContinuous && !dgll_cgll_row_ldofmap[ldof] ) {
m_nTotDofs_Dest++;
dgll_cgll_row_ldofmap[ldof] = true;
}
if ( !isTgtContinuous ) m_nTotDofs_Dest++;
row_dofmap[ idof ] = ldof;
row_ldofmap[ ldof ] = idof;
row_gdofmap[ idof ] = tgt_soln_gdofs[idof] - 1;
if (vprint) std::cout << "Row: " << idof << ", " << ldof << ", " << tgt_soln_gdofs[idof] - 1 << "\n";
}
}
}
}
// Let us also allocate the local representation of the sparse matrix
#ifdef MOAB_HAVE_EIGEN
// if (vprint)
{
std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_dofmap.size() << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_dofmap.size() << "\n";
// std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
}
#endif
return moab::MB_SUCCESS;
}
| moab::ErrorCode moab::TempestOfflineMap::SetDofMapTags | ( | const std::string | srcDofTagName, |
| const std::string | tgtDofTagName | ||
| ) | [private] |
Store the tag names associated with global DoF ids for source and target meshes.
Definition at line 144 of file TempestOfflineMap.cpp.
References MB_CHK_ERR, MB_SUCCESS, MB_TAG_ANY, and MB_TYPE_INTEGER.
{
moab::ErrorCode rval;
rval = mbCore->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);
rval = mbCore->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);
return moab::MB_SUCCESS;
}
| std::vector<unsigned long> moab::TempestOfflineMap::col_dofmap |
Definition at line 411 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::col_gdofmap |
Definition at line 412 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::col_ldofmap |
Definition at line 413 of file TempestOfflineMap.hpp.
| DataMatrix3D<int> moab::TempestOfflineMap::dataGLLNodesDest |
Definition at line 415 of file TempestOfflineMap.hpp.
| DataMatrix3D<int> moab::TempestOfflineMap::dataGLLNodesSrc |
Definition at line 415 of file TempestOfflineMap.hpp.
| DataMatrix3D<int> moab::TempestOfflineMap::dataGLLNodesSrcCov |
Definition at line 415 of file TempestOfflineMap.hpp.
Definition at line 425 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
Definition at line 425 of file TempestOfflineMap.hpp.
Referenced by LinearRemapFVtoFV_Tempest_MOAB(), and TempestOfflineMap().
Definition at line 416 of file TempestOfflineMap.hpp.
Definition at line 410 of file TempestOfflineMap.hpp.
The original tag data and local to global DoF mapping to associate matrix values to solution.
Definition at line 410 of file TempestOfflineMap.hpp.
The boolean flag representing whether the root process has the updated global view.
Definition at line 375 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
Definition at line 420 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
Definition at line 421 of file TempestOfflineMap.hpp.
Referenced by LinearRemapFVtoFV_Tempest_MOAB(), and TempestOfflineMap().
Definition at line 422 of file TempestOfflineMap.hpp.
Referenced by LinearRemapFVtoFV_Tempest_MOAB(), and TempestOfflineMap().
Definition at line 423 of file TempestOfflineMap.hpp.
Referenced by LinearRemapFVtoFV_Tempest_MOAB(), and TempestOfflineMap().
Definition at line 418 of file TempestOfflineMap.hpp.
Definition at line 418 of file TempestOfflineMap.hpp.
Definition at line 417 of file TempestOfflineMap.hpp.
Definition at line 417 of file TempestOfflineMap.hpp.
Definition at line 417 of file TempestOfflineMap.hpp.
Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections.
The fundamental remapping operator object.
Definition at line 364 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
Definition at line 416 of file TempestOfflineMap.hpp.
| OfflineMap* moab::TempestOfflineMap::m_weightMapGlobal |
The SparseMatrix representing this operator.
Definition at line 370 of file TempestOfflineMap.hpp.
Referenced by GetGlobalSourceAreas(), GetGlobalTargetAreas(), and TempestOfflineMap().
The DataVector that stores the global (GID-based) areas of the source mesh.
The DataVector that stores the global (GID-based) areas of the target mesh. The reference to the moab::Core object that contains source/target and overlap sets.
Definition at line 398 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
Definition at line 426 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
| std::vector<unsigned long> moab::TempestOfflineMap::row_dofmap |
Definition at line 411 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::row_gdofmap |
Definition at line 412 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::row_ldofmap |
Definition at line 413 of file TempestOfflineMap.hpp.
Definition at line 426 of file TempestOfflineMap.hpp.
Referenced by TempestOfflineMap().
| std::vector<unsigned long> moab::TempestOfflineMap::srccol_dofmap |
Definition at line 411 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::srccol_gdofmap |
Definition at line 412 of file TempestOfflineMap.hpp.
| std::vector<unsigned long> moab::TempestOfflineMap::srccol_ldofmap |
Definition at line 413 of file TempestOfflineMap.hpp.
1.7.6.1