moab::TempestOnlineMap Class Reference

An offline map between two Meshes. More...

#include <TempestOnlineMap.hpp>

Inheritance diagram for moab::TempestOnlineMap:

Collaboration diagram for moab::TempestOnlineMap:

Public Types
enum	DiscretizationType { DiscretizationType_FV , DiscretizationType_CGLL , DiscretizationType_DGLL , DiscretizationType_PCLOUD }
typedef double(*	sample_function) (double, double)

Public Member Functions
	TempestOnlineMap (moab::TempestRemapper *remapper)
	Generate the metadata associated with the offline map. More...
virtual	~TempestOnlineMap ()
	Define a virtual destructor. More...
moab::ErrorCode	GenerateRemappingWeights (std::string strInputType, std::string strOutputType, const GenerateOfflineMapAlgorithmOptions &mapOptions, const std::string &srcDofTagName="GLOBAL_ID", const std::string &tgtDofTagName="GLOBAL_ID")
	Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix. More...
moab::ErrorCode	ReadParallelMap (const char *strSource, const std::vector< int > &owned_dof_ids, bool row_major_ownership=true)
	Generate the metadata associated with the offline map. More...
moab::ErrorCode	WriteParallelMap (const std::string &strTarget, const std::map< std::string, std::string > &attrMap)
	Write the TempestOnlineMap to a parallel NetCDF file. More...
virtual int	IsConsistent (double dTolerance)
	Determine if the map is first-order accurate. More...
virtual int	IsConservative (double dTolerance)
	Determine if the map is conservative. More...
virtual int	IsMonotone (double dTolerance)
	Determine if the map is monotone. More...
const DataArray1D< double > &	GetGlobalSourceAreas () const
	If we computed the reduction, get the vector representing the source areas for all entities in the mesh More...
const DataArray1D< double > &	GetGlobalTargetAreas () const
	If we computed the reduction, get the vector representing the target areas for all entities in the mesh More...
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 More...
moab::ErrorCode	SetDOFmapAssociation (DiscretizationType srcType, bool isSrcContinuous, DataArray3D< int > *srcdataGLLNodes, DataArray3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, bool isDestContinuous, DataArray3D< int > *tgtdataGLLNodes)
	Compute the association between the solution tag global DoF numbering and the local matrix numbering so that matvec operations can be performed consistently. More...
int	GetSourceGlobalNDofs ()
	Get the number of total Degrees-Of-Freedom defined on the source mesh. More...
int	GetDestinationGlobalNDofs ()
	Get the number of total Degrees-Of-Freedom defined on the destination mesh. More...
int	GetSourceLocalNDofs ()
	Get the number of local Degrees-Of-Freedom defined on the source mesh. More...
int	GetDestinationLocalNDofs ()
	Get the number of local Degrees-Of-Freedom defined on the destination mesh. More...
int	GetSourceNDofsPerElement ()
	Get the number of Degrees-Of-Freedom per element on the source mesh. More...
int	GetDestinationNDofsPerElement ()
	Get the number of Degrees-Of-Freedom per element on the destination mesh. More...
void	SetSourceNDofsPerElement (int ns)
	Set the number of Degrees-Of-Freedom per element on the source mesh. More...
void	SetDestinationNDofsPerElement (int nt)
	Get the number of Degrees-Of-Freedom per element on the destination mesh. More...
int	GetRowGlobalDoF (int localID) const
	Get the global Degrees-Of-Freedom ID on the destination mesh. More...
int	GetIndexOfRowGlobalDoF (int globalRowDoF) const
	Get the index of globaRowDoF. More...
int	GetColGlobalDoF (int localID) const
	Get the global Degrees-Of-Freedom ID on the source mesh. More...
int	GetIndexOfColGlobalDoF (int globalColDoF) const
	Get the index of globaColDoF. More...
moab::ErrorCode	ApplyWeights (std::vector< double > &srcVals, std::vector< double > &tgtVals, bool transpose=false)
	Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More...
moab::ErrorCode	ApplyWeights (moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose=false)
	Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More...
moab::ErrorCode	DefineAnalyticalSolution (moab::Tag &exactSolnTag, const std::string &solnName, Remapper::IntersectionContext ctx, sample_function testFunction, moab::Tag *clonedSolnTag=NULL, std::string cloneSolnName="")
	Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user. More...
moab::ErrorCode	ComputeMetrics (Remapper::IntersectionContext ctx, moab::Tag &exactTag, moab::Tag &approxTag, std::map< std::string, double > &metrics, bool verbose=true)
	Compute the error between a sampled (exact) solution and a projected solution in various error norms. More...
moab::ErrorCode	fill_row_ids (std::vector< int > &ids_of_interest)
moab::ErrorCode	fill_col_ids (std::vector< int > &ids_of_interest)
moab::ErrorCode	set_col_dc_dofs (std::vector< int > &values_entities)
moab::ErrorCode	set_row_dc_dofs (std::vector< int > &values_entities)

Private Member Functions
moab::ErrorCode	LinearRemapNN_MOAB (bool use_GID_matching=false, bool strict_check=false)
	Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh. More...
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. More...
void	LinearRemapSE0_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian)
	Generate the OfflineMap for linear conserative element-average spectral element to element average remapping. More...
void	LinearRemapSE4_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, int nMonotoneType, bool fContinuousIn, bool fNoConservation)
	Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping. More...
void	LinearRemapFVtoGLL_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, const DataArray1D< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation)
	Generate the OfflineMap for remapping from finite volumes to finite elements. More...
void	LinearRemapGLLtoGLL2_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut, bool fNoConservation)
	Generate the OfflineMap for remapping from finite elements to finite elements. More...
void	LinearRemapGLLtoGLL2_Pointwise_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut)
	Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation). More...
moab::ErrorCode	WriteSCRIPMapFile (const std::string &strOutputFile, const std::map< std::string, std::string > &attrMap)
	Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections. More...
moab::ErrorCode	WriteHDF5MapFile (const std::string &filename)
	Parallel I/O with NetCDF to write out the SCRIP file from multiple processors. More...
void	setup_sizes_dimensions ()

Private Attributes
moab::TempestRemapper *	m_remapper
	The fundamental remapping operator object. More...
moab::Interface *	m_interface
	The reference to the moab::Core object that contains source/target and overlap sets. More...
moab::Tag	m_dofTagSrc
	The original tag data and local to global DoF mapping to associate matrix values to solution More...
moab::Tag	m_dofTagDest
std::vector< unsigned >	row_gdofmap
std::vector< unsigned >	col_gdofmap
std::vector< unsigned >	srccol_gdofmap
std::vector< int >	row_dtoc_dofmap
std::vector< int >	col_dtoc_dofmap
std::vector< int >	srccol_dtoc_dofmap
std::map< int, int >	rowMap
std::map< int, int >	colMap
DataArray3D< int >	dataGLLNodesSrc
DataArray3D< int >	dataGLLNodesSrcCov
DataArray3D< int >	dataGLLNodesDest
DiscretizationType	m_srcDiscType
DiscretizationType	m_destDiscType
int	m_nTotDofs_Src
int	m_nTotDofs_SrcCov
int	m_nTotDofs_Dest
int	m_nDofsPEl_Src
int	m_nDofsPEl_Dest
DiscretizationType	m_eInputType
DiscretizationType	m_eOutputType
bool	m_bConserved
int	m_iMonotonicity
Mesh *	m_meshInput
Mesh *	m_meshInputCov
Mesh *	m_meshOutput
Mesh *	m_meshOverlap
bool	is_parallel
bool	is_root
int	rank
int	size

Detailed Description

An offline map between two Meshes.

Definition at line 48 of file TempestOnlineMap.hpp.

Member Typedef Documentation

sample_function

typedef double( * moab::TempestOnlineMap::sample_function) (double, double)

Definition at line 349 of file TempestOnlineMap.hpp.

Member Enumeration Documentation

DiscretizationType

enum moab::TempestOnlineMap::DiscretizationType

Enumerator
DiscretizationType_FV
DiscretizationType_CGLL
DiscretizationType_DGLL
DiscretizationType_PCLOUD

Definition at line 64 of file TempestOnlineMap.hpp.

 {
 DiscretizationType_FV,
 DiscretizationType_CGLL,
 DiscretizationType_DGLL,
 DiscretizationType_PCLOUD
 };

Constructor & Destructor Documentation

TempestOnlineMap()

moab::TempestOnlineMap::TempestOnlineMap

(

moab::TempestRemapper *

remapper

)

Generate the metadata associated with the offline map.

Definition at line 55 of file TempestOnlineMap.cpp.

 : OfflineMap(), m_remapper( remapper )
{
 // Get the references for the MOAB core objects
 m_interface = m_remapper->get_interface();
#ifdef MOAB_HAVE_MPI
 m_pcomm = m_remapper->get_parallel_communicator();
#endif
 
 // Update the references to the meshes
 m_meshInput = remapper->GetMesh( moab::Remapper::SourceMesh );
 m_meshInputCov = remapper->GetCoveringMesh();
 m_meshOutput = remapper->GetMesh( moab::Remapper::TargetMesh );
 m_meshOverlap = remapper->GetMesh( moab::Remapper::OverlapMesh );
 
 is_parallel = remapper->is_parallel;
 is_root = remapper->is_root;
 rank = remapper->rank;
 size = remapper->size;
 
 // Initialize dimension information from file
 this->setup_sizes_dimensions();
}

References moab::Remapper::get_interface(), moab::TempestRemapper::GetCoveringMesh(), moab::TempestRemapper::GetMesh(), is_parallel, moab::TempestRemapper::is_parallel, is_root, moab::TempestRemapper::is_root, m_interface, m_meshInput, m_meshInputCov, m_meshOutput, m_meshOverlap, m_remapper, moab::Remapper::OverlapMesh, rank, moab::TempestRemapper::rank, setup_sizes_dimensions(), size, moab::TempestRemapper::size, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.

~TempestOnlineMap()

moab::TempestOnlineMap::~TempestOnlineMap ( )

virtual

Define a virtual destructor.

Definition at line 103 of file TempestOnlineMap.cpp.

{
 m_interface = NULL;
#ifdef MOAB_HAVE_MPI
 m_pcomm = NULL;
#endif
 m_meshInput = NULL;
 m_meshOutput = NULL;
 m_meshOverlap = NULL;
}

Member Function Documentation

ApplyWeights() [1/2]

moab::ErrorCode moab::TempestOnlineMap::ApplyWeights	(	moab::Tag	srcSolutionTag,
		moab::Tag	tgtSolutionTag,
		bool	transpose = false
	)

Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals

Definition at line 1676 of file TempestOnlineMap.cpp.

{
 moab::ErrorCode rval;
 
 std::vector< double > solSTagVals;
 std::vector< double > solTTagVals;
 
 moab::Range sents, tents;
 if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
 {
 if( m_remapper->point_cloud_source )
 {
 moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
 solSTagVals.resize( covSrcEnts.size(), -1.0 );
 sents = covSrcEnts;
 }
 else
 {
 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
 -1.0 );
 sents = covSrcEnts;
 }
 if( m_remapper->point_cloud_target )
 {
 moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
 solTTagVals.resize( tgtEnts.size(), -1.0 );
 tents = tgtEnts;
 }
 else
 {
 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
 solTTagVals.resize(
 tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
 tents = tgtEnts;
 }
 }
 else
 {
 moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
 moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
 solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
 -1.0 );
 solTTagVals.resize(
 tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
 
 sents = covSrcEnts;
 tents = tgtEnts;
 }
 
 // The tag data is np*np*n_el_src
 rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );
 
 // Compute the application of weights on the suorce solution data and store it in the
 // destination solution vector data Optionally, can also perform the transpose application of
 // the weight matrix. Set the 3rd argument to true if this is needed
 rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );
 
 // The tag data is np*np*n_el_dest
 rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 return moab::MB_SUCCESS;
}

References moab::Remapper::CoveringMesh, ErrorCode, MB_CHK_SET_ERR, MB_SUCCESS, moab::Range::size(), and moab::Remapper::TargetMesh.

ApplyWeights() [2/2]

moab::ErrorCode moab::TempestOnlineMap::ApplyWeights	(	std::vector< double > &	srcVals,
		std::vector< double > &	tgtVals,
		bool	transpose = false
	)

Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals

Referenced by main().

ComputeMetrics()

moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics	(	Remapper::IntersectionContext	ctx,
		moab::Tag &	exactTag,
		moab::Tag &	approxTag,
		std::map< std::string, double > &	metrics,
		bool	verbose = true
	)

Compute the error between a sampled (exact) solution and a projected solution in various error norms.

Definition at line 2104 of file TempestOnlineMap.cpp.

{
 moab::ErrorCode rval;
 const bool outputEnabled = ( is_root );
 int discOrder;
 DiscretizationType discMethod;
 moab::EntityHandle meshset;
 moab::Range entities;
 Mesh* trmesh;
 switch( ctx )
 {
 case Remapper::SourceMesh:
 meshset = m_remapper->m_covering_source_set;
 trmesh = m_remapper->m_covering_source;
 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
 : m_remapper->m_covering_source_entities );
 discOrder = m_nDofsPEl_Src;
 discMethod = m_eInputType;
 break;
 
 case Remapper::TargetMesh:
 meshset = m_remapper->m_target_set;
 trmesh = m_remapper->m_target;
 entities =
 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
 discOrder = m_nDofsPEl_Dest;
 discMethod = m_eOutputType;
 break;
 
 default:
 if( outputEnabled )
 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
 return moab::MB_FAILURE;
 }
 
 // Let us create teh solution tag with appropriate information for name, discretization order
 // (DoF space)
 std::string exactTagName, projTagName;
 const int ntotsize = entities.size() * discOrder * discOrder;
 int ntotsize_glob = 0;
 std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
 rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
 rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
 rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
 rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
 
 std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 ); // L1Err, L2Err, LinfErr
 for( int i = 0; i < ntotsize; ++i )
 {
 const double error = fabs( exactSolution[i] - projSolution[i] );
 errnorms[0] += error;
 errnorms[1] += error * error;
 errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] );
 }
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 {
 MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
 MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
 MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
 }
#else
 ntotsize_glob = ntotsize;
 globerrnorms = errnorms;
#endif
 globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob );
 globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob );
 
 metrics.clear();
 metrics["L1Error"] = globerrnorms[0];
 metrics["L2Error"] = globerrnorms[1];
 metrics["LinfError"] = globerrnorms[2];
 
 if( verbose && is_root )
 {
 std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
 std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
 std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
 std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl;
 }
 
 return moab::MB_SUCCESS;
}

References entities, moab::error(), ErrorCode, MB_CHK_ERR, MB_SUCCESS, moab::Remapper::SourceMesh, moab::Remapper::TargetMesh, and verbose.

Referenced by main().

DefineAnalyticalSolution()

moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution	(	moab::Tag &	exactSolnTag,
		const std::string &	solnName,
		Remapper::IntersectionContext	ctx,
		sample_function	testFunction,
		moab::Tag *	clonedSolnTag = NULL,
		std::string	cloneSolnName = ""
	)

Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user.

Definition at line 1742 of file TempestOnlineMap.cpp.

