Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
moab::TempestOnlineMap Class Reference

An offline map between two Meshes. More...

#include <TempestOnlineMap.hpp>

+ Inheritance diagram for moab::TempestOnlineMap:
+ Collaboration diagram for moab::TempestOnlineMap:

Public Types

enum  DiscretizationType { DiscretizationType_FV = 0 , DiscretizationType_CGLL = 1 , DiscretizationType_DGLL = 2 , DiscretizationType_PCLOUD = 3 }
 
enum  CAASType {
  CAAS_NONE = 0 , CAAS_GLOBAL = 1 , CAAS_LOCAL = 2 , CAAS_LOCAL_ADJACENT = 3 ,
  CAAS_QLT = 4
}
 
typedef double(* sample_function) (double, double)
 

Public Member Functions

 TempestOnlineMap (moab::TempestRemapper *remapper)
 Generate the metadata associated with the offline map. More...
 
virtual ~TempestOnlineMap ()
 Define a virtual destructor. More...
 
moab::ErrorCode GenerateRemappingWeights (std::string strInputType, std::string strOutputType, const GenerateOfflineMapAlgorithmOptions &mapOptions, const std::string &srcDofTagName="GLOBAL_ID", const std::string &tgtDofTagName="GLOBAL_ID")
 Generate the offline map, given the source and target mesh and discretization details. This method generates the mapping between the two meshes based on the overlap and stores the result in the SparseMatrix. More...
 
moab::ErrorCode ReadParallelMap (const char *strSource, const std::vector< int > &owned_dof_ids, bool row_major_ownership=true)
 Generate the metadata associated with the offline map. More...
 
moab::ErrorCode WriteParallelMap (const std::string &strTarget, const std::map< std::string, std::string > &attrMap)
 Write the TempestOnlineMap to a parallel NetCDF file. More...
 
virtual int IsConsistent (double dTolerance)
 Determine if the map is first-order accurate. More...
 
virtual int IsConservative (double dTolerance)
 Determine if the map is conservative. More...
 
virtual int IsMonotone (double dTolerance)
 Determine if the map is monotone. More...
 
const DataArray1D< double > & GetGlobalSourceAreas () const
 If we computed the reduction, get the vector representing the source areas for all entities in the mesh More...
 
const DataArray1D< double > & GetGlobalTargetAreas () const
 If we computed the reduction, get the vector representing the target areas for all entities in the mesh More...
 
void PrintMapStatistics ()
 
moab::ErrorCode SetDOFmapTags (const std::string srcDofTagName, const std::string tgtDofTagName)
 Store the tag names associated with global DoF ids for source and target meshes to be used for mapping.

Parameters
srcDofTagNameThe tag name associated with global DoF ids for the source mesh
tgtDofTagNameThe tag name associated with global DoF ids for the target mesh
More...
 
moab::ErrorCode SetDOFmapAssociation (DiscretizationType srcType, int srcOrder, bool isSrcContinuous, DataArray3D< int > *srcdataGLLNodes, DataArray3D< int > *srcdataGLLNodesSrc, DiscretizationType destType, int destOrder, bool isTgtContinuous, DataArray3D< int > *tgtdataGLLNodes)
 
std::pair< double, double > ApplyBoundsLimiting (std::vector< double > &dataInDouble, std::vector< double > &dataOutDouble, CAASType caasType=CAAS_GLOBAL, int caasIteration=0, double mismatch=0.0)
 
void ComputeAdjacencyRelations (std::vector< std::unordered_set< int > > &vecAdjFaces, int nrings, const Range &entities, bool useMOABAdjacencies=true, Mesh *trMesh=nullptr)
 
int GetSourceGlobalNDofs ()
 Get the number of total Degrees-Of-Freedom defined on the source mesh. More...
 
int GetDestinationGlobalNDofs ()
 Get the number of total Degrees-Of-Freedom defined on the destination mesh. More...
 
int GetSourceLocalNDofs ()
 Get the number of local Degrees-Of-Freedom defined on the source mesh. More...
 
int GetDestinationLocalNDofs ()
 Get the number of local Degrees-Of-Freedom defined on the destination mesh. More...
 
int GetSourceNDofsPerElement ()
 Get the number of Degrees-Of-Freedom per element on the source mesh. More...
 
int GetDestinationNDofsPerElement ()
 Get the number of Degrees-Of-Freedom per element on the destination mesh. More...
 
void SetSourceNDofsPerElement (int ns)
 Set the number of Degrees-Of-Freedom per element on the source mesh. More...
 
void SetDestinationNDofsPerElement (int nt)
 Get the number of Degrees-Of-Freedom per element on the destination mesh. More...
 
int GetRowGlobalDoF (int localID) const
 Get the global Degrees-Of-Freedom ID on the destination mesh. More...
 
int GetIndexOfRowGlobalDoF (int globalRowDoF) const
 Get the index of globaRowDoF. More...
 
int GetColGlobalDoF (int localID) const
 Get the global Degrees-Of-Freedom ID on the source mesh. More...
 
int GetIndexOfColGlobalDoF (int globalColDoF) const
 Get the index of globaColDoF. More...
 
moab::ErrorCode ApplyWeights (moab::Tag srcSolutionTag, moab::Tag tgtSolutionTag, bool transpose=false, CAASType caasType=CAAS_NONE)
 Apply the weight matrix onto the source vector (tag) provided as input, and return the column vector (solution projection) in a tag, after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More...
 
moab::ErrorCode DefineAnalyticalSolution (moab::Tag &exactSolnTag, const std::string &solnName, Remapper::IntersectionContext ctx, sample_function testFunction, moab::Tag *clonedSolnTag=NULL, std::string cloneSolnName="")
 Define an analytical solution over the given (source or target) mesh, as specificed in the context. This routine will define a tag that is compatible with the specified discretization method type and order and sample the solution exactly using the analytical function provided by the user. More...
 
moab::ErrorCode ComputeMetrics (Remapper::IntersectionContext ctx, moab::Tag &exactTag, moab::Tag &approxTag, std::map< std::string, double > &metrics, bool verbose=true)
 Compute the error between a sampled (exact) solution and a projected solution in various error norms. More...
 
moab::ErrorCode fill_row_ids (std::vector< int > &ids_of_interest)
 
moab::ErrorCode fill_col_ids (std::vector< int > &ids_of_interest)
 
moab::ErrorCode set_col_dc_dofs (std::vector< int > &values_entities)
 
moab::ErrorCode set_row_dc_dofs (std::vector< int > &values_entities)
 
void SetMeshInput (Mesh *imesh)
 

Private Member Functions

moab::ErrorCode LinearRemapNN_MOAB (bool use_GID_matching=false, bool strict_check=false)
 Compute the remapping weights as a permutation matrix that relates DoFs on the source mesh to DoFs on the target mesh. More...
 
void LinearRemapFVtoFV_Tempest_MOAB (int nOrder)
 Compute the remapping weights for a FV field defined on the source to a FV field defined on the target mesh. More...
 
void LinearRemapSE0_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian)
 Generate the OfflineMap for linear conserative element-average spectral element to element average remapping. More...
 
void LinearRemapSE4_Tempest_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, int nMonotoneType, bool fContinuousIn, bool fNoConservation)
 Generate the OfflineMap for cubic conserative element-average spectral element to element average remapping. More...
 
void LinearRemapFVtoGLL_MOAB (const DataArray3D< int > &dataGLLNodes, const DataArray3D< double > &dataGLLJacobian, const DataArray1D< double > &dataGLLNodalArea, int nOrder, int nMonotoneType, bool fContinuous, bool fNoConservation)
 Generate the OfflineMap for remapping from finite volumes to finite elements. More...
 
void LinearRemapGLLtoGLL2_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut, bool fNoConservation)
 Generate the OfflineMap for remapping from finite elements to finite elements. More...
 
void LinearRemapGLLtoGLL2_Pointwise_MOAB (const DataArray3D< int > &dataGLLNodesIn, const DataArray3D< double > &dataGLLJacobianIn, const DataArray3D< int > &dataGLLNodesOut, const DataArray3D< double > &dataGLLJacobianOut, const DataArray1D< double > &dataNodalAreaOut, int nPin, int nPout, int nMonotoneType, bool fContinuousIn, bool fContinuousOut)
 Generate the OfflineMap for remapping from finite elements to finite elements (pointwise interpolation). More...
 
moab::ErrorCode WriteSCRIPMapFile (const std::string &strOutputFile, const std::map< std::string, std::string > &attrMap)
 Copy the local matrix from Tempest SparseMatrix representation (ELL) to the parallel CSR Eigen Matrix for scalable application of matvec needed for projections. More...
 
moab::ErrorCode WriteHDF5MapFile (const std::string &filename)
 Parallel I/O with NetCDF to write out the SCRIP file from multiple processors. More...
 
void setup_sizes_dimensions ()
 
void CAASLimiter (std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, double &dMass)
 
double QLTLimiter (int caasIteration, std::vector< double > &dataCorrectedField, std::vector< double > &dataLowerBound, std::vector< double > &dataUpperBound, std::vector< double > &dMassDefect)
 
moab::ErrorCode ApplyWeights (std::vector< double > &srcVals, std::vector< double > &tgtVals, bool transpose=false)
 Apply the weight matrix onto the source vector provided as input, and return the column vector (solution projection) after the map application Compute: tgtVals = A(S->T) * \srcVals, or if (transpose) tgtVals = [A(T->S)]^T * \srcVals More...
 

Private Attributes

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

Detailed Description

An offline map between two Meshes.

Definition at line 52 of file TempestOnlineMap.hpp.

Member Typedef Documentation

◆ sample_function

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

Definition at line 398 of file TempestOnlineMap.hpp.

Member Enumeration Documentation

◆ CAASType

Enumerator
CAAS_NONE 
CAAS_GLOBAL 
CAAS_LOCAL 
CAAS_LOCAL_ADJACENT 
CAAS_QLT 

Definition at line 77 of file TempestOnlineMap.hpp.

78  {
79  CAAS_NONE = 0,
80  CAAS_GLOBAL = 1,
81  CAAS_LOCAL = 2,
83  CAAS_QLT = 4
84  };

◆ DiscretizationType

Enumerator
DiscretizationType_FV 
DiscretizationType_CGLL 
DiscretizationType_DGLL 
DiscretizationType_PCLOUD 

Definition at line 68 of file TempestOnlineMap.hpp.

Constructor & Destructor Documentation

◆ TempestOnlineMap()

moab::TempestOnlineMap::TempestOnlineMap ( moab::TempestRemapper remapper)

Generate the metadata associated with the offline map.

Definition at line 64 of file TempestOnlineMap.cpp.

64  : OfflineMap(), m_remapper( remapper )
65 {
66  // Get the references for the MOAB core objects
68 #ifdef MOAB_HAVE_MPI
69  m_pcomm = m_remapper->get_parallel_communicator();
70 #endif
71 
72  // Update the references to the meshes
74  m_meshInputCov = remapper->GetCoveringMesh();
77 
78  is_parallel = remapper->is_parallel;
79  is_root = remapper->is_root;
80  rank = remapper->rank;
81  size = remapper->size;
82 
83  // set default order
85 
86  // Initialize dimension information from file
87  this->setup_sizes_dimensions();
88 }

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

◆ ~TempestOnlineMap()

moab::TempestOnlineMap::~TempestOnlineMap ( )
virtual

Define a virtual destructor.

Definition at line 115 of file TempestOnlineMap.cpp.

116 {
117  m_interface = nullptr;
118 #ifdef MOAB_HAVE_MPI
119  m_pcomm = nullptr;
120 #endif
121  m_meshInput = nullptr;
122  m_meshOutput = nullptr;
123  m_meshOverlap = nullptr;
124 }

Member Function Documentation

◆ ApplyBoundsLimiting()

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

ApplyBoundsLimiting - Apply bounds limiting to the data field

Parameters
dataInDouble- input data field
dataOutDouble- output data field
caasType- type of limiter
caasIteration- iteration number of limiter
Returns
- pair of mass defect pre and post limiter application

◆ ApplyWeights() [1/2]

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

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

Definition at line 1761 of file TempestOnlineMap.cpp.

1765 {
1766  moab::ErrorCode rval;
1767 
1768  std::vector< double > solSTagVals;
1769  std::vector< double > solTTagVals;
1770 
1771  moab::Range sents, tents;
1773  {
1775  {
1777  solSTagVals.resize( covSrcEnts.size(), -1.0 );
1778  sents = covSrcEnts;
1779  }
1780  else
1781  {
1783  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1784  -1.0 );
1785  sents = covSrcEnts;
1786  }
1788  {
1790  solTTagVals.resize( tgtEnts.size(), -1.0 );
1791  tents = tgtEnts;
1792  }
1793  else
1794  {
1796  solTTagVals.resize(
1797  tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
1798  tents = tgtEnts;
1799  }
1800  }
1801  else
1802  {
1805  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1806  -1.0 );
1807  solTTagVals.resize(
1808  tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
1809 
1810  sents = covSrcEnts;
1811  tents = tgtEnts;
1812  }
1813 
1814  // The tag data is np*np*n_el_src
1815  rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );
1816 
1817  // Compute the application of weights on the suorce solution data and store it in the
1818  // destination solution vector data Optionally, can also perform the transpose application of
1819  // the weight matrix. Set the 3rd argument to true if this is needed
1820  rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );
1821 
1822  if( caasType != CAAS_NONE )
1823  {
1824  std::string tgtSolutionTagName;
1825  rval = m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName );MB_CHK_SET_ERR( rval, "Getting tag name failed" );
1826 
1827  // Perform CAAS iterations iteratively until convergence
1828  constexpr int nmax_caas_iterations = 10;
1829  double mismatch = 1.0;
1830  int caasIteration = 0;
1831  double initialMismatch = 0.0;
1832  while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
1833  caasIteration++ < nmax_caas_iterations ) // iterate until convergence or a maximum of 5 iterations
1834  {
1835  // The tag data is np*np*n_el_dest
1836  rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1837 
1838 // #ifdef MOAB_HAVE_MPI
1839 // rval = m_pcomm->exchange_tags( tgtSolutionTag, tents );MB_CHK_SET_ERR( rval, "Tag exchange failed" );
1840 // #endif
1841 
1842  double dMassDiffPostGlobal;
1843  std::pair< double, double > mDefect =
1844  this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
1845 #ifdef MOAB_HAVE_MPI
1846  double dMassDiffPost = mDefect.second;
1847  MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
1848 #else
1849  dMassDiffPostGlobal = mDefect.second;
1850 #endif
1851  if( caasIteration == 1 ) initialMismatch = mDefect.first;
1852  if( m_remapper->verbose && is_root )
1853  {
1854  printf( "Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
1855  tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
1856  }
1857  mismatch = dMassDiffPostGlobal;
1858  }
1859  }
1860 
1861  // The tag data is np*np*n_el_dest
1862  rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1863 
1864  return moab::MB_SUCCESS;
1865 }

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

Referenced by main().

◆ ApplyWeights() [2/2]

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

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

◆ CAASLimiter()

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

◆ ComputeAdjacencyRelations()

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

Vector storing adjacent Faces.

Definition at line 1712 of file TempestOnlineMap.cpp.

1717 {
1718  assert( nrings > 0 );
1719  assert( useMOABAdjacencies || trMesh != nullptr );
1720 
1721  const size_t nrows = vecAdjFaces.size();
1723  for( size_t index = 0; index < nrows; index++ )
1724  {
1725  vecAdjFaces[index].insert( index ); // add self target face first
1726  {
1727  // Compute the adjacent faces to the target face
1728  if( useMOABAdjacencies )
1729  {
1730  moab::Range ents;
1731  // ents.insert( entities.index( entities[index] ) );
1732  ents.insert( entities[index] );
1733  moab::Range adjEnts;
1734  moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, nrings );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
1735  for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
1736  {
1737  // int adjIndex = m_interface->id_from_handle(*it)-1;
1738  int adjIndex = entities.index( *it );
1739  // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
1740  if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
1741  }
1742  }
1743  else
1744  {
1745  /// Vector storing adjacent Faces.
1746  typedef std::pair< int, int > FaceDistancePair;
1747  typedef std::vector< FaceDistancePair > AdjacentFaceVector;
1748  AdjacentFaceVector adjFaces;
1749  Face& face = trMesh->faces[index];
1750  GetAdjacentFaceVectorByEdge( *trMesh, index, nrings * face.edges.size(), adjFaces );
1751 
1752  // Add the adjacent faces to the target face list
1753  for( auto adjFace : adjFaces )
1754  if( adjFace.first >= 0 )
1755  vecAdjFaces[index].insert( adjFace.first ); // map target face to source face
1756  }
1757  }
1758  }
1759 }

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

◆ ComputeMetrics()

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

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

Definition at line 2229 of file TempestOnlineMap.cpp.

2234 {
2235  moab::ErrorCode rval;
2236  const bool outputEnabled = ( is_root );
2237  int discOrder;
2238  // DiscretizationType discMethod;
2239  // moab::EntityHandle meshset;
2241  // Mesh* trmesh;
2242  switch( ctx )
2243  {
2244  case Remapper::SourceMesh:
2245  // meshset = m_remapper->m_covering_source_set;
2246  // trmesh = m_remapper->m_covering_source;
2249  discOrder = m_nDofsPEl_Src;
2250  // discMethod = m_eInputType;
2251  break;
2252 
2253  case Remapper::TargetMesh:
2254  // meshset = m_remapper->m_target_set;
2255  // trmesh = m_remapper->m_target;
2256  entities =
2258  discOrder = m_nDofsPEl_Dest;
2259  // discMethod = m_eOutputType;
2260  break;
2261 
2262  default:
2263  if( outputEnabled )
2264  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
2265  return moab::MB_FAILURE;
2266  }
2267 
2268  // Let us create teh solution tag with appropriate information for name, discretization order
2269  // (DoF space)
2270  std::string exactTagName, projTagName;
2271  const int ntotsize = entities.size() * discOrder * discOrder;
2272  std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
2273  rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
2274  rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
2275  rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
2276  rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
2277 
2278  const auto& ovents = m_remapper->m_overlap_entities;
2279 
2280  std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 ); // L1Err, L2Err, LinfErr
2281  double sumarea = 0.0;
2282  for( size_t i = 0; i < ovents.size(); ++i )
2283  {
2284  const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
2285  if( srcidx < 0 ) continue; // Skip non-overlapping entities
2286  const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
2287  if( tgtidx < 0 ) continue; // skip ghost target faces
2288  const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
2289  const double error = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
2290  errnorms[0] += ovarea * error;
2291  errnorms[1] += ovarea * error * error;
2292  errnorms[3] = ( error > errnorms[3] ? error : errnorms[3] );
2293  sumarea += ovarea;
2294  }
2295  errnorms[2] = sumarea;
2296 #ifdef MOAB_HAVE_MPI
2297  if( m_pcomm )
2298  {
2299  MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
2300  MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
2301  }
2302 #else
2303  for( int i = 0; i < 4; ++i )
2304  globerrnorms[i] = errnorms[i];
2305 #endif
2306 
2307  globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
2308  globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
2309 
2310  metrics.clear();
2311  metrics["L1Error"] = globerrnorms[0];
2312  metrics["L2Error"] = globerrnorms[1];
2313  metrics["LinfError"] = globerrnorms[3];
2314 
2315  if( verbose && is_root )
2316  {
2317  std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
2318  std::cout << "\t Total Intersection area = " << globerrnorms[2] << std::endl;
2319  std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
2320  std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
2321  std::cout << "\t L_inf error = " << globerrnorms[3] << std::endl;
2322  }
2323 
2324  return moab::MB_SUCCESS;
2325 }

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

Referenced by main().

◆ DefineAnalyticalSolution()

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

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

Definition at line 1867 of file TempestOnlineMap.cpp.

