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 = 0 , DiscretizationType_CGLL = 1 , DiscretizationType_DGLL = 2 , DiscretizationType_PCLOUD = 3 }
enum	CAASType {   CAAS_NONE = 0 , CAAS_GLOBAL = 1 , CAAS_LOCAL = 2 , CAAS_LOCAL_ADJACENT = 3 ,   CAAS_QLT = 4 }
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 > &tgt_dof_ids, std::vector< double > &areaA, int &nA, std::vector< double > &areaB, int &nB)
	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...
void	PrintMapStatistics ()
	Print information and metadata about the remapping weights. 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 to be used for mapping. Parameters srcDofTagName The tag name associated with global DoF ids for the source mesh tgtDofTagName The tag name associated with global DoF ids for the target mesh More...
moab::ErrorCode	SetDOFmapAssociation (DiscretizationType srcType, int srcOrder, bool isSrcContinuous, DataArray3D< int > *srcdataGLLNodes, DataArray3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, int destOrder, bool isTgtContinuous, DataArray3D< int > *tgtdataGLLNodes)
std::pair< double, double >	ApplyBoundsLimiting (std::vector< double > &dataInDouble, std::vector< double > &dataOutDouble, CAASType caasType=CAAS_GLOBAL, int caasIteration=0, double mismatch=0.0)
void	ComputeAdjacencyRelations (std::vector< std::unordered_set< int > > &vecAdjFaces, int nrings, const Range &entities, bool useMOABAdjacencies=true, Mesh *trMesh=nullptr)
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 (moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose=false, CAASType caasType=CAAS_NONE, double default_projection=0.0)
	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_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)
void	SetMeshInput (Mesh *imesh)

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...
template<typename SparseMatrixType >
void	serializeSparseMatrix (const SparseMatrixType &mat, const std::string &filename)
void	setup_sizes_dimensions ()
void	CAASLimiter (std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, double &dMass)
double	QLTLimiter (int caasIteration, std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, std::vector< double > &dMassDefect)
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...

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
int	m_input_order
int	m_output_order
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 61 of file TempestOnlineMap.hpp.

Member Typedef Documentation

sample_function

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

Definition at line 414 of file TempestOnlineMap.hpp.

Member Enumeration Documentation

CAASType

enum moab::TempestOnlineMap::CAASType

Enumerator
CAAS_NONE
CAAS_GLOBAL
CAAS_LOCAL
CAAS_LOCAL_ADJACENT
CAAS_QLT

Definition at line 86 of file TempestOnlineMap.hpp.

    {
        CAAS_NONE           = 0,
        CAAS_GLOBAL         = 1,
        CAAS_LOCAL          = 2,
        CAAS_LOCAL_ADJACENT = 3,
        CAAS_QLT            = 4
    };

DiscretizationType

enum moab::TempestOnlineMap::DiscretizationType

Enumerator
DiscretizationType_FV
DiscretizationType_CGLL
DiscretizationType_DGLL
DiscretizationType_PCLOUD

Definition at line 77 of file TempestOnlineMap.hpp.

    {
        DiscretizationType_FV     = 0,
        DiscretizationType_CGLL   = 1,
        DiscretizationType_DGLL   = 2,
        DiscretizationType_PCLOUD = 3
    };

Constructor & Destructor Documentation

TempestOnlineMap()

moab::TempestOnlineMap::TempestOnlineMap

(

moab::TempestRemapper *

remapper

)

Generate the metadata associated with the offline map.

Definition at line 64 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
 
    // now let us re-update the reference to the input source mesh
    m_meshInput = m_remapper->GetMesh( moab::Remapper::SourceMesh );
    // now let us re-update the reference to the covering mesh
    m_meshInputCov = m_remapper->GetCoveringMesh();
    // now let us re-update the reference to the output target mesh
    m_meshOutput = m_remapper->GetMesh( moab::Remapper::TargetMesh );
    // now let us re-update the reference to the output target mesh
    m_meshOverlap = m_remapper->GetMesh( moab::Remapper::OverlapMesh );
 
    is_parallel = remapper->is_parallel;
    is_root     = remapper->is_root;
    rank        = remapper->rank;
    size        = remapper->size;
 
    // set default order
    m_input_order = m_output_order = 1;
 
    // 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_input_order, m_interface, m_meshInput, m_meshInputCov, m_meshOutput, m_meshOverlap, m_output_order, 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 118 of file TempestOnlineMap.cpp.

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

Member Function Documentation

ApplyBoundsLimiting()

