MOAB: Mesh Oriented datABase  (version 5.3.1)
moab::TempestOnlineMap Class Reference

An offline map between two Meshes. More...

#include <TempestOnlineMap.hpp>

+ Collaboration diagram for moab::TempestOnlineMap:

Public Types

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

Public Member Functions

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

Private Member Functions

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

Private Attributes

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

Detailed Description

An offline map between two Meshes.

Definition at line 48 of file TempestOnlineMap.hpp.


Member Typedef Documentation

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

Definition at line 329 of file TempestOnlineMap.hpp.


Member Enumeration Documentation

Enumerator:
DiscretizationType_FV 
DiscretizationType_CGLL 
DiscretizationType_DGLL 
DiscretizationType_PCLOUD 

Definition at line 64 of file TempestOnlineMap.hpp.


Constructor & Destructor Documentation

Generate the metadata associated with the offline map.

Definition at line 78 of file TempestOnlineMap.cpp.

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

                                                                      : OfflineMap(), m_remapper( remapper )
{
    // Get the references for the MOAB core objects
    m_interface = m_remapper->get_interface();
#ifdef MOAB_HAVE_MPI
    m_pcomm = m_remapper->get_parallel_communicator();
#endif

    // Update the references to the meshes
    m_meshInput    = remapper->GetMesh( moab::Remapper::SourceMesh );
    m_meshInputCov = remapper->GetCoveringMesh();
    m_meshOutput   = remapper->GetMesh( moab::Remapper::TargetMesh );
    m_meshOverlap  = remapper->GetMesh( moab::Remapper::OverlapMesh );

    is_parallel = false;
    is_root     = true;
    rank        = 0;
    root_proc   = rank;
    size        = 1;
#ifdef MOAB_HAVE_MPI
    int flagInit;
    MPI_Initialized( &flagInit );
    if( flagInit )
    {
        is_parallel = true;
        assert( m_pcomm != NULL );
        rank    = m_pcomm->rank();
        size    = m_pcomm->size();
        is_root = ( rank == 0 );
    }
#endif

    // Compute and store the total number of source and target DoFs corresponding
    // to number of rows and columns in the mapping.

    // Initialize dimension information from file
    this->setup_sizes_dimensions();

    // Build a matrix of source and target discretization so that we know how to assign
    // the global DoFs in parallel for the mapping weights
    // For example, FV->FV: rows X cols = faces_source X faces_target
}

Define a virtual destructor.

Definition at line 168 of file TempestOnlineMap.cpp.

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

Member Function Documentation

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

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

Referenced by main().

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

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

Definition at line 1567 of file TempestOnlineMap.cpp.

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

{
    moab::ErrorCode rval;

    std::vector< double > solSTagVals;
    std::vector< double > solTTagVals;

    moab::Range sents, tents;
    if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
    {
        if( m_remapper->point_cloud_source )
        {
            moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
            solSTagVals.resize( covSrcEnts.size(), -1.0 );
            sents = covSrcEnts;
        }
        else
        {
            moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
            solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
                                -1.0 );
            sents = covSrcEnts;
        }
        if( m_remapper->point_cloud_target )
        {
            moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
            solTTagVals.resize( tgtEnts.size(), -1.0 );
            tents = tgtEnts;
        }
        else
        {
            moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
            solTTagVals.resize(
                tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
            tents = tgtEnts;
        }
    }
    else
    {
        moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
        moab::Range& tgtEnts    = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
        solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
                            -1.0 );
        solTTagVals.resize(
            tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );

        sents = covSrcEnts;
        tents = tgtEnts;
    }

    // The tag data is np*np*n_el_src
    rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );

    // Compute the application of weights on the suorce solution data and store it in the
    // destination solution vector data Optionally, can also perform the transpose application of
    // the weight matrix. Set the 3rd argument to true if this is needed
    rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );

    // The tag data is np*np*n_el_dest
    rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

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

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

Definition at line 1995 of file TempestOnlineMap.cpp.

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

Referenced by main().

{
    moab::ErrorCode rval;
    const bool outputEnabled = ( is_root );
    int discOrder;
    DiscretizationType discMethod;
    moab::EntityHandle meshset;
    moab::Range entities;
    Mesh* trmesh;
    switch( ctx )
    {
        case Remapper::SourceMesh:
            meshset    = m_remapper->m_covering_source_set;
            trmesh     = m_remapper->m_covering_source;
            entities   = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
                                                        : m_remapper->m_covering_source_entities );
            discOrder  = m_nDofsPEl_Src;
            discMethod = m_eInputType;
            break;

        case Remapper::TargetMesh:
            meshset = m_remapper->m_target_set;
            trmesh  = m_remapper->m_target;
            entities =
                ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
            discOrder  = m_nDofsPEl_Dest;
            discMethod = m_eOutputType;
            break;

        default:
            if( outputEnabled )
                std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
            return moab::MB_FAILURE;
    }

    // Let us create teh solution tag with appropriate information for name, discretization order
    // (DoF space)
    std::string exactTagName, projTagName;
    const int ntotsize = entities.size() * discOrder * discOrder;
    int ntotsize_glob  = 0;
    std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
    rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
    rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
    rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
    rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );

    std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 );  //  L1Err, L2Err, LinfErr
    for( int i = 0; i < ntotsize; ++i )
    {
        const double error = fabs( exactSolution[i] - projSolution[i] );
        errnorms[0] += error;
        errnorms[1] += error * error;
        errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] );
    }
#ifdef MOAB_HAVE_MPI
    if( m_pcomm )
    {
        MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
        MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
        MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
    }
#else
    ntotsize_glob = ntotsize;
    globerrnorms = errnorms;
#endif
    globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob );
    globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob );

    metrics.clear();
    metrics["L1Error"]   = globerrnorms[0];
    metrics["L2Error"]   = globerrnorms[1];
    metrics["LinfError"] = globerrnorms[2];

    if( verbose && is_root )
    {
        std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
        std::cout << "\t L_1 error   = " << globerrnorms[0] << std::endl;
        std::cout << "\t L_2 error   = " << globerrnorms[1] << std::endl;
        std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl;
    }

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

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

Definition at line 1633 of file TempestOnlineMap.cpp.

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

Referenced by main().

{
    moab::ErrorCode rval;
    const bool outputEnabled = ( is_root );
    int discOrder;
    DiscretizationType discMethod;
    moab::EntityHandle meshset;
    moab::Range entities;
    Mesh* trmesh;
    switch( ctx )
    {
        case Remapper::SourceMesh:
            meshset    = m_remapper->m_covering_source_set;
            trmesh     = m_remapper->m_covering_source;
            entities   = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
                                                        : m_remapper->m_covering_source_entities );
            discOrder  = m_nDofsPEl_Src;
            discMethod = m_eInputType;
            break;

        case Remapper::TargetMesh:
            meshset = m_remapper->m_target_set;
            trmesh  = m_remapper->m_target;
            entities =
                ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
            discOrder  = m_nDofsPEl_Dest;
            discMethod = m_eOutputType;
            break;

        default:
            if( outputEnabled )
                std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
            return moab::MB_FAILURE;
    }

    // Let us create teh solution tag with appropriate information for name, discretization order
    // (DoF space)
    rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
                                        MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
    if( clonedSolnTag != NULL )
    {
        if( cloneSolnName.size() == 0 )
        {
            cloneSolnName = solnName + std::string( "Cloned" );
        }
        rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
                                            *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
    }

    // Triangular quadrature rule
    const int TriQuadratureOrder = 10;

    if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;

    TriangularQuadratureRule triquadrule( TriQuadratureOrder );

    const int TriQuadraturePoints = triquadrule.GetPoints();

    const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
    const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();

    // Output data
    DataArray1D< double > dVar;
    DataArray1D< double > dVarMB;  // re-arranged local MOAB vector

    // Nodal geometric area
    DataArray1D< double > dNodeArea;

    // Calculate element areas
    // trmesh->CalculateFaceAreas(fContainsConcaveFaces);

    if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
    {
        /* Get the spectral points and sample the functionals accordingly */
        const bool fGLL          = true;
        const bool fGLLIntegrate = false;

        // Generate grid metadata
        DataArray3D< int > dataGLLNodes;
        DataArray3D< double > dataGLLJacobian;

        GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );

        // Number of elements
        int nElements = trmesh->faces.size();

        // Verify all elements are quadrilaterals
        for( int k = 0; k < nElements; k++ )
        {
            const Face& face = trmesh->faces[k];

            if( face.edges.size() != 4 )
            {
                _EXCEPTIONT( "Non-quadrilateral face detected; "
                             "incompatible with --gll" );
            }
        }

        // Number of unique nodes
        int iMaxNode = 0;
        for( int i = 0; i < discOrder; i++ )
        {
            for( int j = 0; j < discOrder; j++ )
            {
                for( int k = 0; k < nElements; k++ )
                {
                    if( dataGLLNodes[i][j][k] > iMaxNode )
                    {
                        iMaxNode = dataGLLNodes[i][j][k];
                    }
                }
            }
        }

        // Get Gauss-Lobatto quadrature nodes
        DataArray1D< double > dG;
        DataArray1D< double > dW;

        GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );

        // Get Gauss quadrature nodes
        const int nGaussP = 10;

        DataArray1D< double > dGaussG;
        DataArray1D< double > dGaussW;

        GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );

        // Allocate data
        dVar.Allocate( iMaxNode );
        dVarMB.Allocate( discOrder * discOrder * nElements );
        dNodeArea.Allocate( iMaxNode );

        // Sample data
        for( int k = 0; k < nElements; k++ )
        {
            const Face& face = trmesh->faces[k];

            // Sample data at GLL nodes
            if( fGLL )
            {
                for( int i = 0; i < discOrder; i++ )
                {
                    for( int j = 0; j < discOrder; j++ )
                    {

                        // Apply local map
                        Node node;
                        Node dDx1G;
                        Node dDx2G;

                        ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );

                        // Sample data at this point
                        double dNodeLon = atan2( node.y, node.x );
                        if( dNodeLon < 0.0 )
                        {
                            dNodeLon += 2.0 * M_PI;
                        }
                        double dNodeLat = asin( node.z );

                        double dSample = ( *testFunction )( dNodeLon, dNodeLat );

                        dVar[dataGLLNodes[j][i][k] - 1] = dSample;
                    }
                }
                // High-order Gaussian integration over basis function
            }
            else
            {
                DataArray2D< double > dCoeff( discOrder, discOrder );

                for( int p = 0; p < nGaussP; p++ )
                {
                    for( int q = 0; q < nGaussP; q++ )
                    {

                        // Apply local map
                        Node node;
                        Node dDx1G;
                        Node dDx2G;

                        ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );

                        // Cross product gives local Jacobian
                        Node nodeCross = CrossProduct( dDx1G, dDx2G );

                        double dJacobian =
                            sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );

                        // Find components of quadrature point in basis
                        // of the first Face
                        SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );

                        // Sample data at this point
                        double dNodeLon = atan2( node.y, node.x );
                        if( dNodeLon < 0.0 )
                        {
                            dNodeLon += 2.0 * M_PI;
                        }
                        double dNodeLat = asin( node.z );

                        double dSample = ( *testFunction )( dNodeLon, dNodeLat );

                        // Integrate
                        for( int i = 0; i < discOrder; i++ )
                        {
                            for( int j = 0; j < discOrder; j++ )
                            {

                                double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;

                                dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;

                                dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
                            }
                        }
                    }
                }
            }
        }

        // Divide by area
        if( fGLLIntegrate )
        {
            for( size_t i = 0; i < dVar.GetRows(); i++ )
            {
                dVar[i] /= dNodeArea[i];
            }
        }

        // Let us rearrange the data based on DoF ID specification
        if( ctx == Remapper::SourceMesh )
        {
            for( unsigned j = 0; j < entities.size(); j++ )
                for( int p = 0; p < discOrder; p++ )
                    for( int q = 0; q < discOrder; q++ )
                    {
                        const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
                        dVarMB[offsetDOF]   = dVar[col_dtoc_dofmap[offsetDOF]];
                    }
        }
        else
        {
            for( unsigned j = 0; j < entities.size(); j++ )
                for( int p = 0; p < discOrder; p++ )
                    for( int q = 0; q < discOrder; q++ )
                    {
                        const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
                        dVarMB[offsetDOF]   = dVar[row_dtoc_dofmap[offsetDOF]];
                    }
        }

        // Set the tag data
        rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
    }
    else
    {
        // assert( discOrder == 1 );
        if( discMethod == DiscretizationType_FV )
        {
            /* Compute an element-wise integral to store the sampled solution based on Quadrature
             * rules */
            // Resize the array
            dVar.Allocate( trmesh->faces.size() );

            std::vector< Node >& nodes = trmesh->nodes;

            // Loop through all Faces
            for( size_t i = 0; i < trmesh->faces.size(); i++ )
            {
                const Face& face = trmesh->faces[i];

                // Loop through all sub-triangles
                for( size_t j = 0; j < face.edges.size() - 2; j++ )
                {

                    const Node& node0 = nodes[face[0]];
                    const Node& node1 = nodes[face[j + 1]];
                    const Node& node2 = nodes[face[j + 2]];

                    // Triangle area
                    Face faceTri( 3 );
                    faceTri.SetNode( 0, face[0] );
                    faceTri.SetNode( 1, face[j + 1] );
                    faceTri.SetNode( 2, face[j + 2] );

                    double dTriangleArea = CalculateFaceArea( faceTri, nodes );

                    // Calculate the element average
                    double dTotalSample = 0.0;

                    // Loop through all quadrature points
                    for( int k = 0; k < TriQuadraturePoints; k++ )
                    {
                        Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
                                       TriQuadratureG[k][2] * node2.x,
                                   TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
                                       TriQuadratureG[k][2] * node2.y,
                                   TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
                                       TriQuadratureG[k][2] * node2.z );

                        double dMagnitude = node.Magnitude();
                        node.x /= dMagnitude;
                        node.y /= dMagnitude;
                        node.z /= dMagnitude;

                        double dLon = atan2( node.y, node.x );
                        if( dLon < 0.0 )
                        {
                            dLon += 2.0 * M_PI;
                        }
                        double dLat = asin( node.z );

                        double dSample = ( *testFunction )( dLon, dLat );

                        dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
                    }

                    dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
                }
            }
            rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
        }
        else /* discMethod == DiscretizationType_PCLOUD */
        {
            /* Get the coordinates of the vertices and sample the functionals accordingly */
            std::vector< Node >& nodes = trmesh->nodes;

            // Resize the array
            dVar.Allocate( nodes.size() );

            for( size_t j = 0; j < nodes.size(); j++ )
            {
                Node& node        = nodes[j];
                double dMagnitude = node.Magnitude();
                node.x /= dMagnitude;
                node.y /= dMagnitude;
                node.z /= dMagnitude;
                double dLon = atan2( node.y, node.x );
                if( dLon < 0.0 )
                {
                    dLon += 2.0 * M_PI;
                }
                double dLat = asin( node.z );

                double dSample = ( *testFunction )( dLon, dLat );
                dVar[j]        = dSample;
            }

            rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
        }
    }

    return moab::MB_SUCCESS;
}
moab::ErrorCode moab::TempestOnlineMap::fill_col_ids ( std::vector< int > &  ids_of_interest) [inline]