1873 {
1874  moab::ErrorCode rval;
1875  const bool outputEnabled = ( is_root );
1876  int discOrder;
1877  DiscretizationType discMethod;
1878  // moab::EntityHandle meshset;
1880  Mesh* trmesh;
1881  switch( ctx )
1882  {
1883  case Remapper::SourceMesh:
1884  // meshset = m_remapper->m_covering_source_set;
1885  trmesh = m_remapper->m_covering_source;
1888  discOrder = m_nDofsPEl_Src;
1889  discMethod = m_eInputType;
1890  break;
1891 
1892  case Remapper::TargetMesh:
1893  // meshset = m_remapper->m_target_set;
1894  trmesh = m_remapper->m_target;
1895  entities =
1897  discOrder = m_nDofsPEl_Dest;
1898  discMethod = m_eOutputType;
1899  break;
1900 
1901  default:
1902  if( outputEnabled )
1903  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
1904  return moab::MB_FAILURE;
1905  }
1906 
1907  // Let us create teh solution tag with appropriate information for name, discretization order
1908  // (DoF space)
1909  rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
1910  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1911  if( clonedSolnTag != nullptr )
1912  {
1913  if( cloneSolnName.size() == 0 )
1914  {
1915  cloneSolnName = solnName + std::string( "Cloned" );
1916  }
1917  rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
1918  *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1919  }
1920 
1921  // Triangular quadrature rule
1922  const int TriQuadratureOrder = 10;
1923 
1924  if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
1925 
1926  TriangularQuadratureRule triquadrule( TriQuadratureOrder );
1927 
1928  const int TriQuadraturePoints = triquadrule.GetPoints();
1929 
1930  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1931  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1932 
1933  // Output data
1934  DataArray1D< double > dVar;
1935  DataArray1D< double > dVarMB; // re-arranged local MOAB vector
1936 
1937  // Nodal geometric area
1938  DataArray1D< double > dNodeArea;
1939 
1940  // Calculate element areas
1941  // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
1942 
1943  if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
1944  {
1945  /* Get the spectral points and sample the functionals accordingly */
1946  const bool fGLL = true;
1947  const bool fGLLIntegrate = false;
1948 
1949  // Generate grid metadata
1950  DataArray3D< int > dataGLLNodes;
1951  DataArray3D< double > dataGLLJacobian;
1952 
1953  GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
1954 
1955  // Number of elements
1956  int nElements = trmesh->faces.size();
1957 
1958  // Verify all elements are quadrilaterals
1959  for( int k = 0; k < nElements; k++ )
1960  {
1961  const Face& face = trmesh->faces[k];
1962 
1963  if( face.edges.size() != 4 )
1964  {
1965  _EXCEPTIONT( "Non-quadrilateral face detected; "
1966  "incompatible with --gll" );
1967  }
1968  }
1969 
1970  // Number of unique nodes
1971  int iMaxNode = 0;
1972  for( int i = 0; i < discOrder; i++ )
1973  {
1974  for( int j = 0; j < discOrder; j++ )
1975  {
1976  for( int k = 0; k < nElements; k++ )
1977  {
1978  if( dataGLLNodes[i][j][k] > iMaxNode )
1979  {
1980  iMaxNode = dataGLLNodes[i][j][k];
1981  }
1982  }
1983  }
1984  }
1985 
1986  // Get Gauss-Lobatto quadrature nodes
1987  DataArray1D< double > dG;
1988  DataArray1D< double > dW;
1989 
1990  GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
1991 
1992  // Get Gauss quadrature nodes
1993  const int nGaussP = 10;
1994 
1995  DataArray1D< double > dGaussG;
1996  DataArray1D< double > dGaussW;
1997 
1998  GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
1999 
2000  // Allocate data
2001  dVar.Allocate( iMaxNode );
2002  dVarMB.Allocate( discOrder * discOrder * nElements );
2003  dNodeArea.Allocate( iMaxNode );
2004 
2005  // Sample data
2006  for( int k = 0; k < nElements; k++ )
2007  {
2008  const Face& face = trmesh->faces[k];
2009 
2010  // Sample data at GLL nodes
2011  if( fGLL )
2012  {
2013  for( int i = 0; i < discOrder; i++ )
2014  {
2015  for( int j = 0; j < discOrder; j++ )
2016  {
2017 
2018  // Apply local map
2019  Node node;
2020  Node dDx1G;
2021  Node dDx2G;
2022 
2023  ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
2024 
2025  // Sample data at this point
2026  double dNodeLon = atan2( node.y, node.x );
2027  if( dNodeLon < 0.0 )
2028  {
2029  dNodeLon += 2.0 * M_PI;
2030  }
2031  double dNodeLat = asin( node.z );
2032 
2033  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
2034 
2035  dVar[dataGLLNodes[j][i][k] - 1] = dSample;
2036  }
2037  }
2038  // High-order Gaussian integration over basis function
2039  }
2040  else
2041  {
2042  DataArray2D< double > dCoeff( discOrder, discOrder );
2043 
2044  for( int p = 0; p < nGaussP; p++ )
2045  {
2046  for( int q = 0; q < nGaussP; q++ )
2047  {
2048 
2049  // Apply local map
2050  Node node;
2051  Node dDx1G;
2052  Node dDx2G;
2053 
2054  ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
2055 
2056  // Cross product gives local Jacobian
2057  Node nodeCross = CrossProduct( dDx1G, dDx2G );
2058 
2059  double dJacobian =
2060  sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
2061 
2062  // Find components of quadrature point in basis
2063  // of the first Face
2064  SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
2065 
2066  // Sample data at this point
2067  double dNodeLon = atan2( node.y, node.x );
2068  if( dNodeLon < 0.0 )
2069  {
2070  dNodeLon += 2.0 * M_PI;
2071  }
2072  double dNodeLat = asin( node.z );
2073 
2074  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
2075 
2076  // Integrate
2077  for( int i = 0; i < discOrder; i++ )
2078  {
2079  for( int j = 0; j < discOrder; j++ )
2080  {
2081 
2082  double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
2083 
2084  dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
2085 
2086  dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
2087  }
2088  }
2089  }
2090  }
2091  }
2092  }
2093 
2094  // Divide by area
2095  if( fGLLIntegrate )
2096  {
2097  for( size_t i = 0; i < dVar.GetRows(); i++ )
2098  {
2099  dVar[i] /= dNodeArea[i];
2100  }
2101  }
2102 
2103  // Let us rearrange the data based on DoF ID specification
2104  if( ctx == Remapper::SourceMesh )
2105  {
2106  for( unsigned j = 0; j < entities.size(); j++ )
2107  for( int p = 0; p < discOrder; p++ )
2108  for( int q = 0; q < discOrder; q++ )
2109  {
2110  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2111  dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
2112  }
2113  }
2114  else
2115  {
2116  for( unsigned j = 0; j < entities.size(); j++ )
2117  for( int p = 0; p < discOrder; p++ )
2118  for( int q = 0; q < discOrder; q++ )
2119  {
2120  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2121  dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
2122  }
2123  }
2124 
2125  // Set the tag data
2126  rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
2127  }
2128  else
2129  {
2130  // assert( discOrder == 1 );
2131  if( discMethod == DiscretizationType_FV )
2132  {
2133  /* Compute an element-wise integral to store the sampled solution based on Quadrature
2134  * rules */
2135  // Resize the array
2136  dVar.Allocate( trmesh->faces.size() );
2137 
2138  std::vector< Node >& nodes = trmesh->nodes;
2139 
2140  // Loop through all Faces
2141  for( size_t i = 0; i < trmesh->faces.size(); i++ )
2142  {
2143  const Face& face = trmesh->faces[i];
2144 
2145  // Loop through all sub-triangles
2146  for( size_t j = 0; j < face.edges.size() - 2; j++ )
2147  {
2148 
2149  const Node& node0 = nodes[face[0]];
2150  const Node& node1 = nodes[face[j + 1]];
2151  const Node& node2 = nodes[face[j + 2]];
2152 
2153  // Triangle area
2154  Face faceTri( 3 );
2155  faceTri.SetNode( 0, face[0] );
2156  faceTri.SetNode( 1, face[j + 1] );
2157  faceTri.SetNode( 2, face[j + 2] );
2158 
2159  double dTriangleArea = CalculateFaceArea( faceTri, nodes );
2160 
2161  // Calculate the element average
2162  double dTotalSample = 0.0;
2163 
2164  // Loop through all quadrature points
2165  for( int k = 0; k < TriQuadraturePoints; k++ )
2166  {
2167  Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
2168  TriQuadratureG[k][2] * node2.x,
2169  TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
2170  TriQuadratureG[k][2] * node2.y,
2171  TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
2172  TriQuadratureG[k][2] * node2.z );
2173 
2174  double dMagnitude = node.Magnitude();
2175  node.x /= dMagnitude;
2176  node.y /= dMagnitude;
2177  node.z /= dMagnitude;
2178 
2179  double dLon = atan2( node.y, node.x );
2180  if( dLon < 0.0 )
2181  {
2182  dLon += 2.0 * M_PI;
2183  }
2184  double dLat = asin( node.z );
2185 
2186  double dSample = ( *testFunction )( dLon, dLat );
2187 
2188  dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
2189  }
2190 
2191  dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
2192  }
2193  }
2194  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2195  }
2196  else /* discMethod == DiscretizationType_PCLOUD */
2197  {
2198  /* Get the coordinates of the vertices and sample the functionals accordingly */
2199  std::vector< Node >& nodes = trmesh->nodes;
2200 
2201  // Resize the array
2202  dVar.Allocate( nodes.size() );
2203 
2204  for( size_t j = 0; j < nodes.size(); j++ )
2205  {
2206  Node& node = nodes[j];
2207  double dMagnitude = node.Magnitude();
2208  node.x /= dMagnitude;
2209  node.y /= dMagnitude;
2210  node.z /= dMagnitude;
2211  double dLon = atan2( node.y, node.x );
2212  if( dLon < 0.0 )
2213  {
2214  dLon += 2.0 * M_PI;
2215  }
2216  double dLat = asin( node.z );
2217 
2218  double dSample = ( *testFunction )( dLon, dLat );
2219  dVar[j] = dSample;
2220  }
2221 
2222  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2223  }
2224  }
2225 
2226  return moab::MB_SUCCESS;
2227 }

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

Referenced by main().

◆ fill_col_ids()

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

Definition at line 433 of file TempestOnlineMap.hpp.

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

References col_gdofmap, and MB_SUCCESS.

◆ fill_row_ids()

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

Definition at line 423 of file TempestOnlineMap.hpp.

424  {
425  ids_of_interest.reserve( row_gdofmap.size() );
426  // need to add 1
427  for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
428  ids_of_interest.push_back( *it + 1 );
429 
430  return moab::MB_SUCCESS;
431  }

References MB_SUCCESS, and row_gdofmap.

◆ GenerateRemappingWeights()

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

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

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

Definition at line 682 of file TempestOnlineMap.cpp.

687 {
688  NcError error( NcError::silent_nonfatal );
689 
690  moab::DebugOutput dbgprint( std::cout, rank, 0 );
691  dbgprint.set_prefix( "[TempestOnlineMap]: " );
692  moab::ErrorCode rval;
693 
694  const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
695  const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
696  const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
697 
698  // Build a matrix of source and target discretization so that we know how
699  // to assign the global DoFs in parallel for the mapping weights.
700  // For example,
701  // for FV->FV: the rows represented target DoFs and cols represent source DoFs
702  try
703  {
704  // Check command line parameters (data type arguments)
705  STLStringHelper::ToLower( strInputType );
706  STLStringHelper::ToLower( strOutputType );
707 
708  DiscretizationType eInputType;
709  DiscretizationType eOutputType;
710 
711  if( strInputType == "fv" )
712  {
713  eInputType = DiscretizationType_FV;
714  }
715  else if( strInputType == "cgll" )
716  {
717  eInputType = DiscretizationType_CGLL;
718  }
719  else if( strInputType == "dgll" )
720  {
721  eInputType = DiscretizationType_DGLL;
722  }
723  else if( strInputType == "pcloud" )
724  {
725  eInputType = DiscretizationType_PCLOUD;
726  }
727  else
728  {
729  _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
730  }
731 
732  if( strOutputType == "fv" )
733  {
734  eOutputType = DiscretizationType_FV;
735  }
736  else if( strOutputType == "cgll" )
737  {
738  eOutputType = DiscretizationType_CGLL;
739  }
740  else if( strOutputType == "dgll" )
741  {
742  eOutputType = DiscretizationType_DGLL;
743  }
744  else if( strOutputType == "pcloud" )
745  {
746  eOutputType = DiscretizationType_PCLOUD;
747  }
748  else
749  {
750  _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
751  }
752 
753  // set all required input params
754  m_bConserved = !mapOptions.fNoConservation;
755  m_eInputType = eInputType;
756  m_eOutputType = eOutputType;
757 
758  // Method flags
759  std::string strMapAlgorithm( "" );
760  int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
761 
762  // Make an index of method arguments
763  std::set< std::string > setMethodStrings;
764  {
765  int iLast = 0;
766  for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
767  {
768  if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
769  {
770  std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
771  STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
772  if( strMethodString.length() > 0 )
773  {
774  setMethodStrings.insert( strMethodString );
775  }
776  iLast = i + 1;
777  }
778  }
779  }
780 
781  for( auto it : setMethodStrings )
782  {
783  // Piecewise constant monotonicity
784  if( it == "mono2" )
785  {
786  if( nMonotoneType != 0 )
787  {
788  _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
789  }
791  {
792  _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
793  }
794  nMonotoneType = 2;
795 
796  // Piecewise linear monotonicity
797  }
798  else if( it == "mono3" )
799  {
800  if( nMonotoneType != 0 )
801  {
802  _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
803  }
805  {
806  _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
807  }
808  nMonotoneType = 3;
809 
810  // Volumetric remapping from FV to GLL
811  }
812  else if( it == "volumetric" )
813  {
815  {
816  _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
817  }
818  strMapAlgorithm = "volumetric";
819 
820  // Inverse distance mapping
821  }
822  else if( it == "invdist" )
823  {
825  {
826  _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
827  }
828  strMapAlgorithm = "invdist";
829 
830  // Delaunay triangulation mapping
831  }
832  else if( it == "delaunay" )
833  {
835  {
836  _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
837  }
838  strMapAlgorithm = "delaunay";
839 
840  // Bilinear
841  }
842  else if( it == "bilin" )
843  {
845  {
846  _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
847  }
848  strMapAlgorithm = "fvbilin";
849 
850  // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
851  }
852  else if( it == "intbilin" )
853  {
855  {
856  _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
857  }
859  {
860  strMapAlgorithm = "fvintbilin";
861  }
862  else
863  {
864  strMapAlgorithm = "mono3";
865  }
866 
867  // Integrated bilinear with generalized Barycentric coordinates
868  }
869  else if( it == "intbilingb" )
870  {
872  {
873  _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
874  }
875  strMapAlgorithm = "fvintbilingb";
876  }
877  else
878  {
879  _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
880  }
881  }
882 
885  : mapOptions.nPin );
888  : mapOptions.nPout );
889 
890  rval = SetDOFmapTags( srcDofTagName, tgtDofTagName );MB_CHK_ERR( rval );
891 
892  /// the tag should be created already in the e3sm workflow; if not, create it here
893  Tag areaTag;
894  rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
896  if( MB_ALREADY_ALLOCATED == rval )
897  {
898  if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
899  }
900 
901  double dTotalAreaInput = 0.0, dTotalAreaOutput = 0.0;
902  if( !m_bPointCloudSource )
903  {
904  // Calculate Input Mesh Face areas
905  if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
906  double dTotalAreaInput_loc = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
907  rval = m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea );MB_CHK_ERR( rval );
908  dTotalAreaInput = dTotalAreaInput_loc;
909 #ifdef MOAB_HAVE_MPI
910  if( m_pcomm )
911  MPI_Reduce( &dTotalAreaInput_loc, &dTotalAreaInput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
912 #endif
913  if( is_root ) dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput );
914 
915  // Input mesh areas
916  m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
917  // we do not need to set the area on coverage mesh, only on source and target meshes
918  // rval = m_interface->tag_set_data( areaTag, m_remapper->m_covering_source_entities, m_meshInputCov->vecFaceArea );MB_CHK_ERR( rval );
919  }
920 
921  if( !m_bPointCloudTarget )
922  {
923  // Calculate Output Mesh Face areas
924  if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
925  double dTotalAreaOutput_loc = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
926  dTotalAreaOutput = dTotalAreaOutput_loc;
927 #ifdef MOAB_HAVE_MPI
928  if( m_pcomm )
929  MPI_Reduce( &dTotalAreaOutput_loc, &dTotalAreaOutput, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
930 #endif
931  if( is_root ) dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", dTotalAreaOutput );
932  rval = m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea );MB_CHK_ERR( rval );
933  }
934 
935  if( !m_bPointCloud )
936  {
937  // Verify that overlap mesh is in the correct order, only if size == 1
938  if( 1 == size )
939  {
940  int ixSourceFaceMax = ( -1 );
941  int ixTargetFaceMax = ( -1 );
942 
943  if( m_meshOverlap->vecSourceFaceIx.size() != m_meshOverlap->vecTargetFaceIx.size() )
944  {
945  _EXCEPTIONT( "Invalid overlap mesh:\n"
946  " Possible mesh file corruption?" );
947  }
948 
949  for( unsigned i = 0; i < m_meshOverlap->faces.size(); i++ )
950  {
951  if( m_meshOverlap->vecSourceFaceIx[i] + 1 > ixSourceFaceMax )
952  ixSourceFaceMax = m_meshOverlap->vecSourceFaceIx[i] + 1;
953 
954  if( m_meshOverlap->vecTargetFaceIx[i] + 1 > ixTargetFaceMax )
955  ixTargetFaceMax = m_meshOverlap->vecTargetFaceIx[i] + 1;
956  }
957 
958  // Check for forward correspondence in overlap mesh
959  if( m_meshInputCov->faces.size() - ixSourceFaceMax == 0 )
960  {
961  if( is_root ) dbgprint.printf( 0, "Overlap mesh forward correspondence found\n" );
962  }
963  else if( m_meshOutput->faces.size() - ixTargetFaceMax == 0 )
964  { // Check for reverse correspondence in overlap mesh
965  if( is_root ) dbgprint.printf( 0, "Overlap mesh reverse correspondence found (reversing)\n" );
966 
967  // Reorder overlap mesh
968  m_meshOverlap->ExchangeFirstAndSecondMesh();
969  }
970  else
971  { // No correspondence found
972  _EXCEPTION4( "Invalid overlap mesh:\n No correspondence found with input and output meshes "
973  "(%i,%i) vs (%i,%i)",
974  m_meshInputCov->faces.size(), m_meshOutput->faces.size(), ixSourceFaceMax,
975  ixTargetFaceMax );
976  }
977  }
978 
979  // Calculate Face areas
980  if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
981  double dTotalAreaOverlap_loc = m_meshOverlap->CalculateFaceAreas( false );
982  double dTotalAreaOverlap = dTotalAreaOverlap_loc;
983 #ifdef MOAB_HAVE_MPI
984  if( m_pcomm )
985  MPI_Reduce( &dTotalAreaOverlap_loc, &dTotalAreaOverlap, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
986 #endif
987  if( is_root ) dbgprint.printf( 0, "Overlap Mesh Area: %1.15e\n", dTotalAreaOverlap );
988 
989  // Correct areas to match the areas calculated in the overlap mesh
990  // if (fCorrectAreas) // In MOAB-TempestRemap, we will always keep this to be true
991  {
992  if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
993  DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
994  DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
995 
996  assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
997  assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
998  assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
999 
1000  assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
1001  assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
1002 
1003  for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
1004  {
1005  if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
1006  continue; // skip this cell since it is ghosted
1007 
1008  // let us recompute the source/target areas based on overlap mesh areas
1009  assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
1010  dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
1011  assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
1012  dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
1013  }
1014 
1015  for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
1016  {
1017  if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
1018  {
1019  m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
1020  }
1021  }
1022  for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
1023  {
1024  if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
1025  {
1026  m_meshOutput->vecFaceArea[i] = dTargetArea[i];
1027  }
1028  }
1029  }
1030 
1031  // Set source mesh areas in map
1032  if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
1033  {
1034  this->SetSourceAreas( m_meshInputCov->vecFaceArea );
1035  if( m_meshInputCov->vecMask.size() )
1036  {
1037  this->SetSourceMask( m_meshInputCov->vecMask );
1038  }
1039  }
1040 
1041  // Set target mesh areas in map
1042  if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
1043  {
1044  this->SetTargetAreas( m_meshOutput->vecFaceArea );
1045  if( m_meshOutput->vecMask.size() )
1046  {
1047  this->SetTargetMask( m_meshOutput->vecMask );
1048  }
1049  }
1050 
1051  /*
1052  // Recalculate input mesh area from overlap mesh
1053  if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
1054  dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
1055  dbgprint.printf(0, "Recalculating source mesh areas\n");
1056  dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
1057  dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
1058  }
1059  */
1060  }
1061 
1062  this->m_pmeshSource = m_meshInputCov;
1063  this->m_pmeshOverlap = m_meshOverlap;
1064 
1065  // Finite volume input / Finite volume output
1066  if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1067  {
1068  // Generate reverse node array and edge map
1069  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1070  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1071 
1072  // Initialize coordinates for map
1073  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1074  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1075 
1076  this->m_pdataGLLNodesIn = nullptr;
1077  this->m_pdataGLLNodesOut = nullptr;
1078 
1079  // Finite volume input / Finite element output
1080  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
1081  mapOptions.nPout, false, nullptr );MB_CHK_ERR( rval );
1082 
1083  // Construct remap for FV-FV
1084  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1085 
1086  // Construct OfflineMap
1087  if( strMapAlgorithm == "invdist" )
1088  {
1089  if( is_root ) dbgprint.printf( 0, "Calculating map (invdist)\n" );
1090  if( m_meshInputCov->faces.size() )
1091  LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1092  }
1093  else if( strMapAlgorithm == "delaunay" )
1094  {
1095  if( is_root ) dbgprint.printf( 0, "Calculating map (delaunay)\n" );
1096  if( m_meshInputCov->faces.size() )
1097  LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1098  }
1099  else if( strMapAlgorithm == "fvintbilin" )
1100  {
1101  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilin)\n" );
1102  if( m_meshInputCov->faces.size() )
1103  LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1104  }
1105  else if( strMapAlgorithm == "fvintbilingb" )
1106  {
1107  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilingb)\n" );
1108  if( m_meshInputCov->faces.size() )
1109  LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1110  *this );
1111  }
1112  else if( strMapAlgorithm == "fvbilin" )
1113  {
1114 #ifdef VERBOSE
1115  if( is_root )
1116  {
1117  m_meshInputCov->Write( "SourceMeshMBTR.g" );
1118  m_meshOutput->Write( "TargetMeshMBTR.g" );
1119  }
1120  else
1121  {
1122  m_meshInputCov->Write( "SourceMeshMBTR" + std::to_string( rank ) + ".g" );
1123  m_meshOutput->Write( "TargetMeshMBTR" + std::to_string( rank ) + ".g" );
1124  }
1125 #endif
1126  if( is_root ) dbgprint.printf( 0, "Calculating map (bilin)\n" );
1127  if( m_meshInputCov->faces.size() )
1128  LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
1129  }
1130  else
1131  {
1132  if( is_root ) dbgprint.printf( 0, "Calculating conservative FV-FV map\n" );
1133  if( m_meshInputCov->faces.size() )
1134  {
1135 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1136  LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1137  ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
1138 #else
1139  LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
1140 #endif
1141  }
1142  }
1143  }
1144  else if( eInputType == DiscretizationType_FV )
1145  {
1146  DataArray3D< double > dataGLLJacobian;
1147 
1148  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1149  double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
1150  dataGLLNodesDest, dataGLLJacobian );
1151 
1152  double dNumericalArea = dNumericalArea_loc;
1153 #ifdef MOAB_HAVE_MPI
1154  if( m_pcomm )
1155  MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1156 #endif
1157  if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
1158 
1159  // Initialize coordinates for map
1160  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1161  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1162 
1163  this->m_pdataGLLNodesIn = nullptr;
1164  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1165 
1166  // Generate the continuous Jacobian
1167  bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
1168 
1169  if( eOutputType == DiscretizationType_CGLL )
1170  {
1171  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
1172  }
1173  else
1174  {
1175  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
1176  }
1177 
1178  // Generate reverse node array and edge map
1179  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1180  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1181 
1182  // Finite volume input / Finite element output
1183  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
1184  mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1185  &dataGLLNodesDest );MB_CHK_ERR( rval );
1186 
1187  // Generate remap weights
1188  if( strMapAlgorithm == "volumetric" )
1189  {
1190  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
1191  LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
1192  dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
1193  nMonotoneType, fContinuous, mapOptions.fNoConservation );
1194  }
1195  else
1196  {
1197  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
1198  LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
1199  this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
1200  mapOptions.fNoConservation );
1201  }
1202  }
1203  else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
1204  {
1205  DataArray3D< double > dataGLLJacobian;
1206  if( !m_bPointCloudSource )
1207  {
1208  // Generate reverse node array and edge map
1209  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1210  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
1211 
1212  // Initialize coordinates for map
1213  if( eInputType == DiscretizationType_FV )
1214  {
1215  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1216  }
1217  else
1218  {
1219  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1220  DataArray3D< double > dataGLLJacobianSrc;
1221  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1222  dataGLLJacobian );
1223  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1224  dataGLLJacobianSrc );
1225  }
1226  }
1227  // else { /* Source is a point cloud dataset */ }
1228 
1229  if( !m_bPointCloudTarget )
1230  {
1231  // Generate reverse node array and edge map
1232  if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
1233  if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
1234 
1235  // Initialize coordinates for map
1236  if( eOutputType == DiscretizationType_FV )
1237  {
1238  this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
1239  }
1240  else
1241  {
1242  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1243  GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1244  dataGLLJacobian );
1245  }
1246  }
1247  // else { /* Target is a point cloud dataset */ }
1248 
1249  // Finite volume input / Finite element output
1250  rval = this->SetDOFmapAssociation(
1251  eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1252  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrcCov ),
1253  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrc ),
1254  eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1255  ( m_bPointCloudTarget ? nullptr : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
1256 
1257  // Construct remap
1258  if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
1259  rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
1260  }
1261  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1262  {
1263  DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
1264 
1265  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1266  // double dNumericalAreaCov_loc =
1267  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1268  dataGLLJacobian );
1269 
1270  double dNumericalArea_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
1271  dataGLLNodesSrc, dataGLLJacobianSrc );
1272 
1273  // if ( is_root ) dbgprint.printf ( 0, "Input Mesh: Coverage Area: %1.15e, Output Area:
1274  // %1.15e\n", dNumericalAreaCov_loc, dTotalAreaOutput_loc );
1275  // assert(dNumericalAreaCov_loc >= dTotalAreaOutput_loc);
1276 
1277  double dNumericalArea = dNumericalArea_loc;
1278 #ifdef MOAB_HAVE_MPI
1279  if( m_pcomm )
1280  MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1281 #endif
1282  if( is_root )
1283  {
1284  dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalArea );
1285 
1286  if( fabs( dNumericalArea - dTotalAreaInput ) > 1.0e-12 )
1287  {
1288  dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
1289  "numerical area and geometric area\n" );
1290  }
1291  }
1292 
1293  if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
1294  {
1295  _EXCEPTIONT( "Number of element does not match between metadata and "
1296  "input mesh" );
1297  }
1298 
1299  // Initialize coordinates for map
1300  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1301  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1302 
1303  // Generate the continuous Jacobian for input mesh
1304  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1305 
1306  if( eInputType == DiscretizationType_CGLL )
1307  {
1308  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
1309  }
1310  else
1311  {
1312  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
1313  }
1314 
1315  // Finite element input / Finite volume output
1316  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1317  &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1318  false, nullptr );MB_CHK_ERR( rval );
1319 
1320  // Generate remap
1321  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1322 
1323  if( strMapAlgorithm == "volumetric" )
1324  {
1325  _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
1326  "GLL input mesh" );
1327  }
1328 
1329  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1330  this->m_pdataGLLNodesOut = nullptr;
1331 
1332 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1333  LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
1334  nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
1335  *this );
1336 #else
1337  LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
1338  mapOptions.fNoConservation );
1339 #endif
1340  }
1341  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
1342  {
1343  DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
1344  DataArray3D< double > dataGLLJacobianOut;
1345 
1346  // Input metadata
1347  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1348  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1349  dataGLLJacobianIn );
1350 
1351  double dNumericalAreaSrc_loc = GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble,
1352  dataGLLNodesSrc, dataGLLJacobianSrc );
1353  double dNumericalAreaIn = dNumericalAreaSrc_loc;
1354 #ifdef MOAB_HAVE_MPI
1355  if( m_pcomm )
1356  MPI_Reduce( &dNumericalAreaSrc_loc, &dNumericalAreaIn, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1357 #endif
1358  if( is_root )
1359  {
1360  dbgprint.printf( 0, "Input Mesh Numerical Area: %1.15e\n", dNumericalAreaIn );
1361 
1362  if( fabs( dNumericalAreaIn - dTotalAreaInput ) > 1.0e-12 )
1363  {
1364  dbgprint.printf( 0, "WARNING: Significant mismatch between input mesh "
1365  "numerical area and geometric area\n" );
1366  }
1367  }
1368 
1369  // Output metadata
1370  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1371  double dNumericalAreaOut_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
1372  dataGLLNodesDest, dataGLLJacobianOut );
1373  double dNumericalAreaOut = dNumericalAreaOut_loc;
1374 #ifdef MOAB_HAVE_MPI
1375  if( m_pcomm )
1376  MPI_Reduce( &dNumericalAreaOut_loc, &dNumericalAreaOut, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1377 #endif
1378  if( is_root )
1379  {
1380  dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalAreaOut );
1381 
1382  if( fabs( dNumericalAreaOut - dTotalAreaOutput ) > 1.0e-12 )
1383  {
1384  if( is_root )
1385  dbgprint.printf( 0, "WARNING: Significant mismatch between output mesh "
1386  "numerical area and geometric area\n" );
1387  }
1388  }
1389 
1390  // Initialize coordinates for map
1391  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1392  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1393 
1394  // Generate the continuous Jacobian for input mesh
1395  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1396 
1397  if( eInputType == DiscretizationType_CGLL )
1398  {
1399  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
1400  }
1401  else
1402  {
1403  GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
1404  }
1405 
1406  // Generate the continuous Jacobian for output mesh
1407  bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
1408 
1409  if( eOutputType == DiscretizationType_CGLL )
1410  {
1411  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
1412  }
1413  else
1414  {
1415  GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
1416  }
1417 
1418  // Input Finite Element to Output Finite Element
1419  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1420  &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1421  ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
1422 
1423  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1424  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1425 
1426  // Generate remap
1427  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1428 
1429 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1430  LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
1431  dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1432  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1433  fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *this );
1434 #else
1435  LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1436  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1437  fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
1438 #endif
1439  }
1440  else
1441  {
1442  _EXCEPTIONT( "Not implemented" );
1443  }
1444 
1445 #ifdef MOAB_HAVE_EIGEN3
1446  copy_tempest_sparsemat_to_eigen3();
1447 #endif
1448 
1449 #ifdef MOAB_HAVE_MPI
1450  {
1451  // Remove ghosted entities from overlap set
1452  moab::Range ghostedEnts;
1453  rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
1455  rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
1456  }
1457 #endif
1458  // Verify consistency, conservation and monotonicity, globally
1459  if( !mapOptions.fNoCheck )
1460  {
1461  if( is_root ) dbgprint.printf( 0, "Verifying map" );
1462  this->IsConsistent( 1.0e-8 );
1463  if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
1464 
1465  if( nMonotoneType != 0 )
1466  {
1467  this->IsMonotone( 1.0e-12 );
1468  }
1469  }
1470  }
1471  catch( Exception& e )
1472  {
1473  dbgprint.printf( 0, "%s", e.ToString().c_str() );
1474  return ( moab::MB_FAILURE );
1475  }
1476  catch( ... )
1477  {
1478  return ( moab::MB_FAILURE );
1479  }
1480  return moab::MB_SUCCESS;
1481 }

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