std::pair< double, double > moab::TempestOnlineMap::ApplyBoundsLimiting	(	std::vector< double > &	dataInDouble,
		std::vector< double > &	dataOutDouble,
		CAASType	caasType = CAAS_GLOBAL,
		int	caasIteration = 0,
		double	mismatch = 0.0
	)

ApplyBoundsLimiting - Apply bounds limiting to the data field

Parameters

dataInDouble	- input data field
dataOutDouble	- output data field
caasType	- type of limiter
caasIteration	- iteration number of limiter

Returns

- pair of mass defect pre and post limiter application

ApplyWeights() [1/2]

moab::ErrorCode moab::TempestOnlineMap::ApplyWeights	(	moab::Tag	srcSolutionTag,
		moab::Tag	tgtSolutionTag,
		bool	transpose = false,
		CAASType	caasType = CAAS_NONE,
		double	default_projection = 0.0
	)

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 1415 of file TempestOnlineMap.cpp.

{
    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(), default_projection );
            sents = covSrcEnts;
        }
        else
        {
            moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
            solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
                                default_projection );
            sents = covSrcEnts;
        }
        if( m_remapper->point_cloud_target )
        {
            moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
            solTTagVals.resize( tgtEnts.size(), default_projection );
            tents = tgtEnts;
        }
        else
        {
            moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
            solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
                                    this->GetDestinationNDofsPerElement(),
                                default_projection );
            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(),
                            default_projection );
        solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
                                this->GetDestinationNDofsPerElement(),
                            default_projection );
 
        sents = covSrcEnts;
        tents = tgtEnts;
    }
 
    // The tag data is np*np*n_el_src
    MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
                    "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
    MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, transpose ),
                    "Applying remap operator onto source vector data failed" );
 
    // The tag data is np*np*n_el_dest
    MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
                    "Setting target tag data failed" );
 
    if( caasType != CAAS_NONE )
    {
        std::string tgtSolutionTagName;
        MB_CHK_SET_ERR( m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName ), "Getting tag name failed" );
 
        // Perform CAAS iterations iteratively until convergence
        constexpr int nmax_caas_iterations = 10;
        double mismatch                    = 1.0;
        int caasIteration                  = 0;
        double initialMismatch             = 0.0;
        while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
               caasIteration++ < nmax_caas_iterations )  // iterate until convergence or a maximum of 5 iterations
        {
            double dMassDiffPostGlobal;
            std::pair< double, double > mDefect =
                this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
#ifdef MOAB_HAVE_MPI
            double dMassDiffPost = mDefect.second;
            MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
#else
            dMassDiffPostGlobal = mDefect.second;
#endif
            if( caasIteration == 1 ) initialMismatch = mDefect.first;
            if( m_remapper->verbose && is_root )
            {
                printf( "Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
                        tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
            }
            mismatch = dMassDiffPostGlobal;
 
            // The tag data is np*np*n_el_dest
            MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
                            "Setting local tag data failed" );
        }
    }
 
    return moab::MB_SUCCESS;
}

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

Referenced by main().

ApplyWeights() [2/2]

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

private

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

CAASLimiter()

void moab::TempestOnlineMap::CAASLimiter ( std::vector< double > &  dataCorrectedField, std::vector< double > &  dataLowerBound, std::vector< double > &  dataUpperBound, double &  dMass  )

private

ComputeAdjacencyRelations()

void moab::TempestOnlineMap::ComputeAdjacencyRelations	(	std::vector< std::unordered_set< int > > &	vecAdjFaces,
		int	nrings,
		const Range &	entities,
		bool	useMOABAdjacencies = true,
		Mesh *	trMesh = nullptr
	)

Parameters

vecAdjFaces
nrings
entities
useMOABAdjacencies
trMesh

Returns

Vector storing adjacent Faces.

Definition at line 1366 of file TempestOnlineMap.cpp.