{
 moab::ErrorCode rval;
 const bool outputEnabled = ( is_root );
 int discOrder;
 DiscretizationType discMethod;
 moab::EntityHandle meshset;
 moab::Range entities;
 Mesh* trmesh;
 switch( ctx )
 {
 case Remapper::SourceMesh:
 meshset = m_remapper->m_covering_source_set;
 trmesh = m_remapper->m_covering_source;
 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
 : m_remapper->m_covering_source_entities );
 discOrder = m_nDofsPEl_Src;
 discMethod = m_eInputType;
 break;
 
 case Remapper::TargetMesh:
 meshset = m_remapper->m_target_set;
 trmesh = m_remapper->m_target;
 entities =
 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
 discOrder = m_nDofsPEl_Dest;
 discMethod = m_eOutputType;
 break;
 
 default:
 if( outputEnabled )
 std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
 return moab::MB_FAILURE;
 }
 
 // Let us create teh solution tag with appropriate information for name, discretization order
 // (DoF space)
 rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
 MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 if( clonedSolnTag != NULL )
 {
 if( cloneSolnName.size() == 0 )
 {
 cloneSolnName = solnName + std::string( "Cloned" );
 }
 rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
 *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
 }
 
 // Triangular quadrature rule
 const int TriQuadratureOrder = 10;
 
 if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
 
 TriangularQuadratureRule triquadrule( TriQuadratureOrder );
 
 const int TriQuadraturePoints = triquadrule.GetPoints();
 
 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
 
 // Output data
 DataArray1D< double > dVar;
 DataArray1D< double > dVarMB; // re-arranged local MOAB vector
 
 // Nodal geometric area
 DataArray1D< double > dNodeArea;
 
 // Calculate element areas
 // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
 
 if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
 {
 /* Get the spectral points and sample the functionals accordingly */
 const bool fGLL = true;
 const bool fGLLIntegrate = false;
 
 // Generate grid metadata
 DataArray3D< int > dataGLLNodes;
 DataArray3D< double > dataGLLJacobian;
 
 GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
 
 // Number of elements
 int nElements = trmesh->faces.size();
 
 // Verify all elements are quadrilaterals
 for( int k = 0; k < nElements; k++ )
 {
 const Face& face = trmesh->faces[k];
 
 if( face.edges.size() != 4 )
 {
 _EXCEPTIONT( "Non-quadrilateral face detected; "
 "incompatible with --gll" );
 }
 }
 
 // Number of unique nodes
 int iMaxNode = 0;
 for( int i = 0; i < discOrder; i++ )
 {
 for( int j = 0; j < discOrder; j++ )
 {
 for( int k = 0; k < nElements; k++ )
 {
 if( dataGLLNodes[i][j][k] > iMaxNode )
 {
 iMaxNode = dataGLLNodes[i][j][k];
 }
 }
 }
 }
 
 // Get Gauss-Lobatto quadrature nodes
 DataArray1D< double > dG;
 DataArray1D< double > dW;
 
 GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
 
 // Get Gauss quadrature nodes
 const int nGaussP = 10;
 
 DataArray1D< double > dGaussG;
 DataArray1D< double > dGaussW;
 
 GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
 
 // Allocate data
 dVar.Allocate( iMaxNode );
 dVarMB.Allocate( discOrder * discOrder * nElements );
 dNodeArea.Allocate( iMaxNode );
 
 // Sample data
 for( int k = 0; k < nElements; k++ )
 {
 const Face& face = trmesh->faces[k];
 
 // Sample data at GLL nodes
 if( fGLL )
 {
 for( int i = 0; i < discOrder; i++ )
 {
 for( int j = 0; j < discOrder; j++ )
 {
 
 // Apply local map
 Node node;
 Node dDx1G;
 Node dDx2G;
 
 ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
 
 // Sample data at this point
 double dNodeLon = atan2( node.y, node.x );
 if( dNodeLon < 0.0 )
 {
 dNodeLon += 2.0 * M_PI;
 }
 double dNodeLat = asin( node.z );
 
 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
 
 dVar[dataGLLNodes[j][i][k] - 1] = dSample;
 }
 }
 // High-order Gaussian integration over basis function
 }
 else
 {
 DataArray2D< double > dCoeff( discOrder, discOrder );
 
 for( int p = 0; p < nGaussP; p++ )
 {
 for( int q = 0; q < nGaussP; q++ )
 {
 
 // Apply local map
 Node node;
 Node dDx1G;
 Node dDx2G;
 
 ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
 
 // Cross product gives local Jacobian
 Node nodeCross = CrossProduct( dDx1G, dDx2G );
 
 double dJacobian =
 sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
 
 // Find components of quadrature point in basis
 // of the first Face
 SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
 
 // Sample data at this point
 double dNodeLon = atan2( node.y, node.x );
 if( dNodeLon < 0.0 )
 {
 dNodeLon += 2.0 * M_PI;
 }
 double dNodeLat = asin( node.z );
 
 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
 
 // Integrate
 for( int i = 0; i < discOrder; i++ )
 {
 for( int j = 0; j < discOrder; j++ )
 {
 
 double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
 
 dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
 
 dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
 }
 }
 }
 }
 }
 }
 
 // Divide by area
 if( fGLLIntegrate )
 {
 for( size_t i = 0; i < dVar.GetRows(); i++ )
 {
 dVar[i] /= dNodeArea[i];
 }
 }
 
 // Let us rearrange the data based on DoF ID specification
 if( ctx == Remapper::SourceMesh )
 {
 for( unsigned j = 0; j < entities.size(); j++ )
 for( int p = 0; p < discOrder; p++ )
 for( int q = 0; q < discOrder; q++ )
 {
 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
 dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
 }
 }
 else
 {
 for( unsigned j = 0; j < entities.size(); j++ )
 for( int p = 0; p < discOrder; p++ )
 for( int q = 0; q < discOrder; q++ )
 {
 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
 dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
 }
 }
 
 // Set the tag data
 rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
 }
 else
 {
 // assert( discOrder == 1 );
 if( discMethod == DiscretizationType_FV )
 {
 /* Compute an element-wise integral to store the sampled solution based on Quadrature
 * rules */
 // Resize the array
 dVar.Allocate( trmesh->faces.size() );
 
 std::vector< Node >& nodes = trmesh->nodes;
 
 // Loop through all Faces
 for( size_t i = 0; i < trmesh->faces.size(); i++ )
 {
 const Face& face = trmesh->faces[i];
 
 // Loop through all sub-triangles
 for( size_t j = 0; j < face.edges.size() - 2; j++ )
 {
 
 const Node& node0 = nodes[face[0]];
 const Node& node1 = nodes[face[j + 1]];
 const Node& node2 = nodes[face[j + 2]];
 
 // Triangle area
 Face faceTri( 3 );
 faceTri.SetNode( 0, face[0] );
 faceTri.SetNode( 1, face[j + 1] );
 faceTri.SetNode( 2, face[j + 2] );
 
 double dTriangleArea = CalculateFaceArea( faceTri, nodes );
 
 // Calculate the element average
 double dTotalSample = 0.0;
 
 // Loop through all quadrature points
 for( int k = 0; k < TriQuadraturePoints; k++ )
 {
 Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
 TriQuadratureG[k][2] * node2.x,
 TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
 TriQuadratureG[k][2] * node2.y,
 TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
 TriQuadratureG[k][2] * node2.z );
 
 double dMagnitude = node.Magnitude();
 node.x /= dMagnitude;
 node.y /= dMagnitude;
 node.z /= dMagnitude;
 
 double dLon = atan2( node.y, node.x );
 if( dLon < 0.0 )
 {
 dLon += 2.0 * M_PI;
 }
 double dLat = asin( node.z );
 
 double dSample = ( *testFunction )( dLon, dLat );
 
 dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
 }
 
 dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
 }
 }
 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
 }
 else /* discMethod == DiscretizationType_PCLOUD */
 {
 /* Get the coordinates of the vertices and sample the functionals accordingly */
 std::vector< Node >& nodes = trmesh->nodes;
 
 // Resize the array
 dVar.Allocate( nodes.size() );
 
 for( size_t j = 0; j < nodes.size(); j++ )
 {
 Node& node = nodes[j];
 double dMagnitude = node.Magnitude();
 node.x /= dMagnitude;
 node.y /= dMagnitude;
 node.z /= dMagnitude;
 double dLon = atan2( node.y, node.x );
 if( dLon < 0.0 )
 {
 dLon += 2.0 * M_PI;
 }
 double dLat = asin( node.z );
 
 double dSample = ( *testFunction )( dLon, dLat );
 dVar[j] = dSample;
 }
 
 rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
 }
 }
 
 return moab::MB_SUCCESS;
}

References entities, ErrorCode, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, moab::Remapper::SourceMesh, and moab::Remapper::TargetMesh.

Referenced by main().

fill_col_ids()

moab::ErrorCode moab::TempestOnlineMap::fill_col_ids ( std::vector< int > &  ids_of_interest )

inline

Definition at line 384 of file TempestOnlineMap.hpp.

 {
 ids_of_interest.reserve( col_gdofmap.size() );
 // need to add 1
 for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
 ids_of_interest.push_back( *it + 1 );
 return moab::MB_SUCCESS;
 }

References col_gdofmap, and MB_SUCCESS.

fill_row_ids()

moab::ErrorCode moab::TempestOnlineMap::fill_row_ids ( std::vector< int > &  ids_of_interest )

inline

Definition at line 374 of file TempestOnlineMap.hpp.

 {
 ids_of_interest.reserve( row_gdofmap.size() );
 // need to add 1
 for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
 ids_of_interest.push_back( *it + 1 );
 
 return moab::MB_SUCCESS;
 }

References MB_SUCCESS, and row_gdofmap.

GenerateRemappingWeights()

moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights	(	std::string	strInputType,
		std::string	strOutputType,
		const GenerateOfflineMapAlgorithmOptions &	mapOptions,
		const std::string &	srcDofTagName = "GLOBAL_ID",
		const std::string &	tgtDofTagName = "GLOBAL_ID"
	)

Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix.

the tag should be created already in the e3sm workflow; if not, create it here

Definition at line 691 of file TempestOnlineMap.cpp.

{
 NcError error( NcError::silent_nonfatal );
 
 moab::DebugOutput dbgprint( std::cout, rank, 0 );
 dbgprint.set_prefix( "[TempestOnlineMap]: " );
 moab::ErrorCode rval;
 
 const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
 const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
 const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
 
 // 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,
 // for FV->FV: the rows represented target DoFs and cols represent source DoFs
 try
 {
 // Check command line parameters (data type arguments)
 STLStringHelper::ToLower( strInputType );
 STLStringHelper::ToLower( strOutputType );
 
 DiscretizationType eInputType;
 DiscretizationType eOutputType;
 
 if( strInputType == "fv" )
 {
 eInputType = DiscretizationType_FV;
 }
 else if( strInputType == "cgll" )
 {
 eInputType = DiscretizationType_CGLL;
 }
 else if( strInputType == "dgll" )
 {
 eInputType = DiscretizationType_DGLL;
 }
 else if( strInputType == "pcloud" )
 {
 eInputType = DiscretizationType_PCLOUD;
 }
 else
 {
 _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
 }
 
 if( strOutputType == "fv" )
 {
 eOutputType = DiscretizationType_FV;
 }
 else if( strOutputType == "cgll" )
 {
 eOutputType = DiscretizationType_CGLL;
 }
 else if( strOutputType == "dgll" )
 {
 eOutputType = DiscretizationType_DGLL;
 }
 else if( strOutputType == "pcloud" )
 {
 eOutputType = DiscretizationType_PCLOUD;
 }
 else
 {
 _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
 }
 
 // set all required input params
 m_bConserved = !mapOptions.fNoConservation;
 m_eInputType = eInputType;
 m_eOutputType = eOutputType;
 
 // Method flags
 std::string strMapAlgorithm( "" );
 int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
 
 // Make an index of method arguments
 std::set< std::string > setMethodStrings;
 {
 int iLast = 0;
 for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
 {
 if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
 {
 std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
 STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
 if( strMethodString.length() > 0 )
 {
 setMethodStrings.insert( strMethodString );
 }
 iLast = i + 1;
 }
 }
 }
 
 for( auto it : setMethodStrings )
 {
 // Piecewise constant monotonicity
 if( it == "mono2" )
 {
 if( nMonotoneType != 0 )
 {
 _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
 }
 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
 }
 nMonotoneType = 2;
 
 // Piecewise linear monotonicity
 }
 else if( it == "mono3" )
 {
 if( nMonotoneType != 0 )
 {
 _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
 }
 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
 }
 nMonotoneType = 3;
 
 // Volumetric remapping from FV to GLL
 }
 else if( it == "volumetric" )
 {
 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
 }
 strMapAlgorithm = "volumetric";
 
 // Inverse distance mapping
 }
 else if( it == "invdist" )
 {
 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
 }
 strMapAlgorithm = "invdist";
 
 // Delaunay triangulation mapping
 }
 else if( it == "delaunay" )
 {
 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
 }
 strMapAlgorithm = "delaunay";
 
 // Bilinear
 }
 else if( it == "bilin" )
 {
 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
 }
 strMapAlgorithm = "fvbilin";
 
 // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
 }
 else if( it == "intbilin" )
 {
 if( m_eOutputType != DiscretizationType_FV )
 {
 _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
 }
 if( m_eInputType == DiscretizationType_FV )
 {
 strMapAlgorithm = "fvintbilin";
 }
 else
 {
 strMapAlgorithm = "mono3";
 }
 
 // Integrated bilinear with generalized Barycentric coordinates
 }
 else if( it == "intbilingb" )
 {
 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
 {
 _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
 }
 strMapAlgorithm = "fvintbilingb";
 }
 else
 {
 _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
 }
 }
 
 m_nDofsPEl_Src =
 ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
 : mapOptions.nPin );
 m_nDofsPEl_Dest =
 ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
 : mapOptions.nPout );
 
 rval = SetDOFmapTags( srcDofTagName, tgtDofTagName );MB_CHK_ERR( rval );
 
 /// the tag should be created already in the e3sm workflow; if not, create it here
 Tag areaTag;
 rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
 MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT );
 if( MB_ALREADY_ALLOCATED == rval )
 {
 if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
 }
 
 double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0;
 if( !m_bPointCloudSource )
 {
 // Calculate Input Mesh Face areas
 if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
 double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
 rval = m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea );MB_CHK_ERR( rval );
 dTotalAreaInput = dTotalAreaInput_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
 
 // Input mesh areas
 m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
 // we do not need to set the area on coverage mesh, only on source and target meshes
 // rval = m_interface->tag_set_data( areaTag, m_remapper->m_covering_source_entities, m_meshInputCov->vecFaceArea );MB_CHK_ERR( rval );
 }
 
 if( !m_bPointCloudTarget )
 {
 // Calculate Output Mesh Face areas
 if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
 double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
 dTotalAreaOutput = dTotalAreaOutput_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
 rval = m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea );MB_CHK_ERR( rval );
 }
 
 if( !m_bPointCloud )
 {
 // Verify that overlap mesh is in the correct order, only if size == 1
 if( 1 == size )
 {
 int ixSourceFaceMax = ( -1 );
 int ixTargetFaceMax = ( -1 );
 
 if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() )
 {
 _EXCEPTIONT( "Invalid overlap mesh:\n"
 " Possible mesh file corruption?" );
 }
 
 for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ )
 {
 if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
 ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;
 
 if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
 ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
 }
 
 // Check for forward correspondence in overlap mesh
 if( m_meshInputCov->faces.size() - ixSourceFaceMax == 0 )
 {
 if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" );
 }
 else if( m_meshOutput->faces.size() - ixSourceFaceMax == 0 )
 { // Check for reverse correspondence in overlap mesh
 if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" );
 
 // Reorder overlap mesh
 m_meshOverlap->ExchangeFirstAndSecondMesh();
 }
 else
 { // No correspondence found
 _EXCEPTION4( "Invalid overlap mesh:\n No correspondence found with input and output meshes "
 "(%i,%i) vs (%i,%i)",
 m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax,
 ixTargetFaceMax );
 }
 }
 
 // Calculate Face areas
 if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
 double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false );
 double dTotalAreaOverlap = dTotalAreaOverlap_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
 
 // Correct areas to match the areas calculated in the overlap mesh
 // if (fCorrectAreas) // In MOAB-TempestRemap, we will always keep this to be true
 {
 if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
 DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
 DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
 
 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
 assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
 assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
 
 assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
 assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
 
 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
 {
 if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
 continue; // skip this cell since it is ghosted
 
 // let us recompute the source/target areas based on overlap mesh areas
 assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
 dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
 assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
 dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
 }
 
 for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
 {
 if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
 {
 m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
 }
 }
 for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
 {
 if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
 {
 m_meshOutput->vecFaceArea[i] = dTargetArea[i];
 }
 }
 }
 
 // Set source mesh areas in map
 if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
 {
 this->SetSourceAreas( m_meshInputCov->vecFaceArea );
 if( m_meshInputCov->vecMask.size() )
 {
 this->SetSourceMask( m_meshInputCov->vecMask );
 }
 }
 
 // Set target mesh areas in map
 if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
 {
 this->SetTargetAreas( m_meshOutput->vecFaceArea );
 if( m_meshOutput->vecMask.size() )
 {
 this->SetTargetMask( m_meshOutput->vecMask );
 }
 }
 
 /*
 // Recalculate input mesh area from overlap mesh
 if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
 dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
 dbgprint.printf(0, "Recalculating source mesh areas\n");
 dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
 dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
 }
 */
 }
 
 // Finite volume input / Finite volume output
 if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
 {
 // Generate reverse node array and edge map
 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
 
 // Initialize coordinates for map
 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
 
 // Finite volume input / Finite element output
 rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, false, NULL );MB_CHK_ERR( rval );
 
 // Construct remap for FV-FV
 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
 
 // Construct OfflineMap
 if( strMapAlgorithm == "invdist" )
 {
 if( is_root ) AnnounceStartBlock( "Calculating map (invdist)" );
 LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
 }
 else if( strMapAlgorithm == "delaunay" )
 {
 if( is_root ) AnnounceStartBlock( "Calculating map (delaunay)" );
 LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
 }
 else if( strMapAlgorithm == "fvintbilin" )
 {
 if( is_root ) AnnounceStartBlock( "Calculating map (intbilin)" );
 LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
 }
 else if( strMapAlgorithm == "fvintbilingb" )
 {
 if( is_root ) AnnounceStartBlock( "Calculating map (intbilingb)" );
 LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
 }
 else if( strMapAlgorithm == "fvbilin" )
 {
#ifdef VERBOSE
 if( is_root )
 {
 m_meshInputCov->Write( "SourceMeshMBTR.g" );
 m_meshOutput->Write( "TargetMeshMBTR.g" );
 }
 else
 {
 m_meshInputCov->Write( "SourceMeshMBTR" + std::to_string( rank ) + ".g" );
 m_meshOutput->Write( "TargetMeshMBTR" + std::to_string( rank ) + ".g" );
 }
#endif
 if( is_root ) AnnounceStartBlock( "Calculating map (bilin)" );
 LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
 }
 else
 {
 if( is_root ) AnnounceStartBlock( "Calculating map (default)" );
 // LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
 // ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
 LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
 }
 }
 else if( eInputType == DiscretizationType_FV )
 {
 DataArray3D< double > dataGLLJacobian;
 
 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
 double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
 dataGLLNodesDest, dataGLLJacobian );
 
 double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
 
 // Initialize coordinates for map
 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
 
 // Generate the continuous Jacobian
 bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
 
 if( eOutputType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
 }
 
 // Generate reverse node array and edge map
 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
 
 // Finite volume input / Finite element output
 rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType,
 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
 
 // Generate remap weights
 if( strMapAlgorithm == "volumetric" )
 {
 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
 LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
 dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
 nMonotoneType, fContinuous, mapOptions.fNoConservation );
 }
 else
 {
 if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
 LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
 this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
 mapOptions.fNoConservation );
 }
 }
 else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
 {
 DataArray3D< double > dataGLLJacobian;
 if( !m_bPointCloudSource )
 {
 // Generate reverse node array and edge map
 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
 
 // Initialize coordinates for map
 if( eInputType == DiscretizationType_FV )
 {
 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
 }
 else
 {
 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
 DataArray3D< double > dataGLLJacobianSrc;
 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
 dataGLLJacobian );
 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
 dataGLLJacobianSrc );
 }
 }
 // else { /* Source is a point cloud dataset */ }
 
 if( !m_bPointCloudTarget )
 {
 // Generate reverse node array and edge map
 if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
 if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
 
 // Initialize coordinates for map
 if( eOutputType == DiscretizationType_FV )
 {
 this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
 }
 else
 {
 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
 dataGLLJacobian );
 }
 }
 // else { /* Target is a point cloud dataset */ }
 
 // Finite volume input / Finite element output
 rval = this->SetDOFmapAssociation(
 eInputType, ( eInputType == DiscretizationType_CGLL ),
 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrcCov ),
 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrc ), eOutputType,
 ( eOutputType == DiscretizationType_CGLL ), ( m_bPointCloudTarget ? NULL : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
 
 // Construct remap
 if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
 rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
 }
 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
 {
 DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
 
 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
 // double dNumericalAreaCov_loc =
 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
 dataGLLJacobian );
 
 double dNumericalArea_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
 dataGLLNodesSrc, dataGLLJacobianSrc );
 
 // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area:
 // %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
 // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);
 
 double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root )
 {
 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );
 
 if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
 {
 dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
 "numerical area and geometric area\n" );
 }
 }
 
 if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
 {
 _EXCEPTIONT( "Number of element does not match between metadata and "
 "input mesh" );
 }
 
 // Initialize coordinates for map
 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
 
 // Generate the continuous Jacobian for input mesh
 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
 
 if( eInputType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
 }
 
 // Finite element input / Finite volume output
 rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, false, NULL );MB_CHK_ERR( rval );
 
 // Generate remap
 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
 
 if( strMapAlgorithm == "volumetric" )
 {
 _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
 "GLL input mesh" );
 }
 
 LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
 mapOptions.fNoConservation );
 }
 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
 {
 DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
 DataArray3D< double > dataGLLJacobianOut;
 
 // Input metadata
 if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
 dataGLLJacobianIn );
 
 double dNumericalAreaSrc_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
 dataGLLNodesSrc, dataGLLJacobianSrc );
 double dNumericalAreaIn = dNumericalAreaSrc_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root )
 {
 dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn );
 
 if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
 {
 dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
 "numerical area and geometric area\n" );
 }
 }
 
 // Output metadata
 if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
 double dNumericalAreaOut_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
 dataGLLNodesDest, dataGLLJacobianOut );
 double dNumericalAreaOut = dNumericalAreaOut_loc;