Referenced by main().

◆ GetColGlobalDoF()

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

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

Definition at line 557 of file TempestOnlineMap.hpp.

558 {
559  return col_gdofmap[localColID];
560 }

◆ GetDestinationGlobalNDofs()

int moab::TempestOnlineMap::GetDestinationGlobalNDofs ( )

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

◆ GetDestinationLocalNDofs()

int moab::TempestOnlineMap::GetDestinationLocalNDofs ( )

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

◆ GetDestinationNDofsPerElement()

int moab::TempestOnlineMap::GetDestinationNDofsPerElement ( )
inline

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

Definition at line 575 of file TempestOnlineMap.hpp.

576 {
577  return m_nDofsPEl_Dest;
578 }

◆ GetGlobalSourceAreas()

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

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

◆ GetGlobalTargetAreas()

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

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

◆ GetIndexOfColGlobalDoF()

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

Get the index of globaColDoF.

Definition at line 562 of file TempestOnlineMap.hpp.

563 {
564  return globalColDoF + 1; // temporary
565 }

◆ GetIndexOfRowGlobalDoF()

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

Get the index of globaRowDoF.

Definition at line 551 of file TempestOnlineMap.hpp.

552 {
553  return globalRowDoF + 1;
554 }

◆ GetRowGlobalDoF()

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

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

Definition at line 546 of file TempestOnlineMap.hpp.

547 {
548  return row_gdofmap[localRowID];
549 }

◆ GetSourceGlobalNDofs()

int moab::TempestOnlineMap::GetSourceGlobalNDofs ( )

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

◆ GetSourceLocalNDofs()

int moab::TempestOnlineMap::GetSourceLocalNDofs ( )

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

◆ GetSourceNDofsPerElement()

int moab::TempestOnlineMap::GetSourceNDofsPerElement ( )
inline

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

Definition at line 568 of file TempestOnlineMap.hpp.

569 {
570  return m_nDofsPEl_Src;
571 }

◆ IsConservative()

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

Determine if the map is conservative.

Definition at line 1531 of file TempestOnlineMap.cpp.

1532 {
1533 #ifndef MOAB_HAVE_MPI
1534 
1535  return OfflineMap::IsConservative( dTolerance );
1536 
1537 #else
1538  // return OfflineMap::IsConservative(dTolerance);
1539 
1540  int ierr;
1541  // Get map entries
1542  DataArray1D< int > dataRows;
1543  DataArray1D< int > dataCols;
1544  DataArray1D< double > dataEntries;
1545  const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1546  const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1547 
1548  // Calculate column sums
1549  std::vector< int > dColumnsUnique;
1550  std::vector< double > dColumnSums;
1551 
1552  int nColumns = m_mapRemap.GetColumns();
1553  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1554  dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1555  dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1556 
1557  for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
1558  {
1559  dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1560 
1561  assert( dataCols[i] < m_nTotDofs_SrcCov );
1562 
1563  // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
1564  int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1565  // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1566  dColumnsUnique[dataCols[i]] = colGID;
1567 
1568  // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
1569  // std::endl;
1570  }
1571 
1572  int rootProc = 0;
1573  std::vector< int > nElementsInProc;
1574  const int nDATA = 3;
1575  nElementsInProc.resize( size * nDATA );
1576  int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1577  ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1578  if( ierr != MPI_SUCCESS ) return -1;
1579 
1580  int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
1581  std::vector< int > dColumnIndices;
1582  std::vector< double > dColumnSourceAreas;
1583  std::vector< double > dColumnSumsTotal;
1584  std::vector< int > displs, rcount;
1585  if( rank == rootProc )
1586  {
1587  displs.resize( size + 1, 0 );
1588  rcount.resize( size, 0 );
1589  int gsum = 0;
1590  for( int ir = 0; ir < size; ++ir )
1591  {
1592  nTotVals += nElementsInProc[ir * nDATA];
1593  nTotColumns += nElementsInProc[ir * nDATA + 1];
1594  nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1595 
1596  displs[ir] = gsum;
1597  rcount[ir] = nElementsInProc[ir * nDATA + 1];
1598  gsum += rcount[ir];
1599 
1600  // printf( "%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns, displs[ir], rcount[ir], gsum );
1601  }
1602 
1603  printf( "Total nnz: %d, global source elements = %d\n", nTotVals, gsum );
1604 
1605  dColumnIndices.resize( nTotColumns, -1 );
1606  dColumnSumsTotal.resize( nTotColumns, 0.0 );
1607  // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
1608  }
1609 
1610  // Gather all ColumnSums to root process and accumulate
1611  // We expect that the sums of all columns equate to 1.0 within user specified tolerance
1612  // Need to do a gatherv here since different processes have different number of elements
1613  // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
1614  // MPI_SUM, 0, m_pcomm->comm());
1615  ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
1616  displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1617  if( ierr != MPI_SUCCESS ) return -1;
1618  ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
1619  displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1620  if( ierr != MPI_SUCCESS ) return -1;
1621  // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
1622  // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
1623  // MPI_SUCCESS ) return -1;
1624 
1625  // Clean out unwanted arrays now
1626  dColumnSums.clear();
1627  dColumnsUnique.clear();
1628 
1629  // Verify all column sums equal the input Jacobian
1630  int fConservative = 0;
1631  if( rank == rootProc )
1632  {
1633  displs[size] = ( nTotColumns );
1634  // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
1635  std::map< int, double > dColumnSumsOnRoot;
1636  // std::map<int, double> dColumnSourceAreasOnRoot;
1637  for( int ir = 0; ir < size; ir++ )
1638  {
1639  for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1640  {
1641  if( dColumnIndices[ips] < 0 ) continue;
1642  // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
1643  assert( dColumnIndices[ips] < nTotColumnsUnq );
1644  dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
1645  // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
1646  // dColumnSourceAreas[ dColumnIndices[ips] ]
1647  }
1648  }
1649 
1650  for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1651  {
1652  // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
1653  if( fabs( it->second - 1.0 ) > dTolerance )
1654  {
1655  fConservative++;
1656  Announce( "TempestOnlineMap is not conservative in column "
1657  // "%i (%1.15e)", it->first, it->second );
1658  "%i (%1.15e)",
1659  it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
1660  }
1661  }
1662  }
1663 
1664  // TODO: Just do a broadcast from root instead of a reduction
1665  ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1666  if( ierr != MPI_SUCCESS ) return -1;
1667 
1668  return fConservative;
1669 #endif
1670 }

References size.

◆ IsConsistent()

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

Determine if the map is first-order accurate.

Definition at line 1485 of file TempestOnlineMap.cpp.

1486 {
1487 #ifndef MOAB_HAVE_MPI
1488 
1489  return OfflineMap::IsConsistent( dTolerance );
1490 
1491 #else
1492 
1493  // Get map entries
1494  DataArray1D< int > dataRows;
1495  DataArray1D< int > dataCols;
1496  DataArray1D< double > dataEntries;
1497 
1498  // Calculate row sums
1499  DataArray1D< double > dRowSums;
1500  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1501  dRowSums.Allocate( m_mapRemap.GetRows() );
1502 
1503  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1504  {
1505  dRowSums[dataRows[i]] += dataEntries[i];
1506  }
1507 
1508  // Verify all row sums are equal to 1
1509  int fConsistent = 0;
1510  for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1511  {
1512  if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1513  {
1514  fConsistent++;
1515  int rowGID = row_gdofmap[i];
1516  Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1517  }
1518  }
1519 
1520  int ierr;
1521  int fConsistentGlobal = 0;
1522  ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1523  if( ierr != MPI_SUCCESS ) return -1;
1524 
1525  return fConsistentGlobal;
1526 #endif
1527 }

◆ IsMonotone()

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

Determine if the map is monotone.

Definition at line 1674 of file TempestOnlineMap.cpp.

1675 {
1676 #ifndef MOAB_HAVE_MPI
1677 
1678  return OfflineMap::IsMonotone( dTolerance );
1679 
1680 #else
1681 
1682  // Get map entries
1683  DataArray1D< int > dataRows;
1684  DataArray1D< int > dataCols;
1685  DataArray1D< double > dataEntries;
1686 
1687  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1688 
1689  // Verify all entries are in the range [0,1]
1690  int fMonotone = 0;
1691  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1692  {
1693  if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1694  {
1695  fMonotone++;
1696 
1697  Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1698  }
1699  }
1700 
1701  int ierr;
1702  int fMonotoneGlobal = 0;
1703  ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1704  if( ierr != MPI_SUCCESS ) return -1;
1705 
1706  return fMonotoneGlobal;
1707 #endif
1708 }

◆ LinearRemapFVtoFV_Tempest_MOAB()

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

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

Definition at line 121 of file TempestLinearRemap.cpp.

122 {
123  // Order of triangular quadrature rule
124  const int TriQuadRuleOrder = 4;
125 
126  // Verify ReverseNodeArray has been calculated
127  if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
128  {
129  _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInputCov" );
130  }
131 
132  // Triangular quadrature rule
133  TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
134 
135  // Number of coefficients needed at this order
136 #ifdef RECTANGULAR_TRUNCATION
137  int nCoefficients = nOrder * nOrder;
138 #endif
139 #ifdef TRIANGULAR_TRUNCATION
140  int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
141 #endif
142 
143  // Number of faces you need
144  const int nRequiredFaceSetSize = nCoefficients;
145 
146  // Fit weight exponent
147  const int nFitWeightsExponent = nOrder + 2;
148 
149  // Announcements
150  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
151  dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
152  if( is_root )
153  {
154  dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
155  dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
156  dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
157  dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
158  dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
159  }
160 
161  // Current overlap face
162  int ixOverlap = 0;
163 #ifdef VERBOSE
164  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
165 #endif
166  DataArray2D< double > dIntArray;
167  DataArray1D< double > dConstraint( nCoefficients );
168 
169  // Loop through all faces on m_meshInputCov
170  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
171  {
172  // Output every 1000 elements
173 #ifdef VERBOSE
174  if( ixFirst % outputFrequency == 0 && is_root )
175  {
176  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
177  }
178 #endif
179  // Find the set of Faces that overlap faceFirst
180  int ixOverlapBegin = ixOverlap;
181  unsigned ixOverlapEnd = ixOverlapBegin;
182 
183  for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
184  {
185  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 ) break;
186  }
187 
188  unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
189 
190  if( nOverlapFaces == 0 ) continue;
191 
192  // Build integration array
193  BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
194  nOrder, dIntArray );
195 
196  // Set of Faces to use in building the reconstruction and associated
197  // distance metric.
198  AdjacentFaceVector vecAdjFaces;
199 
200  GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
201 
202  // Number of adjacent Faces
203  int nAdjFaces = vecAdjFaces.size();
204 
205  // Determine the conservative constraint equation
206  double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
207  dConstraint.Zero();
208  for( int p = 0; p < nCoefficients; p++ )
209  {
210  for( unsigned j = 0; j < nOverlapFaces; j++ )
211  {
212  dConstraint[p] += dIntArray[p][j];
213  }
214  dConstraint[p] /= dFirstArea;
215  }
216 
217  // Build the fit array from the integration operator
218  DataArray2D< double > dFitArray;
219  DataArray1D< double > dFitWeights;
220  DataArray2D< double > dFitArrayPlus;
221 
222  BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
223  dFitArray, dFitWeights );
224 
225  // Compute the inverse fit array
226  bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
227 
228  // Multiply integration array and fit array
229  DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
230  if( fSuccess )
231  {
232  // Multiply integration array and inverse fit array
233  for( int i = 0; i < nAdjFaces; i++ )
234  {
235  for( size_t j = 0; j < nOverlapFaces; j++ )
236  {
237  for( int k = 0; k < nCoefficients; k++ )
238  {
239  dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
240  }
241  }
242  }
243 
244  // Unable to invert fit array, drop to 1st order. In this case
245  // dFitArrayPlus(0,0) = 1 and all other entries are zero.
246  }
247  else
248  {
249  dComposedArray.Zero();
250  for( size_t j = 0; j < nOverlapFaces; j++ )
251  {
252  dComposedArray( 0, j ) += dIntArray( 0, j );
253  }
254  }
255 
256  // Put composed array into map
257  for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
258  {
259  for( unsigned j = 0; j < nOverlapFaces; j++ )
260  {
261  int& ixFirstFaceLoc = vecAdjFaces[i].first;
262  int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
263  // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
264  // ixFirstFaceLoc); int ixSecondFaceGlob =
265  // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
266 
267  // signal to not participate, because it is a ghost target
268  if( ixSecondFaceLoc < 0 ) continue; // do not do anything
269 
270  m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
271  dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
272  }
273  }
274 
275  // Increment the current overlap index
276  ixOverlap += nOverlapFaces;
277  }
278 
279  return;
280 }

References dbgprint.

◆ LinearRemapFVtoGLL_MOAB()

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

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

◆ LinearRemapGLLtoGLL2_MOAB()

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

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

Definition at line 1402 of file TempestLinearRemap.cpp.