{
    assert( nrings > 0 );
    assert( useMOABAdjacencies || trMesh != nullptr );
 
    const size_t nrows = vecAdjFaces.size();
    moab::MeshTopoUtil mtu( m_interface );
    for( size_t index = 0; index < nrows; index++ )
    {
        vecAdjFaces[index].insert( index );  // add self target face first
        {
            // Compute the adjacent faces to the target face
            if( useMOABAdjacencies )
            {
                moab::Range ents;
                // ents.insert( entities.index( entities[index] ) );
                ents.insert( entities[index] );
                moab::Range adjEnts;
                moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, nrings );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
                for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
                {
                    // int adjIndex = m_interface->id_from_handle(*it)-1;
                    int adjIndex = entities.index( *it );
                    // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
                    if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
                }
            }
            else
            {
                ///  Vector storing adjacent Faces.
                typedef std::pair< int, int > FaceDistancePair;
                typedef std::vector< FaceDistancePair > AdjacentFaceVector;
                AdjacentFaceVector adjFaces;
                Face& face = trMesh->faces[index];
                GetAdjacentFaceVectorByEdge( *trMesh, index, nrings * face.edges.size(), adjFaces );
 
                // Add the adjacent faces to the target face list
                for( auto adjFace : adjFaces )
                    if( adjFace.first >= 0 )
                        vecAdjFaces[index].insert( adjFace.first );  // map target face to source face
            }
        }
    }
}

References moab::Range::begin(), moab::Range::end(), entities, ErrorCode, moab::MeshTopoUtil::get_bridge_adjacencies(), moab::Range::insert(), and MB_CHK_SET_ERR_CONT.

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 1884 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;
    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 );
 
    const auto& ovents = m_remapper->m_overlap_entities;
 
    std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 );  //  L1Err, L2Err, LinfErr
    double sumarea = 0.0;
    for( size_t i = 0; i < ovents.size(); ++i )
    {
        const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
        if( srcidx < 0 ) continue;  // Skip non-overlapping entities
        const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
        if( tgtidx < 0 ) continue;  // skip ghost target faces
        const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
        const double error  = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
        errnorms[0] += ovarea * error;
        errnorms[1] += ovarea * error * error;
        errnorms[3] = ( error > errnorms[3] ? error : errnorms[3] );
        sumarea += ovarea;
    }
    errnorms[2] = sumarea;
#ifdef MOAB_HAVE_MPI
    if( m_pcomm )
    {
        MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
        MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
    }
#else
    for( int i = 0; i < 4; ++i )
        globerrnorms[i] = errnorms[i];
#endif
 
    globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
    globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
 
    metrics.clear();
    metrics["L1Error"]   = globerrnorms[0];
    metrics["L2Error"]   = globerrnorms[1];
    metrics["LinfError"] = globerrnorms[3];
 
    if( verbose && is_root )
    {
        std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
        std::cout << "\t Total Intersection area = " << globerrnorms[2] << 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[3] << 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 1522 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 != nullptr )
    {
        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 439 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.

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 431 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 );
 
        // Set the source and target mesh objects
        MB_CHK_ERR( SetDOFmapTags( srcDofTagName, tgtDofTagName ) );
 
        ///   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 local_areas[3] = { 0.0, 0.0, 0.0 }, global_areas[3] = { 0.0, 0.0, 0.0 };
        if( !m_bPointCloudSource )
        {
            // Calculate Input Mesh Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
            local_areas[0] = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
            // Set source element areas as tag on the source mesh
            MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea ) );
 
            // Update coverage source mesh areas as well.
            m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
        }
 
        if( !m_bPointCloudTarget )
        {
            // Calculate Output Mesh Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
            local_areas[1] = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
            // Set target element areas as tag on the target mesh
            MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea ) );
        }
 
        if( !m_bPointCloud )
        {
            // Verify that overlap mesh is in the correct order (sanity check)
            assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->vecTargetFaceIx.size() );
 
            // Calculate Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
            local_areas[2] =
                m_meshOverlap->CalculateFaceAreas( mapOptions.fSourceConcave || mapOptions.fTargetConcave );
 
            // store it as global output for now - used later in reduction
            std::copy( local_areas, local_areas + 3, global_areas );