Definition at line 358 of file TempestOnlineMap.hpp.

References col_gdofmap, and MB_SUCCESS.

    {
        ids_of_interest.reserve( col_gdofmap.size() );
        // need to add 1
        for( auto it = col_gdofmap.begin(); it != col_gdofmap.end(); it++ )
            ids_of_interest.push_back( *it + 1 );
        return moab::MB_SUCCESS;
    }
moab::ErrorCode moab::TempestOnlineMap::fill_row_ids ( std::vector< int > &  ids_of_interest) [inline]

Definition at line 348 of file TempestOnlineMap.hpp.

References MB_SUCCESS, and row_gdofmap.

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

        return moab::MB_SUCCESS;
    }
moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights ( std::string  strInputType = "fv",
std::string  strOutputType = "fv",
const int  nPin = 1,
const int  nPout = 1,
bool  fBubble = false,
int  fMonotoneTypeID = 0,
bool  fVolumetric = false,
bool  fNoConservation = false,
bool  fNoCheck = false,
const std::string  srcDofTagName = "GLOBAL_ID",
const std::string  tgtDofTagName = "GLOBAL_ID",
const bool  fInputConcave = false,
const bool  fOutputConcave = false 
)

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

Definition at line 754 of file TempestOnlineMap.cpp.

References dbgprint, moab::error(), ErrorCode, LinearRemapFVtoGLL(), LinearRemapFVtoGLL_Volumetric(), MB_CHK_ERR, MB_CHK_SET_ERR, MB_SUCCESS, moab::Remapper::OverlapMesh, moab::DebugOutput::printf(), rank, and moab::DebugOutput::set_prefix().

Referenced by main().

{
    NcError error( NcError::silent_nonfatal );

    moab::DebugOutput dbgprint( std::cout, rank, 0 );
    dbgprint.set_prefix( "[TempestOnlineMap]: " );
    moab::ErrorCode rval;

    const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
    const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
    const bool m_bPointCloud       = m_bPointCloudSource || m_bPointCloudTarget;

    try
    {
        // Check command line parameters (data type arguments)
        STLStringHelper::ToLower( strInputType );
        STLStringHelper::ToLower( strOutputType );

        DiscretizationType eInputType;
        DiscretizationType eOutputType;
        int fNoCheckGlob = ( fNoCheck ? 1 : 0 );

        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() );
        }

        // Monotonicity flags
        int nMonotoneType = fMonotoneTypeID;
        m_bConserved      = !fNoConservation;
        m_iMonotonicity   = fMonotoneTypeID;
        m_eInputType      = eInputType;
        m_eOutputType     = eOutputType;

        m_nDofsPEl_Src =
            ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1 : nPin );
        m_nDofsPEl_Dest =
            ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1 : nPout );

        rval = SetDOFmapTags( srcDofTagName, tgtDofTagName );MB_CHK_ERR( rval );

        double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0;
        if( !m_bPointCloudSource )
        {
            // Calculate Input Mesh Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
            double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( fInputConcave );
            dTotalAreaInput            = dTotalAreaInput_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
            if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );

            // Input mesh areas
            m_meshInputCov->CalculateFaceAreas( fInputConcave );
        }

        if( !m_bPointCloudTarget )
        {
            // Calculate Output Mesh Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
            double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( fOutputConcave );
            dTotalAreaOutput            = dTotalAreaOutput_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
            if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
        }

        if( !m_bPointCloud )
        {
            // Verify that overlap mesh is in the correct order
            int ixSourceFaceMax = ( -1 );
            int ixTargetFaceMax = ( -1 );

            if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() )
            {
                _EXCEPTIONT( "Invalid overlap mesh:\n"
                             "    Possible mesh file corruption?" );
            }

            for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ )
            {
                if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
                    ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;

                if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
                    ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
            }

            // Check for forward correspondence in overlap mesh
            if( m_meshInput->faces.size() - ixSourceFaceMax == 0 )
            {
                if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" );
            }
            else if( m_meshOutput->faces.size() - ixSourceFaceMax == 0 )
            {  // Check for reverse correspondence in overlap mesh
                if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" );

                // Reorder overlap mesh
                m_meshOverlap->ExchangeFirstAndSecondMesh();
            }
            // else
            // {   // No correspondence found
            //     _EXCEPTION4 ( "Invalid overlap mesh:\n"
            //                   "    No correspondence found with input and output meshes (%i,%i)
            //                   vs (%i,%i)", m_meshInputCov->faces.size(),
            //                   m_meshOutput->faces.size(), ixSourceFaceMax, ixTargetFaceMax );
            // }

            // Calculate Face areas
            if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
            double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false );
            double dTotalAreaOverlap     = dTotalAreaOverlap_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
#endif
            if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );

            // Correct areas to match the areas calculated in the overlap mesh
            // if (fCorrectAreas)
            {
                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.IsAttached() )
                {
                    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.IsAttached() )
                {
                    this->SetTargetMask( m_meshOutput->vecMask );
                }
            }

            // Partial cover
            if( fabs( dTotalAreaOutput - dTotalAreaInput ) > 1.0e-10 &&
                fabs( dTotalAreaOverlap - dTotalAreaInput ) > 1.0e-10 )
            {
                if( !fNoCheckGlob )
                {
                    dbgprint.printf( rank, "WARNING: Significant mismatch between overlap mesh area "
                                           "and input mesh area.\n  Automatically enabling --nocheck\n" );
                    fNoCheckGlob = 1;
                }
            }

            /*
                // Recalculate input mesh area from overlap mesh
                if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
                    dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
                    dbgprint.printf(0, "Recalculating source mesh areas\n");
                    dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
                    dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
                }
            */
        }

        // Finite volume input / Finite volume output
        if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
        {
            // Generate reverse node array and edge map
            if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
            if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );

            // Initialize coordinates for map
            this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
            this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );

            // Finite volume input / Finite element output
            rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType, false, NULL );MB_CHK_ERR( rval );

            // Construct remap
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
            LinearRemapFVtoFV_Tempest_MOAB( nPin );
        }
        else if( eInputType == DiscretizationType_FV )
        {
            DataArray3D< double > dataGLLJacobian;

            if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
            double dNumericalArea_loc =
                GenerateMetaData( *m_meshOutput, nPout, fBubble, dataGLLNodesDest, dataGLLJacobian );

            double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Allreduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 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, nPout, dataGLLNodesDest );

            // Generate the continuous Jacobian
            bool fContinuous = ( eOutputType == DiscretizationType_CGLL );

            if( eOutputType == DiscretizationType_CGLL )
            {
                GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
            }
            else
            {
                GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
            }

            // Generate reverse node array and edge map
            if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
            if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );

            // Finite volume input / Finite element output
            rval = this->SetDOFmapAssociation( eInputType, false, NULL, NULL, eOutputType,
                                               ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );

            // Generate remap weights
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );

            if( fVolumetric )
            {
                LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
                                               dataGLLJacobian, this->GetTargetAreas(), nPin, *this, nMonotoneType,
                                               fContinuous, fNoConservation );
            }
            else
            {
                LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
                                    this->GetTargetAreas(), nPin, *this, nMonotoneType, fContinuous, 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, nPin, fBubble, dataGLLNodesSrcCov, dataGLLJacobian );
                    GenerateMetaData( *m_meshInput, nPin, fBubble, 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, nPout, fBubble, dataGLLNodesDest, dataGLLJacobian );
                }
            }
            // else { /* Target is a point cloud dataset */ }

            // Finite volume input / Finite element output
            rval = this->SetDOFmapAssociation(
                eInputType, ( eInputType == DiscretizationType_CGLL ),
                ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrcCov ),
                ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? NULL : &dataGLLNodesSrc ), eOutputType,
                ( eOutputType == DiscretizationType_CGLL ), ( m_bPointCloudTarget ? NULL : &dataGLLNodesDest ) );MB_CHK_ERR( rval );

            // Construct remap
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
            rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
        }
        else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
        {
            DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;

            if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
            // double dNumericalAreaCov_loc =
            GenerateMetaData( *m_meshInputCov, nPin, fBubble, dataGLLNodesSrcCov, dataGLLJacobian );

            double dNumericalArea_loc =
                GenerateMetaData( *m_meshInput, nPin, fBubble, dataGLLNodesSrc, dataGLLJacobianSrc );

            // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area:
            // %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
            // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);

            double dNumericalArea = dNumericalArea_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Allreduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
#endif
            if( is_root )
            {
                dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );

                if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
                {
                    dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
                                        "numerical area and geometric area\n" );
                }
            }

            if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
            {
                _EXCEPTIONT( "Number of element does not match between metadata and "
                             "input mesh" );
            }

            // Initialize coordinates for map
            this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, nPin, dataGLLNodesSrcCov );
            this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );

            // Generate the continuous Jacobian for input mesh
            bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );

            if( eInputType == DiscretizationType_CGLL )
            {
                GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
            }
            else
            {
                GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
            }

            // Finite element input / Finite volume output
            rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
                                               &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, false, NULL );MB_CHK_ERR( rval );

            // Generate remap
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );

            if( fVolumetric )
            {
                _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
                             "GLL input mesh" );
            }

            LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
                                         fNoConservation );
        }
        else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
        {
            DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
            DataArray3D< double > dataGLLJacobianOut;

            // Input metadata
            if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
            double dNumericalAreaIn_loc =
                GenerateMetaData( *m_meshInputCov, nPin, fBubble, dataGLLNodesSrcCov, dataGLLJacobianIn );

            double dNumericalAreaSrc_loc =
                GenerateMetaData( *m_meshInput, nPin, fBubble, dataGLLNodesSrc, dataGLLJacobianSrc );

            assert( dNumericalAreaIn_loc >= dNumericalAreaSrc_loc );

            double dNumericalAreaIn = dNumericalAreaSrc_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Allreduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
#endif
            if( is_root ) dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn );

            if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
            {
                dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
                                    "numerical area and geometric area\n" );
            }

            // Output metadata
            if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
            double dNumericalAreaOut_loc =
                GenerateMetaData( *m_meshOutput, nPout, fBubble, dataGLLNodesDest, dataGLLJacobianOut );

            double dNumericalAreaOut = dNumericalAreaOut_loc;
#ifdef MOAB_HAVE_MPI
            if( m_pcomm )
                MPI_Allreduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
#endif
            if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut );

            if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
            {
                if( is_root )
                    dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh "
                                        "numerical area and geometric area\n" );
            }

            // Initialize coordinates for map
            this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, nPin, dataGLLNodesSrcCov );
            this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, nPout, dataGLLNodesDest );

            // Generate the continuous Jacobian for input mesh
            bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );

            if( eInputType == DiscretizationType_CGLL )
            {
                GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
            }
            else
            {
                GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
            }

            // Generate the continuous Jacobian for output mesh
            bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );

            if( eOutputType == DiscretizationType_CGLL )
            {
                GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
            }
            else
            {
                GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
            }

            // Input Finite Element to Output Finite Element
            rval = this->SetDOFmapAssociation( eInputType, ( eInputType == DiscretizationType_CGLL ),
                                               &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType,
                                               ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );

            // Generate remap
            if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );

            LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
                                       this->GetTargetAreas(), nPin, nPout, nMonotoneType, fContinuousIn,
                                       fContinuousOut, fNoConservation );
        }
        else
        {
            _EXCEPTIONT( "Not implemented" );
        }

#ifdef MOAB_HAVE_EIGEN3
        copy_tempest_sparsemat_to_eigen3();
#endif

        int fNoCheckLoc = fNoCheckGlob;
#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" );
        }

        // Let us see if we need to perform consistency/conservation checks
        if( m_pcomm ) MPI_Allreduce( &fNoCheckLoc, &fNoCheckGlob, 1, MPI_INT, MPI_MAX, m_pcomm->comm() );
#endif
        // Verify consistency, conservation and monotonicity, globally
        if( !fNoCheckLoc )
        {
            if( is_root ) dbgprint.printf( 0, "Verifying map" );
            this->IsConsistent( 1.0e-8 );
            if( !fNoConservation ) this->IsConservative( 1.0e-8 );

            if( nMonotoneType != 0 )
            {
                this->IsMonotone( 1.0e-12 );
            }
        }
    }
    catch( Exception& e )
    {
        dbgprint.printf( 0, "%s", e.ToString().c_str() );
        return ( moab::MB_FAILURE );
    }
    catch( ... )
    {
        return ( moab::MB_FAILURE );
    }
    return moab::MB_SUCCESS;
}
int moab::TempestOnlineMap::GetColGlobalDoF ( int  localID) const [inline]

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

