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