#ifdef MOAB_HAVE_MPI
            // reduce the local source, target and overlap mesh areas to global areas
            if( m_pcomm && is_parallel )
                MPI_Reduce( local_areas, global_areas, 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
            if( is_root )
            {
                dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", global_areas[0] );
                dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", global_areas[1] );
                dbgprint.printf( 0, "Overlap Mesh Recovered Area: %1.15e\n", global_areas[2] );
            }
 
            // Correct areas to match the areas calculated in the overlap mesh
            constexpr bool fCorrectAreas = true;
            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 );
 
            this->m_pdataGLLNodesIn  = nullptr;
            this->m_pdataGLLNodesOut = nullptr;
 
            // Finite volume input / Finite element output
            rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
                                               mapOptions.nPout, false, nullptr );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 ) dbgprint.printf( 0, "Calculating map (invdist)\n" );
                if( m_meshInputCov->faces.size() )
                    LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
            }
            else if( strMapAlgorithm == "delaunay" )
            {
                if( is_root ) dbgprint.printf( 0, "Calculating map (delaunay)\n" );
                if( m_meshInputCov->faces.size() )
                    LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
            }
            else if( strMapAlgorithm == "fvintbilin" )
            {
                if( is_root ) dbgprint.printf( 0, "Calculating map (intbilin)\n" );
                if( m_meshInputCov->faces.size() )
                    LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
            }
            else if( strMapAlgorithm == "fvintbilingb" )
            {
                if( is_root ) dbgprint.printf( 0, "Calculating map (intbilingb)\n" );
                if( m_meshInputCov->faces.size() )
                    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 ) dbgprint.printf( 0, "Calculating map (bilin)\n" );
                if( m_meshInputCov->faces.size() )
                    LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
            }
            else
            {
                if( is_root ) dbgprint.printf( 0, "Calculating conservative FV-FV map\n" );
                if( m_meshInputCov->faces.size() )
                {
#ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
                    LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
                                       ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
#else
                    LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
#endif
                }
            }
        }
        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 );
 
            this->m_pdataGLLNodesIn  = nullptr;
            this->m_pdataGLLNodesOut = &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, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
                                               mapOptions.nPout, ( 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, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
                ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrcCov ),
                ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrc ),
                eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
                ( m_bPointCloudTarget ? nullptr : &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" );
            // generate metadata for the input meshes (both source and covering source)
            GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
                              dataGLLJacobianSrc );
            GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
                              dataGLLJacobian );
 
            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, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
                                               &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
                                               false, nullptr );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" );
            }
 
            this->m_pdataGLLNodesIn  = &dataGLLNodesSrcCov;
            this->m_pdataGLLNodesOut = nullptr;
 
#ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
            LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
                            nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
                            *this );
#else
            LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
                                         mapOptions.fNoConservation );
#endif
        }
        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" );
            // generate metadata for the input meshes (both source and covering source)
            GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
                              dataGLLJacobianSrc );
            // now coverage
            GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
                              dataGLLJacobianIn );
            // Output metadata
            if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
            GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
                              dataGLLJacobianOut );
 
            // 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, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
                                               &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
                                               ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
 
            this->m_pdataGLLNodesIn  = &dataGLLNodesSrcCov;
            this->m_pdataGLLNodesOut = &dataGLLNodesDest;
 
            // Generate remap
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
 
#ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
            LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
                                            dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
                                            this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
                                            fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *this );
#else
            LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
                                       this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
                                       fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
