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