#ifdef MOAB_HAVE_MPI
 if( m_pcomm )
 MPI_Reduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
 if( is_root )
 {
 dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut );
 
 if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
 {
 if( is_root )
 dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh "
 "numerical area and geometric area\n" );
 }
 }
 
 // Initialize coordinates for map
 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
 
 // Generate the continuous Jacobian for input mesh
 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
 
 if( eInputType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
 }
 
 // Generate the continuous Jacobian for output mesh
 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
 
 if( eOutputType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
 }
 
 // Input Finite Element to Output Finite Element
 rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType,
 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
 
 // Generate remap
 if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
 
 LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
 fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
 }
 else
 {
 _EXCEPTIONT( "Not implemented" );
 }
 
#ifdef MOAB_HAVE_EIGEN3
 copy_tempest_sparsemat_to_eigen3();
#endif
 
#ifdef MOAB_HAVE_MPI
 {
 // Remove ghosted entities from overlap set
 moab::Range ghostedEnts;
 rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
 moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
 rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
 }
#endif
 // Verify consistency, conservation and monotonicity, globally
 if( !mapOptions.fNoCheck )
 {
 if( is_root ) dbgprint.printf( 0, "Verifying map" );
 this->IsConsistent( 1.0e-8 );
 if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
 
 if( nMonotoneType != 0 )
 {
 this->IsMonotone( 1.0e-12 );
 }
 }
 }
 catch( Exception& e )
 {
 dbgprint.printf( 0, "%s", e.ToString().c_str() );
 return ( moab::MB_FAILURE );
 }
 catch( ... )
 {
 return ( moab::MB_FAILURE );
 }
 return moab::MB_SUCCESS;
}

References dbgprint, moab::error(), ErrorCode, MB_ALREADY_ALLOCATED, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_EXCL, MB_TYPE_DOUBLE, moab::Remapper::OverlapMesh, and size.

Referenced by main().

GetColGlobalDoF()

int moab::TempestOnlineMap::GetColGlobalDoF ( int  localID ) const

inline

Get the global Degrees-Of-Freedom ID on the source mesh.

Definition at line 481 of file TempestOnlineMap.hpp.

{
 return col_gdofmap[localColID];
}

GetDestinationGlobalNDofs()

int moab::TempestOnlineMap::GetDestinationGlobalNDofs

(

)

Get the number of total Degrees-Of-Freedom defined on the destination mesh.

GetDestinationLocalNDofs()

int moab::TempestOnlineMap::GetDestinationLocalNDofs

(

)

Get the number of local Degrees-Of-Freedom defined on the destination mesh.

GetDestinationNDofsPerElement()

int moab::TempestOnlineMap::GetDestinationNDofsPerElement ( )

inline

Get the number of Degrees-Of-Freedom per element on the destination mesh.

Definition at line 499 of file TempestOnlineMap.hpp.

{
 return m_nDofsPEl_Dest;
}

GetGlobalSourceAreas()

const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalSourceAreas

(

)

const

If we computed the reduction, get the vector representing the source areas for all entities in the mesh

GetGlobalTargetAreas()

const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalTargetAreas

(

)

const

If we computed the reduction, get the vector representing the target areas for all entities in the mesh

GetIndexOfColGlobalDoF()

int moab::TempestOnlineMap::GetIndexOfColGlobalDoF ( int  globalColDoF ) const

inline

Get the index of globaColDoF.

Definition at line 486 of file TempestOnlineMap.hpp.

{
 return globalColDoF + 1; // temporary
}

GetIndexOfRowGlobalDoF()

int moab::TempestOnlineMap::GetIndexOfRowGlobalDoF ( int  globalRowDoF ) const

inline

Get the index of globaRowDoF.

Definition at line 475 of file TempestOnlineMap.hpp.

{
 return globalRowDoF + 1;
}

GetRowGlobalDoF()

int moab::TempestOnlineMap::GetRowGlobalDoF ( int  localID ) const

inline

Get the global Degrees-Of-Freedom ID on the destination mesh.

Definition at line 470 of file TempestOnlineMap.hpp.

{
 return row_gdofmap[localRowID];
}

GetSourceGlobalNDofs()

int moab::TempestOnlineMap::GetSourceGlobalNDofs

(

)

Get the number of total Degrees-Of-Freedom defined on the source mesh.

GetSourceLocalNDofs()

int moab::TempestOnlineMap::GetSourceLocalNDofs

(

)

Get the number of local Degrees-Of-Freedom defined on the source mesh.

GetSourceNDofsPerElement()

int moab::TempestOnlineMap::GetSourceNDofsPerElement ( )

inline

Get the number of Degrees-Of-Freedom per element on the source mesh.

Definition at line 492 of file TempestOnlineMap.hpp.

{
 return m_nDofsPEl_Src;
}

IsConservative()

int moab::TempestOnlineMap::IsConservative ( double  dTolerance )

virtual

Determine if the map is conservative.

Definition at line 1496 of file TempestOnlineMap.cpp.

{
#ifndef MOAB_HAVE_MPI
 
 return OfflineMap::IsConservative( dTolerance );
 
#else
 // return OfflineMap::IsConservative(dTolerance);
 
 int ierr;
 // Get map entries
 DataArray1D< int > dataRows;
 DataArray1D< int > dataCols;
 DataArray1D< double > dataEntries;
 const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
 const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
 
 // Calculate column sums
 std::vector< int > dColumnsUnique;
 std::vector< double > dColumnSums;
 
 int nColumns = m_mapRemap.GetColumns();
 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
 dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
 dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
 
 for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
 {
 dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
 
 assert( dataCols[i] < m_nTotDofs_SrcCov );
 
 // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
 int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
 // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
 dColumnsUnique[dataCols[i]] = colGID;
 
 // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
 // std::endl;
 }
 
 int rootProc = 0;
 std::vector< int > nElementsInProc;
 const int nDATA = 3;
 if( !rank ) nElementsInProc.resize( size * nDATA );
 int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
 ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 
 int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
 std::vector< int > dColumnIndices;
 std::vector< double > dColumnSourceAreas;
 std::vector< double > dColumnSumsTotal;
 std::vector< int > displs, rcount;
 if( rank == rootProc )
 {
 displs.resize( size + 1, 0 );
 rcount.resize( size, 0 );
 int gsum = 0;
 for( int ir = 0; ir < size; ++ir )
 {
 nTotVals += nElementsInProc[ir * nDATA];
 nTotColumns += nElementsInProc[ir * nDATA + 1];
 nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
 
 displs[ir] = gsum;
 rcount[ir] = nElementsInProc[ir * nDATA + 1];
 gsum += rcount[ir];
 
 // printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns,
 // displs[ir], rcount[ir], gsum);
 }
 
 dColumnIndices.resize( nTotColumns, -1 );
 dColumnSumsTotal.resize( nTotColumns, 0.0 );
 // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
 }
 
 // Gather all ColumnSums to root process and accumulate
 // We expect that the sums of all columns equate to 1.0 within user specified tolerance
 // Need to do a gatherv here since different processes have different number of elements
 // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
 // MPI_SUM, 0, m_pcomm->comm());
 ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
 displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
 displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
 // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
 // MPI_SUCCESS ) return -1;
 
 // Clean out unwanted arrays now
 dColumnSums.clear();
 dColumnsUnique.clear();
 
 // Verify all column sums equal the input Jacobian
 int fConservative = 0;
 if( rank == rootProc )
 {
 displs[size] = ( nTotColumns );
 // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
 std::map< int, double > dColumnSumsOnRoot;
 // std::map<int, double> dColumnSourceAreasOnRoot;
 for( int ir = 0; ir < size; ir++ )
 {
 for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
 {
 if( dColumnIndices[ips] < 0 ) continue;
 // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
 assert( dColumnIndices[ips] < nTotColumnsUnq );
 dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
 // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
 // dColumnSourceAreas[ dColumnIndices[ips] ]
 }
 }
 
 for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
 {
 // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
 if( fabs( it->second - 1.0 ) > dTolerance )
 {
 fConservative++;
 Announce( "TempestOnlineMap is not conservative in column "
 // "%i (%1.15e)", it->first, it->second );
 "%i (%1.15e)",
 it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
 }
 }
 }
 
 // TODO: Just do a broadcast from root instead of a reduction
 ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 
 return fConservative;
#endif
}

References size.

IsConsistent()

int moab::TempestOnlineMap::IsConsistent ( double  dTolerance )

virtual

Determine if the map is first-order accurate.

Definition at line 1450 of file TempestOnlineMap.cpp.

{
#ifndef MOAB_HAVE_MPI
 
 return OfflineMap::IsConsistent( dTolerance );
 
#else
 
 // Get map entries
 DataArray1D< int > dataRows;
 DataArray1D< int > dataCols;
 DataArray1D< double > dataEntries;
 
 // Calculate row sums
 DataArray1D< double > dRowSums;
 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
 dRowSums.Allocate( m_mapRemap.GetRows() );
 
 for( unsigned i = 0; i < dataRows.GetRows(); i++ )
 {
 dRowSums[dataRows[i]] += dataEntries[i];
 }
 
 // Verify all row sums are equal to 1
 int fConsistent = 0;
 for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
 {
 if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
 {
 fConsistent++;
 int rowGID = row_gdofmap[i];
 Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
 }
 }
 
 int ierr;
 int fConsistentGlobal = 0;
 ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 
 return fConsistentGlobal;
#endif
}

IsMonotone()

int moab::TempestOnlineMap::IsMonotone ( double  dTolerance )

virtual

Determine if the map is monotone.

Definition at line 1638 of file TempestOnlineMap.cpp.

{
#ifndef MOAB_HAVE_MPI
 
 return OfflineMap::IsMonotone( dTolerance );
 
#else
 
 // Get map entries
 DataArray1D< int > dataRows;
 DataArray1D< int > dataCols;
 DataArray1D< double > dataEntries;
 
 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
 
 // Verify all entries are in the range [0,1]
 int fMonotone = 0;
 for( unsigned i = 0; i < dataRows.GetRows(); i++ )
 {
 if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
 {
 fMonotone++;
 
 Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
 }
 }
 
 int ierr;
 int fMonotoneGlobal = 0;
 ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
 if( ierr != MPI_SUCCESS ) return -1;
 
 return fMonotoneGlobal;
#endif
}

LinearRemapFVtoFV_Tempest_MOAB()

void moab::TempestOnlineMap::LinearRemapFVtoFV_Tempest_MOAB ( int  nOrder )

private

Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh.

Definition at line 115 of file TempestLinearRemap.cpp.

{
 // Order of triangular quadrature rule
 const int TriQuadRuleOrder = 4;
 
 // Verify ReverseNodeArray has been calculated
 if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
 {
 _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInput" );
 }
 
 // Triangular quadrature rule
 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
 
 // Number of coefficients needed at this order
#ifdef RECTANGULAR_TRUNCATION
 int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
 int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif
 
 // Number of faces you need
 const int nRequiredFaceSetSize = nCoefficients;
 
 // Fit weight exponent
 const int nFitWeightsExponent = nOrder + 2;
 
 // Announcemnets
 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
 dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
 if( is_root )
 {
 dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
 dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
 dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
 dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
 }
 
 // Current overlap face
 int ixOverlap = 0;
#ifdef VERBOSE
 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
#endif
 DataArray2D< double > dIntArray;
 DataArray1D< double > dConstraint( nCoefficients );
 
 // Loop through all faces on m_meshInput
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
 // Output every 1000 elements
#ifdef VERBOSE
 if( ixFirst % outputFrequency == 0 && is_root )
 {
 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
 }
#endif
 // Find the set of Faces that overlap faceFirst
 int ixOverlapBegin = ixOverlap;
 unsigned ixOverlapEnd = ixOverlapBegin;
 
 for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
 {
 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 ) break;
 }
 
 unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
 
 if( nOverlapFaces == 0 ) continue;
 
 // Build integration array
 BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
 nOrder, dIntArray );
 
 // Set of Faces to use in building the reconstruction and associated
 // distance metric.
 AdjacentFaceVector vecAdjFaces;
 
 GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
 
 // Number of adjacent Faces
 int nAdjFaces = vecAdjFaces.size();
 
 // Determine the conservative constraint equation
 double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
 dConstraint.Zero();
 for( int p = 0; p < nCoefficients; p++ )
 {
 for( unsigned j = 0; j < nOverlapFaces; j++ )
 {
 dConstraint[p] += dIntArray[p][j];
 }
 dConstraint[p] /= dFirstArea;
 }
 
 // Build the fit array from the integration operator
 DataArray2D< double > dFitArray;
 DataArray1D< double > dFitWeights;
 DataArray2D< double > dFitArrayPlus;
 
 BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
 dFitArray, dFitWeights );
 
 // Compute the inverse fit array
 bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
 
 // Multiply integration array and fit array
 DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
 if( fSuccess )
 {
 // Multiply integration array and inverse fit array
 for( int i = 0; i < nAdjFaces; i++ )
 {
 for( size_t j = 0; j < nOverlapFaces; j++ )
 {
 for( int k = 0; k < nCoefficients; k++ )
 {
 dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
 }
 }
 }
 
 // Unable to invert fit array, drop to 1st order. In this case
 // dFitArrayPlus(0,0) = 1 and all other entries are zero.
 }
 else
 {
 dComposedArray.Zero();
 for( size_t j = 0; j < nOverlapFaces; j++ )
 {
 dComposedArray( 0, j ) += dIntArray( 0, j );
 }
 }
 
 // Put composed array into map
 for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
 {
 for( unsigned j = 0; j < nOverlapFaces; j++ )
 {
 int& ixFirstFaceLoc = vecAdjFaces[i].first;
 int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
 // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
 // ixFirstFaceLoc); int ixSecondFaceGlob =
 // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
 
 // signal to not participate, because it is a ghost target
 if( ixSecondFaceLoc < 0 ) continue; // do not do anything
 
 m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
 dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
 }
 }
 
 // Increment the current overlap index
 ixOverlap += nOverlapFaces;
 }
 
 return;
}

References dbgprint.

LinearRemapFVtoGLL_MOAB()

void moab::TempestOnlineMap::LinearRemapFVtoGLL_MOAB ( const DataArray3D< int > &  dataGLLNodes, const DataArray3D< double > &  dataGLLJacobian, const DataArray1D< double > &  dataGLLNodalArea, int  nOrder, int  nMonotoneType, bool  fContinuous, bool  fNoConservation  )

private

Generate the OfflineMap for remapping from finite volumes to finite elements.

LinearRemapGLLtoGLL2_MOAB()

void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB ( const DataArray3D< int > &  dataGLLNodesIn, const DataArray3D< double > &  dataGLLJacobianIn, const DataArray3D< int > &  dataGLLNodesOut, const DataArray3D< double > &  dataGLLJacobianOut, const DataArray1D< double > &  dataNodalAreaOut, int  nPin, int  nPout, int  nMonotoneType, bool  fContinuousIn, bool  fContinuousOut, bool  fNoConservation  )