#endif
        }
        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, and moab::Remapper::OverlapMesh.

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 567 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 585 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 572 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 561 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 556 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 578 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 1185 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;
    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 );
        }
 
        printf( "Total nnz: %d, global source elements = %d\n", nTotVals, 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 1139 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 1328 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 121 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_meshInputCov" );
    }
 
    // 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;
 
    // Announcements
    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_meshInputCov
    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 1583 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_meshInputCov
    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, and ForceIntArrayConsistencyConservation().

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 2134 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.

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 58 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_FAILURE;
    }
}

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 1162 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++ )
        {
            // if( m_meshOverlap->vecTargetFaceIx[ixOverlapTemp] < 0 ) continue;  // skip ghost target faces
            // 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().

PrintMapStatistics()

void moab::TempestOnlineMap::PrintMapStatistics

(

)

Print information and metadata about the remapping weights.

Definition at line 284 of file TempestLinearRemap.cpp.

{
    int nrows = m_weightMatrix.rows();      // Number of rows
    int ncols = m_weightMatrix.cols();      // Number of columns
    int NNZ   = m_weightMatrix.nonZeros();  // Number of non zero values
#ifdef MOAB_HAVE_MPI
    // find out min/max for NNZ, ncols, nrows
    // should work on std c++ 11
    int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
    int rarr3[6];
    MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
 
    int total[3];
    MPI_Reduce( arr3, total, 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
    if( !rank )
        std::cout << "-> Rows (min/max/sum): (" << rarr3[1] << " / " << -rarr3[4] << " / " << total[1] << "), "
                  << " Cols (min/max/sum): (" << rarr3[2] << " / " << -rarr3[5] << " / " << total[2] << "), "
                  << " NNZ (min/max/sum): (" << rarr3[0] << " / " << -rarr3[3] << " / " << total[0] << ")\n";
#else
    std::cout << "-> Rows: " << nrows << ", Cols: " << ncols << ", NNZ: " << NNZ << "\n";
#endif
}

Referenced by main().

QLTLimiter()

double moab::TempestOnlineMap::QLTLimiter ( int  caasIteration, std::vector< double > &  dataCorrectedField, std::vector< double > &  dataLowerBound, std::vector< double > &  dataUpperBound, std::vector< double > &  dMassDefect  )

private

ReadParallelMap()

moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap	(	const char *	strSource,
		const std::vector< int > &	tgt_dof_ids,
		std::vector< double > &	areaA,
		int &	nA,
		std::vector< double > &	areaB,
		int &	nB
	)

Generate the metadata associated with the offline map.

Read the OfflineMap from a NetCDF file.

Definition at line 1205 of file TempestOnlineMapIO.cpp.

{
    NcError error( NcError::silent_nonfatal );
 
    NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
    NcVar *varAreaA = NULL, *varAreaB = NULL;
    int nS = 0;
#ifdef MOAB_HAVE_PNETCDF
    // some variables will be used just in the case netcdfpar reader fails
    int ncfile = -1;
    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 == nullptr )                                                                      \
        {                                                                                         \
            _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" );
 
        varAreaA = ncMap.get_var( "area_a" );
        CHECK_EXCEPTION( varAreaA, "variable", "area_a" );
 
        varAreaB = ncMap.get_var( "area_b" );
        CHECK_EXCEPTION( varAreaB, "variable", "area_b" );
 
#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 );
        ncMap.enable_var_par_access( varAreaA, is_independent );
        ncMap.enable_var_par_access( varAreaB, 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 );
        ERR_PARNC(
            ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile ) );  // bail out completely
        ERR_PARNC( ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited ) );
        // find dimension ids for n_S
        int ins;
        ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_s", &ins ) );
        MPI_Offset leng;
        ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
        nS = (int)leng;
        ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_a", &ins ) );
        ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
        nA = (int)leng;
        ERR_PARNC( ncmpi_inq_dimid( ncfile, "n_b", &ins ) );
        ERR_PARNC( ncmpi_inq_dimlen( ncfile, ins, &leng ) );
        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;
    }
 
    int localSizeA   = nA / size;
    long offsetReadA = rank * localSizeA;
    // leftovers on last rank
    if( rank == size - 1 )
    {
        localSizeA += nA % size;
    }
 
    int localSizeB   = nB / size;
    long offsetReadB = rank * localSizeB;
    // leftovers on last rank
    if( rank == size - 1 )
    {
        localSizeB += nB % size;
    }
 
    std::vector< int > vecRow, vecCol;
    std::vector< double > vecS;
    vecRow.resize( localSize );
    vecCol.resize( localSize );
    vecS.resize( localSize );
    vecAreaA.resize( localSizeA );
    vecAreaB.resize( localSizeB );
 
    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 );
 
        varAreaA->set_cur( (long)( offsetReadA ) );
        varAreaA->get( &( vecAreaA[0] ), localSizeA );
 
        varAreaB->set_cur( (long)( offsetReadB ) );
        varAreaB->get( &( vecAreaB[0] ), localSizeB );
 
        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;
        // get the sparse matrix values, row and column indices
        ERR_PARNC( ncmpi_inq_varid( ncfile, "S", &varid ) );
        ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] ) );
        ERR_PARNC( ncmpi_inq_varid( ncfile, "row", &varid ) );
        ERR_PARNC( ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] ) );
        ERR_PARNC( ncmpi_inq_varid( ncfile, "col", &varid ) );
        ERR_PARNC( ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] ) );
 
        ERR_PARNC( ncmpi_inq_varid( ncfile, "area_a", &varid ) );
        MPI_Offset startA = (MPI_Offset)offsetReadA;
        MPI_Offset countA = (MPI_Offset)localSizeA;
        ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &startA, &countA, &vecAreaA[0] ) );
 
        ERR_PARNC( ncmpi_inq_varid( ncfile, "area_b", &varid ) );
        MPI_Offset startB = (MPI_Offset)offsetReadB;
        MPI_Offset countB = (MPI_Offset)localSizeB;
        ERR_PARNC( ncmpi_get_vara_double_all( ncfile, varid, &startB, &countB, &vecAreaB[0] ) );
 
        ERR_PARNC( ncmpi_close( ncfile ) );
#endif
    }
 
    // NOTE: we can check if the following workflow would help here.
    // Initial experiments did not show same index map.
    // Might be useful to understand why that is the case.
    // First, ensure that all local and global indices are computed
    // this->m_remapper->ComputeGlobalLocalMaps();
    // Next, reuse the global to local map
    // const auto& rowMap = m_remapper->gid_to_lid_tgt;
    // const auto& colMap = m_remapper->gid_to_lid_covsrc;