1413 {
1414  // Triangular quadrature rule
1415  TriangularQuadratureRule triquadrule( 8 );
1416 
1417  const DataArray2D< double >& dG = triquadrule.GetG();
1418  const DataArray1D< double >& dW = triquadrule.GetW();
1419 
1420  // Get SparseMatrix represntation of the OfflineMap
1421  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1422 
1423  // Sample coefficients
1424  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1425  DataArray2D< double > dSampleCoeffOut( nPout, nPout );
1426 
1427  // Announcemnets
1428  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1429  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
1430  if( is_root )
1431  {
1432  dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
1433  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1434  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1435  }
1436 
1437  // Build the integration array for each element on m_meshOverlap
1438  DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
1439 
1440  // Number of overlap Faces per source Face
1441  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1442 
1443  int ixOverlap = 0;
1444  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1445  {
1446  // Determine how many overlap Faces and triangles are present
1447  int nOverlapFaces = 0;
1448  size_t ixOverlapTemp = ixOverlap;
1449  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1450  {
1451  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1452  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
1453  {
1454  break;
1455  }
1456 
1457  nOverlapFaces++;
1458  }
1459 
1460  nAllOverlapFaces[ixFirst] = nOverlapFaces;
1461 
1462  // Increment the current overlap index
1463  ixOverlap += nAllOverlapFaces[ixFirst];
1464  }
1465 
1466  // Geometric area of each output node
1467  DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
1468 
1469  // Area of each overlap element in the output basis
1470  DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
1471 
1472  // Loop through all faces on m_meshInputCov
1473  ixOverlap = 0;
1474 #ifdef VERBOSE
1475  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1476 #endif
1477  if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
1478 
1479  // generic triangle used for area computation, for triangles around the center of overlap face;
1480  // used for overlap faces with more than 4 edges;
1481  // nodes array will be set for each triangle;
1482  // these triangles are not part of the mesh structure, they are just temporary during
1483  // aforementioned decomposition.
1484  Face faceTri( 3 );
1485  NodeVector nodes( 3 );
1486  faceTri.SetNode( 0, 0 );
1487  faceTri.SetNode( 1, 1 );
1488  faceTri.SetNode( 2, 2 );
1489 
1490  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1491  {
1492 #ifdef VERBOSE
1493  // Announce computation progress
1494  if( ixFirst % outputFrequency == 0 && is_root )
1495  {
1496  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1497  }
1498 #endif
1499  // Quantities from the First Mesh
1500  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1501 
1502  const NodeVector& nodesFirst = m_meshInputCov->nodes;
1503 
1504  // Number of overlapping Faces and triangles
1505  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1506 
1507  if( !nOverlapFaces ) continue;
1508 
1509  // // Calculate total element Jacobian
1510  // double dTotalJacobian = 0.0;
1511  // for (int s = 0; s < nPin; s++) {
1512  // for (int t = 0; t < nPin; t++) {
1513  // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
1514  // }
1515  // }
1516 
1517  // Loop through all Overlap Faces
1518  for( int i = 0; i < nOverlapFaces; i++ )
1519  {
1520  // Quantities from the overlap Mesh
1521  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
1522 
1523  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1524 
1525  // Quantities from the Second Mesh
1526  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1527 
1528  // signal to not participate, because it is a ghost target
1529  if( ixSecond < 0 ) continue; // do not do anything
1530 
1531  const NodeVector& nodesSecond = m_meshOutput->nodes;
1532 
1533  const Face& faceSecond = m_meshOutput->faces[ixSecond];
1534 
1535  int nbEdges = faceOverlap.edges.size();
1536  int nOverlapTriangles = 1;
1537  Node center; // not used if nbEdges == 3
1538  if( nbEdges > 3 )
1539  { // decompose from center in this case
1540  nOverlapTriangles = nbEdges;
1541  for( int k = 0; k < nbEdges; k++ )
1542  {
1543  const Node& node = nodesOverlap[faceOverlap[k]];
1544  center = center + node;
1545  }
1546  center = center / nbEdges;
1547  center = center.Normalized(); // project back on sphere of radius 1
1548  }
1549 
1550  Node node0, node1, node2;
1551  double dTriArea;
1552 
1553  // Loop over all sub-triangles of this Overlap Face
1554  for( int j = 0; j < nOverlapTriangles; j++ )
1555  {
1556  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1557  {
1558  node0 = nodesOverlap[faceOverlap[0]];
1559  node1 = nodesOverlap[faceOverlap[1]];
1560  node2 = nodesOverlap[faceOverlap[2]];
1561  dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1562  }
1563  else // decompose polygon in triangles around the center
1564  {
1565  node0 = center;
1566  node1 = nodesOverlap[faceOverlap[j]];
1567  int j1 = ( j + 1 ) % nbEdges;
1568  node2 = nodesOverlap[faceOverlap[j1]];
1569  nodes[0] = center;
1570  nodes[1] = node1;
1571  nodes[2] = node2;
1572  dTriArea = CalculateFaceArea( faceTri, nodes );
1573  }
1574 
1575  for( int k = 0; k < triquadrule.GetPoints(); k++ )
1576  {
1577  // Get the nodal location of this point
1578  double dX[3];
1579 
1580  dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
1581  dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
1582  dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
1583 
1584  double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
1585 
1586  dX[0] /= dMag;
1587  dX[1] /= dMag;
1588  dX[2] /= dMag;
1589 
1590  Node nodeQuadrature( dX[0], dX[1], dX[2] );
1591 
1592  // Find the components of this quadrature point in the basis
1593  // of the first Face.
1594  double dAlphaIn;
1595  double dBetaIn;
1596 
1597  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
1598 
1599  // Find the components of this quadrature point in the basis
1600  // of the second Face.
1601  double dAlphaOut;
1602  double dBetaOut;
1603 
1604  ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
1605 
1606  /*
1607  // Check inverse map value
1608  if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
1609  (dBetaIn < 0.0) || (dBetaIn > 1.0)
1610  ) {
1611  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1612  dAlphaIn, dBetaIn);
1613  }
1614 
1615  // Check inverse map value
1616  if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
1617  (dBetaOut < 0.0) || (dBetaOut > 1.0)
1618  ) {
1619  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1620  dAlphaOut, dBetaOut);
1621  }
1622  */
1623  // Sample the First finite element at this point
1624  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1625 
1626  // Sample the Second finite element at this point
1627  SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
1628 
1629  // Overlap output area
1630  for( int s = 0; s < nPout; s++ )
1631  {
1632  for( int t = 0; t < nPout; t++ )
1633  {
1634  double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
1635 
1636  dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
1637 
1638  dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
1639  }
1640  }
1641 
1642  // Compute overlap integral
1643  int ixp = 0;
1644  for( int p = 0; p < nPin; p++ )
1645  {
1646  for( int q = 0; q < nPin; q++ )
1647  {
1648  int ixs = 0;
1649  for( int s = 0; s < nPout; s++ )
1650  {
1651  for( int t = 0; t < nPout; t++ )
1652  {
1653  // Sample the Second finite element at this point
1654  dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
1655  dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
1656 
1657  ixs++;
1658  }
1659  }
1660 
1661  ixp++;
1662  }
1663  }
1664  }
1665  }
1666  }
1667 
1668  // Coefficients
1669  DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
1670 
1671  for( int i = 0; i < nOverlapFaces; i++ )
1672  {
1673  // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1674 
1675  int ixp = 0;
1676  for( int p = 0; p < nPin; p++ )
1677  {
1678  for( int q = 0; q < nPin; q++ )
1679  {
1680  int ixs = 0;
1681  for( int s = 0; s < nPout; s++ )
1682  {
1683  for( int t = 0; t < nPout; t++ )
1684  {
1685  dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
1686  dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1687 
1688  ixs++;
1689  }
1690  }
1691 
1692  ixp++;
1693  }
1694  }
1695  }
1696 
1697  // Source areas
1698  DataArray1D< double > vecSourceArea( nPin * nPin );
1699 
1700  for( int p = 0; p < nPin; p++ )
1701  {
1702  for( int q = 0; q < nPin; q++ )
1703  {
1704  vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
1705  }
1706  }
1707 
1708  // Target areas
1709  DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
1710 
1711  for( int i = 0; i < nOverlapFaces; i++ )
1712  {
1713  // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1714  int ixs = 0;
1715  for( int s = 0; s < nPout; s++ )
1716  {
1717  for( int t = 0; t < nPout; t++ )
1718  {
1719  vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
1720 
1721  ixs++;
1722  }
1723  }
1724  }
1725 
1726  // Force consistency and conservation
1727  if( !fNoConservation )
1728  {
1729  ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
1730  }
1731 
1732  // Update global coefficients
1733  for( int i = 0; i < nOverlapFaces; i++ )
1734  {
1735  int ixp = 0;
1736  for( int p = 0; p < nPin; p++ )
1737  {
1738  for( int q = 0; q < nPin; q++ )
1739  {
1740  int ixs = 0;
1741  for( int s = 0; s < nPout; s++ )
1742  {
1743  for( int t = 0; t < nPout; t++ )
1744  {
1745  dGlobalIntArray[ixp][ixOverlap + i][ixs] =
1746  dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1747 
1748  ixs++;
1749  }
1750  }
1751 
1752  ixp++;
1753  }
1754  }
1755  }
1756 
1757 #ifdef VVERBOSE
1758  // Check column sums (conservation)
1759  for( int i = 0; i < nPin * nPin; i++ )
1760  {
1761  double dColSum = 0.0;
1762  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1763  {
1764  dColSum += dCoeff[j][i] * vecTargetArea[j];
1765  }
1766  printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
1767  }
1768 
1769  // Check row sums (consistency)
1770  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1771  {
1772  double dRowSum = 0.0;
1773  for( int i = 0; i < nPin * nPin; i++ )
1774  {
1775  dRowSum += dCoeff[j][i];
1776  }
1777  printf( "Row %i: %1.15e\n", j, dRowSum );
1778  }
1779 #endif
1780 
1781  // Increment the current overlap index
1782  ixOverlap += nOverlapFaces;
1783  }
1784 
1785  // Build redistribution map within target element
1786  if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
1787  DataArray1D< double > dRedistSourceArea( nPout * nPout );
1788  DataArray1D< double > dRedistTargetArea( nPout * nPout );
1789  std::vector< DataArray2D< double > > dRedistributionMaps;
1790  dRedistributionMaps.resize( m_meshOutput->faces.size() );
1791 
1792  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1793  {
1794  dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
1795 
1796  for( int i = 0; i < nPout * nPout; i++ )
1797  {
1798  dRedistributionMaps[ixSecond][i][i] = 1.0;
1799  }
1800 
1801  for( int s = 0; s < nPout * nPout; s++ )
1802  {
1803  dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
1804  }
1805 
1806  for( int s = 0; s < nPout * nPout; s++ )
1807  {
1808  dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
1809  }
1810 
1811  if( !fNoConservation )
1812  {
1813  ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
1814  ( nMonotoneType != 0 ) );
1815 
1816  for( int s = 0; s < nPout * nPout; s++ )
1817  {
1818  for( int t = 0; t < nPout * nPout; t++ )
1819  {
1820  dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
1821  }
1822  }
1823  }
1824  }
1825 
1826  // Construct the total geometric area
1827  DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
1828  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1829  {
1830  for( int s = 0; s < nPout; s++ )
1831  {
1832  for( int t = 0; t < nPout; t++ )
1833  {
1834  dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
1835  dGeometricOutputArea[ixSecond][s * nPout + t];
1836  }
1837  }
1838  }
1839 
1840  // Compose the integration operator with the output map
1841  ixOverlap = 0;
1842 
1843  if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
1844 
1845  // Map from source DOFs to target DOFs with redistribution applied
1846  DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
1847 
1848  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1849  {
1850 #ifdef VERBOSE
1851  // Announce computation progress
1852  if( ixFirst % outputFrequency == 0 && is_root )
1853  {
1854  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1855  }
1856 #endif
1857  // Number of overlapping Faces and triangles
1858  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1859 
1860  if( !nOverlapFaces ) continue;
1861 
1862  // Put composed array into map
1863  for( int j = 0; j < nOverlapFaces; j++ )
1864  {
1865  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1866 
1867  // signal to not participate, because it is a ghost target
1868  if( ixSecondFace < 0 ) continue; // do not do anything
1869 
1870  dRedistributedOp.Zero();
1871  for( int p = 0; p < nPin * nPin; p++ )
1872  {
1873  for( int s = 0; s < nPout * nPout; s++ )
1874  {
1875  for( int t = 0; t < nPout * nPout; t++ )
1876  {
1877  dRedistributedOp[p][s] +=
1878  dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
1879  }
1880  }
1881  }
1882 
1883  int ixp = 0;
1884  for( int p = 0; p < nPin; p++ )
1885  {
1886  for( int q = 0; q < nPin; q++ )
1887  {
1888  int ixFirstNode;
1889  if( fContinuousIn )
1890  {
1891  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
1892  }
1893  else
1894  {
1895  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
1896  }
1897 
1898  int ixs = 0;
1899  for( int s = 0; s < nPout; s++ )
1900  {
1901  for( int t = 0; t < nPout; t++ )
1902  {
1903  int ixSecondNode;
1904  if( fContinuousOut )
1905  {
1906  ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
1907 
1908  if( !fNoConservation )
1909  {
1910  smatMap( ixSecondNode, ixFirstNode ) +=
1911  dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
1912  }
1913  else
1914  {
1915  smatMap( ixSecondNode, ixFirstNode ) +=
1916  dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
1917  }
1918  }
1919  else
1920  {
1921  ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
1922 
1923  if( !fNoConservation )
1924  {
1925  smatMap( ixSecondNode, ixFirstNode ) +=
1926  dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
1927  }
1928  else
1929  {
1930  smatMap( ixSecondNode, ixFirstNode ) +=
1931  dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
1932  }
1933  }
1934 
1935  ixs++;
1936  }
1937  }
1938 
1939  ixp++;
1940  }
1941  }
1942  }
1943 
1944  // Increment the current overlap index
1945  ixOverlap += nOverlapFaces;
1946  }
1947 
1948  return;
1949 }

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

◆ LinearRemapGLLtoGLL2_Pointwise_MOAB()

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

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

Definition at line 1953 of file TempestLinearRemap.cpp.

1963 {
1964  // Gauss-Lobatto quadrature within Faces
1965  DataArray1D< double > dGL;
1966  DataArray1D< double > dWL;
1967 
1968  GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
1969 
1970  // Get SparseMatrix represntation of the OfflineMap
1971  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1972 
1973  // Sample coefficients
1974  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1975 
1976  // Announcemnets
1977  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1978  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
1979  if( is_root )
1980  {
1981  dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
1982  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1983  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1984  }
1985 
1986  // Number of overlap Faces per source Face
1987  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1988 
1989  int ixOverlap = 0;
1990 
1991  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1992  {
1993  size_t ixOverlapTemp = ixOverlap;
1994  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1995  {
1996  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1997 
1998  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1999 
2000  nAllOverlapFaces[ixFirst]++;
2001  }
2002 
2003  // Increment the current overlap index
2004  ixOverlap += nAllOverlapFaces[ixFirst];
2005  }
2006 
2007  // Number of times this point was found
2008  DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
2009 
2010  ixOverlap = 0;
2011 #ifdef VERBOSE
2012  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
2013 #endif
2014  // Loop through all faces on m_meshInputCov
2015  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2016  {
2017 #ifdef VERBOSE
2018  // Announce computation progress
2019  if( ixFirst % outputFrequency == 0 && is_root )
2020  {
2021  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
2022  }
2023 #endif
2024  // Quantities from the First Mesh
2025  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
2026 
2027  const NodeVector& nodesFirst = m_meshInputCov->nodes;
2028 
2029  // Number of overlapping Faces and triangles
2030  int nOverlapFaces = nAllOverlapFaces[ixFirst];
2031 
2032  // Loop through all Overlap Faces
2033  for( int i = 0; i < nOverlapFaces; i++ )
2034  {
2035  // Quantities from the Second Mesh
2036  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
2037 
2038  // signal to not participate, because it is a ghost target
2039  if( ixSecond < 0 ) continue; // do not do anything
2040 
2041  const NodeVector& nodesSecond = m_meshOutput->nodes;
2042  const Face& faceSecond = m_meshOutput->faces[ixSecond];
2043 
2044  // Loop through all nodes on the second face
2045  for( int s = 0; s < nPout; s++ )
2046  {
2047  for( int t = 0; t < nPout; t++ )
2048  {
2049  size_t ixSecondNode;
2050  if( fContinuousOut )
2051  {
2052  ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
2053  }
2054  else
2055  {
2056  ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
2057  }
2058 
2059  if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
2060 
2061  // Check if this node has been found already
2062  if( fSecondNodeFound[ixSecondNode] ) continue;
2063 
2064  // Check this node
2065  Node node;
2066  Node dDx1G;
2067  Node dDx2G;
2068 
2069  ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
2070 
2071  // Find the components of this quadrature point in the basis
2072  // of the first Face.
2073  double dAlphaIn;
2074  double dBetaIn;
2075 
2076  ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
2077 
2078  // Check if this node is within the first Face
2079  if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
2080  ( dBetaIn > 1.0 + 1.0e-10 ) )
2081  continue;
2082 
2083  // Node is within the overlap region, mark as found
2084  fSecondNodeFound[ixSecondNode] = true;
2085 
2086  // Sample the First finite element at this point
2087  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
2088 
2089  // Add to map
2090  for( int p = 0; p < nPin; p++ )
2091  {
2092  for( int q = 0; q < nPin; q++ )
2093  {
2094  int ixFirstNode;
2095  if( fContinuousIn )
2096  {
2097  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
2098  }
2099  else
2100  {
2101  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
2102  }
2103 
2104  smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
2105  }
2106  }
2107  }
2108  }
2109  }
2110 
2111  // Increment the current overlap index
2112  ixOverlap += nOverlapFaces;
2113  }
2114 
2115  // Check for missing samples
2116  for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
2117  {
2118  if( !fSecondNodeFound[i] )
2119  {
2120  _EXCEPTION1( "Can't sample point %i", i );
2121  }
2122  }
2123 
2124  return;
2125 }

References dbgprint.

◆ LinearRemapNN_MOAB()

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

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

Definition at line 58 of file TempestLinearRemap.cpp.

59 {
60  /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov) */
61 
62 #ifdef VVERBOSE
63  {
64  std::ofstream output_file( "rowcolindices.txt", std::ios::out );
65  output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
66  << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
67  output_file << "Rows \n";
68  for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
69  output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
70  output_file << "Cols \n";
71  for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
72  output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
73  output_file.flush(); // required here
74  output_file.close();
75  }
76 #endif
77 
78  if( use_GID_matching )
79  {
80  std::map< unsigned, unsigned > src_gl;
81  for( unsigned it = 0; it < col_gdofmap.size(); ++it )
82  src_gl[col_gdofmap[it]] = it;
83 
84  std::map< unsigned, unsigned >::iterator iter;
85  for( unsigned it = 0; it < row_gdofmap.size(); ++it )
86  {
87  unsigned row = row_gdofmap[it];
88  iter = src_gl.find( row );
89  if( strict_check && iter == src_gl.end() )
90  {
91  std::cout << "Searching for global target DOF " << row
92  << " but could not find correspondence in source mesh.\n";
93  assert( false );
94  }
95  else if( iter == src_gl.end() )
96  {
97  continue;
98  }
99  else
100  {
101  unsigned icol = src_gl[row];
102  unsigned irow = it;
103 
104  // Set the permutation matrix in local space
105  m_mapRemap( irow, icol ) = 1.0;
106  }
107  }
108 
109  return moab::MB_SUCCESS;
110  }
111  else
112  {
113  /* Create a Kd-tree to perform local queries to find nearest neighbors */
114 
115  return moab::MB_FAILURE;
116  }
117 }

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

◆ LinearRemapSE0_Tempest_MOAB()

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

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

◆ LinearRemapSE4_Tempest_MOAB()

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

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

Definition at line 978 of file TempestLinearRemap.cpp.

983 {
984  // Order of the polynomial interpolant
985  int nP = dataGLLNodes.GetRows();
986 
987  // Order of triangular quadrature rule
988  const int TriQuadRuleOrder = 4;
989 
990  // Triangular quadrature rule
991  TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
992 
993  int TriQuadraturePoints = triquadrule.GetPoints();
994 
995  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
996 
997  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
998 
999  // Sample coefficients
1000  DataArray2D< double > dSampleCoeff( nP, nP );
1001 
1002  // GLL Quadrature nodes on quadrilateral elements
1003  DataArray1D< double > dG;
1004  DataArray1D< double > dW;
1005  GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
1006 
1007  // Announcements
1008  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1009  dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
1010  if( is_root )
1011  {
1012  dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
1013  dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
1014  dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
1015  }
1016 
1017  // Get SparseMatrix represntation of the OfflineMap
1018  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1019 
1020  // NodeVector from m_meshOverlap
1021  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1022  const NodeVector& nodesFirst = m_meshInputCov->nodes;
1023 
1024  // Vector of source areas
1025  DataArray1D< double > vecSourceArea( nP * nP );
1026 
1027  DataArray1D< double > vecTargetArea;
1028  DataArray2D< double > dCoeff;
1029 
1030 #ifdef VERBOSE
1031  std::stringstream sstr;
1032  sstr << "remapdata_" << rank << ".txt";
1033  std::ofstream output_file( sstr.str() );
1034 #endif
1035 
1036  // Current Overlap Face
1037  int ixOverlap = 0;
1038 #ifdef VERBOSE
1039  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1040 #endif
1041  // generic triangle used for area computation, for triangles around the center of overlap face;
1042  // used for overlap faces with more than 4 edges;
1043  // nodes array will be set for each triangle;
1044  // these triangles are not part of the mesh structure, they are just temporary during
1045  // aforementioned decomposition.
1046  Face faceTri( 3 );
1047  NodeVector nodes( 3 );
1048  faceTri.SetNode( 0, 0 );
1049  faceTri.SetNode( 1, 1 );
1050  faceTri.SetNode( 2, 2 );
1051 
1052  // Loop over all input Faces
1053  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1054  {
1055  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1056 
1057  if( faceFirst.edges.size() != 4 )
1058  {
1059  _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
1060  }
1061 #ifdef VERBOSE
1062  // Announce computation progress
1063  if( ixFirst % outputFrequency == 0 && is_root )
1064  {
1065  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1066  }
1067 #endif
1068  // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
1069  // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
1070  // However, the relation with MOAB and Tempest will go out of the roof
1071 
1072  // Determine how many overlap Faces and triangles are present
1073  int nOverlapFaces = 0;
1074  size_t ixOverlapTemp = ixOverlap;
1075  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1076  {
1077  // if( m_meshOverlap->vecTargetFaceIx[ixOverlapTemp] < 0 ) continue; // skip ghost target faces
1078  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1079  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1080 
1081  nOverlapFaces++;
1082  }
1083 
1084  // No overlaps
1085  if( nOverlapFaces == 0 )
1086  {
1087  continue;
1088  }
1089 
1090  // Allocate remap coefficients array for meshFirst Face
1091  DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
1092 
1093  // Find the local remap coefficients
1094  for( int j = 0; j < nOverlapFaces; j++ )
1095  {
1096  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
1097  if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
1098  {
1099  Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
1100  m_meshOverlap->vecFaceArea[ixOverlap + j] );
1101  int n = faceOverlap.edges.size();
1102  Announce( "Number nodes: %d", n );
1103  for( int k = 0; k < n; k++ )
1104  {
1105  Node nd = nodesOverlap[faceOverlap[k]];
1106  Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
1107  }
1108  continue;
1109  }
1110 
1111  // #ifdef VERBOSE
1112  // if ( is_root )
1113  // Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
1114  // m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
1115  // m_meshOverlap->vecFaceArea[ixOverlap + j] );
1116  // #endif
1117 
1118  int nbEdges = faceOverlap.edges.size();
1119  int nOverlapTriangles = 1;
1120  Node center; // not used if nbEdges == 3
1121  if( nbEdges > 3 )
1122  { // decompose from center in this case
1123  nOverlapTriangles = nbEdges;
1124  for( int k = 0; k < nbEdges; k++ )
1125  {
1126  const Node& node = nodesOverlap[faceOverlap[k]];
1127  center = center + node;
1128  }
1129  center = center / nbEdges;
1130  center = center.Normalized(); // project back on sphere of radius 1
1131  }
1132 
1133  Node node0, node1, node2;
1134  double dTriangleArea;
1135 
1136  // Loop over all sub-triangles of this Overlap Face
1137  for( int k = 0; k < nOverlapTriangles; k++ )
1138  {
1139  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1140  {
1141  node0 = nodesOverlap[faceOverlap[0]];
1142  node1 = nodesOverlap[faceOverlap[1]];
1143  node2 = nodesOverlap[faceOverlap[2]];
1144  dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1145  }
1146  else // decompose polygon in triangles around the center
1147  {
1148  node0 = center;
1149  node1 = nodesOverlap[faceOverlap[k]];
1150  int k1 = ( k + 1 ) % nbEdges;
1151  node2 = nodesOverlap[faceOverlap[k1]];
1152  nodes[0] = center;
1153  nodes[1] = node1;
1154  nodes[2] = node2;
1155  dTriangleArea = CalculateFaceArea( faceTri, nodes );
1156  }
1157  // Coordinates of quadrature Node
1158  for( int l = 0; l < TriQuadraturePoints; l++ )
1159  {
1160  Node nodeQuadrature;
1161  nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
1162  TriQuadratureG[l][2] * node2.x;
1163 
1164  nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
1165  TriQuadratureG[l][2] * node2.y;
1166 
1167  nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
1168  TriQuadratureG[l][2] * node2.z;
1169 
1170  nodeQuadrature = nodeQuadrature.Normalized();
1171 
1172  // Find components of quadrature point in basis
1173  // of the first Face
1174  double dAlpha;
1175  double dBeta;
1176 
1177  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
1178 
1179  // Check inverse map value
1180  if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
1181  ( dBeta > 1.0 + 1.0e-13 ) )
1182  {
1183  _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
1184  "(%1.5e %1.5e)",
1185  j, l, dAlpha, dBeta );
1186  }
1187 
1188  // Sample the finite element at this point
1189  SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
1190 
1191  // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
1192  for( int p = 0; p < nP; p++ )
1193  {
1194  for( int q = 0; q < nP; q++ )
1195  {
1196  dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
1197  m_meshOverlap->vecFaceArea[ixOverlap + j];
1198  }
1199  }
1200  }
1201  }
1202  }
1203 
1204 #ifdef VERBOSE
1205  output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1206  for( int j = 0; j < nOverlapFaces; j++ )
1207  {
1208  for( int p = 0; p < nP; p++ )
1209  {
1210  for( int q = 0; q < nP; q++ )
1211  {
1212  output_file << dRemapCoeff[p][q][j] << " ";
1213  }
1214  }
1215  }
1216  output_file << std::endl;
1217 #endif
1218 
1219  // Force consistency and conservation
1220  if( !fNoConservation )
1221  {
1222  double dTargetArea = 0.0;
1223  for( int j = 0; j < nOverlapFaces; j++ )
1224  {
1225  dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
1226  }
1227 
1228  for( int p = 0; p < nP; p++ )
1229  {
1230  for( int q = 0; q < nP; q++ )
1231  {
1232  vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1233  }
1234  }
1235 
1236  const double areaTolerance = 1e-10;
1237  // Source elements are completely covered by target volumes
1238  if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
1239  {
1240  vecTargetArea.Allocate( nOverlapFaces );
1241  for( int j = 0; j < nOverlapFaces; j++ )
1242  {
1243  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1244  }
1245 
1246  dCoeff.Allocate( nOverlapFaces, nP * nP );
1247 
1248  for( int j = 0; j < nOverlapFaces; j++ )
1249  {
1250  for( int p = 0; p < nP; p++ )
1251  {
1252  for( int q = 0; q < nP; q++ )
1253  {
1254  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1255  }
1256  }
1257  }
1258 
1259  // Target volumes only partially cover source elements
1260  }
1261  else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
1262  {
1263  double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
1264 
1265  vecTargetArea.Allocate( nOverlapFaces + 1 );
1266  for( int j = 0; j < nOverlapFaces; j++ )
1267  {
1268  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1269  }
1270  vecTargetArea[nOverlapFaces] = dExtraneousArea;
1271 
1272 #ifdef VERBOSE
1273  Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
1274  m_meshInputCov->vecFaceArea[ixFirst] );
1275 #endif
1276  if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
1277  {
1278  _EXCEPTIONT( "Partial element area exceeds total element area" );
1279  }
1280 
1281  dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
1282 
1283  for( int j = 0; j < nOverlapFaces; j++ )
1284  {
1285  for( int p = 0; p < nP; p++ )
1286  {
1287  for( int q = 0; q < nP; q++ )
1288  {
1289  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1290  }
1291  }
1292  }
1293  for( int p = 0; p < nP; p++ )
1294  {
1295  for( int q = 0; q < nP; q++ )
1296  {
1297  dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1298  }
1299  }
1300  for( int j = 0; j < nOverlapFaces; j++ )
1301  {
1302  for( int p = 0; p < nP; p++ )
1303  {
1304  for( int q = 0; q < nP; q++ )
1305  {
1306  dCoeff[nOverlapFaces][p * nP + q] -=
1307  dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
1308  }
1309  }
1310  }
1311  for( int p = 0; p < nP; p++ )
1312  {
1313  for( int q = 0; q < nP; q++ )
1314  {
1315  dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
1316  }
1317  }
1318 
1319  // Source elements only partially cover target volumes
1320  }
1321  else
1322  {
1323  Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
1324  m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
1325  _EXCEPTIONT( "Target grid must be a subset of source grid" );
1326  }
1327 
1328  ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
1329  /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
1330 
1331  for( int j = 0; j < nOverlapFaces; j++ )
1332  {
1333  for( int p = 0; p < nP; p++ )
1334  {
1335  for( int q = 0; q < nP; q++ )
1336  {
1337  dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
1338  }
1339  }
1340  }
1341  }
1342 
1343 #ifdef VERBOSE
1344  // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1345  // for ( int j = 0; j < nOverlapFaces; j++ )
1346  // {
1347  // for ( int p = 0; p < nP; p++ )
1348  // {
1349  // for ( int q = 0; q < nP; q++ )
1350  // {
1351  // output_file << dRemapCoeff[p][q][j] << " ";
1352  // }
1353  // }
1354  // }
1355  // output_file << std::endl;
1356 #endif
1357 
1358  // Put these remap coefficients into the SparseMatrix map
1359  for( int j = 0; j < nOverlapFaces; j++ )
1360  {
1361  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1362 
1363  // signal to not participate, because it is a ghost target
1364  if( ixSecondFace < 0 ) continue; // do not do anything
1365 
1366  for( int p = 0; p < nP; p++ )
1367  {
1368  for( int q = 0; q < nP; q++ )
1369  {
1370  if( fContinuousIn )
1371  {
1372  int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
1373 
1374  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1375  m_meshOverlap->vecFaceArea[ixOverlap + j] /
1376  m_meshOutput->vecFaceArea[ixSecondFace];
1377  }
1378  else
1379  {
1380  int ixFirstNode = ixFirst * nP * nP + p * nP + q;
1381 
1382  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1383  m_meshOverlap->vecFaceArea[ixOverlap + j] /
1384  m_meshOutput->vecFaceArea[ixSecondFace];
1385  }
1386  }
1387  }
1388  }
1389  // Increment the current overlap index
1390  ixOverlap += nOverlapFaces;
1391  }
1392 #ifdef VERBOSE
1393  output_file.flush(); // required here
1394  output_file.close();
1395 #endif
1396 
1397  return;
1398 }

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