Definition at line 450 of file TempestOnlineMap.hpp.

{
    return col_gdofmap[localColID];
}

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

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

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

Definition at line 468 of file TempestOnlineMap.hpp.

{
    return m_nDofsPEl_Dest;
}
const DataArray1D< double >& moab::TempestOnlineMap::GetGlobalSourceAreas ( ) const

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

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

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

int moab::TempestOnlineMap::GetIndexOfColGlobalDoF ( int  globalColDoF) const [inline]

Get the index of globaColDoF.

Definition at line 455 of file TempestOnlineMap.hpp.

{
    return globalColDoF + 1;  // temporary
}
int moab::TempestOnlineMap::GetIndexOfRowGlobalDoF ( int  globalRowDoF) const [inline]

Get the index of globaRowDoF.

Definition at line 444 of file TempestOnlineMap.hpp.

{
    return globalRowDoF + 1;
}
int moab::TempestOnlineMap::GetRowGlobalDoF ( int  localID) const [inline]

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

Definition at line 439 of file TempestOnlineMap.hpp.

{
    return row_gdofmap[localRowID];
}

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

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

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

Definition at line 461 of file TempestOnlineMap.hpp.

{
    return m_nDofsPEl_Src;
}
int moab::TempestOnlineMap::IsConservative ( double  dTolerance) [virtual]

Determine if the map is conservative.

Definition at line 1387 of file TempestOnlineMap.cpp.

References ierr, rank, and size.

{
#ifndef MOAB_HAVE_MPI

    return OfflineMap::IsConservative( dTolerance );

#else
    // return OfflineMap::IsConservative(dTolerance);

    int ierr;
    // Get map entries
    DataArray1D< int > dataRows;
    DataArray1D< int > dataCols;
    DataArray1D< double > dataEntries;
    const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
    const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();

    // Calculate column sums
    std::vector< int > dColumnsUnique;
    std::vector< double > dColumnSums;

    int nColumns = m_mapRemap.GetColumns();
    m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
    dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
    dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );

    for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
    {
        dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];

        assert( dataCols[i] < m_nTotDofs_SrcCov );

        // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
        int colGID = this->GetColGlobalDoF( dataCols[i] );  // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
        // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
        dColumnsUnique[dataCols[i]] = colGID;

        // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
        // std::endl;
    }

    int rootProc = 0;
    std::vector< int > nElementsInProc;
    const int nDATA = 3;
    if( !rank ) nElementsInProc.resize( size * nDATA );
    int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
    ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;

    int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
    std::vector< int > dColumnIndices;
    std::vector< double > dColumnSourceAreas;
    std::vector< double > dColumnSumsTotal;
    std::vector< int > displs, rcount;
    if( rank == rootProc )
    {
        displs.resize( size + 1, 0 );
        rcount.resize( size, 0 );
        int gsum = 0;
        for( int ir = 0; ir < size; ++ir )
        {
            nTotVals += nElementsInProc[ir * nDATA];
            nTotColumns += nElementsInProc[ir * nDATA + 1];
            nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];

            displs[ir] = gsum;
            rcount[ir] = nElementsInProc[ir * nDATA + 1];
            gsum += rcount[ir];

            // printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns,
            // displs[ir], rcount[ir], gsum);
        }

        dColumnIndices.resize( nTotColumns, -1 );
        dColumnSumsTotal.resize( nTotColumns, 0.0 );
        // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
    }

    // Gather all ColumnSums to root process and accumulate
    // We expect that the sums of all columns equate to 1.0 within user specified tolerance
    // Need to do a gatherv here since different processes have different number of elements
    // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
    // MPI_SUM, 0, m_pcomm->comm());
    ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
                        displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;
    ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
                        displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;
    // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
    // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
    // MPI_SUCCESS ) return -1;

    // Clean out unwanted arrays now
    dColumnSums.clear();
    dColumnsUnique.clear();

    // Verify all column sums equal the input Jacobian
    int fConservative = 0;
    if( rank == rootProc )
    {
        displs[size] = ( nTotColumns );
        // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
        std::map< int, double > dColumnSumsOnRoot;
        // std::map<int, double> dColumnSourceAreasOnRoot;
        for( int ir = 0; ir < size; ir++ )
        {
            for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
            {
                if( dColumnIndices[ips] < 0 ) continue;
                // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
                assert( dColumnIndices[ips] < nTotColumnsUnq );
                dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips];  // / dColumnSourceAreas[ips];
                // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
                // dColumnSourceAreas[ dColumnIndices[ips] ]
            }
        }

        for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
        {
            // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
            if( fabs( it->second - 1.0 ) > dTolerance )
            {
                fConservative++;
                Announce( "TempestOnlineMap is not conservative in column "
                          // "%i (%1.15e)", it->first, it->second );
                          "%i (%1.15e)",
                          it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
            }
        }
    }

    // TODO: Just do a broadcast from root instead of a reduction
    ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;

    return fConservative;
#endif
}
int moab::TempestOnlineMap::IsConsistent ( double  dTolerance) [virtual]

Determine if the map is first-order accurate.

Definition at line 1341 of file TempestOnlineMap.cpp.

References ierr.

Referenced by read_map_from_disk().

{
#ifndef MOAB_HAVE_MPI

    return OfflineMap::IsConsistent( dTolerance );

#else

    // Get map entries
    DataArray1D< int > dataRows;
    DataArray1D< int > dataCols;
    DataArray1D< double > dataEntries;

    // Calculate row sums
    DataArray1D< double > dRowSums;
    m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
    dRowSums.Allocate( m_mapRemap.GetRows() );

    for( unsigned i = 0; i < dataRows.GetRows(); i++ )
    {
        dRowSums[dataRows[i]] += dataEntries[i];
    }

    // Verify all row sums are equal to 1
    int fConsistent = 0;
    for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
    {
        if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
        {
            fConsistent++;
            int rowGID = row_gdofmap[i];
            Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
        }
    }

    int ierr;
    int fConsistentGlobal = 0;
    ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;

    return fConsistentGlobal;
#endif
}
int moab::TempestOnlineMap::IsMonotone ( double  dTolerance) [virtual]

Determine if the map is monotone.

Definition at line 1529 of file TempestOnlineMap.cpp.

References ierr.

Referenced by read_map_from_disk().

{
#ifndef MOAB_HAVE_MPI

    return OfflineMap::IsMonotone( dTolerance );

#else

    // Get map entries
    DataArray1D< int > dataRows;
    DataArray1D< int > dataCols;
    DataArray1D< double > dataEntries;

    m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );

    // Verify all entries are in the range [0,1]
    int fMonotone = 0;
    for( unsigned i = 0; i < dataRows.GetRows(); i++ )
    {
        if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
        {
            fMonotone++;

            Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
        }
    }

    int ierr;
    int fMonotoneGlobal = 0;
    ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
    if( ierr != MPI_SUCCESS ) return -1;

    return fMonotoneGlobal;
#endif
}

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

Definition at line 130 of file TempestLinearRemap.cpp.

References BuildFitArray(), BuildIntegrationArray(), dbgprint, GetAdjacentFaceVectorByEdge(), InvertFitArray_Corrected(), moab::DebugOutput::printf(), rank, and moab::DebugOutput::set_prefix().

{
    // Order of triangular quadrature rule
    const int TriQuadRuleOrder = 4;

    // Verify ReverseNodeArray has been calculated
    if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
    {
        _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInput" );
    }

    // Triangular quadrature rule
    TriangularQuadratureRule triquadrule( TriQuadRuleOrder );

    // Number of coefficients needed at this order
#ifdef RECTANGULAR_TRUNCATION
    int nCoefficients = nOrder * nOrder;
#endif
#ifdef TRIANGULAR_TRUNCATION
    int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
#endif

    // Number of faces you need
    const int nRequiredFaceSetSize = nCoefficients;

    // Fit weight exponent
    const int nFitWeightsExponent = nOrder + 2;

    // Announcemnets
    moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
    dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
    if( is_root )
    {
        dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
        dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
        dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
        dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
        dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
    }

    // Current overlap face
    int ixOverlap                  = 0;
    const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;

    // Loop through all faces on m_meshInput
    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        // Output every 1000 elements
        if( ixFirst % outputFrequency == 0 && is_root )
        {
            dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
        }

        // Find the set of Faces that overlap faceFirst
        int ixOverlapBegin    = ixOverlap;
        unsigned ixOverlapEnd = ixOverlapBegin;

        for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
        {
            if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 )
            {
                break;
            }
        }

        unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
        // if ( is_root ) Announce ( "Element %i / %i :: [%i, %i]", ixFirst,
        // m_meshInputCov->faces.size(), ixOverlapBegin, ixOverlapEnd );

        if( nOverlapFaces == 0 ) continue;

        // Build integration array
        DataArray2D< double > dIntArray;

        BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
                               nOrder, dIntArray );

        // Set of Faces to use in building the reconstruction and associated
        // distance metric.
        AdjacentFaceVector vecAdjFaces;

        GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );

        // Number of adjacent Faces
        int nAdjFaces = vecAdjFaces.size();

        // Determine the conservative constraint equation
        DataArray1D< double > dConstraint( nCoefficients );

        double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];

        for( int p = 0; p < nCoefficients; p++ )
        {
            for( unsigned j = 0; j < nOverlapFaces; j++ )
            {
                dConstraint[p] += dIntArray[p][j];
            }
            dConstraint[p] /= dFirstArea;
        }

        // Build the fit array from the integration operator
        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
        InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );

        // Multiply integration array and fit array
        DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );

        for( int i = 0; i < nAdjFaces; i++ )
        {
            for( unsigned j = 0; j < nOverlapFaces; j++ )
            {
                for( int k = 0; k < nCoefficients; k++ )
                {
                    dComposedArray[i][j] += dIntArray[k][j] * dFitArrayPlus[i][k];
                }
            }
        }

        // Put composed array into map
        for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
        {
            for( unsigned j = 0; j < nOverlapFaces; j++ )
            {
                int& ixFirstFaceLoc  = vecAdjFaces[i].first;
                int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
                // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
                // ixFirstFaceLoc); int ixSecondFaceGlob =
                // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);

                // signal to not participate, because it is a ghost target
                if( ixSecondFaceLoc < 0 ) continue;  // do not do anything

                m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
                    dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
            }
        }

        // Increment the current overlap index
        ixOverlap += nOverlapFaces;
    }

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

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

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

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

Definition at line 891 of file TempestLinearRemap.cpp.

References center(), dbgprint, ForceIntArrayConsistencyConservation(), moab::DebugOutput::printf(), moab::DebugOutput::set_prefix(), and t.

{
    // Triangular quadrature rule
    TriangularQuadratureRule triquadrule( 8 );

    const DataArray2D< double >& dG = triquadrule.GetG();
    const DataArray1D< double >& dW = triquadrule.GetW();

    // Get SparseMatrix represntation of the OfflineMap
    SparseMatrix< double >& smatMap = this->GetSparseMatrix();

    // Sample coefficients
    DataArray2D< double > dSampleCoeffIn( nPin, nPin );
    DataArray2D< double > dSampleCoeffOut( nPout, nPout );

    // Announcemnets
    moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
    dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
    if( is_root )
    {
        dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
        dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
        dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
    }

    // Build the integration array for each element on m_meshOverlap
    DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );

    // Number of overlap Faces per source Face
    DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );

    int ixOverlap = 0;
    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        // Determine how many overlap Faces and triangles are present
        int nOverlapFaces    = 0;
        size_t ixOverlapTemp = ixOverlap;
        for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
        {
            // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
            if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
            {
                break;
            }

            nOverlapFaces++;
        }

        nAllOverlapFaces[ixFirst] = nOverlapFaces;

        // Increment the current overlap index
        ixOverlap += nAllOverlapFaces[ixFirst];
    }

    // Geometric area of each output node
    DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );

    // Area of each overlap element in the output basis
    DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );

    // Loop through all faces on m_meshInput
    ixOverlap                      = 0;
    const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;

    if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );

    // generic triangle used for area computation, for triangles around the center of overlap face;
    // used for overlap faces with more than 4 edges;
    // nodes array will be set for each triangle;
    // these triangles are not part of the mesh structure, they are just temporary during
    //   aforementioned decomposition.
    Face faceTri( 3 );
    NodeVector nodes( 3 );
    faceTri.SetNode( 0, 0 );
    faceTri.SetNode( 1, 1 );
    faceTri.SetNode( 2, 2 );

    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        // Announce computation progress
        if( ixFirst % outputFrequency == 0 && is_root )
        {
            dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
        }

        // 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++ )
    {
        // Announce computation progress
        if( ixFirst % outputFrequency == 0 && is_root )
        {
            dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
        }

        // Number of overlapping Faces and triangles
        int nOverlapFaces = nAllOverlapFaces[ixFirst];

        if( !nOverlapFaces ) continue;

        // Put composed array into map
        for( int j = 0; j < nOverlapFaces; j++ )
        {
            int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];

            // signal to not participate, because it is a ghost target
            if( ixSecondFace < 0 ) continue;  // do not do anything

            dRedistributedOp.Zero();
            for( int p = 0; p < nPin * nPin; p++ )
            {
                for( int s = 0; s < nPout * nPout; s++ )
                {
                    for( int t = 0; t < nPout * nPout; t++ )
                    {
                        dRedistributedOp[p][s] +=
                            dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
                    }
                }
            }

            int ixp = 0;
            for( int p = 0; p < nPin; p++ )
            {
                for( int q = 0; q < nPin; q++ )
                {
                    int ixFirstNode;
                    if( fContinuousIn )
                    {
                        ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
                    }
                    else
                    {
                        ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
                    }

                    int ixs = 0;
                    for( int s = 0; s < nPout; s++ )
                    {
                        for( int t = 0; t < nPout; t++ )
                        {
                            int ixSecondNode;
                            if( fContinuousOut )
                            {
                                ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;

                                if( !fNoConservation )
                                {
                                    smatMap( ixSecondNode, ixFirstNode ) +=
                                        dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
                                }
                                else
                                {
                                    smatMap( ixSecondNode, ixFirstNode ) +=
                                        dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
                                }
                            }
                            else
                            {
                                ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;

                                if( !fNoConservation )
                                {
                                    smatMap( ixSecondNode, ixFirstNode ) +=
                                        dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
                                }
                                else
                                {
                                    smatMap( ixSecondNode, ixFirstNode ) +=
                                        dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
                                }
                            }

                            ixs++;
                        }
                    }

                    ixp++;
                }
            }
        }

        // Increment the current overlap index
        ixOverlap += nOverlapFaces;
    }

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

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

