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"
26 #include "LinearRemapSE0.h"
27 #include "LinearRemapFV.h"
40 #ifdef MOAB_HAVE_NETCDFPAR
43 #include "netcdfcpp.h"
56 #define MPI_CHK_ERR( err ) \
59 std::cout << "MPI Failure. ErrorCode (" << ( err ) << ") "; \
60 std::cout << "\nMPI Aborting... \n"; \
61 return moab::MB_FAILURE; \
69 m_pcomm =
m_remapper->get_parallel_communicator();
97 std::vector< std::string > dimNames;
98 std::vector< int > dimSizes;
99 dimNames.push_back(
"num_elem" );
100 dimSizes.push_back( m_meshInputCov->faces.size() );
102 this->InitializeSourceDimensions( dimNames, dimSizes );
107 std::vector< std::string > dimNames;
108 std::vector< int > dimSizes;
109 dimNames.push_back(
"num_elem" );
110 dimSizes.push_back( m_meshOutput->faces.size() );
112 this->InitializeTargetDimensions( dimNames, dimSizes );
120 m_interface =
nullptr;
124 m_meshInput =
nullptr;
125 m_meshOutput =
nullptr;
126 m_meshOverlap =
nullptr;
132 const std::string tgtDofTagName )
137 tagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
143 MB_CHK_SET_ERR( MB_FAILURE,
"DoF tag is not set correctly for source mesh." );
148 tagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
153 MB_CHK_SET_ERR( MB_FAILURE,
"DoF tag is not set correctly for target mesh." );
165 bool isSrcContinuous,
166 DataArray3D< int >* srcdataGLLNodes,
167 DataArray3D< int >* srcdataGLLNodesSrc,
170 bool isTgtContinuous,
171 DataArray3D< int >* tgtdataGLLNodes )
173 std::vector< bool > dgll_cgll_row_ldofmap, dgll_cgll_col_ldofmap, dgll_cgll_covcol_ldofmap;
174 std::vector< int > src_soln_gdofs, locsrc_soln_gdofs, tgt_soln_gdofs;
177 m_srcDiscType = srcType;
178 m_destDiscType = destType;
179 m_input_order = srcOrder;
180 m_output_order = destOrder;
182 bool vprint = is_root &&
false;
188 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src, -1 );
189 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] );
MB_CHK_ERR( rval );
190 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src );
191 rval = m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] );
MB_CHK_ERR( rval );
192 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * m_nDofsPEl_Dest * m_nDofsPEl_Dest );
193 rval = m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] );
MB_CHK_ERR( rval );
198 std::ofstream output_file(
"sourcecov-gids-0.txt" );
199 output_file <<
"I, GDOF\n";
200 for(
unsigned i = 0; i < src_soln_gdofs.size(); ++i )
201 output_file << i <<
", " << src_soln_gdofs[i] <<
"\n";
203 output_file <<
"ELEMID, IDOF, LDOF, GDOF, NDOF\n";
204 m_nTotDofs_SrcCov = 0;
205 if( isSrcContinuous )
206 dgll_cgll_covcol_ldofmap.resize(
207 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src,
false );
208 for(
unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
210 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
212 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
214 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
215 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
216 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
219 dgll_cgll_covcol_ldofmap[localDOF] =
true;
221 output_file << m_remapper->lid_to_gid_covsrc[j] <<
", " << offsetDOF <<
", " << localDOF
222 <<
", " << src_soln_gdofs[offsetDOF] <<
", " << m_nTotDofs_SrcCov <<
"\n";
228 dgll_cgll_covcol_ldofmap.clear();
232 std::ofstream output_file(
"source-gids-0.txt" );
233 output_file <<
"I, GDOF\n";
234 for(
unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
235 output_file << i <<
", " << locsrc_soln_gdofs[i] <<
"\n";
237 output_file <<
"ELEMID, IDOF, LDOF, GDOF, NDOF\n";
239 if( isSrcContinuous )
240 dgll_cgll_col_ldofmap.resize(
241 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src,
false );
242 for(
unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
244 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
246 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
248 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
249 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
250 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
253 dgll_cgll_col_ldofmap[localDOF] =
true;
255 output_file << m_remapper->lid_to_gid_src[j] <<
", " << offsetDOF <<
", " << localDOF
256 <<
", " << locsrc_soln_gdofs[offsetDOF] <<
", " << m_nTotDofs_Src <<
"\n";
262 dgll_cgll_col_ldofmap.clear();
266 std::ofstream output_file(
"target-gids-0.txt" );
267 output_file <<
"I, GDOF\n";
268 for(
unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
269 output_file << i <<
", " << tgt_soln_gdofs[i] <<
"\n";
271 output_file <<
"ELEMID, IDOF, GDOF, NDOF\n";
274 for(
unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
276 output_file << m_remapper->lid_to_gid_tgt[i] <<
", " << i <<
", " << tgt_soln_gdofs[i] <<
", "
277 << m_nTotDofs_Dest <<
"\n";
288 std::ofstream output_file(
"sourcecov-gids-1.txt" );
289 output_file <<
"I, GDOF\n";
290 for(
unsigned i = 0; i < src_soln_gdofs.size(); ++i )
291 output_file << i <<
", " << src_soln_gdofs[i] <<
"\n";
293 output_file <<
"ELEMID, IDOF, LDOF, GDOF, NDOF\n";
294 m_nTotDofs_SrcCov = 0;
295 if( isSrcContinuous )
296 dgll_cgll_covcol_ldofmap.resize(
297 m_remapper->m_covering_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src,
false );
298 for(
unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
300 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
302 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
304 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
305 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
306 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
309 dgll_cgll_covcol_ldofmap[localDOF] =
true;
311 output_file << m_remapper->lid_to_gid_covsrc[j] <<
", " << offsetDOF <<
", " << localDOF
312 <<
", " << src_soln_gdofs[offsetDOF] <<
", " << m_nTotDofs_SrcCov <<
"\n";
318 dgll_cgll_covcol_ldofmap.clear();
322 std::ofstream output_file(
"source-gids-1.txt" );
323 output_file <<
"I, GDOF\n";
324 for(
unsigned i = 0; i < locsrc_soln_gdofs.size(); ++i )
325 output_file << i <<
", " << locsrc_soln_gdofs[i] <<
"\n";
327 output_file <<
"ELEMID, IDOF, LDOF, GDOF, NDOF\n";
329 if( isSrcContinuous )
330 dgll_cgll_col_ldofmap.resize(
331 m_remapper->m_source_entities.size() * m_nDofsPEl_Src * m_nDofsPEl_Src,
false );
332 for(
unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
334 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
336 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
338 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
339 const int offsetDOF = j * m_nDofsPEl_Src * m_nDofsPEl_Src + p * m_nDofsPEl_Src + q;
340 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
343 dgll_cgll_col_ldofmap[localDOF] =
true;
345 output_file << m_remapper->lid_to_gid_src[j] <<
", " << offsetDOF <<
", " << localDOF
346 <<
", " << locsrc_soln_gdofs[offsetDOF] <<
", " << m_nTotDofs_Src <<
"\n";
352 dgll_cgll_col_ldofmap.clear();
356 std::ofstream output_file(
"target-gids-1.txt" );
357 output_file <<
"I, GDOF\n";
358 for(
unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
359 output_file << i <<
", " << tgt_soln_gdofs[i] <<
"\n";
361 output_file <<
"ELEMID, IDOF, GDOF, NDOF\n";
364 for(
unsigned i = 0; i < tgt_soln_gdofs.size(); ++i )
366 output_file << m_remapper->lid_to_gid_tgt[i] <<
", " << i <<
", " << tgt_soln_gdofs[i] <<
", "
367 << m_nTotDofs_Dest <<
"\n";
379 int srcTagSize = ( m_eInputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Src * m_nDofsPEl_Src );
380 if( m_remapper->point_cloud_source )
382 assert( m_nDofsPEl_Src == 1 );
383 col_gdofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
384 col_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
385 src_soln_gdofs.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
386 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_vertices, &src_soln_gdofs[0] ) );
391 col_gdofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
392 col_dtoc_dofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
393 src_soln_gdofs.resize( m_remapper->m_covering_source_entities.size() * srcTagSize, UINT_MAX );
394 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_covering_source_entities, &src_soln_gdofs[0] ) );
397 #ifdef ALTERNATE_NUMBERING_IMPLEMENTATION
398 unsigned maxSrcIndx = 0;
401 std::vector< int > locdofs( srcTagSize );
402 std::map< Node, moab::EntityHandle > mapLocalMBNodes;
404 for(
unsigned iel = 0; iel < m_remapper->m_covering_source_entities.size(); ++iel )
406 EntityHandle eh = m_remapper->m_covering_source_entities[iel];
407 rval = m_interface->get_coords( &eh, 1, elcoords );
MB_CHK_ERR( rval );
408 Node elCentroid( elcoords[0], elcoords[1], elcoords[2] );
409 mapLocalMBNodes.insert( std::pair< Node, moab::EntityHandle >( elCentroid, eh ) );
412 const NodeVector& nodes = m_remapper->m_covering_source->nodes;
413 for(
unsigned j = 0; j < m_remapper->m_covering_source->faces.size(); j++ )
415 const Face&
face = m_remapper->m_covering_source->faces[j];
418 centroid.x = centroid.y = centroid.z = 0.0;
419 for(
unsigned l = 0; l <
face.edges.size(); ++l )
421 centroid.x += nodes[
face[l]].x;
422 centroid.y += nodes[
face[l]].y;
423 centroid.z += nodes[
face[l]].z;
425 const double factor = 1.0 /
face.edges.size();
426 centroid.x *= factor;
427 centroid.y *= factor;
428 centroid.z *= factor;
431 if( mapLocalMBNodes.find( centroid ) != mapLocalMBNodes.end() )
433 current_eh = mapLocalMBNodes[centroid];
436 rval = m_interface->tag_get_data( m_dofTagSrc, ¤t_eh, 1, &locdofs[0] );
MB_CHK_ERR( rval );
437 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
439 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
441 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
442 const int offsetDOF = p * m_nDofsPEl_Src + q;
443 maxSrcIndx = ( localDOF > maxSrcIndx ? localDOF : maxSrcIndx );
444 std::cout <<
"Col: " << current_eh <<
", " << m_remapper->lid_to_gid_covsrc[j] <<
", " << offsetDOF
445 <<
", " << localDOF <<
", " << locdofs[offsetDOF] - 1 <<
", " << maxSrcIndx <<
"\n";
451 m_nTotDofs_SrcCov = 0;
452 if( srcdataGLLNodes ==
nullptr )
455 std::vector<int> sorted_tmp_dofs( src_soln_gdofs.size() );
456 std::copy( src_soln_gdofs.begin(), src_soln_gdofs.end(), sorted_tmp_dofs.begin() );
457 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
459 for(
unsigned i = 0; i < col_gdofmap.size(); ++i )
461 auto gdof = sorted_tmp_dofs[i];
464 col_gdofmap[i] = gdof - 1;
465 col_dtoc_dofmap[i] = i;
466 if( vprint ) std::cout <<
"Col: " << i <<
", " << col_gdofmap[i] <<
"\n";
472 if( isSrcContinuous )
473 dgll_cgll_covcol_ldofmap.resize( m_remapper->m_covering_source_entities.size() * srcTagSize,
false );
475 for(
unsigned j = 0; j < m_remapper->m_covering_source_entities.size(); j++ )
477 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
479 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
481 const int localDOF = ( *srcdataGLLNodes )[p][q][j] - 1;
482 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
483 if( isSrcContinuous && !dgll_cgll_covcol_ldofmap[localDOF] )
486 dgll_cgll_covcol_ldofmap[localDOF] =
true;
488 if( !isSrcContinuous ) m_nTotDofs_SrcCov++;
489 assert( src_soln_gdofs[offsetDOF] > 0 );
490 col_gdofmap[localDOF] = src_soln_gdofs[offsetDOF] - 1;
491 col_dtoc_dofmap[offsetDOF] = localDOF;
493 std::cout <<
"Col: " << offsetDOF <<
", " << localDOF <<
", " << col_gdofmap[offsetDOF] <<
", "
494 << m_nTotDofs_SrcCov <<
"\n";
500 if( m_remapper->point_cloud_source )
502 assert( m_nDofsPEl_Src == 1 );
503 srccol_gdofmap.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
504 srccol_dtoc_dofmap.resize( m_remapper->m_covering_source_vertices.size(), UINT_MAX );
505 locsrc_soln_gdofs.resize( m_remapper->m_source_vertices.size(), UINT_MAX );
506 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_vertices, &locsrc_soln_gdofs[0] ) );
510 srccol_gdofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
511 srccol_dtoc_dofmap.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
512 locsrc_soln_gdofs.resize( m_remapper->m_source_entities.size() * srcTagSize, UINT_MAX );
513 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagSrc, m_remapper->m_source_entities, &locsrc_soln_gdofs[0] ) );
518 if( srcdataGLLNodesSrc ==
nullptr )
521 std::vector< int > sorted_tmp_dofs( locsrc_soln_gdofs.size() );
522 std::copy( locsrc_soln_gdofs.begin(), locsrc_soln_gdofs.end(), sorted_tmp_dofs.begin() );
523 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
525 for(
unsigned i = 0; i < srccol_gdofmap.size(); ++i )
527 auto gdof = sorted_tmp_dofs[i];
529 srccol_gdofmap[i] = gdof - 1;
530 srccol_dtoc_dofmap[i] = i;
536 if( isSrcContinuous ) dgll_cgll_col_ldofmap.resize( m_remapper->m_source_entities.size() * srcTagSize,
false );
538 for(
unsigned j = 0; j < m_remapper->m_source_entities.size(); j++ )
540 for(
int p = 0; p < m_nDofsPEl_Src; p++ )
542 for(
int q = 0; q < m_nDofsPEl_Src; q++ )
544 const int localDOF = ( *srcdataGLLNodesSrc )[p][q][j] - 1;
545 const int offsetDOF = j * srcTagSize + p * m_nDofsPEl_Src + q;
546 if( isSrcContinuous && !dgll_cgll_col_ldofmap[localDOF] )
549 dgll_cgll_col_ldofmap[localDOF] =
true;
551 if( !isSrcContinuous ) m_nTotDofs_Src++;
552 assert( locsrc_soln_gdofs[offsetDOF] > 0 );
553 srccol_gdofmap[localDOF] = locsrc_soln_gdofs[offsetDOF] - 1;
554 srccol_dtoc_dofmap[offsetDOF] = localDOF;
560 int tgtTagSize = ( m_eOutputType == DiscretizationType_FV ? 1 : m_nDofsPEl_Dest * m_nDofsPEl_Dest );
561 if( m_remapper->point_cloud_target )
563 assert( m_nDofsPEl_Dest == 1 );
564 row_gdofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
565 row_dtoc_dofmap.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
566 tgt_soln_gdofs.resize( m_remapper->m_target_vertices.size(), UINT_MAX );
567 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_vertices, &tgt_soln_gdofs[0] ) );
572 row_gdofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
573 row_dtoc_dofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
574 tgt_soln_gdofs.resize( m_remapper->m_target_entities.size() * tgtTagSize, UINT_MAX );
575 MB_CHK_ERR( m_interface->tag_get_data( m_dofTagDest, m_remapper->m_target_entities, &tgt_soln_gdofs[0] ) );
581 if( tgtdataGLLNodes ==
nullptr )
584 std::vector< int > sorted_tmp_dofs( tgt_soln_gdofs.size() );
585 std::copy( tgt_soln_gdofs.begin(), tgt_soln_gdofs.end(), sorted_tmp_dofs.begin() );
586 std::sort( sorted_tmp_dofs.begin(), sorted_tmp_dofs.end() );
588 for(
unsigned i = 0; i < row_gdofmap.size(); ++i )
590 auto gdof = sorted_tmp_dofs[i];
592 row_gdofmap[i] = gdof - 1;
593 row_dtoc_dofmap[i] = i;
594 if( vprint ) std::cout <<
"Row: " << i <<
", " << row_gdofmap[i] <<
"\n";
600 if( isTgtContinuous ) dgll_cgll_row_ldofmap.resize( m_remapper->m_target_entities.size() * tgtTagSize,
false );
602 for(
unsigned j = 0; j < m_remapper->m_target_entities.size(); j++ )
604 for(
int p = 0; p < m_nDofsPEl_Dest; p++ )
606 for(
int q = 0; q < m_nDofsPEl_Dest; q++ )
608 const int localDOF = ( *tgtdataGLLNodes )[p][q][j] - 1;
609 const int offsetDOF = j * tgtTagSize + p * m_nDofsPEl_Dest + q;
610 if( isTgtContinuous && !dgll_cgll_row_ldofmap[localDOF] )
613 dgll_cgll_row_ldofmap[localDOF] =
true;
615 if( !isTgtContinuous ) m_nTotDofs_Dest++;
616 assert( tgt_soln_gdofs[offsetDOF] > 0 );
617 row_gdofmap[localDOF] = tgt_soln_gdofs[offsetDOF] - 1;
618 row_dtoc_dofmap[offsetDOF] = localDOF;
620 std::cout <<
"Row: " << offsetDOF <<
", " << localDOF <<
", " << row_gdofmap[offsetDOF] <<
", "
621 << m_nTotDofs_Dest <<
"\n";
628 #if defined( MOAB_HAVE_EIGEN3 ) && defined( VERBOSE )
631 std::cout <<
"[" << rank <<
"]" <<
"DoFs: row = " << m_nTotDofs_Dest <<
", " << row_gdofmap.size()
632 <<
", col = " << m_nTotDofs_Src <<
", " << m_nTotDofs_SrcCov <<
", " << col_gdofmap.size() <<
"\n";
638 #ifdef CHECK_INCREASING_DOF
639 for(
size_t i = 0; i < row_gdofmap.size() - 1; i++ )
641 if( row_gdofmap[i] > row_gdofmap[i + 1] )
642 std::cout <<
" on rank " << rank <<
" in row_gdofmap[" << i <<
"]=" << row_gdofmap[i] <<
" > row_gdofmap["
643 << i + 1 <<
"]=" << row_gdofmap[i + 1] <<
" \n";
645 for(
size_t i = 0; i < col_gdofmap.size() - 1; i++ )
647 if( col_gdofmap[i] > col_gdofmap[i + 1] )
648 std::cout <<
" on rank " << rank <<
" in col_gdofmap[" << i <<
"]=" << col_gdofmap[i] <<
" > col_gdofmap["
649 << i + 1 <<
"]=" << col_gdofmap[i + 1] <<
" \n";
665 col_dtoc_dofmap.resize( values_entities.size(), -1 );
666 for(
size_t j = 0; j < values_entities.size(); j++ )
669 const auto it = colMap.find( values_entities[j] - 1 );
670 if( it != colMap.end() ) col_dtoc_dofmap[j] = it->second;
679 row_dtoc_dofmap.resize( values_entities.size(), -1 );
680 for(
size_t j = 0; j < values_entities.size(); j++ )
683 const auto it = rowMap.find( values_entities[j] - 1 );
684 if( it != rowMap.end() ) row_dtoc_dofmap[j] = it->second;
691 std::string strOutputType,
692 const GenerateOfflineMapAlgorithmOptions& mapOptions,
693 const std::string& srcDofTagName,
694 const std::string& tgtDofTagName )
696 NcError
error( NcError::silent_nonfatal );
699 dbgprint.set_prefix(
"[TempestOnlineMap]: " );
702 const bool m_bPointCloudSource = ( m_remapper->point_cloud_source );
703 const bool m_bPointCloudTarget = ( m_remapper->point_cloud_target );
704 const bool m_bPointCloud = m_bPointCloudSource || m_bPointCloudTarget;
713 STLStringHelper::ToLower( strInputType );
714 STLStringHelper::ToLower( strOutputType );
719 if( strInputType ==
"fv" )
721 eInputType = DiscretizationType_FV;
723 else if( strInputType ==
"cgll" )
725 eInputType = DiscretizationType_CGLL;
727 else if( strInputType ==
"dgll" )
729 eInputType = DiscretizationType_DGLL;
731 else if( strInputType ==
"pcloud" )
733 eInputType = DiscretizationType_PCLOUD;
737 _EXCEPTION1(
"Invalid \"in_type\" value (%s), expected [fv|cgll|dgll]", strInputType.c_str() );
740 if( strOutputType ==
"fv" )
742 eOutputType = DiscretizationType_FV;
744 else if( strOutputType ==
"cgll" )
746 eOutputType = DiscretizationType_CGLL;
748 else if( strOutputType ==
"dgll" )
750 eOutputType = DiscretizationType_DGLL;
752 else if( strOutputType ==
"pcloud" )
754 eOutputType = DiscretizationType_PCLOUD;
758 _EXCEPTION1(
"Invalid \"out_type\" value (%s), expected [fv|cgll|dgll]", strOutputType.c_str() );
762 m_bConserved = !mapOptions.fNoConservation;
763 m_eInputType = eInputType;
764 m_eOutputType = eOutputType;
767 std::string strMapAlgorithm(
"" );
768 int nMonotoneType = ( mapOptions.fMonotone ) ? ( 1 ) : ( 0 );
771 std::set< std::string > setMethodStrings;
774 for(
size_t i = 0; i <= mapOptions.strMethod.length(); i++ )
776 if( ( i == mapOptions.strMethod.length() ) || ( mapOptions.strMethod[i] ==
';' ) )
778 std::string strMethodString = mapOptions.strMethod.substr( iLast, i - iLast );
779 STLStringHelper::RemoveWhitespaceInPlace( strMethodString );
780 if( strMethodString.length() > 0 )
782 setMethodStrings.insert( strMethodString );
789 for(
auto it : setMethodStrings )
794 if( nMonotoneType != 0 )
796 _EXCEPTIONT(
"Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
798 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
800 _EXCEPTIONT(
"--method \"mono2\" is only used when remapping to/from CGLL or DGLL grids" );
806 else if( it ==
"mono3" )
808 if( nMonotoneType != 0 )
810 _EXCEPTIONT(
"Multiple monotonicity specifications found (--mono) or (--method \"mono#\")" );
812 if( ( m_eInputType == DiscretizationType_FV ) && ( m_eOutputType == DiscretizationType_FV ) )
814 _EXCEPTIONT(
"--method \"mono3\" is only used when remapping to/from CGLL or DGLL grids" );
820 else if( it ==
"volumetric" )
822 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType == DiscretizationType_FV ) )
824 _EXCEPTIONT(
"--method \"volumetric\" may only be used for FV->CGLL or FV->DGLL remapping" );
826 strMapAlgorithm =
"volumetric";
830 else if( it ==
"invdist" )
832 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
834 _EXCEPTIONT(
"--method \"invdist\" may only be used for FV->FV remapping" );
836 strMapAlgorithm =
"invdist";
840 else if( it ==
"delaunay" )
842 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
844 _EXCEPTIONT(
"--method \"delaunay\" may only be used for FV->FV remapping" );
846 strMapAlgorithm =
"delaunay";
850 else if( it ==
"bilin" )
852 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
854 _EXCEPTIONT(
"--method \"bilin\" may only be used for FV->FV remapping" );
856 strMapAlgorithm =
"fvbilin";
860 else if( it ==
"intbilin" )
862 if( m_eOutputType != DiscretizationType_FV )
864 _EXCEPTIONT(
"--method \"intbilin\" may only be used when mapping to FV." );
866 if( m_eInputType == DiscretizationType_FV )
868 strMapAlgorithm =
"fvintbilin";
872 strMapAlgorithm =
"mono3";
877 else if( it ==
"intbilingb" )
879 if( ( m_eInputType != DiscretizationType_FV ) || ( m_eOutputType != DiscretizationType_FV ) )
881 _EXCEPTIONT(
"--method \"intbilingb\" may only be used for FV->FV remapping" );
883 strMapAlgorithm =
"fvintbilingb";
887 _EXCEPTION1(
"Invalid --method argument \"%s\"", it.c_str() );
892 ( m_eInputType == DiscretizationType_FV || m_eInputType == DiscretizationType_PCLOUD ? 1
895 ( m_eOutputType == DiscretizationType_FV || m_eOutputType == DiscretizationType_PCLOUD ? 1
896 : mapOptions.nPout );
899 MB_CHK_ERR( SetDOFmapTags( srcDofTagName, tgtDofTagName ) );
903 rval = m_interface->tag_get_handle(
"aream", 1,
MB_TYPE_DOUBLE, areaTag,
907 if( is_root )
dbgprint.printf( 0,
"aream tag already defined \n" );
910 double local_areas[3] = { 0.0, 0.0, 0.0 }, global_areas[3] = { 0.0, 0.0, 0.0 };
911 if( !m_bPointCloudSource )
914 if( is_root )
dbgprint.printf( 0,
"Calculating input mesh Face areas\n" );
915 local_areas[0] = m_meshInput->CalculateFaceAreas( mapOptions.fSourceConcave );
917 MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_source_entities, m_meshInput->vecFaceArea ) );
920 m_meshInputCov->CalculateFaceAreas( mapOptions.fSourceConcave );
923 if( !m_bPointCloudTarget )
926 if( is_root )
dbgprint.printf( 0,
"Calculating output mesh Face areas\n" );
927 local_areas[1] = m_meshOutput->CalculateFaceAreas( mapOptions.fTargetConcave );
929 MB_CHK_ERR( m_interface->tag_set_data( areaTag, m_remapper->m_target_entities, m_meshOutput->vecFaceArea ) );
935 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->vecTargetFaceIx.size() );
938 if( is_root )
dbgprint.printf( 0,
"Calculating overlap mesh Face areas\n" );
940 m_meshOverlap->CalculateFaceAreas( mapOptions.fSourceConcave || mapOptions.fTargetConcave );
943 std::copy( local_areas, local_areas + 3, global_areas );
946 if( m_pcomm && is_parallel )
947 MPI_Reduce( local_areas, global_areas, 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
951 dbgprint.printf( 0,
"Input Mesh Geometric Area: %1.15e\n", global_areas[0] );
952 dbgprint.printf( 0,
"Output Mesh Geometric Area: %1.15e\n", global_areas[1] );
953 dbgprint.printf( 0,
"Overlap Mesh Recovered Area: %1.15e\n", global_areas[2] );
957 constexpr
bool fCorrectAreas =
true;
960 if( is_root )
dbgprint.printf( 0,
"Correcting source/target areas to overlap mesh areas\n" );
961 DataArray1D< double > dSourceArea( m_meshInputCov->faces.size() );
962 DataArray1D< double > dTargetArea( m_meshOutput->faces.size() );
964 assert( m_meshOverlap->vecSourceFaceIx.size() == m_meshOverlap->faces.size() );
965 assert( m_meshOverlap->vecTargetFaceIx.size() == m_meshOverlap->faces.size() );
966 assert( m_meshOverlap->vecFaceArea.GetRows() == m_meshOverlap->faces.size() );
968 assert( m_meshInputCov->vecFaceArea.GetRows() == m_meshInputCov->faces.size() );
969 assert( m_meshOutput->vecFaceArea.GetRows() == m_meshOutput->faces.size() );
971 for(
size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
973 if( m_meshOverlap->vecSourceFaceIx[i] < 0 || m_meshOverlap->vecTargetFaceIx[i] < 0 )
977 assert(
static_cast< size_t >( m_meshOverlap->vecSourceFaceIx[i] ) < m_meshInputCov->faces.size() );
978 dSourceArea[m_meshOverlap->vecSourceFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
979 assert(
static_cast< size_t >( m_meshOverlap->vecTargetFaceIx[i] ) < m_meshOutput->faces.size() );
980 dTargetArea[m_meshOverlap->vecTargetFaceIx[i]] += m_meshOverlap->vecFaceArea[i];
983 for(
size_t i = 0; i < m_meshInputCov->faces.size(); i++ )
985 if( fabs( dSourceArea[i] - m_meshInputCov->vecFaceArea[i] ) < 1.0e-10 )
987 m_meshInputCov->vecFaceArea[i] = dSourceArea[i];
990 for(
size_t i = 0; i < m_meshOutput->faces.size(); i++ )
992 if( fabs( dTargetArea[i] - m_meshOutput->vecFaceArea[i] ) < 1.0e-10 )
994 m_meshOutput->vecFaceArea[i] = dTargetArea[i];
1000 if( !m_bPointCloudSource && eInputType == DiscretizationType_FV )
1002 this->SetSourceAreas( m_meshInputCov->vecFaceArea );
1003 if( m_meshInputCov->vecMask.size() )
1005 this->SetSourceMask( m_meshInputCov->vecMask );
1010 if( !m_bPointCloudTarget && eOutputType == DiscretizationType_FV )
1012 this->SetTargetAreas( m_meshOutput->vecFaceArea );
1013 if( m_meshOutput->vecMask.size() )
1015 this->SetTargetMask( m_meshOutput->vecMask );
1031 if( ( eInputType == DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1034 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1035 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap(
false );
1038 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1039 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1041 this->m_pdataGLLNodesIn =
nullptr;
1042 this->m_pdataGLLNodesOut =
nullptr;
1045 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin,
false,
nullptr,
nullptr, eOutputType,
1046 mapOptions.nPout,
false,
nullptr );
MB_CHK_ERR( rval );
1049 if( is_root )
dbgprint.printf( 0,
"Calculating remap weights\n" );
1052 if( strMapAlgorithm ==
"invdist" )
1054 if( is_root )
dbgprint.printf( 0,
"Calculating map (invdist)\n" );
1055 if( m_meshInputCov->faces.size() )
1056 LinearRemapFVtoFVInvDist( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *
this );
1058 else if( strMapAlgorithm ==
"delaunay" )
1060 if( is_root )
dbgprint.printf( 0,
"Calculating map (delaunay)\n" );
1061 if( m_meshInputCov->faces.size() )
1062 LinearRemapTriangulation( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *
this );
1064 else if( strMapAlgorithm ==
"fvintbilin" )
1066 if( is_root )
dbgprint.printf( 0,
"Calculating map (intbilin)\n" );
1067 if( m_meshInputCov->faces.size() )
1068 LinearRemapIntegratedBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *
this );
1070 else if( strMapAlgorithm ==
"fvintbilingb" )
1072 if( is_root )
dbgprint.printf( 0,
"Calculating map (intbilingb)\n" );
1073 if( m_meshInputCov->faces.size() )
1074 LinearRemapIntegratedGeneralizedBarycentric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1077 else if( strMapAlgorithm ==
"fvbilin" )
1082 m_meshInputCov->Write(
"SourceMeshMBTR.g" );
1083 m_meshOutput->Write(
"TargetMeshMBTR.g" );
1087 m_meshInputCov->Write(
"SourceMeshMBTR" + std::to_string( rank ) +
".g" );
1088 m_meshOutput->Write(
"TargetMeshMBTR" + std::to_string( rank ) +
".g" );
1091 if( is_root )
dbgprint.printf( 0,
"Calculating map (bilin)\n" );
1092 if( m_meshInputCov->faces.size() )
1093 LinearRemapBilinear( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, *
this );
1097 if( is_root )
dbgprint.printf( 0,
"Calculating conservative FV-FV map\n" );
1098 if( m_meshInputCov->faces.size() )
1100 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1101 LinearRemapFVtoFV( *m_meshInputCov, *m_meshOutput, *m_meshOverlap,
1102 ( mapOptions.fMonotone ) ? ( 1 ) : ( mapOptions.nPin ), *
this );
1104 LinearRemapFVtoFV_Tempest_MOAB( ( mapOptions.fMonotone ? 1 : mapOptions.nPin ) );
1109 else if( eInputType == DiscretizationType_FV )
1111 DataArray3D< double > dataGLLJacobian;
1113 if( is_root )
dbgprint.printf( 0,
"Generating output mesh meta data\n" );
1114 double dNumericalArea_loc = GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble,
1115 dataGLLNodesDest, dataGLLJacobian );
1117 double dNumericalArea = dNumericalArea_loc;
1118 #ifdef MOAB_HAVE_MPI
1120 MPI_Reduce( &dNumericalArea_loc, &dNumericalArea, 1, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
1122 if( is_root )
dbgprint.printf( 0,
"Output Mesh Numerical Area: %1.15e\n", dNumericalArea );
1125 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1126 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1128 this->m_pdataGLLNodesIn =
nullptr;
1129 this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1132 bool fContinuous = ( eOutputType == DiscretizationType_CGLL );
1134 if( eOutputType == DiscretizationType_CGLL )
1136 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobian, this->GetTargetAreas() );
1140 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetTargetAreas() );
1144 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1145 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap(
false );
1148 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin,
false,
nullptr,
nullptr, eOutputType,
1149 mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1153 if( strMapAlgorithm ==
"volumetric" )
1155 if( is_root )
dbgprint.printf( 0,
"Calculating remapping weights for FV->GLL (volumetric)\n" );
1156 LinearRemapFVtoGLL_Volumetric( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest,
1157 dataGLLJacobian, this->GetTargetAreas(), mapOptions.nPin, *
this,
1158 nMonotoneType, fContinuous, mapOptions.fNoConservation );
1162 if( is_root )
dbgprint.printf( 0,
"Calculating remapping weights for FV->GLL\n" );
1163 LinearRemapFVtoGLL( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesDest, dataGLLJacobian,
1164 this->GetTargetAreas(), mapOptions.nPin, *
this, nMonotoneType, fContinuous,
1165 mapOptions.fNoConservation );
1168 else if( ( eInputType == DiscretizationType_PCLOUD ) || ( eOutputType == DiscretizationType_PCLOUD ) )
1170 DataArray3D< double > dataGLLJacobian;
1171 if( !m_bPointCloudSource )
1174 if( m_meshInputCov->revnodearray.size() == 0 ) m_meshInputCov->ConstructReverseNodeArray();
1175 if( m_meshInputCov->edgemap.size() == 0 ) m_meshInputCov->ConstructEdgeMap(
false );
1178 if( eInputType == DiscretizationType_FV )
1180 this->InitializeSourceCoordinatesFromMeshFV( *m_meshInputCov );
1184 if( is_root )
dbgprint.printf( 0,
"Generating input mesh meta data\n" );
1185 DataArray3D< double > dataGLLJacobianSrc;
1186 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1188 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1189 dataGLLJacobianSrc );
1194 if( !m_bPointCloudTarget )
1197 if( m_meshOutput->revnodearray.size() == 0 ) m_meshOutput->ConstructReverseNodeArray();
1198 if( m_meshOutput->edgemap.size() == 0 ) m_meshOutput->ConstructEdgeMap(
false );
1201 if( eOutputType == DiscretizationType_FV )
1203 this->InitializeSourceCoordinatesFromMeshFV( *m_meshOutput );
1207 if( is_root )
dbgprint.printf( 0,
"Generating output mesh meta data\n" );
1208 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1215 rval = this->SetDOFmapAssociation(
1216 eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1217 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ?
nullptr : &dataGLLNodesSrcCov ),
1218 ( m_bPointCloudSource || eInputType == DiscretizationType_FV ?
nullptr : &dataGLLNodesSrc ),
1219 eOutputType, mapOptions.nPout, ( eOutputType == DiscretizationType_CGLL ),
1220 ( m_bPointCloudTarget ?
nullptr : &dataGLLNodesDest ) );
MB_CHK_ERR( rval );
1223 if( is_root )
dbgprint.printf( 0,
"Calculating remap weights with Nearest-Neighbor method\n" );
1224 rval = LinearRemapNN_MOAB(
true ,
false );
MB_CHK_ERR( rval );
1226 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType == DiscretizationType_FV ) )
1228 DataArray3D< double > dataGLLJacobianSrc, dataGLLJacobian;
1230 if( is_root )
dbgprint.printf( 0,
"Generating input mesh meta data\n" );
1232 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1233 dataGLLJacobianSrc );
1234 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1237 if( dataGLLNodesSrcCov.GetSubColumns() != m_meshInputCov->faces.size() )
1239 _EXCEPTIONT(
"Number of element does not match between metadata and "
1244 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1245 this->InitializeTargetCoordinatesFromMeshFV( *m_meshOutput );
1248 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1250 if( eInputType == DiscretizationType_CGLL )
1252 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobian, this->GetSourceAreas() );
1256 GenerateDiscontinuousJacobian( dataGLLJacobian, this->GetSourceAreas() );
1260 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1261 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1265 if( is_root )
dbgprint.printf( 0,
"Calculating remap weights\n" );
1267 if( strMapAlgorithm ==
"volumetric" )
1269 _EXCEPTIONT(
"Unimplemented: Volumetric currently unavailable for"
1273 this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1274 this->m_pdataGLLNodesOut =
nullptr;
1276 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1277 LinearRemapSE4( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov, dataGLLJacobian,
1278 nMonotoneType, fContinuousIn, mapOptions.fNoConservation, mapOptions.fSparseConstraints,
1281 LinearRemapSE4_Tempest_MOAB( dataGLLNodesSrcCov, dataGLLJacobian, nMonotoneType, fContinuousIn,
1282 mapOptions.fNoConservation );
1285 else if( ( eInputType != DiscretizationType_FV ) && ( eOutputType != DiscretizationType_FV ) )
1287 DataArray3D< double > dataGLLJacobianIn, dataGLLJacobianSrc;
1288 DataArray3D< double > dataGLLJacobianOut;
1291 if( is_root )
dbgprint.printf( 0,
"Generating input mesh meta data\n" );
1293 GenerateMetaData( *m_meshInput, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrc,
1294 dataGLLJacobianSrc );
1296 GenerateMetaData( *m_meshInputCov, mapOptions.nPin, mapOptions.fNoBubble, dataGLLNodesSrcCov,
1297 dataGLLJacobianIn );
1299 if( is_root )
dbgprint.printf( 0,
"Generating output mesh meta data\n" );
1300 GenerateMetaData( *m_meshOutput, mapOptions.nPout, mapOptions.fNoBubble, dataGLLNodesDest,
1301 dataGLLJacobianOut );
1304 this->InitializeSourceCoordinatesFromMeshFE( *m_meshInputCov, mapOptions.nPin, dataGLLNodesSrcCov );
1305 this->InitializeTargetCoordinatesFromMeshFE( *m_meshOutput, mapOptions.nPout, dataGLLNodesDest );
1308 bool fContinuousIn = ( eInputType == DiscretizationType_CGLL );
1310 if( eInputType == DiscretizationType_CGLL )
1312 GenerateUniqueJacobian( dataGLLNodesSrcCov, dataGLLJacobianIn, this->GetSourceAreas() );
1316 GenerateDiscontinuousJacobian( dataGLLJacobianIn, this->GetSourceAreas() );
1320 bool fContinuousOut = ( eOutputType == DiscretizationType_CGLL );
1322 if( eOutputType == DiscretizationType_CGLL )
1324 GenerateUniqueJacobian( dataGLLNodesDest, dataGLLJacobianOut, this->GetTargetAreas() );
1328 GenerateDiscontinuousJacobian( dataGLLJacobianOut, this->GetTargetAreas() );
1332 rval = this->SetDOFmapAssociation( eInputType, mapOptions.nPin, ( eInputType == DiscretizationType_CGLL ),
1333 &dataGLLNodesSrcCov, &dataGLLNodesSrc, eOutputType, mapOptions.nPout,
1334 ( eOutputType == DiscretizationType_CGLL ), &dataGLLNodesDest );
MB_CHK_ERR( rval );
1336 this->m_pdataGLLNodesIn = &dataGLLNodesSrcCov;
1337 this->m_pdataGLLNodesOut = &dataGLLNodesDest;
1340 if( is_root )
dbgprint.printf( 0,
"Calculating remap weights\n" );
1342 #ifdef USE_NATIVE_TEMPESTREMAP_ROUTINES
1343 LinearRemapGLLtoGLL_Integrated( *m_meshInputCov, *m_meshOutput, *m_meshOverlap, dataGLLNodesSrcCov,
1344 dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1345 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1346 fContinuousIn, fContinuousOut, mapOptions.fSparseConstraints, *
this );
1348 LinearRemapGLLtoGLL2_MOAB( dataGLLNodesSrcCov, dataGLLJacobianIn, dataGLLNodesDest, dataGLLJacobianOut,
1349 this->GetTargetAreas(), mapOptions.nPin, mapOptions.nPout, nMonotoneType,
1350 fContinuousIn, fContinuousOut, mapOptions.fNoConservation );
1355 _EXCEPTIONT(
"Not implemented" );
1358 #ifdef MOAB_HAVE_EIGEN3
1359 copy_tempest_sparsemat_to_eigen3();
1362 #ifdef MOAB_HAVE_MPI
1366 rval = m_remapper->GetOverlapAugmentedEntities( ghostedEnts );
MB_CHK_ERR( rval );
1368 rval = m_interface->remove_entities( m_meshOverlapSet, ghostedEnts );
MB_CHK_SET_ERR( rval,
"Deleting ghosted entities failed" );
1372 if( !mapOptions.fNoCheck )
1374 if( is_root )
dbgprint.printf( 0,
"Verifying map" );
1375 this->IsConsistent( 1.0e-8 );
1376 if( !mapOptions.fNoConservation ) this->IsConservative( 1.0e-8 );
1378 if( nMonotoneType != 0 )
1380 this->IsMonotone( 1.0e-12 );
1384 catch( Exception& e )
1386 dbgprint.printf( 0,
"%s", e.ToString().c_str() );
1387 return ( moab::MB_FAILURE );
1391 return ( moab::MB_FAILURE );
1400 #ifndef MOAB_HAVE_MPI
1402 return OfflineMap::IsConsistent( dTolerance );
1407 DataArray1D< int > dataRows;
1408 DataArray1D< int > dataCols;
1409 DataArray1D< double > dataEntries;
1412 DataArray1D< double > dRowSums;
1413 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1414 dRowSums.Allocate( m_mapRemap.GetRows() );
1416 for(
unsigned i = 0; i < dataRows.GetRows(); i++ )
1418 dRowSums[dataRows[i]] += dataEntries[i];
1422 int fConsistent = 0;
1423 for(
unsigned i = 0; i < dRowSums.GetRows(); i++ )
1425 if( fabs( dRowSums[i] - 1.0 ) > dTolerance )
1428 int rowGID = row_gdofmap[i];
1429 Announce(
"TempestOnlineMap is not consistent in row %i (%1.15e)", rowGID, dRowSums[i] );
1434 int fConsistentGlobal = 0;
1435 ierr = MPI_Allreduce( &fConsistent, &fConsistentGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1436 if( ierr != MPI_SUCCESS )
return -1;
1438 return fConsistentGlobal;
1446 #ifndef MOAB_HAVE_MPI
1448 return OfflineMap::IsConservative( dTolerance );
1455 DataArray1D< int > dataRows;
1456 DataArray1D< int > dataCols;
1457 DataArray1D< double > dataEntries;
1458 const DataArray1D< double >& dTargetAreas = this->GetTargetAreas();
1459 const DataArray1D< double >& dSourceAreas = this->GetSourceAreas();
1462 std::vector< int > dColumnsUnique;
1463 std::vector< double > dColumnSums;
1465 int nColumns = m_mapRemap.GetColumns();
1466 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1467 dColumnSums.resize( m_nTotDofs_SrcCov, 0.0 );
1468 dColumnsUnique.resize( m_nTotDofs_SrcCov, -1 );
1470 for(
unsigned i = 0; i < dataEntries.GetRows(); i++ )
1472 dColumnSums[dataCols[i]] += dataEntries[i] * dTargetAreas[dataRows[i]] / dSourceAreas[dataCols[i]];
1474 assert( dataCols[i] < m_nTotDofs_SrcCov );
1477 int colGID = this->GetColGlobalDoF( dataCols[i] );
1479 dColumnsUnique[dataCols[i]] = colGID;
1486 std::vector< int > nElementsInProc;
1487 const int nDATA = 3;
1488 nElementsInProc.resize(
size * nDATA );
1489 int senddata[nDATA] = { nColumns, m_nTotDofs_SrcCov, m_nTotDofs_Src };
1490 ierr = MPI_Gather( senddata, nDATA, MPI_INT, nElementsInProc.data(), nDATA, MPI_INT, rootProc, m_pcomm->comm() );
1491 if( ierr != MPI_SUCCESS )
return -1;
1493 int nTotVals = 0, nTotColumns = 0, nTotColumnsUnq = 0;
1494 std::vector< int > dColumnIndices;
1495 std::vector< double > dColumnSourceAreas;
1496 std::vector< double > dColumnSumsTotal;
1497 std::vector< int > displs, rcount;
1498 if( rank == rootProc )
1500 displs.resize(
size + 1, 0 );
1501 rcount.resize(
size, 0 );
1503 for(
int ir = 0; ir <
size; ++ir )
1505 nTotVals += nElementsInProc[ir * nDATA];
1506 nTotColumns += nElementsInProc[ir * nDATA + 1];
1507 nTotColumnsUnq += nElementsInProc[ir * nDATA + 2];
1510 rcount[ir] = nElementsInProc[ir * nDATA + 1];
1516 printf(
"Total nnz: %d, global source elements = %d\n", nTotVals, gsum );
1518 dColumnIndices.resize( nTotColumns, -1 );
1519 dColumnSumsTotal.resize( nTotColumns, 0.0 );
1528 ierr = MPI_Gatherv( &dColumnsUnique[0], m_nTotDofs_SrcCov, MPI_INT, &dColumnIndices[0], rcount.data(),
1529 displs.data(), MPI_INT, rootProc, m_pcomm->comm() );
1530 if( ierr != MPI_SUCCESS )
return -1;
1531 ierr = MPI_Gatherv( &dColumnSums[0], m_nTotDofs_SrcCov, MPI_DOUBLE, &dColumnSumsTotal[0], rcount.data(),
1532 displs.data(), MPI_DOUBLE, rootProc, m_pcomm->comm() );
1533 if( ierr != MPI_SUCCESS )
return -1;
1539 dColumnSums.clear();
1540 dColumnsUnique.clear();
1543 int fConservative = 0;
1544 if( rank == rootProc )
1546 displs[
size] = ( nTotColumns );
1548 std::map< int, double > dColumnSumsOnRoot;
1550 for(
int ir = 0; ir <
size; ir++ )
1552 for(
int ips = displs[ir]; ips < displs[ir + 1]; ips++ )
1554 if( dColumnIndices[ips] < 0 )
continue;
1556 assert( dColumnIndices[ips] < nTotColumnsUnq );
1557 dColumnSumsOnRoot[dColumnIndices[ips]] += dColumnSumsTotal[ips];
1563 for( std::map< int, double >::iterator it = dColumnSumsOnRoot.begin(); it != dColumnSumsOnRoot.end(); ++it )
1566 if( fabs( it->second - 1.0 ) > dTolerance )
1569 Announce(
"TempestOnlineMap is not conservative in column "
1572 it->first, it->second );
1578 ierr = MPI_Bcast( &fConservative, 1, MPI_INT, rootProc, m_pcomm->comm() );
1579 if( ierr != MPI_SUCCESS )
return -1;
1581 return fConservative;
1589 #ifndef MOAB_HAVE_MPI
1591 return OfflineMap::IsMonotone( dTolerance );
1596 DataArray1D< int > dataRows;
1597 DataArray1D< int > dataCols;
1598 DataArray1D< double > dataEntries;
1600 m_mapRemap.GetEntries( dataRows, dataCols, dataEntries );
1604 for(
unsigned i = 0; i < dataRows.GetRows(); i++ )
1606 if( ( dataEntries[i] < -dTolerance ) || ( dataEntries[i] > 1.0 + dTolerance ) )
1610 Announce(
"TempestOnlineMap is not monotone in entry (%i): %1.15e", i, dataEntries[i] );
1615 int fMonotoneGlobal = 0;
1616 ierr = MPI_Allreduce( &fMonotone, &fMonotoneGlobal, 1, MPI_INT, MPI_SUM, m_pcomm->comm() );
1617 if( ierr != MPI_SUCCESS )
return -1;
1619 return fMonotoneGlobal;
1628 bool useMOABAdjacencies,
1631 assert( nrings > 0 );
1632 assert( useMOABAdjacencies || trMesh !=
nullptr );
1634 const size_t nrows = vecAdjFaces.size();
1636 for(
size_t index = 0; index < nrows; index++ )
1638 vecAdjFaces[index].insert( index );
1641 if( useMOABAdjacencies )
1651 int adjIndex =
entities.index( *it );
1653 if( adjIndex >= 0 ) vecAdjFaces[index].insert( adjIndex );
1662 Face&
face = trMesh->faces[index];
1663 GetAdjacentFaceVectorByEdge( *trMesh, index, nrings *
face.edges.size(), adjFaces );
1666 for(
auto adjFace : adjFaces )
1667 if( adjFace.first >= 0 )
1668 vecAdjFaces[index].insert( adjFace.first );
1678 double default_projection )
1680 std::vector< double > solSTagVals;
1681 std::vector< double > solTTagVals;
1684 if( m_remapper->point_cloud_source || m_remapper->point_cloud_target )
1686 if( m_remapper->point_cloud_source )
1689 solSTagVals.resize( covSrcEnts.
size(), default_projection );
1695 solSTagVals.resize( covSrcEnts.
size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1696 default_projection );
1699 if( m_remapper->point_cloud_target )
1702 solTTagVals.resize( tgtEnts.
size(), default_projection );
1708 solTTagVals.resize( tgtEnts.
size() * this->GetDestinationNDofsPerElement() *
1709 this->GetDestinationNDofsPerElement(),
1710 default_projection );
1718 solSTagVals.resize( covSrcEnts.
size() * this->GetSourceNDofsPerElement() * this->GetSourceNDofsPerElement(),
1719 default_projection );
1720 solTTagVals.resize( tgtEnts.
size() * this->GetDestinationNDofsPerElement() *
1721 this->GetDestinationNDofsPerElement(),
1722 default_projection );
1729 MB_CHK_SET_ERR( m_interface->tag_get_data( srcSolutionTag, sents, &solSTagVals[0] ),
1730 "Getting local tag data failed" );
1735 MB_CHK_SET_ERR( this->ApplyWeights( solSTagVals, solTTagVals, transpose ),
1736 "Applying remap operator onto source vector data failed" );
1739 MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1740 "Setting target tag data failed" );
1742 if( caasType != CAAS_NONE )
1744 std::string tgtSolutionTagName;
1745 MB_CHK_SET_ERR( m_interface->tag_get_name( tgtSolutionTag, tgtSolutionTagName ),
"Getting tag name failed" );
1748 constexpr
int nmax_caas_iterations = 10;
1749 double mismatch = 1.0;
1750 int caasIteration = 0;
1751 double initialMismatch = 0.0;
1752 while( ( fabs( mismatch / initialMismatch ) > 1e-15 && fabs( mismatch ) > 1e-15 ) &&
1753 caasIteration++ < nmax_caas_iterations )
1755 double dMassDiffPostGlobal;
1756 std::pair< double, double > mDefect =
1757 this->ApplyBoundsLimiting( solSTagVals, solTTagVals, caasType, caasIteration, mismatch );
1758 #ifdef MOAB_HAVE_MPI
1759 double dMassDiffPost = mDefect.second;
1760 MPI_Allreduce( &dMassDiffPost, &dMassDiffPostGlobal, 1, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
1762 dMassDiffPostGlobal = mDefect.second;
1764 if( caasIteration == 1 ) initialMismatch = mDefect.first;
1765 if( m_remapper->verbose && is_root )
1767 printf(
"Field {%s} -> CAAS iteration: %d, mass defect: %3.4e, post-CAAS: %3.4e\n",
1768 tgtSolutionTagName.c_str(), caasIteration, mDefect.first, dMassDiffPostGlobal );
1770 mismatch = dMassDiffPostGlobal;
1773 MB_CHK_SET_ERR( m_interface->tag_set_data( tgtSolutionTag, tents, &solTTagVals[0] ),
1774 "Setting local tag data failed" );
1782 const std::string& solnName,
1784 sample_function testFunction,
1786 std::string cloneSolnName )
1789 const bool outputEnabled = ( is_root );
1799 trmesh = m_remapper->m_covering_source;
1800 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
1801 : m_remapper->m_covering_source_entities );
1802 discOrder = m_nDofsPEl_Src;
1803 discMethod = m_eInputType;
1808 trmesh = m_remapper->m_target;
1810 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
1811 discOrder = m_nDofsPEl_Dest;
1812 discMethod = m_eOutputType;
1817 std::cout <<
"Invalid context specified for defining an analytical solution tag" << std::endl;
1818 return moab::MB_FAILURE;
1823 rval = m_interface->tag_get_handle( solnName.c_str(), discOrder * discOrder,
MB_TYPE_DOUBLE, solnTag,
1825 if( clonedSolnTag !=
nullptr )
1827 if( cloneSolnName.size() == 0 )
1829 cloneSolnName = solnName + std::string(
"Cloned" );
1831 rval = m_interface->tag_get_handle( cloneSolnName.c_str(), discOrder * discOrder,
MB_TYPE_DOUBLE,
1836 const int TriQuadratureOrder = 10;
1838 if( outputEnabled ) std::cout <<
"Using triangular quadrature of order " << TriQuadratureOrder << std::endl;
1840 TriangularQuadratureRule triquadrule( TriQuadratureOrder );
1842 const int TriQuadraturePoints = triquadrule.GetPoints();
1844 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1845 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1848 DataArray1D< double > dVar;
1849 DataArray1D< double > dVarMB;
1852 DataArray1D< double > dNodeArea;
1857 if( discMethod == DiscretizationType_CGLL || discMethod == DiscretizationType_DGLL )
1860 const bool fGLL =
true;
1861 const bool fGLLIntegrate =
false;
1864 DataArray3D< int > dataGLLNodes;
1865 DataArray3D< double > dataGLLJacobian;
1867 GenerateMetaData( *trmesh, discOrder,
false, dataGLLNodes, dataGLLJacobian );
1870 int nElements = trmesh->faces.size();
1873 for(
int k = 0; k < nElements; k++ )
1875 const Face&
face = trmesh->faces[k];
1877 if(
face.edges.size() != 4 )
1879 _EXCEPTIONT(
"Non-quadrilateral face detected; "
1880 "incompatible with --gll" );
1886 for(
int i = 0; i < discOrder; i++ )
1888 for(
int j = 0; j < discOrder; j++ )
1890 for(
int k = 0; k < nElements; k++ )
1892 if( dataGLLNodes[i][j][k] > iMaxNode )
1894 iMaxNode = dataGLLNodes[i][j][k];
1901 DataArray1D< double > dG;
1902 DataArray1D< double > dW;
1904 GaussLobattoQuadrature::GetPoints( discOrder, 0.0, 1.0, dG, dW );
1907 const int nGaussP = 10;
1909 DataArray1D< double > dGaussG;
1910 DataArray1D< double > dGaussW;
1912 GaussQuadrature::GetPoints( nGaussP, 0.0, 1.0, dGaussG, dGaussW );
1915 dVar.Allocate( iMaxNode );
1916 dVarMB.Allocate( discOrder * discOrder * nElements );
1917 dNodeArea.Allocate( iMaxNode );
1920 for(
int k = 0; k < nElements; k++ )
1922 const Face&
face = trmesh->faces[k];
1927 for(
int i = 0; i < discOrder; i++ )
1929 for(
int j = 0; j < discOrder; j++ )
1937 ApplyLocalMap(
face, trmesh->nodes, dG[i], dG[j], node, dDx1G, dDx2G );
1940 double dNodeLon = atan2( node.y, node.x );
1941 if( dNodeLon < 0.0 )
1943 dNodeLon += 2.0 * M_PI;
1945 double dNodeLat = asin( node.z );
1947 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1949 dVar[dataGLLNodes[j][i][k] - 1] = dSample;
1956 DataArray2D< double > dCoeff( discOrder, discOrder );
1958 for(
int p = 0; p < nGaussP; p++ )
1960 for(
int q = 0; q < nGaussP; q++ )
1968 ApplyLocalMap(
face, trmesh->nodes, dGaussG[p], dGaussG[q], node, dDx1G, dDx2G );
1971 Node nodeCross = CrossProduct( dDx1G, dDx2G );
1974 sqrt( nodeCross.x * nodeCross.x + nodeCross.y * nodeCross.y + nodeCross.z * nodeCross.z );
1978 SampleGLLFiniteElement( 0, discOrder, dGaussG[p], dGaussG[q], dCoeff );
1981 double dNodeLon = atan2( node.y, node.x );
1982 if( dNodeLon < 0.0 )
1984 dNodeLon += 2.0 * M_PI;
1986 double dNodeLat = asin( node.z );
1988 double dSample = ( *testFunction )( dNodeLon, dNodeLat );
1991 for(
int i = 0; i < discOrder; i++ )
1993 for(
int j = 0; j < discOrder; j++ )
1996 double dNodalArea = dCoeff[i][j] * dGaussW[p] * dGaussW[q] * dJacobian;
1998 dVar[dataGLLNodes[i][j][k] - 1] += dSample * dNodalArea;
2000 dNodeArea[dataGLLNodes[i][j][k] - 1] += dNodalArea;
2011 for(
size_t i = 0; i < dVar.GetRows(); i++ )
2013 dVar[i] /= dNodeArea[i];
2020 for(
unsigned j = 0; j <
entities.size(); j++ )
2021 for(
int p = 0; p < discOrder; p++ )
2022 for(
int q = 0; q < discOrder; q++ )
2024 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2025 dVarMB[offsetDOF] = dVar[col_dtoc_dofmap[offsetDOF]];
2030 for(
unsigned j = 0; j <
entities.size(); j++ )
2031 for(
int p = 0; p < discOrder; p++ )
2032 for(
int q = 0; q < discOrder; q++ )
2034 const int offsetDOF = j * discOrder * discOrder + p * discOrder + q;
2035 dVarMB[offsetDOF] = dVar[row_dtoc_dofmap[offsetDOF]];
2040 rval = m_interface->tag_set_data( solnTag,
entities, &dVarMB[0] );
MB_CHK_ERR( rval );
2045 if( discMethod == DiscretizationType_FV )
2050 dVar.Allocate( trmesh->faces.size() );
2052 std::vector< Node >& nodes = trmesh->nodes;
2055 for(
size_t i = 0; i < trmesh->faces.size(); i++ )
2057 const Face&
face = trmesh->faces[i];
2060 for(
size_t j = 0; j <
face.edges.size() - 2; j++ )
2063 const Node& node0 = nodes[
face[0]];
2064 const Node& node1 = nodes[
face[j + 1]];
2065 const Node& node2 = nodes[
face[j + 2]];
2069 faceTri.SetNode( 0,
face[0] );
2070 faceTri.SetNode( 1,
face[j + 1] );
2071 faceTri.SetNode( 2,
face[j + 2] );
2073 double dTriangleArea = CalculateFaceArea( faceTri, nodes );
2076 double dTotalSample = 0.0;
2079 for(
int k = 0; k < TriQuadraturePoints; k++ )
2081 Node node( TriQuadratureG[k][0] * node0.x + TriQuadratureG[k][1] * node1.x +
2082 TriQuadratureG[k][2] * node2.x,
2083 TriQuadratureG[k][0] * node0.y + TriQuadratureG[k][1] * node1.y +
2084 TriQuadratureG[k][2] * node2.y,
2085 TriQuadratureG[k][0] * node0.z + TriQuadratureG[k][1] * node1.z +
2086 TriQuadratureG[k][2] * node2.z );
2088 double dMagnitude = node.Magnitude();
2089 node.x /= dMagnitude;
2090 node.y /= dMagnitude;
2091 node.z /= dMagnitude;
2093 double dLon = atan2( node.y, node.x );
2098 double dLat = asin( node.z );
2100 double dSample = ( *testFunction )( dLon, dLat );
2102 dTotalSample += dSample * TriQuadratureW[k] * dTriangleArea;
2105 dVar[i] += dTotalSample / trmesh->vecFaceArea[i];
2113 std::vector< Node >& nodes = trmesh->nodes;
2116 dVar.Allocate( nodes.size() );
2118 for(
size_t j = 0; j < nodes.size(); j++ )
2120 Node& node = nodes[j];
2121 double dMagnitude = node.Magnitude();
2122 node.x /= dMagnitude;
2123 node.y /= dMagnitude;
2124 node.z /= dMagnitude;
2125 double dLon = atan2( node.y, node.x );
2130 double dLat = asin( node.z );
2132 double dSample = ( *testFunction )( dLon, dLat );
2146 std::map< std::string, double >& metrics,
2150 const bool outputEnabled = ( is_root );
2161 entities = ( m_remapper->point_cloud_source ? m_remapper->m_covering_source_vertices
2162 : m_remapper->m_covering_source_entities );
2163 discOrder = m_nDofsPEl_Src;
2171 ( m_remapper->point_cloud_target ? m_remapper->m_target_vertices : m_remapper->m_target_entities );
2172 discOrder = m_nDofsPEl_Dest;
2178 std::cout <<
"Invalid context specified for defining an analytical solution tag" << std::endl;
2179 return moab::MB_FAILURE;
2184 std::string exactTagName, projTagName;
2185 const int ntotsize =
entities.size() * discOrder * discOrder;
2186 std::vector< double > exactSolution( ntotsize, 0.0 ), projSolution( ntotsize, 0.0 );
2187 rval = m_interface->tag_get_name( exactTag, exactTagName );
MB_CHK_ERR( rval );
2188 rval = m_interface->tag_get_data( exactTag,
entities, &exactSolution[0] );
MB_CHK_ERR( rval );
2189 rval = m_interface->tag_get_name( approxTag, projTagName );
MB_CHK_ERR( rval );
2190 rval = m_interface->tag_get_data( approxTag,
entities, &projSolution[0] );
MB_CHK_ERR( rval );
2192 const auto& ovents = m_remapper->m_overlap_entities;
2194 std::vector< double > errnorms( 4, 0.0 ), globerrnorms( 4, 0.0 );
2195 double sumarea = 0.0;
2196 for(
size_t i = 0; i < ovents.size(); ++i )
2198 const int srcidx = m_remapper->m_overlap->vecSourceFaceIx[i];
2199 if( srcidx < 0 )
continue;
2200 const int tgtidx = m_remapper->m_overlap->vecTargetFaceIx[i];
2201 if( tgtidx < 0 )
continue;
2202 const double ovarea = m_remapper->m_overlap->vecFaceArea[i];
2203 const double error = fabs( exactSolution[tgtidx] - projSolution[tgtidx] );
2204 errnorms[0] += ovarea *
error;
2206 errnorms[3] = (
error > errnorms[3] ?
error : errnorms[3] );
2209 errnorms[2] = sumarea;
2210 #ifdef MOAB_HAVE_MPI
2213 MPI_Reduce( &errnorms[0], &globerrnorms[0], 3, MPI_DOUBLE, MPI_SUM, 0, m_pcomm->comm() );
2214 MPI_Reduce( &errnorms[3], &globerrnorms[3], 1, MPI_DOUBLE, MPI_MAX, 0, m_pcomm->comm() );
2217 for(
int i = 0; i < 4; ++i )
2218 globerrnorms[i] = errnorms[i];
2221 globerrnorms[0] = ( globerrnorms[0] / globerrnorms[2] );
2222 globerrnorms[1] = std::sqrt( globerrnorms[1] / globerrnorms[2] );
2225 metrics[
"L1Error"] = globerrnorms[0];
2226 metrics[
"L2Error"] = globerrnorms[1];
2227 metrics[
"LinfError"] = globerrnorms[3];
2231 std::cout <<
"Error metrics when comparing " << projTagName <<
" against " << exactTagName << std::endl;
2232 std::cout <<
"\t Total Intersection area = " << globerrnorms[2] << std::endl;
2233 std::cout <<
"\t L_1 error = " << globerrnorms[0] << std::endl;
2234 std::cout <<
"\t L_2 error = " << globerrnorms[1] << std::endl;
2235 std::cout <<
"\t L_inf error = " << globerrnorms[3] << std::endl;