private

Generate the OfflineMap for remapping from finite elements to finite elements.

Definition at line 895 of file TempestLinearRemap.cpp.

{
 // Triangular quadrature rule
 TriangularQuadratureRule triquadrule( 8 );
 
 const DataArray2D< double >& dG = triquadrule.GetG();
 const DataArray1D< double >& dW = triquadrule.GetW();
 
 // Get SparseMatrix represntation of the OfflineMap
 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
 
 // Sample coefficients
 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
 DataArray2D< double > dSampleCoeffOut( nPout, nPout );
 
 // Announcemnets
 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
 if( is_root )
 {
 dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
 }
 
 // Build the integration array for each element on m_meshOverlap
 DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
 
 // Number of overlap Faces per source Face
 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
 
 int ixOverlap = 0;
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
 // Determine how many overlap Faces and triangles are present
 int nOverlapFaces = 0;
 size_t ixOverlapTemp = ixOverlap;
 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
 {
 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
 {
 break;
 }
 
 nOverlapFaces++;
 }
 
 nAllOverlapFaces[ixFirst] = nOverlapFaces;
 
 // Increment the current overlap index
 ixOverlap += nAllOverlapFaces[ixFirst];
 }
 
 // Geometric area of each output node
 DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
 
 // Area of each overlap element in the output basis
 DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
 
 // Loop through all faces on m_meshInput
 ixOverlap = 0;
#ifdef VERBOSE
 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
#endif
 if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
 
 // generic triangle used for area computation, for triangles around the center of overlap face;
 // used for overlap faces with more than 4 edges;
 // nodes array will be set for each triangle;
 // these triangles are not part of the mesh structure, they are just temporary during
 // aforementioned decomposition.
 Face faceTri( 3 );
 NodeVector nodes( 3 );
 faceTri.SetNode( 0, 0 );
 faceTri.SetNode( 1, 1 );
 faceTri.SetNode( 2, 2 );
 
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
#ifdef VERBOSE
 // Announce computation progress
 if( ixFirst % outputFrequency == 0 && is_root )
 {
 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
 }
#endif
 // Quantities from the First Mesh
 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
 
 const NodeVector& nodesFirst = m_meshInputCov->nodes;
 
 // Number of overlapping Faces and triangles
 int nOverlapFaces = nAllOverlapFaces[ixFirst];
 
 if( !nOverlapFaces ) continue;
 
 // // Calculate total element Jacobian
 // double dTotalJacobian = 0.0;
 // for (int s = 0; s < nPin; s++) {
 // for (int t = 0; t < nPin; t++) {
 // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
 // }
 // }
 
 // Loop through all Overlap Faces
 for( int i = 0; i < nOverlapFaces; i++ )
 {
 // Quantities from the overlap Mesh
 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
 
 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
 
 // Quantities from the Second Mesh
 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
 
 // signal to not participate, because it is a ghost target
 if( ixSecond < 0 ) continue; // do not do anything
 
 const NodeVector& nodesSecond = m_meshOutput->nodes;
 
 const Face& faceSecond = m_meshOutput->faces[ixSecond];
 
 int nbEdges = faceOverlap.edges.size();
 int nOverlapTriangles = 1;
 Node center; // not used if nbEdges == 3
 if( nbEdges > 3 )
 { // decompose from center in this case
 nOverlapTriangles = nbEdges;
 for( int k = 0; k < nbEdges; k++ )
 {
 const Node& node = nodesOverlap[faceOverlap[k]];
 center = center + node;
 }
 center = center / nbEdges;
 center = center.Normalized(); // project back on sphere of radius 1
 }
 
 Node node0, node1, node2;
 double dTriArea;
 
 // Loop over all sub-triangles of this Overlap Face
 for( int j = 0; j < nOverlapTriangles; j++ )
 {
 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
 {
 node0 = nodesOverlap[faceOverlap[0]];
 node1 = nodesOverlap[faceOverlap[1]];
 node2 = nodesOverlap[faceOverlap[2]];
 dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
 }
 else // decompose polygon in triangles around the center
 {
 node0 = center;
 node1 = nodesOverlap[faceOverlap[j]];
 int j1 = ( j + 1 ) % nbEdges;
 node2 = nodesOverlap[faceOverlap[j1]];
 nodes[0] = center;
 nodes[1] = node1;
 nodes[2] = node2;
 dTriArea = CalculateFaceArea( faceTri, nodes );
 }
 
 for( int k = 0; k < triquadrule.GetPoints(); k++ )
 {
 // Get the nodal location of this point
 double dX[3];
 
 dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
 dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
 dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
 
 double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
 
 dX[0] /= dMag;
 dX[1] /= dMag;
 dX[2] /= dMag;
 
 Node nodeQuadrature( dX[0], dX[1], dX[2] );
 
 // Find the components of this quadrature point in the basis
 // of the first Face.
 double dAlphaIn;
 double dBetaIn;
 
 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
 
 // Find the components of this quadrature point in the basis
 // of the second Face.
 double dAlphaOut;
 double dBetaOut;
 
 ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
 
 /*
 // Check inverse map value
 if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
 (dBetaIn < 0.0) || (dBetaIn > 1.0)
 ) {
 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
 dAlphaIn, dBetaIn);
 }
 
 // Check inverse map value
 if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
 (dBetaOut < 0.0) || (dBetaOut > 1.0)
 ) {
 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
 dAlphaOut, dBetaOut);
 }
 */
 // Sample the First finite element at this point
 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
 
 // Sample the Second finite element at this point
 SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
 
 // Overlap output area
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
 
 dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
 
 dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
 }
 }
 
 // Compute overlap integral
 int ixp = 0;
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 int ixs = 0;
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 // Sample the Second finite element at this point
 dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
 dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
 
 ixs++;
 }
 }
 
 ixp++;
 }
 }
 }
 }
 }
 
 // Coefficients
 DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
 
 for( int i = 0; i < nOverlapFaces; i++ )
 {
 // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
 
 int ixp = 0;
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 int ixs = 0;
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
 dOverlapOutputArea[ixOverlap + i][s * nPout + t];
 
 ixs++;
 }
 }
 
 ixp++;
 }
 }
 }
 
 // Source areas
 DataArray1D< double > vecSourceArea( nPin * nPin );
 
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
 }
 }
 
 // Target areas
 DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
 
 for( int i = 0; i < nOverlapFaces; i++ )
 {
 // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
 int ixs = 0;
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
 
 ixs++;
 }
 }
 }
 
 // Force consistency and conservation
 if( !fNoConservation )
 {
 ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
 }
 
 // Update global coefficients
 for( int i = 0; i < nOverlapFaces; i++ )
 {
 int ixp = 0;
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 int ixs = 0;
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 dGlobalIntArray[ixp][ixOverlap + i][ixs] =
 dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
 
 ixs++;
 }
 }
 
 ixp++;
 }
 }
 }
 
#ifdef VVERBOSE
 // Check column sums (conservation)
 for( int i = 0; i < nPin * nPin; i++ )
 {
 double dColSum = 0.0;
 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
 {
 dColSum += dCoeff[j][i] * vecTargetArea[j];
 }
 printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
 }
 
 // Check row sums (consistency)
 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
 {
 double dRowSum = 0.0;
 for( int i = 0; i < nPin * nPin; i++ )
 {
 dRowSum += dCoeff[j][i];
 }
 printf( "Row %i: %1.15e\n", j, dRowSum );
 }
#endif
 
 // Increment the current overlap index
 ixOverlap += nOverlapFaces;
 }
 
 // Build redistribution map within target element
 if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
 DataArray1D< double > dRedistSourceArea( nPout * nPout );
 DataArray1D< double > dRedistTargetArea( nPout * nPout );
 std::vector< DataArray2D< double > > dRedistributionMaps;
 dRedistributionMaps.resize( m_meshOutput->faces.size() );
 
 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
 {
 dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
 
 for( int i = 0; i < nPout * nPout; i++ )
 {
 dRedistributionMaps[ixSecond][i][i] = 1.0;
 }
 
 for( int s = 0; s < nPout * nPout; s++ )
 {
 dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
 }
 
 for( int s = 0; s < nPout * nPout; s++ )
 {
 dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
 }
 
 if( !fNoConservation )
 {
 ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
 ( nMonotoneType != 0 ) );
 
 for( int s = 0; s < nPout * nPout; s++ )
 {
 for( int t = 0; t < nPout * nPout; t++ )
 {
 dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
 }
 }
 }
 }
 
 // Construct the total geometric area
 DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
 {
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
 dGeometricOutputArea[ixSecond][s * nPout + t];
 }
 }
 }
 
 // Compose the integration operator with the output map
 ixOverlap = 0;
 
 if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
 
 // Map from source DOFs to target DOFs with redistribution applied
 DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
 
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
#ifdef VERBOSE
 // Announce computation progress
 if( ixFirst % outputFrequency == 0 && is_root )
 {
 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
 }
#endif
 // Number of overlapping Faces and triangles
 int nOverlapFaces = nAllOverlapFaces[ixFirst];
 
 if( !nOverlapFaces ) continue;
 
 // Put composed array into map
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
 
 // signal to not participate, because it is a ghost target
 if( ixSecondFace < 0 ) continue; // do not do anything
 
 dRedistributedOp.Zero();
 for( int p = 0; p < nPin * nPin; p++ )
 {
 for( int s = 0; s < nPout * nPout; s++ )
 {
 for( int t = 0; t < nPout * nPout; t++ )
 {
 dRedistributedOp[p][s] +=
 dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
 }
 }
 }
 
 int ixp = 0;
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 int ixFirstNode;
 if( fContinuousIn )
 {
 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
 }
 else
 {
 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
 }
 
 int ixs = 0;
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 int ixSecondNode;
 if( fContinuousOut )
 {
 ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
 
 if( !fNoConservation )
 {
 smatMap( ixSecondNode, ixFirstNode ) +=
 dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
 }
 else
 {
 smatMap( ixSecondNode, ixFirstNode ) +=
 dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
 }
 }
 else
 {
 ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
 
 if( !fNoConservation )
 {
 smatMap( ixSecondNode, ixFirstNode ) +=
 dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
 }
 else
 {
 smatMap( ixSecondNode, ixFirstNode ) +=
 dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
 }
 }
 
 ixs++;
 }
 }
 
 ixp++;
 }
 }
 }
 
 // Increment the current overlap index
 ixOverlap += nOverlapFaces;
 }
 
 return;
}

References center(), dbgprint, ForceIntArrayConsistencyConservation(), and t.

LinearRemapGLLtoGLL2_Pointwise_MOAB()

void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB ( const DataArray3D< int > &  dataGLLNodesIn, const DataArray3D< double > &  dataGLLJacobianIn, const DataArray3D< int > &  dataGLLNodesOut, const DataArray3D< double > &  dataGLLJacobianOut, const DataArray1D< double > &  dataNodalAreaOut, int  nPin, int  nPout, int  nMonotoneType, bool  fContinuousIn, bool  fContinuousOut  )

private

Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation).

Definition at line 1446 of file TempestLinearRemap.cpp.

{
 // Gauss-Lobatto quadrature within Faces
 DataArray1D< double > dGL;
 DataArray1D< double > dWL;
 
 GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
 
 // Get SparseMatrix represntation of the OfflineMap
 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
 
 // Sample coefficients
 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
 
 // Announcemnets
 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
 if( is_root )
 {
 dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
 }
 
 // Number of overlap Faces per source Face
 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
 
 int ixOverlap = 0;
 
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
 size_t ixOverlapTemp = ixOverlap;
 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
 {
 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
 
 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
 
 nAllOverlapFaces[ixFirst]++;
 }
 
 // Increment the current overlap index
 ixOverlap += nAllOverlapFaces[ixFirst];
 }
 
 // Number of times this point was found
 DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
 
 ixOverlap = 0;
#ifdef VERBOSE
 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
#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 && is_root )
 {
 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
 }
#endif
 // Quantities from the First Mesh
 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
 
 const NodeVector& nodesFirst = m_meshInputCov->nodes;
 
 // Number of overlapping Faces and triangles
 int nOverlapFaces = nAllOverlapFaces[ixFirst];
 
 // Loop through all Overlap Faces
 for( int i = 0; i < nOverlapFaces; i++ )
 {
 // Quantities from the Second Mesh
 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
 
 // signal to not participate, because it is a ghost target
 if( ixSecond < 0 ) continue; // do not do anything
 
 const NodeVector& nodesSecond = m_meshOutput->nodes;
 const Face& faceSecond = m_meshOutput->faces[ixSecond];
 
 // Loop through all nodes on the second face
 for( int s = 0; s < nPout; s++ )
 {
 for( int t = 0; t < nPout; t++ )
 {
 size_t ixSecondNode;
 if( fContinuousOut )
 {
 ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
 }
 else
 {
 ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
 }
 
 if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
 
 // Check if this node has been found already
 if( fSecondNodeFound[ixSecondNode] ) continue;
 
 // Check this node
 Node node;
 Node dDx1G;
 Node dDx2G;
 
 ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
 
 // Find the components of this quadrature point in the basis
 // of the first Face.
 double dAlphaIn;
 double dBetaIn;
 
 ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
 
 // Check if this node is within the first Face
 if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
 ( dBetaIn > 1.0 + 1.0e-10 ) )
 continue;
 
 // Node is within the overlap region, mark as found
 fSecondNodeFound[ixSecondNode] = true;
 
 // Sample the First finite element at this point
 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
 
 // Add to map
 for( int p = 0; p < nPin; p++ )
 {
 for( int q = 0; q < nPin; q++ )
 {
 int ixFirstNode;
 if( fContinuousIn )
 {
 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
 }
 else
 {
 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
 }
 
 smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
 }
 }
 }
 }
 }
 
 // Increment the current overlap index
 ixOverlap += nOverlapFaces;
 }
 
 // Check for missing samples
 for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
 {
 if( !fSecondNodeFound[i] )
 {
 _EXCEPTION1( "Can't sample point %i", i );
 }
 }
 
 return;
}

References dbgprint, and t.

LinearRemapNN_MOAB()

moab::ErrorCode moab::TempestOnlineMap::LinearRemapNN_MOAB ( bool  use_GID_matching = false, bool  strict_check = false  )

private

Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh.

Definition at line 52 of file TempestLinearRemap.cpp.

{
 /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov) */
 
#ifdef VVERBOSE
 {
 std::ofstream output_file( "rowcolindices.txt", std::ios::out );
 output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
 << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
 output_file << "Rows \n";
 for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
 output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
 output_file << "Cols \n";
 for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
 output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
 output_file.flush(); // required here
 output_file.close();
 }
#endif
 
 if( use_GID_matching )
 {
 std::map< unsigned, unsigned > src_gl;
 for( unsigned it = 0; it < col_gdofmap.size(); ++it )
 src_gl[col_gdofmap[it]] = it;
 
 std::map< unsigned, unsigned >::iterator iter;
 for( unsigned it = 0; it < row_gdofmap.size(); ++it )
 {
 unsigned row = row_gdofmap[it];
 iter = src_gl.find( row );
 if( strict_check && iter == src_gl.end() )
 {
 std::cout << "Searching for global target DOF " << row
 << " but could not find correspondence in source mesh.\n";
 assert( false );
 }
 else if( iter == src_gl.end() )
 {
 continue;
 }
 else
 {
 unsigned icol = src_gl[row];
 unsigned irow = it;
 
 // Set the permutation matrix in local space
 m_mapRemap( irow, icol ) = 1.0;
 }
 }
 
 return moab::MB_SUCCESS;
 }
 else
 {
 /* Create a Kd-tree to perform local queries to find nearest neighbors */
 
 return moab::MB_SUCCESS;
 }
}

References col_gdofmap, m_nTotDofs_Dest, m_nTotDofs_SrcCov, MB_SUCCESS, and row_gdofmap.

LinearRemapSE0_Tempest_MOAB()

void moab::TempestOnlineMap::LinearRemapSE0_Tempest_MOAB ( const DataArray3D< int > &  dataGLLNodes, const DataArray3D< double > &  dataGLLJacobian  )

private

Generate the OfflineMap for linear conserative element-average spectral element to element average remapping.

LinearRemapSE4_Tempest_MOAB()

void moab::TempestOnlineMap::LinearRemapSE4_Tempest_MOAB ( const DataArray3D< int > &  dataGLLNodes, const DataArray3D< double > &  dataGLLJacobian, int  nMonotoneType, bool  fContinuousIn, bool  fNoConservation  )

private

Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping.

Definition at line 469 of file TempestLinearRemap.cpp.