Definition at line 1433 of file TempestLinearRemap.cpp.

References dbgprint, moab::DebugOutput::printf(), moab::DebugOutput::set_prefix(), and t.

{
    // Gauss-Lobatto quadrature within Faces
    DataArray1D< double > dGL;
    DataArray1D< double > dWL;

    GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );

    // Get SparseMatrix represntation of the OfflineMap
    SparseMatrix< double >& smatMap = this->GetSparseMatrix();

    // Sample coefficients
    DataArray2D< double > dSampleCoeffIn( nPin, nPin );

    // Announcemnets
    moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
    dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
    if( is_root )
    {
        dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
        dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
        dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
    }

    // Number of overlap Faces per source Face
    DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );

    int ixOverlap = 0;

    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        size_t ixOverlapTemp = ixOverlap;
        for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
        {
            // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];

            if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;

            nAllOverlapFaces[ixFirst]++;
        }

        // Increment the current overlap index
        ixOverlap += nAllOverlapFaces[ixFirst];
    }

    // Number of times this point was found
    DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );

    ixOverlap                      = 0;
    const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;

    // Loop through all faces on m_meshInputCov
    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        // Announce computation progress
        if( ixFirst % outputFrequency == 0 && is_root )
        {
            dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
        }

        // Quantities from the First Mesh
        const Face& faceFirst = m_meshInputCov->faces[ixFirst];

        const NodeVector& nodesFirst = m_meshInputCov->nodes;

        // Number of overlapping Faces and triangles
        int nOverlapFaces = nAllOverlapFaces[ixFirst];

        // Loop through all Overlap Faces
        for( int i = 0; i < nOverlapFaces; i++ )
        {

            // Quantities from the Second Mesh
            int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];

            // signal to not participate, because it is a ghost target
            if( ixSecond < 0 ) continue;  // do not do anything

            const NodeVector& nodesSecond = m_meshOutput->nodes;

            const Face& faceSecond = m_meshOutput->faces[ixSecond];

            // Loop through all nodes on the second face
            for( int s = 0; s < nPout; s++ )
            {
                for( int t = 0; t < nPout; t++ )
                {
                    size_t ixSecondNode;
                    if( fContinuousOut )
                    {
                        ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
                    }
                    else
                    {
                        ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
                    }

                    if( ixSecondNode >= fSecondNodeFound.GetRows() )
                    {
                        _EXCEPTIONT( "Logic error" );
                    }

                    // Check if this node has been found already
                    if( fSecondNodeFound[ixSecondNode] )
                    {
                        continue;
                    }

                    // Check this node
                    Node node;
                    Node dDx1G;
                    Node dDx2G;

                    ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );

                    // Find the components of this quadrature point in the basis
                    // of the first Face.
                    double dAlphaIn;
                    double dBetaIn;

                    ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );

                    // Check if this node is within the first Face
                    if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
                        ( dBetaIn > 1.0 + 1.0e-10 ) )
                    {
                        continue;
                    }

                    // Node is within the overlap region, mark as found
                    fSecondNodeFound[ixSecondNode] = true;

                    // Sample the First finite element at this point
                    SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );

                    // Add to map
                    for( int p = 0; p < nPin; p++ )
                    {
                        for( int q = 0; q < nPin; q++ )
                        {
                            int ixFirstNode;
                            if( fContinuousIn )
                            {
                                ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
                            }
                            else
                            {
                                ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
                            }

                            smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
                        }
                    }
                }
            }
        }

        // Increment the current overlap index
        ixOverlap += nOverlapFaces;
    }

    // Check for missing samples
    for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
    {
        if( !fSecondNodeFound[i] )
        {
            _EXCEPTION1( "Can't sample point %i", i );
        }
    }

    return;
}
moab::ErrorCode moab::TempestOnlineMap::LinearRemapNN_MOAB ( bool  use_GID_matching = false,
bool  strict_check = false 
) [private]

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

Definition at line 67 of file TempestLinearRemap.cpp.

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

{
    /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov)  */

#ifdef VVERBOSE
    {
        std::ofstream output_file( "rowcolindices.txt", std::ios::out );
        output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
                    << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
        output_file << "Rows \n";
        for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
            output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
        output_file << "Cols \n";
        for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
            output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
        output_file.flush();  // required here
        output_file.close();
    }
#endif

    if( use_GID_matching )
    {
        std::map< unsigned, unsigned > src_gl;
        for( unsigned it = 0; it < col_gdofmap.size(); ++it )
            src_gl[col_gdofmap[it]] = it;

        std::map< unsigned, unsigned >::iterator iter;
        for( unsigned it = 0; it < row_gdofmap.size(); ++it )
        {
            unsigned row = row_gdofmap[it];
            iter         = src_gl.find( row );
            if( strict_check && iter == src_gl.end() )
            {
                std::cout << "Searching for global target DOF " << row
                          << " but could not find correspondence in source mesh.\n";
                assert( false );
            }
            else if( iter == src_gl.end() )
            {
                continue;
            }
            else
            {
                unsigned icol = src_gl[row];
                unsigned irow = it;

                // Set the permutation matrix in local space
                m_mapRemap( irow, icol ) = 1.0;
            }
        }

        return moab::MB_SUCCESS;
    }
    else
    {
        /* Create a Kd-tree to perform local queries to find nearest neighbors */

        return moab::MB_SUCCESS;
    }
}
void moab::TempestOnlineMap::LinearRemapSE0_Tempest_MOAB ( const DataArray3D< int > &  dataGLLNodes,
const DataArray3D< double > &  dataGLLJacobian 
) [private]

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

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

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

Definition at line 470 of file TempestLinearRemap.cpp.

References center(), dbgprint, ForceConsistencyConservation3(), l, n, moab::DebugOutput::printf(), and moab::DebugOutput::set_prefix().

{
    // Order of the polynomial interpolant
    int nP = dataGLLNodes.GetRows();

    // Order of triangular quadrature rule
    const int TriQuadRuleOrder = 4;

    // Triangular quadrature rule
    TriangularQuadratureRule triquadrule( TriQuadRuleOrder );

    int TriQuadraturePoints = triquadrule.GetPoints();

    const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();

    const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();

    // Sample coefficients
    DataArray2D< double > dSampleCoeff( nP, nP );

    // GLL Quadrature nodes on quadrilateral elements
    DataArray1D< double > dG;
    DataArray1D< double > dW;
    GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );

    // Announcemnets
    moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
    dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
    if( is_root )
    {
        dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
        dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
        dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
    }

    // Get SparseMatrix represntation of the OfflineMap
    SparseMatrix< double >& smatMap = this->GetSparseMatrix();

    // NodeVector from m_meshOverlap
    const NodeVector& nodesOverlap = m_meshOverlap->nodes;
    const NodeVector& nodesFirst   = m_meshInputCov->nodes;

    // Vector of source areas
    DataArray1D< double > vecSourceArea( nP * nP );

    DataArray1D< double > vecTargetArea;
    DataArray2D< double > dCoeff;

#ifdef VERBOSE
    std::stringstream sstr;
    sstr << "remapdata_" << rank << ".txt";
    std::ofstream output_file( sstr.str() );
#endif

    // Current Overlap Face
    int ixOverlap                  = 0;
    const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;

    // generic triangle used for area computation, for triangles around the center of overlap face;
    // used for overlap faces with more than 4 edges;
    // nodes array will be set for each triangle;
    // these triangles are not part of the mesh structure, they are just temporary during
    //   aforementioned decomposition.
    Face faceTri( 3 );
    NodeVector nodes( 3 );
    faceTri.SetNode( 0, 0 );
    faceTri.SetNode( 1, 1 );
    faceTri.SetNode( 2, 2 );

    // Loop over all input Faces
    for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
    {
        const Face& faceFirst = m_meshInputCov->faces[ixFirst];

        if( faceFirst.edges.size() != 4 )
        {
            _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
        }

        // Announce computation progress
        if( ixFirst % outputFrequency == 0 && is_root )
        {
            dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
        }

        // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
        // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
        // However, the relation with MOAB and Tempest will go out of the roof

        // Determine how many overlap Faces and triangles are present
        int nOverlapFaces    = 0;
        size_t ixOverlapTemp = ixOverlap;
        for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
        {
            // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
            if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
            {
                break;
            }

            nOverlapFaces++;
        }

        // No overlaps
        if( nOverlapFaces == 0 )
        {
            continue;
        }

        // Allocate remap coefficients array for meshFirst Face
        DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );

        // Find the local remap coefficients
        for( int j = 0; j < nOverlapFaces; j++ )
        {
            const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
            if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 )  // machine precision
            {
                Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
                          m_meshOverlap->vecFaceArea[ixOverlap + j] );
                int n = faceOverlap.edges.size();
                Announce( "Number nodes: %d", n );
                for( int k = 0; k < n; k++ )
                {
                    Node nd = nodesOverlap[faceOverlap[k]];
                    Announce( "Node %d  %d  : %1.10e  %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
                }
                continue;
            }

            // #ifdef VERBOSE
            // if ( is_root )
            //     Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
            //     m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
            //     m_meshOverlap->vecFaceArea[ixOverlap + j] );
            // #endif

            int nbEdges = faceOverlap.edges.size();
            int nOverlapTriangles = 1;
            Node center; // not used if nbEdges == 3
            if (nbEdges > 3) { // decompose from center in this case
                nOverlapTriangles = nbEdges;
                for (int k = 0; k < nbEdges; k++) {
                    const Node &node = nodesOverlap[faceOverlap[k]];
                    center = center + node;
                }
                center = center / nbEdges;
                center = center.Normalized();// project back on sphere of radius 1
            }

            Node node0, node1, node2;
            double dTriangleArea;

            // Loop over all sub-triangles of this Overlap Face
            for( int k = 0; k < nOverlapTriangles; k++ )
            {
                if (nbEdges == 3) // will come here only once, nOverlapTriangles == 1 in this case
                {
                    node0 = nodesOverlap[faceOverlap[0]];
                    node1 = nodesOverlap[faceOverlap[1]];
                    node2 = nodesOverlap[faceOverlap[2]];
                    dTriangleArea = CalculateFaceArea(faceOverlap, nodesOverlap);
                }
                else // decompose polygon in triangles around the center
                {
                    node0 = center;
                    node1 = nodesOverlap[faceOverlap[k]];
                    int k1 = (k + 1) % nbEdges;
                    node2 = nodesOverlap[faceOverlap[k1]];
                    nodes[0] = center;
                    nodes[1] = node1;
                    nodes[2] = node2;
                    dTriangleArea = CalculateFaceArea(faceTri, nodes);
                }
                // Coordinates of quadrature Node
                for( int l = 0; l < TriQuadraturePoints; l++ )
                {
                    Node nodeQuadrature;
                    nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
                                       TriQuadratureG[l][2] * node2.x;

                    nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
                                       TriQuadratureG[l][2] * node2.y;

                    nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
                                       TriQuadratureG[l][2] * node2.z;

                    nodeQuadrature = nodeQuadrature.Normalized();

                    // Find components of quadrature point in basis
                    // of the first Face
                    double dAlpha;
                    double dBeta;

                    ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );

                    // Check inverse map value
                    if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
                        ( dBeta > 1.0 + 1.0e-13 ) )
                    {
                        _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
                                     "(%1.5e %1.5e)",
                                     j, l, dAlpha, dBeta );
                    }

                    // Sample the finite element at this point
                    SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );

                    // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
                    for( int p = 0; p < nP; p++ )
                    {
                        for( int q = 0; q < nP; q++ )
                        {
                            dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
                                                    m_meshOverlap->vecFaceArea[ixOverlap + j];
                        }
                    }
                }
            }
        }

#ifdef VERBOSE
        output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
        for( int j = 0; j < nOverlapFaces; j++ )
        {
            for( int p = 0; p < nP; p++ )
            {
                for( int q = 0; q < nP; q++ )
                {
                    output_file << dRemapCoeff[p][q][j] << " ";
                }
            }
        }
        output_file << std::endl;
