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