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