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