#ifdef MOAB_HAVE_EIGEN3
 
    typedef Eigen::Triplet< double > Triplet;
    std::vector< Triplet > tripletList;
 
#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
 
        int nPerPart        = nDofs / size;
        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; sparse matrix is 0 based
            int colval  = vecCol[i] - 1;
            int to_proc = -1;
 
            to_proc = rowval / nPerPart;
            if( to_proc == size ) to_proc = size - 1;
 
            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 rendezvous 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 ?
                to_proc     = dof_val / nPerPart;
                if( to_proc == size ) to_proc = size - 1;
 
                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() );
 
            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  tl_re.vi_rd[2 * 0 + 1];
 
                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 + 1];  // 1 for row, 2 for column // first value, it should be
                if (startDofIndex.find(valDof)==startDofIndex.end())
                    continue;
                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;
        }
 
        // set of row and col used on this task
        std::set< int > rowSet;
        std::set< int > colSet;
        // populate the sparsematrix, using rowMap and colMap
        int n = tl->get_n();
        for( int i = 0; i < n; i++ )
        {
            const int vecRowValue = tl->vi_wr[3 * i + 1];
            const int vecColValue = tl->vi_wr[3 * i + 2];
            rowSet.insert( vecRowValue );
            colSet.insert( vecColValue );
        }
        int index = 0;
        row_gdofmap.resize( rowSet.size() );
        for( auto setIt : rowSet )
        {
            row_gdofmap[index] = setIt;
            rowMap[setIt]      = index++;
        }
        m_nTotDofs_Dest = index;
        index           = 0;
        col_gdofmap.resize( colSet.size() );
        for( auto setIt : colSet )
        {
            col_gdofmap[index] = setIt;
            colMap[setIt]      = index++;
        }
        m_nTotDofs_SrcCov = index;
 
        tripletList.reserve( n );
        for( int i = 0; i < n; i++ )
        {
            const int vecRowValue = tl->vi_wr[3 * i + 1];
            const int vecColValue = tl->vi_wr[3 * i + 2];
            double value          = tl->vr_wr[i];
            tripletList.emplace_back( rowMap[vecRowValue], colMap[vecColValue], value );
        }
        tl->reset();
    }
    else
#endif
    {
        // set of row and col used on this task
        std::set< int > rowSet;
        std::set< int > colSet;
        // populate the sparsematrix, using rowMap and colMap
        for( int i = 0; i < nS; i++ )
        {
            const int vecRowValue = vecRow[i] - 1;
            const int vecColValue = vecCol[i] - 1;
            rowSet.insert( vecRowValue );
            colSet.insert( vecColValue );
        }
 
        int index = 0;
        row_gdofmap.resize( rowSet.size() );
        for( auto setIt : rowSet )
        {
            row_gdofmap[index] = setIt;
            rowMap[setIt]      = index++;
        }
        m_nTotDofs_Dest = index;
        index           = 0;
        col_gdofmap.resize( colSet.size() );
        for( auto setIt : colSet )
        {
            col_gdofmap[index] = setIt;
            colMap[setIt]      = index++;
        }
        m_nTotDofs_SrcCov = index;
 
        tripletList.reserve( nS );
        for( int i = 0; i < nS; i++ )
        {
            const int vecRowValue = vecRow[i] - 1;  // the rows, cols are 1 based in the file
            const int vecColValue = vecCol[i] - 1;  // sparse matrix will be 0 based
            double value          = vecS[i];
            tripletList.emplace_back( rowMap[vecRowValue], colMap[vecColValue], value );
        }
    }
 
    m_weightMatrix.resize( m_nTotDofs_Dest, m_nTotDofs_SrcCov );
    m_rowVector.resize( m_nTotDofs_Dest );
    m_colVector.resize( m_nTotDofs_SrcCov );
    m_nTotDofs_Src = m_nTotDofs_SrcCov;  // do we need both?
    m_weightMatrix.setFromTriplets( tripletList.begin(), tripletList.end() );
    // Reset the source and target data first
    m_rowVector.setZero();
    m_colVector.setZero();
#ifdef VERBOSE
    serializeSparseMatrix( m_weightMatrix, "map_operator_" + std::to_string( rank ) + ".txt" );