#endif

        // Force consistency and conservation
        if( !fNoConservation )
        {
            double dTargetArea = 0.0;
            for( int j = 0; j < nOverlapFaces; j++ )
            {
                dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
            }

            for( int p = 0; p < nP; p++ )
            {
                for( int q = 0; q < nP; q++ )
                {
                    vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
                }
            }

            const double areaTolerance = 1e-10;
            // Source elements are completely covered by target volumes
            if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
            {
                vecTargetArea.Allocate( nOverlapFaces );
                for( int j = 0; j < nOverlapFaces; j++ )
                {
                    vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
                }

                dCoeff.Allocate( nOverlapFaces, nP * nP );

                for( int j = 0; j < nOverlapFaces; j++ )
                {
                    for( int p = 0; p < nP; p++ )
                    {
                        for( int q = 0; q < nP; q++ )
                        {
                            dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
                        }
                    }
                }

                // Target volumes only partially cover source elements
            }
            else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
            {
                double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;

                vecTargetArea.Allocate( nOverlapFaces + 1 );
                for( int j = 0; j < nOverlapFaces; j++ )
                {
                    vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
                }
                vecTargetArea[nOverlapFaces] = dExtraneousArea;

#ifdef VERBOSE
                Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
                          m_meshInputCov->vecFaceArea[ixFirst] );
#endif
                if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
                {
                    _EXCEPTIONT( "Partial element area exceeds total element area" );
                }

                dCoeff.Allocate( nOverlapFaces + 1, nP * nP );

                for( int j = 0; j < nOverlapFaces; j++ )
                {
                    for( int p = 0; p < nP; p++ )
                    {
                        for( int q = 0; q < nP; q++ )
                        {
                            dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
                        }
                    }
                }
                for( int p = 0; p < nP; p++ )
                {
                    for( int q = 0; q < nP; q++ )
                    {
                        dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
                    }
                }
                for( int j = 0; j < nOverlapFaces; j++ )
                {
                    for( int p = 0; p < nP; p++ )
                    {
                        for( int q = 0; q < nP; q++ )
                        {
                            dCoeff[nOverlapFaces][p * nP + q] -=
                                dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
                        }
                    }
                }
                for( int p = 0; p < nP; p++ )
                {
                    for( int q = 0; q < nP; q++ )
                    {
                        dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
                    }
                }

                // Source elements only partially cover target volumes
            }
            else
            {
                Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
                          m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
                _EXCEPTIONT( "Target grid must be a subset of source grid" );
            }

            ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
                                           /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );

            for( int j = 0; j < nOverlapFaces; j++ )
            {
                for( int p = 0; p < nP; p++ )
                {
                    for( int q = 0; q < nP; q++ )
                    {
                        dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
                    }
                }
            }
        }

#ifdef VERBOSE
        // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
        // for ( int j = 0; j < nOverlapFaces; j++ )
        // {
        //     for ( int p = 0; p < nP; p++ )
        //     {
        //         for ( int q = 0; q < nP; q++ )
        //         {
        //             output_file << dRemapCoeff[p][q][j] << " ";
        //         }
        //     }
        // }
        // output_file << std::endl;
#endif

        // Put these remap coefficients into the SparseMatrix map
        for( int j = 0; j < nOverlapFaces; j++ )
        {
            int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];

            // signal to not participate, because it is a ghost target
            if( ixSecondFace < 0 ) continue;  // do not do anything

            for( int p = 0; p < nP; p++ )
            {
                for( int q = 0; q < nP; q++ )
                {
                    if( fContinuousIn )
                    {
                        int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;

                        smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
                                                                m_meshOverlap->vecFaceArea[ixOverlap + j] /
                                                                m_meshOutput->vecFaceArea[ixSecondFace];
                    }
                    else
                    {
                        int ixFirstNode = ixFirst * nP * nP + p * nP + q;

                        smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
                                                                m_meshOverlap->vecFaceArea[ixOverlap + j] /
                                                                m_meshOutput->vecFaceArea[ixSecondFace];
                    }
                }
            }
        }
        // Increment the current overlap index
        ixOverlap += nOverlapFaces;
    }
#ifdef VERBOSE
    output_file.flush();  // required here
    output_file.close();
#endif

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

Generate the metadata associated with the offline map.

Read the OfflineMap from a NetCDF file.

Definition at line 1145 of file TempestOnlineMapIO.cpp.

References CHECK_EXCEPTION, moab::TupleList::enableWriteAccess(), moab::error(), moab::TupleList::get_n(), moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_SUCCESS, n, rank, size, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.

Referenced by read_map_from_disk().

{
    NcError error( NcError::silent_nonfatal );

    NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
    int nS = 0, nA = 0, nB = 0;
#ifdef MOAB_HAVE_PNETCDF
    // some variables will be used just in the case netcdfpar reader fails
    int ncfile = -1, ret = 0;
    int ndims, nvars, ngatts, unlimited;
#endif
#ifdef MOAB_HAVE_NETCDFPAR
    bool is_independent = true;
    ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
    // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
    NcFile ncMap( strSource, NcFile::ReadOnly );
#endif

#define CHECK_EXCEPTION( obj, type, varstr )                                                      \
    {                                                                                             \
        if( obj == NULL )                                                                         \
        {                                                                                         \
            _EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
        }                                                                                         \
    }

    // Read SparseMatrix entries

    if( ncMap.is_valid() )
    {
        NcDim* dimNS = ncMap.get_dim( "n_s" );
        CHECK_EXCEPTION( dimNS, "dimension", "n_s" );

        NcDim* dimNA = ncMap.get_dim( "n_a" );
        CHECK_EXCEPTION( dimNA, "dimension", "n_a" );

        NcDim* dimNB = ncMap.get_dim( "n_b" );
        CHECK_EXCEPTION( dimNB, "dimension", "n_b" );

        // store total number of nonzeros
        nS = dimNS->size();
        nA = dimNA->size();
        nB = dimNB->size();

        varRow = ncMap.get_var( "row" );
        CHECK_EXCEPTION( varRow, "variable", "row" );

        varCol = ncMap.get_var( "col" );
        CHECK_EXCEPTION( varCol, "variable", "col" );

        varS = ncMap.get_var( "S" );
        CHECK_EXCEPTION( varS, "variable", "S" );

#ifdef MOAB_HAVE_NETCDFPAR
        ncMap.enable_var_par_access( varRow, is_independent );
        ncMap.enable_var_par_access( varCol, is_independent );
        ncMap.enable_var_par_access( varS, is_independent );
#endif
    }
    else
    {
#ifdef MOAB_HAVE_PNETCDF
        // read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
        // why build wth pnetcdf without MPI ?
        // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
        ret = ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile ); ERR_PARNC( ret );  // bail out completely
        ret = ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited ); ERR_PARNC( ret );
        // find dimension ids for n_S
        int ins;
        ret = ncmpi_inq_dimid( ncfile, "n_s", &ins ); ERR_PARNC( ret );
        MPI_Offset leng;
        ret = ncmpi_inq_dimlen( ncfile, ins, &leng ); ERR_PARNC( ret );
        nS  = (int)leng;
        ret = ncmpi_inq_dimid( ncfile, "n_b", &ins ); ERR_PARNC( ret );
        ret = ncmpi_inq_dimlen( ncfile, ins, &leng ); ERR_PARNC( ret );
        nB = (int)leng;
#else
        _EXCEPTION1( "cannot read the file %s", strSource );
#endif
    }

    std::vector< int > rowOwnership;
    // if owned_dof_ids = NULL, use the default trivial partitioning scheme
    if( owned_dof_ids.size() == 0 )
    {
        // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
        rowOwnership.resize( size );
        int nGRowPerPart   = nB / size;
        int nGRowRemainder = nB % size;  // Keep the remainder in root
        rowOwnership[0]    = nGRowPerPart + nGRowRemainder;
        for( int ip = 1, roffset = rowOwnership[0]; ip < size; ++ip )
        {
            roffset += nGRowPerPart;
            rowOwnership[ip] = roffset;
        }
    }

    // Let us declare the map object for every process
    SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();

    int localSize   = nS / size;
    long offsetRead = rank * localSize;
    // leftovers on last rank
    if( rank == size - 1 )
    {
        localSize += nS % size;
    }

    std::vector< int > vecRow, vecCol;
    std::vector< double > vecS;
    vecRow.resize( localSize );
    vecCol.resize( localSize );
    vecS.resize( localSize );

    if( ncMap.is_valid() )
    {
        varRow->set_cur( (long)( offsetRead ) );
        varRow->get( &( vecRow[0] ), localSize );

        varCol->set_cur( (long)( offsetRead ) );
        varCol->get( &( vecCol[0] ), localSize );

        varS->set_cur( (long)( offsetRead ) );
        varS->get( &( vecS[0] ), localSize );

        ncMap.close();
    }
    else
    {
#ifdef MOAB_HAVE_PNETCDF
        // fill the local vectors with the variables from pnetcdf file; first inquire, then fill
        MPI_Offset start = (MPI_Offset)offsetRead;
        MPI_Offset count = (MPI_Offset)localSize;
        int varid;
        ret = ncmpi_inq_varid( ncfile, "S", &varid );  ERR_PARNC( ret );
        ret = ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] ); ERR_PARNC( ret );
        ret = ncmpi_inq_varid( ncfile, "row", &varid ); ERR_PARNC( ret );
        ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] ); ERR_PARNC( ret );
        ret = ncmpi_inq_varid( ncfile, "col", &varid ); ERR_PARNC( ret );
        ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] ); ERR_PARNC( ret );
    ret = ncmpi_close(ncfile); ERR_PARNC( ret );
#endif
    }

    // Now let us set the necessary global-to-local ID maps so that A*x operations
    // can be performed cleanly as if map was computed online
    row_dtoc_dofmap.clear();
    // row_dtoc_dofmap.reserve( nB / size );
    col_dtoc_dofmap.clear();
    rowMap.clear();
    colMap.clear();
    // col_dtoc_dofmap.reserve( 2 * nA / size );
    // row_dtoc_dofmap.resize( m_nTotDofs_Dest, UINT_MAX );
    // col_dtoc_dofmap.resize( m_nTotDofs_SrcCov, UINT_MAX );

#ifdef MOAB_HAVE_MPI
    // bother with tuple list only if size > 1
    // otherwise, just fill the sparse matrix
    if( size > 1 )
    {
        // send to
        moab::TupleList tl;
        unsigned numr = 1;                          //
        tl.initialize( 3, 0, 0, numr, localSize );  // to proc, row, col, value
        tl.enableWriteAccess();
        // populate
        for( int i = 0; i < localSize; i++ )
        {
            int rowval  = vecRow[i] - 1;  // dofs are 1 based in the file
            int colval  = vecCol[i] - 1;
            int to_proc = -1;
            //

            if( rowOwnership[0] > rowval )
                to_proc = 0;
            else
            {
                for( int ip = 1; ip < size; ++ip )
                {
                    if( rowOwnership[ip - 1] <= rowval && rowOwnership[ip] > rowval )
                    {
                        to_proc = ip;
                        break;
                    }
                }
            }

            int n               = tl.get_n();
            tl.vi_wr[3 * n]     = to_proc;
            tl.vi_wr[3 * n + 1] = rowval;
            tl.vi_wr[3 * n + 2] = colval;
            tl.vr_wr[n]         = vecS[i];

            tl.inc_n();
        }

        // now do the heavy communication
        ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl, 0 );

        int rindexMax = 0, cindexMax = 0;
        // populate the sparsematrix, using rowMap and colMap
        int n = tl.get_n();
        for( int i = 0; i < n; i++ )
        {
            int rindex, cindex;
            const int& vecRowValue = tl.vi_wr[3 * i + 1];
            const int& vecColValue = tl.vi_wr[3 * i + 2];

            std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
            if( riter == rowMap.end() )
            {
                rowMap[vecRowValue] = rindexMax;
                rindex              = rindexMax;
                row_gdofmap.push_back( vecRowValue );
                // row_dtoc_dofmap.push_back( vecRowValue );
                rindexMax++;
            }
            else
                rindex = riter->second;

            std::map< int, int >::iterator citer = colMap.find( vecColValue );
            if( citer == colMap.end() )
            {
                colMap[vecColValue] = cindexMax;
                cindex              = cindexMax;
                col_gdofmap.push_back( vecColValue );
                // col_dtoc_dofmap.push_back( vecColValue );
                cindexMax++;
            }
            else
                cindex = citer->second;

            sparseMatrix( rindex, cindex ) = tl.vr_wr[i];
        }
    }
    else
#endif
    {
        int rindexMax = 0, cindexMax = 0;

        for( int i = 0; i < nS; i++ )
        {
            int rindex, cindex;
            const int& vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
            const int& vecColValue = vecCol[i] - 1;

            std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
            if( riter == rowMap.end() )
            {
                rowMap[vecRowValue] = rindexMax;
                rindex              = rindexMax;
                row_gdofmap.push_back( vecRowValue );
                // row_dtoc_dofmap.push_back( vecRowValue );
                rindexMax++;
            }
            else
                rindex = riter->second;

            std::map< int, int >::iterator citer = colMap.find( vecColValue );
            if( citer == colMap.end() )
            {
                colMap[vecColValue] = cindexMax;
                cindex              = cindexMax;
                col_gdofmap.push_back( vecColValue );
                // col_dtoc_dofmap.push_back( vecColValue );
                cindexMax++;
            }
            else
                cindex = citer->second;

            sparseMatrix( rindex, cindex ) = vecS[i];
        }
    }

    m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
    m_nTotDofs_Dest   = sparseMatrix.GetRows();

#ifdef MOAB_HAVE_EIGEN3
    this->copy_tempest_sparsemat_to_eigen3();
#endif

    // Reset the source and target data first
    m_rowVector.setZero();
    m_colVector.setZero();

    return moab::MB_SUCCESS;
}
moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs ( std::vector< int > &  values_entities)

Definition at line 710 of file TempestOnlineMap.cpp.

References MB_SUCCESS.

{
    // col_gdofmap has global dofs , that should be in the list of values, such that
    // row_dtoc_dofmap[offsetDOF] = localDOF;
    // //  we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
    // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
    // form first inverse

    col_dtoc_dofmap.resize( values_entities.size() );
    for( int j = 0; j < (int)values_entities.size(); j++ )
    {
        if( colMap.find( values_entities[j] - 1 ) != colMap.end() )
            col_dtoc_dofmap[j] = colMap[values_entities[j] - 1];
        else
        {
            col_dtoc_dofmap[j] = -1;  // signal that this value should not be used in
            // std::cout <<"values_entities[j] -  1: " << values_entities[j] -  1 <<" at index j = " << j <<  " not
            // found in colMap \n";
        }
    }
    return moab::MB_SUCCESS;
}
moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs ( std::vector< int > &  values_entities)