{
 // Order of the polynomial interpolant
 int nP = dataGLLNodes.GetRows();
 
 // Order of triangular quadrature rule
 const int TriQuadRuleOrder = 4;
 
 // Triangular quadrature rule
 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
 
 int TriQuadraturePoints = triquadrule.GetPoints();
 
 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
 
 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
 
 // Sample coefficients
 DataArray2D< double > dSampleCoeff( nP, nP );
 
 // GLL Quadrature nodes on quadrilateral elements
 DataArray1D< double > dG;
 DataArray1D< double > dW;
 GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
 
 // Announcements
 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
 dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
 if( is_root )
 {
 dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
 dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
 }
 
 // Get SparseMatrix represntation of the OfflineMap
 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
 
 // NodeVector from m_meshOverlap
 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
 const NodeVector& nodesFirst = m_meshInputCov->nodes;
 
 // Vector of source areas
 DataArray1D< double > vecSourceArea( nP * nP );
 
 DataArray1D< double > vecTargetArea;
 DataArray2D< double > dCoeff;
 
#ifdef VERBOSE
 std::stringstream sstr;
 sstr << "remapdata_" << rank << ".txt";
 std::ofstream output_file( sstr.str() );
#endif
 
 // Current Overlap Face
 int ixOverlap = 0;
#ifdef VERBOSE
 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
#endif
 // generic triangle used for area computation, for triangles around the center of overlap face;
 // used for overlap faces with more than 4 edges;
 // nodes array will be set for each triangle;
 // these triangles are not part of the mesh structure, they are just temporary during
 // aforementioned decomposition.
 Face faceTri( 3 );
 NodeVector nodes( 3 );
 faceTri.SetNode( 0, 0 );
 faceTri.SetNode( 1, 1 );
 faceTri.SetNode( 2, 2 );
 
 // Loop over all input Faces
 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
 {
 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
 
 if( faceFirst.edges.size() != 4 )
 {
 _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
 }
#ifdef VERBOSE
 // Announce computation progress
 if( ixFirst % outputFrequency == 0 && is_root )
 {
 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
 }
#endif
 // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
 // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
 // However, the relation with MOAB and Tempest will go out of the roof
 
 // Determine how many overlap Faces and triangles are present
 int nOverlapFaces = 0;
 size_t ixOverlapTemp = ixOverlap;
 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
 {
 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
 {
 break;
 }
 
 nOverlapFaces++;
 }
 
 // No overlaps
 if( nOverlapFaces == 0 )
 {
 continue;
 }
 
 // Allocate remap coefficients array for meshFirst Face
 DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
 
 // Find the local remap coefficients
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
 if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
 {
 Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
 m_meshOverlap->vecFaceArea[ixOverlap + j] );
 int n = faceOverlap.edges.size();
 Announce( "Number nodes: %d", n );
 for( int k = 0; k < n; k++ )
 {
 Node nd = nodesOverlap[faceOverlap[k]];
 Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
 }
 continue;
 }
 
 // #ifdef VERBOSE
 // if ( is_root )
 // Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
 // m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
 // m_meshOverlap->vecFaceArea[ixOverlap + j] );
 // #endif
 
 int nbEdges = faceOverlap.edges.size();
 int nOverlapTriangles = 1;
 Node center; // not used if nbEdges == 3
 if( nbEdges > 3 )
 { // decompose from center in this case
 nOverlapTriangles = nbEdges;
 for( int k = 0; k < nbEdges; k++ )
 {
 const Node& node = nodesOverlap[faceOverlap[k]];
 center = center + node;
 }
 center = center / nbEdges;
 center = center.Normalized(); // project back on sphere of radius 1
 }
 
 Node node0, node1, node2;
 double dTriangleArea;
 
 // Loop over all sub-triangles of this Overlap Face
 for( int k = 0; k < nOverlapTriangles; k++ )
 {
 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
 {
 node0 = nodesOverlap[faceOverlap[0]];
 node1 = nodesOverlap[faceOverlap[1]];
 node2 = nodesOverlap[faceOverlap[2]];
 dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
 }
 else // decompose polygon in triangles around the center
 {
 node0 = center;
 node1 = nodesOverlap[faceOverlap[k]];
 int k1 = ( k + 1 ) % nbEdges;
 node2 = nodesOverlap[faceOverlap[k1]];
 nodes[0] = center;
 nodes[1] = node1;
 nodes[2] = node2;
 dTriangleArea = CalculateFaceArea( faceTri, nodes );
 }
 // Coordinates of quadrature Node
 for( int l = 0; l < TriQuadraturePoints; l++ )
 {
 Node nodeQuadrature;
 nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
 TriQuadratureG[l][2] * node2.x;
 
 nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
 TriQuadratureG[l][2] * node2.y;
 
 nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
 TriQuadratureG[l][2] * node2.z;
 
 nodeQuadrature = nodeQuadrature.Normalized();
 
 // Find components of quadrature point in basis
 // of the first Face
 double dAlpha;
 double dBeta;
 
 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
 
 // Check inverse map value
 if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
 ( dBeta > 1.0 + 1.0e-13 ) )
 {
 _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
 "(%1.5e %1.5e)",
 j, l, dAlpha, dBeta );
 }
 
 // Sample the finite element at this point
 SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
 
 // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
 m_meshOverlap->vecFaceArea[ixOverlap + j];
 }
 }
 }
 }
 }
 
#ifdef VERBOSE
 output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 output_file << dRemapCoeff[p][q][j] << " ";
 }
 }
 }
 output_file << std::endl;
#endif
 
 // Force consistency and conservation
 if( !fNoConservation )
 {
 double dTargetArea = 0.0;
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
 }
 
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
 }
 }
 
 const double areaTolerance = 1e-10;
 // Source elements are completely covered by target volumes
 if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
 {
 vecTargetArea.Allocate( nOverlapFaces );
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
 }
 
 dCoeff.Allocate( nOverlapFaces, nP * nP );
 
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
 }
 }
 }
 
 // Target volumes only partially cover source elements
 }
 else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
 {
 double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
 
 vecTargetArea.Allocate( nOverlapFaces + 1 );
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
 }
 vecTargetArea[nOverlapFaces] = dExtraneousArea;
 
#ifdef VERBOSE
 Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
 m_meshInputCov->vecFaceArea[ixFirst] );
#endif
 if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
 {
 _EXCEPTIONT( "Partial element area exceeds total element area" );
 }
 
 dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
 
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
 }
 }
 }
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
 }
 }
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dCoeff[nOverlapFaces][p * nP + q] -=
 dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
 }
 }
 }
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
 }
 }
 
 // Source elements only partially cover target volumes
 }
 else
 {
 Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
 m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
 _EXCEPTIONT( "Target grid must be a subset of source grid" );
 }
 
 ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
 /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
 
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
 }
 }
 }
 }
 
#ifdef VERBOSE
 // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
 // for ( int j = 0; j < nOverlapFaces; j++ )
 // {
 // for ( int p = 0; p < nP; p++ )
 // {
 // for ( int q = 0; q < nP; q++ )
 // {
 // output_file << dRemapCoeff[p][q][j] << " ";
 // }
 // }
 // }
 // output_file << std::endl;
#endif
 
 // Put these remap coefficients into the SparseMatrix map
 for( int j = 0; j < nOverlapFaces; j++ )
 {
 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
 
 // signal to not participate, because it is a ghost target
 if( ixSecondFace < 0 ) continue; // do not do anything
 
 for( int p = 0; p < nP; p++ )
 {
 for( int q = 0; q < nP; q++ )
 {
 if( fContinuousIn )
 {
 int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
 
 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
 m_meshOverlap->vecFaceArea[ixOverlap + j] /
 m_meshOutput->vecFaceArea[ixSecondFace];
 }
 else
 {
 int ixFirstNode = ixFirst * nP * nP + p * nP + q;
 
 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
 m_meshOverlap->vecFaceArea[ixOverlap + j] /
 m_meshOutput->vecFaceArea[ixSecondFace];
 }
 }
 }
 }
 // Increment the current overlap index
 ixOverlap += nOverlapFaces;
 }
#ifdef VERBOSE
 output_file.flush(); // required here
 output_file.close();
#endif
 
 return;
}

References center(), dbgprint, and ForceConsistencyConservation3().

ReadParallelMap()

moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap	(	const char *	strSource,
		const std::vector< int > &	owned_dof_ids,
		bool	row_major_ownership = true
	)

Generate the metadata associated with the offline map.

Read the OfflineMap from a NetCDF file.

Definition at line 1195 of file TempestOnlineMapIO.cpp.

{
 NcError error( NcError::silent_nonfatal );
 
 NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
 int nS = 0, nA = 0, nB = 0;
#ifdef MOAB_HAVE_PNETCDF
 // some variables will be used just in the case netcdfpar reader fails
 int ncfile = -1, ret = 0;
 int ndims, nvars, ngatts, unlimited;
#endif
#ifdef MOAB_HAVE_NETCDFPAR
 bool is_independent = true;
 ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
 NcFile ncMap( strSource, NcFile::ReadOnly );
#endif
 
#define CHECK_EXCEPTION( obj, type, varstr ) \
 { \
 if( obj == NULL ) \
 { \
 _EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
 } \
 }
 
 // Read SparseMatrix entries
 
 if( ncMap.is_valid() )
 {
 NcDim* dimNS = ncMap.get_dim( "n_s" );
 CHECK_EXCEPTION( dimNS, "dimension", "n_s" );
 
 NcDim* dimNA = ncMap.get_dim( "n_a" );
 CHECK_EXCEPTION( dimNA, "dimension", "n_a" );
 
 NcDim* dimNB = ncMap.get_dim( "n_b" );
 CHECK_EXCEPTION( dimNB, "dimension", "n_b" );
 
 // store total number of nonzeros
 nS = dimNS->size();
 nA = dimNA->size();
 nB = dimNB->size();
 
 varRow = ncMap.get_var( "row" );
 CHECK_EXCEPTION( varRow, "variable", "row" );
 
 varCol = ncMap.get_var( "col" );
 CHECK_EXCEPTION( varCol, "variable", "col" );
 
 varS = ncMap.get_var( "S" );
 CHECK_EXCEPTION( varS, "variable", "S" );
 
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varRow, is_independent );
 ncMap.enable_var_par_access( varCol, is_independent );
 ncMap.enable_var_par_access( varS, is_independent );
#endif
 }
 else
 {
#ifdef MOAB_HAVE_PNETCDF
 // read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
 // why build wth pnetcdf without MPI ?
 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
 ret = ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile );
 ERR_PARNC( ret ); // bail out completely
 ret = ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited );
 ERR_PARNC( ret );
 // find dimension ids for n_S
 int ins;
 ret = ncmpi_inq_dimid( ncfile, "n_s", &ins );
 ERR_PARNC( ret );
 MPI_Offset leng;
 ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
 ERR_PARNC( ret );
 nS = (int)leng;
 ret = ncmpi_inq_dimid( ncfile, "n_b", &ins );
 ERR_PARNC( ret );
 ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
 ERR_PARNC( ret );
 nB = (int)leng;
#else
 _EXCEPTION1( "cannot read the file %s", strSource );
#endif
 }
 
 // Let us declare the map object for every process
 SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
 
 int localSize = nS / size;
 long offsetRead = rank * localSize;
 // leftovers on last rank
 if( rank == size - 1 )
 {
 localSize += nS % size;
 }
 
 std::vector< int > vecRow, vecCol;
 std::vector< double > vecS;
 vecRow.resize( localSize );
 vecCol.resize( localSize );
 vecS.resize( localSize );
 
 if( ncMap.is_valid() )
 {
 varRow->set_cur( (long)( offsetRead ) );
 varRow->get( &( vecRow[0] ), localSize );
 
 varCol->set_cur( (long)( offsetRead ) );
 varCol->get( &( vecCol[0] ), localSize );
 
 varS->set_cur( (long)( offsetRead ) );
 varS->get( &( vecS[0] ), localSize );
 
 ncMap.close();
 }
 else
 {
#ifdef MOAB_HAVE_PNETCDF
 // fill the local vectors with the variables from pnetcdf file; first inquire, then fill
 MPI_Offset start = (MPI_Offset)offsetRead;
 MPI_Offset count = (MPI_Offset)localSize;
 int varid;
 ret = ncmpi_inq_varid( ncfile, "S", &varid );
 ERR_PARNC( ret );
 ret = ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] );
 ERR_PARNC( ret );
 ret = ncmpi_inq_varid( ncfile, "row", &varid );
 ERR_PARNC( ret );
 ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] );
 ERR_PARNC( ret );
 ret = ncmpi_inq_varid( ncfile, "col", &varid );
 ERR_PARNC( ret );
 ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] );
 ERR_PARNC( ret );
 ret = ncmpi_close( ncfile );
 ERR_PARNC( ret );
#endif
 }
 
 // Now let us set the necessary global-to-local ID maps so that A*x operations
 // can be performed cleanly as if map was computed online
 row_dtoc_dofmap.clear();
 // row_dtoc_dofmap.reserve( nB / size );
 col_dtoc_dofmap.clear();
 rowMap.clear();
 colMap.clear();
 // col_dtoc_dofmap.reserve( 2 * nA / size );
 // row_dtoc_dofmap.resize( m_nTotDofs_Dest, UINT_MAX );
 // col_dtoc_dofmap.resize( m_nTotDofs_SrcCov, UINT_MAX );
 
#ifdef MOAB_HAVE_MPI
 // bother with tuple list only if size > 1
 // otherwise, just fill the sparse matrix
 if( size > 1 )
 {
 std::vector< int > ownership;
 // the default trivial partitioning scheme
 int nDofs = nB; // this is for row partitioning
 if( !row_partition ) nDofs = nA; // column partitioning
 
 // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
 ownership.resize( size );
 int nPerPart = nDofs / size;
 int nRemainder = nDofs % size; // Keep the remainder in root
 ownership[0] = nPerPart + nRemainder;
 for( int ip = 1, roffset = ownership[0]; ip < size; ++ip )
 {
 roffset += nPerPart;
 ownership[ip] = roffset;
 }
 moab::TupleList* tl = new moab::TupleList;
 unsigned numr = 1; //
 tl->initialize( 3, 0, 0, numr, localSize ); // to proc, row, col, value
 tl->enableWriteAccess();
 // populate
 for( int i = 0; i < localSize; i++ )
 {
 int rowval = vecRow[i] - 1; // dofs are 1 based in the file
 int colval = vecCol[i] - 1;
 int to_proc = -1;
 int dof_val = colval;
 if( row_partition ) dof_val = rowval;
 
 if( ownership[0] > dof_val )
 to_proc = 0;
 else
 {
 for( int ip = 1; ip < size; ++ip )
 {
 if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
 {
 to_proc = ip;
 break;
 }
 }
 }
 
 int n = tl->get_n();
 tl->vi_wr[3 * n] = to_proc;
 tl->vi_wr[3 * n + 1] = rowval;
 tl->vi_wr[3 * n + 2] = colval;
 tl->vr_wr[n] = vecS[i];
 tl->inc_n();
 }
 // heavy communication
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl, 0 );
 
 if( owned_dof_ids.size() > 0 )
 {
 // we need to send desired dof to the rendez_vous point
 moab::TupleList tl_re; //
 tl_re.initialize( 2, 0, 0, 0, owned_dof_ids.size() ); // to proc, value
 tl_re.enableWriteAccess();
 // send first to rendez_vous point, decided by trivial partitioning
 
 for( size_t i = 0; i < owned_dof_ids.size(); i++ )
 {
 int to_proc = -1;
 int dof_val = owned_dof_ids[i] - 1; // dofs are 1 based in the file, partition from 0 ?
 
 if( ownership[0] > dof_val )
 to_proc = 0;
 else
 {
 for( int ip = 1; ip < size; ++ip )
 {
 if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
 {
 to_proc = ip;
 break;
 }
 }
 }
 
 int n = tl_re.get_n();
 tl_re.vi_wr[2 * n] = to_proc;
 tl_re.vi_wr[2 * n + 1] = dof_val;
 
 tl_re.inc_n();
 }
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl_re, 0 );
 // now we know in tl_re where do we need to send back dof_val
 moab::TupleList::buffer sort_buffer;
 sort_buffer.buffer_init( tl_re.get_n() );
 tl_re.sort( 1, &sort_buffer ); // so now we order by value
 
 sort_buffer.buffer_init( tl->get_n() );
 int indexOrder = 2; // colVal
 if( row_partition ) indexOrder = 1; // rowVal
 //tl->sort( indexOrder, &sort_buffer );
 
 std::map< int, int > startDofIndex, endDofIndex; // indices in tl_re for values we want
 int dofVal = -1;
 if( tl_re.get_n() > 0 ) dofVal = tl_re.vi_rd[1]; // first dof val on this rank
 startDofIndex[dofVal] = 0;
 endDofIndex[dofVal] = 0; // start and end
 for( unsigned k = 1; k < tl_re.get_n(); k++ )
 {
 int newDof = tl_re.vi_rd[2 * k + 1];
 if( dofVal == newDof )
 {
 endDofIndex[dofVal] = k; // increment by 1 actually
 }
 else
 {
 dofVal = newDof;
 startDofIndex[dofVal] = k;
 endDofIndex[dofVal] = k;
 }
 }
 
 // basically, for each value we are interested in, index in tl_re with those values are
 // tl_re.vi_rd[2*startDofIndex+1] == valDof == tl_re.vi_rd[2*endDofIndex+1]
 // so now we have ordered
 // tl_re shows to what proc do we need to send the tuple (row, col, val)
 moab::TupleList* tl_back = new moab::TupleList;
 unsigned numr = 1; //
 // localSize is a good guess, but maybe it should be bigger ?
 // this could be bigger for repeated dofs
 tl_back->initialize( 3, 0, 0, numr, tl->get_n() ); // to proc, row, col, value
 tl_back->enableWriteAccess();
 // now loop over tl and tl_re to see where to send
 // form the new tuple, which will contain the desired dofs per task, per row or column distribution
 
 for( unsigned k = 0; k < tl->get_n(); k++ )
 {
 int valDof = tl->vi_rd[3 * k + indexOrder]; // 1 for row, 2 for column // first value, it should be
 for( int ire = startDofIndex[valDof]; ire <= endDofIndex[valDof]; ire++ )
 {
 int to_proc = tl_re.vi_rd[2 * ire];
 int n = tl_back->get_n();
 tl_back->vi_wr[3 * n] = to_proc;
 tl_back->vi_wr[3 * n + 1] = tl->vi_rd[3 * k + 1]; // row
 tl_back->vi_wr[3 * n + 2] = tl->vi_rd[3 * k + 2]; // col
 tl_back->vr_wr[n] = tl->vr_rd[k];
 tl_back->inc_n();
 }
 }
 
 // now communicate to the desired tasks:
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl_back, 0 );
 
 tl_re.reset(); // clear memory, although this will go out of scope
 tl->reset();
 tl = tl_back;
 }
 
 int rindexMax = 0, cindexMax = 0;
 // populate the sparsematrix, using rowMap and colMap
 int n = tl->get_n();
 for( int i = 0; i < n; i++ )
 {
 int rindex, cindex;
 const int& vecRowValue = tl->vi_wr[3 * i + 1];
 const int& vecColValue = tl->vi_wr[3 * i + 2];
 
 std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
 if( riter == rowMap.end() )
 {
 rowMap[vecRowValue] = rindexMax;
 rindex = rindexMax;
 row_gdofmap.push_back( vecRowValue );
 // row_dtoc_dofmap.push_back( vecRowValue );
 rindexMax++;
 }
 else
 rindex = riter->second;
 
 std::map< int, int >::iterator citer = colMap.find( vecColValue );
 if( citer == colMap.end() )
 {
 colMap[vecColValue] = cindexMax;
 cindex = cindexMax;
 col_gdofmap.push_back( vecColValue );
 // col_dtoc_dofmap.push_back( vecColValue );
 cindexMax++;
 }
 else
 cindex = citer->second;
 
 sparseMatrix( rindex, cindex ) = tl->vr_wr[i];
 }
 tl->reset();
 }
 else
