Mesh Oriented datABase  (version 5.6.0)
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"
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  << col_gdofmap.size() << "\n";
67  output_file << "Rows \n";
68  for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
69  output_file << iv << " " << row_gdofmap[iv] << "\n";
70  output_file << "Cols \n";
71  for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
72  output_file << iv << " " << col_gdofmap[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 
264  // signal to not participate, because it is a ghost target
265  if( ixSecondFaceLoc < 0 ) continue; // do not do anything
266 
267  m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
268  dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
269  }
270  }
271 
272  // Increment the current overlap index
273  ixOverlap += nOverlapFaces;
274  }
275 
276  return;
277 }
278 
279 ///////////////////////////////////////////////////////////////////////////////
280 
281 void moab::TempestOnlineMap::PrintMapStatistics()
282 {
283  int nrows = m_weightMatrix.rows(); // Number of rows
284  int ncols = m_weightMatrix.cols(); // Number of columns
285  int NNZ = m_weightMatrix.nonZeros(); // Number of non zero values
286 #ifdef MOAB_HAVE_MPI
287  // find out min/max for NNZ, ncols, nrows
288  // should work on std c++ 11
289  int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
290  int rarr3[6] = {0, 0, 0, 0, 0, 0};
291  MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
292 
293  int total[3] = {0, 0, 0};
294  MPI_Reduce( arr3, total, 3, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
295  if( !rank )
296  std::cout << "-> Rows (min/max/sum): (" << rarr3[1] << " / " << -rarr3[4] << " / " << total[1] << "), "
297  << " Cols (min/max/sum): (" << rarr3[2] << " / " << -rarr3[5] << " / " << total[2] << "), "
298  << " NNZ (min/max/sum): (" << rarr3[0] << " / " << -rarr3[3] << " / " << total[0] << ")\n";
299 #else
300  std::cout << "-> Rows: " << nrows << ", Cols: " << ncols << ", NNZ: " << NNZ << "\n";
301 #endif
302 }
303 
304 #ifdef MOAB_HAVE_EIGEN3
305 void moab::TempestOnlineMap::copy_tempest_sparsemat_to_eigen3()
306 {
307 #ifndef VERBOSE
308 #define VERBOSE_ACTIVATED
309 // #define VERBOSE
310 #endif
311 
312  /* Should the columns be the global size of the matrix ? */
313  m_weightMatrix.resize( m_nTotDofs_Dest, m_nTotDofs_SrcCov );
314  m_rowVector.resize( m_weightMatrix.rows() );
315  m_colVector.resize( m_weightMatrix.cols() );
316 
317 #ifdef VERBOSE
318  int locrows = std::max( m_mapRemap.GetRows(), m_nTotDofs_Dest );
319  int loccols = std::max( m_mapRemap.GetColumns(), m_nTotDofs_SrcCov );
320 
321  std::cout << m_weightMatrix.rows() << ", " << locrows << ", " << m_weightMatrix.cols() << ", " << loccols << "\n";
322  // assert(m_weightMatrix.rows() == locrows && m_weightMatrix.cols() == loccols);
323 #endif
324 
325  DataArray1D< int > lrows;
326  DataArray1D< int > lcols;
327  DataArray1D< double > lvals;
328  m_mapRemap.GetEntries( lrows, lcols, lvals );
329  size_t locvals = lvals.GetRows();
330 
331  // first matrix
332  typedef Eigen::Triplet< double > Triplet;
333  std::vector< Triplet > tripletList;
334  tripletList.reserve( locvals );
335  for( size_t iv = 0; iv < locvals; iv++ )
336  {
337  tripletList.push_back( Triplet( lrows[iv], lcols[iv], lvals[iv] ) );
338  }
339  m_weightMatrix.setFromTriplets( tripletList.begin(), tripletList.end() );
340  m_weightMatrix.makeCompressed();
341 
342 #ifdef VERBOSE
343  std::stringstream sstr;
344  sstr << "tempestmatrix.txt.0000" << rank;
345  std::ofstream output_file( sstr.str(), std::ios::out );
346  output_file << "0 " << locrows << " 0 " << loccols << "\n";
347  for( unsigned iv = 0; iv < locvals; iv++ )
348  {
349  output_file << GetRowGlobalDoF( lrows[iv] ) << " " << GetColGlobalDoF( lcols[iv] ) << " " << lvals[iv]
350  << "\n";
351  }
352  output_file.flush(); // required here
353  output_file.close();
354 #endif
355 
356 #ifdef VERBOSE_ACTIVATED
357 #undef VERBOSE_ACTIVATED
358 #undef VERBOSE
359 #endif
360  return;
361 }
362 #endif
363 
364 ///////////////////////////////////////////////////////////////////////////////
365 
366 template < typename T >
367 static std::vector< size_t > sort_indexes( const std::vector< T >& v )
368 {
369  // initialize original index locations
370  std::vector< size_t > idx( v.size() );
371  std::iota( idx.begin(), idx.end(), 0 );
372 
373  // sort indexes based on comparing values in v
374  // using std::stable_sort instead of std::sort
375  // to avoid unnecessary index re-orderings
376  // when v contains elements of equal values
377  std::stable_sort( idx.begin(), idx.end(), [&v]( size_t i1, size_t i2 ) { return fabs( v[i1] ) > fabs( v[i2] ); } );
378 
379  return idx;
380 }
381 
382 double moab::TempestOnlineMap::QLTLimiter( int caasIteration,
383  std::vector< double >& dataCorrectedField,
384  std::vector< double >& dataLowerBound,
385  std::vector< double >& dataUpperBound,
386  std::vector< double >& dMassDefect )
387 {
388  const size_t nrows = dataCorrectedField.size();
389  double dMassL = 0.0;
390  double dMassU = 0.0;
391  std::vector< double > dataCorrection( nrows );
392  double dMassDiffCum = 0.0;
393  double dLMinusU = fabs( dataUpperBound[0] - dataLowerBound[0] );
394  const DataArray1D< double >& dTargetAreas = this->m_remapper->m_target->vecFaceArea;
395 
396  // std::vector< size_t > sortedIdx = sort_indexes( dMassDefect );
397  std::vector< std::unordered_set< int > > vecAdjTargetFaces( nrows );
398  constexpr bool useMOABAdjacencies = true;
399 #ifdef USE_ComputeAdjacencyRelations
400  if( useMOABAdjacencies )
401  ComputeAdjacencyRelations( vecAdjTargetFaces, caasIteration, m_remapper->m_target_entities,
402  useMOABAdjacencies );
403  else
404  ComputeAdjacencyRelations( vecAdjTargetFaces, caasIteration, m_remapper->m_target_entities, useMOABAdjacencies,
405  this->m_remapper->m_target );
406 #else
407  moab::MeshTopoUtil mtu( m_interface );
408  ;
409 #endif
410 
411  for( size_t i = 0; i < nrows; i++ )
412  {
413  // size_t index = sortedIdx[i];
414  size_t index = i;
415  dataCorrection[index] = fmax( dataLowerBound[index], fmin( dataUpperBound[index], 0.0 ) );
416  // dMassDiff[index] = dMassDefect[index] - dTargetAreas[index] * dataCorrection[index];
417  // dMassDiff[index] = dMassDefect[index];
418 
419  dMassL += dTargetAreas[index] * dataLowerBound[index];
420  dMassU += dTargetAreas[index] * dataUpperBound[index];
421  dLMinusU = fmax( dLMinusU, fabs( dataUpperBound[index] - dataLowerBound[index] ) );
422  dMassDiffCum += dMassDefect[index] - dTargetAreas[index] * dataCorrection[index];
423 
424 #ifndef USE_ComputeAdjacencyRelations
425  vecAdjTargetFaces[index].insert( index ); // add self target face first
426  {
427  // Compute the adjacent faces to the target face
428  if( useMOABAdjacencies )
429  {
430  moab::Range ents;
431  // ents.insert( m_remapper->m_target_entities.index( m_remapper->m_target_entities[index] ) );
432  ents.insert( m_remapper->m_target_entities[index] );
433  moab::Range adjEnts;
434  moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, caasIteration );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
435  for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
436  {
437  // int adjIndex = m_interface->id_from_handle(*it)-1;
438  int adjIndex = m_remapper->m_target_entities.index( *it );
439  // printf("rank: %d, Element %lu, entity: %lu, adjIndex %d\n", rank, index, *it, adjIndex);
440  if( adjIndex >= 0 ) vecAdjTargetFaces[index].insert( adjIndex );
441  }
442  }
443  else
444  {
445  AdjacentFaceVector vecAdjFaces;
446  GetAdjacentFaceVectorByEdge( *this->m_remapper->m_target, index,
447  ( m_output_order + 1 ) * ( m_output_order + 1 ) * ( m_output_order + 1 ),
448  // ( m_output_order + 1 ) * ( m_output_order + 1 ),
449  // ( 4 ) * ( m_output_order + 1 ) * ( m_output_order + 1 ),
450  vecAdjFaces );
451 
452  // Add the adjacent faces to the target face list
453  for( auto adjFace : vecAdjFaces )
454  if( adjFace.first >= 0 )
455  vecAdjTargetFaces[index].insert( adjFace.first ); // map target face to source face
456  }
457  }
458 #endif
459  }
460 
461 #ifdef MOAB_HAVE_MPI
462  std::vector< double > localDefects( 5, 0.0 ), globalDefects( 5, 0.0 );
463  localDefects[0] = dMassL;
464  localDefects[1] = dMassU;
465  localDefects[2] = dMassDiffCum;
466  localDefects[3] = dLMinusU;
467  // localDefects[4] = dMassCorrectU;
468 
469  MPI_Allreduce( localDefects.data(), globalDefects.data(), 4, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
470 
471  dMassL = globalDefects[0];
472  dMassU = globalDefects[1];
473  dMassDiffCum = globalDefects[2];
474  dLMinusU = globalDefects[3];
475  // dMassCorrectU = globalDefects[4];
476 #endif
477 
478  //If the upper and lower bounds are too close together, just clip
479  if( fabs( dMassDiffCum ) < 1e-15 || dLMinusU < 1e-15 )
480  {
481  for( size_t i = 0; i < nrows; i++ )
482  dataCorrectedField[i] += dataCorrection[i];
483  return dMassDiffCum;
484  }
485  else
486  {
487  if( dMassL > dMassDiffCum )
488  {
489  Announce( "Lower bound mass exceeds target mass by %1.15e: CAAS will need another iteration",
490  dMassL - dMassDiffCum );
491  dMassDiffCum = dMassL;
492  // dMass -= dMassL;
493  }
494  else if( dMassU < dMassDiffCum )
495  {
496  Announce( "Target mass exceeds upper bound mass by %1.15e: CAAS will need another iteration",
497  dMassDiffCum - dMassU );
498  dMassDiffCum = dMassU;
499  // dMass -= dMassU;
500  }
501 
502  // TODO: optimize away dataMassVec by a simple transient double within the loop
503  // DataArray1D< double > dataMassVec( nrows ); //vector of mass redistribution
504  for( size_t i = 0; i < nrows; i++ )
505  {
506  // size_t index = sortedIdx[i];
507  size_t index = i;
508  const std::unordered_set< int >& neighbors = vecAdjTargetFaces[index];
509  if( dMassDefect[index] > 0.0 )
510  {
511  double dMassCorrectU = 0.0;
512  for( auto it : neighbors )
513  dMassCorrectU += dTargetAreas[it] * ( dataUpperBound[it] - dataCorrection[it] );
514 
515  // double dMassDiffCumOld = dMassDefect[index];
516  for( auto it : neighbors )
517  dataCorrection[it] +=
518  dMassDefect[index] * ( dataUpperBound[it] - dataCorrection[it] ) / dMassCorrectU;
519  }
520  else
521  {
522  double dMassCorrectL = 0.0;
523  for( auto it : neighbors )
524  dMassCorrectL += dTargetAreas[it] * ( dataCorrection[it] - dataLowerBound[it] );
525 
526  // double dMassDiffCumOld = dMassDefect[index];
527  for( auto it : neighbors )
528  dataCorrection[it] +=
529  dMassDefect[index] * ( dataCorrection[it] - dataLowerBound[it] ) / dMassCorrectL;
530  }
531  }
532 
533  for( size_t i = 0; i < nrows; i++ )
534  dataCorrectedField[i] += dataCorrection[i];
535  }
536 
537  return dMassDiffCum;
538 }
539 
540 void moab::TempestOnlineMap::CAASLimiter( std::vector< double >& dataCorrectedField,
541  std::vector< double >& dataLowerBound,
542  std::vector< double >& dataUpperBound,
543  double& dMass )
544 {
545  const size_t nrows = dataCorrectedField.size();
546  double dMassL = 0.0;
547  double dMassU = 0.0;
548  std::vector< double > dataCorrection( nrows );
549  const DataArray1D< double >& dTargetAreas = this->m_remapper->m_target->vecFaceArea;
550  double dMassDiff = dMass;
551  double dLMinusU = fabs( dataUpperBound[0] - dataLowerBound[0] );
552  double dMassCorrectU = 0.0;
553  double dMassCorrectL = 0.0;
554  for( size_t i = 0; i < nrows; i++ )
555  {
556  dataCorrection[i] = fmax( dataLowerBound[i], fmin( dataUpperBound[i], 0.0 ) );
557  dMassL += dTargetAreas[i] * dataLowerBound[i];
558  dMassU += dTargetAreas[i] * dataUpperBound[i];
559  dMassDiff -= dTargetAreas[i] * dataCorrection[i];
560  dLMinusU = fmax( dLMinusU, fabs( dataUpperBound[i] - dataLowerBound[i] ) );
561  dMassCorrectL += dTargetAreas[i] * ( dataCorrection[i] - dataLowerBound[i] );
562  dMassCorrectU += dTargetAreas[i] * ( dataUpperBound[i] - dataCorrection[i] );
563  }
564 
565 #ifdef MOAB_HAVE_MPI
566  std::vector< double > localDefects( 5, 0.0 ), globalDefects( 5, 0.0 );
567  localDefects[0] = dMassL;
568  localDefects[1] = dMassU;
569  localDefects[2] = dMassDiff;
570  localDefects[3] = dMassCorrectL;
571  localDefects[4] = dMassCorrectU;
572 
573  MPI_Allreduce( localDefects.data(), globalDefects.data(), 5, MPI_DOUBLE, MPI_SUM, m_pcomm->comm() );
574 
575  dMassL = globalDefects[0];
576  dMassU = globalDefects[1];
577  dMassDiff = globalDefects[2];
578  dMassCorrectL = globalDefects[3];
579  dMassCorrectU = globalDefects[4];
580 #endif
581 
582  //If the upper and lower bounds are too close together, just clip
583  if( fabs( dMassDiff ) < 1e-15 || fabs( dLMinusU ) < 1e-15 )
584  {
585  for( size_t i = 0; i < nrows; i++ )
586  dataCorrectedField[i] += dataCorrection[i];
587  return;
588  }
589  else
590  {
591  if( dMassL > dMassDiff )
592  {
593  Announce( "%d: Lower bound mass exceeds target mass by %1.15e: CAAS will need another iteration", rank,
594  dMassL - dMassDiff );
595  dMassDiff = dMassL;
596  dMass -= dMassL;
597  }
598  else if( dMassU < dMassDiff )
599  {
600  Announce( "%d: Target mass exceeds upper bound mass by %1.15e: CAAS will need another iteration", rank,
601  dMassDiff - dMassU );
602  dMassDiff = dMassU;
603  dMass -= dMassU;
604  }
605 
606  // TODO: optimize away dataMassVec by a simple transient double within the loop
607  DataArray1D< double > dataMassVec( nrows ); //vector of mass redistribution
608  if( dMassDiff > 0.0 )
609  {
610  for( size_t i = 0; i < nrows; i++ )
611  {
612  dataMassVec[i] = ( dataUpperBound[i] - dataCorrection[i] ) / dMassCorrectU;
613  dataCorrection[i] += dMassDiff * dataMassVec[i];
614  }
615  }
616  else
617  {
618  for( size_t i = 0; i < nrows; i++ )
619  {
620  dataMassVec[i] = ( dataCorrection[i] - dataLowerBound[i] ) / dMassCorrectL;
621  dataCorrection[i] += dMassDiff * dataMassVec[i];
622  }
623  }
624 
625  for( size_t i = 0; i < nrows; i++ )
626  dataCorrectedField[i] += dataCorrection[i];
627  }
628 
629  return;
630 }
631 
632 std::pair< double, double > moab::TempestOnlineMap::ApplyBoundsLimiting( std::vector< double >& dataInDouble,
633  std::vector< double >& dataOutDouble,
634  CAASType caasType,
635  int caasIteration,
636  double mismatch )
637 {
638  // Currently only implemented for FV to FV remapping
639  // We should generalize this to other types of remapping
640  assert( !dataGLLNodesSrcCov.IsAttached() && !dataGLLNodesDest.IsAttached() );
641 
642  std::pair< double, double > massDefect( 0.0, 0.0 );
643 
644  // Check if the source and target data are of the same size
645  const size_t nTargetCount = dataOutDouble.size();
646  const DataArray1D< double >& m_dOverlapAreas = this->m_remapper->m_overlap->vecFaceArea;
647 
648  // Apply the offline map to the data
649  double dMassDiff = 0.0;
650  std::vector< double > x( nTargetCount );
651  std::vector< double > dataLowerBound( nTargetCount );
652  std::vector< double > dataUpperBound( nTargetCount );
653  std::vector< double > massVector( nTargetCount );
654  std::vector< std::unordered_set< int > > vecSourceOvTarget( nTargetCount );
655 
656 #undef USE_ComputeAdjacencyRelations
657  constexpr bool useMOABAdjacencies = true;
658 #ifdef USE_ComputeAdjacencyRelations
659  // Compute the adjacent faces to the source face
660  // However, calling MOAB to do this does not work correctly as we need ixS to be the index
661  // Cannot just iterate over all entities in the source covering mesh
662  if( caasType == CAAS_QLT || caasType == CAAS_LOCAL_ADJACENT )
663  {
664  if( useMOABAdjacencies )
665  {
666  moab::ErrorCode rval =
667  ComputeAdjacencyRelations( vecSourceOvTarget, caasIteration, m_remapper->m_covering_source_entities,
668  useMOABAdjacencies );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
669  }
670  else
671  {
672  moab::ErrorCode rval =
673  ComputeAdjacencyRelations( vecSourceOvTarget, caasIteration, m_remapper->m_covering_source_entities,
674  useMOABAdjacencies, m_meshInputCov );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
675  }
676  }
677 #else
678  moab::MeshTopoUtil mtu( m_interface );
679 #endif
680 
681  // Initialize the bounds on the given source and target data
682  double dSourceMin = dataInDouble[0];
683  double dSourceMax = dataInDouble[0];
684  double dTargetMin = dataOutDouble[0];
685  double dTargetMax = dataOutDouble[0];
686  for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
687  {
688  const int ixS = m_meshOverlap->vecSourceFaceIx[i];
689  const int ixT = m_meshOverlap->vecTargetFaceIx[i];
690 
691  if( ixT < 0 ) continue; // skip ghost target faces
692 
693  assert( m_dOverlapAreas[i] > 0.0 );
694  assert( ixS >= 0 );
695  assert( ixT >= 0 );
696 
697 #ifndef USE_ComputeAdjacencyRelations
698  // Compute the adjacent faces to the target face
699  vecSourceOvTarget[ixT].insert( ixS ); // map target face to source face
700  if( ( caasType == CAAS_QLT || caasType == CAAS_LOCAL_ADJACENT ) )
701  {
702  if( useMOABAdjacencies )
703  {
704  moab::Range ents;
705  ents.insert( m_remapper->m_covering_source_entities[ixS] );
706  moab::Range adjEnts;
707  moab::ErrorCode rval = mtu.get_bridge_adjacencies( ents, 0, 2, adjEnts, caasIteration );MB_CHK_SET_ERR_CONT( rval, "Failed to get adjacent faces" );
708  for( moab::Range::iterator it = adjEnts.begin(); it != adjEnts.end(); ++it )
709  {
710  int adjIndex = m_remapper->m_covering_source_entities.index( *it );
711  if( adjIndex >= 0 ) vecSourceOvTarget[ixT].insert( adjIndex );
712  }
713  }
714  else
715  {
716  // Compute the adjacent faces to the target face
717  AdjacentFaceVector vecAdjFaces;
718  GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixS,
719  ( caasIteration ) * ( m_input_order + 1 ) * ( m_input_order + 1 ),
720  vecAdjFaces );
721 
722  //Compute min/max over neighboring faces
723  for( size_t iadj = 0; iadj < vecAdjFaces.size(); iadj++ )
724  vecSourceOvTarget[ixT].insert( vecAdjFaces[iadj].first ); // map target face to source face
725  }
726  }
727 #endif
728 
729  // Update the min and max values of the source data
730  dSourceMax = fmax( dSourceMax, dataInDouble[ixS] );
731  dSourceMin = fmin( dSourceMin, dataInDouble[ixS] );
732 
733  // Update the min and max values of the target data
734  dTargetMin = fmin( dTargetMin, dataOutDouble[ixT] );
735  dTargetMax = fmax( dTargetMax, dataOutDouble[ixT] );
736 
737  const double locMassDiff = ( dataInDouble[ixS] * m_dOverlapAreas[i] ) - // source mass
738  ( dataOutDouble[ixT] * m_dOverlapAreas[i] ); // target mass
739 
740  // Update the mass difference between source and target faces
741  // linked to the overlap mesh element
742  dMassDiff += locMassDiff; // target mass
743  massVector[ixT] += locMassDiff;
744  }
745 
746 #ifdef MOAB_HAVE_MPI
747  std::vector< double > localMinMaxDefects( 5, 0.0 ), globalMinMaxDefects( 5, 0.0 );
748  localMinMaxDefects[0] = dSourceMin;
749  localMinMaxDefects[1] = dTargetMin;
750  localMinMaxDefects[2] = dSourceMax;
751  localMinMaxDefects[3] = dTargetMax;
752  localMinMaxDefects[4] = dMassDiff;
753 
754  if( caasType == CAAS_GLOBAL )
755  {
756  MPI_Allreduce( localMinMaxDefects.data(), globalMinMaxDefects.data(), 2, MPI_DOUBLE, MPI_MIN, m_pcomm->comm() );
757  MPI_Allreduce( localMinMaxDefects.data() + 2, globalMinMaxDefects.data() + 2, 2, MPI_DOUBLE, MPI_MAX,
758  m_pcomm->comm() );
759  dSourceMin = globalMinMaxDefects[0];
760  dSourceMax = globalMinMaxDefects[2];
761  dTargetMin = globalMinMaxDefects[1];
762  dTargetMax = globalMinMaxDefects[3];
763  }
764  if( caasIteration == 1 )
765  MPI_Allreduce( localMinMaxDefects.data() + 4, globalMinMaxDefects.data() + 4, 1, MPI_DOUBLE, MPI_SUM,
766  m_pcomm->comm() );
767  else
768  globalMinMaxDefects[4] = mismatch;
769 
770  dMassDiff = localMinMaxDefects[4];
771  // massDefect.first = localMinMaxDefects[4];
772  massDefect.first = globalMinMaxDefects[4];
773 #else
774 
775  // massDefect.first = fabs( dMassDiff / ( dSourceMax - dSourceMin ) );
776  massDefect.first = dMassDiff;
777 #endif
778 
779  // Early exit if the values are monotone already.
780  // if( ( dTargetMax <= dSourceMax && dTargetMin <= dSourceMin ) || fabs( massDefect.first ) < 1e-16 )
781  if( fabs( massDefect.first ) > 1e-20 )
782  {
783  if( caasType == CAAS_GLOBAL )
784  {
785  for( size_t i = 0; i < nTargetCount; i++ )
786  {
787  dataLowerBound[i] = dSourceMin - dataOutDouble[i];
788  dataUpperBound[i] = dSourceMax - dataOutDouble[i];
789  }
790  } // if( caasType == CAAS_GLOBAL )
791  else // caasType == CAAS_LOCAL
792  {
793  // Compute the local min and max values of the target data
794  std::vector< double > vecLocalUpperBound( nTargetCount );
795  std::vector< double > vecLocalLowerBound( nTargetCount );
796  // Loop over the target faces and compute the min and max values
797  // of the source data linked to the target faces
798  for( size_t i = 0; i < nTargetCount; i++ )
799  {
800  assert( vecSourceOvTarget[i].size() );
801 
802  double dMinI = 1E10; // dataInDouble[vecSourceOvTarget[i][0]];
803  double dMaxI = -1E10; // dataInDouble[vecSourceOvTarget[i][0]];
804 
805  // Compute max over intersecting source faces
806  for( const auto& srcElem : vecSourceOvTarget[i] )
807  {
808  dMinI = fmin( dMinI, dataInDouble[srcElem] ); // min over intersecting source faces
809  dMaxI = fmax( dMaxI, dataInDouble[srcElem] ); // max over intersecting source faces
810  }
811 
812  // Update the min and max values of the target data
813  vecLocalLowerBound[i] = dMinI;
814  vecLocalUpperBound[i] = dMaxI;
815  }
816 
817  for( size_t i = 0; i < nTargetCount; i++ )
818  {
819  dataLowerBound[i] = vecLocalLowerBound[i] - dataOutDouble[i];
820  dataUpperBound[i] = vecLocalUpperBound[i] - dataOutDouble[i];
821  }
822  } // caasType == CAAS_LOCAL
823 
824  // Invoke CAAS or QLT application on the map
825  if( fabs( dMassDiff ) > 1e-20 )
826  {
827  if( caasType == CAAS_QLT )
828  dMassDiff = QLTLimiter( caasIteration, dataOutDouble, dataLowerBound, dataUpperBound, massVector );
829  else
830  CAASLimiter( dataOutDouble, dataLowerBound, dataUpperBound, dMassDiff );
831  }
832 
833  // Announce output mass
834  double dMassDiffPost = 0.0;
835  for( size_t i = 0; i < m_meshOverlap->faces.size(); i++ )
836  {
837  const int ixS = m_meshOverlap->vecSourceFaceIx[i];
838  const int ixT = m_meshOverlap->vecTargetFaceIx[i];
839 
840  if( ixT < 0 ) continue; // skip ghost target faces
841 
842  // Update the mass difference between source and target faces
843  // linked to the overlap mesh element
844  dMassDiffPost += ( dataInDouble[ixS] * m_dOverlapAreas[i] ) - // source mass
845  ( dataOutDouble[ixT] * m_dOverlapAreas[i] ); // target mass
846  }
847  // massDefect.second = fabs( dMassDiffPost / ( dSourceMax - dSourceMin ) );
848  massDefect.second = dMassDiffPost;
849  }
850 
851  // Ideally should perform an AllReduce here to get the global mass difference across all processors
852  // But if we satisfy the constraint on every task, essentially, the global mass difference should be zero!
853  return massDefect;
854 }
855 
856 ///////////////////////////////////////////////////////////////////////////////
857 
858 extern void ForceConsistencyConservation3( const DataArray1D< double >& vecSourceArea,
859  const DataArray1D< double >& vecTargetArea,
860  DataArray2D< double >& dCoeff,
861  bool fMonotone,
862  bool fSparseConstraints = false );
863 
864 ///////////////////////////////////////////////////////////////////////////////
865 
866 extern void ForceIntArrayConsistencyConservation( const DataArray1D< double >& vecSourceArea,
867  const DataArray1D< double >& vecTargetArea,
868  DataArray2D< double >& dCoeff,
869  bool fMonotone );
870 
871 ///////////////////////////////////////////////////////////////////////////////
872 
873 void moab::TempestOnlineMap::LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
874  const DataArray3D< double >& dataGLLJacobian,
875  int nMonotoneType,
876  bool fContinuousIn,
877  bool fNoConservation,
878  bool fSparseConstraints )
879 {
880  // Order of the polynomial interpolant
881  int nP = dataGLLNodes.GetRows();
882 
883  // Order of triangular quadrature rule
884  const int TriQuadRuleOrder = 4;
885 
886  // Triangular quadrature rule
887  TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
888 
889  int TriQuadraturePoints = triquadrule.GetPoints();
890 
891  const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
892 
893  const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
894 
895  // Sample coefficients
896  DataArray2D< double > dSampleCoeff( nP, nP );
897 
898  // GLL Quadrature nodes on quadrilateral elements
899  DataArray1D< double > dG;
900  DataArray1D< double > dW;
901  GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
902 
903  // Announcements
904  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
905  dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
906  if( is_root )
907  {
908  dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
909  dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
910  dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
911  }
912 
913  // Get SparseMatrix represntation of the OfflineMap
914  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
915 
916  // NodeVector from m_meshOverlap
917  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
918  const NodeVector& nodesFirst = m_meshInputCov->nodes;
919 
920  // Vector of source areas
921  DataArray1D< double > vecSourceArea( nP * nP );
922 
923  DataArray1D< double > vecTargetArea;
924  DataArray2D< double > dCoeff;
925 
926 #ifdef VERBOSE
927  std::stringstream sstr;
928  sstr << "remapdata_" << rank << ".txt";
929  std::ofstream output_file( sstr.str() );
930 #endif
931 
932  // Current Overlap Face
933  int ixOverlap = 0;
934 #ifdef VERBOSE
935  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
936 #endif
937  // generic triangle used for area computation, for triangles around the center of overlap face;
938  // used for overlap faces with more than 4 edges;
939  // nodes array will be set for each triangle;
940  // these triangles are not part of the mesh structure, they are just temporary during
941  // aforementioned decomposition.
942  Face faceTri( 3 );
943  NodeVector nodes( 3 );
944  faceTri.SetNode( 0, 0 );
945  faceTri.SetNode( 1, 1 );
946  faceTri.SetNode( 2, 2 );
947 
948  // Loop over all input Faces
949  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
950  {
951  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
952 
953  if( faceFirst.edges.size() != 4 )
954  {
955  _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
956  }
957 #ifdef VERBOSE
958  // Announce computation progress
959  if( ixFirst % outputFrequency == 0 && is_root )
960  {
961  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
962  }
963 #endif
964  // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
965  // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
966  // However, the relation with MOAB and Tempest will go out of the roof
967 
968  // Determine how many overlap Faces and triangles are present
969  int nOverlapFaces = 0;
970  size_t ixOverlapTemp = ixOverlap;
971  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
972  {
973  // if( m_meshOverlap->vecTargetFaceIx[ixOverlapTemp] < 0 ) continue; // skip ghost target faces
974  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
975  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
976 
977  nOverlapFaces++;
978  }
979 
980  // No overlaps
981  if( nOverlapFaces == 0 ) continue;
982 
983  // Allocate remap coefficients array for meshFirst Face
984  DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
985 
986  // Find the local remap coefficients
987  for( int j = 0; j < nOverlapFaces; j++ )
988  {
989  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
990  if( m_meshOverlap->vecFaceArea[ixOverlap + j] < std::numeric_limits<double>::epsilon() ) // machine precision
991  {
992  if (false) { // verbose detailed output about small overlap elements (near machine precision area)
993  Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
994  m_meshOverlap->vecFaceArea[ixOverlap + j] );
995  int n = faceOverlap.edges.size();
996  Announce( "Number nodes: %d", n );
997  for( int k = 0; k < n; k++ )
998  {
999  Node nd = nodesOverlap[faceOverlap[k]];
1000  Announce( "Node %d %d : %1.10e %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
1001  }
1002  }
1003  continue;
1004  }
1005 
1006  int nbEdges = faceOverlap.edges.size();
1007  int nOverlapTriangles = 1;
1008  Node center; // not used if nbEdges == 3
1009  if( nbEdges > 3 )
1010  { // decompose from center in this case
1011  nOverlapTriangles = nbEdges;
1012  for( int k = 0; k < nbEdges; k++ )
1013  {
1014  const Node& node = nodesOverlap[faceOverlap[k]];
1015  center = center + node;
1016  }
1017  center = center / nbEdges;
1018  center = center.Normalized(); // project back on sphere of radius 1
1019  }
1020 
1021  Node node0, node1, node2;
1022  double dTriangleArea;
1023 
1024  // Loop over all sub-triangles of this Overlap Face
1025  for( int k = 0; k < nOverlapTriangles; k++ )
1026  {
1027  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1028  {
1029  node0 = nodesOverlap[faceOverlap[0]];
1030  node1 = nodesOverlap[faceOverlap[1]];
1031  node2 = nodesOverlap[faceOverlap[2]];
1032  dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1033  }
1034  else // decompose polygon in triangles around the center
1035  {
1036  node0 = center;
1037  node1 = nodesOverlap[faceOverlap[k]];
1038  int k1 = ( k + 1 ) % nbEdges;
1039  node2 = nodesOverlap[faceOverlap[k1]];
1040  nodes[0] = center;
1041  nodes[1] = node1;
1042  nodes[2] = node2;
1043  dTriangleArea = CalculateFaceArea( faceTri, nodes );
1044  }
1045  // Coordinates of quadrature Node
1046  for( int l = 0; l < TriQuadraturePoints; l++ )
1047  {
1048  Node nodeQuadrature;
1049  nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
1050  TriQuadratureG[l][2] * node2.x;
1051 
1052  nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
1053  TriQuadratureG[l][2] * node2.y;
1054 
1055  nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
1056  TriQuadratureG[l][2] * node2.z;
1057 
1058  nodeQuadrature = nodeQuadrature.Normalized();
1059 
1060  // Find components of quadrature point in basis
1061  // of the first Face
1062  double dAlpha;
1063  double dBeta;
1064 
1065  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
1066 
1067  // Check inverse map value
1068  if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
1069  ( dBeta > 1.0 + 1.0e-13 ) )
1070  {
1071  _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
1072  "(%1.5e %1.5e)",
1073  j, l, dAlpha, dBeta );
1074  }
1075 
1076  // Sample the finite element at this point
1077  SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
1078 
1079  // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
1080  for( int p = 0; p < nP; p++ )
1081  {
1082  for( int q = 0; q < nP; q++ )
1083  {
1084  dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
1085  m_meshOverlap->vecFaceArea[ixOverlap + j];
1086  }
1087  }
1088  }
1089  }
1090  }
1091 
1092 #ifdef VERBOSE
1093  output_file << "[" << ( ixFirst < (int)col_gdofmap.size() ? (int)col_gdofmap[ixFirst] : -1 ) << "] \t";
1094  for( int j = 0; j < nOverlapFaces; j++ )
1095  {
1096  for( int p = 0; p < nP; p++ )
1097  {
1098  for( int q = 0; q < nP; q++ )
1099  {
1100  output_file << dRemapCoeff[p][q][j] << " ";
1101  }
1102  }
1103  }
1104  output_file << std::endl;
1105 #endif
1106 
1107  // Force consistency and conservation
1108  if( !fNoConservation )
1109  {
1110  double dTargetArea = 0.0;
1111  for( int j = 0; j < nOverlapFaces; j++ )
1112  {
1113  dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
1114  }
1115 
1116  for( int p = 0; p < nP; p++ )
1117  {
1118  for( int q = 0; q < nP; q++ )
1119  {
1120  vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1121  }
1122  }
1123 
1124  const double areaTolerance = 1e-10;
1125  // Source elements are completely covered by target volumes
1126  if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
1127  {
1128  vecTargetArea.Allocate( nOverlapFaces );
1129  for( int j = 0; j < nOverlapFaces; j++ )
1130  {
1131  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1132  }
1133 
1134  dCoeff.Allocate( nOverlapFaces, nP * nP );
1135 
1136  for( int j = 0; j < nOverlapFaces; j++ )
1137  {
1138  for( int p = 0; p < nP; p++ )
1139  {
1140  for( int q = 0; q < nP; q++ )
1141  {
1142  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1143  }
1144  }
1145  }
1146 
1147  // Target volumes only partially cover source elements
1148  }
1149  else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
1150  {
1151  double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
1152 
1153  vecTargetArea.Allocate( nOverlapFaces + 1 );
1154  for( int j = 0; j < nOverlapFaces; j++ )
1155  {
1156  vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
1157  }
1158  vecTargetArea[nOverlapFaces] = dExtraneousArea;
1159 
1160 #ifdef VERBOSE
1161  Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
1162  m_meshInputCov->vecFaceArea[ixFirst] );
1163 #endif
1164  if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
1165  {
1166  _EXCEPTIONT( "Partial element area exceeds total element area" );
1167  }
1168 
1169  dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
1170 
1171  for( int j = 0; j < nOverlapFaces; j++ )
1172  {
1173  for( int p = 0; p < nP; p++ )
1174  {
1175  for( int q = 0; q < nP; q++ )
1176  {
1177  dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
1178  }
1179  }
1180  }
1181  for( int p = 0; p < nP; p++ )
1182  {
1183  for( int q = 0; q < nP; q++ )
1184  {
1185  dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
1186  }
1187  }
1188  for( int j = 0; j < nOverlapFaces; j++ )
1189  {
1190  for( int p = 0; p < nP; p++ )
1191  {
1192  for( int q = 0; q < nP; q++ )
1193  {
1194  dCoeff[nOverlapFaces][p * nP + q] -=
1195  dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
1196  }
1197  }
1198  }
1199  for( int p = 0; p < nP; p++ )
1200  {
1201  for( int q = 0; q < nP; q++ )
1202  {
1203  dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
1204  }
1205  }
1206 
1207  // Source elements only partially cover target volumes
1208  }
1209  else
1210  {
1211  Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
1212  m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
1213  _EXCEPTIONT( "Target grid must be a subset of source grid" );
1214  }
1215 
1216  ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 ),
1217  fSparseConstraints );
1218 
1219  for( int j = 0; j < nOverlapFaces; j++ )
1220  {
1221  for( int p = 0; p < nP; p++ )
1222  {
1223  for( int q = 0; q < nP; q++ )
1224  {
1225  dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
1226  }
1227  }
1228  }
1229  }
1230 
1231 #ifdef VERBOSE
1232  // output_file << "[" << col_gdofmap[ixFirst] << "] \t";
1233  // for ( int j = 0; j < nOverlapFaces; j++ )
1234  // {
1235  // for ( int p = 0; p < nP; p++ )
1236  // {
1237  // for ( int q = 0; q < nP; q++ )
1238  // {
1239  // output_file << dRemapCoeff[p][q][j] << " ";
1240  // }
1241  // }
1242  // }
1243  // output_file << std::endl;
1244 #endif
1245 
1246  // Put these remap coefficients into the SparseMatrix map
1247  for( int j = 0; j < nOverlapFaces; j++ )
1248  {
1249  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1250 
1251  // signal to not participate, because it is a ghost target
1252  if( ixSecondFace < 0 ) continue; // do not do anything
1253 
1254  for( int p = 0; p < nP; p++ )
1255  {
1256  for( int q = 0; q < nP; q++ )
1257  {
1258  if( fContinuousIn )
1259  {
1260  int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
1261 
1262  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1263  m_meshOverlap->vecFaceArea[ixOverlap + j] /
1264  m_meshOutput->vecFaceArea[ixSecondFace];
1265  }
1266  else
1267  {
1268  int ixFirstNode = ixFirst * nP * nP + p * nP + q;
1269 
1270  smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
1271  m_meshOverlap->vecFaceArea[ixOverlap + j] /
1272  m_meshOutput->vecFaceArea[ixSecondFace];
1273  }
1274  }
1275  }
1276  }
1277  // Increment the current overlap index
1278  ixOverlap += nOverlapFaces;
1279  }
1280 #ifdef VERBOSE
1281  output_file.flush(); // required here
1282  output_file.close();
1283 #endif
1284 
1285  return;
1286 }
1287 
1288 ///////////////////////////////////////////////////////////////////////////////
1289 
1290 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
1291  const DataArray3D< double >& dataGLLJacobianIn,
1292  const DataArray3D< int >& dataGLLNodesOut,
1293  const DataArray3D< double >& dataGLLJacobianOut,
1294  const DataArray1D< double >& dataNodalAreaOut,
1295  int nPin,
1296  int nPout,
1297  int nMonotoneType,
1298  bool fContinuousIn,
1299  bool fContinuousOut,
1300  bool fNoConservation )
1301 {
1302  // Triangular quadrature rule
1303  TriangularQuadratureRule triquadrule( 8 );
1304 
1305  const DataArray2D< double >& dG = triquadrule.GetG();
1306  const DataArray1D< double >& dW = triquadrule.GetW();
1307 
1308  // Get SparseMatrix represntation of the OfflineMap
1309  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1310 
1311  // Sample coefficients
1312  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1313  DataArray2D< double > dSampleCoeffOut( nPout, nPout );
1314 
1315  // Announcemnets
1316  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1317  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
1318  if( is_root )
1319  {
1320  dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
1321  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1322  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1323  }
1324 
1325  // Build the integration array for each element on m_meshOverlap
1326  DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
1327 
1328  // Number of overlap Faces per source Face
1329  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1330 
1331  int ixOverlap = 0;
1332  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1333  {
1334  // Determine how many overlap Faces and triangles are present
1335  int nOverlapFaces = 0;
1336  size_t ixOverlapTemp = ixOverlap;
1337  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1338  {
1339  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1340  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
1341  {
1342  break;
1343  }
1344 
1345  nOverlapFaces++;
1346  }
1347 
1348  nAllOverlapFaces[ixFirst] = nOverlapFaces;
1349 
1350  // Increment the current overlap index
1351  ixOverlap += nAllOverlapFaces[ixFirst];
1352  }
1353 
1354  // Geometric area of each output node
1355  DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
1356 
1357  // Area of each overlap element in the output basis
1358  DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
1359 
1360  // Loop through all faces on m_meshInputCov
1361  ixOverlap = 0;
1362 #ifdef VERBOSE
1363  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1364 #endif
1365  if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
1366 
1367  // generic triangle used for area computation, for triangles around the center of overlap face;
1368  // used for overlap faces with more than 4 edges;
1369  // nodes array will be set for each triangle;
1370  // these triangles are not part of the mesh structure, they are just temporary during
1371  // aforementioned decomposition.
1372  Face faceTri( 3 );
1373  NodeVector nodes( 3 );
1374  faceTri.SetNode( 0, 0 );
1375  faceTri.SetNode( 1, 1 );
1376  faceTri.SetNode( 2, 2 );
1377 
1378  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1379  {
1380 #ifdef VERBOSE
1381  // Announce computation progress
1382  if( ixFirst % outputFrequency == 0 && is_root )
1383  {
1384  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1385  }
1386 #endif
1387  // Quantities from the First Mesh
1388  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1389 
1390  const NodeVector& nodesFirst = m_meshInputCov->nodes;
1391 
1392  // Number of overlapping Faces and triangles
1393  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1394 
1395  if( !nOverlapFaces ) continue;
1396 
1397  // // Calculate total element Jacobian
1398  // double dTotalJacobian = 0.0;
1399  // for (int s = 0; s < nPin; s++) {
1400  // for (int t = 0; t < nPin; t++) {
1401  // dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
1402  // }
1403  // }
1404 
1405  // Loop through all Overlap Faces
1406  for( int i = 0; i < nOverlapFaces; i++ )
1407  {
1408  // Quantities from the overlap Mesh
1409  const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
1410 
1411  const NodeVector& nodesOverlap = m_meshOverlap->nodes;
1412 
1413  // Quantities from the Second Mesh
1414  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1415 
1416  // signal to not participate, because it is a ghost target
1417  if( ixSecond < 0 ) continue; // do not do anything
1418 
1419  const NodeVector& nodesSecond = m_meshOutput->nodes;
1420 
1421  const Face& faceSecond = m_meshOutput->faces[ixSecond];
1422 
1423  int nbEdges = faceOverlap.edges.size();
1424  int nOverlapTriangles = 1;
1425  Node center; // not used if nbEdges == 3
1426  if( nbEdges > 3 )
1427  { // decompose from center in this case
1428  nOverlapTriangles = nbEdges;
1429  for( int k = 0; k < nbEdges; k++ )
1430  {
1431  const Node& node = nodesOverlap[faceOverlap[k]];
1432  center = center + node;
1433  }
1434  center = center / nbEdges;
1435  center = center.Normalized(); // project back on sphere of radius 1
1436  }
1437 
1438  Node node0, node1, node2;
1439  double dTriArea;
1440 
1441  // Loop over all sub-triangles of this Overlap Face
1442  for( int j = 0; j < nOverlapTriangles; j++ )
1443  {
1444  if( nbEdges == 3 ) // will come here only once, nOverlapTriangles == 1 in this case
1445  {
1446  node0 = nodesOverlap[faceOverlap[0]];
1447  node1 = nodesOverlap[faceOverlap[1]];
1448  node2 = nodesOverlap[faceOverlap[2]];
1449  dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
1450  }
1451  else // decompose polygon in triangles around the center
1452  {
1453  node0 = center;
1454  node1 = nodesOverlap[faceOverlap[j]];
1455  int j1 = ( j + 1 ) % nbEdges;
1456  node2 = nodesOverlap[faceOverlap[j1]];
1457  nodes[0] = center;
1458  nodes[1] = node1;
1459  nodes[2] = node2;
1460  dTriArea = CalculateFaceArea( faceTri, nodes );
1461  }
1462 
1463  for( int k = 0; k < triquadrule.GetPoints(); k++ )
1464  {
1465  // Get the nodal location of this point
1466  double dX[3];
1467 
1468  dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
1469  dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
1470  dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
1471 
1472  double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
1473 
1474  dX[0] /= dMag;
1475  dX[1] /= dMag;
1476  dX[2] /= dMag;
1477 
1478  Node nodeQuadrature( dX[0], dX[1], dX[2] );
1479 
1480  // Find the components of this quadrature point in the basis
1481  // of the first Face.
1482  double dAlphaIn;
1483  double dBetaIn;
1484 
1485  ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
1486 
1487  // Find the components of this quadrature point in the basis
1488  // of the second Face.
1489  double dAlphaOut;
1490  double dBetaOut;
1491 
1492  ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
1493 
1494  /*
1495  // Check inverse map value
1496  if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
1497  (dBetaIn < 0.0) || (dBetaIn > 1.0)
1498  ) {
1499  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1500  dAlphaIn, dBetaIn);
1501  }
1502 
1503  // Check inverse map value
1504  if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
1505  (dBetaOut < 0.0) || (dBetaOut > 1.0)
1506  ) {
1507  _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
1508  dAlphaOut, dBetaOut);
1509  }
1510  */
1511  // Sample the First finite element at this point
1512  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1513 
1514  // Sample the Second finite element at this point
1515  SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
1516 
1517  // Overlap output area
1518  for( int s = 0; s < nPout; s++ )
1519  {
1520  for( int t = 0; t < nPout; t++ )
1521  {
1522  double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
1523 
1524  dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
1525 
1526  dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
1527  }
1528  }
1529 
1530  // Compute overlap integral
1531  int ixp = 0;
1532  for( int p = 0; p < nPin; p++ )
1533  {
1534  for( int q = 0; q < nPin; q++ )
1535  {
1536  int ixs = 0;
1537  for( int s = 0; s < nPout; s++ )
1538  {
1539  for( int t = 0; t < nPout; t++ )
1540  {
1541  // Sample the Second finite element at this point
1542  dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
1543  dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
1544 
1545  ixs++;
1546  }
1547  }
1548 
1549  ixp++;
1550  }
1551  }
1552  }
1553  }
1554  }
1555 
1556  // Coefficients
1557  DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
1558 
1559  for( int i = 0; i < nOverlapFaces; i++ )
1560  {
1561  // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1562 
1563  int ixp = 0;
1564  for( int p = 0; p < nPin; p++ )
1565  {
1566  for( int q = 0; q < nPin; q++ )
1567  {
1568  int ixs = 0;
1569  for( int s = 0; s < nPout; s++ )
1570  {
1571  for( int t = 0; t < nPout; t++ )
1572  {
1573  dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
1574  dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1575 
1576  ixs++;
1577  }
1578  }
1579 
1580  ixp++;
1581  }
1582  }
1583  }
1584 
1585  // Source areas
1586  DataArray1D< double > vecSourceArea( nPin * nPin );
1587 
1588  for( int p = 0; p < nPin; p++ )
1589  {
1590  for( int q = 0; q < nPin; q++ )
1591  {
1592  vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
1593  }
1594  }
1595 
1596  // Target areas
1597  DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
1598 
1599  for( int i = 0; i < nOverlapFaces; i++ )
1600  {
1601  // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1602  int ixs = 0;
1603  for( int s = 0; s < nPout; s++ )
1604  {
1605  for( int t = 0; t < nPout; t++ )
1606  {
1607  vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
1608 
1609  ixs++;
1610  }
1611  }
1612  }
1613 
1614  // Force consistency and conservation
1615  if( !fNoConservation )
1616  {
1617  ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
1618  }
1619 
1620  // Update global coefficients
1621  for( int i = 0; i < nOverlapFaces; i++ )
1622  {
1623  int ixp = 0;
1624  for( int p = 0; p < nPin; p++ )
1625  {
1626  for( int q = 0; q < nPin; q++ )
1627  {
1628  int ixs = 0;
1629  for( int s = 0; s < nPout; s++ )
1630  {
1631  for( int t = 0; t < nPout; t++ )
1632  {
1633  dGlobalIntArray[ixp][ixOverlap + i][ixs] =
1634  dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
1635 
1636  ixs++;
1637  }
1638  }
1639 
1640  ixp++;
1641  }
1642  }
1643  }
1644 
1645 #ifdef VVERBOSE
1646  // Check column sums (conservation)
1647  for( int i = 0; i < nPin * nPin; i++ )
1648  {
1649  double dColSum = 0.0;
1650  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1651  {
1652  dColSum += dCoeff[j][i] * vecTargetArea[j];
1653  }
1654  printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
1655  }
1656 
1657  // Check row sums (consistency)
1658  for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
1659  {
1660  double dRowSum = 0.0;
1661  for( int i = 0; i < nPin * nPin; i++ )
1662  {
1663  dRowSum += dCoeff[j][i];
1664  }
1665  printf( "Row %i: %1.15e\n", j, dRowSum );
1666  }
1667 #endif
1668 
1669  // Increment the current overlap index
1670  ixOverlap += nOverlapFaces;
1671  }
1672 
1673  // Build redistribution map within target element
1674  if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
1675  DataArray1D< double > dRedistSourceArea( nPout * nPout );
1676  DataArray1D< double > dRedistTargetArea( nPout * nPout );
1677  std::vector< DataArray2D< double > > dRedistributionMaps;
1678  dRedistributionMaps.resize( m_meshOutput->faces.size() );
1679 
1680  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1681  {
1682  dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
1683 
1684  for( int i = 0; i < nPout * nPout; i++ )
1685  {
1686  dRedistributionMaps[ixSecond][i][i] = 1.0;
1687  }
1688 
1689  for( int s = 0; s < nPout * nPout; s++ )
1690  {
1691  dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
1692  }
1693 
1694  for( int s = 0; s < nPout * nPout; s++ )
1695  {
1696  dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
1697  }
1698 
1699  if( !fNoConservation )
1700  {
1701  ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
1702  ( nMonotoneType != 0 ) );
1703 
1704  for( int s = 0; s < nPout * nPout; s++ )
1705  {
1706  for( int t = 0; t < nPout * nPout; t++ )
1707  {
1708  dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
1709  }
1710  }
1711  }
1712  }
1713 
1714  // Construct the total geometric area
1715  DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
1716  for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
1717  {
1718  for( int s = 0; s < nPout; s++ )
1719  {
1720  for( int t = 0; t < nPout; t++ )
1721  {
1722  dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
1723  dGeometricOutputArea[ixSecond][s * nPout + t];
1724  }
1725  }
1726  }
1727 
1728  // Compose the integration operator with the output map
1729  ixOverlap = 0;
1730 
1731  if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
1732 
1733  // Map from source DOFs to target DOFs with redistribution applied
1734  DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
1735 
1736  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1737  {
1738 #ifdef VERBOSE
1739  // Announce computation progress
1740  if( ixFirst % outputFrequency == 0 && is_root )
1741  {
1742  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1743  }
1744 #endif
1745  // Number of overlapping Faces and triangles
1746  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1747 
1748  if( !nOverlapFaces ) continue;
1749 
1750  // Put composed array into map
1751  for( int j = 0; j < nOverlapFaces; j++ )
1752  {
1753  int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
1754 
1755  // signal to not participate, because it is a ghost target
1756  if( ixSecondFace < 0 ) continue; // do not do anything
1757 
1758  dRedistributedOp.Zero();
1759  for( int p = 0; p < nPin * nPin; p++ )
1760  {
1761  for( int s = 0; s < nPout * nPout; s++ )
1762  {
1763  for( int t = 0; t < nPout * nPout; t++ )
1764  {
1765  dRedistributedOp[p][s] +=
1766  dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
1767  }
1768  }
1769  }
1770 
1771  int ixp = 0;
1772  for( int p = 0; p < nPin; p++ )
1773  {
1774  for( int q = 0; q < nPin; q++ )
1775  {
1776  int ixFirstNode;
1777  if( fContinuousIn )
1778  {
1779  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
1780  }
1781  else
1782  {
1783  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
1784  }
1785 
1786  int ixs = 0;
1787  for( int s = 0; s < nPout; s++ )
1788  {
1789  for( int t = 0; t < nPout; t++ )
1790  {
1791  int ixSecondNode;
1792  if( fContinuousOut )
1793  {
1794  ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
1795 
1796  if( !fNoConservation )
1797  {
1798  smatMap( ixSecondNode, ixFirstNode ) +=
1799  dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
1800  }
1801  else
1802  {
1803  smatMap( ixSecondNode, ixFirstNode ) +=
1804  dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
1805  }
1806  }
1807  else
1808  {
1809  ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
1810 
1811  if( !fNoConservation )
1812  {
1813  smatMap( ixSecondNode, ixFirstNode ) +=
1814  dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
1815  }
1816  else
1817  {
1818  smatMap( ixSecondNode, ixFirstNode ) +=
1819  dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
1820  }
1821  }
1822 
1823  ixs++;
1824  }
1825  }
1826 
1827  ixp++;
1828  }
1829  }
1830  }
1831 
1832  // Increment the current overlap index
1833  ixOverlap += nOverlapFaces;
1834  }
1835 
1836  return;
1837 }
1838 
1839 ///////////////////////////////////////////////////////////////////////////////
1840 
1841 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
1842  const DataArray3D< double >& /*dataGLLJacobianIn*/,
1843  const DataArray3D< int >& dataGLLNodesOut,
1844  const DataArray3D< double >& /*dataGLLJacobianOut*/,
1845  const DataArray1D< double >& dataNodalAreaOut,
1846  int nPin,
1847  int nPout,
1848  int nMonotoneType,
1849  bool fContinuousIn,
1850  bool fContinuousOut )
1851 {
1852  // Gauss-Lobatto quadrature within Faces
1853  DataArray1D< double > dGL;
1854  DataArray1D< double > dWL;
1855 
1856  GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
1857 
1858  // Get SparseMatrix represntation of the OfflineMap
1859  SparseMatrix< double >& smatMap = this->GetSparseMatrix();
1860 
1861  // Sample coefficients
1862  DataArray2D< double > dSampleCoeffIn( nPin, nPin );
1863 
1864  // Announcemnets
1865  moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
1866  dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
1867  if( is_root )
1868  {
1869  dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
1870  dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
1871  dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
1872  }
1873 
1874  // Number of overlap Faces per source Face
1875  DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
1876 
1877  int ixOverlap = 0;
1878 
1879  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1880  {
1881  size_t ixOverlapTemp = ixOverlap;
1882  for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
1883  {
1884  // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
1885 
1886  if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
1887 
1888  nAllOverlapFaces[ixFirst]++;
1889  }
1890 
1891  // Increment the current overlap index
1892  ixOverlap += nAllOverlapFaces[ixFirst];
1893  }
1894 
1895  // Number of times this point was found
1896  DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
1897 
1898  ixOverlap = 0;
1899 #ifdef VERBOSE
1900  const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
1901 #endif
1902  // Loop through all faces on m_meshInputCov
1903  for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
1904  {
1905 #ifdef VERBOSE
1906  // Announce computation progress
1907  if( ixFirst % outputFrequency == 0 && is_root )
1908  {
1909  dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
1910  }
1911 #endif
1912  // Quantities from the First Mesh
1913  const Face& faceFirst = m_meshInputCov->faces[ixFirst];
1914 
1915  const NodeVector& nodesFirst = m_meshInputCov->nodes;
1916 
1917  // Number of overlapping Faces and triangles
1918  int nOverlapFaces = nAllOverlapFaces[ixFirst];
1919 
1920  // Loop through all Overlap Faces
1921  for( int i = 0; i < nOverlapFaces; i++ )
1922  {
1923  // Quantities from the Second Mesh
1924  int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
1925 
1926  // signal to not participate, because it is a ghost target
1927  if( ixSecond < 0 ) continue; // do not do anything
1928 
1929  const NodeVector& nodesSecond = m_meshOutput->nodes;
1930  const Face& faceSecond = m_meshOutput->faces[ixSecond];
1931 
1932  // Loop through all nodes on the second face
1933  for( int s = 0; s < nPout; s++ )
1934  {
1935  for( int t = 0; t < nPout; t++ )
1936  {
1937  size_t ixSecondNode;
1938  if( fContinuousOut )
1939  {
1940  ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
1941  }
1942  else
1943  {
1944  ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
1945  }
1946 
1947  if( ixSecondNode >= fSecondNodeFound.GetRows() ) _EXCEPTIONT( "Logic error" );
1948 
1949  // Check if this node has been found already
1950  if( fSecondNodeFound[ixSecondNode] ) continue;
1951 
1952  // Check this node
1953  Node node;
1954  Node dDx1G;
1955  Node dDx2G;
1956 
1957  ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
1958 
1959  // Find the components of this quadrature point in the basis
1960  // of the first Face.
1961  double dAlphaIn;
1962  double dBetaIn;
1963 
1964  ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
1965 
1966  // Check if this node is within the first Face
1967  if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
1968  ( dBetaIn > 1.0 + 1.0e-10 ) )
1969  continue;
1970 
1971  // Node is within the overlap region, mark as found
1972  fSecondNodeFound[ixSecondNode] = true;
1973 
1974  // Sample the First finite element at this point
1975  SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
1976 
1977  // Add to map
1978  for( int p = 0; p < nPin; p++ )
1979  {
1980  for( int q = 0; q < nPin; q++ )
1981  {
1982  int ixFirstNode;
1983  if( fContinuousIn )
1984  {
1985  ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
1986  }
1987  else
1988  {
1989  ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
1990  }
1991 
1992  smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
1993  }
1994  }
1995  }
1996  }
1997  }
1998 
1999  // Increment the current overlap index
2000  ixOverlap += nOverlapFaces;
2001  }
2002 
2003  // Check for missing samples
2004  for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
2005  {
2006  if( !fSecondNodeFound[i] )
2007  {
2008  _EXCEPTION1( "Can't sample point %i", i );
2009  }
2010  }
2011 
2012  return;
2013 }
2014 
2015 ///////////////////////////////////////////////////////////////////////////////