◆ PrintMapStatistics()

void moab::TempestOnlineMap::PrintMapStatistics ( )

Definition at line 284 of file TempestLinearRemap.cpp.

285 {
286  int nrows = m_weightMatrix.rows(); // Number of rows
287  int ncols = m_weightMatrix.cols(); // Number of columns
288  int NNZ = m_weightMatrix.nonZeros(); // Number of non zero values
289 #ifdef MOAB_HAVE_MPI
290  // find out min/max for NNZ, ncols, nrows
291  // should work on std c++ 11
292  int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
293  int rarr3[6];
294  MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
295 
296  int total[3];
297  MPI_Reduce( arr3, total, 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
298  if( !rank )
299  std::cout << "-> Rows (min/max/sum): (" << rarr3[1] << " / " << -rarr3[4] << " / " << total[1] << "), "
300  << " Cols (min/max/sum): (" << rarr3[2] << " / " << -rarr3[5] << " / " << total[2] << "), "
301  << " NNZ (min/max/sum): (" << rarr3[0] << " / " << -rarr3[3] << " / " << total[0] << ")\n";
302 #else
303  std::cout << "-> Rows: " << nrows << ", Cols: " << ncols << ", NNZ: " << NNZ << "\n";
304 #endif
305 }

Referenced by main().

◆ QLTLimiter()

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

◆ ReadParallelMap()

moab::ErrorCode moab::TempestOnlineMap::ReadParallelMap ( const char *  strSource,
const std::vector< int > &  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 1165 of file TempestOnlineMapIO.cpp.

1168 {
1169  NcError error( NcError::silent_nonfatal );
1170 
1171  NcVar *varRow = NULL, *varCol = NULL, *varS = NULL;
1172  int nS = 0, nA = 0, nB = 0;
1173 #ifdef MOAB_HAVE_PNETCDF
1174  // some variables will be used just in the case netcdfpar reader fails
1175  int ncfile = -1, ret = 0;
1176  int ndims, nvars, ngatts, unlimited;
1177 #endif
1178 #ifdef MOAB_HAVE_NETCDFPAR
1179  bool is_independent = true;
1180  ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1181  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
1182 #else
1183  NcFile ncMap( strSource, NcFile::ReadOnly );
1184 #endif
1185 
1186 #define CHECK_EXCEPTION( obj, type, varstr ) \
1187  { \
1188  if( obj == NULL ) \
1189  { \
1190  _EXCEPTION3( "Map file \"%s\" does not contain %s \"%s\"", strSource, type, varstr ); \
1191  } \
1192  }
1193 
1194  // Read SparseMatrix entries
1195 
1196  if( ncMap.is_valid() )
1197  {
1198  NcDim* dimNS = ncMap.get_dim( "n_s" );
1199  CHECK_EXCEPTION( dimNS, "dimension", "n_s" );
1200 
1201  NcDim* dimNA = ncMap.get_dim( "n_a" );
1202  CHECK_EXCEPTION( dimNA, "dimension", "n_a" );
1203 
1204  NcDim* dimNB = ncMap.get_dim( "n_b" );
1205  CHECK_EXCEPTION( dimNB, "dimension", "n_b" );
1206 
1207  // store total number of nonzeros
1208  nS = dimNS->size();
1209  nA = dimNA->size();
1210  nB = dimNB->size();
1211 
1212  varRow = ncMap.get_var( "row" );
1213  CHECK_EXCEPTION( varRow, "variable", "row" );
1214 
1215  varCol = ncMap.get_var( "col" );
1216  CHECK_EXCEPTION( varCol, "variable", "col" );
1217 
1218  varS = ncMap.get_var( "S" );
1219  CHECK_EXCEPTION( varS, "variable", "S" );
1220 
1221 #ifdef MOAB_HAVE_NETCDFPAR
1222  ncMap.enable_var_par_access( varRow, is_independent );
1223  ncMap.enable_var_par_access( varCol, is_independent );
1224  ncMap.enable_var_par_access( varS, is_independent );
1225 #endif
1226  }
1227  else
1228  {
1229 #ifdef MOAB_HAVE_PNETCDF
1230  // read the file using pnetcdf directly, in parallel; need to have MPI, we do not check that anymore
1231  // why build wth pnetcdf without MPI ?
1232  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strSource, NcFile::ReadOnly, NcFile::Netcdf4 );
1233  ret = ncmpi_open( m_pcomm->comm(), strSource, NC_NOWRITE, MPI_INFO_NULL, &ncfile );
1234  ERR_PARNC( ret ); // bail out completely
1235  ret = ncmpi_inq( ncfile, &ndims, &nvars, &ngatts, &unlimited );
1236  ERR_PARNC( ret );
1237  // find dimension ids for n_S
1238  int ins;
1239  ret = ncmpi_inq_dimid( ncfile, "n_s", &ins );
1240  ERR_PARNC( ret );
1241  MPI_Offset leng;
1242  ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
1243  ERR_PARNC( ret );
1244  nS = (int)leng;
1245  ret = ncmpi_inq_dimid( ncfile, "n_b", &ins );
1246  ERR_PARNC( ret );
1247  ret = ncmpi_inq_dimlen( ncfile, ins, &leng );
1248  ERR_PARNC( ret );
1249  nB = (int)leng;
1250 #else
1251  _EXCEPTION1( "cannot read the file %s", strSource );
1252 #endif
1253  }
1254 
1255  // Let us declare the map object for every process
1256  SparseMatrix< double >& sparseMatrix = this->GetSparseMatrix();
1257 
1258  int localSize = nS / size;
1259  long offsetRead = rank * localSize;
1260  // leftovers on last rank
1261  if( rank == size - 1 )
1262  {
1263  localSize += nS % size;
1264  }
1265 
1266  std::vector< int > vecRow, vecCol;
1267  std::vector< double > vecS;
1268  vecRow.resize( localSize );
1269  vecCol.resize( localSize );
1270  vecS.resize( localSize );
1271 
1272  if( ncMap.is_valid() )
1273  {
1274  varRow->set_cur( (long)( offsetRead ) );
1275  varRow->get( &( vecRow[0] ), localSize );
1276 
1277  varCol->set_cur( (long)( offsetRead ) );
1278  varCol->get( &( vecCol[0] ), localSize );
1279 
1280  varS->set_cur( (long)( offsetRead ) );
1281  varS->get( &( vecS[0] ), localSize );
1282 
1283  ncMap.close();
1284  }
1285  else
1286  {
1287 #ifdef MOAB_HAVE_PNETCDF
1288  // fill the local vectors with the variables from pnetcdf file; first inquire, then fill
1289  MPI_Offset start = (MPI_Offset)offsetRead;
1290  MPI_Offset count = (MPI_Offset)localSize;
1291  int varid;
1292  ret = ncmpi_inq_varid( ncfile, "S", &varid );
1293  ERR_PARNC( ret );
1294  ret = ncmpi_get_vara_double_all( ncfile, varid, &start, &count, &vecS[0] );
1295  ERR_PARNC( ret );
1296  ret = ncmpi_inq_varid( ncfile, "row", &varid );
1297  ERR_PARNC( ret );
1298  ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecRow[0] );
1299  ERR_PARNC( ret );
1300  ret = ncmpi_inq_varid( ncfile, "col", &varid );
1301  ERR_PARNC( ret );
1302  ret = ncmpi_get_vara_int_all( ncfile, varid, &start, &count, &vecCol[0] );
1303  ERR_PARNC( ret );
1304  ret = ncmpi_close( ncfile );
1305  ERR_PARNC( ret );
1306 #endif
1307  }
1308 
1309  // Now let us set the necessary global-to-local ID maps so that A*x operations
1310  // can be performed cleanly as if map was computed online
1311  row_dtoc_dofmap.clear();
1312  // row_dtoc_dofmap.reserve( nB / size );
1313  col_dtoc_dofmap.clear();
1314  rowMap.clear();
1315  colMap.clear();
1316  // col_dtoc_dofmap.reserve( 2 * nA / size );
1317  // row_dtoc_dofmap.resize( m_nTotDofs_Dest, UINT_MAX );
1318  // col_dtoc_dofmap.resize( m_nTotDofs_SrcCov, UINT_MAX );
1319 
1320 #ifdef MOAB_HAVE_MPI
1321  // bother with tuple list only if size > 1
1322  // otherwise, just fill the sparse matrix
1323  if( size > 1 )
1324  {
1325  std::vector< int > ownership;
1326  // the default trivial partitioning scheme
1327  int nDofs = nB; // this is for row partitioning
1328  if( !row_partition ) nDofs = nA; // column partitioning
1329 
1330  // assert(row_major_ownership == true); // this block is valid only for row-based partitioning
1331  ownership.resize( size );
1332  int nPerPart = nDofs / size;
1333  int nRemainder = nDofs % size; // Keep the remainder in root
1334  ownership[0] = nPerPart + nRemainder;
1335  for( int ip = 1, roffset = ownership[0]; ip < size; ++ip )
1336  {
1337  roffset += nPerPart;
1338  ownership[ip] = roffset;
1339  }
1340  moab::TupleList* tl = new moab::TupleList;
1341  unsigned numr = 1; //
1342  tl->initialize( 3, 0, 0, numr, localSize ); // to proc, row, col, value
1343  tl->enableWriteAccess();
1344  // populate
1345  for( int i = 0; i < localSize; i++ )
1346  {
1347  int rowval = vecRow[i] - 1; // dofs are 1 based in the file
1348  int colval = vecCol[i] - 1;
1349  int to_proc = -1;
1350  int dof_val = colval;
1351  if( row_partition ) dof_val = rowval;
1352 
1353  if( ownership[0] > dof_val )
1354  to_proc = 0;
1355  else
1356  {
1357  for( int ip = 1; ip < size; ++ip )
1358  {
1359  if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1360  {
1361  to_proc = ip;
1362  break;
1363  }
1364  }
1365  }
1366 
1367  int n = tl->get_n();
1368  tl->vi_wr[3 * n] = to_proc;
1369  tl->vi_wr[3 * n + 1] = rowval;
1370  tl->vi_wr[3 * n + 2] = colval;
1371  tl->vr_wr[n] = vecS[i];
1372  tl->inc_n();
1373  }
1374  // heavy communication
1375  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl, 0 );
1376 
1377  if( owned_dof_ids.size() > 0 )
1378  {
1379  // we need to send desired dof to the rendez_vous point
1380  moab::TupleList tl_re; //
1381  tl_re.initialize( 2, 0, 0, 0, owned_dof_ids.size() ); // to proc, value
1382  tl_re.enableWriteAccess();
1383  // send first to rendez_vous point, decided by trivial partitioning
1384 
1385  for( size_t i = 0; i < owned_dof_ids.size(); i++ )
1386  {
1387  int to_proc = -1;
1388  int dof_val = owned_dof_ids[i] - 1; // dofs are 1 based in the file, partition from 0 ?
1389 
1390  if( ownership[0] > dof_val )
1391  to_proc = 0;
1392  else
1393  {
1394  for( int ip = 1; ip < size; ++ip )
1395  {
1396  if( ownership[ip - 1] <= dof_val && ownership[ip] > dof_val )
1397  {
1398  to_proc = ip;
1399  break;
1400  }
1401  }
1402  }
1403 
1404  int n = tl_re.get_n();
1405  tl_re.vi_wr[2 * n] = to_proc;
1406  tl_re.vi_wr[2 * n + 1] = dof_val;
1407 
1408  tl_re.inc_n();
1409  }
1410  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tl_re, 0 );
1411  // now we know in tl_re where do we need to send back dof_val
1412  moab::TupleList::buffer sort_buffer;
1413  sort_buffer.buffer_init( tl_re.get_n() );
1414  tl_re.sort( 1, &sort_buffer ); // so now we order by value
1415 
1416  sort_buffer.buffer_init( tl->get_n() );
1417  int indexOrder = 2; // colVal
1418  if( row_partition ) indexOrder = 1; // rowVal
1419  //tl->sort( indexOrder, &sort_buffer );
1420 
1421  std::map< int, int > startDofIndex, endDofIndex; // indices in tl_re for values we want
1422  int dofVal = -1;
1423  if( tl_re.get_n() > 0 ) dofVal = tl_re.vi_rd[1]; // first dof val on this rank
1424  startDofIndex[dofVal] = 0;
1425  endDofIndex[dofVal] = 0; // start and end
1426  for( unsigned k = 1; k < tl_re.get_n(); k++ )
1427  {
1428  int newDof = tl_re.vi_rd[2 * k + 1];
1429  if( dofVal == newDof )
1430  {
1431  endDofIndex[dofVal] = k; // increment by 1 actually
1432  }
1433  else
1434  {
1435  dofVal = newDof;
1436  startDofIndex[dofVal] = k;
1437  endDofIndex[dofVal] = k;
1438  }
1439  }
1440 
1441  // basically, for each value we are interested in, index in tl_re with those values are
1442  // tl_re.vi_rd[2*startDofIndex+1] == valDof == tl_re.vi_rd[2*endDofIndex+1]
1443  // so now we have ordered
1444  // tl_re shows to what proc do we need to send the tuple (row, col, val)
1445  moab::TupleList* tl_back = new moab::TupleList;
1446  unsigned numr = 1; //
1447  // localSize is a good guess, but maybe it should be bigger ?
1448  // this could be bigger for repeated dofs
1449  tl_back->initialize( 3, 0, 0, numr, tl->get_n() ); // to proc, row, col, value
1450  tl_back->enableWriteAccess();
1451  // now loop over tl and tl_re to see where to send
1452  // form the new tuple, which will contain the desired dofs per task, per row or column distribution
1453 
1454  for( unsigned k = 0; k < tl->get_n(); k++ )
1455  {
1456  int valDof = tl->vi_rd[3 * k + indexOrder]; // 1 for row, 2 for column // first value, it should be
1457  for( int ire = startDofIndex[valDof]; ire <= endDofIndex[valDof]; ire++ )
1458  {
1459  int to_proc = tl_re.vi_rd[2 * ire];
1460  int n = tl_back->get_n();
1461  tl_back->vi_wr[3 * n] = to_proc;
1462  tl_back->vi_wr[3 * n + 1] = tl->vi_rd[3 * k + 1]; // row
1463  tl_back->vi_wr[3 * n + 2] = tl->vi_rd[3 * k + 2]; // col
1464  tl_back->vr_wr[n] = tl->vr_rd[k];
1465  tl_back->inc_n();
1466  }
1467  }
1468 
1469  // now communicate to the desired tasks:
1470  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, *tl_back, 0 );
1471 
1472  tl_re.reset(); // clear memory, although this will go out of scope
1473  tl->reset();
1474  tl = tl_back;
1475  }
1476 
1477  int rindexMax = 0, cindexMax = 0;
1478  // populate the sparsematrix, using rowMap and colMap
1479  int n = tl->get_n();
1480  for( int i = 0; i < n; i++ )
1481  {
1482  int rindex, cindex;
1483  const int& vecRowValue = tl->vi_wr[3 * i + 1];
1484  const int& vecColValue = tl->vi_wr[3 * i + 2];
1485 
1486  std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
1487  if( riter == rowMap.end() )
1488  {
1489  rowMap[vecRowValue] = rindexMax;
1490  rindex = rindexMax;
1491  row_gdofmap.push_back( vecRowValue );
1492  // row_dtoc_dofmap.push_back( vecRowValue );
1493  rindexMax++;
1494  }
1495  else
1496  rindex = riter->second;
1497 
1498  std::map< int, int >::iterator citer = colMap.find( vecColValue );
1499  if( citer == colMap.end() )
1500  {
1501  colMap[vecColValue] = cindexMax;
1502  cindex = cindexMax;
1503  col_gdofmap.push_back( vecColValue );
1504  // col_dtoc_dofmap.push_back( vecColValue );
1505  cindexMax++;
1506  }
1507  else
1508  cindex = citer->second;
1509 
1510  sparseMatrix( rindex, cindex ) = tl->vr_wr[i];
1511  }
1512  tl->reset();
1513  }
1514  else
1515 #endif
1516  {
1517  int rindexMax = 0, cindexMax = 0;
1518 
1519  for( int i = 0; i < nS; i++ )
1520  {
1521  int rindex, cindex;
1522  const int& vecRowValue = vecRow[i] - 1; // the rows, cols are 1 based in the file
1523  const int& vecColValue = vecCol[i] - 1;
1524 
1525  std::map< int, int >::iterator riter = rowMap.find( vecRowValue );
1526  if( riter == rowMap.end() )
1527  {
1528  rowMap[vecRowValue] = rindexMax;
1529  rindex = rindexMax;
1530  row_gdofmap.push_back( vecRowValue );
1531  // row_dtoc_dofmap.push_back( vecRowValue );
1532  rindexMax++;
1533  }
1534  else
1535  rindex = riter->second;
1536 
1537  std::map< int, int >::iterator citer = colMap.find( vecColValue );
1538  if( citer == colMap.end() )
1539  {
1540  colMap[vecColValue] = cindexMax;
1541  cindex = cindexMax;
1542  col_gdofmap.push_back( vecColValue );
1543  // col_dtoc_dofmap.push_back( vecColValue );
1544  cindexMax++;
1545  }
1546  else
1547  cindex = citer->second;
1548 
1549  sparseMatrix( rindex, cindex ) = vecS[i];
1550  }
1551  }
1552 
1553  m_nTotDofs_SrcCov = sparseMatrix.GetColumns();
1554  m_nTotDofs_Dest = sparseMatrix.GetRows();
1555 
1556 #ifdef MOAB_HAVE_EIGEN3
1557  this->copy_tempest_sparsemat_to_eigen3();
1558 #endif
1559 
1560  // Reset the source and target data first
1561  m_rowVector.setZero();
1562  m_colVector.setZero();
1563 
1564  return moab::MB_SUCCESS;
1565 }

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

◆ set_col_dc_dofs()

moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs ( std::vector< int > &  values_entities)

Definition at line 638 of file TempestOnlineMap.cpp.

639 {
640  // col_gdofmap has global dofs , that should be in the list of values, such that
641  // row_dtoc_dofmap[offsetDOF] = localDOF;
642  // // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
643  // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
644  // form first inverse
645 
646  col_dtoc_dofmap.resize( values_entities.size() );
647  for( int j = 0; j < (int)values_entities.size(); j++ )
648  {
649  if( colMap.find( values_entities[j] - 1 ) != colMap.end() )
650  col_dtoc_dofmap[j] = colMap[values_entities[j] - 1];
651  else
652  {
653  col_dtoc_dofmap[j] = -1; // signal that this value should not be used in
654  // std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
655  // found in colMap \n";
656  }
657  }
658  return moab::MB_SUCCESS;
659 }