#endif
 {
 int rindexMax = 0, cindexMax = 0;
 
 for( int i = 0; i < nS; i++ )
 {
 int rindex, cindex;
 const int& vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
 const int& vecColValue = vecCol[i] - 1;
 
 std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
 if( riter == rowMap.end() )
 {
 rowMap[vecRowValue] = rindexMax;
 rindex = rindexMax;
 row_gdofmap.push_back( vecRowValue );
 // row_dtoc_dofmap.push_back( vecRowValue );
 rindexMax++;
 }
 else
 rindex = riter->second;
 
 std::map< int, int >::iterator citer = colMap.find( vecColValue );
 if( citer == colMap.end() )
 {
 colMap[vecColValue] = cindexMax;
 cindex = cindexMax;
 col_gdofmap.push_back( vecColValue );
 // col_dtoc_dofmap.push_back( vecColValue );
 cindexMax++;
 }
 else
 cindex = citer->second;
 
 sparseMatrix( rindex, cindex ) = vecS[i];
 }
 }
 
 m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
 m_nTotDofs_Dest = sparseMatrix.GetRows();
 
#ifdef MOAB_HAVE_EIGEN3
 this->copy_tempest_sparsemat_to_eigen3();
#endif
 
 // Reset the source and target data first
 m_rowVector.setZero();
 m_colVector.setZero();
 
 return moab::MB_SUCCESS;
}

References CHECK_EXCEPTION, moab::TupleList::enableWriteAccess(), moab::error(), moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_SUCCESS, moab::TupleList::reset(), size, moab::TupleList::sort(), moab::TupleList::vi_rd, moab::TupleList::vi_wr, moab::TupleList::vr_rd, and moab::TupleList::vr_wr.

set_col_dc_dofs()

moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs

(

std::vector< int > &

values_entities

)

Definition at line 647 of file TempestOnlineMap.cpp.

{
 // col_gdofmap has global dofs , that should be in the list of values, such that
 // row_dtoc_dofmap[offsetDOF] = localDOF;
 // // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
 // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
 // form first inverse
 
 col_dtoc_dofmap.resize( values_entities.size() );
 for( int j = 0; j < (int)values_entities.size(); j++ )
 {
 if( colMap.find( values_entities[j] - 1 ) != colMap.end() )
 col_dtoc_dofmap[j] = colMap[values_entities[j] - 1];
 else
 {
 col_dtoc_dofmap[j] = -1; // signal that this value should not be used in
 // std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
 // found in colMap \n";
 }
 }
 return moab::MB_SUCCESS;
}

References MB_SUCCESS.

set_row_dc_dofs()

moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs

(

std::vector< int > &

values_entities

)

Definition at line 670 of file TempestOnlineMap.cpp.

{
 // row_dtoc_dofmap = values_entities; // needs to point to local
 // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
 
 row_dtoc_dofmap.resize( values_entities.size() );
 for( int j = 0; j < (int)values_entities.size(); j++ )
 {
 if( rowMap.find( values_entities[j] - 1 ) != rowMap.end() )
 row_dtoc_dofmap[j] = rowMap[values_entities[j] - 1]; // values are 1 based, but rowMap, colMap are not
 else
 {
 row_dtoc_dofmap[j] = -1; // not all values are used
 // std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
 // found in rowMap \n";
 }
 }
 return moab::MB_SUCCESS;
}

References MB_SUCCESS.

SetDestinationNDofsPerElement()

void moab::TempestOnlineMap::SetDestinationNDofsPerElement ( int  nt )

inline

Get the number of Degrees-Of-Freedom per element on the destination mesh.

Definition at line 510 of file TempestOnlineMap.hpp.

{
 m_nDofsPEl_Dest = nt;
}

SetDOFmapAssociation()

moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation	(	DiscretizationType	srcType,
		bool	isSrcContinuous,
		DataArray3D< int > *	srcdataGLLNodes,
		DataArray3D< int > *	srcdataGLLNodesSrc,
		DiscretizationType	destType,
		bool	isDestContinuous,
		DataArray3D< int > *	tgtdataGLLNodes
	)

Compute the association between the solution tag global DoF numbering and the local matrix numbering so that matvec operations can be performed consistently.

Definition at line 172 of file TempestOnlineMap.cpp.

{
 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;
 
 // Compute and store the total number of source and target DoFs corresponding
 // to number of rows and columns in the mapping.
#ifdef VVERBOSE
 {
 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
 rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
 
 if( is_root )
 {
 {
 std::ofstream output_file( "sourcecov-gids-0.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
 output_file << i << ", " << src_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
 m_nTotDofs_SrcCov = 0;
 if( isSrcContinuous )
 dgll_cgll_covcol_ldofmap.resize(
 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
 {
 m_nTotDofs_SrcCov++;
 dgll_cgll_covcol_ldofmap[localDOF] = true;
 }
 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
 << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
 }
 }
 }
 output_file.flush(); // required here
 output_file.close();
 dgll_cgll_covcol_ldofmap.clear();
 }
 
 {
 std::ofstream output_file( "source-gids-0.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
 m_nTotDofs_Src = 0;
 if( isSrcContinuous )
 dgll_cgll_col_ldofmap.resize(
 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
 {
 m_nTotDofs_Src++;
 dgll_cgll_col_ldofmap[localDOF] = true;
 }
 output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
 << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
 }
 }
 }
 output_file.flush(); // required here
 output_file.close();
 dgll_cgll_col_ldofmap.clear();
 }
 
 {
 std::ofstream output_file( "target-gids-0.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
 m_nTotDofs_Dest = 0;
 
 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
 {
 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
 << m_nTotDofs_Dest << "\n";
 m_nTotDofs_Dest++;
 }
 
 output_file.flush(); // required here
 output_file.close();
 }
 }
 else
 {
 {
 std::ofstream output_file( "sourcecov-gids-1.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
 output_file << i << ", " << src_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
 m_nTotDofs_SrcCov = 0;
 if( isSrcContinuous )
 dgll_cgll_covcol_ldofmap.resize(
 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
 {
 m_nTotDofs_SrcCov++;
 dgll_cgll_covcol_ldofmap[localDOF] = true;
 }
 output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
 << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
 }
 }
 }
 output_file.flush(); // required here
 output_file.close();
 dgll_cgll_covcol_ldofmap.clear();
 }
 
 {
 std::ofstream output_file( "source-gids-1.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
 output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
 m_nTotDofs_Src = 0;
 if( isSrcContinuous )
 dgll_cgll_col_ldofmap.resize(
 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
 {
 m_nTotDofs_Src++;
 dgll_cgll_col_ldofmap[localDOF] = true;
 }
 output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
 << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
 }
 }
 }
 output_file.flush(); // required here
 output_file.close();
 dgll_cgll_col_ldofmap.clear();
 }
 
 {
 std::ofstream output_file( "target-gids-1.txt" );
 output_file << "I, GDOF\n";
 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
 output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
 
 output_file << "ELEMID, IDOF, GDOF, NDOF\n";
 m_nTotDofs_Dest = 0;
 
 for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
 {
 output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
 << m_nTotDofs_Dest << "\n";
 m_nTotDofs_Dest++;
 }
 
 output_file.flush(); // required here
 output_file.close();
 }
 }
 }
#endif
 
 // Now compute the mapping and store it for the covering mesh
 int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
 if( m_remapper->point_cloud_source )
 {
 assert( m_nDofsPEl_Src == 1 );
 col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
 col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
 src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
 srcTagSize = 1;
 }
 else
 {
 col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
 col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
 }
 
 // std::cout << "TOnlineMap: Process: " << rank << " and covering entities = [" <<
 // col_dofmap.size() << ", " << src_soln_gdofs.size() << "]\n"; MPI_Barrier(MPI_COMM_WORLD);
 
#ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
 unsigned maxSrcIndx = 0;
 
 // for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
 std::vector< int > locdofs( srcTagSize );
 std::map< Node, moab::EntityHandle > mapLocalMBNodes;
 double elcoords[3];
 for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
 {
 EntityHandle eh = m_remapper->m_covering_source_entities[iel];
 rval = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
 Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
 mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
 }
 
 const NodeVector& nodes = m_remapper->m_covering_source->nodes;
 for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
 {
 const Face& face = m_remapper->m_covering_source->faces[j];
 
 Node centroid;
 centroid.x = centroid.y = centroid.z = 0.0;
 for( unsigned l = 0; l < face.edges.size(); ++l )
 {
 centroid.x += nodes[face[l]].x;
 centroid.y += nodes[face[l]].y;
 centroid.z += nodes[face[l]].z;
 }
 const double factor = 1.0 / face.edges.size();
 centroid.x *= factor;
 centroid.y *= factor;
 centroid.z *= factor;
 
 EntityHandle current_eh;
 if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
 {
 current_eh = mapLocalMBNodes[centroid];
 }
 
 rval = m_interface->tag_get_data( m_dofTagSrc, &current_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
 const int offsetDOF = p * m_nDofsPEl_Src + q;
 maxSrcIndx = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
 std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
 << ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
 }
 }
 }
#endif
 
 m_nTotDofs_SrcCov = 0;
 if( srcdataGLLNodes == NULL )
 { /* we only have a mapping for elements as DoFs */
 for( unsigned i = 0; i < col_gdofmap.size(); ++i )
 {
 assert( src_soln_gdofs[i] > 0 );
 col_gdofmap[i] = src_soln_gdofs[i] - 1;
 col_dtoc_dofmap[i] = i;
 if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
 m_nTotDofs_SrcCov++;
 }
 }
 else
 {
 if( isSrcContinuous )
 dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
 // Put these remap coefficients into the SparseMatrix map
 for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
 {
 m_nTotDofs_SrcCov++;
 dgll_cgll_covcol_ldofmap[localDOF] = true;
 }
 if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
 assert( src_soln_gdofs[offsetDOF] > 0 );
 col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
 col_dtoc_dofmap[offsetDOF] = localDOF;
 if( vprint )
 std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", "
 << localDOF << ", " << col_gdofmap[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
 }
 }
 }
 }
 
 if( m_remapper->point_cloud_source )
 {
 assert( m_nDofsPEl_Src == 1 );
 srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
 srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
 locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
 }
 else
 {
 srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
 srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
 }
 
 // Now compute the mapping and store it for the original source mesh
 m_nTotDofs_Src = 0;
 if( srcdataGLLNodesSrc == NULL )
 { /* we only have a mapping for elements as DoFs */
 for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
 {
 assert( locsrc_soln_gdofs[i] > 0 );
 srccol_gdofmap[i] = locsrc_soln_gdofs[i] - 1;
 srccol_dtoc_dofmap[i] = i;
 m_nTotDofs_Src++;
 }
 }
 else
 {
 if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
 // Put these remap coefficients into the SparseMatrix map
 for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Src; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Src; q++ )
 {
 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
 {
 m_nTotDofs_Src++;
 dgll_cgll_col_ldofmap[localDOF] = true;
 }
 if( !isSrcContinuous ) m_nTotDofs_Src++;
 assert( locsrc_soln_gdofs[offsetDOF] > 0 );
 srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
 srccol_dtoc_dofmap[offsetDOF] = localDOF;
 }
 }
 }
 }
 
 int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
 if( m_remapper->point_cloud_target )
 {
 assert( m_nDofsPEl_Dest == 1 );
 row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
 row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
 tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
 tgtTagSize = 1;
 }
 else
 {
 row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
 row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
 rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
 }
 
 // Now compute the mapping and store it for the target mesh
 // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
 m_nTotDofs_Dest = 0;
 if( tgtdataGLLNodes == NULL )
 { /* we only have a mapping for elements as DoFs */
 for( unsigned i = 0; i < row_gdofmap.size(); ++i )
 {
 assert( tgt_soln_gdofs[i] > 0 );
 row_gdofmap[i] = tgt_soln_gdofs[i] - 1;
 row_dtoc_dofmap[i] = i;
 if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
 m_nTotDofs_Dest++;
 }
 }
 else
 {
 if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
 // Put these remap coefficients into the SparseMatrix map
 for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
 {
 for( int p = 0; p < m_nDofsPEl_Dest; p++ )
 {
 for( int q = 0; q < m_nDofsPEl_Dest; q++ )
 {
 const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
 const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
 if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
 {
 m_nTotDofs_Dest++;
 dgll_cgll_row_ldofmap[localDOF] = true;
 }
 if( !isTgtContinuous ) m_nTotDofs_Dest++;
 assert( tgt_soln_gdofs[offsetDOF] > 0 );
 row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
 row_dtoc_dofmap[offsetDOF] = localDOF;
 if( vprint )
 std::cout << "Row: " << m_remapper->lid_to_gid_tgt[j] << ", " << offsetDOF << ", " << localDOF
 << ", " << row_gdofmap[offsetDOF] << ", " << m_nTotDofs_Dest << "\n";
 }
 }
 }
 }
 
 // Let us also allocate the local representation of the sparse matrix
#if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
 if( vprint )
 {
 std::cout << "[" << rank << "]"
 << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size() << ", col = " << m_nTotDofs_Src
 << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
 // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
 }
#endif
 
 // check monotonicity of row_gdofmap and col_gdofmap
#ifdef CHECK_INCREASING_DOF
 for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
 {
 if( row_gdofmap[i] > row_gdofmap[i + 1] )
 std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
 << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
 }
 for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
 {
 if( col_gdofmap[i] > col_gdofmap[i + 1] )
 std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
 << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
 }
#endif
 
 return moab::MB_SUCCESS;
}

References ErrorCode, MB_CHK_ERR, and MB_SUCCESS.

SetDOFmapTags()

moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags	(	const std::string	srcDofTagName,
		const std::string	tgtDofTagName
	)

Store the tag names associated with global DoF ids for source and target meshes

Definition at line 116 of file TempestOnlineMap.cpp.

{
 moab::ErrorCode rval;
 
 int tagSize = 0;
 tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
 rval =
 m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
 
 if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
 {
 int ntot_elements = 0, nelements = m_remapper->m_source_entities.size();
#ifdef MOAB_HAVE_MPI
 int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
 if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
#else
 ntot_elements = nelements;
#endif
 
 rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_covering_source_entities,
 &m_remapper->m_source_entities, srcDofTagName, m_nDofsPEl_Src );MB_CHK_ERR( rval );
 
 rval = m_interface->tag_get_handle( srcDofTagName.c_str(), m_nDofsPEl_Src * m_nDofsPEl_Src, MB_TYPE_INTEGER,
 this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR( rval );
 }
 else
 MB_CHK_ERR( rval );
 
 tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
 rval =
 m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
 if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
 {
 int ntot_elements = 0, nelements = m_remapper->m_target_entities.size();
#ifdef MOAB_HAVE_MPI
 int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
 if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
#else
 ntot_elements = nelements;
#endif
 
 rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_target_entities, NULL, tgtDofTagName,
 m_nDofsPEl_Dest );MB_CHK_ERR( rval );
 
 rval = m_interface->tag_get_handle( tgtDofTagName.c_str(), m_nDofsPEl_Dest * m_nDofsPEl_Dest, MB_TYPE_INTEGER,
 this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR( rval );
 }
 else
 MB_CHK_ERR( rval );
 
 return moab::MB_SUCCESS;
}

References ErrorCode, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_ANY, MB_TAG_NOT_FOUND, and MB_TYPE_INTEGER.

SetSourceNDofsPerElement()

void moab::TempestOnlineMap::SetSourceNDofsPerElement ( int  ns )

inline

Set the number of Degrees-Of-Freedom per element on the source mesh.

Definition at line 505 of file TempestOnlineMap.hpp.

{
 m_nDofsPEl_Src = ns;
}

setup_sizes_dimensions()

void moab::TempestOnlineMap::setup_sizes_dimensions ( )

private

Definition at line 78 of file TempestOnlineMap.cpp.

{
 if( m_meshInputCov )
 {
 std::vector< std::string > dimNames;
 std::vector< int > dimSizes;
 dimNames.push_back( "num_elem" );
 dimSizes.push_back( m_meshInputCov->faces.size() );
 
 this->InitializeSourceDimensions( dimNames, dimSizes );
 }
 
 if( m_meshOutput )
 {
 std::vector< std::string > dimNames;
 std::vector< int > dimSizes;
 dimNames.push_back( "num_elem" );
 dimSizes.push_back( m_meshOutput->faces.size() );
 
 this->InitializeTargetDimensions( dimNames, dimSizes );
 }
}