#endif
// #ifdef MOAB_HAVE_EIGEN3
#endif
    // TODO: make this flexible and read the order from map with help of metadata
    m_nDofsPEl_Src  = 1;  // always assume FV-FV maps are read from file
    m_nDofsPEl_Dest = 1;  // always assume FV-FV maps are read from file
 
    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.

serializeSparseMatrix()

template<typename SparseMatrixType >

void moab::TempestOnlineMap::serializeSparseMatrix ( const SparseMatrixType &  mat, const std::string &  filename  )

private

set_col_dc_dofs()

moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs

(

std::vector< int > &

values_entities

)

Definition at line 397 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
    //
    // resize and initialize to -1 to signal that this value should not be used, if not set below
    col_dtoc_dofmap.resize( values_entities.size(), -1 );
    for( size_t j = 0; j < values_entities.size(); j++ )
    {
        // values are 1 based, but rowMap, colMap are not
        const auto it = colMap.find( values_entities[j] - 1 );
        if( it != colMap.end() ) col_dtoc_dofmap[j] = it->second;
    }
    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 416 of file TempestOnlineMap.cpp.

{
    //  we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
    // resize and initialize to -1 to signal that this value should not be used, if not set below
    row_dtoc_dofmap.resize( values_entities.size(), -1 );
    for( size_t j = 0; j < values_entities.size(); j++ )
    {
        // values are 1 based, but rowMap, colMap are not
        const auto it = rowMap.find( values_entities[j] - 1 );
        if( it != rowMap.end() ) row_dtoc_dofmap[j] = it->second;
    }
    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 596 of file TempestOnlineMap.hpp.

{
    m_nDofsPEl_Dest = nt;
}

SetDOFmapAssociation()

moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation	(	DiscretizationType	srcType,
		int	srcOrder,
		bool	isSrcContinuous,
		DataArray3D< int > *	srcdataGLLNodes,
		DataArray3D< int > *	srcdataGLLNodesSrc,
		DiscretizationType	destType,
		int	destOrder,
		bool	isTgtContinuous,
		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.

Parameters

srcType	The discretization type of the source mesh
srcOrder	The order of the discretization on the source mesh
isSrcContinuous	The continuity of the discretization on the source mesh
srcdataGLLNodes	The GLL nodes on the source mesh
srcdataGLLNodesSrc	The GLL nodes on the source mesh
destType	The discretization type of the destination mesh
destOrder	The order of the discretization on the destination mesh
isTgtContinuous	The continuity of the discretization on the destination mesh
tgtdataGLLNodes	The GLL nodes on the destination mesh

Definition at line 163 of file TempestOnlineMap.cpp.

{
    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;
    m_input_order  = srcOrder;
    m_output_order = destOrder;
 
    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.
    // 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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] ) );
        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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] ) );
    }
 
    m_nTotDofs_SrcCov = 0;
    if( srcdataGLLNodes == nullptr )
    {
        /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < col_gdofmap.size(); ++i )
        {
            auto gdof = src_soln_gdofs[i];
            assert( gdof > 0 );
            col_gdofmap[i]     = gdof - 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: " << 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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] ) );
    }
    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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] ) );
    }
 
    // Now compute the mapping and store it for the original source mesh
    m_nTotDofs_Src = 0;
    if( srcdataGLLNodesSrc == nullptr )
    {
        /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
        {
            auto gdof = locsrc_soln_gdofs[i];
            assert( gdof > 0 );
            srccol_gdofmap[i]     = gdof - 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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] ) );
        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 );
        MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] ) );
    }
 
    // 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 == nullptr )
    {
        /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < row_gdofmap.size(); ++i )
        {
            auto gdof = tgt_soln_gdofs[i];
            assert( gdof > 0 );
            row_gdofmap[i]     = gdof - 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: " << 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 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 to be used for mapping.

Parameters

srcDofTagName	The tag name associated with global DoF ids for the source mesh
tgtDofTagName	The tag name associated with global DoF ids for the target mesh

Definition at line 131 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 )
    {
        MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for source mesh." );
    }
    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 )
    {
        MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for target mesh." );
    }
    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.

SetMeshInput()

void moab::TempestOnlineMap::SetMeshInput ( Mesh *  imesh )

inline

Definition at line 453 of file TempestOnlineMap.hpp.

    {
        m_meshInput = imesh;
    };

References m_meshInput.

SetSourceNDofsPerElement()

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

inline

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

Definition at line 591 of file TempestOnlineMap.hpp.

{
    m_nDofsPEl_Src = ns;
}

setup_sizes_dimensions()