References MB_SUCCESS.

◆ set_row_dc_dofs()

moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs ( std::vector< int > &  values_entities)

Definition at line 661 of file TempestOnlineMap.cpp.

662 {
663  // row_dtoc_dofmap = values_entities; // needs to point to local
664  // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
665 
666  row_dtoc_dofmap.resize( values_entities.size() );
667  for( int j = 0; j < (int)values_entities.size(); j++ )
668  {
669  if( rowMap.find( values_entities[j] - 1 ) != rowMap.end() )
670  row_dtoc_dofmap[j] = rowMap[values_entities[j] - 1]; // values are 1 based, but rowMap, colMap are not
671  else
672  {
673  row_dtoc_dofmap[j] = -1; // not all values are used
674  // std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
675  // found in rowMap \n";
676  }
677  }
678  return moab::MB_SUCCESS;
679 }

References MB_SUCCESS.

◆ SetDestinationNDofsPerElement()

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

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

Definition at line 586 of file TempestOnlineMap.hpp.

587 {
588  m_nDofsPEl_Dest = nt;
589 }

◆ SetDOFmapAssociation()

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

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

Parameters
srcTypeThe discretization type of the source mesh
srcOrderThe order of the discretization on the source mesh
isSrcContinuousThe continuity of the discretization on the source mesh
srcdataGLLNodesThe GLL nodes on the source mesh
srcdataGLLNodesSrcThe GLL nodes on the source mesh
destTypeThe discretization type of the destination mesh
destOrderThe order of the discretization on the destination mesh
isTgtContinuousThe continuity of the discretization on the destination mesh
tgtdataGLLNodesThe GLL nodes on the destination mesh

Definition at line 160 of file TempestOnlineMap.cpp.

169 {
170  moab::ErrorCode rval;
171  std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
172  std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
173 
174  // We are assuming that these are element based tags that are sized: np * np
175  m_srcDiscType = srcType;
176  m_destDiscType = destType;
177  m_input_order = srcOrder;
178  m_output_order = destOrder;
179 
180  bool vprint = is_root && false;
181 
182  // Compute and store the total number of source and target DoFs corresponding
183  // to number of rows and columns in the mapping.
184 #ifdef VVERBOSE
185  {
186  src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
187  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
188  locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
189  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
190  tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
191  rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
192 
193  if( is_root )
194  {
195  {
196  std::ofstream output_file( "sourcecov-gids-0.txt" );
197  output_file << "I, GDOF\n";
198  for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
199  output_file << i << ", " << src_soln_gdofs[i] << "\n";
200 
201  output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
202  m_nTotDofs_SrcCov = 0;
203  if( isSrcContinuous )
204  dgll_cgll_covcol_ldofmap.resize(
206  for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
207  {
208  for( int p = 0; p < m_nDofsPEl_Src; p++ )
209  {
210  for( int q = 0; q < m_nDofsPEl_Src; q++ )
211  {
212  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
213  const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
214  if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
215  {
217  dgll_cgll_covcol_ldofmap[localDOF] = true;
218  }
219  output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
220  << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
221  }
222  }
223  }
224  output_file.flush(); // required here
225  output_file.close();
226  dgll_cgll_covcol_ldofmap.clear();
227  }
228 
229  {
230  std::ofstream output_file( "source-gids-0.txt" );
231  output_file << "I, GDOF\n";
232  for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
233  output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
234 
235  output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
236  m_nTotDofs_Src = 0;
237  if( isSrcContinuous )
238  dgll_cgll_col_ldofmap.resize(
240  for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
241  {
242  for( int p = 0; p < m_nDofsPEl_Src; p++ )
243  {
244  for( int q = 0; q < m_nDofsPEl_Src; q++ )
245  {
246  const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
247  const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
248  if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
249  {
250  m_nTotDofs_Src++;
251  dgll_cgll_col_ldofmap[localDOF] = true;
252  }
253  output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
254  << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
255  }
256  }
257  }
258  output_file.flush(); // required here
259  output_file.close();
260  dgll_cgll_col_ldofmap.clear();
261  }
262 
263  {
264  std::ofstream output_file( "target-gids-0.txt" );
265  output_file << "I, GDOF\n";
266  for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
267  output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
268 
269  output_file << "ELEMID, IDOF, GDOF, NDOF\n";
270  m_nTotDofs_Dest = 0;
271 
272  for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
273  {
274  output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
275  << m_nTotDofs_Dest << "\n";
276  m_nTotDofs_Dest++;
277  }
278 
279  output_file.flush(); // required here
280  output_file.close();
281  }
282  }
283  else
284  {
285  {
286  std::ofstream output_file( "sourcecov-gids-1.txt" );
287  output_file << "I, GDOF\n";
288  for( unsigned i = 0; i < src_soln_gdofs.size(); ++i )
289  output_file << i << ", " << src_soln_gdofs[i] << "\n";
290 
291  output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
292  m_nTotDofs_SrcCov = 0;
293  if( isSrcContinuous )
294  dgll_cgll_covcol_ldofmap.resize(
296  for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
297  {
298  for( int p = 0; p < m_nDofsPEl_Src; p++ )
299  {
300  for( int q = 0; q < m_nDofsPEl_Src; q++ )
301  {
302  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
303  const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
304  if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
305  {
307  dgll_cgll_covcol_ldofmap[localDOF] = true;
308  }
309  output_file << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", " << localDOF
310  << ", " << src_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
311  }
312  }
313  }
314  output_file.flush(); // required here
315  output_file.close();
316  dgll_cgll_covcol_ldofmap.clear();
317  }
318 
319  {
320  std::ofstream output_file( "source-gids-1.txt" );
321  output_file << "I, GDOF\n";
322  for( unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
323  output_file << i << ", " << locsrc_soln_gdofs[i] << "\n";
324 
325  output_file << "ELEMID, IDOF, LDOF, GDOF, NDOF\n";
326  m_nTotDofs_Src = 0;
327  if( isSrcContinuous )
328  dgll_cgll_col_ldofmap.resize(
330  for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
331  {
332  for( int p = 0; p < m_nDofsPEl_Src; p++ )
333  {
334  for( int q = 0; q < m_nDofsPEl_Src; q++ )
335  {
336  const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
337  const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
338  if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
339  {
340  m_nTotDofs_Src++;
341  dgll_cgll_col_ldofmap[localDOF] = true;
342  }
343  output_file << m_remapper->lid_to_gid_src[j] << ", " << offsetDOF << ", " << localDOF
344  << ", " << locsrc_soln_gdofs[offsetDOF] << ", " << m_nTotDofs_Src << "\n";
345  }
346  }
347  }
348  output_file.flush(); // required here
349  output_file.close();
350  dgll_cgll_col_ldofmap.clear();
351  }
352 
353  {
354  std::ofstream output_file( "target-gids-1.txt" );
355  output_file << "I, GDOF\n";
356  for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
357  output_file << i << ", " << tgt_soln_gdofs[i] << "\n";
358 
359  output_file << "ELEMID, IDOF, GDOF, NDOF\n";
360  m_nTotDofs_Dest = 0;
361 
362  for( unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
363  {
364  output_file << m_remapper->lid_to_gid_tgt[i] << ", " << i << ", " << tgt_soln_gdofs[i] << ", "
365  << m_nTotDofs_Dest << "\n";
366  m_nTotDofs_Dest++;
367  }
368 
369  output_file.flush(); // required here
370  output_file.close();
371  }
372  }
373  }
374 #endif
375 
376  // Now compute the mapping and store it for the covering mesh
377  int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
379  {
380  assert( m_nDofsPEl_Src == 1 );
381  col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
383  src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
384  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
385  srcTagSize = 1;
386  }
387  else
388  {
389  col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
390  col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
391  src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
392  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );MB_CHK_ERR( rval );
393  }
394 
395  // std::cout << "TOnlineMap: Process: " << rank << " and covering entities = [" <<
396  // col_dofmap.size() << ", " << src_soln_gdofs.size() << "]\n"; MPI_Barrier(MPI_COMM_WORLD);
397 
398 #ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
399  unsigned maxSrcIndx = 0;
400 
401  // for ( unsigned j = 0; j < m_covering_source_entities.size(); j++ )
402  std::vector< int > locdofs( srcTagSize );
403  std::map< Node, moab::EntityHandle > mapLocalMBNodes;
404  double elcoords[3];
405  for( unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
406  {
408  rval = m_interface->get_coords( &eh, 1, elcoords );MB_CHK_ERR( rval );
409  Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
410  mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
411  }
412 
413  const NodeVector& nodes = m_remapper->m_covering_source->nodes;
414  for( unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
415  {
416  const Face& face = m_remapper->m_covering_source->faces[j];
417 
418  Node centroid;
419  centroid.x = centroid.y = centroid.z = 0.0;
420  for( unsigned l = 0; l < face.edges.size(); ++l )
421  {
422  centroid.x += nodes[face[l]].x;
423  centroid.y += nodes[face[l]].y;
424  centroid.z += nodes[face[l]].z;
425  }
426  const double factor = 1.0 / face.edges.size();
427  centroid.x *= factor;
428  centroid.y *= factor;
429  centroid.z *= factor;
430 
431  EntityHandle current_eh;
432  if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
433  {
434  current_eh = mapLocalMBNodes[centroid];
435  }
436 
437  rval = m_interface->tag_get_data( m_dofTagSrc, &current_eh, 1, &locdofs[0] );MB_CHK_ERR( rval );
438  for( int p = 0; p < m_nDofsPEl_Src; p++ )
439  {
440  for( int q = 0; q < m_nDofsPEl_Src; q++ )
441  {
442  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
443  const int offsetDOF = p * m_nDofsPEl_Src + q;
444  maxSrcIndx = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
445  std::cout << "Col: " << current_eh << ", " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF
446  << ", " << localDOF << ", " << locdofs[offsetDOF] - 1 << ", " << maxSrcIndx << "\n";
447  }
448  }
449  }
450 #endif
451 
452  m_nTotDofs_SrcCov = 0;
453  if( srcdataGLLNodes == nullptr )
454  { /* we only have a mapping for elements as DoFs */
455  for( unsigned i = 0; i < col_gdofmap.size(); ++i )
456  {
457  assert( src_soln_gdofs[i] > 0 );
458  col_gdofmap[i] = src_soln_gdofs[i] - 1;
459  col_dtoc_dofmap[i] = i;
460  if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
462  }
463  }
464  else
465  {
466  if( isSrcContinuous )
467  dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
468  // Put these remap coefficients into the SparseMatrix map
469  for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
470  {
471  for( int p = 0; p < m_nDofsPEl_Src; p++ )
472  {
473  for( int q = 0; q < m_nDofsPEl_Src; q++ )
474  {
475  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
476  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
477  if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
478  {
480  dgll_cgll_covcol_ldofmap[localDOF] = true;
481  }
482  if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
483  assert( src_soln_gdofs[offsetDOF] > 0 );
484  col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
485  col_dtoc_dofmap[offsetDOF] = localDOF;
486  if( vprint )
487  std::cout << "Col: " << m_remapper->lid_to_gid_covsrc[j] << ", " << offsetDOF << ", "
488  << localDOF << ", " << col_gdofmap[offsetDOF] << ", " << m_nTotDofs_SrcCov << "\n";
489  }
490  }
491  }
492  }
493 
495  {
496  assert( m_nDofsPEl_Src == 1 );
497  srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
499  locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
500  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
501  }
502  else
503  {
504  srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
505  srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
506  locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
507  rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );MB_CHK_ERR( rval );
508  }
509 
510  // Now compute the mapping and store it for the original source mesh
511  m_nTotDofs_Src = 0;
512  if( srcdataGLLNodesSrc == nullptr )
513  { /* we only have a mapping for elements as DoFs */
514  for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
515  {
516  assert( locsrc_soln_gdofs[i] > 0 );
517  srccol_gdofmap[i] = locsrc_soln_gdofs[i] - 1;
518  srccol_dtoc_dofmap[i] = i;
519  m_nTotDofs_Src++;
520  }
521  }
522  else
523  {
524  if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
525  // Put these remap coefficients into the SparseMatrix map
526  for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
527  {
528  for( int p = 0; p < m_nDofsPEl_Src; p++ )
529  {
530  for( int q = 0; q < m_nDofsPEl_Src; q++ )
531  {
532  const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
533  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
534  if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
535  {
536  m_nTotDofs_Src++;
537  dgll_cgll_col_ldofmap[localDOF] = true;
538  }
539  if( !isSrcContinuous ) m_nTotDofs_Src++;
540  assert( locsrc_soln_gdofs[offsetDOF] > 0 );
541  srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
542  srccol_dtoc_dofmap[offsetDOF] = localDOF;
543  }
544  }
545  }
546  }
547 
548  int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
550  {
551  assert( m_nDofsPEl_Dest == 1 );
552  row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
553  row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
554  tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
555  rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
556  tgtTagSize = 1;
557  }
558  else
559  {
560  row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
561  row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
562  tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
563  rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );MB_CHK_ERR( rval );
564  }
565 
566  // Now compute the mapping and store it for the target mesh
567  // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
568  m_nTotDofs_Dest = 0;
569  if( tgtdataGLLNodes == nullptr )
570  { /* we only have a mapping for elements as DoFs */
571  for( unsigned i = 0; i < row_gdofmap.size(); ++i )
572  {
573  assert( tgt_soln_gdofs[i] > 0 );
574  row_gdofmap[i] = tgt_soln_gdofs[i] - 1;
575  row_dtoc_dofmap[i] = i;
576  if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
577  m_nTotDofs_Dest++;
578  }
579  }
580  else
581  {
582  if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
583  // Put these remap coefficients into the SparseMatrix map
584  for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
585  {
586  for( int p = 0; p < m_nDofsPEl_Dest; p++ )
587  {
588  for( int q = 0; q < m_nDofsPEl_Dest; q++ )
589  {
590  const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
591  const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
592  if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
593  {
594  m_nTotDofs_Dest++;
595  dgll_cgll_row_ldofmap[localDOF] = true;
596  }
597  if( !isTgtContinuous ) m_nTotDofs_Dest++;
598  assert( tgt_soln_gdofs[offsetDOF] > 0 );
599  row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
600  row_dtoc_dofmap[offsetDOF] = localDOF;
601  if( vprint )
602  std::cout << "Row: " << m_remapper->lid_to_gid_tgt[j] << ", " << offsetDOF << ", " << localDOF
603  << ", " << row_gdofmap[offsetDOF] << ", " << m_nTotDofs_Dest << "\n";
604  }
605  }
606  }
607  }
608 
609  // Let us also allocate the local representation of the sparse matrix
610 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
611  if( vprint )
612  {
613  std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size()
614  << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
615  // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
616  }
617 #endif
618 
619  // check monotonicity of row_gdofmap and col_gdofmap
620 #ifdef CHECK_INCREASING_DOF
621  for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
622  {
623  if( row_gdofmap[i] > row_gdofmap[i + 1] )
624  std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
625  << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
626  }
627  for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
628  {
629  if( col_gdofmap[i] > col_gdofmap[i + 1] )
630  std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
631  << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
632  }
633 #endif
634 
635  return moab::MB_SUCCESS;
636 }

References ErrorCode, MB_CHK_ERR, and MB_SUCCESS.

◆ SetDOFmapTags()

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

Store the tag names associated with global DoF ids for source and target meshes to be used for mapping.

Parameters
srcDofTagNameThe tag name associated with global DoF ids for the source mesh
tgtDofTagNameThe tag name associated with global DoF ids for the target mesh

Definition at line 128 of file TempestOnlineMap.cpp.

130 {
131  moab::ErrorCode rval;
132 
133  int tagSize = 0;
135  rval =
136  m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
137 
139  {
140  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for source mesh." );
141  }
142  else
143  MB_CHK_ERR( rval );
144 
146  rval =
147  m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
149  {
150  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for target mesh." );
151  }
152  else
153  MB_CHK_ERR( rval );
154 
155  return moab::MB_SUCCESS;
156 }

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

◆ SetMeshInput()

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

Definition at line 447 of file TempestOnlineMap.hpp.

448  {
449  m_meshInput = imesh;
450  };

References m_meshInput.

◆ SetSourceNDofsPerElement()

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

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

Definition at line 581 of file TempestOnlineMap.hpp.

582 {
583  m_nDofsPEl_Src = ns;
584 }

◆ setup_sizes_dimensions()

void moab::TempestOnlineMap::setup_sizes_dimensions ( )
private

Definition at line 90 of file TempestOnlineMap.cpp.

91 {
92  if( m_meshInputCov )
93  {
94  std::vector< std::string > dimNames;
95  std::vector< int > dimSizes;
96  dimNames.push_back( "num_elem" );
97  dimSizes.push_back( m_meshInputCov->faces.size() );
98 
99  this->InitializeSourceDimensions( dimNames, dimSizes );
100  }
101 
102  if( m_meshOutput )
103  {
104  std::vector< std::string > dimNames;
105  std::vector< int > dimSizes;
106  dimNames.push_back( "num_elem" );
107  dimSizes.push_back( m_meshOutput->faces.size() );
108 
109  this->InitializeTargetDimensions( dimNames, dimSizes );
110  }
111 }

Referenced by TempestOnlineMap().

◆ WriteHDF5MapFile()

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

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

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

Definition at line 790 of file TempestOnlineMapIO.cpp.

791 {
792  moab::ErrorCode rval;
793 
794  /**
795  * Need to get the global maximum of number of vertices per element
796  * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
797  *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
798  *when writing this out, other processes may have a different size for the same array. This is
799  *hence a mess to consolidate in h5mtoscrip eventually.
800  **/
801 
802  /* Let us compute all relevant data for the current original source mesh on the process */
803  DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
804  DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
805  DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
807  {
808  this->InitializeCoordinatesFromMeshFV(
809  *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
810  ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
812 
813  vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
814  for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
815  vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
816  }
817  else
818  {
819  DataArray3D< double > dataGLLJacobianSrc;
820  this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
821  dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
822 
823  // Generate the continuous Jacobian for input mesh
824  GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
825 
827  {
828  GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
829  }
830  else
831  {
832  GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
833  }
834  }
835 
837  {
838  this->InitializeCoordinatesFromMeshFV(
839  *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
840  ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ) /* fLatLon = false */,
842 
843  vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
844  for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
845  vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
846  }
847  else
848  {
849  DataArray3D< double > dataGLLJacobianDest;
850  this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
851  dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
852 
853  // Generate the continuous Jacobian for input mesh
854  GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
855 
857  {
858  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
859  }
860  else
861  {
862  GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
863  }
864  }
865 
866  moab::EntityHandle& m_meshOverlapSet = m_remapper->m_overlap_set;
867  int tot_src_ents = m_remapper->m_source_entities.size();
868  int tot_tgt_ents = m_remapper->m_target_entities.size();
869  int tot_src_size = dSourceCenterLon.GetRows();
870  int tot_tgt_size = m_dTargetCenterLon.GetRows();
871  int tot_vsrc_size = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
872  int tot_vtgt_size = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
873 
874  const int weightMatNNZ = m_weightMatrix.nonZeros();
875  moab::Tag tagMapMetaData, tagMapIndexRow, tagMapIndexCol, tagMapValues, srcEleIDs, tgtEleIDs;
876  rval = m_interface->tag_get_handle( "SMAT_DATA", 13, moab::MB_TYPE_INTEGER, tagMapMetaData,
877  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
878  rval = m_interface->tag_get_handle( "SMAT_ROWS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexRow,
879  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
880  rval = m_interface->tag_get_handle( "SMAT_COLS", weightMatNNZ, moab::MB_TYPE_INTEGER, tagMapIndexCol,
881  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
882  rval = m_interface->tag_get_handle( "SMAT_VALS", weightMatNNZ, moab::MB_TYPE_DOUBLE, tagMapValues,
883  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
884  rval = m_interface->tag_get_handle( "SourceGIDS", tot_src_size, moab::MB_TYPE_INTEGER, srcEleIDs,
885  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
886  rval = m_interface->tag_get_handle( "TargetGIDS", tot_tgt_size, moab::MB_TYPE_INTEGER, tgtEleIDs,
887  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
888  moab::Tag srcAreaValues, tgtAreaValues;
889  rval = m_interface->tag_get_handle( "SourceAreas", tot_src_size, moab::MB_TYPE_DOUBLE, srcAreaValues,
890  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
891  rval = m_interface->tag_get_handle( "TargetAreas", tot_tgt_size, moab::MB_TYPE_DOUBLE, tgtAreaValues,
892  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
893  moab::Tag tagSrcCoordsCLon, tagSrcCoordsCLat, tagTgtCoordsCLon, tagTgtCoordsCLat;
894  rval = m_interface->tag_get_handle( "SourceCoordCenterLon", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLon,
895  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
896  rval = m_interface->tag_get_handle( "SourceCoordCenterLat", tot_src_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsCLat,
897  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
898  rval = m_interface->tag_get_handle( "TargetCoordCenterLon", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLon,
899  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
900  rval = m_interface->tag_get_handle( "TargetCoordCenterLat", tot_tgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsCLat,
901  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
902  moab::Tag tagSrcCoordsVLon, tagSrcCoordsVLat, tagTgtCoordsVLon, tagTgtCoordsVLat;
903  rval = m_interface->tag_get_handle( "SourceCoordVertexLon", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLon,
904  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
905  rval = m_interface->tag_get_handle( "SourceCoordVertexLat", tot_vsrc_size, moab::MB_TYPE_DOUBLE, tagSrcCoordsVLat,
906  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
907  rval = m_interface->tag_get_handle( "TargetCoordVertexLon", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLon,
908  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
909  rval = m_interface->tag_get_handle( "TargetCoordVertexLat", tot_vtgt_size, moab::MB_TYPE_DOUBLE, tagTgtCoordsVLat,
910  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
911  moab::Tag srcMaskValues, tgtMaskValues;
912  if( m_iSourceMask.IsAttached() )
913  {
914  rval = m_interface->tag_get_handle( "SourceMask", m_iSourceMask.GetRows(), moab::MB_TYPE_INTEGER, srcMaskValues,
915  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
916  }
917  if( m_iTargetMask.IsAttached() )
918  {
919  rval = m_interface->tag_get_handle( "TargetMask", m_iTargetMask.GetRows(), moab::MB_TYPE_INTEGER, tgtMaskValues,
920  moab::MB_TAG_CREAT | moab::MB_TAG_SPARSE | moab::MB_TAG_VARLEN );MB_CHK_SET_ERR( rval, "Retrieving tag handles failed" );
921  }
922 
923  std::vector< int > smatrowvals( weightMatNNZ ), smatcolvals( weightMatNNZ );
924  std::vector< double > smatvals( weightMatNNZ );
925  // const double* smatvals = m_weightMatrix.valuePtr();
926  // Loop over the matrix entries and find the max global ID for rows and columns
927  for( int k = 0, offset = 0; k < m_weightMatrix.outerSize(); ++k )
928  {
929  for( moab::TempestOnlineMap::WeightMatrix::InnerIterator it( m_weightMatrix, k ); it; ++it, ++offset )
930  {
931  smatrowvals[offset] = this->GetRowGlobalDoF( it.row() );
932  smatcolvals[offset] = this->GetColGlobalDoF( it.col() );
933  smatvals[offset] = it.value();
934  }
935  }
936 
937  /* Set the global IDs for the DoFs */
938  ////
939  // col_gdofmap [ col_ldofmap [ 0 : local_ndofs ] ] = GDOF
940  // row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
941  ////
942  int maxrow = 0, maxcol = 0;
943  std::vector< int > src_global_dofs( tot_src_size ), tgt_global_dofs( tot_tgt_size );
944  for( int i = 0; i < tot_src_size; ++i )
945  {
946  src_global_dofs[i] = srccol_gdofmap[i];
947  maxcol = ( src_global_dofs[i] > maxcol ) ? src_global_dofs[i] : maxcol;
948  }
949 
950  for( int i = 0; i < tot_tgt_size; ++i )
951  {
952  tgt_global_dofs[i] = row_gdofmap[i];
953  maxrow = ( tgt_global_dofs[i] > maxrow ) ? tgt_global_dofs[i] : maxrow;
954  }
955 
956  ///////////////////////////////////////////////////////////////////////////
957  // The metadata in H5M file contains the following data:
958  //
959  // 1. n_a: Total source entities: (number of elements in source mesh)
960  // 2. n_b: Total target entities: (number of elements in target mesh)
961  // 3. nv_a: Max edge size of elements in source mesh
962  // 4. nv_b: Max edge size of elements in target mesh
963  // 5. maxrows: Number of rows in remap weight matrix
964  // 6. maxcols: Number of cols in remap weight matrix
965  // 7. nnz: Number of total nnz in sparse remap weight matrix
966  // 8. np_a: The order of the field description on the source mesh: >= 1
967  // 9. np_b: The order of the field description on the target mesh: >= 1
968  // 10. method_a: The type of discretization for field on source mesh: [0 = FV, 1 = cGLL, 2 =
969  // dGLL]
970  // 11. method_b: The type of discretization for field on target mesh: [0 = FV, 1 = cGLL, 2 =
971  // dGLL]
972  // 12. conserved: Flag to specify whether the remap operator has conservation constraints: [0,
973  // 1]
974  // 13. monotonicity: Flags to specify whether the remap operator has monotonicity constraints:
975  // [0, 1, 2]
976  //
977  ///////////////////////////////////////////////////////////////////////////
978  int map_disc_details[6];
979  map_disc_details[0] = m_nDofsPEl_Src;
980  map_disc_details[1] = m_nDofsPEl_Dest;
982  ? 0
983  : ( m_srcDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
985  ? 0
986  : ( m_destDiscType == DiscretizationType_CGLL ? 1 : 2 ) );
987  map_disc_details[4] = ( m_bConserved ? 1 : 0 );
988  map_disc_details[5] = m_iMonotonicity;
989 
990 #ifdef MOAB_HAVE_MPI
991  int loc_smatmetadata[13] = { tot_src_ents,
992  tot_tgt_ents,
995  maxrow + 1,
996  maxcol + 1,
997  weightMatNNZ,
998  map_disc_details[0],
999  map_disc_details[1],
1000  map_disc_details[2],
1001  map_disc_details[3],
1002  map_disc_details[4],
1003  map_disc_details[5] };
1004  rval = m_interface->tag_set_data( tagMapMetaData, &m_meshOverlapSet, 1, &loc_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1005  int glb_smatmetadata[13] = { 0,
1006  0,
1007  0,
1008  0,
1009  0,
1010  0,
1011  0,
1012  map_disc_details[0],
1013  map_disc_details[1],
1014  map_disc_details[2],
1015  map_disc_details[3],
1016  map_disc_details[4],
1017  map_disc_details[5] };
1018  int loc_buf[7] = {
1019  tot_src_ents, tot_tgt_ents, weightMatNNZ, m_remapper->max_source_edges, m_remapper->max_target_edges,
1020  maxrow, maxcol };
1021  int glb_buf[4] = { 0, 0, 0, 0 };
1022  MPI_Reduce( &loc_buf[0], &glb_buf[0], 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
1023  glb_smatmetadata[0] = glb_buf[0];
1024  glb_smatmetadata[1] = glb_buf[1];
1025  glb_smatmetadata[6] = glb_buf[2];
1026  MPI_Reduce( &loc_buf[3], &glb_buf[0], 4, MPI_INT, MPI_MAX, 0, m_pcomm->comm() );
1027  glb_smatmetadata[2] = glb_buf[0];
1028  glb_smatmetadata[3] = glb_buf[1];
1029  glb_smatmetadata[4] = glb_buf[2];
1030  glb_smatmetadata[5] = glb_buf[3];
1031 #else
1032  int glb_smatmetadata[13] = { tot_src_ents,
1033  tot_tgt_ents,
1036  maxrow,
1037  maxcol,
1038  weightMatNNZ,
1039  map_disc_details[0],
1040  map_disc_details[1],
1041  map_disc_details[2],
1042  map_disc_details[3],
1043  map_disc_details[4],
1044  map_disc_details[5] };
1045 #endif
1046  // These values represent number of rows and columns. So should be 1-based.
1047  glb_smatmetadata[4]++;
1048  glb_smatmetadata[5]++;
1049 
1050  if( this->is_root )
1051  {
1052  std::cout << " " << this->rank << " Writing remap weights with size [" << glb_smatmetadata[4] << " X "
1053  << glb_smatmetadata[5] << "] and NNZ = " << glb_smatmetadata[6] << std::endl;
1054  EntityHandle root_set = 0;
1055  rval = m_interface->tag_set_data( tagMapMetaData, &root_set, 1, &glb_smatmetadata[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1056  }
1057 
1058  int dsize;
1059  const int numval = weightMatNNZ;
1060  const void* smatrowvals_d = smatrowvals.data();
1061  const void* smatcolvals_d = smatcolvals.data();
1062  const void* smatvals_d = smatvals.data();
1063  rval = m_interface->tag_set_by_ptr( tagMapIndexRow, &m_meshOverlapSet, 1, &smatrowvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1064  rval = m_interface->tag_set_by_ptr( tagMapIndexCol, &m_meshOverlapSet, 1, &smatcolvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1065  rval = m_interface->tag_set_by_ptr( tagMapValues, &m_meshOverlapSet, 1, &smatvals_d, &numval );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1066 
1067  /* Set the global IDs for the DoFs */
1068  const void* srceleidvals_d = src_global_dofs.data();
1069  const void* tgteleidvals_d = tgt_global_dofs.data();
1070  dsize = src_global_dofs.size();
1071  rval = m_interface->tag_set_by_ptr( srcEleIDs, &m_meshOverlapSet, 1, &srceleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1072  dsize = tgt_global_dofs.size();
1073  rval = m_interface->tag_set_by_ptr( tgtEleIDs, &m_meshOverlapSet, 1, &tgteleidvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1074 
1075  /* Set the source and target areas */
1076  const void* srcareavals_d = vecSourceFaceArea;
1077  const void* tgtareavals_d = vecTargetFaceArea;
1078  dsize = tot_src_size;
1079  rval = m_interface->tag_set_by_ptr( srcAreaValues, &m_meshOverlapSet, 1, &srcareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1080  dsize = tot_tgt_size;
1081  rval = m_interface->tag_set_by_ptr( tgtAreaValues, &m_meshOverlapSet, 1, &tgtareavals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1082 
1083  /* Set the coordinates for source and target center vertices */
1084  const void* srccoordsclonvals_d = &dSourceCenterLon[0];
1085  const void* srccoordsclatvals_d = &dSourceCenterLat[0];
1086  dsize = dSourceCenterLon.GetRows();
1087  rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLon, &m_meshOverlapSet, 1, &srccoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1088  rval = m_interface->tag_set_by_ptr( tagSrcCoordsCLat, &m_meshOverlapSet, 1, &srccoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1089  const void* tgtcoordsclonvals_d = &m_dTargetCenterLon[0];
1090  const void* tgtcoordsclatvals_d = &m_dTargetCenterLat[0];
1091  dsize = vecTargetFaceArea.GetRows();
1092  rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLon, &m_meshOverlapSet, 1, &tgtcoordsclonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1093  rval = m_interface->tag_set_by_ptr( tagTgtCoordsCLat, &m_meshOverlapSet, 1, &tgtcoordsclatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1094 
1095  /* Set the coordinates for source and target element vertices */
1096  const void* srccoordsvlonvals_d = &( dSourceVertexLon[0][0] );
1097  const void* srccoordsvlatvals_d = &( dSourceVertexLat[0][0] );
1098  dsize = dSourceVertexLon.GetRows() * dSourceVertexLon.GetColumns();
1099  rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLon, &m_meshOverlapSet, 1, &srccoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1100  rval = m_interface->tag_set_by_ptr( tagSrcCoordsVLat, &m_meshOverlapSet, 1, &srccoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1101  const void* tgtcoordsvlonvals_d = &( m_dTargetVertexLon[0][0] );
1102  const void* tgtcoordsvlatvals_d = &( m_dTargetVertexLat[0][0] );
1103  dsize = m_dTargetVertexLon.GetRows() * m_dTargetVertexLon.GetColumns();
1104  rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLon, &m_meshOverlapSet, 1, &tgtcoordsvlonvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1105  rval = m_interface->tag_set_by_ptr( tagTgtCoordsVLat, &m_meshOverlapSet, 1, &tgtcoordsvlatvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1106 
1107  /* Set the masks for source and target meshes if available */
1108  if( m_iSourceMask.IsAttached() )
1109  {
1110  const void* srcmaskvals_d = m_iSourceMask;
1111  dsize = m_iSourceMask.GetRows();
1112  rval = m_interface->tag_set_by_ptr( srcMaskValues, &m_meshOverlapSet, 1, &srcmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1113  }
1114 
1115  if( m_iTargetMask.IsAttached() )
1116  {
1117  const void* tgtmaskvals_d = m_iTargetMask;
1118  dsize = m_iTargetMask.GetRows();
1119  rval = m_interface->tag_set_by_ptr( tgtMaskValues, &m_meshOverlapSet, 1, &tgtmaskvals_d, &dsize );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1120  }
1121 
1122 #ifdef MOAB_HAVE_MPI
1123  const char* writeOptions = ( this->size > 1 ? "PARALLEL=WRITE_PART" : "" );
1124 #else
1125  const char* writeOptions = "";
1126 #endif
1127 
1128  // EntityHandle sets[3] = {m_remapper->m_source_set, m_remapper->m_target_set, m_remapper->m_overlap_set};
1129  EntityHandle sets[1] = { m_remapper->m_overlap_set };
1130  rval = m_interface->write_file( strOutputFile.c_str(), NULL, writeOptions, sets, 1 );MB_CHK_ERR( rval );
1131 
1132 #ifdef WRITE_SCRIP_FILE
1133  sstr.str( "" );
1134  sstr << ctx.outFilename.substr( 0, lastindex ) << "_" << proc_id << ".nc";
1135  std::map< std::string, std::string > mapAttributes;
1136  mapAttributes["Creator"] = "MOAB mbtempest workflow";
1137  if( !ctx.proc_id ) std::cout << "Writing offline map to file: " << sstr.str() << std::endl;
1138  this->Write( strOutputFile.c_str(), mapAttributes, NcFile::Netcdf4 );
1139  sstr.str( "" );
1140 #endif
1141 
1142  return moab::MB_SUCCESS;
1143 }

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

◆ WriteParallelMap()

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

Write the TempestOnlineMap to a parallel NetCDF file.

Definition at line 156 of file TempestOnlineMapIO.cpp.

158 {
159  moab::ErrorCode rval;
160 
161  size_t lastindex = strFilename.find_last_of( "." );
162  std::string extension = strFilename.substr( lastindex + 1, strFilename.size() );
163 
164  // Write the map file to disk in parallel
165  if( extension == "nc" )
166  {
167  /* Invoke the actual call to write the parallel map to disk in SCRIP format */
168  rval = this->WriteSCRIPMapFile( strFilename.c_str(), attrMap );MB_CHK_ERR( rval );
169  }
170  else
171  {
172  /* Write to the parallel H5M format */
173  rval = this->WriteHDF5MapFile( strFilename.c_str() );MB_CHK_ERR( rval );
174  }
175 
176  return rval;
177 }

References ErrorCode, and MB_CHK_ERR.

Referenced by main().

◆ WriteSCRIPMapFile()

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

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

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

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

Definition at line 181 of file TempestOnlineMapIO.cpp.

183 {
184  NcError error( NcError::silent_nonfatal );
185 
186 #ifdef MOAB_HAVE_NETCDFPAR
187  bool is_independent = true;
188  ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcFile::Replace, NcFile::Netcdf4 );
189  // ParNcFile ncMap( m_pcomm->comm(), MPI_INFO_NULL, strFilename.c_str(), NcmpiFile::replace, NcmpiFile::classic5 );
190 #else
191  NcFile ncMap( strFilename.c_str(), NcFile::Replace );
192 #endif
193 
194  if( !ncMap.is_valid() )
195  {
196  _EXCEPTION1( "Unable to open output map file \"%s\"", strFilename.c_str() );
197  }
198 
199  // Attributes
200  // ncMap.add_att( "Title", "MOAB-TempestRemap Online Regridding Weight Generator" );
201  auto it = attrMap.begin();
202  while( it != attrMap.end() )
203  {
204  // set the map attributes
205  ncMap.add_att( it->first.c_str(), it->second.c_str() );
206  // increment iterator
207  it++;
208  }
209 
210  /**
211  * Need to get the global maximum of number of vertices per element
212  * Key issue is that when calling InitializeCoordinatesFromMeshFV, the allocation for
213  *dVertexLon/dVertexLat are made based on the maximum vertices in the current process. However,
214  *when writing this out, other processes may have a different size for the same array. This is
215  *hence a mess to consolidate in h5mtoscrip eventually.
216  **/
217 
218  /* Let us compute all relevant data for the current original source mesh on the process */
219  DataArray1D< double > vecSourceFaceArea, vecTargetFaceArea;
220  DataArray1D< double > dSourceCenterLon, dSourceCenterLat, dTargetCenterLon, dTargetCenterLat;
221  DataArray2D< double > dSourceVertexLon, dSourceVertexLat, dTargetVertexLon, dTargetVertexLat;
223  {
224  this->InitializeCoordinatesFromMeshFV(
225  *m_meshInput, dSourceCenterLon, dSourceCenterLat, dSourceVertexLon, dSourceVertexLat,
226  ( this->m_remapper->m_source_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
228 
229  vecSourceFaceArea.Allocate( m_meshInput->vecFaceArea.GetRows() );
230  for( unsigned i = 0; i < m_meshInput->vecFaceArea.GetRows(); ++i )
231  vecSourceFaceArea[i] = m_meshInput->vecFaceArea[i];
232  }
233  else
234  {
235  DataArray3D< double > dataGLLJacobianSrc;
236  this->InitializeCoordinatesFromMeshFE( *m_meshInput, m_nDofsPEl_Src, dataGLLNodesSrc, dSourceCenterLon,
237  dSourceCenterLat, dSourceVertexLon, dSourceVertexLat );
238 
239  // Generate the continuous Jacobian for input mesh
240  GenerateMetaData( *m_meshInput, m_nDofsPEl_Src, false /* fBubble */, dataGLLNodesSrc, dataGLLJacobianSrc );
241 
243  {
244  GenerateUniqueJacobian( dataGLLNodesSrc, dataGLLJacobianSrc, vecSourceFaceArea );
245  }
246  else
247  {
248  GenerateDiscontinuousJacobian( dataGLLJacobianSrc, vecSourceFaceArea );
249  }
250  }
251 
253  {
254  this->InitializeCoordinatesFromMeshFV(
255  *m_meshOutput, dTargetCenterLon, dTargetCenterLat, dTargetVertexLon, dTargetVertexLat,
256  ( this->m_remapper->m_target_type == moab::TempestRemapper::RLL ), /* fLatLon = false */
258 
259  vecTargetFaceArea.Allocate( m_meshOutput->vecFaceArea.GetRows() );
260  for( unsigned i = 0; i < m_meshOutput->vecFaceArea.GetRows(); ++i )
261  {
262  vecTargetFaceArea[i] = m_meshOutput->vecFaceArea[i];
263  }
264  }
265  else
266  {
267  DataArray3D< double > dataGLLJacobianDest;
268  this->InitializeCoordinatesFromMeshFE( *m_meshOutput, m_nDofsPEl_Dest, dataGLLNodesDest, dTargetCenterLon,
269  dTargetCenterLat, dTargetVertexLon, dTargetVertexLat );
270 
271  // Generate the continuous Jacobian for input mesh
272  GenerateMetaData( *m_meshOutput, m_nDofsPEl_Dest, false /* fBubble */, dataGLLNodesDest, dataGLLJacobianDest );
273 
275  {
276  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianDest, vecTargetFaceArea );
277  }
278  else
279  {
280  GenerateDiscontinuousJacobian( dataGLLJacobianDest, vecTargetFaceArea );
281  }
282  }
283 
284  // Map dimensions
285  unsigned nA = ( vecSourceFaceArea.GetRows() );
286  unsigned nB = ( vecTargetFaceArea.GetRows() );
287 
288  std::vector< int > masksA, masksB;
289  ErrorCode rval = m_remapper->GetIMasks( moab::Remapper::SourceMesh, masksA );MB_CHK_SET_ERR( rval, "Trouble getting masks for source" );
290  rval = m_remapper->GetIMasks( moab::Remapper::TargetMesh, masksB );MB_CHK_SET_ERR( rval, "Trouble getting masks for target" );
291 
292  // Number of nodes per Face
293  int nSourceNodesPerFace = dSourceVertexLon.GetColumns();
294  int nTargetNodesPerFace = dTargetVertexLon.GetColumns();
295 
296  // if source or target cells have triangles at poles, center of those triangles need to come from
297  // the original quad, not from center in 3d, converted to 2d again
298  // start copy OnlineMap.cpp tempestremap
299  // right now, do this only for source mesh; copy the logic for target mesh
300  for( unsigned i = 0; i < nA; i++ )
301  {
302  const Face& face = m_meshInput->faces[i];
303 
304  int nNodes = face.edges.size();
305  int indexNodeAtPole = -1;
306  if( 3 == nNodes ) // check if one node at the poles
307  {
308  for( int j = 0; j < nNodes; j++ )
309  if( fabs( fabs( dSourceVertexLat[i][j] ) - 90.0 ) < 1.0e-12 )
310  {
311  indexNodeAtPole = j;
312  break;
313  }
314  }
315  if( indexNodeAtPole < 0 ) continue; // continue i loop, do nothing
316  // recompute center of cell, from 3d data; add one 2 nodes at pole, and average
317  int nodeAtPole = face[indexNodeAtPole]; // use the overloaded operator
318  Node nodePole = m_meshInput->nodes[nodeAtPole];
319  Node newCenter = nodePole * 2;
320  for( int j = 1; j < nNodes; j++ )
321  {
322  int indexi = ( indexNodeAtPole + j ) % nNodes; // nNodes is 3 !
323  const Node& node = m_meshInput->nodes[face[indexi]];
324  newCenter = newCenter + node;
325  }
326  newCenter = newCenter * 0.25;
327  newCenter = newCenter.Normalized();
328 
329 #ifdef VERBOSE
330  double iniLon = dSourceCenterLon[i], iniLat = dSourceCenterLat[i];
331 #endif
332  // dSourceCenterLon, dSourceCenterLat
333  XYZtoRLL_Deg( newCenter.x, newCenter.y, newCenter.z, dSourceCenterLon[i], dSourceCenterLat[i] );
334 #ifdef VERBOSE
335  std::cout << " modify center of triangle from " << iniLon << " " << iniLat << " to " << dSourceCenterLon[i]
336  << " " << dSourceCenterLat[i] << "\n";
337 #endif
338  }
339 
340  // first move data if in parallel
341 #if defined( MOAB_HAVE_MPI )
342  int max_row_dof, max_col_dof; // output; arrays will be re-distributed in chunks [maxdof/size]
343  // if (size > 1)
344  {
345  int ierr = rearrange_arrays_by_dofs( srccol_gdofmap, vecSourceFaceArea, dSourceCenterLon, dSourceCenterLat,
346  dSourceVertexLon, dSourceVertexLat, masksA, nA, nSourceNodesPerFace,
347  max_col_dof ); // now nA will be close to maxdof/size
348  if( ierr != 0 )
349  {
350  _EXCEPTION1( "Unable to arrange source data %d ", nA );
351  }
352  // rearrange target data: (nB)
353  //
354  ierr = rearrange_arrays_by_dofs( row_gdofmap, vecTargetFaceArea, dTargetCenterLon, dTargetCenterLat,
355  dTargetVertexLon, dTargetVertexLat, masksB, nB, nTargetNodesPerFace,
356  max_row_dof ); // now nA will be close to maxdof/size
357  if( ierr != 0 )
358  {
359  _EXCEPTION1( "Unable to arrange target data %d ", nB );
360  }
361  }
362 #endif
363 
364  // Number of non-zeros in the remap matrix operator
365  int nS = m_weightMatrix.nonZeros();
366 
367 #if defined( MOAB_HAVE_MPI ) && defined( MOAB_HAVE_NETCDFPAR )
368  int locbuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
369  int offbuf[3] = { 0, 0, 0 };
370  int globuf[5] = { 0, 0, 0, 0, 0 };
371  MPI_Scan( locbuf, offbuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
372  MPI_Allreduce( locbuf, globuf, 3, MPI_INT, MPI_SUM, m_pcomm->comm() );
373  MPI_Allreduce( &locbuf[3], &globuf[3], 2, MPI_INT, MPI_MAX, m_pcomm->comm() );
374 
375  // MPI_Scan is inclusive of data in current rank; modify accordingly.
376  offbuf[0] -= nA;
377  offbuf[1] -= nB;
378  offbuf[2] -= nS;
379 
380 #else
381  int offbuf[3] = { 0, 0, 0 };
382  int globuf[5] = { (int)nA, (int)nB, nS, nSourceNodesPerFace, nTargetNodesPerFace };
383 #endif
384 
385  std::vector< std::string > srcdimNames, tgtdimNames;
386  std::vector< int > srcdimSizes, tgtdimSizes;
387  {
389  {
390  srcdimNames.push_back( "lat" );
391  srcdimNames.push_back( "lon" );
392  srcdimSizes.resize( 2, 0 );
393  srcdimSizes[0] = m_remapper->m_source_metadata[0];
394  srcdimSizes[1] = m_remapper->m_source_metadata[1];
395  }
396  else
397  {
398  srcdimNames.push_back( "num_elem" );
399  srcdimSizes.push_back( globuf[0] );
400  }
401 
403  {
404  tgtdimNames.push_back( "lat" );
405  tgtdimNames.push_back( "lon" );
406  tgtdimSizes.resize( 2, 0 );
407  tgtdimSizes[0] = m_remapper->m_target_metadata[0];
408  tgtdimSizes[1] = m_remapper->m_target_metadata[1];
409  }
410  else
411  {
412  tgtdimNames.push_back( "num_elem" );
413  tgtdimSizes.push_back( globuf[1] );
414  }
415  }
416 
417  // Write output dimensions entries
418  unsigned nSrcGridDims = ( srcdimSizes.size() );
419  unsigned nDstGridDims = ( tgtdimSizes.size() );
420 
421  NcDim* dimSrcGridRank = ncMap.add_dim( "src_grid_rank", nSrcGridDims );
422  NcDim* dimDstGridRank = ncMap.add_dim( "dst_grid_rank", nDstGridDims );
423 
424  NcVar* varSrcGridDims = ncMap.add_var( "src_grid_dims", ncInt, dimSrcGridRank );
425  NcVar* varDstGridDims = ncMap.add_var( "dst_grid_dims", ncInt, dimDstGridRank );
426 
427 #ifdef MOAB_HAVE_NETCDFPAR
428  ncMap.enable_var_par_access( varSrcGridDims, is_independent );
429  ncMap.enable_var_par_access( varDstGridDims, is_independent );
430 #endif
431 
432  // write dimension names
433  {
434  char szDim[64];
435  for( unsigned i = 0; i < srcdimSizes.size(); i++ )
436  {
437  varSrcGridDims->set_cur( nSrcGridDims - i - 1 );
438  varSrcGridDims->put( &( srcdimSizes[nSrcGridDims - i - 1] ), 1 );
439  }
440 
441  for( unsigned i = 0; i < srcdimSizes.size(); i++ )
442  {
443  snprintf( szDim, 64, "name%i", i );
444  varSrcGridDims->add_att( szDim, srcdimNames[nSrcGridDims - i - 1].c_str() );
445  }
446 
447  for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
448  {
449  varDstGridDims->set_cur( nDstGridDims - i - 1 );
450  varDstGridDims->put( &( tgtdimSizes[nDstGridDims - i - 1] ), 1 );
451  }
452 
453  for( unsigned i = 0; i < tgtdimSizes.size(); i++ )
454  {
455  snprintf( szDim, 64, "name%i", i );
456  varDstGridDims->add_att( szDim, tgtdimNames[nDstGridDims - i - 1].c_str() );
457  }
458  }
459 
460  // Source and Target mesh resolutions
461  NcDim* dimNA = ncMap.add_dim( "n_a", globuf[0] );
462  NcDim* dimNB = ncMap.add_dim( "n_b", globuf[1] );
463 
464  // Number of nodes per Face
465  NcDim* dimNVA = ncMap.add_dim( "nv_a", globuf[3] );
466  NcDim* dimNVB = ncMap.add_dim( "nv_b", globuf[4] );
467 
468  // Write coordinates
469  NcVar* varYCA = ncMap.add_var( "yc_a", ncDouble, dimNA );
470  NcVar* varYCB = ncMap.add_var( "yc_b", ncDouble, dimNB );
471 
472  NcVar* varXCA = ncMap.add_var( "xc_a", ncDouble, dimNA );
473  NcVar* varXCB = ncMap.add_var( "xc_b", ncDouble, dimNB );
474 
475  NcVar* varYVA = ncMap.add_var( "yv_a", ncDouble, dimNA, dimNVA );
476  NcVar* varYVB = ncMap.add_var( "yv_b", ncDouble, dimNB, dimNVB );
477 
478  NcVar* varXVA = ncMap.add_var( "xv_a", ncDouble, dimNA, dimNVA );
479  NcVar* varXVB = ncMap.add_var( "xv_b", ncDouble, dimNB, dimNVB );
480 
481  // Write masks
482  NcVar* varMaskA = ncMap.add_var( "mask_a", ncInt, dimNA );
483  NcVar* varMaskB = ncMap.add_var( "mask_b", ncInt, dimNB );
484 
485 #ifdef MOAB_HAVE_NETCDFPAR
486  ncMap.enable_var_par_access( varYCA, is_independent );
487  ncMap.enable_var_par_access( varYCB, is_independent );
488  ncMap.enable_var_par_access( varXCA, is_independent );
489  ncMap.enable_var_par_access( varXCB, is_independent );
490  ncMap.enable_var_par_access( varYVA, is_independent );
491  ncMap.enable_var_par_access( varYVB, is_independent );
492  ncMap.enable_var_par_access( varXVA, is_independent );
493  ncMap.enable_var_par_access( varXVB, is_independent );
494  ncMap.enable_var_par_access( varMaskA, is_independent );
495  ncMap.enable_var_par_access( varMaskB, is_independent );
496 #endif
497 
498  varYCA->add_att( "units", "degrees" );
499  varYCB->add_att( "units", "degrees" );
500 
501  varXCA->add_att( "units", "degrees" );
502  varXCB->add_att( "units", "degrees" );
503 
504  varYVA->add_att( "units", "degrees" );
505  varYVB->add_att( "units", "degrees" );
506 
507  varXVA->add_att( "units", "degrees" );
508  varXVB->add_att( "units", "degrees" );
509 
510  // Verify dimensionality
511  if( dSourceCenterLon.GetRows() != nA )
512  {
513  _EXCEPTIONT( "Mismatch between dSourceCenterLon and nA" );
514  }
515  if( dSourceCenterLat.GetRows() != nA )
516  {
517  _EXCEPTIONT( "Mismatch between dSourceCenterLat and nA" );
518  }
519  if( dTargetCenterLon.GetRows() != nB )
520  {
521  _EXCEPTIONT( "Mismatch between dTargetCenterLon and nB" );
522  }
523  if( dTargetCenterLat.GetRows() != nB )
524  {
525  _EXCEPTIONT( "Mismatch between dTargetCenterLat and nB" );
526  }
527  if( dSourceVertexLon.GetRows() != nA )
528  {
529  _EXCEPTIONT( "Mismatch between dSourceVertexLon and nA" );
530  }
531  if( dSourceVertexLat.GetRows() != nA )
532  {
533  _EXCEPTIONT( "Mismatch between dSourceVertexLat and nA" );
534  }
535  if( dTargetVertexLon.GetRows() != nB )
536  {
537  _EXCEPTIONT( "Mismatch between dTargetVertexLon and nB" );
538  }
539  if( dTargetVertexLat.GetRows() != nB )
540  {
541  _EXCEPTIONT( "Mismatch between dTargetVertexLat and nB" );
542  }
543 
544  varYCA->set_cur( (long)offbuf[0] );
545  varYCA->put( &( dSourceCenterLat[0] ), nA );
546  varYCB->set_cur( (long)offbuf[1] );
547  varYCB->put( &( dTargetCenterLat[0] ), nB );
548 
549  varXCA->set_cur( (long)offbuf[0] );
550  varXCA->put( &( dSourceCenterLon[0] ), nA );
551  varXCB->set_cur( (long)offbuf[1] );
552  varXCB->put( &( dTargetCenterLon[0] ), nB );
553 
554  varYVA->set_cur( (long)offbuf[0] );
555  varYVA->put( &( dSourceVertexLat[0][0] ), nA, nSourceNodesPerFace );
556  varYVB->set_cur( (long)offbuf[1] );
557  varYVB->put( &( dTargetVertexLat[0][0] ), nB, nTargetNodesPerFace );
558 
559  varXVA->set_cur( (long)offbuf[0] );
560  varXVA->put( &( dSourceVertexLon[0][0] ), nA, nSourceNodesPerFace );
561  varXVB->set_cur( (long)offbuf[1] );
562  varXVB->put( &( dTargetVertexLon[0][0] ), nB, nTargetNodesPerFace );
563 
564  varMaskA->set_cur( (long)offbuf[0] );
565  varMaskA->put( &( masksA[0] ), nA );
566  varMaskB->set_cur( (long)offbuf[1] );
567  varMaskB->put( &( masksB[0] ), nB );
568 
569  // Write areas
570  NcVar* varAreaA = ncMap.add_var( "area_a", ncDouble, dimNA );
571 #ifdef MOAB_HAVE_NETCDFPAR
572  ncMap.enable_var_par_access( varAreaA, is_independent );
573 #endif
574  varAreaA->set_cur( (long)offbuf[0] );
575  varAreaA->put( &( vecSourceFaceArea[0] ), nA );
576 
577  NcVar* varAreaB = ncMap.add_var( "area_b", ncDouble, dimNB );
578 #ifdef MOAB_HAVE_NETCDFPAR
579  ncMap.enable_var_par_access( varAreaB, is_independent );
580 #endif
581  varAreaB->set_cur( (long)offbuf[1] );
582  varAreaB->put( &( vecTargetFaceArea[0] ), nB );
583 
584  // Write SparseMatrix entries
585  DataArray1D< int > vecRow( nS );
586  DataArray1D< int > vecCol( nS );
587  DataArray1D< double > vecS( nS );
588  DataArray1D< double > dFracA( nA );
589  DataArray1D< double > dFracB( nB );
590 
591  moab::TupleList tlValRow, tlValCol;
592  unsigned numr = 1; //
593  // value has to be sent to processor row/nB for for fracA and col/nA for fracB
594  // vecTargetArea (indexRow ) has to be sent for fracA (index col?)
595  // vecTargetFaceArea will have to be sent to col index, with its index !
596  tlValRow.initialize( 2, 0, 0, numr, nS ); // to proc(row), global row , value
597  tlValCol.initialize( 3, 0, 0, numr, nS ); // to proc(col), global row / col, value
598  tlValRow.enableWriteAccess();
599  tlValCol.enableWriteAccess();
600  /*
601  *
602  dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
603  dFracB[ row ] += val ;
604  */
605  int offset = 0;
606 #if defined( MOAB_HAVE_MPI )
607  int nAbase = ( max_col_dof + 1 ) / size; // it is nA, except last rank ( == size - 1 )
608  int nBbase = ( max_row_dof + 1 ) / size; // it is nB, except last rank ( == size - 1 )
609 #endif
610  for( int i = 0; i < m_weightMatrix.outerSize(); ++i )
611  {
612  for( WeightMatrix::InnerIterator it( m_weightMatrix, i ); it; ++it )
613  {
614  vecRow[offset] = 1 + this->GetRowGlobalDoF( it.row() ); // row index
615  vecCol[offset] = 1 + this->GetColGlobalDoF( it.col() ); // col index
616  vecS[offset] = it.value(); // value
617 
618 #if defined( MOAB_HAVE_MPI )
619  {
620  // value M(row, col) will contribute to procRow and procCol values for fracA and fracB
621  int procRow = ( vecRow[offset] - 1 ) / nBbase;
622  if( procRow >= size ) procRow = size - 1;
623  int procCol = ( vecCol[offset] - 1 ) / nAbase;
624  if( procCol >= size ) procCol = size - 1;
625  int nrInd = tlValRow.get_n();
626  tlValRow.vi_wr[2 * nrInd] = procRow;
627  tlValRow.vi_wr[2 * nrInd + 1] = vecRow[offset] - 1;
628  tlValRow.vr_wr[nrInd] = vecS[offset];
629  tlValRow.inc_n();
630  int ncInd = tlValCol.get_n();
631  tlValCol.vi_wr[3 * ncInd] = procCol;
632  tlValCol.vi_wr[3 * ncInd + 1] = vecRow[offset] - 1;
633  tlValCol.vi_wr[3 * ncInd + 2] = vecCol[offset] - 1; // this is column
634  tlValCol.vr_wr[ncInd] = vecS[offset];
635  tlValCol.inc_n();
636  }
637 
638 #endif
639  offset++;
640  }
641  }
642 #if defined( MOAB_HAVE_MPI )
643  // need to send values for their row and col processors, to compute fractions there
644  // now do the heavy communication
645  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValCol, 0 );
646  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tlValRow, 0 );
647 
648  // we have now, for example, dFracB[ row ] += val ;
649  // so we know that on current task, we received tlValRow
650  // reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
651  // dFracB[ row ] += val ;
652  for( unsigned i = 0; i < tlValRow.get_n(); i++ )
653  {
654  // int fromProc = tlValRow.vi_wr[2 * i];
655  int gRowInd = tlValRow.vi_wr[2 * i + 1];
656  int localIndexRow = gRowInd - nBbase * rank; // modulo nBbase rank is from 0 to size - 1;
657  double wgt = tlValRow.vr_wr[i];
658  assert( localIndexRow >= 0 );
659  assert( nB - localIndexRow > 0 );
660  dFracB[localIndexRow] += wgt;
661  }
662  // to compute dFracA we need vecTargetFaceArea[ row ]; we know the row, and we can get the proc we need it from
663 
664  std::set< int > neededRows;
665  for( unsigned i = 0; i < tlValCol.get_n(); i++ )
666  {
667  int rRowInd = tlValCol.vi_wr[3 * i + 1];
668  neededRows.insert( rRowInd );
669  // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
670  }
671  moab::TupleList tgtAreaReq;
672  tgtAreaReq.initialize( 2, 0, 0, 0, neededRows.size() );
673  tgtAreaReq.enableWriteAccess();
674  for( std::set< int >::iterator sit = neededRows.begin(); sit != neededRows.end(); sit++ )
675  {
676  int neededRow = *sit;
677  int procRow = neededRow / nBbase;
678  if( procRow >= size ) procRow = size - 1;
679  int nr = tgtAreaReq.get_n();
680  tgtAreaReq.vi_wr[2 * nr] = procRow;
681  tgtAreaReq.vi_wr[2 * nr + 1] = neededRow;
682  tgtAreaReq.inc_n();
683  }
684 
685  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaReq, 0 );
686  // we need to send back the tgtArea corresponding to row
687  moab::TupleList tgtAreaInfo; // load it with tgtArea at row
688  tgtAreaInfo.initialize( 2, 0, 0, 1, tgtAreaReq.get_n() );
689  tgtAreaInfo.enableWriteAccess();
690  for( unsigned i = 0; i < tgtAreaReq.get_n(); i++ )
691  {
692  int from_proc = tgtAreaReq.vi_wr[2 * i];
693  int row = tgtAreaReq.vi_wr[2 * i + 1];
694  int locaIndexRow = row - rank * nBbase;
695  double areaToSend = vecTargetFaceArea[locaIndexRow];
696  // int remoteIndex = tgtAreaReq.vi_wr[3*i + 2] ;
697 
698  tgtAreaInfo.vi_wr[2 * i] = from_proc; // send back requested info
699  tgtAreaInfo.vi_wr[2 * i + 1] = row;
700  tgtAreaInfo.vr_wr[i] = areaToSend; // this will be tgt area at row
701  tgtAreaInfo.inc_n();
702  }
703  ( m_pcomm->proc_config().crystal_router() )->gs_transfer( 1, tgtAreaInfo, 0 );
704 
705  std::map< int, double > areaAtRow;
706  for( unsigned i = 0; i < tgtAreaInfo.get_n(); i++ )
707  {
708  // we have received from proc, value for row !
709  int row = tgtAreaInfo.vi_wr[2 * i + 1];
710  areaAtRow[row] = tgtAreaInfo.vr_wr[i];
711  }
712 
713  // we have now for rows the
714  // it is ordered by index, so:
715  // now compute reminder dFracA[ col ] += val / vecSourceFaceArea[ col ] * vecTargetFaceArea[ row ];
716  // tgtAreaInfo will have at index i the area we need (from row)
717  // there should be an easier way :(
718  for( unsigned i = 0; i < tlValCol.get_n(); i++ )
719  {
720  int rRowInd = tlValCol.vi_wr[3 * i + 1];
721  int colInd = tlValCol.vi_wr[3 * i + 2];
722  double val = tlValCol.vr_wr[i];
723  int localColInd = colInd - rank * nAbase; // < local nA
724  // we need vecTargetFaceAreaGlobal[ rRowInd ]; this exists on proc procRow
725  auto itMap = areaAtRow.find( rRowInd ); // it should be different from end
726  if( itMap != areaAtRow.end() )
727  {
728  double areaRow = itMap->second; // we fished a lot for this !
729  dFracA[localColInd] += val / vecSourceFaceArea[localColInd] * areaRow;
730  }
731  }
732 
733 #endif
734  // Load in data
735  NcDim* dimNS = ncMap.add_dim( "n_s", globuf[2] );
736 
737  NcVar* varRow = ncMap.add_var( "row", ncInt, dimNS );
738  NcVar* varCol = ncMap.add_var( "col", ncInt, dimNS );
739  NcVar* varS = ncMap.add_var( "S", ncDouble, dimNS );
740 #ifdef MOAB_HAVE_NETCDFPAR
741  ncMap.enable_var_par_access( varRow, is_independent );
742  ncMap.enable_var_par_access( varCol, is_independent );
743  ncMap.enable_var_par_access( varS, is_independent );
744 #endif
745 
746  varRow->set_cur( (long)offbuf[2] );
747  varRow->put( vecRow, nS );
748 
749  varCol->set_cur( (long)offbuf[2] );
750  varCol->put( vecCol, nS );
751 
752  varS->set_cur( (long)offbuf[2] );
753  varS->put( &( vecS[0] ), nS );
754 
755  // Calculate and write fractional coverage arrays
756  NcVar* varFracA = ncMap.add_var( "frac_a", ncDouble, dimNA );
757 #ifdef MOAB_HAVE_NETCDFPAR
758  ncMap.enable_var_par_access( varFracA, is_independent );
759 #endif
760  varFracA->add_att( "name", "fraction of target coverage of source dof" );
761  varFracA->add_att( "units", "unitless" );
762  varFracA->set_cur( (long)offbuf[0] );
763  varFracA->put( &( dFracA[0] ), nA );
764 
765  NcVar* varFracB = ncMap.add_var( "frac_b", ncDouble, dimNB );
766 #ifdef MOAB_HAVE_NETCDFPAR
767  ncMap.enable_var_par_access( varFracB, is_independent );
768 #endif
769  varFracB->add_att( "name", "fraction of source coverage of target dof" );
770  varFracB->add_att( "units", "unitless" );
771  varFracB->set_cur( (long)offbuf[1] );
772  varFracB->put( &( dFracB[0] ), nB );
773 
774  // Add global attributes
775  // std::map<std::string, std::string>::const_iterator iterAttributes =
776  // mapAttributes.begin();
777  // for (; iterAttributes != mapAttributes.end(); iterAttributes++) {
778  // ncMap.add_att(
779  // iterAttributes->first.c_str(),
780  // iterAttributes->second.c_str());
781  // }
782 
783  ncMap.close();
784 
785  return moab::MB_SUCCESS;
786 }

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

Member Data Documentation

◆ col_dtoc_dofmap

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

Definition at line 520 of file TempestOnlineMap.hpp.

◆ col_gdofmap

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

Definition at line 517 of file TempestOnlineMap.hpp.

Referenced by fill_col_ids(), and LinearRemapNN_MOAB().

◆ colMap

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

Definition at line 522 of file TempestOnlineMap.hpp.

◆ dataGLLNodesDest

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

Definition at line 525 of file TempestOnlineMap.hpp.

◆ dataGLLNodesSrc

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

Definition at line 525 of file TempestOnlineMap.hpp.

◆ dataGLLNodesSrcCov

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

Definition at line 525 of file TempestOnlineMap.hpp.

◆ is_parallel

bool moab::TempestOnlineMap::is_parallel
private

Definition at line 540 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ is_root

bool moab::TempestOnlineMap::is_root
private

Definition at line 540 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_bConserved

bool moab::TempestOnlineMap::m_bConserved
private

Definition at line 532 of file TempestOnlineMap.hpp.

◆ m_destDiscType

DiscretizationType moab::TempestOnlineMap::m_destDiscType
private

Definition at line 526 of file TempestOnlineMap.hpp.

◆ m_dofTagDest

moab::Tag moab::TempestOnlineMap::m_dofTagDest
private

Definition at line 516 of file TempestOnlineMap.hpp.

◆ m_dofTagSrc

moab::Tag moab::TempestOnlineMap::m_dofTagSrc
private

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

Definition at line 516 of file TempestOnlineMap.hpp.

◆ m_eInputType

DiscretizationType moab::TempestOnlineMap::m_eInputType
private

Definition at line 531 of file TempestOnlineMap.hpp.

◆ m_eOutputType

DiscretizationType moab::TempestOnlineMap::m_eOutputType
private

Definition at line 531 of file TempestOnlineMap.hpp.

◆ m_iMonotonicity

int moab::TempestOnlineMap::m_iMonotonicity
private

Definition at line 533 of file TempestOnlineMap.hpp.

◆ m_input_order

int moab::TempestOnlineMap::m_input_order
private

Definition at line 523 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_interface

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

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

Definition at line 504 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshInput

Mesh* moab::TempestOnlineMap::m_meshInput
private

Definition at line 535 of file TempestOnlineMap.hpp.

Referenced by SetMeshInput(), and TempestOnlineMap().

◆ m_meshInputCov

Mesh* moab::TempestOnlineMap::m_meshInputCov
private

Definition at line 536 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshOutput

Mesh* moab::TempestOnlineMap::m_meshOutput
private

Definition at line 537 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshOverlap

Mesh* moab::TempestOnlineMap::m_meshOverlap
private

Definition at line 538 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_nDofsPEl_Dest

int moab::TempestOnlineMap::m_nDofsPEl_Dest
private

Definition at line 530 of file TempestOnlineMap.hpp.

◆ m_nDofsPEl_Src

int moab::TempestOnlineMap::m_nDofsPEl_Src
private

Definition at line 530 of file TempestOnlineMap.hpp.

◆ m_nTotDofs_Dest

int moab::TempestOnlineMap::m_nTotDofs_Dest
private

Definition at line 527 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

◆ m_nTotDofs_Src

int moab::TempestOnlineMap::m_nTotDofs_Src
private

Definition at line 527 of file TempestOnlineMap.hpp.

◆ m_nTotDofs_SrcCov

int moab::TempestOnlineMap::m_nTotDofs_SrcCov
private

Definition at line 527 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

◆ m_output_order

int moab::TempestOnlineMap::m_output_order
private

Definition at line 523 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_remapper

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

The fundamental remapping operator object.

Definition at line 491 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_srcDiscType

DiscretizationType moab::TempestOnlineMap::m_srcDiscType
private

Definition at line 526 of file TempestOnlineMap.hpp.

◆ rank

int moab::TempestOnlineMap::rank
private

Definition at line 541 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ row_dtoc_dofmap

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

Definition at line 520 of file TempestOnlineMap.hpp.

◆ row_gdofmap

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

Definition at line 517 of file TempestOnlineMap.hpp.

Referenced by fill_row_ids(), and LinearRemapNN_MOAB().

◆ rowMap

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

Definition at line 522 of file TempestOnlineMap.hpp.

◆ size

int moab::TempestOnlineMap::size
private

Definition at line 541 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ srccol_dtoc_dofmap

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

Definition at line 520 of file TempestOnlineMap.hpp.

◆ srccol_gdofmap

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

Definition at line 517 of file TempestOnlineMap.hpp.


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