Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
TempestOnlineMap.cpp
Go to the documentation of this file.
1 /*
2  * =====================================================================================
3  *
4  * Filename: TempestOnlineMap.hpp
5  *
6  * Description: Interface to the TempestRemap library to compute the consistent,
7  * and accurate high-order conservative remapping weights for overlap
8  * grids on the sphere in climate simulations.
9  *
10  * Author: Vijay S. Mahadevan (vijaysm), [email protected]
11  *
12  * =====================================================================================
13  */
14 
15 #include "Announce.h"
16 #include "DataArray3D.h"
17 #include "FiniteVolumeTools.h"
18 #include "FiniteElementTools.h"
19 #include "TriangularQuadrature.h"
20 #include "GaussQuadrature.h"
21 #include "GaussLobattoQuadrature.h"
22 #include "SparseMatrix.h"
23 #include "STLStringHelper.h"
24 #include "LinearRemapFV.h"
25 
26 #include "LinearRemapSE0.h"
27 #include "LinearRemapFV.h"
28 
29 #include "moab/Remapping/TempestOnlineMap.hpp"
30 #include "DebugOutput.hpp"
31 #include "moab/TupleList.hpp"
32 #include "moab/MeshTopoUtil.hpp"
33 
34 #include <fstream>
35 #include <cmath>
36 #include <cstdlib>
37 #include <numeric>
38 #include <algorithm>
39 
40 #ifdef MOAB_HAVE_NETCDFPAR
41 #include "netcdfcpp_par.hpp"
42 #else
43 #include "netcdfcpp.h"
44 #endif
45 
46 // #define USE_NATIVE_TEMPESTREMAP_ROUTINES
47 
48 ///////////////////////////////////////////////////////////////////////////////
49 
50 // #define VERBOSE
51 // #define VVERBOSE
52 // #define CHECK_INCREASING_DOF
53 
54 ///////////////////////////////////////////////////////////////////////////////
55 
56 #define MPI_CHK_ERR( err ) \
57  if( err ) \
58  { \
59  std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
60  std::cout << "\nMPI Aborting... \n"; \
61  return moab::MB_FAILURE; \
62  }
63 
64 moab::TempestOnlineMap::TempestOnlineMap( moab::TempestRemapper* remapper ) : OfflineMap(), m_remapper( remapper )
65 {
66  // Get the references for the MOAB core objects
67  m_interface = m_remapper->get_interface();
68 #ifdef MOAB_HAVE_MPI
69  m_pcomm = m_remapper->get_parallel_communicator();
70 #endif
71 
72  // now let us re-update the reference to the input source mesh
73  m_meshInput = m_remapper->GetMesh( moab::Remapper::SourceMesh );
74  // now let us re-update the reference to the covering mesh
75  m_meshInputCov = m_remapper->GetCoveringMesh();
76  // now let us re-update the reference to the output target mesh
77  m_meshOutput = m_remapper->GetMesh( moab::Remapper::TargetMesh );
78  // now let us re-update the reference to the output target mesh
79  m_meshOverlap = m_remapper->GetMesh( moab::Remapper::OverlapMesh );
80 
81  is_parallel = remapper->is_parallel;
82  is_root = remapper->is_root;
83  rank = remapper->rank;
84  size = remapper->size;
85 
86  // set default order
87  m_input_order = m_output_order = 1;
88 
89  // Initialize dimension information from file
90  this->setup_sizes_dimensions();
91 }
92 
93 void moab::TempestOnlineMap::setup_sizes_dimensions()
94 {
95  if( m_meshInputCov )
96  {
97  std::vector< std::string > dimNames;
98  std::vector< int > dimSizes;
99  dimNames.push_back( "num_elem" );
100  dimSizes.push_back( m_meshInputCov->faces.size() );
101 
102  this->InitializeSourceDimensions( dimNames, dimSizes );
103  }
104 
105  if( m_meshOutput )
106  {
107  std::vector< std::string > dimNames;
108  std::vector< int > dimSizes;
109  dimNames.push_back( "num_elem" );
110  dimSizes.push_back( m_meshOutput->faces.size() );
111 
112  this->InitializeTargetDimensions( dimNames, dimSizes );
113  }
114 }
115 
116 ///////////////////////////////////////////////////////////////////////////////
117 
118 moab::TempestOnlineMap::~TempestOnlineMap()
119 {
120  m_interface = nullptr;
121 #ifdef MOAB_HAVE_MPI
122  m_pcomm = nullptr;
123 #endif
124  m_meshInput = nullptr;
125  m_meshOutput = nullptr;
126  m_meshOverlap = nullptr;
127 }
128 
129 ///////////////////////////////////////////////////////////////////////////////
130 
131 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapTags( const std::string srcDofTagName,
132  const std::string tgtDofTagName )
133 {
134  moab::ErrorCode rval;
135 
136  int tagSize = 0;
137  tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
138  rval =
139  m_interface->tag_get_handle( srcDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagSrc, MB_TAG_ANY );
140 
141  if( rval == moab::MB_TAG_NOT_FOUND && m_eInputType != DiscretizationType_FV )
142  {
143  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for source mesh." );
144  }
145  else
146  MB_CHK_ERR( rval );
147 
148  tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
149  rval =
150  m_interface->tag_get_handle( tgtDofTagName.c_str(), tagSize, MB_TYPE_INTEGER, this->m_dofTagDest, MB_TAG_ANY );
151  if( rval == moab::MB_TAG_NOT_FOUND && m_eOutputType != DiscretizationType_FV )
152  {
153  MB_CHK_SET_ERR( MB_FAILURE, "DoF tag is not set correctly for target mesh." );
154  }
155  else
156  MB_CHK_ERR( rval );
157 
158  return moab::MB_SUCCESS;
159 }
160 
161 ///////////////////////////////////////////////////////////////////////////////
162 
163 moab::ErrorCode moab::TempestOnlineMap::SetDOFmapAssociation( DiscretizationType srcType,
164  int srcOrder,
165  bool isSrcContinuous,
166  DataArray3D< int >* srcdataGLLNodes,
167  DataArray3D< int >* srcdataGLLNodesSrc,
168  DiscretizationType destType,
169  int destOrder,
170  bool isTgtContinuous,
171  DataArray3D< int >* tgtdataGLLNodes )
172 {
173  std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
174  std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
175 
176  // We are assuming that these are element based tags that are sized: np * np
177  m_srcDiscType = srcType;
178  m_destDiscType = destType;
179  m_input_order = srcOrder;
180  m_output_order = destOrder;
181 
182  bool vprint = is_root && false;
183 
184  // Compute and store the total number of source and target DoFs corresponding
185  // to number of rows and columns in the mapping.
186  // Now compute the mapping and store it for the covering mesh
187  int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
188  if( m_remapper->point_cloud_source )
189  {
190  assert( m_nDofsPEl_Src == 1 );
191  col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
192  col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
193  src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
194  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] ) );
195  srcTagSize = 1;
196  }
197  else
198  {
199  col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
200  col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
201  src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
202  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] ) );
203  }
204 
205  m_nTotDofs_SrcCov = 0;
206  if( srcdataGLLNodes == nullptr )
207  {
208  /* we only have a mapping for elements as DoFs */
209  for( unsigned i = 0; i < col_gdofmap.size(); ++i )
210  {
211  auto gdof = src_soln_gdofs[i];
212  assert( gdof > 0 );
213  col_gdofmap[i] = gdof - 1;
214  col_dtoc_dofmap[i] = i;
215  if( vprint ) std::cout << "Col: " << i << ", " << col_gdofmap[i] << "\n";
216  m_nTotDofs_SrcCov++;
217  }
218  }
219  else
220  {
221  if( isSrcContinuous )
222  dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, false );
223  // Put these remap coefficients into the SparseMatrix map
224  for( unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
225  {
226  for( int p = 0; p < m_nDofsPEl_Src; p++ )
227  {
228  for( int q = 0; q < m_nDofsPEl_Src; q++ )
229  {
230  const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
231  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
232  if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
233  {
234  m_nTotDofs_SrcCov++;
235  dgll_cgll_covcol_ldofmap[localDOF] = true;
236  }
237  if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
238  assert( src_soln_gdofs[offsetDOF] > 0 );
239  col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
240  col_dtoc_dofmap[offsetDOF] = localDOF;
241  if( vprint )
242  std::cout << "Col: " << offsetDOF << ", " << localDOF << ", " << col_gdofmap[offsetDOF] << ", "
243  << m_nTotDofs_SrcCov << "\n";
244  }
245  }
246  }
247  }
248 
249  if( m_remapper->point_cloud_source )
250  {
251  assert( m_nDofsPEl_Src == 1 );
252  srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
253  srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
254  locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
255  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] ) );
256  }
257  else
258  {
259  srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
260  srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
261  locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
262  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] ) );
263  }
264 
265  // Now compute the mapping and store it for the original source mesh
266  m_nTotDofs_Src = 0;
267  if( srcdataGLLNodesSrc == nullptr )
268  {
269  /* we only have a mapping for elements as DoFs */
270  for( unsigned i = 0; i < srccol_gdofmap.size(); ++i )
271  {
272  auto gdof = locsrc_soln_gdofs[i];
273  assert( gdof > 0 );
274  srccol_gdofmap[i] = gdof - 1;
275  srccol_dtoc_dofmap[i] = i;
276  m_nTotDofs_Src++;
277  }
278  }
279  else
280  {
281  if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, false );
282  // Put these remap coefficients into the SparseMatrix map
283  for( unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
284  {
285  for( int p = 0; p < m_nDofsPEl_Src; p++ )
286  {
287  for( int q = 0; q < m_nDofsPEl_Src; q++ )
288  {
289  const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
290  const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
291  if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
292  {
293  m_nTotDofs_Src++;
294  dgll_cgll_col_ldofmap[localDOF] = true;
295  }
296  if( !isSrcContinuous ) m_nTotDofs_Src++;
297  assert( locsrc_soln_gdofs[offsetDOF] > 0 );
298  srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
299  srccol_dtoc_dofmap[offsetDOF] = localDOF;
300  }
301  }
302  }
303  }
304 
305  int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
306  if( m_remapper->point_cloud_target )
307  {
308  assert( m_nDofsPEl_Dest == 1 );
309  row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
310  row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
311  tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
312  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] ) );
313  tgtTagSize = 1;
314  }
315  else
316  {
317  row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
318  row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
319  tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
320  MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] ) );
321  }
322 
323  // Now compute the mapping and store it for the target mesh
324  // To access the GID for each row: row_gdofmap [ row_ldofmap [ 0 : local_ndofs ] ] = GDOF
325  m_nTotDofs_Dest = 0;
326  if( tgtdataGLLNodes == nullptr )
327  {
328  /* we only have a mapping for elements as DoFs */
329  for( unsigned i = 0; i < row_gdofmap.size(); ++i )
330  {
331  auto gdof = tgt_soln_gdofs[i];
332  assert( gdof > 0 );
333  row_gdofmap[i] = gdof - 1;
334  row_dtoc_dofmap[i] = i;
335  if( vprint ) std::cout << "Row: " << i << ", " << row_gdofmap[i] << "\n";
336  m_nTotDofs_Dest++;
337  }
338  }
339  else
340  {
341  if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, false );
342  // Put these remap coefficients into the SparseMatrix map
343  for( unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
344  {
345  for( int p = 0; p < m_nDofsPEl_Dest; p++ )
346  {
347  for( int q = 0; q < m_nDofsPEl_Dest; q++ )
348  {
349  const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
350  const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
351  if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
352  {
353  m_nTotDofs_Dest++;
354  dgll_cgll_row_ldofmap[localDOF] = true;
355  }
356  if( !isTgtContinuous ) m_nTotDofs_Dest++;
357  assert( tgt_soln_gdofs[offsetDOF] > 0 );
358  row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
359  row_dtoc_dofmap[offsetDOF] = localDOF;
360  if( vprint )
361  std::cout << "Row: " << offsetDOF << ", " << localDOF << ", " << row_gdofmap[offsetDOF] << ", "
362  << m_nTotDofs_Dest << "\n";
363  }
364  }
365  }
366  }
367 
368  // Let us also allocate the local representation of the sparse matrix
369 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
370  if( vprint )
371  {
372  std::cout << "[" << rank << "]" << "DoFs: row = " << m_nTotDofs_Dest << ", " << row_gdofmap.size()
373  << ", col = " << m_nTotDofs_Src << ", " << m_nTotDofs_SrcCov << ", " << col_gdofmap.size() << "\n";
374  // std::cout << "Max col_dofmap: " << maxcol << ", Min col_dofmap" << mincol << "\n";
375  }
376 #endif
377 
378  // check monotonicity of row_gdofmap and col_gdofmap
379 #ifdef CHECK_INCREASING_DOF
380  for( size_t i = 0; i < row_gdofmap.size() - 1; i++ )
381  {
382  if( row_gdofmap[i] > row_gdofmap[i + 1] )
383  std::cout << " on rank " << rank << " in row_gdofmap[" << i << "]=" << row_gdofmap[i] << " > row_gdofmap["
384  << i + 1 << "]=" << row_gdofmap[i + 1] << " \n";
385  }
386  for( size_t i = 0; i < col_gdofmap.size() - 1; i++ )
387  {
388  if( col_gdofmap[i] > col_gdofmap[i + 1] )
389  std::cout << " on rank " << rank << " in col_gdofmap[" << i << "]=" << col_gdofmap[i] << " > col_gdofmap["
390  << i + 1 << "]=" << col_gdofmap[i + 1] << " \n";
391  }
392 #endif
393 
394  return moab::MB_SUCCESS;
395 }
396 
397 moab::ErrorCode moab::TempestOnlineMap::set_col_dc_dofs( std::vector< int >& values_entities )
398 {
399  // col_gdofmap has global dofs , that should be in the list of values, such that
400  // row_dtoc_dofmap[offsetDOF] = localDOF;
401  // we need to find col_dtoc_dofmap such that: col_gdofmap[ col_dtoc_dofmap[i] ] == values_entities [i];
402  // we know that col_gdofmap[0..(nbcols-1)] = global_col_dofs -> in values_entities
403  // form first inverse
404  //
405  // resize and initialize to -1 to signal that this value should not be used, if not set below
406  col_dtoc_dofmap.resize( values_entities.size(), -1 );
407  for( size_t j = 0; j < values_entities.size(); j++ )
408  {
409  // values are 1 based, but rowMap, colMap are not
410  const auto it = colMap.find( values_entities[j] - 1 );
411  if( it != colMap.end() ) col_dtoc_dofmap[j] = it->second;
412  }
413  return moab::MB_SUCCESS;
414 }
415 
416 moab::ErrorCode moab::TempestOnlineMap::set_row_dc_dofs( std::vector< int >& values_entities )
417 {
418  // we need to find row_dtoc_dofmap such that: row_gdofmap[ row_dtoc_dofmap[i] ] == values_entities [i];
419  // resize and initialize to -1 to signal that this value should not be used, if not set below
420  row_dtoc_dofmap.resize( values_entities.size(), -1 );
421  for( size_t j = 0; j < values_entities.size(); j++ )
422  {
423  // values are 1 based, but rowMap, colMap are not
424  const auto it = rowMap.find( values_entities[j] - 1 );
425  if( it != rowMap.end() ) row_dtoc_dofmap[j] = it->second;
426  }
427  return moab::MB_SUCCESS;
428 }
429 ///////////////////////////////////////////////////////////////////////////////
430 
431 moab::ErrorCode moab::TempestOnlineMap::GenerateRemappingWeights( std::string strInputType,
432  std::string strOutputType,
433  const GenerateOfflineMapAlgorithmOptions& mapOptions,
434  const std::string& srcDofTagName,
435  const std::string& tgtDofTagName )
436 {
437  NcError error( NcError::silent_nonfatal );
438 
439  moab::DebugOutput dbgprint( std::cout, rank, 0 );
440  dbgprint.set_prefix( "[TempestOnlineMap]: " );
441  moab::ErrorCode rval;
442 
443  const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
444  const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
445  const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
446 
447  // Build a matrix of source and target discretization so that we know how
448  // to assign the global DoFs in parallel for the mapping weights.
449  // For example,
450  // for FV->FV: the rows represented target DoFs and cols represent source DoFs
451  try
452  {
453  // Check command line parameters (data type arguments)
454  STLStringHelper::ToLower( strInputType );
455  STLStringHelper::ToLower( strOutputType );
456 
457  DiscretizationType eInputType;
458  DiscretizationType eOutputType;
459 
460  if( strInputType == "fv" )
461  {
462  eInputType = DiscretizationType_FV;
463  }
464  else if( strInputType == "cgll" )
465  {
466  eInputType = DiscretizationType_CGLL;
467  }
468  else if( strInputType == "dgll" )
469  {
470  eInputType = DiscretizationType_DGLL;
471  }
472  else if( strInputType == "pcloud" )
473  {
474  eInputType = DiscretizationType_PCLOUD;
475  }
476  else
477  {
478  _EXCEPTION1( "Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
479  }
480 
481  if( strOutputType == "fv" )
482  {
483  eOutputType = DiscretizationType_FV;
484  }
485  else if( strOutputType == "cgll" )
486  {
487  eOutputType = DiscretizationType_CGLL;
488  }
489  else if( strOutputType == "dgll" )
490  {
491  eOutputType = DiscretizationType_DGLL;
492  }
493  else if( strOutputType == "pcloud" )
494  {
495  eOutputType = DiscretizationType_PCLOUD;
496  }
497  else
498  {
499  _EXCEPTION1( "Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
500  }
501 
502  // set all required input params
503  m_bConserved = !mapOptions.fNoConservation;
504  m_eInputType = eInputType;
505  m_eOutputType = eOutputType;
506 
507  // Method flags
508  std::string strMapAlgorithm( "" );
509  int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
510 
511  // Make an index of method arguments
512  std::set< std::string > setMethodStrings;
513  {
514  int iLast = 0;
515  for( size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
516  {
517  if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] == ';' ) )
518  {
519  std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
520  STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
521  if( strMethodString.length() > 0 )
522  {
523  setMethodStrings.insert( strMethodString );
524  }
525  iLast = i + 1;
526  }
527  }
528  }
529 
530  for( auto it : setMethodStrings )
531  {
532  // Piecewise constant monotonicity
533  if( it == "mono2" )
534  {
535  if( nMonotoneType != 0 )
536  {
537  _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
538  }
539  if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
540  {
541  _EXCEPTIONT( "--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
542  }
543  nMonotoneType = 2;
544 
545  // Piecewise linear monotonicity
546  }
547  else if( it == "mono3" )
548  {
549  if( nMonotoneType != 0 )
550  {
551  _EXCEPTIONT( "Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
552  }
553  if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
554  {
555  _EXCEPTIONT( "--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
556  }
557  nMonotoneType = 3;
558 
559  // Volumetric remapping from FV to GLL
560  }
561  else if( it == "volumetric" )
562  {
563  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
564  {
565  _EXCEPTIONT( "--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
566  }
567  strMapAlgorithm = "volumetric";
568 
569  // Inverse distance mapping
570  }
571  else if( it == "invdist" )
572  {
573  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
574  {
575  _EXCEPTIONT( "--method \"invdist\" may only be used for FV->FV remapping" );
576  }
577  strMapAlgorithm = "invdist";
578 
579  // Delaunay triangulation mapping
580  }
581  else if( it == "delaunay" )
582  {
583  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
584  {
585  _EXCEPTIONT( "--method \"delaunay\" may only be used for FV->FV remapping" );
586  }
587  strMapAlgorithm = "delaunay";
588 
589  // Bilinear
590  }
591  else if( it == "bilin" )
592  {
593  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
594  {
595  _EXCEPTIONT( "--method \"bilin\" may only be used for FV->FV remapping" );
596  }
597  strMapAlgorithm = "fvbilin";
598 
599  // Integrated bilinear (same as mono3 when source grid is CGLL/DGLL)
600  }
601  else if( it == "intbilin" )
602  {
603  if( m_eOutputType != DiscretizationType_FV )
604  {
605  _EXCEPTIONT( "--method \"intbilin\" may only be used when mapping to FV." );
606  }
607  if( m_eInputType == DiscretizationType_FV )
608  {
609  strMapAlgorithm = "fvintbilin";
610  }
611  else
612  {
613  strMapAlgorithm = "mono3";
614  }
615 
616  // Integrated bilinear with generalized Barycentric coordinates
617  }
618  else if( it == "intbilingb" )
619  {
620  if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
621  {
622  _EXCEPTIONT( "--method \"intbilingb\" may only be used for FV->FV remapping" );
623  }
624  strMapAlgorithm = "fvintbilingb";
625  }
626  else
627  {
628  _EXCEPTION1( "Invalid --method argument \"%s\"", it.c_str() );
629  }
630  }
631 
632  m_nDofsPEl_Src =
633  ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
634  : mapOptions.nPin );
635  m_nDofsPEl_Dest =
636  ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
637  : mapOptions.nPout );
638 
639  // Set the source and target mesh objects
640  MB_CHK_ERR( SetDOFmapTags( srcDofTagName, tgtDofTagName ) );
641 
642  /// the tag should be created already in the e3sm workflow; if not, create it here
643  Tag areaTag;
644  rval = m_interface->tag_get_handle( "aream", 1, MB_TYPE_DOUBLE, areaTag,
645  MB_TAG_DENSE | MB_TAG_EXCL | MB_TAG_CREAT );
646  if( MB_ALREADY_ALLOCATED == rval )
647  {
648  if( is_root ) dbgprint.printf( 0, "aream tag already defined \n" );
649  }
650 
651  double local_areas[3] = { 0.0, 0.0, 0.0 }, global_areas[3] = { 0.0, 0.0, 0.0 };
652  if( !m_bPointCloudSource )
653  {
654  // Calculate Input Mesh Face areas
655  if( is_root ) dbgprint.printf( 0, "Calculating input mesh Face areas\n" );
656  local_areas[0] = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
657  // Set source element areas as tag on the source mesh
658  MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea ) );
659 
660  // Update coverage source mesh areas as well.
661  m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
662  }
663 
664  if( !m_bPointCloudTarget )
665  {
666  // Calculate Output Mesh Face areas
667  if( is_root ) dbgprint.printf( 0, "Calculating output mesh Face areas\n" );
668  local_areas[1] = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
669  // Set target element areas as tag on the target mesh
670  MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea ) );
671  }
672 
673  if( !m_bPointCloud )
674  {
675  // Verify that overlap mesh is in the correct order (sanity check)
676  assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->vecTargetFaceIx.size() );
677 
678  // Calculate Face areas
679  if( is_root ) dbgprint.printf( 0, "Calculating overlap mesh Face areas\n" );
680  local_areas[2] =
681  m_meshOverlap->CalculateFaceAreas( mapOptions.fSourceConcave || mapOptions.fTargetConcave );
682 
683  // store it as global output for now - used later in reduction
684  std::copy( local_areas, local_areas + 3, global_areas );
685 #ifdef MOAB_HAVE_MPI
686  // reduce the local source, target and overlap mesh areas to global areas
687  if( m_pcomm && is_parallel )
688  MPI_Reduce( local_areas, global_areas, 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
689 #endif
690  if( is_root )
691  {
692  dbgprint.printf( 0, "Input Mesh Geometric Area: %1.15e\n", global_areas[0] );
693  dbgprint.printf( 0, "Output Mesh Geometric Area: %1.15e\n", global_areas[1] );
694  dbgprint.printf( 0, "Overlap Mesh Recovered Area: %1.15e\n", global_areas[2] );
695  }
696 
697  // Correct areas to match the areas calculated in the overlap mesh
698  constexpr bool fCorrectAreas = true;
699  if( fCorrectAreas ) // In MOAB-TempestRemap, we will always keep this to be true
700  {
701  if( is_root ) dbgprint.printf( 0, "Correcting source/target areas to overlap mesh areas\n" );
702  DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
703  DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
704 
705  assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
706  assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
707  assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
708 
709  assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
710  assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
711 
712  for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
713  {
714  if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
715  continue; // skip this cell since it is ghosted
716 
717  // let us recompute the source/target areas based on overlap mesh areas
718  assert( static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
719  dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
720  assert( static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
721  dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
722  }
723 
724  for( size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
725  {
726  if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
727  {
728  m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
729  }
730  }
731  for( size_t i = 0; i < m_meshOutput->faces.size(); i++ )
732  {
733  if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
734  {
735  m_meshOutput->vecFaceArea[i] = dTargetArea[i];
736  }
737  }
738  }
739 
740  // Set source mesh areas in map
741  if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
742  {
743  this->SetSourceAreas( m_meshInputCov->vecFaceArea );
744  if( m_meshInputCov->vecMask.size() )
745  {
746  this->SetSourceMask( m_meshInputCov->vecMask );
747  }
748  }
749 
750  // Set target mesh areas in map
751  if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
752  {
753  this->SetTargetAreas( m_meshOutput->vecFaceArea );
754  if( m_meshOutput->vecMask.size() )
755  {
756  this->SetTargetMask( m_meshOutput->vecMask );
757  }
758  }
759 
760  /*
761  // Recalculate input mesh area from overlap mesh
762  if (fabs(dTotalAreaOverlap - dTotalAreaInput) > 1.0e-10) {
763  dbgprint.printf(0, "Overlap mesh only covers a sub-area of the sphere\n");
764  dbgprint.printf(0, "Recalculating source mesh areas\n");
765  dTotalAreaInput = m_meshInput->CalculateFaceAreasFromOverlap(m_meshOverlap);
766  dbgprint.printf(0, "New Input Mesh Geometric Area: %1.15e\n", dTotalAreaInput);
767  }
768  */
769  }
770 
771  // Finite volume input / Finite volume output
772  if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
773  {
774  // Generate reverse node array and edge map
775  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
776  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
777 
778  // Initialize coordinates for map
779  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
780  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
781 
782  this->m_pdataGLLNodesIn = nullptr;
783  this->m_pdataGLLNodesOut = nullptr;
784 
785  // Finite volume input / Finite element output
786  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
787  mapOptions.nPout, false, nullptr );MB_CHK_ERR( rval );
788 
789  // Construct remap for FV-FV
790  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
791 
792  // Construct OfflineMap
793  if( strMapAlgorithm == "invdist" )
794  {
795  if( is_root ) dbgprint.printf( 0, "Calculating map (invdist)\n" );
796  if( m_meshInputCov->faces.size() )
797  LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
798  }
799  else if( strMapAlgorithm == "delaunay" )
800  {
801  if( is_root ) dbgprint.printf( 0, "Calculating map (delaunay)\n" );
802  if( m_meshInputCov->faces.size() )
803  LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
804  }
805  else if( strMapAlgorithm == "fvintbilin" )
806  {
807  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilin)\n" );
808  if( m_meshInputCov->faces.size() )
809  LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
810  }
811  else if( strMapAlgorithm == "fvintbilingb" )
812  {
813  if( is_root ) dbgprint.printf( 0, "Calculating map (intbilingb)\n" );
814  if( m_meshInputCov->faces.size() )
815  LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
816  *this );
817  }
818  else if( strMapAlgorithm == "fvbilin" )
819  {
820 #ifdef VERBOSE
821  if( is_root )
822  {
823  m_meshInputCov->Write( "SourceMeshMBTR.g" );
824  m_meshOutput->Write( "TargetMeshMBTR.g" );
825  }
826  else
827  {
828  m_meshInputCov->Write( "SourceMeshMBTR" + std::to_string( rank ) + ".g" );
829  m_meshOutput->Write( "TargetMeshMBTR" + std::to_string( rank ) + ".g" );
830  }
831 #endif
832  if( is_root ) dbgprint.printf( 0, "Calculating map (bilin)\n" );
833  if( m_meshInputCov->faces.size() )
834  LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *this );
835  }
836  else
837  {
838  if( is_root ) dbgprint.printf( 0, "Calculating conservative FV-FV map\n" );
839  if( m_meshInputCov->faces.size() )
840  {
841 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
842  LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
843  ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *this );
844 #else
845  LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
846 #endif
847  }
848  }
849  }
850  else if( eInputType == DiscretizationType_FV )
851  {
852  DataArray3D< double > dataGLLJacobian;
853 
854  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
855  double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
856  dataGLLNodesDest, dataGLLJacobian );
857 
858  double dNumericalArea = dNumericalArea_loc;
859 #ifdef MOAB_HAVE_MPI
860  if( m_pcomm )
861  MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
862 #endif
863  if( is_root ) dbgprint.printf( 0, "Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
864 
865  // Initialize coordinates for map
866  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
867  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
868 
869  this->m_pdataGLLNodesIn = nullptr;
870  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
871 
872  // Generate the continuous Jacobian
873  bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
874 
875  if( eOutputType == DiscretizationType_CGLL )
876  {
877  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
878  }
879  else
880  {
881  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
882  }
883 
884  // Generate reverse node array and edge map
885  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
886  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
887 
888  // Finite volume input / Finite element output
889  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, false, nullptr, nullptr, eOutputType,
890  mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
891  &dataGLLNodesDest );MB_CHK_ERR( rval );
892 
893  // Generate remap weights
894  if( strMapAlgorithm == "volumetric" )
895  {
896  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL (volumetric)\n" );
897  LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
898  dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *this,
899  nMonotoneType, fContinuous, mapOptions.fNoConservation );
900  }
901  else
902  {
903  if( is_root ) dbgprint.printf( 0, "Calculating remapping weights for FV->GLL\n" );
904  LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
905  this->GetTargetAreas(), mapOptions.nPin, *this, nMonotoneType, fContinuous,
906  mapOptions.fNoConservation );
907  }
908  }
909  else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
910  {
911  DataArray3D< double > dataGLLJacobian;
912  if( !m_bPointCloudSource )
913  {
914  // Generate reverse node array and edge map
915  if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
916  if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap( false );
917 
918  // Initialize coordinates for map
919  if( eInputType == DiscretizationType_FV )
920  {
921  this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
922  }
923  else
924  {
925  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
926  DataArray3D< double > dataGLLJacobianSrc;
927  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
928  dataGLLJacobian );
929  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
930  dataGLLJacobianSrc );
931  }
932  }
933  // else { /* Source is a point cloud dataset */ }
934 
935  if( !m_bPointCloudTarget )
936  {
937  // Generate reverse node array and edge map
938  if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
939  if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap( false );
940 
941  // Initialize coordinates for map
942  if( eOutputType == DiscretizationType_FV )
943  {
944  this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
945  }
946  else
947  {
948  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
949  GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
950  dataGLLJacobian );
951  }
952  }
953  // else { /* Target is a point cloud dataset */ }
954 
955  // Finite volume input / Finite element output
956  rval = this->SetDOFmapAssociation(
957  eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
958  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrcCov ),
959  ( m_bPointCloudSource || eInputType == DiscretizationType_FV ? nullptr : &dataGLLNodesSrc ),
960  eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
961  ( m_bPointCloudTarget ? nullptr : &dataGLLNodesDest ) );MB_CHK_ERR( rval );
962 
963  // Construct remap
964  if( is_root ) dbgprint.printf( 0, "Calculating remap weights with Nearest-Neighbor method\n" );
965  rval = LinearRemapNN_MOAB( true /*use_GID_matching*/, false /*strict_check*/ );MB_CHK_ERR( rval );
966  }
967  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
968  {
969  DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
970 
971  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
972  // generate metadata for the input meshes (both source and covering source)
973  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
974  dataGLLJacobianSrc );
975  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
976  dataGLLJacobian );
977 
978  if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
979  {
980  _EXCEPTIONT( "Number of element does not match between metadata and "
981  "input mesh" );
982  }
983 
984  // Initialize coordinates for map
985  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
986  this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
987 
988  // Generate the continuous Jacobian for input mesh
989  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
990 
991  if( eInputType == DiscretizationType_CGLL )
992  {
993  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
994  }
995  else
996  {
997  GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
998  }
999 
1000  // Finite element input / Finite volume output
1001  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1002  &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1003  false, nullptr );MB_CHK_ERR( rval );
1004 
1005  // Generate remap
1006  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1007 
1008  if( strMapAlgorithm == "volumetric" )
1009  {
1010  _EXCEPTIONT( "Unimplemented: Volumetric currently unavailable for"
1011  "GLL input mesh" );
1012  }
1013 
1014  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1015  this->m_pdataGLLNodesOut = nullptr;
1016 
1017 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1018  LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
1019  nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
1020  *this );
1021 #else
1022  LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
1023  mapOptions.fNoConservation );
1024 #endif
1025  }
1026  else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
1027  {
1028  DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
1029  DataArray3D< double > dataGLLJacobianOut;
1030 
1031  // Input metadata
1032  if( is_root ) dbgprint.printf( 0, "Generating input mesh meta data\n" );
1033  // generate metadata for the input meshes (both source and covering source)
1034  GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1035  dataGLLJacobianSrc );
1036  // now coverage
1037  GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1038  dataGLLJacobianIn );
1039  // Output metadata
1040  if( is_root ) dbgprint.printf( 0, "Generating output mesh meta data\n" );
1041  GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1042  dataGLLJacobianOut );
1043 
1044  // Initialize coordinates for map
1045  this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1046  this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1047 
1048  // Generate the continuous Jacobian for input mesh
1049  bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1050 
1051  if( eInputType == DiscretizationType_CGLL )
1052  {
1053  GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
1054  }
1055  else
1056  {
1057  GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
1058  }
1059 
1060  // Generate the continuous Jacobian for output mesh
1061  bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
1062 
1063  if( eOutputType == DiscretizationType_CGLL )
1064  {
1065  GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
1066  }
1067  else
1068  {
1069  GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
1070  }
1071 
1072  // Input Finite Element to Output Finite Element
1073  rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1074  &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1075  ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );MB_CHK_ERR( rval );
1076 
1077  this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1078  this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1079 
1080  // Generate remap
1081  if( is_root ) dbgprint.printf( 0, "Calculating remap weights\n" );
1082 
1083 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1084  LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
1085  dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1086  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1087  fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *this );
1088 #else
1089  LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1090  this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1091  fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
1092 #endif
1093  }
1094  else
1095  {
1096  _EXCEPTIONT( "Not implemented" );
1097  }
1098 
1099 #ifdef MOAB_HAVE_EIGEN3
1100  copy_tempest_sparsemat_to_eigen3();
1101 #endif
1102 
1103 #ifdef MOAB_HAVE_MPI
1104  {
1105  // Remove ghosted entities from overlap set
1106  moab::Range ghostedEnts;
1107  rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );MB_CHK_ERR( rval );
1108  moab::EntityHandle m_meshOverlapSet = m_remapper->GetMeshSet( moab::Remapper::OverlapMesh );
1109  rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );MB_CHK_SET_ERR( rval, "Deleting ghosted entities failed" );
1110  }
1111 #endif
1112  // Verify consistency, conservation and monotonicity, globally
1113  if( !mapOptions.fNoCheck )
1114  {
1115  if( is_root ) dbgprint.printf( 0, "Verifying map" );
1116  this->IsConsistent( 1.0e-8 );
1117  if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
1118 
1119  if( nMonotoneType != 0 )
1120  {
1121  this->IsMonotone( 1.0e-12 );
1122  }
1123  }
1124  }
1125  catch( Exception& e )
1126  {
1127  dbgprint.printf( 0, "%s", e.ToString().c_str() );
1128  return ( moab::MB_FAILURE );
1129  }
1130  catch( ... )
1131  {
1132  return ( moab::MB_FAILURE );
1133  }
1134  return moab::MB_SUCCESS;
1135 }
1136 
1137 ///////////////////////////////////////////////////////////////////////////////
1138 
1139 int moab::TempestOnlineMap::IsConsistent( double dTolerance )
1140 {
1141 #ifndef MOAB_HAVE_MPI
1142 
1143  return OfflineMap::IsConsistent( dTolerance );
1144 
1145 #else
1146 
1147  // Get map entries
1148  DataArray1D< int > dataRows;
1149  DataArray1D< int > dataCols;
1150  DataArray1D< double > dataEntries;
1151 
1152  // Calculate row sums
1153  DataArray1D< double > dRowSums;
1154  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1155  dRowSums.Allocate( m_mapRemap.GetRows() );
1156 
1157  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1158  {
1159  dRowSums[dataRows[i]] += dataEntries[i];
1160  }
1161 
1162  // Verify all row sums are equal to 1
1163  int fConsistent = 0;
1164  for( unsigned i = 0; i < dRowSums.GetRows(); i++ )
1165  {
1166  if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1167  {
1168  fConsistent++;
1169  int rowGID = row_gdofmap[i];
1170  Announce( "TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1171  }
1172  }
1173 
1174  int ierr;
1175  int fConsistentGlobal = 0;
1176  ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1177  if( ierr != MPI_SUCCESS ) return -1;
1178 
1179  return fConsistentGlobal;
1180 #endif
1181 }
1182 
1183 ///////////////////////////////////////////////////////////////////////////////
1184 
1185 int moab::TempestOnlineMap::IsConservative( double dTolerance )
1186 {
1187 #ifndef MOAB_HAVE_MPI
1188 
1189  return OfflineMap::IsConservative( dTolerance );
1190 
1191 #else
1192  // return OfflineMap::IsConservative(dTolerance);
1193 
1194  int ierr;
1195  // Get map entries
1196  DataArray1D< int > dataRows;
1197  DataArray1D< int > dataCols;
1198  DataArray1D< double > dataEntries;
1199  const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1200  const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1201 
1202  // Calculate column sums
1203  std::vector< int > dColumnsUnique;
1204  std::vector< double > dColumnSums;
1205 
1206  int nColumns = m_mapRemap.GetColumns();
1207  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1208  dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1209  dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1210 
1211  for( unsigned i = 0; i < dataEntries.GetRows(); i++ )
1212  {
1213  dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1214 
1215  assert( dataCols[i] < m_nTotDofs_SrcCov );
1216 
1217  // GID for column DoFs: col_gdofmap[ col_ldofmap [ dataCols[i] ] ]
1218  int colGID = this->GetColGlobalDoF( dataCols[i] ); // col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1219  // int colGID = col_gdofmap[ col_ldofmap [ dataCols[i] ] ];
1220  dColumnsUnique[dataCols[i]] = colGID;
1221 
1222  // std::cout << "Column dataCols[i]=" << dataCols[i] << " with GID = " << colGID <<
1223  // std::endl;
1224  }
1225 
1226  int rootProc = 0;
1227  std::vector< int > nElementsInProc;
1228  const int nDATA = 3;
1229  nElementsInProc.resize( size * nDATA );
1230  int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1231  ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1232  if( ierr != MPI_SUCCESS ) return -1;
1233 
1234  int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
1235  std::vector< int > dColumnIndices;
1236  std::vector< double > dColumnSourceAreas;
1237  std::vector< double > dColumnSumsTotal;
1238  std::vector< int > displs, rcount;
1239  if( rank == rootProc )
1240  {
1241  displs.resize( size + 1, 0 );
1242  rcount.resize( size, 0 );
1243  int gsum = 0;
1244  for( int ir = 0; ir < size; ++ir )
1245  {
1246  nTotVals += nElementsInProc[ir * nDATA];
1247  nTotColumns += nElementsInProc[ir * nDATA + 1];
1248  nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1249 
1250  displs[ir] = gsum;
1251  rcount[ir] = nElementsInProc[ir * nDATA + 1];
1252  gsum += rcount[ir];
1253 
1254  // printf( "%d: nTotColumns: %d, Displs: %d, rcount: %d, gsum = %d\n", ir, nTotColumns, displs[ir], rcount[ir], gsum );
1255  }
1256 
1257  printf( "Total nnz: %d, global source elements = %d\n", nTotVals, gsum );
1258 
1259  dColumnIndices.resize( nTotColumns, -1 );
1260  dColumnSumsTotal.resize( nTotColumns, 0.0 );
1261  // dColumnSourceAreas.resize ( nTotColumns, 0.0 );
1262  }
1263 
1264  // Gather all ColumnSums to root process and accumulate
1265  // We expect that the sums of all columns equate to 1.0 within user specified tolerance
1266  // Need to do a gatherv here since different processes have different number of elements
1267  // MPI_Reduce(&dColumnSums[0], &dColumnSumsTotal[0], m_mapRemap.GetColumns(), MPI_DOUBLE,
1268  // MPI_SUM, 0, m_pcomm->comm());
1269  ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
1270  displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1271  if( ierr != MPI_SUCCESS ) return -1;
1272  ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
1273  displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1274  if( ierr != MPI_SUCCESS ) return -1;
1275  // ierr = MPI_Gatherv ( &dSourceAreas[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSourceAreas[0],
1276  // rcount.data(), displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() ); if ( ierr !=
1277  // MPI_SUCCESS ) return -1;
1278 
1279  // Clean out unwanted arrays now
1280  dColumnSums.clear();
1281  dColumnsUnique.clear();
1282 
1283  // Verify all column sums equal the input Jacobian
1284  int fConservative = 0;
1285  if( rank == rootProc )
1286  {
1287  displs[size] = ( nTotColumns );
1288  // std::vector<double> dColumnSumsOnRoot(nTotColumnsUnq, 0.0);
1289  std::map< int, double > dColumnSumsOnRoot;
1290  // std::map<int, double> dColumnSourceAreasOnRoot;
1291  for( int ir = 0; ir < size; ir++ )
1292  {
1293  for( int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1294  {
1295  if( dColumnIndices[ips] < 0 ) continue;
1296  // printf("%d, %d: dColumnIndices[ips]: %d\n", ir, ips, dColumnIndices[ips]);
1297  assert( dColumnIndices[ips] < nTotColumnsUnq );
1298  dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips]; // / dColumnSourceAreas[ips];
1299  // dColumnSourceAreasOnRoot[ dColumnIndices[ips] ] = dColumnSourceAreas[ips];
1300  // dColumnSourceAreas[ dColumnIndices[ips] ]
1301  }
1302  }
1303 
1304  for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1305  {
1306  // if ( fabs ( it->second - dColumnSourceAreasOnRoot[it->first] ) > dTolerance )
1307  if( fabs( it->second - 1.0 ) > dTolerance )
1308  {
1309  fConservative++;
1310  Announce( "TempestOnlineMap is not conservative in column "
1311  // "%i (%1.15e)", it->first, it->second );
1312  "%i (%1.15e)",
1313  it->first, it->second /* / dColumnSourceAreasOnRoot[it->first] */ );
1314  }
1315  }
1316  }
1317 
1318  // TODO: Just do a broadcast from root instead of a reduction
1319  ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1320  if( ierr != MPI_SUCCESS ) return -1;
1321 
1322  return fConservative;
1323 #endif
1324 }
1325 
1326 ///////////////////////////////////////////////////////////////////////////////
1327 
1328 int moab::TempestOnlineMap::IsMonotone( double dTolerance )
1329 {
1330 #ifndef MOAB_HAVE_MPI
1331 
1332  return OfflineMap::IsMonotone( dTolerance );
1333 
1334 #else
1335 
1336  // Get map entries
1337  DataArray1D< int > dataRows;
1338  DataArray1D< int > dataCols;
1339  DataArray1D< double > dataEntries;
1340 
1341  m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1342 
1343  // Verify all entries are in the range [0,1]
1344  int fMonotone = 0;
1345  for( unsigned i = 0; i < dataRows.GetRows(); i++ )
1346  {
1347  if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1348  {
1349  fMonotone++;
1350 
1351  Announce( "TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1352  }
1353  }
1354 
1355  int ierr;
1356  int fMonotoneGlobal = 0;
1357  ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1358  if( ierr != MPI_SUCCESS ) return -1;
1359 
1360  return fMonotoneGlobal;
1361 #endif
1362 }
1363 
1364 ///////////////////////////////////////////////////////////////////////////////
1365 
1366 void moab::TempestOnlineMap::ComputeAdjacencyRelations( std::vector< std::unordered_set< int > >& vecAdjFaces,
1367  int nrings,
1368  const Range& entities,
1369  bool useMOABAdjacencies,
1370  Mesh* trMesh )
1371 {
1372  assert( nrings > 0 );
1373  assert( useMOABAdjacencies || trMesh != nullptr );
1374 
1375  const size_t nrows = vecAdjFaces.size();
1376  moab::MeshTopoUtil mtu( m_interface );
1377  for( size_t index = 0; index < nrows; index++ )
1378  {
1379  vecAdjFaces[index].insert( index ); // add self target face first
1380  {
1381  // Compute the adjacent faces to the target face
1382  if( useMOABAdjacencies )
1383  {
1384  moab::Range ents;
1385  // ents.insert( entities.index( entities[index] ) );
1386  ents.insert( entities[index] );
1387  moab::Range adjEnts;
1388  moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, nrings );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
1389  for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
1390  {
1391  // int adjIndex = m_interface->id_from_handle(*it)-1;
1392  int adjIndex = entities.index( *it );
1393  // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
1394  if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
1395  }
1396  }
1397  else
1398  {
1399  /// Vector storing adjacent Faces.
1400  typedef std::pair< int, int > FaceDistancePair;
1401  typedef std::vector< FaceDistancePair > AdjacentFaceVector;
1402  AdjacentFaceVector adjFaces;
1403  Face& face = trMesh->faces[index];
1404  GetAdjacentFaceVectorByEdge( *trMesh, index, nrings * face.edges.size(), adjFaces );
1405 
1406  // Add the adjacent faces to the target face list
1407  for( auto adjFace : adjFaces )
1408  if( adjFace.first >= 0 )
1409  vecAdjFaces[index].insert( adjFace.first ); // map target face to source face
1410  }
1411  }
1412  }
1413 }
1414 
1415 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( moab::Tag srcSolutionTag,
1416  moab::Tag tgtSolutionTag,
1417  bool transpose,
1418  CAASType caasType,
1419  double default_projection )
1420 {
1421  std::vector< double > solSTagVals;
1422  std::vector< double > solTTagVals;
1423 
1424  moab::Range sents, tents;
1425  if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
1426  {
1427  if( m_remapper->point_cloud_source )
1428  {
1429  moab::Range& covSrcEnts = m_remapper->GetMeshVertices( moab::Remapper::CoveringMesh );
1430  solSTagVals.resize( covSrcEnts.size(), default_projection );
1431  sents = covSrcEnts;
1432  }
1433  else
1434  {
1435  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1436  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1437  default_projection );
1438  sents = covSrcEnts;
1439  }
1440  if( m_remapper->point_cloud_target )
1441  {
1442  moab::Range& tgtEnts = m_remapper->GetMeshVertices( moab::Remapper::TargetMesh );
1443  solTTagVals.resize( tgtEnts.size(), default_projection );
1444  tents = tgtEnts;
1445  }
1446  else
1447  {
1448  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1449  solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1450  this->GetDestinationNDofsPerElement(),
1451  default_projection );
1452  tents = tgtEnts;
1453  }
1454  }
1455  else
1456  {
1457  moab::Range& covSrcEnts = m_remapper->GetMeshEntities( moab::Remapper::CoveringMesh );
1458  moab::Range& tgtEnts = m_remapper->GetMeshEntities( moab::Remapper::TargetMesh );
1459  solSTagVals.resize( covSrcEnts.size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1460  default_projection );
1461  solTTagVals.resize( tgtEnts.size() * this->GetDestinationNDofsPerElement() *
1462  this->GetDestinationNDofsPerElement(),
1463  default_projection );
1464 
1465  sents = covSrcEnts;
1466  tents = tgtEnts;
1467  }
1468 
1469  // The tag data is np*np*n_el_src
1470  MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
1471  "Getting local tag data failed" );
1472 
1473  // Compute the application of weights on the suorce solution data and store it in the
1474  // destination solution vector data Optionally, can also perform the transpose application of
1475  // the weight matrix. Set the 3rd argument to true if this is needed
1476  MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, transpose ),
1477  "Applying remap operator onto source vector data failed" );
1478 
1479  // The tag data is np*np*n_el_dest
1480  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1481  "Setting target tag data failed" );
1482 
1483  if( caasType != CAAS_NONE )
1484  {
1485  std::string tgtSolutionTagName;
1486  MB_CHK_SET_ERR( m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName ), "Getting tag name failed" );
1487 
1488  // Perform CAAS iterations iteratively until convergence
1489  constexpr int nmax_caas_iterations = 10;
1490  double mismatch = 1.0;
1491  int caasIteration = 0;
1492  double initialMismatch = 0.0;
1493  while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
1494  caasIteration++ < nmax_caas_iterations ) // iterate until convergence or a maximum of 5 iterations
1495  {
1496  double dMassDiffPostGlobal;
1497  std::pair< double, double > mDefect =
1498  this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
1499 #ifdef MOAB_HAVE_MPI
1500  double dMassDiffPost = mDefect.second;
1501  MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
1502 #else
1503  dMassDiffPostGlobal = mDefect.second;
1504 #endif
1505  if( caasIteration == 1 ) initialMismatch = mDefect.first;
1506  if( m_remapper->verbose && is_root )
1507  {
1508  printf( "Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
1509  tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
1510  }
1511  mismatch = dMassDiffPostGlobal;
1512 
1513  // The tag data is np*np*n_el_dest
1514  MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1515  "Setting local tag data failed" );
1516  }
1517  }
1518 
1519  return moab::MB_SUCCESS;
1520 }
1521 
1522 moab::ErrorCode moab::TempestOnlineMap::DefineAnalyticalSolution( moab::Tag& solnTag,
1523  const std::string& solnName,
1524  moab::Remapper::IntersectionContext ctx,
1525  sample_function testFunction,
1526  moab::Tag* clonedSolnTag,
1527  std::string cloneSolnName )
1528 {
1529  moab::ErrorCode rval;
1530  const bool outputEnabled = ( is_root );
1531  int discOrder;
1532  DiscretizationType discMethod;
1533  // moab::EntityHandle meshset;
1534  moab::Range entities;
1535  Mesh* trmesh;
1536  switch( ctx )
1537  {
1538  case Remapper::SourceMesh:
1539  // meshset = m_remapper->m_covering_source_set;
1540  trmesh = m_remapper->m_covering_source;
1541  entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
1542  : m_remapper->m_covering_source_entities );
1543  discOrder = m_nDofsPEl_Src;
1544  discMethod = m_eInputType;
1545  break;
1546 
1547  case Remapper::TargetMesh:
1548  // meshset = m_remapper->m_target_set;
1549  trmesh = m_remapper->m_target;
1550  entities =
1551  ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
1552  discOrder = m_nDofsPEl_Dest;
1553  discMethod = m_eOutputType;
1554  break;
1555 
1556  default:
1557  if( outputEnabled )
1558  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
1559  return moab::MB_FAILURE;
1560  }
1561 
1562  // Let us create teh solution tag with appropriate information for name, discretization order
1563  // (DoF space)
1564  rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE, solnTag,
1565  MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1566  if( clonedSolnTag != nullptr )
1567  {
1568  if( cloneSolnName.size() == 0 )
1569  {
1570  cloneSolnName = solnName + std::string( "Cloned" );
1571  }
1572  rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder, MB_TYPE_DOUBLE,
1573  *clonedSolnTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_ERR( rval );
1574  }
1575 
1576  // Triangular quadrature rule
1577  const int TriQuadratureOrder = 10;
1578 
1579  if( outputEnabled ) std::cout << "Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
1580 
1581  TriangularQuadratureRule triquadrule( TriQuadratureOrder );
1582 
1583  const int TriQuadraturePoints = triquadrule.GetPoints();
1584 
1585  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1586  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1587 
1588  // Output data
1589  DataArray1D< double > dVar;
1590  DataArray1D< double > dVarMB; // re-arranged local MOAB vector
1591 
1592  // Nodal geometric area
1593  DataArray1D< double > dNodeArea;
1594 
1595  // Calculate element areas
1596  // trmesh->CalculateFaceAreas(fContainsConcaveFaces);
1597 
1598  if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
1599  {
1600  /* Get the spectral points and sample the functionals accordingly */
1601  const bool fGLL = true;
1602  const bool fGLLIntegrate = false;
1603 
1604  // Generate grid metadata
1605  DataArray3D< int > dataGLLNodes;
1606  DataArray3D< double > dataGLLJacobian;
1607 
1608  GenerateMetaData( *trmesh, discOrder, false, dataGLLNodes, dataGLLJacobian );
1609 
1610  // Number of elements
1611  int nElements = trmesh->faces.size();
1612 
1613  // Verify all elements are quadrilaterals
1614  for( int k = 0; k < nElements; k++ )
1615  {
1616  const Face& face = trmesh->faces[k];
1617 
1618  if( face.edges.size() != 4 )
1619  {
1620  _EXCEPTIONT( "Non-quadrilateral face detected; "
1621  "incompatible with --gll" );
1622  }
1623  }
1624 
1625  // Number of unique nodes
1626  int iMaxNode = 0;
1627  for( int i = 0; i < discOrder; i++ )
1628  {
1629  for( int j = 0; j < discOrder; j++ )
1630  {
1631  for( int k = 0; k < nElements; k++ )
1632  {
1633  if( dataGLLNodes[i][j][k] > iMaxNode )
1634  {
1635  iMaxNode = dataGLLNodes[i][j][k];
1636  }
1637  }
1638  }
1639  }
1640 
1641  // Get Gauss-Lobatto quadrature nodes
1642  DataArray1D< double > dG;
1643  DataArray1D< double > dW;
1644 
1645  GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
1646 
1647  // Get Gauss quadrature nodes
1648  const int nGaussP = 10;
1649 
1650  DataArray1D< double > dGaussG;
1651  DataArray1D< double > dGaussW;
1652 
1653  GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
1654 
1655  // Allocate data
1656  dVar.Allocate( iMaxNode );
1657  dVarMB.Allocate( discOrder * discOrder * nElements );
1658  dNodeArea.Allocate( iMaxNode );
1659 
1660  // Sample data
1661  for( int k = 0; k < nElements; k++ )
1662  {
1663  const Face& face = trmesh->faces[k];
1664 
1665  // Sample data at GLL nodes
1666  if( fGLL )
1667  {
1668  for( int i = 0; i < discOrder; i++ )
1669  {
1670  for( int j = 0; j < discOrder; j++ )
1671  {
1672 
1673  // Apply local map
1674  Node node;
1675  Node dDx1G;
1676  Node dDx2G;
1677 
1678  ApplyLocalMap( face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
1679 
1680  // Sample data at this point
1681  double dNodeLon = atan2( node.y, node.x );
1682  if( dNodeLon < 0.0 )
1683  {
1684  dNodeLon += 2.0 * M_PI;
1685  }
1686  double dNodeLat = asin( node.z );
1687 
1688  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1689 
1690  dVar[dataGLLNodes[j][i][k] - 1] = dSample;
1691  }
1692  }
1693  // High-order Gaussian integration over basis function
1694  }
1695  else
1696  {
1697  DataArray2D< double > dCoeff( discOrder, discOrder );
1698 
1699  for( int p = 0; p < nGaussP; p++ )
1700  {
1701  for( int q = 0; q < nGaussP; q++ )
1702  {
1703 
1704  // Apply local map
1705  Node node;
1706  Node dDx1G;
1707  Node dDx2G;
1708 
1709  ApplyLocalMap( face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
1710 
1711  // Cross product gives local Jacobian
1712  Node nodeCross = CrossProduct( dDx1G, dDx2G );
1713 
1714  double dJacobian =
1715  sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
1716 
1717  // Find components of quadrature point in basis
1718  // of the first Face
1719  SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
1720 
1721  // Sample data at this point
1722  double dNodeLon = atan2( node.y, node.x );
1723  if( dNodeLon < 0.0 )
1724  {
1725  dNodeLon += 2.0 * M_PI;
1726  }
1727  double dNodeLat = asin( node.z );
1728 
1729  double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1730 
1731  // Integrate
1732  for( int i = 0; i < discOrder; i++ )
1733  {
1734  for( int j = 0; j < discOrder; j++ )
1735  {
1736 
1737  double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
1738 
1739  dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
1740 
1741  dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
1742  }
1743  }
1744  }
1745  }
1746  }
1747  }
1748 
1749  // Divide by area
1750  if( fGLLIntegrate )
1751  {
1752  for( size_t i = 0; i < dVar.GetRows(); i++ )
1753  {
1754  dVar[i] /= dNodeArea[i];
1755  }
1756  }
1757 
1758  // Let us rearrange the data based on DoF ID specification
1759  if( ctx == Remapper::SourceMesh )
1760  {
1761  for( unsigned j = 0; j < entities.size(); j++ )
1762  for( int p = 0; p < discOrder; p++ )
1763  for( int q = 0; q < discOrder; q++ )
1764  {
1765  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
1766  dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
1767  }
1768  }
1769  else
1770  {
1771  for( unsigned j = 0; j < entities.size(); j++ )
1772  for( int p = 0; p < discOrder; p++ )
1773  for( int q = 0; q < discOrder; q++ )
1774  {
1775  const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
1776  dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
1777  }
1778  }
1779 
1780  // Set the tag data
1781  rval = m_interface->tag_set_data( solnTag, entities, &dVarMB[0] );MB_CHK_ERR( rval );
1782  }
1783  else
1784  {
1785  // assert( discOrder == 1 );
1786  if( discMethod == DiscretizationType_FV )
1787  {
1788  /* Compute an element-wise integral to store the sampled solution based on Quadrature
1789  * rules */
1790  // Resize the array
1791  dVar.Allocate( trmesh->faces.size() );
1792 
1793  std::vector< Node >& nodes = trmesh->nodes;
1794 
1795  // Loop through all Faces
1796  for( size_t i = 0; i < trmesh->faces.size(); i++ )
1797  {
1798  const Face& face = trmesh->faces[i];
1799 
1800  // Loop through all sub-triangles
1801  for( size_t j = 0; j < face.edges.size() - 2; j++ )
1802  {
1803 
1804  const Node& node0 = nodes[face[0]];
1805  const Node& node1 = nodes[face[j + 1]];
1806  const Node& node2 = nodes[face[j + 2]];
1807 
1808  // Triangle area
1809  Face faceTri( 3 );
1810  faceTri.SetNode( 0, face[0] );
1811  faceTri.SetNode( 1, face[j + 1] );
1812  faceTri.SetNode( 2, face[j + 2] );
1813 
1814  double dTriangleArea = CalculateFaceArea( faceTri, nodes );
1815 
1816  // Calculate the element average
1817  double dTotalSample = 0.0;
1818 
1819  // Loop through all quadrature points
1820  for( int k = 0; k < TriQuadraturePoints; k++ )
1821  {
1822  Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
1823  TriQuadratureG[k][2] * node2.x,
1824  TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
1825  TriQuadratureG[k][2] * node2.y,
1826  TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
1827  TriQuadratureG[k][2] * node2.z );
1828 
1829  double dMagnitude = node.Magnitude();
1830  node.x /= dMagnitude;
1831  node.y /= dMagnitude;
1832  node.z /= dMagnitude;
1833 
1834  double dLon = atan2( node.y, node.x );
1835  if( dLon < 0.0 )
1836  {
1837  dLon += 2.0 * M_PI;
1838  }
1839  double dLat = asin( node.z );
1840 
1841  double dSample = ( *testFunction )( dLon, dLat );
1842 
1843  dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
1844  }
1845 
1846  dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
1847  }
1848  }
1849  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
1850  }
1851  else /* discMethod == DiscretizationType_PCLOUD */
1852  {
1853  /* Get the coordinates of the vertices and sample the functionals accordingly */
1854  std::vector< Node >& nodes = trmesh->nodes;
1855 
1856  // Resize the array
1857  dVar.Allocate( nodes.size() );
1858 
1859  for( size_t j = 0; j < nodes.size(); j++ )
1860  {
1861  Node& node = nodes[j];
1862  double dMagnitude = node.Magnitude();
1863  node.x /= dMagnitude;
1864  node.y /= dMagnitude;
1865  node.z /= dMagnitude;
1866  double dLon = atan2( node.y, node.x );
1867  if( dLon < 0.0 )
1868  {
1869  dLon += 2.0 * M_PI;
1870  }
1871  double dLat = asin( node.z );
1872 
1873  double dSample = ( *testFunction )( dLon, dLat );
1874  dVar[j] = dSample;
1875  }
1876 
1877  rval = m_interface->tag_set_data( solnTag, entities, &dVar[0] );MB_CHK_ERR( rval );
1878  }
1879  }
1880 
1881  return moab::MB_SUCCESS;
1882 }
1883 
1884 moab::ErrorCode moab::TempestOnlineMap::ComputeMetrics( moab::Remapper::IntersectionContext ctx,
1885  moab::Tag& exactTag,
1886  moab::Tag& approxTag,
1887  std::map< std::string, double >& metrics,
1888  bool verbose )
1889 {
1890  moab::ErrorCode rval;
1891  const bool outputEnabled = ( is_root );
1892  int discOrder;
1893  // DiscretizationType discMethod;
1894  // moab::EntityHandle meshset;
1895  moab::Range entities;
1896  // Mesh* trmesh;
1897  switch( ctx )
1898  {
1899  case Remapper::SourceMesh:
1900  // meshset = m_remapper->m_covering_source_set;
1901  // trmesh = m_remapper->m_covering_source;
1902  entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
1903  : m_remapper->m_covering_source_entities );
1904  discOrder = m_nDofsPEl_Src;
1905  // discMethod = m_eInputType;
1906  break;
1907 
1908  case Remapper::TargetMesh:
1909  // meshset = m_remapper->m_target_set;
1910  // trmesh = m_remapper->m_target;
1911  entities =
1912  ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
1913  discOrder = m_nDofsPEl_Dest;
1914  // discMethod = m_eOutputType;
1915  break;
1916 
1917  default:
1918  if( outputEnabled )
1919  std::cout << "Invalid context specified for defining an analytical solution tag" << std::endl;
1920  return moab::MB_FAILURE;
1921  }
1922 
1923  // Let us create teh solution tag with appropriate information for name, discretization order
1924  // (DoF space)
1925  std::string exactTagName, projTagName;
1926  const int ntotsize = entities.size() * discOrder * discOrder;
1927  std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
1928  rval = m_interface->tag_get_name( exactTag, exactTagName );MB_CHK_ERR( rval );
1929  rval = m_interface->tag_get_data( exactTag, entities, &exactSolution[0] );MB_CHK_ERR( rval );
1930  rval = m_interface->tag_get_name( approxTag, projTagName );MB_CHK_ERR( rval );
1931  rval = m_interface->tag_get_data( approxTag, entities, &projSolution[0] );MB_CHK_ERR( rval );
1932 
1933  const auto& ovents = m_remapper->m_overlap_entities;
1934 
1935  std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 ); // L1Err, L2Err, LinfErr
1936  double sumarea = 0.0;
1937  for( size_t i = 0; i < ovents.size(); ++i )
1938  {
1939  const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
1940  if( srcidx < 0 ) continue; // Skip non-overlapping entities
1941  const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
1942  if( tgtidx < 0 ) continue; // skip ghost target faces
1943  const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
1944  const double error = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
1945  errnorms[0] += ovarea * error;
1946  errnorms[1] += ovarea * error * error;
1947  errnorms[3] = ( error > errnorms[3] ? error : errnorms[3] );
1948  sumarea += ovarea;
1949  }
1950  errnorms[2] = sumarea;
1951 #ifdef MOAB_HAVE_MPI
1952  if( m_pcomm )
1953  {
1954  MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1955  MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
1956  }
1957 #else
1958  for( int i = 0; i < 4; ++i )
1959  globerrnorms[i] = errnorms[i];
1960 #endif
1961 
1962  globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
1963  globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
1964 
1965  metrics.clear();
1966  metrics["L1Error"] = globerrnorms[0];
1967  metrics["L2Error"] = globerrnorms[1];
1968  metrics["LinfError"] = globerrnorms[3];
1969 
1970  if( verbose && is_root )
1971  {
1972  std::cout << "Error metrics when comparing " << projTagName << " against " << exactTagName << std::endl;
1973  std::cout << "\t Total Intersection area = " << globerrnorms[2] << std::endl;
1974  std::cout << "\t L_1 error = " << globerrnorms[0] << std::endl;
1975  std::cout << "\t L_2 error = " << globerrnorms[1] << std::endl;
1976  std::cout << "\t L_inf error = " << globerrnorms[3] << std::endl;
1977  }
1978 
1979  return moab::MB_SUCCESS;
1980 }