Definition at line 733 of file TempestOnlineMap.cpp.

References MB_SUCCESS.

{
    // row_dtoc_dofmap = values_entities; // needs to point to local
    //  we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];

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

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

Definition at line 479 of file TempestOnlineMap.hpp.

{
    m_nDofsPEl_Dest = nt;
}
moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation ( DiscretizationType  srcType,
bool  isSrcContinuous,
DataArray3D< int > *  srcdataGLLNodes,
DataArray3D< int > *  srcdataGLLNodesSrc,
DiscretizationType  destType,
bool  isDestContinuous,
DataArray3D< int > *  tgtdataGLLNodes 
)

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

Definition at line 237 of file TempestOnlineMap.cpp.

References ErrorCode, l, MB_CHK_ERR, MB_SUCCESS, and rank.

{
    moab::ErrorCode rval;
    std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
    std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;

    // We are assuming that these are element based tags that are sized: np * np
    m_srcDiscType  = srcType;
    m_destDiscType = destType;

    bool vprint = is_root && false;

#ifdef VVERBOSE
    {
        src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
        locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
        tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
        rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );

        if( is_root )
        {
            {
                std::ofstream output_file( "sourcecov-gids-0.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
                    output_file << i << ", " << src_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
                m_nTotDofs_SrcCov = 0;
                if( isSrcContinuous )
                    dgll_cgll_covcol_ldofmap.resize(
                        m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
                for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
                {
                    for( int p = 0; p < m_nDofsPEl_Src; p++ )
                    {
                        for( int q = 0; q < m_nDofsPEl_Src; q++ )
                        {
                            const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
                            const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
                            if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
                            {
                                m_nTotDofs_SrcCov++;
                                dgll_cgll_covcol_ldofmap[localDOF] = true;
                            }
                            output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
                                        << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
                        }
                    }
                }
                output_file.flush();  // required here
                output_file.close();
                dgll_cgll_covcol_ldofmap.clear();
            }

            {
                std::ofstream output_file( "source-gids-0.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
                    output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
                m_nTotDofs_Src = 0;
                if( isSrcContinuous )
                    dgll_cgll_col_ldofmap.resize(
                        m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
                for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
                {
                    for( int p = 0; p < m_nDofsPEl_Src; p++ )
                    {
                        for( int q = 0; q < m_nDofsPEl_Src; q++ )
                        {
                            const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
                            const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
                            if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
                            {
                                m_nTotDofs_Src++;
                                dgll_cgll_col_ldofmap[localDOF] = true;
                            }
                            output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
                                        << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
                        }
                    }
                }
                output_file.flush();  // required here
                output_file.close();
                dgll_cgll_col_ldofmap.clear();
            }

            {
                std::ofstream output_file( "target-gids-0.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
                    output_file << i << ", " << tgt_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, GDOF, NDOF\n";
                m_nTotDofs_Dest = 0;

                for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
                {
                    output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
                                << m_nTotDofs_Dest << "\n";
                    m_nTotDofs_Dest++;
                }

                output_file.flush();  // required here
                output_file.close();
            }
        }
        else
        {
            {
                std::ofstream output_file( "sourcecov-gids-1.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
                    output_file << i << ", " << src_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
                m_nTotDofs_SrcCov = 0;
                if( isSrcContinuous )
                    dgll_cgll_covcol_ldofmap.resize(
                        m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
                for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
                {
                    for( int p = 0; p < m_nDofsPEl_Src; p++ )
                    {
                        for( int q = 0; q < m_nDofsPEl_Src; q++ )
                        {
                            const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
                            const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
                            if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
                            {
                                m_nTotDofs_SrcCov++;
                                dgll_cgll_covcol_ldofmap[localDOF] = true;
                            }
                            output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
                                        << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
                        }
                    }
                }
                output_file.flush();  // required here
                output_file.close();
                dgll_cgll_covcol_ldofmap.clear();
            }

            {
                std::ofstream output_file( "source-gids-1.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
                    output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
                m_nTotDofs_Src = 0;
                if( isSrcContinuous )
                    dgll_cgll_col_ldofmap.resize(
                        m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, false );
                for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
                {
                    for( int p = 0; p < m_nDofsPEl_Src; p++ )
                    {
                        for( int q = 0; q < m_nDofsPEl_Src; q++ )
                        {
                            const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
                            const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
                            if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
                            {
                                m_nTotDofs_Src++;
                                dgll_cgll_col_ldofmap[localDOF] = true;
                            }
                            output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
                                        << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
                        }
                    }
                }
                output_file.flush();  // required here
                output_file.close();
                dgll_cgll_col_ldofmap.clear();
            }

            {
                std::ofstream output_file( "target-gids-1.txt" );
                output_file << "I, GDOF\n";
                for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
                    output_file << i << ", " << tgt_soln_gdofs[i] << "\n";

                output_file << "ELEMID, IDOF, GDOF, NDOF\n";
                m_nTotDofs_Dest = 0;

                for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
                {
                    output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
                                << m_nTotDofs_Dest << "\n";
                    m_nTotDofs_Dest++;
                }

                output_file.flush();  // required here
                output_file.close();
            }
        }
    }
#endif

    // Now compute the mapping and store it for the covering mesh
    int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
    if( m_remapper->point_cloud_source )
    {
        assert( m_nDofsPEl_Src == 1 );
        col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
        col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
        src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
        srcTagSize = 1;
    }
    else
    {
        col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
        col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
        src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
    }

    // std::cout << "TOnlineMap: Process: " << rank << " and covering entities = [" <<
    // col_dofmap.size() << ", " << src_soln_gdofs.size() << "]\n"; MPI_Barrier(MPI_COMM_WORLD);

#ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
    unsigned maxSrcIndx = 0;

    // for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
    std::vector< int > locdofs( srcTagSize );
    std::map< Node, moab::EntityHandle > mapLocalMBNodes;
    double elcoords[3];
    for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
    {
        EntityHandle eh = m_remapper->m_covering_source_entities[iel];
        rval            = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
        Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
        mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
    }

    const NodeVector& nodes = m_remapper->m_covering_source->nodes;
    for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
    {
        const Face& face = m_remapper->m_covering_source->faces[j];

        Node centroid;
        centroid.x = centroid.y = centroid.z = 0.0;
        for( unsigned l = 0; l < face.edges.size(); ++l )
        {
            centroid.x += nodes[face[l]].x;
            centroid.y += nodes[face[l]].y;
            centroid.z += nodes[face[l]].z;
        }
        const double factor = 1.0 / face.edges.size();
        centroid.x *= factor;
        centroid.y *= factor;
        centroid.z *= factor;

        EntityHandle current_eh;
        if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
        {
            current_eh = mapLocalMBNodes[centroid];
        }

        rval = m_interface->tag_get_data( m_dofTagSrc, &current_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
        for( int p = 0; p < m_nDofsPEl_Src; p++ )
        {
            for( int q = 0; q < m_nDofsPEl_Src; q++ )
            {
                const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
                const int offsetDOF = p * m_nDofsPEl_Src + q;
                maxSrcIndx          = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
                std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
                          << ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
            }
        }
    }
#endif

    m_nTotDofs_SrcCov = 0;
    if( srcdataGLLNodes == NULL )
    { /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < col_gdofmap.size(); ++i )
        {
            assert( src_soln_gdofs[i] > 0 );
            col_gdofmap[i]     = src_soln_gdofs[i] - 1;
            col_dtoc_dofmap[i] = i;
            if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
            m_nTotDofs_SrcCov++;
        }
    }
    else
    {
        if( isSrcContinuous )
            dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
        // Put these remap coefficients into the SparseMatrix map
        for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
        {
            for( int p = 0; p < m_nDofsPEl_Src; p++ )
            {
                for( int q = 0; q < m_nDofsPEl_Src; q++ )
                {
                    const int localDOF  = ( *srcdataGLLNodes )[p][q][j] - 1;
                    const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
                    if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
                    {
                        m_nTotDofs_SrcCov++;
                        dgll_cgll_covcol_ldofmap[localDOF] = true;
                    }
                    if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
                    assert( src_soln_gdofs[offsetDOF] > 0 );
                    col_gdofmap[localDOF]      = src_soln_gdofs[offsetDOF] - 1;
                    col_dtoc_dofmap[offsetDOF] = localDOF;
                    if( vprint )
                        std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", "
                                  << localDOF << ", " << col_gdofmap[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
                }
            }
        }
    }

    if( m_remapper->point_cloud_source )
    {
        assert( m_nDofsPEl_Src == 1 );
        srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
        srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
        locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
    }
    else
    {
        srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
        srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
        locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
    }

    // Now compute the mapping and store it for the original source mesh
    m_nTotDofs_Src = 0;
    if( srcdataGLLNodesSrc == NULL )
    { /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
        {
            assert( locsrc_soln_gdofs[i] > 0 );
            srccol_gdofmap[i]     = locsrc_soln_gdofs[i] - 1;
            srccol_dtoc_dofmap[i] = i;
            m_nTotDofs_Src++;
        }
    }
    else
    {
        if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
        // Put these remap coefficients into the SparseMatrix map
        for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
        {
            for( int p = 0; p < m_nDofsPEl_Src; p++ )
            {
                for( int q = 0; q < m_nDofsPEl_Src; q++ )
                {
                    const int localDOF  = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
                    const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
                    if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
                    {
                        m_nTotDofs_Src++;
                        dgll_cgll_col_ldofmap[localDOF] = true;
                    }
                    if( !isSrcContinuous ) m_nTotDofs_Src++;
                    assert( locsrc_soln_gdofs[offsetDOF] > 0 );
                    srccol_gdofmap[localDOF]      = locsrc_soln_gdofs[offsetDOF] - 1;
                    srccol_dtoc_dofmap[offsetDOF] = localDOF;
                }
            }
        }
    }

    int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
    if( m_remapper->point_cloud_target )
    {
        assert( m_nDofsPEl_Dest == 1 );
        row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
        row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
        tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
        tgtTagSize = 1;
    }
    else
    {
        row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
        row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
        tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
        rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
    }

    // Now compute the mapping and store it for the target mesh
    // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
    m_nTotDofs_Dest = 0;
    if( tgtdataGLLNodes == NULL )
    { /* we only have a mapping for elements as DoFs */
        for( unsigned i = 0; i < row_gdofmap.size(); ++i )
        {
            assert( tgt_soln_gdofs[i] > 0 );
            row_gdofmap[i]     = tgt_soln_gdofs[i] - 1;
            row_dtoc_dofmap[i] = i;
            if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
            m_nTotDofs_Dest++;
        }
    }
    else
    {
        if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
        // Put these remap coefficients into the SparseMatrix map
        for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
        {
            for( int p = 0; p < m_nDofsPEl_Dest; p++ )
            {
                for( int q = 0; q < m_nDofsPEl_Dest; q++ )
                {
                    const int localDOF  = ( *tgtdataGLLNodes )[p][q][j] - 1;
                    const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
                    if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
                    {
                        m_nTotDofs_Dest++;
                        dgll_cgll_row_ldofmap[localDOF] = true;
                    }
                    if( !isTgtContinuous ) m_nTotDofs_Dest++;
                    assert( tgt_soln_gdofs[offsetDOF] > 0 );
                    row_gdofmap[localDOF]      = tgt_soln_gdofs[offsetDOF] - 1;
                    row_dtoc_dofmap[offsetDOF] = localDOF;
                    if( vprint )
                        std::cout << "Row: " << m_remapper->lid_to_gid_tgt[j] << ", " << offsetDOF << ", " << localDOF
                                  << ", " << row_gdofmap[offsetDOF] << ", " << m_nTotDofs_Dest << "\n";
                }
            }
        }
    }

    // Let us also allocate the local representation of the sparse matrix
#if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
    if( vprint )
    {
        std::cout << "[" << rank << "]"
                  << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size() << ", col = " << m_nTotDofs_Src
                  << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
        // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
    }
#endif

    // check monotonicity of row_gdofmap and col_gdofmap
#ifdef CHECK_INCREASING_DOF
    for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
    {
        if( row_gdofmap[i] > row_gdofmap[i + 1] )
            std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
                      << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
    }
    for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
    {
        if( col_gdofmap[i] > col_gdofmap[i + 1] )
            std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
                      << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
    }
#endif

    return moab::MB_SUCCESS;
}
moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags ( const std::string  srcDofTagName,
const std::string  tgtDofTagName 
)

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

Definition at line 181 of file TempestOnlineMap.cpp.

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

{
    moab::ErrorCode rval;

    int tagSize = 0;
    tagSize     = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
    rval =
        m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );

    if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
    {
        int ntot_elements = 0, nelements = m_remapper->m_source_entities.size();
#ifdef MOAB_HAVE_MPI
        int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
        if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
#else
        ntot_elements = nelements;
#endif

        rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_covering_source_entities,
                                                   &m_remapper->m_source_entities, srcDofTagName, m_nDofsPEl_Src );MB_CHK_ERR( rval );

        rval = m_interface->tag_get_handle( srcDofTagName.c_str(), m_nDofsPEl_Src * m_nDofsPEl_Src, MB_TYPE_INTEGER,
                                            this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR( rval );
    }
    else
        MB_CHK_ERR( rval );

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

        rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_target_entities, NULL, tgtDofTagName,
                                                   m_nDofsPEl_Dest );MB_CHK_ERR( rval );

        rval = m_interface->tag_get_handle( tgtDofTagName.c_str(), m_nDofsPEl_Dest * m_nDofsPEl_Dest, MB_TYPE_INTEGER,
                                            this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR( rval );
    }
    else
        MB_CHK_ERR( rval );

    return moab::MB_SUCCESS;
}

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

Definition at line 474 of file TempestOnlineMap.hpp.

{
    m_nDofsPEl_Src = ns;
}

Definition at line 121 of file TempestOnlineMap.cpp.

References moab::TempestRemapper::RLL.

Referenced by TempestOnlineMap().

{
    if( m_meshInputCov )
    {
        std::vector< std::string > dimNames;
        std::vector< int > dimSizes;
        if( m_remapper->m_source_type == moab::TempestRemapper::RLL && m_remapper->m_source_metadata.size() )
        {
            dimNames.push_back( "lat" );
            dimNames.push_back( "lon" );
            dimSizes.resize( 2, 0 );
            dimSizes[0] = m_remapper->m_source_metadata[1];
            dimSizes[1] = m_remapper->m_source_metadata[2];
        }
        else
        {
            dimNames.push_back( "num_elem" );
            dimSizes.push_back( m_meshInputCov->faces.size() );
        }

        this->InitializeSourceDimensions( dimNames, dimSizes );
    }

    if( m_meshOutput )
    {
        std::vector< std::string > dimNames;
        std::vector< int > dimSizes;
        if( m_remapper->m_target_type == moab::TempestRemapper::RLL && m_remapper->m_target_metadata.size() )
        {
            dimNames.push_back( "lat" );
            dimNames.push_back( "lon" );
            dimSizes.resize( 2, 0 );
            dimSizes[0] = m_remapper->m_target_metadata[1];
            dimSizes[1] = m_remapper->m_target_metadata[2];
        }
        else
        {
            dimNames.push_back( "num_elem" );
            dimSizes.push_back( m_meshOutput->faces.size() );
        }

        this->InitializeTargetDimensions( dimNames, dimSizes );
    }
}
moab::ErrorCode moab::TempestOnlineMap::WriteHDF5MapFile ( const std::string &  filename) [private]

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

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

Definition at line 739 of file TempestOnlineMapIO.cpp.

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

{
    moab::ErrorCode rval;

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

    /* Let us compute all relevant data for the current original source mesh on the process */
    DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
    DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
    DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
    if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
    {
        this->InitializeCoordinatesFromMeshFV(
            *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
            ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
            m_remapper->max_source_edges );

        vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
        for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
            vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
    }
    else
    {
        DataArray3D< double > dataGLLJacobianSrc;
        this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
                                               dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );

        // Generate the continuous Jacobian for input mesh
        GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );

        if( m_srcDiscType == DiscretizationType_CGLL )
        {
            GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, m_meshInput->vecFaceArea );
        }
        else
        {
            GenerateDiscontinuousJacobian( dataGLLJacobianSrc, m_meshInput->vecFaceArea );
        }

        vecSourceFaceArea.Allocate( m_meshInput->faces.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
        int offset = 0;
        for( size_t e = 0; e < m_meshInput->faces.size(); e++ )
        {
            for( int s = 0; s < m_nDofsPEl_Src; s++ )
            {
                for( int t = 0; t < m_nDofsPEl_Src; t++ )
                {
                    vecSourceFaceArea[srccol_dtoc_dofmap[offset + s * m_nDofsPEl_Src + t]] =
                        dataGLLJacobianSrc[s][t][e];
                }
            }
            offset += m_nDofsPEl_Src * m_nDofsPEl_Src;
        }
    }

    if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
    {
        this->InitializeCoordinatesFromMeshFV(
            *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
            ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
            m_remapper->max_target_edges );

        vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
        for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
            vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
    }
    else
    {
        DataArray3D< double > dataGLLJacobianDest;
        this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
                                               dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );

        // Generate the continuous Jacobian for input mesh
        GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );

        if( m_destDiscType == DiscretizationType_CGLL )
        {
            GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, m_meshOutput->vecFaceArea );
        }
        else
        {
            GenerateDiscontinuousJacobian( dataGLLJacobianDest, m_meshOutput->vecFaceArea );
        }

        vecTargetFaceArea.Allocate( m_meshOutput->faces.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
        int offset = 0;
        for( size_t e = 0; e < m_meshOutput->faces.size(); e++ )
        {
            for( int s = 0; s < m_nDofsPEl_Dest; s++ )
            {
                for( int t = 0; t < m_nDofsPEl_Dest; t++ )
                {
                    vecTargetFaceArea[row_dtoc_dofmap[offset + s * m_nDofsPEl_Dest + t]] = dataGLLJacobianDest[s][t][e];
                }
            }
            offset += m_nDofsPEl_Dest * m_nDofsPEl_Dest;
        }
    }

    moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
    int tot_src_ents                     = m_remapper->m_source_entities.size();
    int tot_tgt_ents                     = m_remapper->m_target_entities.size();
    int tot_src_size                     = dSourceCenterLon.GetRows();
    int tot_tgt_size                     = m_dTargetCenterLon.GetRows();
    int tot_vsrc_size                    = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
    int tot_vtgt_size                    = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();

    const int weightMatNNZ = m_weightMatrix.nonZeros();
    moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
    rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    moab::Tag srcAreaValues, tgtAreaValues;
    rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
    rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
    rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
                                        moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    moab::Tag srcMaskValues, tgtMaskValues;
    if( m_iSourceMask.IsAttached() )
    {
        rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
                                            moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    }
    if( m_iTargetMask.IsAttached() )
    {
        rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
                                            moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
    }

    std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
    std::vector< double > smatvals( weightMatNNZ );
    // const double* smatvals = m_weightMatrix.valuePtr();
    // Loop over the matrix entries and find the max global ID for rows and columns
    for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
    {
        for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
        {
            smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
            smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
            smatvals[offset]    = it.value();
        }
    }

    /* Set the global IDs for the DoFs */
    ////
    // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
    // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
    ////
    int maxrow = 0, maxcol = 0;
    std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
    for( int i = 0; i < tot_src_size; ++i )
    {
        src_global_dofs[i] = srccol_gdofmap[i];
        maxcol             = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
    }

    for( int i = 0; i < tot_tgt_size; ++i )
    {
        tgt_global_dofs[i] = row_gdofmap[i];
        maxrow             = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
    }

    ///////////////////////////////////////////////////////////////////////////
    // The metadata in H5M file contains the following data:
    //
    //   1. n_a: Total source entities: (number of elements in source mesh)
    //   2. n_b: Total target entities: (number of elements in target mesh)
    //   3. nv_a: Max edge size of elements in source mesh
    //   4. nv_b: Max edge size of elements in target mesh
    //   5. maxrows: Number of rows in remap weight matrix
    //   6. maxcols: Number of cols in remap weight matrix
    //   7. nnz: Number of total nnz in sparse remap weight matrix
    //   8. np_a: The order of the field description on the source mesh: >= 1
    //   9. np_b: The order of the field description on the target mesh: >= 1
    //   10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
    //   dGLL]
    //   11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
    //   dGLL]
    //   12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
    //   1]
    //   13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
    //   [0, 1, 2]
    //
    ///////////////////////////////////////////////////////////////////////////
    int map_disc_details[6];
    map_disc_details[0] = m_nDofsPEl_Src;
    map_disc_details[1] = m_nDofsPEl_Dest;
    map_disc_details[2] = ( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD
                                ? 0
                                : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
    map_disc_details[3] = ( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD
                                ? 0
                                : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
    map_disc_details[4] = ( m_bConserved ? 1 : 0 );
    map_disc_details[5] = m_iMonotonicity;

#ifdef MOAB_HAVE_MPI
    int loc_smatmetadata[13] = { tot_src_ents,
                                 tot_tgt_ents,
                                 m_remapper->max_source_edges,
                                 m_remapper->max_target_edges,
                                 maxrow + 1,
                                 maxcol + 1,
                                 weightMatNNZ,
                                 map_disc_details[0],
                                 map_disc_details[1],
                                 map_disc_details[2],
                                 map_disc_details[3],
                                 map_disc_details[4],
                                 map_disc_details[5] };
    rval                     = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    int glb_smatmetadata[13] = { 0,
                                 0,
                                 0,
                                 0,
                                 0,
                                 0,
                                 0,
                                 map_disc_details[0],
                                 map_disc_details[1],
                                 map_disc_details[2],
                                 map_disc_details[3],
                                 map_disc_details[4],
                                 map_disc_details[5] };
    int loc_buf[7]           = {
        tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
        maxrow,       maxcol
    };
    int glb_buf[4] = { 0, 0, 0, 0 };
    MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
    glb_smatmetadata[0] = glb_buf[0];
    glb_smatmetadata[1] = glb_buf[1];
    glb_smatmetadata[6] = glb_buf[2];
    MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
    glb_smatmetadata[2] = glb_buf[0];
    glb_smatmetadata[3] = glb_buf[1];
    glb_smatmetadata[4] = glb_buf[2];
    glb_smatmetadata[5] = glb_buf[3];
#else
    int glb_smatmetadata[13] = { tot_src_ents,
                                 tot_tgt_ents,
                                 m_remapper->max_source_edges,
                                 m_remapper->max_target_edges,
                                 maxrow,
                                 maxcol,
                                 weightMatNNZ,
                                 map_disc_details[0],
                                 map_disc_details[1],
                                 map_disc_details[2],
                                 map_disc_details[3],
                                 map_disc_details[4],
                                 map_disc_details[5] };
#endif
    // These values represent number of rows and columns. So should be 1-based.
    glb_smatmetadata[4]++;
    glb_smatmetadata[5]++;

    if( this->is_root )
    {
        std::cout << "  " << this->rank << "  Writing remap weights with size [" << glb_smatmetadata[4] << " X "
                  << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
        EntityHandle root_set = 0;
        rval                  = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    }

    int dsize;
    const int numval          = weightMatNNZ;
    const void* smatrowvals_d = smatrowvals.data();
    const void* smatcolvals_d = smatcolvals.data();
    const void* smatvals_d    = smatvals.data();
    rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

    /* Set the global IDs for the DoFs */
    const void* srceleidvals_d = src_global_dofs.data();
    const void* tgteleidvals_d = tgt_global_dofs.data();
    dsize                      = src_global_dofs.size();
    rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    dsize = tgt_global_dofs.size();
    rval  = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

    /* Set the source and target areas */
    const void* srcareavals_d = vecSourceFaceArea;
    const void* tgtareavals_d = vecTargetFaceArea;
    dsize                     = tot_src_size;
    rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    dsize = tot_tgt_size;
    rval  = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

    /* Set the coordinates for source and target center vertices */
    const void* srccoordsclonvals_d = &dSourceCenterLon[0];
    const void* srccoordsclatvals_d = &dSourceCenterLat[0];
    dsize                           = dSourceCenterLon.GetRows();
    rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
    const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
    dsize                           = vecTargetFaceArea.GetRows();
    rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

    /* Set the coordinates for source and target element vertices */
    const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
    const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
    dsize                           = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
    rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
    const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
    dsize                           = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
    rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );

    /* Set the masks for source and target meshes if available */
    if( m_iSourceMask.IsAttached() )
    {
        const void* srcmaskvals_d = m_iSourceMask;
        dsize                     = m_iSourceMask.GetRows();
        rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    }

    if( m_iTargetMask.IsAttached() )
    {
        const void* tgtmaskvals_d = m_iTargetMask;
        dsize                     = m_iTargetMask.GetRows();
        rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
    }

#ifdef MOAB_HAVE_MPI
    const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
#else
    const char* writeOptions = "";
#endif

    // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
    EntityHandle sets[1] = { m_remapper->m_overlap_set };
    rval                 = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );

#ifdef WRITE_SCRIP_FILE
    sstr.str( "" );
    sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
    std::map< std::string, std::string > mapAttributes;
    mapAttributes["Creator"] = "MOAB mbtempest workflow";
    if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
    this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
    sstr.str( "" );
#endif

    return moab::MB_SUCCESS;
}
moab::ErrorCode moab::TempestOnlineMap::WriteParallelMap ( const std::string &  strTarget)

Write the TempestOnlineMap to a parallel NetCDF file.

Definition at line 146 of file TempestOnlineMapIO.cpp.

References ErrorCode, and MB_CHK_ERR.

Referenced by main().

{
    moab::ErrorCode rval;

    size_t lastindex      = strFilename.find_last_of( "." );
    std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );

    // Write the map file to disk in parallel
    if( extension == "nc" )
    {
        /* Invoke the actual call to write the parallel map to disk in SCRIP format */
        rval = this->WriteSCRIPMapFile( strFilename.c_str() );MB_CHK_ERR( rval );
    }
    else
    {
        /* Write to the parallel H5M format */
        rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
    }

    return rval;
}
moab::ErrorCode moab::TempestOnlineMap::WriteSCRIPMapFile ( const std::string &  strOutputFile) [private]

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

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

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

Definition at line 170 of file TempestOnlineMapIO.cpp.

References moab::TupleList::enableWriteAccess(), moab::error(), moab::TupleList::get_n(), ierr, moab::TupleList::inc_n(), moab::TupleList::initialize(), MB_SUCCESS, nr, rank, moab::TempestRemapper::RLL, size, moab::TupleList::vi_wr, and moab::TupleList::vr_wr.

{
    NcError error( NcError::silent_nonfatal );

#ifdef MOAB_HAVE_NETCDFPAR
    bool is_independent = true;
    ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
    // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
#else
    NcFile ncMap( strFilename.c_str(), NcFile::Replace );
#endif

    if( !ncMap.is_valid() )
    {
        _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
    }

    // Attributes
    ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );

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

    /* Let us compute all relevant data for the current original source mesh on the process */
    DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
    DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
    DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
    if( m_srcDiscType == DiscretizationType_FV || m_srcDiscType == DiscretizationType_PCLOUD )
    {
        this->InitializeCoordinatesFromMeshFV(
            *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
            ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
            m_remapper->max_source_edges );

        vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
        for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
            vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
    }
    else
    {
        DataArray3D< double > dataGLLJacobianSrc;
        this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
                                               dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );

        // Generate the continuous Jacobian for input mesh
        GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );

        if( m_srcDiscType == DiscretizationType_CGLL )
        {
            GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
        }
        else
        {
            GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
        }
    }

    if( m_destDiscType == DiscretizationType_FV || m_destDiscType == DiscretizationType_PCLOUD )
    {
        this->InitializeCoordinatesFromMeshFV(
            *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
            ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
            m_remapper->max_target_edges );

        vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
        for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
        {
            vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
        }
    }
    else
    {
        DataArray3D< double > dataGLLJacobianDest;
        this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
                                               dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );

        // Generate the continuous Jacobian for input mesh
        GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );

        if( m_destDiscType == DiscretizationType_CGLL )
        {
            GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
        }
        else
        {
            GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
        }
    }

    // Map dimensions
    unsigned nA = ( vecSourceFaceArea.GetRows() );
    unsigned nB = ( vecTargetFaceArea.GetRows() );

    // 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( int 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, 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, nB, nTargetNodesPerFace,
                                         max_row_dof );  // now nA will be close to maxdof/size
        if( ierr != 0 )
        {
            _EXCEPTION1( "Unable to arrange target data %d ", nB );
        }
    }