Referenced by TempestOnlineMap().

WriteHDF5MapFile()

moab::ErrorCode moab::TempestOnlineMap::WriteHDF5MapFile ( const std::string &  filename )

private

Parallel I/O with NetCDF to write out the SCRIP file from multiple processors.

Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.

Definition at line 790 of file TempestOnlineMapIO.cpp.

{
 moab::ErrorCode rval;
 
 /**
 * Need to get the global maximum of number of vertices per element
 * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
 *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
 *when writing this out, other processes may have a different size for the same array. This is
 *hence a mess to consolidate in h5mtoscrip eventually.
 **/
 
 /* Let us compute all relevant data for the current original source mesh on the process */
 DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
 DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
 DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
 if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
 {
 this->InitializeCoordinatesFromMeshFV(
 *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
 ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
 m_remapper->max_source_edges );
 
 vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
 for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
 vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
 }
 else
 {
 DataArray3D< double > dataGLLJacobianSrc;
 this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
 dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
 
 // Generate the continuous Jacobian for input mesh
 GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
 
 if( m_srcDiscType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, m_meshInput->vecFaceArea );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianSrc, m_meshInput->vecFaceArea );
 }
 
 vecSourceFaceArea.Allocate( m_meshInput->faces.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
 int offset = 0;
 for( size_t e = 0; e < m_meshInput->faces.size(); e++ )
 {
 for( int s = 0; s < m_nDofsPEl_Src; s++ )
 {
 for( int t = 0; t < m_nDofsPEl_Src; t++ )
 {
 vecSourceFaceArea[srccol_dtoc_dofmap[offset + s * m_nDofsPEl_Src + t]] =
 dataGLLJacobianSrc[s][t][e];
 }
 }
 offset += m_nDofsPEl_Src * m_nDofsPEl_Src;
 }
 }
 
 if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
 {
 this->InitializeCoordinatesFromMeshFV(
 *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
 ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
 m_remapper->max_target_edges );
 
 vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
 for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
 vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
 }
 else
 {
 DataArray3D< double > dataGLLJacobianDest;
 this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
 dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
 
 // Generate the continuous Jacobian for input mesh
 GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
 
 if( m_destDiscType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, m_meshOutput->vecFaceArea );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianDest, m_meshOutput->vecFaceArea );
 }
 
 vecTargetFaceArea.Allocate( m_meshOutput->faces.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
 int offset = 0;
 for( size_t e = 0; e < m_meshOutput->faces.size(); e++ )
 {
 for( int s = 0; s < m_nDofsPEl_Dest; s++ )
 {
 for( int t = 0; t < m_nDofsPEl_Dest; t++ )
 {
 vecTargetFaceArea[row_dtoc_dofmap[offset + s * m_nDofsPEl_Dest + t]] = dataGLLJacobianDest[s][t][e];
 }
 }
 offset += m_nDofsPEl_Dest * m_nDofsPEl_Dest;
 }
 }
 
 moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
 int tot_src_ents = m_remapper->m_source_entities.size();
 int tot_tgt_ents = m_remapper->m_target_entities.size();
 int tot_src_size = dSourceCenterLon.GetRows();
 int tot_tgt_size = m_dTargetCenterLon.GetRows();
 int tot_vsrc_size = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
 int tot_vtgt_size = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
 
 const int weightMatNNZ = m_weightMatrix.nonZeros();
 moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
 rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 moab::Tag srcAreaValues, tgtAreaValues;
 rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
 rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
 rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 moab::Tag srcMaskValues, tgtMaskValues;
 if( m_iSourceMask.IsAttached() )
 {
 rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 }
 if( m_iTargetMask.IsAttached() )
 {
 rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
 moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
 }
 
 std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
 std::vector< double > smatvals( weightMatNNZ );
 // const double* smatvals = m_weightMatrix.valuePtr();
 // Loop over the matrix entries and find the max global ID for rows and columns
 for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
 {
 for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
 {
 smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
 smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
 smatvals[offset] = it.value();
 }
 }
 
 /* Set the global IDs for the DoFs */
 ////
 // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
 // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
 ////
 int maxrow = 0, maxcol = 0;
 std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
 for( int i = 0; i < tot_src_size; ++i )
 {
 src_global_dofs[i] = srccol_gdofmap[i];
 maxcol = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
 }
 
 for( int i = 0; i < tot_tgt_size; ++i )
 {
 tgt_global_dofs[i] = row_gdofmap[i];
 maxrow = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
 }
 
 ///////////////////////////////////////////////////////////////////////////
 // The metadata in H5M file contains the following data:
 //
 // 1. n_a: Total source entities: (number of elements in source mesh)
 // 2. n_b: Total target entities: (number of elements in target mesh)
 // 3. nv_a: Max edge size of elements in source mesh
 // 4. nv_b: Max edge size of elements in target mesh
 // 5. maxrows: Number of rows in remap weight matrix
 // 6. maxcols: Number of cols in remap weight matrix
 // 7. nnz: Number of total nnz in sparse remap weight matrix
 // 8. np_a: The order of the field description on the source mesh: >= 1
 // 9. np_b: The order of the field description on the target mesh: >= 1
 // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
 // dGLL]
 // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
 // dGLL]
 // 12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
 // 1]
 // 13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
 // [0, 1, 2]
 //
 ///////////////////////////////////////////////////////////////////////////
 int map_disc_details[6];
 map_disc_details[0] = m_nDofsPEl_Src;
 map_disc_details[1] = m_nDofsPEl_Dest;
 map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
 ? 0
 : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
 map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
 ? 0
 : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
 map_disc_details[4] = ( m_bConserved ? 1 : 0 );
 map_disc_details[5] = m_iMonotonicity;
 
#ifdef MOAB_HAVE_MPI
 int loc_smatmetadata[13] = { tot_src_ents,
 tot_tgt_ents,
 m_remapper->max_source_edges,
 m_remapper->max_target_edges,
 maxrow + 1,
 maxcol + 1,
 weightMatNNZ,
 map_disc_details[0],
 map_disc_details[1],
 map_disc_details[2],
 map_disc_details[3],
 map_disc_details[4],
 map_disc_details[5] };
 rval = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 int glb_smatmetadata[13] = { 0,
 0,
 0,
 0,
 0,
 0,
 0,
 map_disc_details[0],
 map_disc_details[1],
 map_disc_details[2],
 map_disc_details[3],
 map_disc_details[4],
 map_disc_details[5] };
 int loc_buf[7] = {
 tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
 maxrow, maxcol };
 int glb_buf[4] = { 0, 0, 0, 0 };
 MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
 glb_smatmetadata[0] = glb_buf[0];
 glb_smatmetadata[1] = glb_buf[1];
 glb_smatmetadata[6] = glb_buf[2];
 MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
 glb_smatmetadata[2] = glb_buf[0];
 glb_smatmetadata[3] = glb_buf[1];
 glb_smatmetadata[4] = glb_buf[2];
 glb_smatmetadata[5] = glb_buf[3];
#else
 int glb_smatmetadata[13] = { tot_src_ents,
 tot_tgt_ents,
 m_remapper->max_source_edges,
 m_remapper->max_target_edges,
 maxrow,
 maxcol,
 weightMatNNZ,
 map_disc_details[0],
 map_disc_details[1],
 map_disc_details[2],
 map_disc_details[3],
 map_disc_details[4],
 map_disc_details[5] };
#endif
 // These values represent number of rows and columns. So should be 1-based.
 glb_smatmetadata[4]++;
 glb_smatmetadata[5]++;
 
 if( this->is_root )
 {
 std::cout << " " << this->rank << " Writing remap weights with size [" << glb_smatmetadata[4] << " X "
 << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
 EntityHandle root_set = 0;
 rval = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 }
 
 int dsize;
 const int numval = weightMatNNZ;
 const void* smatrowvals_d = smatrowvals.data();
 const void* smatcolvals_d = smatcolvals.data();
 const void* smatvals_d = smatvals.data();
 rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 /* Set the global IDs for the DoFs */
 const void* srceleidvals_d = src_global_dofs.data();
 const void* tgteleidvals_d = tgt_global_dofs.data();
 dsize = src_global_dofs.size();
 rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 dsize = tgt_global_dofs.size();
 rval = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 /* Set the source and target areas */
 const void* srcareavals_d = vecSourceFaceArea;
 const void* tgtareavals_d = vecTargetFaceArea;
 dsize = tot_src_size;
 rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 dsize = tot_tgt_size;
 rval = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 /* Set the coordinates for source and target center vertices */
 const void* srccoordsclonvals_d = &dSourceCenterLon[0];
 const void* srccoordsclatvals_d = &dSourceCenterLat[0];
 dsize = dSourceCenterLon.GetRows();
 rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
 const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
 dsize = vecTargetFaceArea.GetRows();
 rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 /* Set the coordinates for source and target element vertices */
 const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
 const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
 dsize = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
 rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
 const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
 dsize = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
 rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 
 /* Set the masks for source and target meshes if available */
 if( m_iSourceMask.IsAttached() )
 {
 const void* srcmaskvals_d = m_iSourceMask;
 dsize = m_iSourceMask.GetRows();
 rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 }
 
 if( m_iTargetMask.IsAttached() )
 {
 const void* tgtmaskvals_d = m_iTargetMask;
 dsize = m_iTargetMask.GetRows();
 rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
 }
 
#ifdef MOAB_HAVE_MPI
 const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
 const char* writeOptions = "";
#endif
 
 // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
 EntityHandle sets[1] = { m_remapper->m_overlap_set };
 rval = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
 
#ifdef WRITE_SCRIP_FILE
 sstr.str( "" );
 sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
 std::map< std::string, std::string > mapAttributes;
 mapAttributes["Creator"] = "MOAB mbtempest workflow";
 if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
 this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
 sstr.str( "" );
#endif
 
 return moab::MB_SUCCESS;
}

References ErrorCode, MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_SPARSE, MB_TAG_VARLEN, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, moab::TempestRemapper::RLL, size, and t.

WriteParallelMap()

moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap	(	const std::string &	strTarget,
		const std::map< std::string, std::string > &	attrMap
	)

Write the TempestOnlineMap to a parallel NetCDF file.

Definition at line 156 of file TempestOnlineMapIO.cpp.

{
 moab::ErrorCode rval;
 
 size_t lastindex = strFilename.find_last_of( "." );
 std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
 
 // Write the map file to disk in parallel
 if( extension == "nc" )
 {
 /* Invoke the actual call to write the parallel map to disk in SCRIP format */
 rval = this->WriteSCRIPMapFile( strFilename.c_str(), attrMap );MB_CHK_ERR( rval );
 }
 else
 {
 /* Write to the parallel H5M format */
 rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
 }
 
 return rval;
}

References ErrorCode, and MB_CHK_ERR.

Referenced by main().

WriteSCRIPMapFile()

moab::ErrorCode moab::TempestOnlineMap::WriteSCRIPMapFile ( const std::string &  strOutputFile, const std::map< std::string, std::string > &  attrMap  )

private

Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections.

Parallel I/O with HDF5 to write out the remapping weights from multiple processors.

Need to get the global maximum of number of vertices per element Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However, when writing this out, other processes may have a different size for the same array. This is hence a mess to consolidate in h5mtoscrip eventually.

Definition at line 181 of file TempestOnlineMapIO.cpp.

{
 NcError error( NcError::silent_nonfatal );
 
#ifdef MOAB_HAVE_NETCDFPAR
 bool is_independent = true;
 ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
 // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
 NcFile ncMap( strFilename.c_str(), NcFile::Replace );
#endif
 
 if( !ncMap.is_valid() )
 {
 _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
 }
 
 // Attributes
 // ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
 auto it = attrMap.begin();
 while( it != attrMap.end() )
 {
 // set the map attributes
 ncMap.add_att( it->first.c_str(), it->second.c_str() );
 // increment iterator
 it++;
 }
 
 /**
 * Need to get the global maximum of number of vertices per element
 * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
 *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
 *when writing this out, other processes may have a different size for the same array. This is
 *hence a mess to consolidate in h5mtoscrip eventually.
 **/
 
 /* Let us compute all relevant data for the current original source mesh on the process */
 DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
 DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
 DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
 if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
 {
 this->InitializeCoordinatesFromMeshFV(
 *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
 ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
 m_remapper->max_source_edges );
 
 vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
 for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
 vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
 }
 else
 {
 DataArray3D< double > dataGLLJacobianSrc;
 this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
 dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
 
 // Generate the continuous Jacobian for input mesh
 GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
 
 if( m_srcDiscType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
 }
 }
 
 if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
 {
 this->InitializeCoordinatesFromMeshFV(
 *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
 ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
 m_remapper->max_target_edges );
 
 vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
 for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
 {
 vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
 }
 }
 else
 {
 DataArray3D< double > dataGLLJacobianDest;
 this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
 dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
 
 // Generate the continuous Jacobian for input mesh
 GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
 
 if( m_destDiscType == DiscretizationType_CGLL )
 {
 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
 }
 else
 {
 GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
 }
 }
 
 // Map dimensions
 unsigned nA = ( vecSourceFaceArea.GetRows() );
 unsigned nB = ( vecTargetFaceArea.GetRows() );
 
 std::vector< int > masksA, masksB;
 ErrorCode rval = m_remapper->GetIMasks( moab::Remapper::SourceMesh, masksA );MB_CHK_SET_ERR( rval, "Trouble getting masks for source" );
 rval = m_remapper->GetIMasks( moab::Remapper::TargetMesh, masksB );MB_CHK_SET_ERR( rval, "Trouble getting masks for target" );
 
 // Number of nodes per Face
 int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
 int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
 
 // if source or target cells have triangles at poles, center of those triangles need to come from
 // the original quad, not from center in 3d, converted to 2d again
 // start copy OnlineMap.cpp tempestremap
 // right now, do this only for source mesh; copy the logic for target mesh
 for( unsigned i = 0; i < nA; i++ )
 {
 const Face& face = m_meshInput->faces[i];
 
 int nNodes = face.edges.size();
 int indexNodeAtPole = -1;
 if( 3 == nNodes ) // check if one node at the poles
 {
 for( int j = 0; j < nNodes; j++ )
 if( fabs( fabs( dSourceVertexLat[i][j] ) - 90.0 ) < 1.0e-12 )
 {
 indexNodeAtPole = j;
 break;
 }
 }
 if( indexNodeAtPole < 0 ) continue; // continue i loop, do nothing
 // recompute center of cell, from 3d data; add one 2 nodes at pole, and average
 int nodeAtPole = face[indexNodeAtPole]; // use the overloaded operator
 Node nodePole = m_meshInput->nodes[nodeAtPole];
 Node newCenter = nodePole * 2;
 for( int j = 1; j < nNodes; j++ )
 {
 int indexi = ( indexNodeAtPole + j ) % nNodes; // nNodes is 3 !
 const Node& node = m_meshInput->nodes[face[indexi]];
 newCenter = newCenter + node;
 }
 newCenter = newCenter * 0.25;
 newCenter = newCenter.Normalized();
 
#ifdef VERBOSE
 double iniLon = dSourceCenterLon[i], iniLat = dSourceCenterLat[i];
#endif
 // dSourceCenterLon, dSourceCenterLat
 XYZtoRLL_Deg( newCenter.x, newCenter.y, newCenter.z, dSourceCenterLon[i], dSourceCenterLat[i] );
#ifdef VERBOSE
 std::cout << " modify center of triangle from " << iniLon << " " << iniLat << " to " << dSourceCenterLon[i]
 << " " << dSourceCenterLat[i] << "\n";
#endif
 }
 
 // first move data if in parallel
#if defined( MOAB_HAVE_MPI )
 int max_row_dof, max_col_dof; // output; arrays will be re-distributed in chunks [maxdof/size]
 // if (size > 1)
 {
 int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
 dSourceVertexLon, dSourceVertexLat, masksA, nA, nSourceNodesPerFace,
 max_col_dof ); // now nA will be close to maxdof/size
 if( ierr != 0 )
 {
 _EXCEPTION1( "Unable to arrange source data %d ", nA );
 }
 // rearrange target data: (nB)
 //
 ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
 dTargetVertexLon, dTargetVertexLat, masksB, nB, nTargetNodesPerFace,
 max_row_dof ); // now nA will be close to maxdof/size
 if( ierr != 0 )
 {
 _EXCEPTION1( "Unable to arrange target data %d ", nB );
 }
 }
#endif
 
 // Number of non-zeros in the remap matrix operator
 int nS = m_weightMatrix.nonZeros();
 
#if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
 int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
 int offbuf[3] = { 0, 0, 0 };
 int globuf[5] = { 0, 0, 0, 0, 0 };
 MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
 MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
 MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
 
 // MPI_Scan is inclusive of data in current rank; modify accordingly.
 offbuf[0] -= nA;
 offbuf[1] -= nB;
 offbuf[2] -= nS;
 
#else
 int offbuf[3] = { 0, 0, 0 };
 int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
#endif
 
 std::vector< std::string > srcdimNames, tgtdimNames;
 std::vector< int > srcdimSizes, tgtdimSizes;
 {
 if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
 {
 srcdimNames.push_back( "lat" );
 srcdimNames.push_back( "lon" );
 srcdimSizes.resize( 2, 0 );
 srcdimSizes[0] = m_remapper->m_source_metadata[0];
 srcdimSizes[1] = m_remapper->m_source_metadata[1];
 }
 else
 {
 srcdimNames.push_back( "num_elem" );
 srcdimSizes.push_back( globuf[0] );
 }
 
 if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
 {
 tgtdimNames.push_back( "lat" );
 tgtdimNames.push_back( "lon" );
 tgtdimSizes.resize( 2, 0 );
 tgtdimSizes[0] = m_remapper->m_target_metadata[0];
 tgtdimSizes[1] = m_remapper->m_target_metadata[1];
 }
 else
 {
 tgtdimNames.push_back( "num_elem" );
 tgtdimSizes.push_back( globuf[1] );
 }
 }
 
 // Write output dimensions entries
 unsigned nSrcGridDims = ( srcdimSizes.size() );
 unsigned nDstGridDims = ( tgtdimSizes.size() );
 
 NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
 NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
 
 NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
 NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
 
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varSrcGridDims, is_independent );
 ncMap.enable_var_par_access( varDstGridDims, is_independent );
#endif
 
 // write dimension names
 {
 char szDim[64];
 for( unsigned i = 0; i < srcdimSizes.size(); i++ )
 {
 varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
 varSrcGridDims->put( &( srcdimSizes[nSrcGridDims - i - 1] ), 1 );
 }
 
 for( unsigned i = 0; i < srcdimSizes.size(); i++ )
 {
 snprintf( szDim, 64, "name%i", i );
 varSrcGridDims->add_att( szDim, srcdimNames[nSrcGridDims - i - 1].c_str() );
 }
 
 for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
 {
 varDstGridDims->set_cur( nDstGridDims - i - 1 );
 varDstGridDims->put( &( tgtdimSizes[nDstGridDims - i - 1] ), 1 );
 }
 
 for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
 {
 snprintf( szDim, 64, "name%i", i );
 varDstGridDims->add_att( szDim, tgtdimNames[nDstGridDims - i - 1].c_str() );
 }
 }
 
 // Source and Target mesh resolutions
 NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
 NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
 
 // Number of nodes per Face
 NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
 NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
 
 // Write coordinates
 NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
 NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
 
 NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
 NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
 
 NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
 NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
 
 NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
 NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
 
 // Write masks
 NcVar* varMaskA = ncMap.add_var( "mask_a", ncInt, dimNA );
 NcVar* varMaskB = ncMap.add_var( "mask_b", ncInt, dimNB );
 
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varYCA, is_independent );
 ncMap.enable_var_par_access( varYCB, is_independent );
 ncMap.enable_var_par_access( varXCA, is_independent );
 ncMap.enable_var_par_access( varXCB, is_independent );
 ncMap.enable_var_par_access( varYVA, is_independent );
 ncMap.enable_var_par_access( varYVB, is_independent );
 ncMap.enable_var_par_access( varXVA, is_independent );
 ncMap.enable_var_par_access( varXVB, is_independent );
 ncMap.enable_var_par_access( varMaskA, is_independent );
 ncMap.enable_var_par_access( varMaskB, is_independent );
#endif
 
 varYCA->add_att( "units", "degrees" );
 varYCB->add_att( "units", "degrees" );
 
 varXCA->add_att( "units", "degrees" );
 varXCB->add_att( "units", "degrees" );
 
 varYVA->add_att( "units", "degrees" );
 varYVB->add_att( "units", "degrees" );
 
 varXVA->add_att( "units", "degrees" );
 varXVB->add_att( "units", "degrees" );
 
 // Verify dimensionality
 if( dSourceCenterLon.GetRows() != nA )
 {
 _EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" );
 }
 if( dSourceCenterLat.GetRows() != nA )
 {
 _EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" );
 }
 if( dTargetCenterLon.GetRows() != nB )
 {
 _EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" );
 }
 if( dTargetCenterLat.GetRows() != nB )
 {
 _EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" );
 }
 if( dSourceVertexLon.GetRows() != nA )
 {
 _EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" );
 }
 if( dSourceVertexLat.GetRows() != nA )
 {
 _EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" );
 }
 if( dTargetVertexLon.GetRows() != nB )
 {
 _EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" );
 }
 if( dTargetVertexLat.GetRows() != nB )
 {
 _EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" );
 }
 
 varYCA->set_cur( (long)offbuf[0] );
 varYCA->put( &( dSourceCenterLat[0] ), nA );
 varYCB->set_cur( (long)offbuf[1] );
 varYCB->put( &( dTargetCenterLat[0] ), nB );
 
 varXCA->set_cur( (long)offbuf[0] );
 varXCA->put( &( dSourceCenterLon[0] ), nA );
 varXCB->set_cur( (long)offbuf[1] );
 varXCB->put( &( dTargetCenterLon[0] ), nB );
 
 varYVA->set_cur( (long)offbuf[0] );
 varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
 varYVB->set_cur( (long)offbuf[1] );
 varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
 
 varXVA->set_cur( (long)offbuf[0] );
 varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
 varXVB->set_cur( (long)offbuf[1] );
 varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
 
 varMaskA->set_cur( (long)offbuf[0] );
 varMaskA->put( &( masksA[0] ), nA );
 varMaskB->set_cur( (long)offbuf[1] );
 varMaskB->put( &( masksB[0] ), nB );
 
 // Write areas
 NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varAreaA, is_independent );
#endif
 varAreaA->set_cur( (long)offbuf[0] );
 varAreaA->put( &( vecSourceFaceArea[0] ), nA );
 
 NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varAreaB, is_independent );
#endif
 varAreaB->set_cur( (long)offbuf[1] );
 varAreaB->put( &( vecTargetFaceArea[0] ), nB );
 
 // Write SparseMatrix entries
 DataArray1D< int > vecRow( nS );
 DataArray1D< int > vecCol( nS );
 DataArray1D< double > vecS( nS );
 DataArray1D< double > dFracA( nA );
 DataArray1D< double > dFracB( nB );
 
 moab::TupleList tlValRow, tlValCol;
 unsigned numr = 1; //
 // value has to be sent to processor row/nB for for fracA and col/nA for fracB
 // vecTargetArea (indexRow ) has to be sent for fracA (index col?)
 // vecTargetFaceArea will have to be sent to col index, with its index !
 tlValRow.initialize( 2, 0, 0, numr, nS ); // to proc(row), global row , value
 tlValCol.initialize( 3, 0, 0, numr, nS ); // to proc(col), global row / col, value
 tlValRow.enableWriteAccess();
 tlValCol.enableWriteAccess();
 /*
 *
 dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
 dFracB[ row ] += val ;
 */
 int offset = 0;
#if defined( MOAB_HAVE_MPI )
 int nAbase = ( max_col_dof + 1 ) / size; // it is nA, except last rank ( == size - 1 )
 int nBbase = ( max_row_dof + 1 ) / size; // it is nB, except last rank ( == size - 1 )
#endif
 for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
 {
 for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
 {
 vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() ); // row index
 vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() ); // col index
 vecS[offset] = it.value(); // value
 
#if defined( MOAB_HAVE_MPI )
 {
 // value M(row, col) will contribute to procRow and procCol values for fracA and fracB
 int procRow = ( vecRow[offset] - 1 ) / nBbase;
 if( procRow >= size ) procRow = size - 1;
 int procCol = ( vecCol[offset] - 1 ) / nAbase;
 if( procCol >= size ) procCol = size - 1;
 int nrInd = tlValRow.get_n();
 tlValRow.vi_wr[2 * nrInd] = procRow;
 tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
 tlValRow.vr_wr[nrInd] = vecS[offset];
 tlValRow.inc_n();
 int ncInd = tlValCol.get_n();
 tlValCol.vi_wr[3 * ncInd] = procCol;
 tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
 tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1; // this is column
 tlValCol.vr_wr[ncInd] = vecS[offset];
 tlValCol.inc_n();
 }
 
#endif
 offset++;
 }
 }
#if defined( MOAB_HAVE_MPI )
 // need to send values for their row and col processors, to compute fractions there
 // now do the heavy communication
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );
 
 // we have now, for example, dFracB[ row ] += val ;
 // so we know that on current task, we received tlValRow
 // reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
 // dFracB[ row ] += val ;
 for( unsigned i = 0; i < tlValRow.get_n(); i++ )
 {
 // int fromProc = tlValRow.vi_wr[2 * i];
 int gRowInd = tlValRow.vi_wr[2 * i + 1];
 int localIndexRow = gRowInd - nBbase * rank; // modulo nBbase rank is from 0 to size - 1;
 double wgt = tlValRow.vr_wr[i];
 assert( localIndexRow >= 0 );
 assert( nB - localIndexRow > 0 );
 dFracB[localIndexRow] += wgt;
 }
 // to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from
 
 std::set< int > neededRows;
 for( unsigned i = 0; i < tlValCol.get_n(); i++ )
 {
 int rRowInd = tlValCol.vi_wr[3 * i + 1];
 neededRows.insert( rRowInd );
 // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
 }
 moab::TupleList tgtAreaReq;
 tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
 tgtAreaReq.enableWriteAccess();
 for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
 {
 int neededRow = *sit;
 int procRow = neededRow / nBbase;
 if( procRow >= size ) procRow = size - 1;
 int nr = tgtAreaReq.get_n();
 tgtAreaReq.vi_wr[2 * nr] = procRow;
 tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
 tgtAreaReq.inc_n();
 }
 
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
 // we need to send back the tgtArea corresponding to row
 moab::TupleList tgtAreaInfo; // load it with tgtArea at row
 tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
 tgtAreaInfo.enableWriteAccess();
 for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
 {
 int from_proc = tgtAreaReq.vi_wr[2 * i];
 int row = tgtAreaReq.vi_wr[2 * i + 1];
 int locaIndexRow = row - rank * nBbase;
 double areaToSend = vecTargetFaceArea[locaIndexRow];
 // int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;
 
 tgtAreaInfo.vi_wr[2 * i] = from_proc; // send back requested info
 tgtAreaInfo.vi_wr[2 * i + 1] = row;
 tgtAreaInfo.vr_wr[i] = areaToSend; // this will be tgt area at row
 tgtAreaInfo.inc_n();
 }
 ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );
 
 std::map< int, double > areaAtRow;
 for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
 {
 // we have received from proc, value for row !
 int row = tgtAreaInfo.vi_wr[2 * i + 1];
 areaAtRow[row] = tgtAreaInfo.vr_wr[i];
 }
 
 // we have now for rows the
 // it is ordered by index, so:
 // now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
 // tgtAreaInfo will have at index i the area we need (from row)
 // there should be an easier way :(
 for( unsigned i = 0; i < tlValCol.get_n(); i++ )
 {
 int rRowInd = tlValCol.vi_wr[3 * i + 1];
 int colInd = tlValCol.vi_wr[3 * i + 2];
 double val = tlValCol.vr_wr[i];
 int localColInd = colInd - rank * nAbase; // < local nA
 // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
 auto itMap = areaAtRow.find( rRowInd ); // it should be different from end
 if( itMap != areaAtRow.end() )
 {
 double areaRow = itMap->second; // we fished a lot for this !
 dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
 }
 }
 
#endif
 // Load in data
 NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
 
 NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
 NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
 NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varRow, is_independent );
 ncMap.enable_var_par_access( varCol, is_independent );
 ncMap.enable_var_par_access( varS, is_independent );
#endif
 
 varRow->set_cur( (long)offbuf[2] );
 varRow->put( vecRow, nS );
 
 varCol->set_cur( (long)offbuf[2] );
 varCol->put( vecCol, nS );
 
 varS->set_cur( (long)offbuf[2] );
 varS->put( &( vecS[0] ), nS );
 
 // Calculate and write fractional coverage arrays
 NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varFracA, is_independent );
#endif
 varFracA->add_att( "name", "fraction of target coverage of source dof" );
 varFracA->add_att( "units", "unitless" );
 varFracA->set_cur( (long)offbuf[0] );
 varFracA->put( &( dFracA[0] ), nA );
 
 NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
 ncMap.enable_var_par_access( varFracB, is_independent );
#endif
 varFracB->add_att( "name", "fraction of source coverage of target dof" );
 varFracB->add_att( "units", "unitless" );
 varFracB->set_cur( (long)offbuf[1] );
 varFracB->put( &( dFracB[0] ), nB );
 
 // Add global attributes
 // std::map<std::string, std::string>::const_iterator iterAttributes =
 // mapAttributes.begin();
 // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
 // ncMap.add_att(
 // iterAttributes->first.c_str(),
 // iterAttributes->second.c_str());
 // }
 
 ncMap.close();
 
 return moab::MB_SUCCESS;
}

References moab::TupleList::enableWriteAccess(), moab::error(), ErrorCode, moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_CHK_SET_ERR, MB_SUCCESS, nr, moab::TempestRemapper::RLL, size, moab::Remapper::SourceMesh, moab::Remapper::TargetMesh, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.

Member Data Documentation

col_dtoc_dofmap

std::vector< int > moab::TempestOnlineMap::col_dtoc_dofmap

private

Definition at line 445 of file TempestOnlineMap.hpp.

col_gdofmap

std::vector< unsigned > moab::TempestOnlineMap::col_gdofmap

private

Definition at line 442 of file TempestOnlineMap.hpp.

Referenced by fill_col_ids(), and LinearRemapNN_MOAB().

colMap

std::map< int, int > moab::TempestOnlineMap::colMap

private

Definition at line 447 of file TempestOnlineMap.hpp.

dataGLLNodesDest

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesDest

private

Definition at line 449 of file TempestOnlineMap.hpp.

dataGLLNodesSrc

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrc

private

Definition at line 449 of file TempestOnlineMap.hpp.

dataGLLNodesSrcCov

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrcCov

private

Definition at line 449 of file TempestOnlineMap.hpp.

is_parallel

bool moab::TempestOnlineMap::is_parallel

private

Definition at line 464 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

is_root

bool moab::TempestOnlineMap::is_root

private

Definition at line 464 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_bConserved

bool moab::TempestOnlineMap::m_bConserved

private

Definition at line 456 of file TempestOnlineMap.hpp.

m_destDiscType

DiscretizationType moab::TempestOnlineMap::m_destDiscType

private

Definition at line 450 of file TempestOnlineMap.hpp.

m_dofTagDest

moab::Tag moab::TempestOnlineMap::m_dofTagDest

private

Definition at line 441 of file TempestOnlineMap.hpp.

m_dofTagSrc

moab::Tag moab::TempestOnlineMap::m_dofTagSrc

private

The original tag data and local to global DoF mapping to associate matrix values to solution

Definition at line 441 of file TempestOnlineMap.hpp.

m_eInputType

DiscretizationType moab::TempestOnlineMap::m_eInputType

private

Definition at line 455 of file TempestOnlineMap.hpp.

m_eOutputType

DiscretizationType moab::TempestOnlineMap::m_eOutputType

private

Definition at line 455 of file TempestOnlineMap.hpp.

m_iMonotonicity

int moab::TempestOnlineMap::m_iMonotonicity

private

Definition at line 457 of file TempestOnlineMap.hpp.

m_interface

moab::Interface* moab::TempestOnlineMap::m_interface

private

The reference to the moab::Core object that contains source/target and overlap sets.

Definition at line 429 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshInput

Mesh* moab::TempestOnlineMap::m_meshInput

private

Definition at line 459 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshInputCov

Mesh* moab::TempestOnlineMap::m_meshInputCov

private

Definition at line 460 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshOutput

Mesh* moab::TempestOnlineMap::m_meshOutput

private

Definition at line 461 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshOverlap

Mesh* moab::TempestOnlineMap::m_meshOverlap

private

Definition at line 462 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_nDofsPEl_Dest

int moab::TempestOnlineMap::m_nDofsPEl_Dest

private

Definition at line 454 of file TempestOnlineMap.hpp.

m_nDofsPEl_Src

int moab::TempestOnlineMap::m_nDofsPEl_Src

private

Definition at line 454 of file TempestOnlineMap.hpp.

m_nTotDofs_Dest

int moab::TempestOnlineMap::m_nTotDofs_Dest

private

Definition at line 451 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

m_nTotDofs_Src

int moab::TempestOnlineMap::m_nTotDofs_Src

private

Definition at line 451 of file TempestOnlineMap.hpp.

m_nTotDofs_SrcCov

int moab::TempestOnlineMap::m_nTotDofs_SrcCov

private

Definition at line 451 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

m_remapper

moab::TempestRemapper* moab::TempestOnlineMap::m_remapper

private

The fundamental remapping operator object.

Definition at line 416 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_srcDiscType

DiscretizationType moab::TempestOnlineMap::m_srcDiscType

private

Definition at line 450 of file TempestOnlineMap.hpp.

rank

int moab::TempestOnlineMap::rank

private

Definition at line 465 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

row_dtoc_dofmap

std::vector< int > moab::TempestOnlineMap::row_dtoc_dofmap

private

Definition at line 445 of file TempestOnlineMap.hpp.

row_gdofmap

std::vector< unsigned > moab::TempestOnlineMap::row_gdofmap

private

Definition at line 442 of file TempestOnlineMap.hpp.

Referenced by fill_row_ids(), and LinearRemapNN_MOAB().

rowMap

std::map< int, int > moab::TempestOnlineMap::rowMap

private

Definition at line 447 of file TempestOnlineMap.hpp.

size

int moab::TempestOnlineMap::size

private

Definition at line 465 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

srccol_dtoc_dofmap

std::vector< int > moab::TempestOnlineMap::srccol_dtoc_dofmap

private

Definition at line 445 of file TempestOnlineMap.hpp.

srccol_gdofmap

std::vector< unsigned > moab::TempestOnlineMap::srccol_gdofmap

private

Definition at line 442 of file TempestOnlineMap.hpp.

---

The documentation for this class was generated from the following files:

• TempestOnlineMap.hpp
• TempestLinearRemap.cpp
• TempestOnlineMap.cpp
• TempestOnlineMapIO.cpp