void moab::TempestOnlineMap::setup_sizes_dimensions ( )

private

Definition at line 93 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 829 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, 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 );
        }
    }
 
    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, and size.

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 193 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 218 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();
 
#ifdef VERBOSE
    serializeSparseMatrix( m_weightMatrix, "map_operator_" + std::to_string( rank ) + ".txt" );
#endif
    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 530 of file TempestOnlineMap.hpp.

col_gdofmap

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

private

Definition at line 527 of file TempestOnlineMap.hpp.

Referenced by fill_col_ids(), and LinearRemapNN_MOAB().

colMap

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

private

Definition at line 532 of file TempestOnlineMap.hpp.

dataGLLNodesDest

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesDest

private

Definition at line 535 of file TempestOnlineMap.hpp.

dataGLLNodesSrc

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrc

private

Definition at line 535 of file TempestOnlineMap.hpp.

dataGLLNodesSrcCov

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrcCov

private

Definition at line 535 of file TempestOnlineMap.hpp.

is_parallel

bool moab::TempestOnlineMap::is_parallel

private

Definition at line 550 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

is_root

bool moab::TempestOnlineMap::is_root

private

Definition at line 550 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_bConserved

bool moab::TempestOnlineMap::m_bConserved

private

Definition at line 542 of file TempestOnlineMap.hpp.

m_destDiscType

DiscretizationType moab::TempestOnlineMap::m_destDiscType

private

Definition at line 536 of file TempestOnlineMap.hpp.

m_dofTagDest

moab::Tag moab::TempestOnlineMap::m_dofTagDest

private

Definition at line 526 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 526 of file TempestOnlineMap.hpp.

m_eInputType

DiscretizationType moab::TempestOnlineMap::m_eInputType

private

Definition at line 541 of file TempestOnlineMap.hpp.

m_eOutputType

DiscretizationType moab::TempestOnlineMap::m_eOutputType

private

Definition at line 541 of file TempestOnlineMap.hpp.

m_iMonotonicity

int moab::TempestOnlineMap::m_iMonotonicity

private

Definition at line 543 of file TempestOnlineMap.hpp.

m_input_order

int moab::TempestOnlineMap::m_input_order

private

Definition at line 533 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

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 514 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshInput

Mesh* moab::TempestOnlineMap::m_meshInput

private

Definition at line 545 of file TempestOnlineMap.hpp.

Referenced by SetMeshInput(), and TempestOnlineMap().

m_meshInputCov

Mesh* moab::TempestOnlineMap::m_meshInputCov

private

Definition at line 546 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshOutput

Mesh* moab::TempestOnlineMap::m_meshOutput

private

Definition at line 547 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_meshOverlap

Mesh* moab::TempestOnlineMap::m_meshOverlap

private

Definition at line 548 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_nDofsPEl_Dest

int moab::TempestOnlineMap::m_nDofsPEl_Dest

private

Definition at line 540 of file TempestOnlineMap.hpp.

m_nDofsPEl_Src

int moab::TempestOnlineMap::m_nDofsPEl_Src

private

Definition at line 540 of file TempestOnlineMap.hpp.

m_nTotDofs_Dest

int moab::TempestOnlineMap::m_nTotDofs_Dest

private

Definition at line 537 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

m_nTotDofs_Src

int moab::TempestOnlineMap::m_nTotDofs_Src

private

Definition at line 537 of file TempestOnlineMap.hpp.

m_nTotDofs_SrcCov

int moab::TempestOnlineMap::m_nTotDofs_SrcCov

private

Definition at line 537 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

m_output_order

int moab::TempestOnlineMap::m_output_order

private

Definition at line 533 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_remapper

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

private

The fundamental remapping operator object.

Definition at line 500 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

m_srcDiscType

DiscretizationType moab::TempestOnlineMap::m_srcDiscType

private

Definition at line 536 of file TempestOnlineMap.hpp.

rank

int moab::TempestOnlineMap::rank

private

Definition at line 551 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

row_dtoc_dofmap

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

private

Definition at line 530 of file TempestOnlineMap.hpp.

row_gdofmap

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

private

Definition at line 527 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

rowMap

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

private

Definition at line 532 of file TempestOnlineMap.hpp.

size

int moab::TempestOnlineMap::size

private

Definition at line 551 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

srccol_dtoc_dofmap

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

private

Definition at line 530 of file TempestOnlineMap.hpp.

srccol_gdofmap

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

private

Definition at line 527 of file TempestOnlineMap.hpp.

---

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

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