 |
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 /// \file TempestLinearRemap.cpp
4 /// \author Vijay Mahadevan
5 /// \version Mar 08, 2017
6 ///
7
8 #ifdef WIN32 /* windows */
9 #define _USE_MATH_DEFINES // For M_PI
10 #endif
11
12 #pragma GCC diagnostic push
13 #pragma GCC diagnostic ignored "-Wdeprecated-copy-with-user-provided-copy"
14 #pragma GCC diagnostic ignored "-Wsign-compare"
15
16 #include "Announce.h"
17 #include "DataArray3D.h"
18 #include "FiniteElementTools.h"
19 #include "FiniteVolumeTools.h"
20 #include "GaussLobattoQuadrature.h"
21 #include "TriangularQuadrature.h"
22 #include "MathHelper.h"
23 #include "SparseMatrix.h"
24 #include "OverlapMesh.h"
25
26 #include "DebugOutput.hpp"
27 #include "moab/AdaptiveKDTree.hpp"
28
29 #include "moab/Remapping/TempestOnlineMap.hpp"
30 #include "moab/TupleList.hpp"
31 #include "moab/MeshTopoUtil.hpp"
32
33 #pragma GCC diagnostic pop
34
35 #include <fstream>
36 #include <cmath>
37 #include <cstdlib>
38 #include <sstream>
39 #include <numeric>
40 #include <algorithm>
41 #include <unordered_set>
42
43 // #define VERBOSE
44 #define USE_ComputeAdjacencyRelations
45
46 /// <summary>
47 /// Face index and distance metric pair.
48 /// </summary>
49 typedef std::pair< int, int > FaceDistancePair;
50
51 /// <summary>
52 /// Vector storing adjacent Faces.
53 /// </summary>
54 typedef std::vector< FaceDistancePair > AdjacentFaceVector;
55
56 ///////////////////////////////////////////////////////////////////////////////
57
58 moab::ErrorCode moab::TempestOnlineMap::LinearRemapNN_MOAB( bool use_GID_matching, bool strict_check )
59 {
60 /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov) */
61
62 #ifdef VVERBOSE
63 {
64 std::ofstream output_file( "rowcolindices.txt", std::ios::out );
65 output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
66 << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
67 output_file << "Rows \n";
68 for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
69 output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
70 output_file << "Cols \n";
71 for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
72 output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
73 output_file.flush(); // required here
74 output_file.close();
75 }
76 #endif
77
78 if( use_GID_matching )
79 {
80 std::map< unsigned, unsigned > src_gl;
81 for( unsigned it = 0; it < col_gdofmap.size(); ++it )
82 src_gl[col_gdofmap[it]] = it;
83
84 std::map< unsigned, unsigned >::iterator iter;
85 for( unsigned it = 0; it < row_gdofmap.size(); ++it )
86 {
87 unsigned row = row_gdofmap[it];
88 iter = src_gl.find( row );
89 if( strict_check && iter == src_gl.end() )
90 {
91 std::cout << "Searching for global target DOF " << row
92 << " but could not find correspondence in source mesh.\n";
93 assert( false );
94 }
95 else if( iter == src_gl.end() )
96 {
97 continue;
98 }
99 else
100 {
101 unsigned icol = src_gl[row];
102 unsigned irow = it;
103
104 // Set the permutation matrix in local space
105 m_mapRemap( irow, icol ) = 1.0;
106 }
107 }
108
109 return moab::MB_SUCCESS;
110 }
111 else
112 {
113 /* Create a Kd-tree to perform local queries to find nearest neighbors */
114
115 return moab::MB_FAILURE;
116 }
117 }
118
119 ///////////////////////////////////////////////////////////////////////////////
120
121 void moab::TempestOnlineMap::LinearRemapFVtoFV_Tempest_MOAB( int nOrder )
122 {
123 // Order of triangular quadrature rule
124 const int TriQuadRuleOrder = 4;
125
126 // Verify ReverseNodeArray has been calculated
127 if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
128 {
129 _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInputCov" );
130 }
131
132 // Triangular quadrature rule
133 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
134
135 // Number of coefficients needed at this order
136 #ifdef RECTANGULAR_TRUNCATION
137 int nCoefficients = nOrder * nOrder;
138 #endif
139 #ifdef TRIANGULAR_TRUNCATION
140 int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
141 #endif
142
143 // Number of faces you need
144 const int nRequiredFaceSetSize = nCoefficients;
145
146 // Fit weight exponent
147 const int nFitWeightsExponent = nOrder + 2;
148
149 // Announcements
150 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
151 dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
152 if( is_root )
153 {
154 dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
155 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
156 dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
157 dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
158 dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
159 }
160
161 // Current overlap face
162 int ixOverlap = 0;
163 #ifdef VERBOSE
164 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
165 #endif
166 DataArray2D< double > dIntArray;
167 DataArray1D< double > dConstraint( nCoefficients );
168
169 // Loop through all faces on m_meshInputCov
170 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
171 {
172 // Output every 1000 elements
173 #ifdef VERBOSE
174 if( ixFirst % outputFrequency == 0 && is_root )
175 {
176 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
177 }
178 #endif
179 // Find the set of Faces that overlap faceFirst
180 int ixOverlapBegin = ixOverlap;
181 unsigned ixOverlapEnd = ixOverlapBegin;
182
183 for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
184 {
185 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 ) break;
186 }
187
188 unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
189
190 if( nOverlapFaces == 0 ) continue;
191
192 // Build integration array
193 BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
194 nOrder, dIntArray );
195
196 // Set of Faces to use in building the reconstruction and associated
197 // distance metric.
198 AdjacentFaceVector vecAdjFaces;
199
200 GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
201
202 // Number of adjacent Faces
203 int nAdjFaces = vecAdjFaces.size();
204
205 // Determine the conservative constraint equation
206 double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
207 dConstraint.Zero();
208 for( int p = 0; p < nCoefficients; p++ )
209 {
210 for( unsigned j = 0; j < nOverlapFaces; j++ )
211 {
212 dConstraint[p] += dIntArray[p][j];
213 }
214 dConstraint[p] /= dFirstArea;
215 }
216
217 // Build the fit array from the integration operator
218 DataArray2D< double > dFitArray;
219 DataArray1D< double > dFitWeights;
220 DataArray2D< double > dFitArrayPlus;
221
222 BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
223 dFitArray, dFitWeights );
224
225 // Compute the inverse fit array
226 bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
227
228 // Multiply integration array and fit array
229 DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
230 if( fSuccess )
231 {
232 // Multiply integration array and inverse fit array
233 for( int i = 0; i < nAdjFaces; i++ )
234 {
235 for( size_t j = 0; j < nOverlapFaces; j++ )
236 {
237 for( int k = 0; k < nCoefficients; k++ )
238 {
239 dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
240 }
241 }
242 }
243
244 // Unable to invert fit array, drop to 1st order. In this case
245 // dFitArrayPlus(0,0) = 1 and all other entries are zero.
246 }
247 else
248 {
249 dComposedArray.Zero();
250 for( size_t j = 0; j < nOverlapFaces; j++ )
251 {
252 dComposedArray( 0, j ) += dIntArray( 0, j );
253 }
254 }
255
256 // Put composed array into map
257 for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
258 {
259 for( unsigned j = 0; j < nOverlapFaces; j++ )
260 {
261 int& ixFirstFaceLoc = vecAdjFaces[i].first;
262 int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
263 // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
264 // ixFirstFaceLoc); int ixSecondFaceGlob =
265 // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
266
267 // signal to not participate, because it is a ghost target
268 if( ixSecondFaceLoc < 0 ) continue; // do not do anything
269
270 m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
271 dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
272 }
273 }
274
275 // Increment the current overlap index
276 ixOverlap += nOverlapFaces;
277 }
278
279 return;
280 }
281
282 ///////////////////////////////////////////////////////////////////////////////
283
284 void moab::TempestOnlineMap::PrintMapStatistics()
285 {
286 int nrows = m_weightMatrix.rows(); // Number of rows
287 int ncols = m_weightMatrix.cols(); // Number of columns
288 int NNZ = m_weightMatrix.nonZeros(); // Number of non zero values
289 #ifdef MOAB_HAVE_MPI
290 // find out min/max for NNZ, ncols, nrows
291 // should work on std c++ 11
292 int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
293 int rarr3[6];
294 MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
295
296 int total[3];
297 MPI_Reduce( arr3, total, 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
298 if( !rank )
299 std::cout << "-> Rows (min/max/sum): (" << rarr3[1] << " / " << -rarr3[4] << " / " << total[1] << "), "
300 << " Cols (min/max/sum): (" << rarr3[2] << " / " << -rarr3[5] << " / " << total[2] << "), "
301 << " NNZ (min/max/sum): (" << rarr3[0] << " / " << -rarr3[3] << " / " << total[0] << ")\n";
302 #else
303 std::cout << "-> Rows: " << nrows << ", Cols: " << ncols << ", NNZ: " << NNZ << "\n";
304 #endif
305 }
306
307 #ifdef MOAB_HAVE_EIGEN3
308 void moab::TempestOnlineMap::copy_tempest_sparsemat_to_eigen3()
309 {
310 #ifndef VERBOSE
311 #define VERBOSE_ACTIVATED
312 // #define VERBOSE
313 #endif
314
315 /* Should the columns be the global size of the matrix ? */
316 m_weightMatrix.resize( m_nTotDofs_Dest, m_nTotDofs_SrcCov );
317 m_rowVector.resize( m_weightMatrix.rows() );
318 m_colVector.resize( m_weightMatrix.cols() );
319
320 #ifdef VERBOSE
321 int locrows = std::max( m_mapRemap.GetRows(), m_nTotDofs_Dest );
322 int loccols = std::max( m_mapRemap.GetColumns(), m_nTotDofs_SrcCov );
323
324 std::cout << m_weightMatrix.rows() << ", " << locrows << ", " << m_weightMatrix.cols() << ", " << loccols << "\n";
325 // assert(m_weightMatrix.rows() == locrows && m_weightMatrix.cols() == loccols);
326 #endif
327
328 DataArray1D< int > lrows;
329 DataArray1D< int > lcols;
330 DataArray1D< double > lvals;
331 m_mapRemap.GetEntries( lrows, lcols, lvals );
332 size_t locvals = lvals.GetRows();
333
334 // first matrix
335 typedef Eigen::Triplet< double > Triplet;
336 std::vector< Triplet > tripletList;
337 tripletList.reserve( locvals );
338 for( size_t iv = 0; iv < locvals; iv++ )
339 {
340 tripletList.push_back( Triplet( lrows[iv], lcols[iv], lvals[iv] ) );
341 }
342 m_weightMatrix.setFromTriplets( tripletList.begin(), tripletList.end() );
343 m_weightMatrix.makeCompressed();
344
345 #ifdef VERBOSE
346 std::stringstream sstr;
347 sstr << "tempestmatrix.txt.0000" << rank;
348 std::ofstream output_file( sstr.str(), std::ios::out );
349 output_file << "0 " << locrows << " 0 " << loccols << "\n";
350 for( unsigned iv = 0; iv < locvals; iv++ )
351 {
352 // output_file << lrows[iv] << " " << row_ldofmap[lrows[iv]] << " " <<
353 // row_gdofmap[row_ldofmap[lrows[iv]]] << " " << col_gdofmap[col_ldofmap[lcols[iv]]] << " "
354 // << lvals[iv] << "\n";
355 output_file << row_gdofmap[row_ldofmap[lrows[iv]]] << " " << col_gdofmap[col_ldofmap[lcols[iv]]] << " "
356 << lvals[iv] << "\n";
357 }
358 output_file.flush(); // required here
359 output_file.close();
360 #endif
361
362 #ifdef VERBOSE_ACTIVATED
363 #undef VERBOSE_ACTIVATED
364 #undef VERBOSE
365 #endif
366 return;
367 }
368
369 ///////////////////////////////////////////////////////////////////////////////
370
371 template < typename T >
372 static std::vector< size_t > sort_indexes( const std::vector< T >& v )
373 {
374 // initialize original index locations
375 std::vector< size_t > idx( v.size() );
376 std::iota( idx.begin(), idx.end(), 0 );
377
378 // sort indexes based on comparing values in v
379 // using std::stable_sort instead of std::sort
380 // to avoid unnecessary index re-orderings
381 // when v contains elements of equal values
382 std::stable_sort( idx.begin(), idx.end(), [&v]( size_t i1, size_t i2 ) { return fabs( v[i1] ) > fabs( v[i2] ); } );
383
384 return idx;
385 }
386
387 double moab::TempestOnlineMap::QLTLimiter( int caasIteration,
388 std::vector< double >& dataCorrectedField,
389 std::vector< double >& dataLowerBound,
390 std::vector< double >& dataUpperBound,
391 std::vector< double >& dMassDefect )
392 {
393 const size_t nrows = dataCorrectedField.size();
394 double dMassL = 0.0;
395 double dMassU = 0.0;
396 std::vector< double > dataCorrection( nrows );
397 double dMassDiffCum = 0.0;
398 double dLMinusU = fabs( dataUpperBound[0] - dataLowerBound[0] );
399 const DataArray1D< double >& dTargetAreas = this->m_remapper->m_target->vecFaceArea;
400
401 // std::vector< size_t > sortedIdx = sort_indexes( dMassDefect );
402 std::vector< std::unordered_set< int > > vecAdjTargetFaces( nrows );
403 constexpr bool useMOABAdjacencies = true;
404 #ifdef USE_ComputeAdjacencyRelations
405 if( useMOABAdjacencies )
406 ComputeAdjacencyRelations( vecAdjTargetFaces, caasIteration, m_remapper->m_target_entities,
407 useMOABAdjacencies );
408 else
409 ComputeAdjacencyRelations( vecAdjTargetFaces, caasIteration, m_remapper->m_target_entities, useMOABAdjacencies,
410 this->m_remapper->m_target );
411 #else
412 moab::MeshTopoUtil mtu( m_interface );
413 ;
414 #endif
415
416 for( size_t i = 0; i < nrows; i++ )
417 {
418 // size_t index = sortedIdx[i];
419 size_t index = i;
420 dataCorrection[index] = fmax( dataLowerBound[index], fmin( dataUpperBound[index], 0.0 ) );
421 // dMassDiff[index] = dMassDefect[index] - dTargetAreas[index] * dataCorrection[index];
422 // dMassDiff[index] = dMassDefect[index];
423
424 dMassL += dTargetAreas[index] * dataLowerBound[index];
425 dMassU += dTargetAreas[index] * dataUpperBound[index];
426 dLMinusU = fmax( dLMinusU, fabs( dataUpperBound[index] - dataLowerBound[index] ) );
427 dMassDiffCum += dMassDefect[index] - dTargetAreas[index] * dataCorrection[index];
428
429 #ifndef USE_ComputeAdjacencyRelations
430 vecAdjTargetFaces[index].insert( index ); // add self target face first
431 {
432 // Compute the adjacent faces to the target face
433 if( useMOABAdjacencies )
434 {
435 moab::Range ents;
436 // ents.insert( m_remapper->m_target_entities.index( m_remapper->m_target_entities[index] ) );
437 ents.insert( m_remapper->m_target_entities[index] );
438 moab::Range adjEnts;
439 moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, caasIteration );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
440 for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
441 {
442 // int adjIndex = m_interface->id_from_handle(*it)-1;
443 int adjIndex = m_remapper->m_target_entities.index( *it );
444 // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
445 if( adjIndex >= 0 ) vecAdjTargetFaces[index].insert( adjIndex );
446 }
447 }
448 else
449 {
450 AdjacentFaceVector vecAdjFaces;
451 GetAdjacentFaceVectorByEdge( *this->m_remapper->m_target, index,
452 ( m_output_order + 1 ) * ( m_output_order + 1 ) * ( m_output_order + 1 ),
453 // ( m_output_order + 1 ) * ( m_output_order + 1 ),
454 // ( 4 ) * ( m_output_order + 1 ) * ( m_output_order + 1 ),
455 vecAdjFaces );
456
457 // Add the adjacent faces to the target face list
458 for( auto adjFace : vecAdjFaces )
459 if( adjFace.first >= 0 )
460 vecAdjTargetFaces[index].insert( adjFace.first ); // map target face to source face
461 }
462 }
463 #endif
464 }
465
466 #ifdef MOAB_HAVE_MPI
467 std::vector< double > localDefects( 5, 0.0 ), globalDefects( 5, 0.0 );
468 localDefects[0] = dMassL;
469 localDefects[1] = dMassU;
470 localDefects[2] = dMassDiffCum;
471 localDefects[3] = dLMinusU;
472 // localDefects[4] = dMassCorrectU;
473
474 MPI_Allreduce( localDefects.data(), globalDefects.data(), 4, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
475
476 dMassL = globalDefects[0];
477 dMassU = globalDefects[1];
478 dMassDiffCum = globalDefects[2];
479 dLMinusU = globalDefects[3];
480 // dMassCorrectU = globalDefects[4];
481 #endif
482
483 //If the upper and lower bounds are too close together, just clip
484 if( fabs( dMassDiffCum ) < 1e-15 || dLMinusU < 1e-15 )
485 {
486 for( size_t i = 0; i < nrows; i++ )
487 dataCorrectedField[i] += dataCorrection[i];
488 return dMassDiffCum;
489 }
490 else
491 {
492 if( dMassL > dMassDiffCum )
493 {
494 Announce( "Lower bound mass exceeds target mass by %1.15e: CAAS will need another iteration",
495 dMassL - dMassDiffCum );
496 dMassDiffCum = dMassL;
497 // dMass -= dMassL;
498 }
499 else if( dMassU < dMassDiffCum )
500 {
501 Announce( "Target mass exceeds upper bound mass by %1.15e: CAAS will need another iteration",
502 dMassDiffCum - dMassU );
503 dMassDiffCum = dMassU;
504 // dMass -= dMassU;
505 }
506
507 // TODO: optimize away dataMassVec by a simple transient double within the loop
508 // DataArray1D< double > dataMassVec( nrows ); //vector of mass redistribution
509 for( size_t i = 0; i < nrows; i++ )
510 {
511 // size_t index = sortedIdx[i];
512 size_t index = i;
513 const std::unordered_set< int >& neighbors = vecAdjTargetFaces[index];
514 if( dMassDefect[index] > 0.0 )
515 {
516 double dMassCorrectU = 0.0;
517 for( auto it : neighbors )
518 dMassCorrectU += dTargetAreas[it] * ( dataUpperBound[it] - dataCorrection[it] );
519
520 // double dMassDiffCumOld = dMassDefect[index];
521 for( auto it : neighbors )
522 dataCorrection[it] +=
523 dMassDefect[index] * ( dataUpperBound[it] - dataCorrection[it] ) / dMassCorrectU;
524 }
525 else
526 {
527 double dMassCorrectL = 0.0;
528 for( auto it : neighbors )
529 dMassCorrectL += dTargetAreas[it] * ( dataCorrection[it] - dataLowerBound[it] );
530
531 // double dMassDiffCumOld = dMassDefect[index];
532 for( auto it : neighbors )
533 dataCorrection[it] +=
534 dMassDefect[index] * ( dataCorrection[it] - dataLowerBound[it] ) / dMassCorrectL;
535 }
536 }
537
538 for( size_t i = 0; i < nrows; i++ )
539 dataCorrectedField[i] += dataCorrection[i];
540 }
541
542 return dMassDiffCum;
543 }
544
545 void moab::TempestOnlineMap::CAASLimiter( std::vector< double >& dataCorrectedField,
546 std::vector< double >& dataLowerBound,
547 std::vector< double >& dataUpperBound,
548 double& dMass )
549 {
550 const size_t nrows = dataCorrectedField.size();
551 double dMassL = 0.0;
552 double dMassU = 0.0;
553 std::vector< double > dataCorrection( nrows );
554 const DataArray1D< double >& dTargetAreas = this->m_remapper->m_target->vecFaceArea;
555 double dMassDiff = dMass;
556 double dLMinusU = fabs( dataUpperBound[0] - dataLowerBound[0] );
557 double dMassCorrectU = 0.0;
558 double dMassCorrectL = 0.0;
559 for( size_t i = 0; i < nrows; i++ )
560 {
561 dataCorrection[i] = fmax( dataLowerBound[i], fmin( dataUpperBound[i], 0.0 ) );
562 dMassL += dTargetAreas[i] * dataLowerBound[i];
563 dMassU += dTargetAreas[i] * dataUpperBound[i];
564 dMassDiff -= dTargetAreas[i] * dataCorrection[i];
565 dLMinusU = fmax( dLMinusU, fabs( dataUpperBound[i] - dataLowerBound[i] ) );
566 dMassCorrectL += dTargetAreas[i] * ( dataCorrection[i] - dataLowerBound[i] );
567 dMassCorrectU += dTargetAreas[i] * ( dataUpperBound[i] - dataCorrection[i] );
568 }
569
570 #ifdef MOAB_HAVE_MPI
571 std::vector< double > localDefects( 5, 0.0 ), globalDefects( 5, 0.0 );
572 localDefects[0] = dMassL;
573 localDefects[1] = dMassU;
574 localDefects[2] = dMassDiff;
575 localDefects[3] = dMassCorrectL;
576 localDefects[4] = dMassCorrectU;
577
578 MPI_Allreduce( localDefects.data(), globalDefects.data(), 5, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
579
580 dMassL = globalDefects[0];
581 dMassU = globalDefects[1];
582 dMassDiff = globalDefects[2];
583 dMassCorrectL = globalDefects[3];
584 dMassCorrectU = globalDefects[4];
585 #endif
586
587 //If the upper and lower bounds are too close together, just clip
588 if( fabs( dMassDiff ) < 1e-15 || fabs( dLMinusU ) < 1e-15 )
589 {
590 for( size_t i = 0; i < nrows; i++ )
591 dataCorrectedField[i] += dataCorrection[i];
592 return;
593 }
594 else
595 {
596 if( dMassL > dMassDiff )
597 {
598 Announce( "%d: Lower bound mass exceeds target mass by %1.15e: CAAS will need another iteration", rank,
599 dMassL - dMassDiff );
600 dMassDiff = dMassL;
601 dMass -= dMassL;
602 }
603 else if( dMassU < dMassDiff )
604 {
605 Announce( "%d: Target mass exceeds upper bound mass by %1.15e: CAAS will need another iteration", rank,
606 dMassDiff - dMassU );
607 dMassDiff = dMassU;
608 dMass -= dMassU;
609 }
610
611 // TODO: optimize away dataMassVec by a simple transient double within the loop
612 DataArray1D< double > dataMassVec( nrows ); //vector of mass redistribution
613 if( dMassDiff > 0.0 )
614 {
615 for( size_t i = 0; i < nrows; i++ )
616 {
617 dataMassVec[i] = ( dataUpperBound[i] - dataCorrection[i] ) / dMassCorrectU;
618 dataCorrection[i] += dMassDiff * dataMassVec[i];
619 }
620 }
621 else
622 {
623 for( size_t i = 0; i < nrows; i++ )
624 {
625 dataMassVec[i] = ( dataCorrection[i] - dataLowerBound[i] ) / dMassCorrectL;
626 dataCorrection[i] += dMassDiff * dataMassVec[i];
627 }
628 }
629
630 for( size_t i = 0; i < nrows; i++ )
631 dataCorrectedField[i] += dataCorrection[i];
632 }
633
634 return;
635 }
636
637 std::pair< double, double > moab::TempestOnlineMap::ApplyBoundsLimiting( std::vector< double >& dataInDouble,
638 std::vector< double >& dataOutDouble,
639 CAASType caasType,
640 int caasIteration,
641 double mismatch )
642 {
643 // Currently only implemented for FV to FV remapping
644 // We should generalize this to other types of remapping
645 assert( !dataGLLNodesSrcCov.IsAttached() && !dataGLLNodesDest.IsAttached() );
646
647 std::pair< double, double > massDefect( 0.0, 0.0 );
648
649 // Check if the source and target data are of the same size
650 const size_t nTargetCount = dataOutDouble.size();
651 const DataArray1D< double >& m_dOverlapAreas = this->m_remapper->m_overlap->vecFaceArea;
652
653 // Apply the offline map to the data
654 double dMassDiff = 0.0;
655 std::vector< double > x( nTargetCount );
656 std::vector< double > dataLowerBound( nTargetCount );
657 std::vector< double > dataUpperBound( nTargetCount );
658 std::vector< double > massVector( nTargetCount );
659 std::vector< std::unordered_set< int > > vecSourceOvTarget( nTargetCount );
660
661 #undef USE_ComputeAdjacencyRelations
662 constexpr bool useMOABAdjacencies = true;
663 #ifdef USE_ComputeAdjacencyRelations
664 // Compute the adjacent faces to the source face
665 // However, calling MOAB to do this does not work correctly as we need ixS to be the index
666 // Cannot just iterate over all entities in the source covering mesh
667 if( caasType == CAAS_QLT || caasType == CAAS_LOCAL_ADJACENT )
668 {
669 if( useMOABAdjacencies )
670 {
671 moab::ErrorCode rval =
672 ComputeAdjacencyRelations( vecSourceOvTarget, caasIteration, m_remapper->m_covering_source_entities,
673 useMOABAdjacencies );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
674 }
675 else
676 {
677 moab::ErrorCode rval =
678 ComputeAdjacencyRelations( vecSourceOvTarget, caasIteration, m_remapper->m_covering_source_entities,
679 useMOABAdjacencies, m_meshInputCov );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
680 }
681 }
682 #else
683 moab::MeshTopoUtil mtu( m_interface );
684 #endif
685
686 // Initialize the bounds on the given source and target data
687 double dSourceMin = dataInDouble[0];
688 double dSourceMax = dataInDouble[0];
689 double dTargetMin = dataOutDouble[0];
690 double dTargetMax = dataOutDouble[0];
691 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
692 {
693 const int ixS = m_meshOverlap->vecSourceFaceIx[i];
694 const int ixT = m_meshOverlap->vecTargetFaceIx[i];
695
696 if( ixT < 0 ) continue; // skip ghost target faces
697
698 assert( m_dOverlapAreas[i] > 0.0 );
699 assert( ixS >= 0 );
700 assert( ixT >= 0 );
701
702 #ifndef USE_ComputeAdjacencyRelations
703 // Compute the adjacent faces to the target face
704 vecSourceOvTarget[ixT].insert( ixS ); // map target face to source face
705 if( ( caasType == CAAS_QLT || caasType == CAAS_LOCAL_ADJACENT ) )
706 {
707 if( useMOABAdjacencies )
708 {
709 moab::Range ents;
710 ents.insert( m_remapper->m_covering_source_entities[ixS] );
711 moab::Range adjEnts;
712 moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, caasIteration );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
713 for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
714 {
715 int adjIndex = m_remapper->m_covering_source_entities.index( *it );
716 if( adjIndex >= 0 ) vecSourceOvTarget[ixT].insert( adjIndex );
717 }
718 }
719 else
720 {
721 // Compute the adjacent faces to the target face
722 AdjacentFaceVector vecAdjFaces;
723 GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixS,
724 ( caasIteration ) * ( m_input_order + 1 ) * ( m_input_order + 1 ),
725 vecAdjFaces );
726
727 //Compute min/max over neighboring faces
728 for( size_t iadj = 0; iadj < vecAdjFaces.size(); iadj++ )
729 vecSourceOvTarget[ixT].insert( vecAdjFaces[iadj].first ); // map target face to source face
730 }
731 }
732 #endif
733
734 // Update the min and max values of the source data
735 dSourceMax = fmax( dSourceMax, dataInDouble[ixS] );
736 dSourceMin = fmin( dSourceMin, dataInDouble[ixS] );
737
738 // Update the min and max values of the target data
739 dTargetMin = fmin( dTargetMin, dataOutDouble[ixT] );
740 dTargetMax = fmax( dTargetMax, dataOutDouble[ixT] );
741
742 const double locMassDiff = ( dataInDouble[ixS] * m_dOverlapAreas[i] ) - // source mass
743 ( dataOutDouble[ixT] * m_dOverlapAreas[i] ); // target mass
744
745 // Update the mass difference between source and target faces
746 // linked to the overlap mesh element
747 dMassDiff += locMassDiff; // target mass
748 massVector[ixT] += locMassDiff;
749 }
750
751 #ifdef MOAB_HAVE_MPI
752 std::vector< double > localMinMaxDefects( 5, 0.0 ), globalMinMaxDefects( 5, 0.0 );
753 localMinMaxDefects[0] = dSourceMin;
754 localMinMaxDefects[1] = dTargetMin;
755 localMinMaxDefects[2] = dSourceMax;
756 localMinMaxDefects[3] = dTargetMax;
757 localMinMaxDefects[4] = dMassDiff;
758
759 if( caasType == CAAS_GLOBAL )
760 {
761 MPI_Allreduce( localMinMaxDefects.data(), globalMinMaxDefects.data(), 2, MPI_DOUBLE, MPI_MIN, m_pcomm->comm() );
762 MPI_Allreduce( localMinMaxDefects.data() + 2, globalMinMaxDefects.data() + 2, 2, MPI_DOUBLE, MPI_MAX,
763 m_pcomm->comm() );
764 dSourceMin = globalMinMaxDefects[0];
765 dSourceMax = globalMinMaxDefects[2];
766 dTargetMin = globalMinMaxDefects[1];
767 dTargetMax = globalMinMaxDefects[3];
768 }
769 if( caasIteration == 1 )
770 MPI_Allreduce( localMinMaxDefects.data() + 4, globalMinMaxDefects.data() + 4, 1, MPI_DOUBLE, MPI_SUM,
771 m_pcomm->comm() );
772 else
773 globalMinMaxDefects[4] = mismatch;
774
775 dMassDiff = localMinMaxDefects[4];
776 // massDefect.first = localMinMaxDefects[4];
777 massDefect.first = globalMinMaxDefects[4];
778 #else
779
780 // massDefect.first = fabs( dMassDiff / ( dSourceMax - dSourceMin ) );
781 massDefect.first = dMassDiff;
782 #endif
783
784 // Early exit if the values are monotone already.
785 // if( ( dTargetMax <= dSourceMax && dTargetMin <= dSourceMin ) || fabs( massDefect.first ) < 1e-16 )
786 if( fabs( massDefect.first ) > 1e-20 )
787 {
788 if( caasType == CAAS_GLOBAL )
789 {
790 for( size_t i = 0; i < nTargetCount; i++ )
791 {
792 dataLowerBound[i] = dSourceMin - dataOutDouble[i];
793 dataUpperBound[i] = dSourceMax - dataOutDouble[i];
794 }
795 } // if( caasType == CAAS_GLOBAL )
796 else // caasType == CAAS_LOCAL
797 {
798 // Compute the local min and max values of the target data
799 std::vector< double > vecLocalUpperBound( nTargetCount );
800 std::vector< double > vecLocalLowerBound( nTargetCount );
801 // Loop over the target faces and compute the min and max values
802 // of the source data linked to the target faces
803 for( size_t i = 0; i < nTargetCount; i++ )
804 {
805 assert( vecSourceOvTarget[i].size() );
806
807 double dMinI = 1E10; // dataInDouble[vecSourceOvTarget[i][0]];
808 double dMaxI = -1E10; // dataInDouble[vecSourceOvTarget[i][0]];
809
810 // Compute max over intersecting source faces
811 for( const auto& srcElem : vecSourceOvTarget[i] )
812 {
813 dMinI = fmin( dMinI, dataInDouble[srcElem] ); // min over intersecting source faces
814 dMaxI = fmax( dMaxI, dataInDouble[srcElem] ); // max over intersecting source faces
815 }
816
817 // Update the min and max values of the target data
818 vecLocalLowerBound[i] = dMinI;
819 vecLocalUpperBound[i] = dMaxI;
820 }
821
822 for( size_t i = 0; i < nTargetCount; i++ )
823 {
824 dataLowerBound[i] = vecLocalLowerBound[i] - dataOutDouble[i];
825 dataUpperBound[i] = vecLocalUpperBound[i] - dataOutDouble[i];
826 }
827 } // caasType == CAAS_LOCAL
828
829 // Invoke CAAS or QLT application on the map
830 if( fabs( dMassDiff ) > 1e-20 )
831 {
832 if( caasType == CAAS_QLT )
833 dMassDiff = QLTLimiter( caasIteration, dataOutDouble, dataLowerBound, dataUpperBound, massVector );
834 else
835 CAASLimiter( dataOutDouble, dataLowerBound, dataUpperBound, dMassDiff );
836 }
837
838 // Announce output mass
839 double dMassDiffPost = 0.0;
840 for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
841 {
842 const int ixS = m_meshOverlap->vecSourceFaceIx[i];
843 const int ixT = m_meshOverlap->vecTargetFaceIx[i];
844
845 if( ixT < 0 ) continue; // skip ghost target faces
846
847 // Update the mass difference between source and target faces
848 // linked to the overlap mesh element
849 dMassDiffPost += ( dataInDouble[ixS] * m_dOverlapAreas[i] ) - // source mass
850 ( dataOutDouble[ixT] * m_dOverlapAreas[i] ); // target mass
851 }
852 // massDefect.second = fabs( dMassDiffPost / ( dSourceMax - dSourceMin ) );
853 massDefect.second = dMassDiffPost;
854 }
855
856 // Ideally should perform an AllReduce here to get the global mass difference across all processors
857 // But if we satisfy the constraint on every task, essentially, the global mass difference should be zero!
858 return massDefect;
859 }
860
861 ///////////////////////////////////////////////////////////////////////////////
862
863 // ** Kahan Summation Algorithm for improved numerical accuracy **
864 struct KahanSum
865 {
866 double sum = 0.0;
867 double correction = 0.0;
868
869 void add( double value )
870 {
871 double y = value - correction; // Correct the input
872 double t = sum + y; // Perform the sum
873 correction = ( t - sum ) - y; // Update correction
874 sum = t; // Store the new sum
875 }
876
877 double result() const
878 {
879 return sum;
880 }
881 };
882
883 // Pairwise summation helper function
884 inline double pairwiseSum( const std::set< double >& sorted )
885 {
886 if( sorted.empty() ) return 0.0;
887 if( sorted.size() == 1 ) return *sorted.begin();
888
889 // Accumulate pairwise to minimize rounding error
890 double sum = 0.0;
891 for( double val : sorted )
892 sum += val;
893 return sum;
894 }
895
896 // Pairwise summation helper function
897 inline double pairwiseKahanSum( const std::set< double >& sorted )
898 {
899 if( sorted.empty() ) return 0.0;
900 if( sorted.size() == 1 ) return *sorted.begin();
901
902 // Accumulate pairwise to minimize rounding error
903 // Apply pairwise summation with Kahan correction
904 KahanSum kahan;
905 for( double val : sorted )
906 kahan.add( val );
907 return kahan.result();
908 }
909
910 // Sparse matrix-vector multiplication using pairwise summation
911 inline void deterministicSparseMatVecMul( const typename moab::TempestOnlineMap::WeightMatrix& A,
912 const typename moab::TempestOnlineMap::WeightColVector& x,
913 typename moab::TempestOnlineMap::WeightRowVector& result )
914 {
915 constexpr bool useKahanSum = false;
916 constexpr bool usePairwiseSum = false;
917
918 result.setZero(); // Ensure no uninitialized memory issues
919
920 // Iterate row-wise to enforce a fixed summation order
921 for( int row = 0; row < A.outerSize(); ++row )
922 {
923 std::set< double > accumulators;
924 for( typename moab::TempestOnlineMap::WeightMatrix::InnerIterator it( A, row ); it; ++it )
925 {
926 // accumulators contains the sorted values of the product: A(row, col) * x(col)
927 accumulators.insert( it.value() * x( it.col() ) );
928 }
929 if( usePairwiseSum ) result( row ) = pairwiseSum( accumulators );
930 if( useKahanSum ) result( row ) = pairwiseKahanSum( accumulators );
931
932 if( !usePairwiseSum && !useKahanSum )
933 {
934 double sum = 0.0;
935 for( double val : accumulators )
936 sum += val;
937 result( row ) = sum;
938 }
939 }
940 }
941
942 //
943 // Perform a deterministic sparse matrix-vector multiplication
944 inline void deterministicSparseMatVecMulKahan( const typename moab::TempestOnlineMap::WeightMatrix& A,
945 const typename moab::TempestOnlineMap::WeightColVector& x,
946 typename moab::TempestOnlineMap::WeightRowVector& result )
947 {
948 result.setZero(); // Ensure no uninitialized memory issues
949
950 // Iterate row-wise to enforce a fixed summation order
951 for( int row = 0; row < A.outerSize(); ++row )
952 {
953 KahanSum kahan;
954 for( typename moab::TempestOnlineMap::WeightMatrix::InnerIterator it( A, row ); it; ++it )
955 {
956 double product = it.value() * x( it.col() ); // Compute product
957 kahan.add( product );
958 }
959
960 result( row ) = kahan.result();
961 }
962 }
963
964 // Perform a deterministic sparse matrix-vector multiplication
965 inline void deterministicSparseMatVecMulClean( const typename moab::TempestOnlineMap::WeightMatrix& A,
966 const typename moab::TempestOnlineMap::WeightColVector& x,
967 typename moab::TempestOnlineMap::WeightRowVector& result )
968 {
969 result.setZero(); // Ensure no uninitialized memory issues
970
971 // Iterate row-wise to enforce a fixed summation order
972 for( int row = 0; row < A.outerSize(); ++row )
973 {
974 for( typename moab::TempestOnlineMap::WeightMatrix::InnerIterator it( A, row ); it; ++it )
975 {
976 const double product = it.value() * x( it.col() ); // Compute product
977 result( it.row() ) += product;
978 }
979 }
980 }
981
982 inline void deterministicSparseMatVecMulNative( const typename moab::TempestOnlineMap::WeightMatrix& A,
983 const typename moab::TempestOnlineMap::WeightColVector& x,
984 typename moab::TempestOnlineMap::WeightRowVector& result )
985 {
986 result = A * x; // Perform the matrix-vector multiplication using Eigen3
987 }
988
989 // Deterministic sparse matrix-vector multiplication with A^T * x using pairwise summation
990 inline void deterministicSparseMatTransposeVecMul( const typename moab::TempestOnlineMap::WeightMatrix& A,
991 const typename moab::TempestOnlineMap::WeightRowVector& x,
992 typename moab::TempestOnlineMap::WeightColVector& result )
993 {
994 result.setZero(); // Ensure no uninitialized memory issues
995
996 // Temporary storage for pairwise summation
997 std::vector< std::set< double > > accumulators( A.cols() );
998
999 // Iterate over A row-wise, but accumulate into result as if computing A^T * x
1000 for( int row = 0; row < A.outerSize(); ++row )
1001 {
1002 for( typename moab::TempestOnlineMap::WeightMatrix::InnerIterator it( A, row ); it; ++it )
1003 {
1004 accumulators[it.col()].insert( it.value() * x( row ) );
1005 }
1006 }
1007
1008 // Compute final sum using pairwise summation for each entry
1009 for( int col = 0; col < A.cols(); ++col )
1010 {
1011 // result( col ) = pairwiseSum( accumulators[col] );
1012 result( col ) = pairwiseKahanSum( accumulators[col] );
1013 }
1014 }
1015
1016 // Perform a deterministic sparse matrix-vector multiplication
1017 inline void deterministicSparseMatTransposeVecMulClean( const typename moab::TempestOnlineMap::WeightMatrix& A,
1018 const typename moab::TempestOnlineMap::WeightRowVector& x,
1019 typename moab::TempestOnlineMap::WeightColVector& result )
1020 {
1021 result.setZero(); // Ensure no uninitialized memory issues
1022
1023 // Iterate over A row-wise, but accumulate into result as if computing A^T * x
1024 for( int row = 0; row < A.outerSize(); ++row )
1025 {
1026 for( typename moab::TempestOnlineMap::WeightMatrix::InnerIterator it( A, row ); it; ++it )
1027 {
1028 const double product = it.value() * x( row ); // Compute product
1029 result( it.col() ) += product; // Accumulate contributions to the corresponding row in A^T
1030 }
1031 }
1032 }
1033
1034 // Perform a deterministic sparse matrix-vector multiplication
1035 inline void deterministicSparseMatTransposeVecMulNative( const typename moab::TempestOnlineMap::WeightMatrix& A,
1036 const typename moab::TempestOnlineMap::WeightRowVector& x,
1037 typename moab::TempestOnlineMap::WeightColVector& result )
1038 {
1039 result = A.adjoint() * x; // Perform the adjoint.matrix-vector multiplication using Eigen3
1040 }
1041 ///////////////////////////////////////////////////////////////////////////////
1042
1043 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( std::vector< double >& srcVals,
1044 std::vector< double >& tgtVals,
1045 bool transpose )
1046 {
1047 // Reset the source and target data first
1048 m_rowVector.setZero();
1049 m_colVector.setZero();
1050 #ifdef VERBOSE
1051 std::stringstream sstr;
1052 static int callId = 0;
1053 callId++;
1054 sstr << "projection_id_" << callId << "_s_" << size << "_rk_" << rank << ".txt";
1055 std::ofstream output_file( sstr.str() );
1056 #endif
1057 // Perform the actual projection of weights: application of weight matrix onto the source
1058 // solution vector
1059 if( transpose )
1060 {
1061 // Permute the source data first
1062 for( unsigned i = 0; i < srcVals.size(); ++i )
1063 {
1064 if( row_dtoc_dofmap[i] >= 0 && row_dtoc_dofmap[i] < m_rowVector.size() )
1065 m_rowVector( row_dtoc_dofmap[i] ) = srcVals[i]; // permute and set the row (source) vector properly
1066 }
1067
1068 // deterministicSparseMatTransposeVecMulClean( m_weightMatrix, m_rowVector, m_colVector );
1069 deterministicSparseMatTransposeVecMul( m_weightMatrix, m_rowVector, m_colVector );
1070 // deterministicSparseMatTransposeVecMulNative( m_weightMatrix, m_rowVector, m_colVector );
1071
1072 // Permute the resulting target data back
1073 for( unsigned i = 0; i < tgtVals.size(); ++i )
1074 {
1075 if( col_dtoc_dofmap[i] >= 0 && col_dtoc_dofmap[i] < m_colVector.size() )
1076 tgtVals[i] = m_colVector( col_dtoc_dofmap[i] ); // permute and set the row (source) vector properly
1077 }
1078 }
1079 else
1080 {
1081 // Permute the source data first
1082 #ifdef VERBOSE
1083 output_file << "ColVector: " << m_colVector.size() << ", SrcVals: " << srcVals.size()
1084 << ", Sizes: " << m_nTotDofs_SrcCov << ", " << col_dtoc_dofmap.size() << "\n";
1085 #endif
1086 for( unsigned i = 0; i < srcVals.size(); ++i )
1087 {
1088 if( col_dtoc_dofmap[i] >= 0 && col_dtoc_dofmap[i] < m_colVector.size() )
1089 {
1090 m_colVector( col_dtoc_dofmap[i] ) = srcVals[i]; // permute and set the row (source) vector properly
1091 #ifdef VERBOSE
1092 output_file << i << " " << col_gdofmap[col_dtoc_dofmap[i]] + 1 << " " << srcVals[i] << "\n";
1093 #endif
1094 }
1095 }
1096
1097 // deterministicSparseMatVecMulClean( m_weightMatrix, m_colVector, m_rowVector );
1098 deterministicSparseMatVecMul( m_weightMatrix, m_colVector, m_rowVector );
1099 // deterministicSparseMatVecMulNative( m_weightMatrix, m_colVector, m_rowVector );
1100 // deterministicSparseMatVecMulKahan( m_weightMatrix, m_colVector, m_rowVector );
1101
1102 // Permute the resulting target data back
1103 #ifdef VERBOSE
1104 output_file
1105 << "RowVector: " << m_rowVector.size() << ", TgtVals:" << tgtVals.size() << ", Sizes: " << m_nTotDofs_Dest
1106 << ", " << row_gdofmap.size() << "\n";
1107 #endif
1108 for( unsigned i = 0; i < tgtVals.size(); ++i )
1109 {
1110 if( row_dtoc_dofmap[i] >= 0 && row_dtoc_dofmap[i] < m_rowVector.size() )
1111 {
1112 tgtVals[i] = m_rowVector( row_dtoc_dofmap[i] ); // permute and set the row (source) vector properly
1113 #ifdef VERBOSE
1114 output_file << i << " " << row_gdofmap[row_dtoc_dofmap[i]] + 1 << " " << tgtVals[i] << "\n";
1115 #endif
1116 }
1117 }
1118 }
1119
1120 // if( caasType != CAAS_NONE )
1121 // {
1122 // constexpr int nmax_caas_iterations = 5;
1123 // double mismatch = 1.0;
1124 // int caasIteration = 0;
1125 // while( mismatch > 1e-15 &&
1126 // caasIteration++ < nmax_caas_iterations ) // iterate until convergence or a maximum of 5 iterations
1127 // {
1128 // std::pair< double, double > mDefect = this->ApplyCAASLimiting( srcVals, tgtVals, caasType );
1129 // if( m_remapper->verbose )
1130 // printf( "Rank %d: -- Iteration: %d, Net original mass defect: %3.4e, mass defect post-CAAS: %3.4e\n",
1131 // m_remapper->rank, caasIteration, mDefect.first, mDefect.second );
1132 // mismatch = mDefect.second;
1133 // }
1134 // }
1135
1136 #ifdef VERBOSE
1137 output_file.flush(); // required here
1138 output_file.close();
1139 #endif
1140
1141 // All done with matvec application
1142 return moab::MB_SUCCESS;
1143 }
1144
1145 #endif
1146
1147 ///////////////////////////////////////////////////////////////////////////////
1148
1149 extern void ForceConsistencyConservation3( const DataArray1D< double >& vecSourceArea,
1150 const DataArray1D< double >& vecTargetArea,
1151 DataArray2D< double >& dCoeff,
1152 bool fMonotone,
1153 bool fSparseConstraints = false );
1154
1155 ///////////////////////////////////////////////////////////////////////////////
1156
1157 extern void ForceIntArrayConsistencyConservation( const DataArray1D< double >& vecSourceArea,
1158 const DataArray1D< double >& vecTargetArea,
1159 DataArray2D< double >& dCoeff,
1160 bool fMonotone );
1161
1162 ///////////////////////////////////////////////////////////////////////////////
1163
1164 void moab::TempestOnlineMap::LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
1165 const DataArray3D< double >& dataGLLJacobian,
1166 int nMonotoneType,
1167 bool fContinuousIn,
1168 bool fNoConservation )
1169 {
1170 // Order of the polynomial interpolant
1171 int nP = dataGLLNodes.GetRows();
1172
1173 // Order of triangular quadrature rule
1174 const int TriQuadRuleOrder = 4;
1175
1176 // Triangular quadrature rule
1177 TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
1178
1179 int TriQuadraturePoints = triquadrule.GetPoints();
1180
1181 const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
1182
1183 const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
1184
1185 // Sample coefficients
1186 DataArray2D< double > dSampleCoeff( nP, nP );
1187
1188 // GLL Quadrature nodes on quadrilateral elements
1189 DataArray1D< double > dG;
1190 DataArray1D< double > dW;
1191 GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
1192
1193 // Announcements
1194 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1195 dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
1196 if( is_root )
1197 {
1198 dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
1199 dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
1200 dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
1201 }
1202
1203 // Get SparseMatrix represntation of the OfflineMap
1204 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1205
1206 // NodeVector from m_meshOverlap
1207 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1208 const NodeVector& nodesFirst = m_meshInputCov->nodes;
1209
1210 // Vector of source areas
1211 DataArray1D< double > vecSourceArea( nP * nP );
1212
1213 DataArray1D< double > vecTargetArea;
1214 DataArray2D< double > dCoeff;
1215
1216 #ifdef VERBOSE
1217 std::stringstream sstr;
1218 sstr << "remapdata_" << rank << ".txt";
1219 std::ofstream output_file( sstr.str() );
1220 #endif
1221
1222 // Current Overlap Face
1223 int ixOverlap = 0;
1224 #ifdef VERBOSE
1225 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1226 #endif
1227 // generic triangle used for area computation, for triangles around the center of overlap face;
1228 // used for overlap faces with more than 4 edges;
1229 // nodes array will be set for each triangle;
1230 // these triangles are not part of the mesh structure, they are just temporary during
1231 // aforementioned decomposition.
1232 Face faceTri( 3 );
1233 NodeVector nodes( 3 );
1234 faceTri.SetNode( 0, 0 );
1235 faceTri.SetNode( 1, 1 );
1236 faceTri.SetNode( 2, 2 );
1237
1238 // Loop over all input Faces
1239 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1240 {
1241 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1242
1243 if( faceFirst.edges.size() != 4 )
1244 {
1245 _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
1246 }
1247 #ifdef VERBOSE
1248 // Announce computation progress
1249 if( ixFirst % outputFrequency == 0 && is_root )
1250 {
1251 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1252 }
1253 #endif
1254 // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
1255 // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
1256 // However, the relation with MOAB and Tempest will go out of the roof
1257
1258 // Determine how many overlap Faces and triangles are present
1259 int nOverlapFaces = 0;
1260 size_t ixOverlapTemp = ixOverlap;
1261 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1262 {
1263 // if( m_meshOverlap->vecTargetFaceIx[ixOverlapTemp] < 0 ) continue; // skip ghost target faces
1264 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1265 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1266
1267 nOverlapFaces++;
1268 }
1269
1270 // No overlaps
1271 if( nOverlapFaces == 0 )
1272 {
1273 continue;
1274 }
1275
1276 // Allocate remap coefficients array for meshFirst Face
1277 DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
1278
1279 // Find the local remap coefficients
1280 for( int j = 0; j < nOverlapFaces; j++ )
1281 {
1282 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
1283 if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 ) // machine precision
1284 {
1285 Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
1286 m_meshOverlap->vecFaceArea[ixOverlap + j] );
1287 int n = faceOverlap.edges.size();
1288 Announce( "Number nodes: %d", n );
1289 for( int k = 0; k < n; k++ )
1290 {
1291 Node nd = nodesOverlap[faceOverlap[k]];
1292 Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
1293 }
1294 continue;
1295 }
1296
1297 // #ifdef VERBOSE
1298 // if ( is_root )
1299 // Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
1300 // m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
1301 // m_meshOverlap->vecFaceArea[ixOverlap + j] );
1302 // #endif
1303
1304 int nbEdges = faceOverlap.edges.size();
1305 int nOverlapTriangles = 1;
1306 Node center; // not used if nbEdges == 3
1307 if( nbEdges > 3 )
1308 { // decompose from center in this case
1309 nOverlapTriangles = nbEdges;
1310 for( int k = 0; k < nbEdges; k++ )
1311 {
1312 const Node& node = nodesOverlap[faceOverlap[k]];
1313 center = center + node;
1314 }
1315 center = center / nbEdges;
1316 center = center.Normalized(); // project back on sphere of radius 1
1317 }
1318
1319 Node node0, node1, node2;
1320 double dTriangleArea;
1321
1322 // Loop over all sub-triangles of this Overlap Face
1323 for( int k = 0; k < nOverlapTriangles; k++ )
1324 {
1325 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1326 {
1327 node0 = nodesOverlap[faceOverlap[0]];
1328 node1 = nodesOverlap[faceOverlap[1]];
1329 node2 = nodesOverlap[faceOverlap[2]];
1330 dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1331 }
1332 else // decompose polygon in triangles around the center
1333 {
1334 node0 = center;
1335 node1 = nodesOverlap[faceOverlap[k]];
1336 int k1 = ( k + 1 ) % nbEdges;
1337 node2 = nodesOverlap[faceOverlap[k1]];
1338 nodes[0] = center;
1339 nodes[1] = node1;
1340 nodes[2] = node2;
1341 dTriangleArea = CalculateFaceArea( faceTri, nodes );
1342 }
1343 // Coordinates of quadrature Node
1344 for( int l = 0; l < TriQuadraturePoints; l++ )
1345 {
1346 Node nodeQuadrature;
1347 nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
1348 TriQuadratureG[l][2] * node2.x;
1349
1350 nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
1351 TriQuadratureG[l][2] * node2.y;
1352
1353 nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
1354 TriQuadratureG[l][2] * node2.z;
1355
1356 nodeQuadrature = nodeQuadrature.Normalized();
1357
1358 // Find components of quadrature point in basis
1359 // of the first Face
1360 double dAlpha;
1361 double dBeta;
1362
1363 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
1364
1365 // Check inverse map value
1366 if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
1367 ( dBeta > 1.0 + 1.0e-13 ) )
1368 {
1369 _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
1370 "(%1.5e %1.5e)",
1371 j, l, dAlpha, dBeta );
1372 }
1373
1374 // Sample the finite element at this point
1375 SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
1376
1377 // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
1378 for( int p = 0; p < nP; p++ )
1379 {
1380 for( int q = 0; q < nP; q++ )
1381 {
1382 dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
1383 m_meshOverlap->vecFaceArea[ixOverlap + j];
1384 }
1385 }
1386 }
1387 }
1388 }
1389
1390 #ifdef VERBOSE
1391 output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1392 for( int j = 0; j < nOverlapFaces; j++ )
1393 {
1394 for( int p = 0; p < nP; p++ )
1395 {
1396 for( int q = 0; q < nP; q++ )
1397 {
1398 output_file << dRemapCoeff[p][q][j] << " ";
1399 }
1400 }
1401 }
1402 output_file << std::endl;
1403 #endif
1404
1405 // Force consistency and conservation
1406 if( !fNoConservation )
1407 {
1408 double dTargetArea = 0.0;
1409 for( int j = 0; j < nOverlapFaces; j++ )
1410 {
1411 dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
1412 }
1413
1414 for( int p = 0; p < nP; p++ )
1415 {
1416 for( int q = 0; q < nP; q++ )
1417 {
1418 vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1419 }
1420 }
1421
1422 const double areaTolerance = 1e-10;
1423 // Source elements are completely covered by target volumes
1424 if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
1425 {
1426 vecTargetArea.Allocate( nOverlapFaces );
1427 for( int j = 0; j < nOverlapFaces; j++ )
1428 {
1429 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1430 }
1431
1432 dCoeff.Allocate( nOverlapFaces, nP * nP );
1433
1434 for( int j = 0; j < nOverlapFaces; j++ )
1435 {
1436 for( int p = 0; p < nP; p++ )
1437 {
1438 for( int q = 0; q < nP; q++ )
1439 {
1440 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1441 }
1442 }
1443 }
1444
1445 // Target volumes only partially cover source elements
1446 }
1447 else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
1448 {
1449 double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
1450
1451 vecTargetArea.Allocate( nOverlapFaces + 1 );
1452 for( int j = 0; j < nOverlapFaces; j++ )
1453 {
1454 vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1455 }
1456 vecTargetArea[nOverlapFaces] = dExtraneousArea;
1457
1458 #ifdef VERBOSE
1459 Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
1460 m_meshInputCov->vecFaceArea[ixFirst] );
1461 #endif
1462 if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
1463 {
1464 _EXCEPTIONT( "Partial element area exceeds total element area" );
1465 }
1466
1467 dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
1468
1469 for( int j = 0; j < nOverlapFaces; j++ )
1470 {
1471 for( int p = 0; p < nP; p++ )
1472 {
1473 for( int q = 0; q < nP; q++ )
1474 {
1475 dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1476 }
1477 }
1478 }
1479 for( int p = 0; p < nP; p++ )
1480 {
1481 for( int q = 0; q < nP; q++ )
1482 {
1483 dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1484 }
1485 }
1486 for( int j = 0; j < nOverlapFaces; j++ )
1487 {
1488 for( int p = 0; p < nP; p++ )
1489 {
1490 for( int q = 0; q < nP; q++ )
1491 {
1492 dCoeff[nOverlapFaces][p * nP + q] -=
1493 dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
1494 }
1495 }
1496 }
1497 for( int p = 0; p < nP; p++ )
1498 {
1499 for( int q = 0; q < nP; q++ )
1500 {
1501 dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
1502 }
1503 }
1504
1505 // Source elements only partially cover target volumes
1506 }
1507 else
1508 {
1509 Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
1510 m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
1511 _EXCEPTIONT( "Target grid must be a subset of source grid" );
1512 }
1513
1514 ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
1515 /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
1516
1517 for( int j = 0; j < nOverlapFaces; j++ )
1518 {
1519 for( int p = 0; p < nP; p++ )
1520 {
1521 for( int q = 0; q < nP; q++ )
1522 {
1523 dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
1524 }
1525 }
1526 }
1527 }
1528
1529 #ifdef VERBOSE
1530 // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
1531 // for ( int j = 0; j < nOverlapFaces; j++ )
1532 // {
1533 // for ( int p = 0; p < nP; p++ )
1534 // {
1535 // for ( int q = 0; q < nP; q++ )
1536 // {
1537 // output_file << dRemapCoeff[p][q][j] << " ";
1538 // }
1539 // }
1540 // }
1541 // output_file << std::endl;
1542 #endif
1543
1544 // Put these remap coefficients into the SparseMatrix map
1545 for( int j = 0; j < nOverlapFaces; j++ )
1546 {
1547 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1548
1549 // signal to not participate, because it is a ghost target
1550 if( ixSecondFace < 0 ) continue; // do not do anything
1551
1552 for( int p = 0; p < nP; p++ )
1553 {
1554 for( int q = 0; q < nP; q++ )
1555 {
1556 if( fContinuousIn )
1557 {
1558 int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
1559
1560 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1561 m_meshOverlap->vecFaceArea[ixOverlap + j] /
1562 m_meshOutput->vecFaceArea[ixSecondFace];
1563 }
1564 else
1565 {
1566 int ixFirstNode = ixFirst * nP * nP + p * nP + q;
1567
1568 smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1569 m_meshOverlap->vecFaceArea[ixOverlap + j] /
1570 m_meshOutput->vecFaceArea[ixSecondFace];
1571 }
1572 }
1573 }
1574 }
1575 // Increment the current overlap index
1576 ixOverlap += nOverlapFaces;
1577 }
1578 #ifdef VERBOSE
1579 output_file.flush(); // required here
1580 output_file.close();
1581 #endif
1582
1583 return;
1584 }
1585
1586 ///////////////////////////////////////////////////////////////////////////////
1587
1588 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
1589 const DataArray3D< double >& dataGLLJacobianIn,
1590 const DataArray3D< int >& dataGLLNodesOut,
1591 const DataArray3D< double >& dataGLLJacobianOut,
1592 const DataArray1D< double >& dataNodalAreaOut,
1593 int nPin,
1594 int nPout,
1595 int nMonotoneType,
1596 bool fContinuousIn,
1597 bool fContinuousOut,
1598 bool fNoConservation )
1599 {
1600 // Triangular quadrature rule
1601 TriangularQuadratureRule triquadrule( 8 );
1602
1603 const DataArray2D< double >& dG = triquadrule.GetG();
1604 const DataArray1D< double >& dW = triquadrule.GetW();
1605
1606 // Get SparseMatrix represntation of the OfflineMap
1607 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1608
1609 // Sample coefficients
1610 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1611 DataArray2D< double > dSampleCoeffOut( nPout, nPout );
1612
1613 // Announcemnets
1614 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1615 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
1616 if( is_root )
1617 {
1618 dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
1619 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1620 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1621 }
1622
1623 // Build the integration array for each element on m_meshOverlap
1624 DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
1625
1626 // Number of overlap Faces per source Face
1627 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1628
1629 int ixOverlap = 0;
1630 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1631 {
1632 // Determine how many overlap Faces and triangles are present
1633 int nOverlapFaces = 0;
1634 size_t ixOverlapTemp = ixOverlap;
1635 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1636 {
1637 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1638 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
1639 {
1640 break;
1641 }
1642
1643 nOverlapFaces++;
1644 }
1645
1646 nAllOverlapFaces[ixFirst] = nOverlapFaces;
1647
1648 // Increment the current overlap index
1649 ixOverlap += nAllOverlapFaces[ixFirst];
1650 }
1651
1652 // Geometric area of each output node
1653 DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
1654
1655 // Area of each overlap element in the output basis
1656 DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
1657
1658 // Loop through all faces on m_meshInputCov
1659 ixOverlap = 0;
1660 #ifdef VERBOSE
1661 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1662 #endif
1663 if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
1664
1665 // generic triangle used for area computation, for triangles around the center of overlap face;
1666 // used for overlap faces with more than 4 edges;
1667 // nodes array will be set for each triangle;
1668 // these triangles are not part of the mesh structure, they are just temporary during
1669 // aforementioned decomposition.
1670 Face faceTri( 3 );
1671 NodeVector nodes( 3 );
1672 faceTri.SetNode( 0, 0 );
1673 faceTri.SetNode( 1, 1 );
1674 faceTri.SetNode( 2, 2 );
1675
1676 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1677 {
1678 #ifdef VERBOSE
1679 // Announce computation progress
1680 if( ixFirst % outputFrequency == 0 && is_root )
1681 {
1682 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1683 }
1684 #endif
1685 // Quantities from the First Mesh
1686 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1687
1688 const NodeVector& nodesFirst = m_meshInputCov->nodes;
1689
1690 // Number of overlapping Faces and triangles
1691 int nOverlapFaces = nAllOverlapFaces[ixFirst];
1692
1693 if( !nOverlapFaces ) continue;
1694
1695 // // Calculate total element Jacobian
1696 // double dTotalJacobian = 0.0;
1697 // for (int s = 0; s < nPin; s++) {
1698 // for (int t = 0; t < nPin; t++) {
1699 // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
1700 // }
1701 // }
1702
1703 // Loop through all Overlap Faces
1704 for( int i = 0; i < nOverlapFaces; i++ )
1705 {
1706 // Quantities from the overlap Mesh
1707 const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
1708
1709 const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1710
1711 // Quantities from the Second Mesh
1712 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1713
1714 // signal to not participate, because it is a ghost target
1715 if( ixSecond < 0 ) continue; // do not do anything
1716
1717 const NodeVector& nodesSecond = m_meshOutput->nodes;
1718
1719 const Face& faceSecond = m_meshOutput->faces[ixSecond];
1720
1721 int nbEdges = faceOverlap.edges.size();
1722 int nOverlapTriangles = 1;
1723 Node center; // not used if nbEdges == 3
1724 if( nbEdges > 3 )
1725 { // decompose from center in this case
1726 nOverlapTriangles = nbEdges;
1727 for( int k = 0; k < nbEdges; k++ )
1728 {
1729 const Node& node = nodesOverlap[faceOverlap[k]];
1730 center = center + node;
1731 }
1732 center = center / nbEdges;
1733 center = center.Normalized(); // project back on sphere of radius 1
1734 }
1735
1736 Node node0, node1, node2;
1737 double dTriArea;
1738
1739 // Loop over all sub-triangles of this Overlap Face
1740 for( int j = 0; j < nOverlapTriangles; j++ )
1741 {
1742 if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1743 {
1744 node0 = nodesOverlap[faceOverlap[0]];
1745 node1 = nodesOverlap[faceOverlap[1]];
1746 node2 = nodesOverlap[faceOverlap[2]];
1747 dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1748 }
1749 else // decompose polygon in triangles around the center
1750 {
1751 node0 = center;
1752 node1 = nodesOverlap[faceOverlap[j]];
1753 int j1 = ( j + 1 ) % nbEdges;
1754 node2 = nodesOverlap[faceOverlap[j1]];
1755 nodes[0] = center;
1756 nodes[1] = node1;
1757 nodes[2] = node2;
1758 dTriArea = CalculateFaceArea( faceTri, nodes );
1759 }
1760
1761 for( int k = 0; k < triquadrule.GetPoints(); k++ )
1762 {
1763 // Get the nodal location of this point
1764 double dX[3];
1765
1766 dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
1767 dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
1768 dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
1769
1770 double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
1771
1772 dX[0] /= dMag;
1773 dX[1] /= dMag;
1774 dX[2] /= dMag;
1775
1776 Node nodeQuadrature( dX[0], dX[1], dX[2] );
1777
1778 // Find the components of this quadrature point in the basis
1779 // of the first Face.
1780 double dAlphaIn;
1781 double dBetaIn;
1782
1783 ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
1784
1785 // Find the components of this quadrature point in the basis
1786 // of the second Face.
1787 double dAlphaOut;
1788 double dBetaOut;
1789
1790 ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
1791
1792 /*
1793 // Check inverse map value
1794 if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
1795 (dBetaIn < 0.0) || (dBetaIn > 1.0)
1796 ) {
1797 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1798 dAlphaIn, dBetaIn);
1799 }
1800
1801 // Check inverse map value
1802 if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
1803 (dBetaOut < 0.0) || (dBetaOut > 1.0)
1804 ) {
1805 _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1806 dAlphaOut, dBetaOut);
1807 }
1808 */
1809 // Sample the First finite element at this point
1810 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1811
1812 // Sample the Second finite element at this point
1813 SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
1814
1815 // Overlap output area
1816 for( int s = 0; s < nPout; s++ )
1817 {
1818 for( int t = 0; t < nPout; t++ )
1819 {
1820 double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
1821
1822 dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
1823
1824 dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
1825 }
1826 }
1827
1828 // Compute overlap integral
1829 int ixp = 0;
1830 for( int p = 0; p < nPin; p++ )
1831 {
1832 for( int q = 0; q < nPin; q++ )
1833 {
1834 int ixs = 0;
1835 for( int s = 0; s < nPout; s++ )
1836 {
1837 for( int t = 0; t < nPout; t++ )
1838 {
1839 // Sample the Second finite element at this point
1840 dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
1841 dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
1842
1843 ixs++;
1844 }
1845 }
1846
1847 ixp++;
1848 }
1849 }
1850 }
1851 }
1852 }
1853
1854 // Coefficients
1855 DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
1856
1857 for( int i = 0; i < nOverlapFaces; i++ )
1858 {
1859 // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1860
1861 int ixp = 0;
1862 for( int p = 0; p < nPin; p++ )
1863 {
1864 for( int q = 0; q < nPin; q++ )
1865 {
1866 int ixs = 0;
1867 for( int s = 0; s < nPout; s++ )
1868 {
1869 for( int t = 0; t < nPout; t++ )
1870 {
1871 dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
1872 dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1873
1874 ixs++;
1875 }
1876 }
1877
1878 ixp++;
1879 }
1880 }
1881 }
1882
1883 // Source areas
1884 DataArray1D< double > vecSourceArea( nPin * nPin );
1885
1886 for( int p = 0; p < nPin; p++ )
1887 {
1888 for( int q = 0; q < nPin; q++ )
1889 {
1890 vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
1891 }
1892 }
1893
1894 // Target areas
1895 DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
1896
1897 for( int i = 0; i < nOverlapFaces; i++ )
1898 {
1899 // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1900 int ixs = 0;
1901 for( int s = 0; s < nPout; s++ )
1902 {
1903 for( int t = 0; t < nPout; t++ )
1904 {
1905 vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
1906
1907 ixs++;
1908 }
1909 }
1910 }
1911
1912 // Force consistency and conservation
1913 if( !fNoConservation )
1914 {
1915 ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
1916 }
1917
1918 // Update global coefficients
1919 for( int i = 0; i < nOverlapFaces; i++ )
1920 {
1921 int ixp = 0;
1922 for( int p = 0; p < nPin; p++ )
1923 {
1924 for( int q = 0; q < nPin; q++ )
1925 {
1926 int ixs = 0;
1927 for( int s = 0; s < nPout; s++ )
1928 {
1929 for( int t = 0; t < nPout; t++ )
1930 {
1931 dGlobalIntArray[ixp][ixOverlap + i][ixs] =
1932 dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1933
1934 ixs++;
1935 }
1936 }
1937
1938 ixp++;
1939 }
1940 }
1941 }
1942
1943 #ifdef VVERBOSE
1944 // Check column sums (conservation)
1945 for( int i = 0; i < nPin * nPin; i++ )
1946 {
1947 double dColSum = 0.0;
1948 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1949 {
1950 dColSum += dCoeff[j][i] * vecTargetArea[j];
1951 }
1952 printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
1953 }
1954
1955 // Check row sums (consistency)
1956 for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1957 {
1958 double dRowSum = 0.0;
1959 for( int i = 0; i < nPin * nPin; i++ )
1960 {
1961 dRowSum += dCoeff[j][i];
1962 }
1963 printf( "Row %i: %1.15e\n", j, dRowSum );
1964 }
1965 #endif
1966
1967 // Increment the current overlap index
1968 ixOverlap += nOverlapFaces;
1969 }
1970
1971 // Build redistribution map within target element
1972 if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
1973 DataArray1D< double > dRedistSourceArea( nPout * nPout );
1974 DataArray1D< double > dRedistTargetArea( nPout * nPout );
1975 std::vector< DataArray2D< double > > dRedistributionMaps;
1976 dRedistributionMaps.resize( m_meshOutput->faces.size() );
1977
1978 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1979 {
1980 dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
1981
1982 for( int i = 0; i < nPout * nPout; i++ )
1983 {
1984 dRedistributionMaps[ixSecond][i][i] = 1.0;
1985 }
1986
1987 for( int s = 0; s < nPout * nPout; s++ )
1988 {
1989 dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
1990 }
1991
1992 for( int s = 0; s < nPout * nPout; s++ )
1993 {
1994 dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
1995 }
1996
1997 if( !fNoConservation )
1998 {
1999 ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
2000 ( nMonotoneType != 0 ) );
2001
2002 for( int s = 0; s < nPout * nPout; s++ )
2003 {
2004 for( int t = 0; t < nPout * nPout; t++ )
2005 {
2006 dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
2007 }
2008 }
2009 }
2010 }
2011
2012 // Construct the total geometric area
2013 DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
2014 for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
2015 {
2016 for( int s = 0; s < nPout; s++ )
2017 {
2018 for( int t = 0; t < nPout; t++ )
2019 {
2020 dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
2021 dGeometricOutputArea[ixSecond][s * nPout + t];
2022 }
2023 }
2024 }
2025
2026 // Compose the integration operator with the output map
2027 ixOverlap = 0;
2028
2029 if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
2030
2031 // Map from source DOFs to target DOFs with redistribution applied
2032 DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
2033
2034 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2035 {
2036 #ifdef VERBOSE
2037 // Announce computation progress
2038 if( ixFirst % outputFrequency == 0 && is_root )
2039 {
2040 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
2041 }
2042 #endif
2043 // Number of overlapping Faces and triangles
2044 int nOverlapFaces = nAllOverlapFaces[ixFirst];
2045
2046 if( !nOverlapFaces ) continue;
2047
2048 // Put composed array into map
2049 for( int j = 0; j < nOverlapFaces; j++ )
2050 {
2051 int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
2052
2053 // signal to not participate, because it is a ghost target
2054 if( ixSecondFace < 0 ) continue; // do not do anything
2055
2056 dRedistributedOp.Zero();
2057 for( int p = 0; p < nPin * nPin; p++ )
2058 {
2059 for( int s = 0; s < nPout * nPout; s++ )
2060 {
2061 for( int t = 0; t < nPout * nPout; t++ )
2062 {
2063 dRedistributedOp[p][s] +=
2064 dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
2065 }
2066 }
2067 }
2068
2069 int ixp = 0;
2070 for( int p = 0; p < nPin; p++ )
2071 {
2072 for( int q = 0; q < nPin; q++ )
2073 {
2074 int ixFirstNode;
2075 if( fContinuousIn )
2076 {
2077 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
2078 }
2079 else
2080 {
2081 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
2082 }
2083
2084 int ixs = 0;
2085 for( int s = 0; s < nPout; s++ )
2086 {
2087 for( int t = 0; t < nPout; t++ )
2088 {
2089 int ixSecondNode;
2090 if( fContinuousOut )
2091 {
2092 ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
2093
2094 if( !fNoConservation )
2095 {
2096 smatMap( ixSecondNode, ixFirstNode ) +=
2097 dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
2098 }
2099 else
2100 {
2101 smatMap( ixSecondNode, ixFirstNode ) +=
2102 dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
2103 }
2104 }
2105 else
2106 {
2107 ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
2108
2109 if( !fNoConservation )
2110 {
2111 smatMap( ixSecondNode, ixFirstNode ) +=
2112 dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
2113 }
2114 else
2115 {
2116 smatMap( ixSecondNode, ixFirstNode ) +=
2117 dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
2118 }
2119 }
2120
2121 ixs++;
2122 }
2123 }
2124
2125 ixp++;
2126 }
2127 }
2128 }
2129
2130 // Increment the current overlap index
2131 ixOverlap += nOverlapFaces;
2132 }
2133
2134 return;
2135 }
2136
2137 ///////////////////////////////////////////////////////////////////////////////
2138
2139 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
2140 const DataArray3D< double >& /*dataGLLJacobianIn*/,
2141 const DataArray3D< int >& dataGLLNodesOut,
2142 const DataArray3D< double >& /*dataGLLJacobianOut*/,
2143 const DataArray1D< double >& dataNodalAreaOut,
2144 int nPin,
2145 int nPout,
2146 int nMonotoneType,
2147 bool fContinuousIn,
2148 bool fContinuousOut )
2149 {
2150 // Gauss-Lobatto quadrature within Faces
2151 DataArray1D< double > dGL;
2152 DataArray1D< double > dWL;
2153
2154 GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
2155
2156 // Get SparseMatrix represntation of the OfflineMap
2157 SparseMatrix< double >& smatMap = this->GetSparseMatrix();
2158
2159 // Sample coefficients
2160 DataArray2D< double > dSampleCoeffIn( nPin, nPin );
2161
2162 // Announcemnets
2163 moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
2164 dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
2165 if( is_root )
2166 {
2167 dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
2168 dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
2169 dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
2170 }
2171
2172 // Number of overlap Faces per source Face
2173 DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
2174
2175 int ixOverlap = 0;
2176
2177 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2178 {
2179 size_t ixOverlapTemp = ixOverlap;
2180 for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
2181 {
2182 // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
2183
2184 if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
2185
2186 nAllOverlapFaces[ixFirst]++;
2187 }
2188
2189 // Increment the current overlap index
2190 ixOverlap += nAllOverlapFaces[ixFirst];
2191 }
2192
2193 // Number of times this point was found
2194 DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
2195
2196 ixOverlap = 0;
2197 #ifdef VERBOSE
2198 const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
2199 #endif
2200 // Loop through all faces on m_meshInputCov
2201 for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
2202 {
2203 #ifdef VERBOSE
2204 // Announce computation progress
2205 if( ixFirst % outputFrequency == 0 && is_root )
2206 {
2207 dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
2208 }
2209 #endif
2210 // Quantities from the First Mesh
2211 const Face& faceFirst = m_meshInputCov->faces[ixFirst];
2212
2213 const NodeVector& nodesFirst = m_meshInputCov->nodes;
2214
2215 // Number of overlapping Faces and triangles
2216 int nOverlapFaces = nAllOverlapFaces[ixFirst];
2217
2218 // Loop through all Overlap Faces
2219 for( int i = 0; i < nOverlapFaces; i++ )
2220 {
2221 // Quantities from the Second Mesh
2222 int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
2223
2224 // signal to not participate, because it is a ghost target
2225 if( ixSecond < 0 ) continue; // do not do anything
2226
2227 const NodeVector& nodesSecond = m_meshOutput->nodes;
2228 const Face& faceSecond = m_meshOutput->faces[ixSecond];
2229
2230 // Loop through all nodes on the second face
2231 for( int s = 0; s < nPout; s++ )
2232 {
2233 for( int t = 0; t < nPout; t++ )
2234 {
2235 size_t ixSecondNode;
2236 if( fContinuousOut )
2237 {
2238 ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
2239 }
2240 else
2241 {
2242 ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
2243 }
2244
2245 if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
2246
2247 // Check if this node has been found already
2248 if( fSecondNodeFound[ixSecondNode] ) continue;
2249
2250 // Check this node
2251 Node node;
2252 Node dDx1G;
2253 Node dDx2G;
2254
2255 ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
2256
2257 // Find the components of this quadrature point in the basis
2258 // of the first Face.
2259 double dAlphaIn;
2260 double dBetaIn;
2261
2262 ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
2263
2264 // Check if this node is within the first Face
2265 if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
2266 ( dBetaIn > 1.0 + 1.0e-10 ) )
2267 continue;
2268
2269 // Node is within the overlap region, mark as found
2270 fSecondNodeFound[ixSecondNode] = true;
2271
2272 // Sample the First finite element at this point
2273 SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
2274
2275 // Add to map
2276 for( int p = 0; p < nPin; p++ )
2277 {
2278 for( int q = 0; q < nPin; q++ )
2279 {
2280 int ixFirstNode;
2281 if( fContinuousIn )
2282 {
2283 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
2284 }
2285 else
2286 {
2287 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
2288 }
2289
2290 smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
2291 }
2292 }
2293 }
2294 }
2295 }
2296
2297 // Increment the current overlap index
2298 ixOverlap += nOverlapFaces;
2299 }
2300
2301 // Check for missing samples
2302 for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
2303 {
2304 if( !fSecondNodeFound[i] )
2305 {
2306 _EXCEPTION1( "Can't sample point %i", i );
2307 }
2308 }
2309
2310 return;
2311 }
2312
2313 ///////////////////////////////////////////////////////////////////////////////
1.9.1.