Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
TempestLinearRemap.cpp
Go to the documentation of this file.
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 ///////////////////////////////////////////////////////////////////////////////