#endif

    // Number of non-zeros in the remap matrix operator
    int nS = m_weightMatrix.nonZeros();

#if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
    int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
    int offbuf[3] = { 0, 0, 0 };
    int globuf[5] = { 0, 0, 0, 0, 0 };
    MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
    MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
    MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );

    // MPI_Scan is inclusive of data in current rank; modify accordingly.
    offbuf[0] -= nA;
    offbuf[1] -= nB;
    offbuf[2] -= nS;

#else
    int offbuf[3]            = { 0, 0, 0 };
    int globuf[5]            = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
#endif

    // Write output dimensions entries
    unsigned nSrcGridDims = ( m_vecSourceDimSizes.size() );
    unsigned nDstGridDims = ( m_vecTargetDimSizes.size() );

    NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
    NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );

    NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
    NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );

#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varSrcGridDims, is_independent );
    ncMap.enable_var_par_access( varDstGridDims, is_independent );
#endif

    char szDim[64];
    if( ( nSrcGridDims == 1 ) && ( m_vecSourceDimSizes[0] != (int)nA ) )
    {
        varSrcGridDims->put( &globuf[0], 1 );
        varSrcGridDims->add_att( "name0", "num_dof" );
    }
    else
    {
        for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
        {
            int tmp = ( i == 0 ? globuf[0] : m_vecSourceDimSizes[i] );
            varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
            varSrcGridDims->put( &( tmp ), 1 );
        }

        for( unsigned i = 0; i < m_vecSourceDimSizes.size(); i++ )
        {
            sprintf( szDim, "name%i", i );
            varSrcGridDims->add_att( szDim, m_vecSourceDimNames[nSrcGridDims - i - 1].c_str() );
        }
    }

    if( ( nDstGridDims == 1 ) && ( m_vecTargetDimSizes[0] != (int)nB ) )
    {
        varDstGridDims->put( &globuf[1], 1 );
        varDstGridDims->add_att( "name0", "num_dof" );
    }
    else
    {
        for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
        {
            int tmp = ( i == 0 ? globuf[1] : m_vecTargetDimSizes[i] );
            varDstGridDims->set_cur( nDstGridDims - i - 1 );
            varDstGridDims->put( &( tmp ), 1 );
        }

        for( unsigned i = 0; i < m_vecTargetDimSizes.size(); i++ )
        {
            sprintf( szDim, "name%i", i );
            varDstGridDims->add_att( szDim, m_vecTargetDimNames[nDstGridDims - i - 1].c_str() );
        }
    }

    // Source and Target mesh resolutions
    NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
    NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );

    // Number of nodes per Face
    NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
    NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );

    // Write coordinates
    NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
    NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );

    NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
    NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );

    NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
    NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );

    NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
    NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );

#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 );
#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 );

    // Write areas
    NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varAreaA, is_independent );
#endif
    varAreaA->set_cur( (long)offbuf[0] );
    varAreaA->put( &( vecSourceFaceArea[0] ), nA );

    NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varAreaB, is_independent );
#endif
    varAreaB->set_cur( (long)offbuf[1] );
    varAreaB->put( &( vecTargetFaceArea[0] ), nB );

    // Write SparseMatrix entries
    DataArray1D< int > vecRow( nS );
    DataArray1D< int > vecCol( nS );
    DataArray1D< double > vecS( nS );
    DataArray1D< double > dFracA( nA );
    DataArray1D< double > dFracB( nB );

    moab::TupleList tlValRow, tlValCol;
    unsigned numr = 1;  //
    // value has to be sent to processor row/nB for for fracA and col/nA for fracB
    // vecTargetArea (indexRow ) has to be sent for fracA (index col?)
    // vecTargetFaceArea will have to be sent to col index, with its index !
    tlValRow.initialize( 2, 0, 0, numr, nS );  // to proc(row),  global row , value
    tlValCol.initialize( 3, 0, 0, numr, nS );  // to proc(col),  global row / col, value
    tlValRow.enableWriteAccess();
    tlValCol.enableWriteAccess();
    /*
     *
             dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
             dFracB[ row ] += val ;
     */
    int offset = 0;
#if defined( MOAB_HAVE_MPI )
    int nAbase = ( max_col_dof + 1 ) / size;  // it is nA, except last rank ( == size - 1 )
    int nBbase = ( max_row_dof + 1 ) / size;  // it is nB, except last rank ( == size - 1 )
#endif
    for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
    {
        for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
        {
            vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() );  // row index
            vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() );  // col index
            vecS[offset]   = it.value();                             // value

#if defined( MOAB_HAVE_MPI )
            {
                // value M(row, col) will contribute to procRow and procCol values for fracA and fracB
                int procRow = ( vecRow[offset] - 1 ) / nBbase;
                if( procRow >= size ) procRow = size - 1;
                int procCol = ( vecCol[offset] - 1 ) / nAbase;
                if( procCol >= size ) procCol = size - 1;
                int nrInd                     = tlValRow.get_n();
                tlValRow.vi_wr[2 * nrInd]     = procRow;
                tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
                tlValRow.vr_wr[nrInd]         = vecS[offset];
                tlValRow.inc_n();
                int ncInd                     = tlValCol.get_n();
                tlValCol.vi_wr[3 * ncInd]     = procCol;
                tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
                tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1;  // this is column
                tlValCol.vr_wr[ncInd]         = vecS[offset];
                tlValCol.inc_n();
            }

#endif
            offset++;
        }
    }
#if defined( MOAB_HAVE_MPI )
    // need to send values for their row and col processors, to compute fractions there
    // now do the heavy communication
    ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
    ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );

    // we have now, for example,  dFracB[ row ] += val ;
    // so we know that on current task, we received tlValRow
    // reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
    //          dFracB[ row ] += val ;
    for( unsigned i = 0; i < tlValRow.get_n(); i++ )
    {
        // int fromProc = tlValRow.vi_wr[2 * i];
        int gRowInd       = tlValRow.vi_wr[2 * i + 1];
        int localIndexRow = gRowInd - nBbase * rank;  // modulo nBbase rank is from 0 to size - 1;
        double wgt        = tlValRow.vr_wr[i];
        assert( localIndexRow >= 0 );
        assert( nB - localIndexRow > 0 );
        dFracB[localIndexRow] += wgt;
    }
    // to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from

    std::set< int > neededRows;
    for( unsigned i = 0; i < tlValCol.get_n(); i++ )
    {
        int rRowInd = tlValCol.vi_wr[3 * i + 1];
        neededRows.insert( rRowInd );
        // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
    }
    moab::TupleList tgtAreaReq;
    tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
    tgtAreaReq.enableWriteAccess();
    for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
    {
        int neededRow = *sit;
        int procRow   = neededRow / nBbase;
        if( procRow >= size ) procRow = size - 1;
        int nr                       = tgtAreaReq.get_n();
        tgtAreaReq.vi_wr[2 * nr]     = procRow;
        tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
        tgtAreaReq.inc_n();
    }

    ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
    // we need to send back the tgtArea corresponding to row
    moab::TupleList tgtAreaInfo;  // load it with tgtArea at row
    tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
    tgtAreaInfo.enableWriteAccess();
    for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
    {
        int from_proc     = tgtAreaReq.vi_wr[2 * i];
        int row           = tgtAreaReq.vi_wr[2 * i + 1];
        int locaIndexRow  = row - rank * nBbase;
        double areaToSend = vecTargetFaceArea[locaIndexRow];
        // int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;

        tgtAreaInfo.vi_wr[2 * i]     = from_proc;  // send back requested info
        tgtAreaInfo.vi_wr[2 * i + 1] = row;
        tgtAreaInfo.vr_wr[i]         = areaToSend;  // this will be tgt area at row
        tgtAreaInfo.inc_n();
    }
    ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );

    std::map< int, double > areaAtRow;
    for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
    {
        // we have received from proc, value for row !
        int row        = tgtAreaInfo.vi_wr[2 * i + 1];
        areaAtRow[row] = tgtAreaInfo.vr_wr[i];
    }

    // we have now for rows the
    // it is ordered by index, so:
    // now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
    // tgtAreaInfo will have at index i the area we need (from row)
    // there should be an easier way :(
    for( unsigned i = 0; i < tlValCol.get_n(); i++ )
    {
        int rRowInd     = tlValCol.vi_wr[3 * i + 1];
        int colInd      = tlValCol.vi_wr[3 * i + 2];
        double val      = tlValCol.vr_wr[i];
        int localColInd = colInd - rank * nAbase;  // < local nA
        // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
        auto itMap = areaAtRow.find( rRowInd );  // it should be different from end
        if( itMap != areaAtRow.end() )
        {
            double areaRow = itMap->second;  // we fished a lot for this !
            dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
        }
    }

#endif
    // Load in data
    NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );

    NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
    NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
    NcVar* varS   = ncMap.add_var( "S", ncDouble, dimNS );
#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varRow, is_independent );
    ncMap.enable_var_par_access( varCol, is_independent );
    ncMap.enable_var_par_access( varS, is_independent );
#endif

    varRow->set_cur( (long)offbuf[2] );
    varRow->put( vecRow, nS );

    varCol->set_cur( (long)offbuf[2] );
    varCol->put( vecCol, nS );

    varS->set_cur( (long)offbuf[2] );
    varS->put( &( vecS[0] ), nS );

    // Calculate and write fractional coverage arrays
    NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varFracA, is_independent );
#endif
    varFracA->add_att( "name", "fraction of target coverage of source dof" );
    varFracA->add_att( "units", "unitless" );
    varFracA->set_cur( (long)offbuf[0] );
    varFracA->put( &( dFracA[0] ), nA );

    NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
#ifdef MOAB_HAVE_NETCDFPAR
    ncMap.enable_var_par_access( varFracB, is_independent );
#endif
    varFracB->add_att( "name", "fraction of source coverage of target dof" );
    varFracB->add_att( "units", "unitless" );
    varFracB->set_cur( (long)offbuf[1] );
    varFracB->put( &( dFracB[0] ), nB );

    // Add global attributes
    // std::map<std::string, std::string>::const_iterator iterAttributes =
    //     mapAttributes.begin();
    // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
    //     ncMap.add_att(
    //         iterAttributes->first.c_str(),
    //         iterAttributes->second.c_str());
    // }

    ncMap.close();

    return moab::MB_SUCCESS;
}

Member Data Documentation

std::vector< int > moab::TempestOnlineMap::col_dtoc_dofmap [private]

Definition at line 414 of file TempestOnlineMap.hpp.

std::vector< unsigned > moab::TempestOnlineMap::col_gdofmap [private]

Definition at line 411 of file TempestOnlineMap.hpp.

Referenced by fill_col_ids(), and LinearRemapNN_MOAB().

std::map< int, int > moab::TempestOnlineMap::colMap [private]

Definition at line 416 of file TempestOnlineMap.hpp.

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesDest [private]

Definition at line 418 of file TempestOnlineMap.hpp.

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrc [private]

Definition at line 418 of file TempestOnlineMap.hpp.

DataArray3D< int > moab::TempestOnlineMap::dataGLLNodesSrcCov [private]

Definition at line 418 of file TempestOnlineMap.hpp.

Definition at line 433 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 433 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 425 of file TempestOnlineMap.hpp.

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

Definition at line 410 of file TempestOnlineMap.hpp.

Definition at line 426 of file TempestOnlineMap.hpp.

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

Definition at line 398 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 428 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 429 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 430 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 431 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 423 of file TempestOnlineMap.hpp.

Definition at line 423 of file TempestOnlineMap.hpp.

Definition at line 420 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

Definition at line 420 of file TempestOnlineMap.hpp.

Definition at line 420 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

The fundamental remapping operator object.

Definition at line 385 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 434 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

Definition at line 434 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

std::vector< int > moab::TempestOnlineMap::row_dtoc_dofmap [private]

Definition at line 414 of file TempestOnlineMap.hpp.

std::vector< unsigned > moab::TempestOnlineMap::row_gdofmap [private]

Definition at line 411 of file TempestOnlineMap.hpp.

Referenced by fill_row_ids(), and LinearRemapNN_MOAB().

std::map< int, int > moab::TempestOnlineMap::rowMap [private]

Definition at line 416 of file TempestOnlineMap.hpp.

Definition at line 434 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

std::vector< int > moab::TempestOnlineMap::srccol_dtoc_dofmap [private]

Definition at line 414 of file TempestOnlineMap.hpp.

std::vector< unsigned > moab::TempestOnlineMap::srccol_gdofmap [private]

Definition at line 411 of file TempestOnlineMap.hpp.

List of all members.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines