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 , DiscretizationType_CGLL , DiscretizationType_DGLL , DiscretizationType_PCLOUD }
 
typedef double(* sample_function) (double, double)
 

Public Member Functions

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

Private Member Functions

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

Private Attributes

moab::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
 
DataArray3D< int > dataGLLNodesSrc
 
DataArray3D< int > dataGLLNodesSrcCov
 
DataArray3D< int > dataGLLNodesDest
 
DiscretizationType m_srcDiscType
 
DiscretizationType m_destDiscType
 
int m_nTotDofs_Src
 
int m_nTotDofs_SrcCov
 
int m_nTotDofs_Dest
 
int m_nDofsPEl_Src
 
int m_nDofsPEl_Dest
 
DiscretizationType m_eInputType
 
DiscretizationType m_eOutputType
 
bool m_bConserved
 
int m_iMonotonicity
 
Mesh * m_meshInput
 
Mesh * m_meshInputCov
 
Mesh * m_meshOutput
 
Mesh * m_meshOverlap
 
bool is_parallel
 
bool is_root
 
int rank
 
int size
 

Detailed Description

An offline map between two Meshes.

Definition at line 48 of file TempestOnlineMap.hpp.

Member Typedef Documentation

◆ sample_function

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

Definition at line 349 of file TempestOnlineMap.hpp.

Member Enumeration Documentation

◆ DiscretizationType

Enumerator
DiscretizationType_FV 
DiscretizationType_CGLL 
DiscretizationType_DGLL 
DiscretizationType_PCLOUD 

Definition at line 64 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 55 of file TempestOnlineMap.cpp.

55  : OfflineMap(), m_remapper( remapper )
56 {
57  // Get the references for the MOAB core objects
59 #ifdef MOAB_HAVE_MPI
60  m_pcomm = m_remapper->get_parallel_communicator();
61 #endif
62 
63  // Update the references to the meshes
65  m_meshInputCov = remapper->GetCoveringMesh();
68 
69  is_parallel = remapper->is_parallel;
70  is_root = remapper->is_root;
71  rank = remapper->rank;
72  size = remapper->size;
73 
74  // Initialize dimension information from file
75  this->setup_sizes_dimensions();
76 }

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

◆ ~TempestOnlineMap()

moab::TempestOnlineMap::~TempestOnlineMap ( )
virtual

Define a virtual destructor.

Definition at line 103 of file TempestOnlineMap.cpp.

104 {
105  m_interface = NULL;
106 #ifdef MOAB_HAVE_MPI
107  m_pcomm = NULL;
108 #endif
109  m_meshInput = NULL;
110  m_meshOutput = NULL;
111  m_meshOverlap = NULL;
112 }

Member Function Documentation

◆ ApplyWeights() [1/2]

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

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

Definition at line 1676 of file TempestOnlineMap.cpp.

1679 {
1680  moab::ErrorCode rval;
1681 
1682  std::vector< double > solSTagVals;
1683  std::vector< double > solTTagVals;
1684 
1685  moab::Range sents, tents;
1687  {
1689  {
1691  solSTagVals.resize( covSrcEnts.size(), -1.0 );
1692  sents = covSrcEnts;
1693  }
1694  else
1695  {
1697  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1698  -1.0 );
1699  sents = covSrcEnts;
1700  }
1702  {
1704  solTTagVals.resize( tgtEnts.size(), -1.0 );
1705  tents = tgtEnts;
1706  }
1707  else
1708  {
1710  solTTagVals.resize(
1711  tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
1712  tents = tgtEnts;
1713  }
1714  }
1715  else
1716  {
1719  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1720  -1.0 );
1721  solTTagVals.resize(
1722  tgtEnts.size() * this->GetDestinationNDofsPerElement() * this->GetDestinationNDofsPerElement(), -1.0 );
1723 
1724  sents = covSrcEnts;
1725  tents = tgtEnts;
1726  }
1727 
1728  // The tag data is np*np*n_el_src
1729  rval = m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] );MB_CHK_SET_ERR( rval, "Getting local tag data failed" );
1730 
1731  // Compute the application of weights on the suorce solution data and store it in the
1732  // destination solution vector data Optionally, can also perform the transpose application of
1733  // the weight matrix. Set the 3rd argument to true if this is needed
1734  rval = this->ApplyWeights( solSTagVals, solTTagVals, transpose );MB_CHK_SET_ERR( rval, "Applying remap operator onto source vector data failed" );
1735 
1736  // The tag data is np*np*n_el_dest
1737  rval = m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] );MB_CHK_SET_ERR( rval, "Setting local tag data failed" );
1738 
1739  return moab::MB_SUCCESS;
1740 }

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

◆ ApplyWeights() [2/2]

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

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

Referenced by main().

◆ ComputeMetrics()

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

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

Definition at line 2104 of file TempestOnlineMap.cpp.

2109 {
2110  moab::ErrorCode rval;
2111  const bool outputEnabled = ( is_root );
2112  int discOrder;
2113  DiscretizationType discMethod;
2114  moab::EntityHandle meshset;
2116  Mesh* trmesh;
2117  switch( ctx )
2118  {
2119  case Remapper::SourceMesh:
2120  meshset = m_remapper->m_covering_source_set;
2121  trmesh = m_remapper->m_covering_source;
2124  discOrder = m_nDofsPEl_Src;
2125  discMethod = m_eInputType;
2126  break;
2127 
2128  case Remapper::TargetMesh:
2129  meshset = m_remapper->m_target_set;
2130  trmesh = m_remapper->m_target;
2131  entities =
2133  discOrder = m_nDofsPEl_Dest;
2134  discMethod = m_eOutputType;
2135  break;
2136 
2137  default:
2138  if( outputEnabled )
2139  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
2140  return moab::MB_FAILURE;
2141  }
2142 
2143  // Let us create teh solution tag with appropriate information for name, discretization order
2144  // (DoF space)
2145  std::string exactTagName, projTagName;
2146  const int ntotsize = entities.size() * discOrder * discOrder;
2147  int ntotsize_glob = 0;
2148  std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
2149  rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
2150  rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
2151  rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
2152  rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
2153 
2154  std::vector< double > errnorms( 3, 0.0 ), globerrnorms( 3, 0.0 ); // L1Err, L2Err, LinfErr
2155  for( int i = 0; i < ntotsize; ++i )
2156  {
2157  const double error = fabs( exactSolution[i] - projSolution[i] );
2158  errnorms[0] += error;
2159  errnorms[1] += error * error;
2160  errnorms[2] = ( error > errnorms[2] ? error : errnorms[2] );
2161  }
2162 #ifdef MOAB_HAVE_MPI
2163  if( m_pcomm )
2164  {
2165  MPI_Reduce( &ntotsize, &ntotsize_glob, 1, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
2166  MPI_Reduce( &errnorms[0], &globerrnorms[0], 2, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
2167  MPI_Reduce( &errnorms[2], &globerrnorms[2], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
2168  }
2169 #else
2170  ntotsize_glob = ntotsize;
2171  globerrnorms = errnorms;
2172 #endif
2173  globerrnorms[0] = ( globerrnorms[0] / ntotsize_glob );
2174  globerrnorms[1] = std::sqrt( globerrnorms[1] / ntotsize_glob );
2175 
2176  metrics.clear();
2177  metrics["L1Error"] = globerrnorms[0];
2178  metrics["L2Error"] = globerrnorms[1];
2179  metrics["LinfError"] = globerrnorms[2];
2180 
2181  if( verbose && is_root )
2182  {
2183  std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
2184  std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
2185  std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
2186  std::cout << "\t L_inf error = " << globerrnorms[2] << std::endl;
2187  }
2188 
2189  return moab::MB_SUCCESS;
2190 }

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

Referenced by main().

◆ DefineAnalyticalSolution()

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

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

Definition at line 1742 of file TempestOnlineMap.cpp.

1748 {
1749  moab::ErrorCode rval;
1750  const bool outputEnabled = ( is_root );
1751  int discOrder;
1752  DiscretizationType discMethod;
1753  moab::EntityHandle meshset;
1755  Mesh* trmesh;
1756  switch( ctx )
1757  {
1758  case Remapper::SourceMesh:
1759  meshset = m_remapper->m_covering_source_set;
1760  trmesh = m_remapper->m_covering_source;
1763  discOrder = m_nDofsPEl_Src;
1764  discMethod = m_eInputType;
1765  break;
1766 
1767  case Remapper::TargetMesh:
1768  meshset = m_remapper->m_target_set;
1769  trmesh = m_remapper->m_target;
1770  entities =
1772  discOrder = m_nDofsPEl_Dest;
1773  discMethod = m_eOutputType;
1774  break;
1775 
1776  default:
1777  if( outputEnabled )
1778  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
1779  return moab::MB_FAILURE;
1780  }
1781 
1782  // Let us create teh solution tag with appropriate information for name, discretization order
1783  // (DoF space)
1784  rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
1785  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1786  if( clonedSolnTag != NULL )
1787  {
1788  if( cloneSolnName.size() == 0 )
1789  {
1790  cloneSolnName = solnName + std::string( "Cloned" );
1791  }
1792  rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
1793  *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1794  }
1795 
1796  // Triangular quadrature rule
1797  const int TriQuadratureOrder = 10;
1798 
1799  if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
1800 
1801  TriangularQuadratureRule triquadrule( TriQuadratureOrder );
1802 
1803  const int TriQuadraturePoints = triquadrule.GetPoints();
1804 
1805  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1806  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1807 
1808  // Output data
1809  DataArray1D< double > dVar;
1810  DataArray1D< double > dVarMB; // re-arranged local MOAB vector
1811 
1812  // Nodal geometric area
1813  DataArray1D< double > dNodeArea;
1814 
1815  // Calculate element areas
1816  // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
1817 
1818  if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
1819  {
1820  /* Get the spectral points and sample the functionals accordingly */
1821  const bool fGLL = true;
1822  const bool fGLLIntegrate = false;
1823 
1824  // Generate grid metadata
1825  DataArray3D< int > dataGLLNodes;
1826  DataArray3D< double > dataGLLJacobian;
1827 
1828  GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
1829 
1830  // Number of elements
1831  int nElements = trmesh->faces.size();
1832 
1833  // Verify all elements are quadrilaterals
1834  for( int k = 0; k < nElements; k++ )
1835  {
1836  const Face& face = trmesh->faces[k];
1837 
1838  if( face.edges.size() != 4 )
1839  {
1840  _EXCEPTIONT( "Non-quadrilateral face detected; "
1841  "incompatible with --gll" );
1842  }
1843  }
1844 
1845  // Number of unique nodes
1846  int iMaxNode = 0;
1847  for( int i = 0; i < discOrder; i++ )
1848  {
1849  for( int j = 0; j < discOrder; j++ )
1850  {
1851  for( int k = 0; k < nElements; k++ )
1852  {
1853  if( dataGLLNodes[i][j][k] > iMaxNode )
1854  {
1855  iMaxNode = dataGLLNodes[i][j][k];
1856  }
1857  }
1858  }
1859  }
1860 
1861  // Get Gauss-Lobatto quadrature nodes
1862  DataArray1D< double > dG;
1863  DataArray1D< double > dW;
1864 
1865  GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
1866 
1867  // Get Gauss quadrature nodes
1868  const int nGaussP = 10;
1869 
1870  DataArray1D< double > dGaussG;
1871  DataArray1D< double > dGaussW;
1872 
1873  GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
1874 
1875  // Allocate data
1876  dVar.Allocate( iMaxNode );
1877  dVarMB.Allocate( discOrder * discOrder * nElements );
1878  dNodeArea.Allocate( iMaxNode );
1879 
1880  // Sample data
1881  for( int k = 0; k < nElements; k++ )
1882  {
1883  const Face& face = trmesh->faces[k];
1884 
1885  // Sample data at GLL nodes
1886  if( fGLL )
1887  {
1888  for( int i = 0; i < discOrder; i++ )
1889  {
1890  for( int j = 0; j < discOrder; j++ )
1891  {
1892 
1893  // Apply local map
1894  Node node;
1895  Node dDx1G;
1896  Node dDx2G;
1897 
1898  ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
1899 
1900  // Sample data at this point
1901  double dNodeLon = atan2( node.y, node.x );
1902  if( dNodeLon < 0.0 )
1903  {
1904  dNodeLon += 2.0 * M_PI;
1905  }
1906  double dNodeLat = asin( node.z );
1907 
1908  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1909 
1910  dVar[dataGLLNodes[j][i][k] - 1] = dSample;
1911  }
1912  }
1913  // High-order Gaussian integration over basis function
1914  }
1915  else
1916  {
1917  DataArray2D< double > dCoeff( discOrder, discOrder );
1918 
1919  for( int p = 0; p < nGaussP; p++ )
1920  {
1921  for( int q = 0; q < nGaussP; q++ )
1922  {
1923 
1924  // Apply local map
1925  Node node;
1926  Node dDx1G;
1927  Node dDx2G;
1928 
1929  ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
1930 
1931  // Cross product gives local Jacobian
1932  Node nodeCross = CrossProduct( dDx1G, dDx2G );
1933 
1934  double dJacobian =
1935  sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
1936 
1937  // Find components of quadrature point in basis
1938  // of the first Face
1939  SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
1940 
1941  // Sample data at this point
1942  double dNodeLon = atan2( node.y, node.x );
1943  if( dNodeLon < 0.0 )
1944  {
1945  dNodeLon += 2.0 * M_PI;
1946  }
1947  double dNodeLat = asin( node.z );
1948 
1949  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1950 
1951  // Integrate
1952  for( int i = 0; i < discOrder; i++ )
1953  {
1954  for( int j = 0; j < discOrder; j++ )
1955  {
1956 
1957  double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
1958 
1959  dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
1960 
1961  dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
1962  }
1963  }
1964  }
1965  }
1966  }
1967  }
1968 
1969  // Divide by area
1970  if( fGLLIntegrate )
1971  {
1972  for( size_t i = 0; i < dVar.GetRows(); i++ )
1973  {
1974  dVar[i] /= dNodeArea[i];
1975  }
1976  }
1977 
1978  // Let us rearrange the data based on DoF ID specification
1979  if( ctx == Remapper::SourceMesh )
1980  {
1981  for( unsigned j = 0; j < entities.size(); j++ )
1982  for( int p = 0; p < discOrder; p++ )
1983  for( int q = 0; q < discOrder; q++ )
1984  {
1985  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
1986  dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
1987  }
1988  }
1989  else
1990  {
1991  for( unsigned j = 0; j < entities.size(); j++ )
1992  for( int p = 0; p < discOrder; p++ )
1993  for( int q = 0; q < discOrder; q++ )
1994  {
1995  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
1996  dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
1997  }
1998  }
1999 
2000  // Set the tag data
2001  rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
2002  }
2003  else
2004  {
2005  // assert( discOrder == 1 );
2006  if( discMethod == DiscretizationType_FV )
2007  {
2008  /* Compute an element-wise integral to store the sampled solution based on Quadrature
2009  * rules */
2010  // Resize the array
2011  dVar.Allocate( trmesh->faces.size() );
2012 
2013  std::vector< Node >& nodes = trmesh->nodes;
2014 
2015  // Loop through all Faces
2016  for( size_t i = 0; i < trmesh->faces.size(); i++ )
2017  {
2018  const Face& face = trmesh->faces[i];
2019 
2020  // Loop through all sub-triangles
2021  for( size_t j = 0; j < face.edges.size() - 2; j++ )
2022  {
2023 
2024  const Node& node0 = nodes[face[0]];
2025  const Node& node1 = nodes[face[j + 1]];
2026  const Node& node2 = nodes[face[j + 2]];
2027 
2028  // Triangle area
2029  Face faceTri( 3 );
2030  faceTri.SetNode( 0, face[0] );
2031  faceTri.SetNode( 1, face[j + 1] );
2032  faceTri.SetNode( 2, face[j + 2] );
2033 
2034  double dTriangleArea = CalculateFaceArea( faceTri, nodes );
2035 
2036  // Calculate the element average
2037  double dTotalSample = 0.0;
2038 
2039  // Loop through all quadrature points
2040  for( int k = 0; k < TriQuadraturePoints; k++ )
2041  {
2042  Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
2043  TriQuadratureG[k][2] * node2.x,
2044  TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
2045  TriQuadratureG[k][2] * node2.y,
2046  TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
2047  TriQuadratureG[k][2] * node2.z );
2048 
2049  double dMagnitude = node.Magnitude();
2050  node.x /= dMagnitude;
2051  node.y /= dMagnitude;
2052  node.z /= dMagnitude;
2053 
2054  double dLon = atan2( node.y, node.x );
2055  if( dLon < 0.0 )
2056  {
2057  dLon += 2.0 * M_PI;
2058  }
2059  double dLat = asin( node.z );
2060 
2061  double dSample = ( *testFunction )( dLon, dLat );
2062 
2063  dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
2064  }
2065 
2066  dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
2067  }
2068  }
2069  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2070  }
2071  else /* discMethod == DiscretizationType_PCLOUD */
2072  {
2073  /* Get the coordinates of the vertices and sample the functionals accordingly */
2074  std::vector< Node >& nodes = trmesh->nodes;
2075 
2076  // Resize the array
2077  dVar.Allocate( nodes.size() );
2078 
2079  for( size_t j = 0; j < nodes.size(); j++ )
2080  {
2081  Node& node = nodes[j];
2082  double dMagnitude = node.Magnitude();
2083  node.x /= dMagnitude;
2084  node.y /= dMagnitude;
2085  node.z /= dMagnitude;
2086  double dLon = atan2( node.y, node.x );
2087  if( dLon < 0.0 )
2088  {
2089  dLon += 2.0 * M_PI;
2090  }
2091  double dLat = asin( node.z );
2092 
2093  double dSample = ( *testFunction )( dLon, dLat );
2094  dVar[j] = dSample;
2095  }
2096 
2097  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
2098  }
2099  }
2100 
2101  return moab::MB_SUCCESS;
2102 }

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

Referenced by main().

◆ fill_col_ids()

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

Definition at line 384 of file TempestOnlineMap.hpp.

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

References col_gdofmap, and MB_SUCCESS.

◆ fill_row_ids()

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

Definition at line 374 of file TempestOnlineMap.hpp.

375  {
376  ids_of_interest.reserve( row_gdofmap.size() );
377  // need to add 1
378  for( auto it = row_gdofmap.begin(); it != row_gdofmap.end(); it++ )
379  ids_of_interest.push_back( *it + 1 );
380 
381  return moab::MB_SUCCESS;
382  }

References MB_SUCCESS, and row_gdofmap.

◆ GenerateRemappingWeights()

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

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

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

Definition at line 691 of file TempestOnlineMap.cpp.

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

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

Referenced by main().

◆ GetColGlobalDoF()

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

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

Definition at line 481 of file TempestOnlineMap.hpp.

482 {
483  return col_gdofmap[localColID];
484 }

◆ GetDestinationGlobalNDofs()

int moab::TempestOnlineMap::GetDestinationGlobalNDofs ( )

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

◆ GetDestinationLocalNDofs()

int moab::TempestOnlineMap::GetDestinationLocalNDofs ( )

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

◆ GetDestinationNDofsPerElement()

int moab::TempestOnlineMap::GetDestinationNDofsPerElement ( )
inline

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

Definition at line 499 of file TempestOnlineMap.hpp.

500 {
501  return m_nDofsPEl_Dest;
502 }

◆ GetGlobalSourceAreas()

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

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

◆ GetGlobalTargetAreas()

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

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

◆ GetIndexOfColGlobalDoF()

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

Get the index of globaColDoF.

Definition at line 486 of file TempestOnlineMap.hpp.

487 {
488  return globalColDoF + 1; // temporary
489 }

◆ GetIndexOfRowGlobalDoF()

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

Get the index of globaRowDoF.

Definition at line 475 of file TempestOnlineMap.hpp.

476 {
477  return globalRowDoF + 1;
478 }

◆ GetRowGlobalDoF()

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

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

Definition at line 470 of file TempestOnlineMap.hpp.

471 {
472  return row_gdofmap[localRowID];
473 }

◆ GetSourceGlobalNDofs()

int moab::TempestOnlineMap::GetSourceGlobalNDofs ( )

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

◆ GetSourceLocalNDofs()

int moab::TempestOnlineMap::GetSourceLocalNDofs ( )

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

◆ GetSourceNDofsPerElement()

int moab::TempestOnlineMap::GetSourceNDofsPerElement ( )
inline

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

Definition at line 492 of file TempestOnlineMap.hpp.

493 {
494  return m_nDofsPEl_Src;
495 }

◆ IsConservative()

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

Determine if the map is conservative.

Definition at line 1496 of file TempestOnlineMap.cpp.

1497 {
1498 #ifndef MOAB_HAVE_MPI
1499 
1500  return OfflineMap::IsConservative( dTolerance );
1501 
1502 #else
1503  // return OfflineMap::IsConservative(dTolerance);
1504 
1505  int ierr;
1506  // Get map entries
1507  DataArray1D< int > dataRows;
1508  DataArray1D< int > dataCols;
1509  DataArray1D< double > dataEntries;
1510  const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1511  const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1512 
1513  // Calculate column sums
1514  std::vector< int > dColumnsUnique;
1515  std::vector< double > dColumnSums;
1516 
1517  int nColumns = m_mapRemap.GetColumns();
1518  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1519  dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1520  dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1521 
1522  for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
1523  {
1524  dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1525 
1526  assert( dataCols[i] < m_nTotDofs_SrcCov );
1527 
1528  // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
1529  int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1530  // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1531  dColumnsUnique[dataCols[i]] = colGID;
1532 
1533  // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
1534  // std::endl;
1535  }
1536 
1537  int rootProc = 0;
1538  std::vector< int > nElementsInProc;
1539  const int nDATA = 3;
1540  if( !rank ) nElementsInProc.resize( size * nDATA );
1541  int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1542  ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1543  if( ierr != MPI_SUCCESS ) return -1;
1544 
1545  int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
1546  std::vector< int > dColumnIndices;
1547  std::vector< double > dColumnSourceAreas;
1548  std::vector< double > dColumnSumsTotal;
1549  std::vector< int > displs, rcount;
1550  if( rank == rootProc )
1551  {
1552  displs.resize( size + 1, 0 );
1553  rcount.resize( size, 0 );
1554  int gsum = 0;
1555  for( int ir = 0; ir < size; ++ir )
1556  {
1557  nTotVals += nElementsInProc[ir * nDATA];
1558  nTotColumns += nElementsInProc[ir * nDATA + 1];
1559  nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1560 
1561  displs[ir] = gsum;
1562  rcount[ir] = nElementsInProc[ir * nDATA + 1];
1563  gsum += rcount[ir];
1564 
1565  // printf("%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns,
1566  // displs[ir], rcount[ir], gsum);
1567  }
1568 
1569  dColumnIndices.resize( nTotColumns, -1 );
1570  dColumnSumsTotal.resize( nTotColumns, 0.0 );
1571  // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
1572  }
1573 
1574  // Gather all ColumnSums to root process and accumulate
1575  // We expect that the sums of all columns equate to 1.0 within user specified tolerance
1576  // Need to do a gatherv here since different processes have different number of elements
1577  // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
1578  // MPI_SUM, 0, m_pcomm->comm());
1579  ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
1580  displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1581  if( ierr != MPI_SUCCESS ) return -1;
1582  ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
1583  displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1584  if( ierr != MPI_SUCCESS ) return -1;
1585  // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
1586  // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
1587  // MPI_SUCCESS ) return -1;
1588 
1589  // Clean out unwanted arrays now
1590  dColumnSums.clear();
1591  dColumnsUnique.clear();
1592 
1593  // Verify all column sums equal the input Jacobian
1594  int fConservative = 0;
1595  if( rank == rootProc )
1596  {
1597  displs[size] = ( nTotColumns );
1598  // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
1599  std::map< int, double > dColumnSumsOnRoot;
1600  // std::map<int, double> dColumnSourceAreasOnRoot;
1601  for( int ir = 0; ir < size; ir++ )
1602  {
1603  for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1604  {
1605  if( dColumnIndices[ips] < 0 ) continue;
1606  // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
1607  assert( dColumnIndices[ips] < nTotColumnsUnq );
1608  dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
1609  // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
1610  // dColumnSourceAreas[ dColumnIndices[ips] ]
1611  }
1612  }
1613 
1614  for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1615  {
1616  // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
1617  if( fabs( it->second - 1.0 ) > dTolerance )
1618  {
1619  fConservative++;
1620  Announce( "TempestOnlineMap is not conservative in column "
1621  // "%i (%1.15e)", it->first, it->second );
1622  "%i (%1.15e)",
1623  it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
1624  }
1625  }
1626  }
1627 
1628  // TODO: Just do a broadcast from root instead of a reduction
1629  ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1630  if( ierr != MPI_SUCCESS ) return -1;
1631 
1632  return fConservative;
1633 #endif
1634 }

References size.

◆ IsConsistent()

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

Determine if the map is first-order accurate.

Definition at line 1450 of file TempestOnlineMap.cpp.

1451 {
1452 #ifndef MOAB_HAVE_MPI
1453 
1454  return OfflineMap::IsConsistent( dTolerance );
1455 
1456 #else
1457 
1458  // Get map entries
1459  DataArray1D< int > dataRows;
1460  DataArray1D< int > dataCols;
1461  DataArray1D< double > dataEntries;
1462 
1463  // Calculate row sums
1464  DataArray1D< double > dRowSums;
1465  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1466  dRowSums.Allocate( m_mapRemap.GetRows() );
1467 
1468  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1469  {
1470  dRowSums[dataRows[i]] += dataEntries[i];
1471  }
1472 
1473  // Verify all row sums are equal to 1
1474  int fConsistent = 0;
1475  for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1476  {
1477  if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1478  {
1479  fConsistent++;
1480  int rowGID = row_gdofmap[i];
1481  Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1482  }
1483  }
1484 
1485  int ierr;
1486  int fConsistentGlobal = 0;
1487  ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1488  if( ierr != MPI_SUCCESS ) return -1;
1489 
1490  return fConsistentGlobal;
1491 #endif
1492 }

◆ IsMonotone()

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

Determine if the map is monotone.

Definition at line 1638 of file TempestOnlineMap.cpp.

1639 {
1640 #ifndef MOAB_HAVE_MPI
1641 
1642  return OfflineMap::IsMonotone( dTolerance );
1643 
1644 #else
1645 
1646  // Get map entries
1647  DataArray1D< int > dataRows;
1648  DataArray1D< int > dataCols;
1649  DataArray1D< double > dataEntries;
1650 
1651  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1652 
1653  // Verify all entries are in the range [0,1]
1654  int fMonotone = 0;
1655  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1656  {
1657  if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1658  {
1659  fMonotone++;
1660 
1661  Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1662  }
1663  }
1664 
1665  int ierr;
1666  int fMonotoneGlobal = 0;
1667  ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1668  if( ierr != MPI_SUCCESS ) return -1;
1669 
1670  return fMonotoneGlobal;
1671 #endif
1672 }

◆ LinearRemapFVtoFV_Tempest_MOAB()

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

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

Definition at line 115 of file TempestLinearRemap.cpp.

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

References dbgprint.

◆ LinearRemapFVtoGLL_MOAB()

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

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

◆ LinearRemapGLLtoGLL2_MOAB()

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

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

Definition at line 895 of file TempestLinearRemap.cpp.

906 {
907  // Triangular quadrature rule
908  TriangularQuadratureRule triquadrule( 8 );
909 
910  const DataArray2D< double >& dG = triquadrule.GetG();
911  const DataArray1D< double >& dW = triquadrule.GetW();
912 
913  // Get SparseMatrix represntation of the OfflineMap
914  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
915 
916  // Sample coefficients
917  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
918  DataArray2D< double > dSampleCoeffOut( nPout, nPout );
919 
920  // Announcemnets
921  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
922  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
923  if( is_root )
924  {
925  dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
926  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
927  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
928  }
929 
930  // Build the integration array for each element on m_meshOverlap
931  DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
932 
933  // Number of overlap Faces per source Face
934  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
935 
936  int ixOverlap = 0;
937  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
938  {
939  // Determine how many overlap Faces and triangles are present
940  int nOverlapFaces = 0;
941  size_t ixOverlapTemp = ixOverlap;
942  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
943  {
944  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
945  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
946  {
947  break;
948  }
949 
950  nOverlapFaces++;
951  }
952 
953  nAllOverlapFaces[ixFirst] = nOverlapFaces;
954 
955  // Increment the current overlap index
956  ixOverlap += nAllOverlapFaces[ixFirst];
957  }
958 
959  // Geometric area of each output node
960  DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
961 
962  // Area of each overlap element in the output basis
963  DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
964 
965  // Loop through all faces on m_meshInput
966  ixOverlap = 0;
967 #ifdef VERBOSE
968  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
969 #endif
970  if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
971 
972  // generic triangle used for area computation, for triangles around the center of overlap face;
973  // used for overlap faces with more than 4 edges;
974  // nodes array will be set for each triangle;
975  // these triangles are not part of the mesh structure, they are just temporary during
976  // aforementioned decomposition.
977  Face faceTri( 3 );
978  NodeVector nodes( 3 );
979  faceTri.SetNode( 0, 0 );
980  faceTri.SetNode( 1, 1 );
981  faceTri.SetNode( 2, 2 );
982 
983  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
984  {
985 #ifdef VERBOSE
986  // Announce computation progress
987  if( ixFirst % outputFrequency == 0 && is_root )
988  {
989  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
990  }
991 #endif
992  // Quantities from the First Mesh
993  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
994 
995  const NodeVector& nodesFirst = m_meshInputCov->nodes;
996 
997  // Number of overlapping Faces and triangles
998  int nOverlapFaces = nAllOverlapFaces[ixFirst];
999 
1000  if( !nOverlapFaces ) continue;
1001 
1002  // // Calculate total element Jacobian
1003  // double dTotalJacobian = 0.0;
1004  // for (int s = 0; s < nPin; s++) {
1005  // for (int t = 0; t < nPin; t++) {
1006  // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
1007  // }
1008  // }
1009 
1010  // Loop through all Overlap Faces
1011  for( int i = 0; i < nOverlapFaces; i++ )
1012  {
1013  // Quantities from the overlap Mesh
1014  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
1015 
1016  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1017 
1018  // Quantities from the Second Mesh
1019  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1020 
1021  // signal to not participate, because it is a ghost target
1022  if( ixSecond < 0 ) continue; // do not do anything
1023 
1024  const NodeVector& nodesSecond = m_meshOutput->nodes;
1025 
1026  const Face& faceSecond = m_meshOutput->faces[ixSecond];
1027 
1028  int nbEdges = faceOverlap.edges.size();
1029  int nOverlapTriangles = 1;
1030  Node center; // not used if nbEdges == 3
1031  if( nbEdges > 3 )
1032  { // decompose from center in this case
1033  nOverlapTriangles = nbEdges;
1034  for( int k = 0; k < nbEdges; k++ )
1035  {
1036  const Node& node = nodesOverlap[faceOverlap[k]];
1037  center = center + node;
1038  }
1039  center = center / nbEdges;
1040  center = center.Normalized(); // project back on sphere of radius 1
1041  }
1042 
1043  Node node0, node1, node2;
1044  double dTriArea;
1045 
1046  // Loop over all sub-triangles of this Overlap Face
1047  for( int j = 0; j < nOverlapTriangles; j++ )
1048  {
1049  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1050  {
1051  node0 = nodesOverlap[faceOverlap[0]];
1052  node1 = nodesOverlap[faceOverlap[1]];
1053  node2 = nodesOverlap[faceOverlap[2]];
1054  dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1055  }
1056  else // decompose polygon in triangles around the center
1057  {
1058  node0 = center;
1059  node1 = nodesOverlap[faceOverlap[j]];
1060  int j1 = ( j + 1 ) % nbEdges;
1061  node2 = nodesOverlap[faceOverlap[j1]];
1062  nodes[0] = center;
1063  nodes[1] = node1;
1064  nodes[2] = node2;
1065  dTriArea = CalculateFaceArea( faceTri, nodes );
1066  }
1067 
1068  for( int k = 0; k < triquadrule.GetPoints(); k++ )
1069  {
1070  // Get the nodal location of this point
1071  double dX[3];
1072 
1073  dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
1074  dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
1075  dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
1076 
1077  double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
1078 
1079  dX[0] /= dMag;
1080  dX[1] /= dMag;
1081  dX[2] /= dMag;
1082 
1083  Node nodeQuadrature( dX[0], dX[1], dX[2] );
1084 
1085  // Find the components of this quadrature point in the basis
1086  // of the first Face.
1087  double dAlphaIn;
1088  double dBetaIn;
1089 
1090  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
1091 
1092  // Find the components of this quadrature point in the basis
1093  // of the second Face.
1094  double dAlphaOut;
1095  double dBetaOut;
1096 
1097  ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
1098 
1099  /*
1100  // Check inverse map value
1101  if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
1102  (dBetaIn < 0.0) || (dBetaIn > 1.0)
1103  ) {
1104  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1105  dAlphaIn, dBetaIn);
1106  }
1107 
1108  // Check inverse map value
1109  if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
1110  (dBetaOut < 0.0) || (dBetaOut > 1.0)
1111  ) {
1112  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1113  dAlphaOut, dBetaOut);
1114  }
1115  */
1116  // Sample the First finite element at this point
1117  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1118 
1119  // Sample the Second finite element at this point
1120  SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
1121 
1122  // Overlap output area
1123  for( int s = 0; s < nPout; s++ )
1124  {
1125  for( int t = 0; t < nPout; t++ )
1126  {
1127  double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
1128 
1129  dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
1130 
1131  dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
1132  }
1133  }
1134 
1135  // Compute overlap integral
1136  int ixp = 0;
1137  for( int p = 0; p < nPin; p++ )
1138  {
1139  for( int q = 0; q < nPin; q++ )
1140  {
1141  int ixs = 0;
1142  for( int s = 0; s < nPout; s++ )
1143  {
1144  for( int t = 0; t < nPout; t++ )
1145  {
1146  // Sample the Second finite element at this point
1147  dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
1148  dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
1149 
1150  ixs++;
1151  }
1152  }
1153 
1154  ixp++;
1155  }
1156  }
1157  }
1158  }
1159  }
1160 
1161  // Coefficients
1162  DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
1163 
1164  for( int i = 0; i < nOverlapFaces; i++ )
1165  {
1166  // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1167 
1168  int ixp = 0;
1169  for( int p = 0; p < nPin; p++ )
1170  {
1171  for( int q = 0; q < nPin; q++ )
1172  {
1173  int ixs = 0;
1174  for( int s = 0; s < nPout; s++ )
1175  {
1176  for( int t = 0; t < nPout; t++ )
1177  {
1178  dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
1179  dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1180 
1181  ixs++;
1182  }
1183  }
1184 
1185  ixp++;
1186  }
1187  }
1188  }
1189 
1190  // Source areas
1191  DataArray1D< double > vecSourceArea( nPin * nPin );
1192 
1193  for( int p = 0; p < nPin; p++ )
1194  {
1195  for( int q = 0; q < nPin; q++ )
1196  {
1197  vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
1198  }
1199  }
1200 
1201  // Target areas
1202  DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
1203 
1204  for( int i = 0; i < nOverlapFaces; i++ )
1205  {
1206  // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1207  int ixs = 0;
1208  for( int s = 0; s < nPout; s++ )
1209  {
1210  for( int t = 0; t < nPout; t++ )
1211  {
1212  vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
1213 
1214  ixs++;
1215  }
1216  }
1217  }
1218 
1219  // Force consistency and conservation
1220  if( !fNoConservation )
1221  {
1222  ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
1223  }
1224 
1225  // Update global coefficients
1226  for( int i = 0; i < nOverlapFaces; i++ )
1227  {
1228  int ixp = 0;
1229  for( int p = 0; p < nPin; p++ )
1230  {
1231  for( int q = 0; q < nPin; q++ )
1232  {
1233  int ixs = 0;
1234  for( int s = 0; s < nPout; s++ )
1235  {
1236  for( int t = 0; t < nPout; t++ )
1237  {
1238  dGlobalIntArray[ixp][ixOverlap + i][ixs] =
1239  dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1240 
1241  ixs++;
1242  }
1243  }
1244 
1245  ixp++;
1246  }
1247  }
1248  }
1249 
1250 #ifdef VVERBOSE
1251  // Check column sums (conservation)
1252  for( int i = 0; i < nPin * nPin; i++ )
1253  {
1254  double dColSum = 0.0;
1255  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1256  {
1257  dColSum += dCoeff[j][i] * vecTargetArea[j];
1258  }
1259  printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
1260  }
1261 
1262  // Check row sums (consistency)
1263  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1264  {
1265  double dRowSum = 0.0;
1266  for( int i = 0; i < nPin * nPin; i++ )
1267  {
1268  dRowSum += dCoeff[j][i];
1269  }
1270  printf( "Row %i: %1.15e\n", j, dRowSum );
1271  }
1272 #endif
1273 
1274  // Increment the current overlap index
1275  ixOverlap += nOverlapFaces;
1276  }
1277 
1278  // Build redistribution map within target element
1279  if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
1280  DataArray1D< double > dRedistSourceArea( nPout * nPout );
1281  DataArray1D< double > dRedistTargetArea( nPout * nPout );
1282  std::vector< DataArray2D< double > > dRedistributionMaps;
1283  dRedistributionMaps.resize( m_meshOutput->faces.size() );
1284 
1285  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1286  {
1287  dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
1288 
1289  for( int i = 0; i < nPout * nPout; i++ )
1290  {
1291  dRedistributionMaps[ixSecond][i][i] = 1.0;
1292  }
1293 
1294  for( int s = 0; s < nPout * nPout; s++ )
1295  {
1296  dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
1297  }
1298 
1299  for( int s = 0; s < nPout * nPout; s++ )
1300  {
1301  dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
1302  }
1303 
1304  if( !fNoConservation )
1305  {
1306  ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
1307  ( nMonotoneType != 0 ) );
1308 
1309  for( int s = 0; s < nPout * nPout; s++ )
1310  {
1311  for( int t = 0; t < nPout * nPout; t++ )
1312  {
1313  dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
1314  }
1315  }
1316  }
1317  }
1318 
1319  // Construct the total geometric area
1320  DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
1321  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1322  {
1323  for( int s = 0; s < nPout; s++ )
1324  {
1325  for( int t = 0; t < nPout; t++ )
1326  {
1327  dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
1328  dGeometricOutputArea[ixSecond][s * nPout + t];
1329  }
1330  }
1331  }
1332 
1333  // Compose the integration operator with the output map
1334  ixOverlap = 0;
1335 
1336  if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
1337 
1338  // Map from source DOFs to target DOFs with redistribution applied
1339  DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
1340 
1341  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1342  {
1343 #ifdef VERBOSE
1344  // Announce computation progress
1345  if( ixFirst % outputFrequency == 0 && is_root )
1346  {
1347  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1348  }
1349 #endif
1350  // Number of overlapping Faces and triangles
1351  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1352 
1353  if( !nOverlapFaces ) continue;
1354 
1355  // Put composed array into map
1356  for( int j = 0; j < nOverlapFaces; j++ )
1357  {
1358  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1359 
1360  // signal to not participate, because it is a ghost target
1361  if( ixSecondFace < 0 ) continue; // do not do anything
1362 
1363  dRedistributedOp.Zero();
1364  for( int p = 0; p < nPin * nPin; p++ )
1365  {
1366  for( int s = 0; s < nPout * nPout; s++ )
1367  {
1368  for( int t = 0; t < nPout * nPout; t++ )
1369  {
1370  dRedistributedOp[p][s] +=
1371  dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
1372  }
1373  }
1374  }
1375 
1376  int ixp = 0;
1377  for( int p = 0; p < nPin; p++ )
1378  {
1379  for( int q = 0; q < nPin; q++ )
1380  {
1381  int ixFirstNode;
1382  if( fContinuousIn )
1383  {
1384  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
1385  }
1386  else
1387  {
1388  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
1389  }
1390 
1391  int ixs = 0;
1392  for( int s = 0; s < nPout; s++ )
1393  {
1394  for( int t = 0; t < nPout; t++ )
1395  {
1396  int ixSecondNode;
1397  if( fContinuousOut )
1398  {
1399  ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
1400 
1401  if( !fNoConservation )
1402  {
1403  smatMap( ixSecondNode, ixFirstNode ) +=
1404  dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
1405  }
1406  else
1407  {
1408  smatMap( ixSecondNode, ixFirstNode ) +=
1409  dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
1410  }
1411  }
1412  else
1413  {
1414  ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
1415 
1416  if( !fNoConservation )
1417  {
1418  smatMap( ixSecondNode, ixFirstNode ) +=
1419  dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
1420  }
1421  else
1422  {
1423  smatMap( ixSecondNode, ixFirstNode ) +=
1424  dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
1425  }
1426  }
1427 
1428  ixs++;
1429  }
1430  }
1431 
1432  ixp++;
1433  }
1434  }
1435  }
1436 
1437  // Increment the current overlap index
1438  ixOverlap += nOverlapFaces;
1439  }
1440 
1441  return;
1442 }

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 1446 of file TempestLinearRemap.cpp.

1456 {
1457  // Gauss-Lobatto quadrature within Faces
1458  DataArray1D< double > dGL;
1459  DataArray1D< double > dWL;
1460 
1461  GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
1462 
1463  // Get SparseMatrix represntation of the OfflineMap
1464  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1465 
1466  // Sample coefficients
1467  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1468 
1469  // Announcemnets
1470  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1471  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
1472  if( is_root )
1473  {
1474  dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
1475  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1476  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1477  }
1478 
1479  // Number of overlap Faces per source Face
1480  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1481 
1482  int ixOverlap = 0;
1483 
1484  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1485  {
1486  size_t ixOverlapTemp = ixOverlap;
1487  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1488  {
1489  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1490 
1491  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1492 
1493  nAllOverlapFaces[ixFirst]++;
1494  }
1495 
1496  // Increment the current overlap index
1497  ixOverlap += nAllOverlapFaces[ixFirst];
1498  }
1499 
1500  // Number of times this point was found
1501  DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
1502 
1503  ixOverlap = 0;
1504 #ifdef VERBOSE
1505  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1506 #endif
1507  // Loop through all faces on m_meshInputCov
1508  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1509  {
1510 #ifdef VERBOSE
1511  // Announce computation progress
1512  if( ixFirst % outputFrequency == 0 && is_root )
1513  {
1514  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1515  }
1516 #endif
1517  // Quantities from the First Mesh
1518  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1519 
1520  const NodeVector& nodesFirst = m_meshInputCov->nodes;
1521 
1522  // Number of overlapping Faces and triangles
1523  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1524 
1525  // Loop through all Overlap Faces
1526  for( int i = 0; i < nOverlapFaces; i++ )
1527  {
1528  // Quantities from the Second Mesh
1529  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1530 
1531  // signal to not participate, because it is a ghost target
1532  if( ixSecond < 0 ) continue; // do not do anything
1533 
1534  const NodeVector& nodesSecond = m_meshOutput->nodes;
1535  const Face& faceSecond = m_meshOutput->faces[ixSecond];
1536 
1537  // Loop through all nodes on the second face
1538  for( int s = 0; s < nPout; s++ )
1539  {
1540  for( int t = 0; t < nPout; t++ )
1541  {
1542  size_t ixSecondNode;
1543  if( fContinuousOut )
1544  {
1545  ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
1546  }
1547  else
1548  {
1549  ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
1550  }
1551 
1552  if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
1553 
1554  // Check if this node has been found already
1555  if( fSecondNodeFound[ixSecondNode] ) continue;
1556 
1557  // Check this node
1558  Node node;
1559  Node dDx1G;
1560  Node dDx2G;
1561 
1562  ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
1563 
1564  // Find the components of this quadrature point in the basis
1565  // of the first Face.
1566  double dAlphaIn;
1567  double dBetaIn;
1568 
1569  ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
1570 
1571  // Check if this node is within the first Face
1572  if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
1573  ( dBetaIn > 1.0 + 1.0e-10 ) )
1574  continue;
1575 
1576  // Node is within the overlap region, mark as found
1577  fSecondNodeFound[ixSecondNode] = true;
1578 
1579  // Sample the First finite element at this point
1580  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1581 
1582  // Add to map
1583  for( int p = 0; p < nPin; p++ )
1584  {
1585  for( int q = 0; q < nPin; q++ )
1586  {
1587  int ixFirstNode;
1588  if( fContinuousIn )
1589  {
1590  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
1591  }
1592  else
1593  {
1594  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
1595  }
1596 
1597  smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
1598  }
1599  }
1600  }
1601  }
1602  }
1603 
1604  // Increment the current overlap index
1605  ixOverlap += nOverlapFaces;
1606  }
1607 
1608  // Check for missing samples
1609  for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
1610  {
1611  if( !fSecondNodeFound[i] )
1612  {
1613  _EXCEPTION1( "Can't sample point %i", i );
1614  }
1615  }
1616 
1617  return;
1618 }

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 52 of file TempestLinearRemap.cpp.

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

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

◆ LinearRemapSE0_Tempest_MOAB()

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

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

◆ LinearRemapSE4_Tempest_MOAB()

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

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

Definition at line 469 of file TempestLinearRemap.cpp.

474 {
475  // Order of the polynomial interpolant
476  int nP = dataGLLNodes.GetRows();
477 
478  // Order of triangular quadrature rule
479  const int TriQuadRuleOrder = 4;
480 
481  // Triangular quadrature rule
482  TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
483 
484  int TriQuadraturePoints = triquadrule.GetPoints();
485 
486  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
487 
488  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
489 
490  // Sample coefficients
491  DataArray2D< double > dSampleCoeff( nP, nP );
492 
493  // GLL Quadrature nodes on quadrilateral elements
494  DataArray1D< double > dG;
495  DataArray1D< double > dW;
496  GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
497 
498  // Announcements
499  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
500  dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
501  if( is_root )
502  {
503  dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
504  dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
505  dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
506  }
507 
508  // Get SparseMatrix represntation of the OfflineMap
509  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
510 
511  // NodeVector from m_meshOverlap
512  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
513  const NodeVector& nodesFirst = m_meshInputCov->nodes;
514 
515  // Vector of source areas
516  DataArray1D< double > vecSourceArea( nP * nP );
517 
518  DataArray1D< double > vecTargetArea;
519  DataArray2D< double > dCoeff;
520 
521 #ifdef VERBOSE
522  std::stringstream sstr;
523  sstr << "remapdata_" << rank << ".txt";
524  std::ofstream output_file( sstr.str() );
525 #endif
526 
527  // Current Overlap Face
528  int ixOverlap = 0;
529 #ifdef VERBOSE
530  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
531 #endif
532  // generic triangle used for area computation, for triangles around the center of overlap face;
533  // used for overlap faces with more than 4 edges;
534  // nodes array will be set for each triangle;
535  // these triangles are not part of the mesh structure, they are just temporary during
536  // aforementioned decomposition.
537  Face faceTri( 3 );
538  NodeVector nodes( 3 );
539  faceTri.SetNode( 0, 0 );
540  faceTri.SetNode( 1, 1 );
541  faceTri.SetNode( 2, 2 );
542 
543  // Loop over all input Faces
544  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
545  {
546  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
547 
548  if( faceFirst.edges.size() != 4 )
549  {
550  _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
551  }
552 #ifdef VERBOSE
553  // Announce computation progress
554  if( ixFirst % outputFrequency == 0 && is_root )
555  {
556  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
557  }
558 #endif
559  // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
560  // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
561  // However, the relation with MOAB and Tempest will go out of the roof
562 
563  // Determine how many overlap Faces and triangles are present
564  int nOverlapFaces = 0;
565  size_t ixOverlapTemp = ixOverlap;
566  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
567  {
568  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
569  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
570  {
571  break;
572  }
573 
574  nOverlapFaces++;
575  }
576 
577  // No overlaps
578  if( nOverlapFaces == 0 )
579  {
580  continue;
581  }
582 
583  // Allocate remap coefficients array for meshFirst Face
584  DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
585 
586  // Find the local remap coefficients
587  for( int j = 0; j < nOverlapFaces; j++ )
588  {
589  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
590  if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
591  {
592  Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
593  m_meshOverlap->vecFaceArea[ixOverlap + j] );
594  int n = faceOverlap.edges.size();
595  Announce( "Number nodes: %d", n );
596  for( int k = 0; k < n; k++ )
597  {
598  Node nd = nodesOverlap[faceOverlap[k]];
599  Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
600  }
601  continue;
602  }
603 
604  // #ifdef VERBOSE
605  // if ( is_root )
606  // Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
607  // m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
608  // m_meshOverlap->vecFaceArea[ixOverlap + j] );
609  // #endif
610 
611  int nbEdges = faceOverlap.edges.size();
612  int nOverlapTriangles = 1;
613  Node center; // not used if nbEdges == 3
614  if( nbEdges > 3 )
615  { // decompose from center in this case
616  nOverlapTriangles = nbEdges;
617  for( int k = 0; k < nbEdges; k++ )
618  {
619  const Node& node = nodesOverlap[faceOverlap[k]];
620  center = center + node;
621  }
622  center = center / nbEdges;
623  center = center.Normalized(); // project back on sphere of radius 1
624  }
625 
626  Node node0, node1, node2;
627  double dTriangleArea;
628 
629  // Loop over all sub-triangles of this Overlap Face
630  for( int k = 0; k < nOverlapTriangles; k++ )
631  {
632  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
633  {
634  node0 = nodesOverlap[faceOverlap[0]];
635  node1 = nodesOverlap[faceOverlap[1]];
636  node2 = nodesOverlap[faceOverlap[2]];
637  dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
638  }
639  else // decompose polygon in triangles around the center
640  {
641  node0 = center;
642  node1 = nodesOverlap[faceOverlap[k]];
643  int k1 = ( k + 1 ) % nbEdges;
644  node2 = nodesOverlap[faceOverlap[k1]];
645  nodes[0] = center;
646  nodes[1] = node1;
647  nodes[2] = node2;
648  dTriangleArea = CalculateFaceArea( faceTri, nodes );
649  }
650  // Coordinates of quadrature Node
651  for( int l = 0; l < TriQuadraturePoints; l++ )
652  {
653  Node nodeQuadrature;
654  nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
655  TriQuadratureG[l][2] * node2.x;
656 
657  nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
658  TriQuadratureG[l][2] * node2.y;
659 
660  nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
661  TriQuadratureG[l][2] * node2.z;
662 
663  nodeQuadrature = nodeQuadrature.Normalized();
664 
665  // Find components of quadrature point in basis
666  // of the first Face
667  double dAlpha;
668  double dBeta;
669 
670  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
671 
672  // Check inverse map value
673  if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
674  ( dBeta > 1.0 + 1.0e-13 ) )
675  {
676  _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
677  "(%1.5e %1.5e)",
678  j, l, dAlpha, dBeta );
679  }
680 
681  // Sample the finite element at this point
682  SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
683 
684  // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
685  for( int p = 0; p < nP; p++ )
686  {
687  for( int q = 0; q < nP; q++ )
688  {
689  dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
690  m_meshOverlap->vecFaceArea[ixOverlap + j];
691  }
692  }
693  }
694  }
695  }
696 
697 #ifdef VERBOSE
698  output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
699  for( int j = 0; j < nOverlapFaces; j++ )
700  {
701  for( int p = 0; p < nP; p++ )
702  {
703  for( int q = 0; q < nP; q++ )
704  {
705  output_file << dRemapCoeff[p][q][j] << " ";
706  }
707  }
708  }
709  output_file << std::endl;
710 #endif
711 
712  // Force consistency and conservation
713  if( !fNoConservation )
714  {
715  double dTargetArea = 0.0;
716  for( int j = 0; j < nOverlapFaces; j++ )
717  {
718  dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
719  }
720 
721  for( int p = 0; p < nP; p++ )
722  {
723  for( int q = 0; q < nP; q++ )
724  {
725  vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
726  }
727  }
728 
729  const double areaTolerance = 1e-10;
730  // Source elements are completely covered by target volumes
731  if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
732  {
733  vecTargetArea.Allocate( nOverlapFaces );
734  for( int j = 0; j < nOverlapFaces; j++ )
735  {
736  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
737  }
738 
739  dCoeff.Allocate( nOverlapFaces, nP * nP );
740 
741  for( int j = 0; j < nOverlapFaces; j++ )
742  {
743  for( int p = 0; p < nP; p++ )
744  {
745  for( int q = 0; q < nP; q++ )
746  {
747  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
748  }
749  }
750  }
751 
752  // Target volumes only partially cover source elements
753  }
754  else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
755  {
756  double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
757 
758  vecTargetArea.Allocate( nOverlapFaces + 1 );
759  for( int j = 0; j < nOverlapFaces; j++ )
760  {
761  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
762  }
763  vecTargetArea[nOverlapFaces] = dExtraneousArea;
764 
765 #ifdef VERBOSE
766  Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
767  m_meshInputCov->vecFaceArea[ixFirst] );
768 #endif
769  if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
770  {
771  _EXCEPTIONT( "Partial element area exceeds total element area" );
772  }
773 
774  dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
775 
776  for( int j = 0; j < nOverlapFaces; j++ )
777  {
778  for( int p = 0; p < nP; p++ )
779  {
780  for( int q = 0; q < nP; q++ )
781  {
782  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
783  }
784  }
785  }
786  for( int p = 0; p < nP; p++ )
787  {
788  for( int q = 0; q < nP; q++ )
789  {
790  dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
791  }
792  }
793  for( int j = 0; j < nOverlapFaces; j++ )
794  {
795  for( int p = 0; p < nP; p++ )
796  {
797  for( int q = 0; q < nP; q++ )
798  {
799  dCoeff[nOverlapFaces][p * nP + q] -=
800  dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
801  }
802  }
803  }
804  for( int p = 0; p < nP; p++ )
805  {
806  for( int q = 0; q < nP; q++ )
807  {
808  dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
809  }
810  }
811 
812  // Source elements only partially cover target volumes
813  }
814  else
815  {
816  Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
817  m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
818  _EXCEPTIONT( "Target grid must be a subset of source grid" );
819  }
820 
821  ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
822  /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
823 
824  for( int j = 0; j < nOverlapFaces; j++ )
825  {
826  for( int p = 0; p < nP; p++ )
827  {
828  for( int q = 0; q < nP; q++ )
829  {
830  dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
831  }
832  }
833  }
834  }
835 
836 #ifdef VERBOSE
837  // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
838  // for ( int j = 0; j < nOverlapFaces; j++ )
839  // {
840  // for ( int p = 0; p < nP; p++ )
841  // {
842  // for ( int q = 0; q < nP; q++ )
843  // {
844  // output_file << dRemapCoeff[p][q][j] << " ";
845  // }
846  // }
847  // }
848  // output_file << std::endl;
849 #endif
850 
851  // Put these remap coefficients into the SparseMatrix map
852  for( int j = 0; j < nOverlapFaces; j++ )
853  {
854  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
855 
856  // signal to not participate, because it is a ghost target
857  if( ixSecondFace < 0 ) continue; // do not do anything
858 
859  for( int p = 0; p < nP; p++ )
860  {
861  for( int q = 0; q < nP; q++ )
862  {
863  if( fContinuousIn )
864  {
865  int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
866 
867  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
868  m_meshOverlap->vecFaceArea[ixOverlap + j] /
869  m_meshOutput->vecFaceArea[ixSecondFace];
870  }
871  else
872  {
873  int ixFirstNode = ixFirst * nP * nP + p * nP + q;
874 
875  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
876  m_meshOverlap->vecFaceArea[ixOverlap + j] /
877  m_meshOutput->vecFaceArea[ixSecondFace];
878  }
879  }
880  }
881  }
882  // Increment the current overlap index
883  ixOverlap += nOverlapFaces;
884  }
885 #ifdef VERBOSE
886  output_file.flush(); // required here
887  output_file.close();
888 #endif
889 
890  return;
891 }

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

◆ ReadParallelMap()

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

Generate the metadata associated with the offline map.

Read the OfflineMap from a NetCDF file.

Definition at line 1194 of file TempestOnlineMapIO.cpp.

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

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

◆ set_col_dc_dofs()

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

Definition at line 647 of file TempestOnlineMap.cpp.

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

References MB_SUCCESS.

◆ set_row_dc_dofs()

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

Definition at line 670 of file TempestOnlineMap.cpp.

671 {
672  // row_dtoc_dofmap = values_entities; // needs to point to local
673  // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
674 
675  row_dtoc_dofmap.resize( values_entities.size() );
676  for( int j = 0; j < (int)values_entities.size(); j++ )
677  {
678  if( rowMap.find( values_entities[j] - 1 ) != rowMap.end() )
679  row_dtoc_dofmap[j] = rowMap[values_entities[j] - 1]; // values are 1 based, but rowMap, colMap are not
680  else
681  {
682  row_dtoc_dofmap[j] = -1; // not all values are used
683  // std::cout <<"values_entities[j] - 1: " << values_entities[j] - 1 <<" at index j = " << j << " not
684  // found in rowMap \n";
685  }
686  }
687  return moab::MB_SUCCESS;
688 }

References MB_SUCCESS.

◆ SetDestinationNDofsPerElement()

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

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

Definition at line 510 of file TempestOnlineMap.hpp.

511 {
512  m_nDofsPEl_Dest = nt;
513 }

◆ SetDOFmapAssociation()

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

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

Definition at line 172 of file TempestOnlineMap.cpp.

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

References ErrorCode, MB_CHK_ERR, and MB_SUCCESS.

◆ SetDOFmapTags()

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

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

Definition at line 116 of file TempestOnlineMap.cpp.

118 {
119  moab::ErrorCode rval;
120 
121  int tagSize = 0;
123  rval =
124  m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
125 
127  {
128  int ntot_elements = 0, nelements = m_remapper->m_source_entities.size();
129 #ifdef MOAB_HAVE_MPI
130  int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
131  if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
132 #else
133  ntot_elements = nelements;
134 #endif
135 
137  &m_remapper->m_source_entities, srcDofTagName, m_nDofsPEl_Src );MB_CHK_ERR( rval );
138 
139  rval = m_interface->tag_get_handle( srcDofTagName.c_str(), m_nDofsPEl_Src * m_nDofsPEl_Src, MB_TYPE_INTEGER,
140  this->m_dofTagSrc, MB_TAG_ANY );MB_CHK_ERR( rval );
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  int ntot_elements = 0, nelements = m_remapper->m_target_entities.size();
151 #ifdef MOAB_HAVE_MPI
152  int ierr = MPI_Allreduce( &nelements, &ntot_elements, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
153  if( ierr != 0 ) MB_CHK_SET_ERR( MB_FAILURE, "MPI_Allreduce failed to get total source elements" );
154 #else
155  ntot_elements = nelements;
156 #endif
157 
158  rval = m_remapper->GenerateCSMeshMetadata( ntot_elements, m_remapper->m_target_entities, NULL, tgtDofTagName,
159  m_nDofsPEl_Dest );MB_CHK_ERR( rval );
160 
161  rval = m_interface->tag_get_handle( tgtDofTagName.c_str(), m_nDofsPEl_Dest * m_nDofsPEl_Dest, MB_TYPE_INTEGER,
162  this->m_dofTagDest, MB_TAG_ANY );MB_CHK_ERR( rval );
163  }
164  else
165  MB_CHK_ERR( rval );
166 
167  return moab::MB_SUCCESS;
168 }

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

◆ SetSourceNDofsPerElement()

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

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

Definition at line 505 of file TempestOnlineMap.hpp.

506 {
507  m_nDofsPEl_Src = ns;
508 }

◆ setup_sizes_dimensions()

void moab::TempestOnlineMap::setup_sizes_dimensions ( )
private

Definition at line 78 of file TempestOnlineMap.cpp.

79 {
80  if( m_meshInputCov )
81  {
82  std::vector< std::string > dimNames;
83  std::vector< int > dimSizes;
84  dimNames.push_back( "num_elem" );
85  dimSizes.push_back( m_meshInputCov->faces.size() );
86 
87  this->InitializeSourceDimensions( dimNames, dimSizes );
88  }
89 
90  if( m_meshOutput )
91  {
92  std::vector< std::string > dimNames;
93  std::vector< int > dimSizes;
94  dimNames.push_back( "num_elem" );
95  dimSizes.push_back( m_meshOutput->faces.size() );
96 
97  this->InitializeTargetDimensions( dimNames, dimSizes );
98  }
99 }

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

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

◆ col_gdofmap

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

Definition at line 442 of file TempestOnlineMap.hpp.

Referenced by fill_col_ids(), and LinearRemapNN_MOAB().

◆ colMap

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

Definition at line 447 of file TempestOnlineMap.hpp.

◆ dataGLLNodesDest

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

Definition at line 449 of file TempestOnlineMap.hpp.

◆ dataGLLNodesSrc

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

Definition at line 449 of file TempestOnlineMap.hpp.

◆ dataGLLNodesSrcCov

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

Definition at line 449 of file TempestOnlineMap.hpp.

◆ is_parallel

bool moab::TempestOnlineMap::is_parallel
private

Definition at line 464 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ is_root

bool moab::TempestOnlineMap::is_root
private

Definition at line 464 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_bConserved

bool moab::TempestOnlineMap::m_bConserved
private

Definition at line 456 of file TempestOnlineMap.hpp.

◆ m_destDiscType

DiscretizationType moab::TempestOnlineMap::m_destDiscType
private

Definition at line 450 of file TempestOnlineMap.hpp.

◆ m_dofTagDest

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

Definition at line 441 of file TempestOnlineMap.hpp.

◆ m_dofTagSrc

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

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

Definition at line 441 of file TempestOnlineMap.hpp.

◆ m_eInputType

DiscretizationType moab::TempestOnlineMap::m_eInputType
private

Definition at line 455 of file TempestOnlineMap.hpp.

◆ m_eOutputType

DiscretizationType moab::TempestOnlineMap::m_eOutputType
private

Definition at line 455 of file TempestOnlineMap.hpp.

◆ m_iMonotonicity

int moab::TempestOnlineMap::m_iMonotonicity
private

Definition at line 457 of file TempestOnlineMap.hpp.

◆ m_interface

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

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

Definition at line 429 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshInput

Mesh* moab::TempestOnlineMap::m_meshInput
private

Definition at line 459 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshInputCov

Mesh* moab::TempestOnlineMap::m_meshInputCov
private

Definition at line 460 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshOutput

Mesh* moab::TempestOnlineMap::m_meshOutput
private

Definition at line 461 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_meshOverlap

Mesh* moab::TempestOnlineMap::m_meshOverlap
private

Definition at line 462 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_nDofsPEl_Dest

int moab::TempestOnlineMap::m_nDofsPEl_Dest
private

Definition at line 454 of file TempestOnlineMap.hpp.

◆ m_nDofsPEl_Src

int moab::TempestOnlineMap::m_nDofsPEl_Src
private

Definition at line 454 of file TempestOnlineMap.hpp.

◆ m_nTotDofs_Dest

int moab::TempestOnlineMap::m_nTotDofs_Dest
private

Definition at line 451 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

◆ m_nTotDofs_Src

int moab::TempestOnlineMap::m_nTotDofs_Src
private

Definition at line 451 of file TempestOnlineMap.hpp.

◆ m_nTotDofs_SrcCov

int moab::TempestOnlineMap::m_nTotDofs_SrcCov
private

Definition at line 451 of file TempestOnlineMap.hpp.

Referenced by LinearRemapNN_MOAB().

◆ m_remapper

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

The fundamental remapping operator object.

Definition at line 416 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ m_srcDiscType

DiscretizationType moab::TempestOnlineMap::m_srcDiscType
private

Definition at line 450 of file TempestOnlineMap.hpp.

◆ rank

int moab::TempestOnlineMap::rank
private

Definition at line 465 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ row_dtoc_dofmap

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

Definition at line 445 of file TempestOnlineMap.hpp.

◆ row_gdofmap

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

Definition at line 442 of file TempestOnlineMap.hpp.

Referenced by fill_row_ids(), and LinearRemapNN_MOAB().

◆ rowMap

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

Definition at line 447 of file TempestOnlineMap.hpp.

◆ size

int moab::TempestOnlineMap::size
private

Definition at line 465 of file TempestOnlineMap.hpp.

Referenced by TempestOnlineMap().

◆ srccol_dtoc_dofmap

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

Definition at line 445 of file TempestOnlineMap.hpp.

◆ srccol_gdofmap

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

Definition at line 442 of file TempestOnlineMap.hpp.


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