MOAB: Mesh Oriented datABase  (version 5.4.1)
TempestLinearRemap.cpp
Go to the documentation of this file.
00001 ///////////////////////////////////////////////////////////////////////////////
00002 ///
00003 /// \file    TempestLinearRemap.cpp
00004 /// \author  Vijay Mahadevan
00005 /// \version Mar 08, 2017
00006 ///
00007 
00008 #ifdef WIN32               /* windows */
00009 #define _USE_MATH_DEFINES  // For M_PI
00010 #endif
00011 #include "Announce.h"
00012 #include "DataArray3D.h"
00013 #include "FiniteElementTools.h"
00014 #include "FiniteVolumeTools.h"
00015 #include "GaussLobattoQuadrature.h"
00016 #include "TriangularQuadrature.h"
00017 #include "MeshUtilitiesFuzzy.h"
00018 #include "MeshUtilitiesExact.h"
00019 #include "MathHelper.h"
00020 #include "SparseMatrix.h"
00021 #include "OverlapMesh.h"
00022 
00023 #include "DebugOutput.hpp"
00024 #include "moab/AdaptiveKDTree.hpp"
00025 
00026 #include "moab/Remapping/TempestOnlineMap.hpp"
00027 #include "moab/TupleList.hpp"
00028 
00029 #ifdef MOAB_HAVE_EIGEN3
00030 #include <Eigen/Dense>
00031 #endif
00032 
00033 #include <fstream>
00034 #include <cmath>
00035 #include <cstdlib>
00036 #include <sstream>
00037 
00038 // #define VERBOSE
00039 
00040 /// <summary>
00041 ///     Face index and distance metric pair.
00042 /// </summary>
00043 typedef std::pair< int, int > FaceDistancePair;
00044 
00045 /// <summary>
00046 ///     Vector storing adjacent Faces.
00047 /// </summary>
00048 typedef std::vector< FaceDistancePair > AdjacentFaceVector;
00049 
00050 ///////////////////////////////////////////////////////////////////////////////
00051 
00052 moab::ErrorCode moab::TempestOnlineMap::LinearRemapNN_MOAB( bool use_GID_matching, bool strict_check )
00053 {
00054     /* m_mapRemap size = (m_nTotDofs_Dest X m_nTotDofs_SrcCov)  */
00055 
00056 #ifdef VVERBOSE
00057     {
00058         std::ofstream output_file( "rowcolindices.txt", std::ios::out );
00059         output_file << m_nTotDofs_Dest << " " << m_nTotDofs_SrcCov << " " << row_gdofmap.size() << " "
00060                     << row_ldofmap.size() << " " << col_gdofmap.size() << " " << col_ldofmap.size() << "\n";
00061         output_file << "Rows \n";
00062         for( unsigned iv = 0; iv < row_gdofmap.size(); iv++ )
00063             output_file << row_gdofmap[iv] << " " << row_dofmap[iv] << "\n";
00064         output_file << "Cols \n";
00065         for( unsigned iv = 0; iv < col_gdofmap.size(); iv++ )
00066             output_file << col_gdofmap[iv] << " " << col_dofmap[iv] << "\n";
00067         output_file.flush();  // required here
00068         output_file.close();
00069     }
00070 #endif
00071 
00072     if( use_GID_matching )
00073     {
00074         std::map< unsigned, unsigned > src_gl;
00075         for( unsigned it = 0; it < col_gdofmap.size(); ++it )
00076             src_gl[col_gdofmap[it]] = it;
00077 
00078         std::map< unsigned, unsigned >::iterator iter;
00079         for( unsigned it = 0; it < row_gdofmap.size(); ++it )
00080         {
00081             unsigned row = row_gdofmap[it];
00082             iter         = src_gl.find( row );
00083             if( strict_check && iter == src_gl.end() )
00084             {
00085                 std::cout << "Searching for global target DOF " << row
00086                           << " but could not find correspondence in source mesh.\n";
00087                 assert( false );
00088             }
00089             else if( iter == src_gl.end() )
00090             {
00091                 continue;
00092             }
00093             else
00094             {
00095                 unsigned icol = src_gl[row];
00096                 unsigned irow = it;
00097 
00098                 // Set the permutation matrix in local space
00099                 m_mapRemap( irow, icol ) = 1.0;
00100             }
00101         }
00102 
00103         return moab::MB_SUCCESS;
00104     }
00105     else
00106     {
00107         /* Create a Kd-tree to perform local queries to find nearest neighbors */
00108 
00109         return moab::MB_SUCCESS;
00110     }
00111 }
00112 
00113 ///////////////////////////////////////////////////////////////////////////////
00114 
00115 void moab::TempestOnlineMap::LinearRemapFVtoFV_Tempest_MOAB( int nOrder )
00116 {
00117     // Order of triangular quadrature rule
00118     const int TriQuadRuleOrder = 4;
00119 
00120     // Verify ReverseNodeArray has been calculated
00121     if( m_meshInputCov->faces.size() > 0 && m_meshInputCov->revnodearray.size() == 0 )
00122     {
00123         _EXCEPTIONT( "ReverseNodeArray has not been calculated for m_meshInput" );
00124     }
00125 
00126     // Triangular quadrature rule
00127     TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
00128 
00129     // Number of coefficients needed at this order
00130 #ifdef RECTANGULAR_TRUNCATION
00131     int nCoefficients = nOrder * nOrder;
00132 #endif
00133 #ifdef TRIANGULAR_TRUNCATION
00134     int nCoefficients = nOrder * ( nOrder + 1 ) / 2;
00135 #endif
00136 
00137     // Number of faces you need
00138     const int nRequiredFaceSetSize = nCoefficients;
00139 
00140     // Fit weight exponent
00141     const int nFitWeightsExponent = nOrder + 2;
00142 
00143     // Announcemnets
00144     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
00145     dbgprint.set_prefix( "[LinearRemapFVtoFV_Tempest_MOAB]: " );
00146     if( is_root )
00147     {
00148         dbgprint.printf( 0, "Finite Volume to Finite Volume Projection\n" );
00149         dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
00150         dbgprint.printf( 0, "Number of coefficients: %i\n", nCoefficients );
00151         dbgprint.printf( 0, "Required adjacency set size: %i\n", nRequiredFaceSetSize );
00152         dbgprint.printf( 0, "Fit weights exponent: %i\n", nFitWeightsExponent );
00153     }
00154 
00155     // Current overlap face
00156     int ixOverlap                  = 0;
00157     const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
00158 
00159     DataArray2D< double > dIntArray;
00160     DataArray1D< double > dConstraint( nCoefficients );
00161 
00162     // Loop through all faces on m_meshInput
00163     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
00164     {
00165         // Output every 1000 elements
00166         if( ixFirst % outputFrequency == 0 && is_root )
00167         {
00168             dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
00169         }
00170 
00171         // Find the set of Faces that overlap faceFirst
00172         int ixOverlapBegin    = ixOverlap;
00173         unsigned ixOverlapEnd = ixOverlapBegin;
00174 
00175         for( ; ixOverlapEnd < m_meshOverlap->faces.size(); ixOverlapEnd++ )
00176         {
00177             if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapEnd] != 0 )
00178                 break;
00179         }
00180 
00181         unsigned nOverlapFaces = ixOverlapEnd - ixOverlapBegin;
00182 
00183         if( nOverlapFaces == 0 ) continue;
00184 
00185         // Build integration array
00186         BuildIntegrationArray( *m_meshInputCov, *m_meshOverlap, triquadrule, ixFirst, ixOverlapBegin, ixOverlapEnd,
00187                                nOrder, dIntArray );
00188 
00189         // Set of Faces to use in building the reconstruction and associated
00190         // distance metric.
00191         AdjacentFaceVector vecAdjFaces;
00192 
00193         GetAdjacentFaceVectorByEdge( *m_meshInputCov, ixFirst, nRequiredFaceSetSize, vecAdjFaces );
00194 
00195         // Number of adjacent Faces
00196         int nAdjFaces = vecAdjFaces.size();
00197 
00198         // Determine the conservative constraint equation
00199         double dFirstArea = m_meshInputCov->vecFaceArea[ixFirst];
00200         dConstraint.Zero();
00201         for( int p = 0; p < nCoefficients; p++ )
00202         {
00203             for( unsigned j = 0; j < nOverlapFaces; j++ )
00204             {
00205                 dConstraint[p] += dIntArray[p][j];
00206             }
00207             dConstraint[p] /= dFirstArea;
00208         }
00209 
00210         // Build the fit array from the integration operator
00211         DataArray2D< double > dFitArray;
00212         DataArray1D< double > dFitWeights;
00213         DataArray2D< double > dFitArrayPlus;
00214 
00215         BuildFitArray( *m_meshInputCov, triquadrule, ixFirst, vecAdjFaces, nOrder, nFitWeightsExponent, dConstraint,
00216                        dFitArray, dFitWeights );
00217 
00218         // Compute the inverse fit array
00219         bool fSuccess = InvertFitArray_Corrected( dConstraint, dFitArray, dFitWeights, dFitArrayPlus );
00220 
00221         // Multiply integration array and fit array
00222         DataArray2D< double > dComposedArray( nAdjFaces, nOverlapFaces );
00223         if( fSuccess )
00224         {
00225             // Multiply integration array and inverse fit array
00226             for( int i = 0; i < nAdjFaces; i++ )
00227             {
00228                 for( size_t j = 0; j < nOverlapFaces; j++ )
00229                 {
00230                     for( int k = 0; k < nCoefficients; k++ )
00231                     {
00232                         dComposedArray( i, j ) += dIntArray( k, j ) * dFitArrayPlus( i, k );
00233                     }
00234                 }
00235             }
00236 
00237             // Unable to invert fit array, drop to 1st order.  In this case
00238             // dFitArrayPlus(0,0) = 1 and all other entries are zero.
00239         }
00240         else
00241         {
00242             dComposedArray.Zero();
00243             for( size_t j = 0; j < nOverlapFaces; j++ )
00244             {
00245                 dComposedArray( 0, j ) += dIntArray( 0, j );
00246             }
00247         }
00248 
00249         // Put composed array into map
00250         for( unsigned i = 0; i < vecAdjFaces.size(); i++ )
00251         {
00252             for( unsigned j = 0; j < nOverlapFaces; j++ )
00253             {
00254                 int& ixFirstFaceLoc  = vecAdjFaces[i].first;
00255                 int& ixSecondFaceLoc = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
00256                 // int ixFirstFaceGlob = m_remapper->GetGlobalID(moab::Remapper::SourceMesh,
00257                 // ixFirstFaceLoc); int ixSecondFaceGlob =
00258                 // m_remapper->GetGlobalID(moab::Remapper::TargetMesh, ixSecondFaceLoc);
00259 
00260                 // signal to not participate, because it is a ghost target
00261                 if( ixSecondFaceLoc < 0 ) continue;  // do not do anything
00262 
00263                 m_mapRemap( ixSecondFaceLoc, ixFirstFaceLoc ) +=
00264                     dComposedArray[i][j] / m_meshOutput->vecFaceArea[ixSecondFaceLoc];
00265             }
00266         }
00267 
00268         // Increment the current overlap index
00269         ixOverlap += nOverlapFaces;
00270     }
00271 
00272     return;
00273 }
00274 
00275 ///////////////////////////////////////////////////////////////////////////////
00276 
00277 #ifdef MOAB_HAVE_EIGEN3
00278 void moab::TempestOnlineMap::copy_tempest_sparsemat_to_eigen3()
00279 {
00280 #ifndef VERBOSE
00281 #define VERBOSE_ACTIVATED
00282 // #define VERBOSE
00283 #endif
00284     if( m_nTotDofs_Dest <= 0 || m_nTotDofs_SrcCov <= 0 )
00285     {
00286         // std::cout << rank << ": rowsize = " <<  m_nTotDofs_Dest << ", colsize = " <<
00287         // m_nTotDofs_SrcCov << "\n";
00288         // return;  // No need to allocate if either rows or cols size are zero
00289     }
00290 
00291     /* Should the columns be the global size of the matrix ? */
00292     m_weightMatrix.resize( m_nTotDofs_Dest, m_nTotDofs_SrcCov );
00293     m_rowVector.resize( m_weightMatrix.rows() );
00294     m_colVector.resize( m_weightMatrix.cols() );
00295 
00296 #ifdef VERBOSE
00297     int locrows = std::max( m_mapRemap.GetRows(), m_nTotDofs_Dest );
00298     int loccols = std::max( m_mapRemap.GetColumns(), m_nTotDofs_SrcCov );
00299 
00300     std::cout << m_weightMatrix.rows() << ", " << locrows << ", " << m_weightMatrix.cols() << ", " << loccols << "\n";
00301     // assert(m_weightMatrix.rows() == locrows && m_weightMatrix.cols() == loccols);
00302 #endif
00303 
00304     DataArray1D< int > lrows;
00305     DataArray1D< int > lcols;
00306     DataArray1D< double > lvals;
00307     m_mapRemap.GetEntries( lrows, lcols, lvals );
00308     size_t locvals = lvals.GetRows();
00309 
00310     // first matrix
00311     typedef Eigen::Triplet< double > Triplet;
00312     std::vector< Triplet > tripletList;
00313     tripletList.reserve( locvals );
00314     for( size_t iv = 0; iv < locvals; iv++ )
00315     {
00316         tripletList.push_back( Triplet( lrows[iv], lcols[iv], lvals[iv] ) );
00317     }
00318 
00319     m_weightMatrix.setFromTriplets( tripletList.begin(), tripletList.end() );
00320     m_weightMatrix.makeCompressed();
00321 
00322     int nrows = m_weightMatrix.rows();      // Number of rows
00323     int ncols = m_weightMatrix.cols();      // Number of columns
00324     int NNZ   = m_weightMatrix.nonZeros();  // Number of non zero values
00325 #ifdef MOAB_HAVE_MPI
00326     // find out min/max for NNZ, ncols, nrows
00327     // should work on std c++ 11
00328     int arr3[6] = { NNZ, nrows, ncols, -NNZ, -nrows, -ncols };
00329     int rarr3[6];
00330     MPI_Reduce( arr3, rarr3, 6, MPI_INT, MPI_MIN, 0, m_pcomm->comm() );
00331 
00332     int total[2];
00333     MPI_Reduce( arr3, total, 2, MPI_INT, MPI_SUM, 0, m_pcomm->comm() );
00334     if( !rank )
00335         std::cout << " Rows:(" << rarr3[1] << ", " << -rarr3[4] << "), Cols:(" << rarr3[2] << ", " << -rarr3[5]
00336                   << "), NNZ:(" << rarr3[0] << ", " << -rarr3[3] << "),  total NNZ:" << total[0]
00337                   << " total rows:" << total[1] << "\n";
00338 #else
00339     std::cout << "nr rows: " << nrows << " cols: " << ncols << " non-zeros: " << NNZ << "\n";
00340 #endif
00341 
00342 #ifdef VERBOSE
00343     std::stringstream sstr;
00344     sstr << "tempestmatrix.txt.0000" << rank;
00345     std::ofstream output_file( sstr.str(), std::ios::out );
00346     output_file << "0 " << locrows << " 0 " << loccols << "\n";
00347     for( unsigned iv = 0; iv < locvals; iv++ )
00348     {
00349         // output_file << lrows[iv] << " " << row_ldofmap[lrows[iv]] << " " <<
00350         // row_gdofmap[row_ldofmap[lrows[iv]]] << " " << col_gdofmap[col_ldofmap[lcols[iv]]] << " "
00351         // << lvals[iv] << "\n";
00352         output_file << row_gdofmap[row_ldofmap[lrows[iv]]] << " " << col_gdofmap[col_ldofmap[lcols[iv]]] << " "
00353                     << lvals[iv] << "\n";
00354     }
00355     output_file.flush();  // required here
00356     output_file.close();
00357 #endif
00358 
00359 #ifdef VERBOSE_ACTIVATED
00360 #undef VERBOSE_ACTIVATED
00361 #undef VERBOSE
00362 #endif
00363     return;
00364 }
00365 
00366 ///////////////////////////////////////////////////////////////////////////////
00367 //#define VERBOSE
00368 moab::ErrorCode moab::TempestOnlineMap::ApplyWeights( std::vector< double >& srcVals,
00369                                                       std::vector< double >& tgtVals,
00370                                                       bool transpose )
00371 {
00372     // Reset the source and target data first
00373     m_rowVector.setZero();
00374     m_colVector.setZero();
00375 
00376 #ifdef VERBOSE
00377     std::stringstream sstr;
00378     static int callId = 0;
00379     callId++;
00380     sstr << "projection_id_" << callId << "_s_" << size << "_rk_" << rank << ".txt";
00381     std::ofstream output_file( sstr.str() );
00382 #endif
00383     // Perform the actual projection of weights: application of weight matrix onto the source
00384     // solution vector
00385     if( transpose )
00386     {
00387         // Permute the source data first
00388         for( unsigned i = 0; i < srcVals.size(); ++i )
00389         {
00390             if( row_dtoc_dofmap[i] >= 0 )
00391                 m_rowVector( row_dtoc_dofmap[i] ) = srcVals[i];  // permute and set the row (source) vector properly
00392         }
00393 
00394         m_colVector = m_weightMatrix.adjoint() * m_rowVector;
00395 
00396         // Permute the resulting target data back
00397         for( unsigned i = 0; i < tgtVals.size(); ++i )
00398         {
00399             if( col_dtoc_dofmap[i] >= 0 )
00400                 tgtVals[i] = m_colVector( col_dtoc_dofmap[i] );  // permute and set the row (source) vector properly
00401         }
00402     }
00403     else
00404     {
00405         // Permute the source data first
00406 #ifdef VERBOSE
00407         output_file << "ColVector: " << m_colVector.size() << ", SrcVals: " << srcVals.size()
00408                     << ", Sizes: " << m_nTotDofs_SrcCov << ", " << col_dtoc_dofmap.size() << "\n";
00409 #endif
00410         for( unsigned i = 0; i < srcVals.size(); ++i )
00411         {
00412             if( col_dtoc_dofmap[i] >= 0 )
00413             {
00414                 m_colVector( col_dtoc_dofmap[i] ) = srcVals[i];  // permute and set the row (source) vector properly
00415 #ifdef VERBOSE
00416                 output_file << i << " " << col_gdofmap[col_dtoc_dofmap[i]] + 1 << "  " << srcVals[i] << "\n";
00417 #endif
00418             }
00419         }
00420 
00421         m_rowVector = m_weightMatrix * m_colVector;
00422 
00423         // Permute the resulting target data back
00424 #ifdef VERBOSE
00425         output_file << "RowVector: " << m_rowVector.size() << ", TgtVals:" << tgtVals.size()
00426                     << ", Sizes: " << m_nTotDofs_Dest << ", " << row_gdofmap.size() << "\n";
00427 #endif
00428         for( unsigned i = 0; i < tgtVals.size(); ++i )
00429         {
00430             if( row_dtoc_dofmap[i] >= 0 )
00431             {
00432                 tgtVals[i] = m_rowVector( row_dtoc_dofmap[i] );  // permute and set the row (source) vector properly
00433 #ifdef VERBOSE
00434                 output_file << i << " " << row_gdofmap[row_dtoc_dofmap[i]] + 1 << "  " << tgtVals[i] << "\n";
00435 #endif
00436             }
00437         }
00438     }
00439 
00440 #ifdef VERBOSE
00441     output_file.flush();  // required here
00442     output_file.close();
00443 #endif
00444 
00445     // All done with matvec application
00446     return moab::MB_SUCCESS;
00447 }
00448 
00449 #endif
00450 //#undef VERBOSE
00451 ///////////////////////////////////////////////////////////////////////////////
00452 
00453 extern void ForceConsistencyConservation3( const DataArray1D< double >& vecSourceArea,
00454                                            const DataArray1D< double >& vecTargetArea,
00455                                            DataArray2D< double >& dCoeff,
00456                                            bool fMonotone,
00457                                            bool fSparseConstraints = false );
00458 
00459 ///////////////////////////////////////////////////////////////////////////////
00460 
00461 extern void ForceIntArrayConsistencyConservation( const DataArray1D< double >& vecSourceArea,
00462                                                   const DataArray1D< double >& vecTargetArea,
00463                                                   DataArray2D< double >& dCoeff,
00464                                                   bool fMonotone );
00465 
00466 ///////////////////////////////////////////////////////////////////////////////
00467 
00468 void moab::TempestOnlineMap::LinearRemapSE4_Tempest_MOAB( const DataArray3D< int >& dataGLLNodes,
00469                                                           const DataArray3D< double >& dataGLLJacobian,
00470                                                           int nMonotoneType,
00471                                                           bool fContinuousIn,
00472                                                           bool fNoConservation )
00473 {
00474     // Order of the polynomial interpolant
00475     int nP = dataGLLNodes.GetRows();
00476 
00477     // Order of triangular quadrature rule
00478     const int TriQuadRuleOrder = 4;
00479 
00480     // Triangular quadrature rule
00481     TriangularQuadratureRule triquadrule( TriQuadRuleOrder );
00482 
00483     int TriQuadraturePoints = triquadrule.GetPoints();
00484 
00485     const DataArray2D< double >& TriQuadratureG = triquadrule.GetG();
00486 
00487     const DataArray1D< double >& TriQuadratureW = triquadrule.GetW();
00488 
00489     // Sample coefficients
00490     DataArray2D< double > dSampleCoeff( nP, nP );
00491 
00492     // GLL Quadrature nodes on quadrilateral elements
00493     DataArray1D< double > dG;
00494     DataArray1D< double > dW;
00495     GaussLobattoQuadrature::GetPoints( nP, 0.0, 1.0, dG, dW );
00496 
00497     // Announcemnets
00498     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
00499     dbgprint.set_prefix( "[LinearRemapSE4_Tempest_MOAB]: " );
00500     if( is_root )
00501     {
00502         dbgprint.printf( 0, "Finite Element to Finite Volume Projection\n" );
00503         dbgprint.printf( 0, "Triangular quadrature rule order %i\n", TriQuadRuleOrder );
00504         dbgprint.printf( 0, "Order of the FE polynomial interpolant: %i\n", nP );
00505     }
00506 
00507     // Get SparseMatrix represntation of the OfflineMap
00508     SparseMatrix< double >& smatMap = this->GetSparseMatrix();
00509 
00510     // NodeVector from m_meshOverlap
00511     const NodeVector& nodesOverlap = m_meshOverlap->nodes;
00512     const NodeVector& nodesFirst   = m_meshInputCov->nodes;
00513 
00514     // Vector of source areas
00515     DataArray1D< double > vecSourceArea( nP * nP );
00516 
00517     DataArray1D< double > vecTargetArea;
00518     DataArray2D< double > dCoeff;
00519 
00520 #ifdef VERBOSE
00521     std::stringstream sstr;
00522     sstr << "remapdata_" << rank << ".txt";
00523     std::ofstream output_file( sstr.str() );
00524 #endif
00525 
00526     // Current Overlap Face
00527     int ixOverlap                  = 0;
00528     const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
00529 
00530     // generic triangle used for area computation, for triangles around the center of overlap face;
00531     // used for overlap faces with more than 4 edges;
00532     // nodes array will be set for each triangle;
00533     // these triangles are not part of the mesh structure, they are just temporary during
00534     //   aforementioned decomposition.
00535     Face faceTri( 3 );
00536     NodeVector nodes( 3 );
00537     faceTri.SetNode( 0, 0 );
00538     faceTri.SetNode( 1, 1 );
00539     faceTri.SetNode( 2, 2 );
00540 
00541     // Loop over all input Faces
00542     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
00543     {
00544         const Face& faceFirst = m_meshInputCov->faces[ixFirst];
00545 
00546         if( faceFirst.edges.size() != 4 )
00547         {
00548             _EXCEPTIONT( "Only quadrilateral elements allowed for SE remapping" );
00549         }
00550 
00551         // Announce computation progress
00552         if( ixFirst % outputFrequency == 0 && is_root )
00553         {
00554             dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
00555         }
00556 
00557         // Need to re-number the overlap elements such that vecSourceFaceIx[a:b] = 0, then 1 and so
00558         // on wrt the input mesh data Then the overlap_end and overlap_begin will be correct.
00559         // However, the relation with MOAB and Tempest will go out of the roof
00560 
00561         // Determine how many overlap Faces and triangles are present
00562         int nOverlapFaces    = 0;
00563         size_t ixOverlapTemp = ixOverlap;
00564         for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
00565         {
00566             // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
00567             if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
00568             {
00569                 break;
00570             }
00571 
00572             nOverlapFaces++;
00573         }
00574 
00575         // No overlaps
00576         if( nOverlapFaces == 0 )
00577         {
00578             continue;
00579         }
00580 
00581         // Allocate remap coefficients array for meshFirst Face
00582         DataArray3D< double > dRemapCoeff( nP, nP, nOverlapFaces );
00583 
00584         // Find the local remap coefficients
00585         for( int j = 0; j < nOverlapFaces; j++ )
00586         {
00587             const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + j];
00588             if( m_meshOverlap->vecFaceArea[ixOverlap + j] < 1.e-16 )  // machine precision
00589             {
00590                 Announce( "Very small overlap at index %i area polygon: (%1.10e )", ixOverlap + j,
00591                           m_meshOverlap->vecFaceArea[ixOverlap + j] );
00592                 int n = faceOverlap.edges.size();
00593                 Announce( "Number nodes: %d", n );
00594                 for( int k = 0; k < n; k++ )
00595                 {
00596                     Node nd = nodesOverlap[faceOverlap[k]];
00597                     Announce( "Node %d  %d  : %1.10e  %1.10e %1.10e ", k, faceOverlap[k], nd.x, nd.y, nd.z );
00598                 }
00599                 continue;
00600             }
00601 
00602             // #ifdef VERBOSE
00603             // if ( is_root )
00604             //     Announce ( "\tLocal ID: %i/%i = %i, areas = %2.8e", j + ixOverlap, nOverlapFaces,
00605             //     m_remapper->lid_to_gid_covsrc[m_meshOverlap->vecSourceFaceIx[ixOverlap + j]],
00606             //     m_meshOverlap->vecFaceArea[ixOverlap + j] );
00607             // #endif
00608 
00609             int nbEdges           = faceOverlap.edges.size();
00610             int nOverlapTriangles = 1;
00611             Node center;  // not used if nbEdges == 3
00612             if( nbEdges > 3 )
00613             {  // decompose from center in this case
00614                 nOverlapTriangles = nbEdges;
00615                 for( int k = 0; k < nbEdges; k++ )
00616                 {
00617                     const Node& node = nodesOverlap[faceOverlap[k]];
00618                     center           = center + node;
00619                 }
00620                 center = center / nbEdges;
00621                 center = center.Normalized();  // project back on sphere of radius 1
00622             }
00623 
00624             Node node0, node1, node2;
00625             double dTriangleArea;
00626 
00627             // Loop over all sub-triangles of this Overlap Face
00628             for( int k = 0; k < nOverlapTriangles; k++ )
00629             {
00630                 if( nbEdges == 3 )  // will come here only once, nOverlapTriangles == 1 in this case
00631                 {
00632                     node0         = nodesOverlap[faceOverlap[0]];
00633                     node1         = nodesOverlap[faceOverlap[1]];
00634                     node2         = nodesOverlap[faceOverlap[2]];
00635                     dTriangleArea = CalculateFaceArea( faceOverlap, nodesOverlap );
00636                 }
00637                 else  // decompose polygon in triangles around the center
00638                 {
00639                     node0         = center;
00640                     node1         = nodesOverlap[faceOverlap[k]];
00641                     int k1        = ( k + 1 ) % nbEdges;
00642                     node2         = nodesOverlap[faceOverlap[k1]];
00643                     nodes[0]      = center;
00644                     nodes[1]      = node1;
00645                     nodes[2]      = node2;
00646                     dTriangleArea = CalculateFaceArea( faceTri, nodes );
00647                 }
00648                 // Coordinates of quadrature Node
00649                 for( int l = 0; l < TriQuadraturePoints; l++ )
00650                 {
00651                     Node nodeQuadrature;
00652                     nodeQuadrature.x = TriQuadratureG[l][0] * node0.x + TriQuadratureG[l][1] * node1.x +
00653                                        TriQuadratureG[l][2] * node2.x;
00654 
00655                     nodeQuadrature.y = TriQuadratureG[l][0] * node0.y + TriQuadratureG[l][1] * node1.y +
00656                                        TriQuadratureG[l][2] * node2.y;
00657 
00658                     nodeQuadrature.z = TriQuadratureG[l][0] * node0.z + TriQuadratureG[l][1] * node1.z +
00659                                        TriQuadratureG[l][2] * node2.z;
00660 
00661                     nodeQuadrature = nodeQuadrature.Normalized();
00662 
00663                     // Find components of quadrature point in basis
00664                     // of the first Face
00665                     double dAlpha;
00666                     double dBeta;
00667 
00668                     ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlpha, dBeta );
00669 
00670                     // Check inverse map value
00671                     if( ( dAlpha < -1.0e-13 ) || ( dAlpha > 1.0 + 1.0e-13 ) || ( dBeta < -1.0e-13 ) ||
00672                         ( dBeta > 1.0 + 1.0e-13 ) )
00673                     {
00674                         _EXCEPTION4( "Inverse Map for element %d and subtriangle %d out of range "
00675                                      "(%1.5e %1.5e)",
00676                                      j, l, dAlpha, dBeta );
00677                     }
00678 
00679                     // Sample the finite element at this point
00680                     SampleGLLFiniteElement( nMonotoneType, nP, dAlpha, dBeta, dSampleCoeff );
00681 
00682                     // Add sample coefficients to the map if m_meshOverlap->vecFaceArea[ixOverlap + j] > 0
00683                     for( int p = 0; p < nP; p++ )
00684                     {
00685                         for( int q = 0; q < nP; q++ )
00686                         {
00687                             dRemapCoeff[p][q][j] += TriQuadratureW[l] * dTriangleArea * dSampleCoeff[p][q] /
00688                                                     m_meshOverlap->vecFaceArea[ixOverlap + j];
00689                         }
00690                     }
00691                 }
00692             }
00693         }
00694 
00695 #ifdef VERBOSE
00696         output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
00697         for( int j = 0; j < nOverlapFaces; j++ )
00698         {
00699             for( int p = 0; p < nP; p++ )
00700             {
00701                 for( int q = 0; q < nP; q++ )
00702                 {
00703                     output_file << dRemapCoeff[p][q][j] << " ";
00704                 }
00705             }
00706         }
00707         output_file << std::endl;
00708 #endif
00709 
00710         // Force consistency and conservation
00711         if( !fNoConservation )
00712         {
00713             double dTargetArea = 0.0;
00714             for( int j = 0; j < nOverlapFaces; j++ )
00715             {
00716                 dTargetArea += m_meshOverlap->vecFaceArea[ixOverlap + j];
00717             }
00718 
00719             for( int p = 0; p < nP; p++ )
00720             {
00721                 for( int q = 0; q < nP; q++ )
00722                 {
00723                     vecSourceArea[p * nP + q] = dataGLLJacobian[p][q][ixFirst];
00724                 }
00725             }
00726 
00727             const double areaTolerance = 1e-10;
00728             // Source elements are completely covered by target volumes
00729             if( fabs( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea ) <= areaTolerance )
00730             {
00731                 vecTargetArea.Allocate( nOverlapFaces );
00732                 for( int j = 0; j < nOverlapFaces; j++ )
00733                 {
00734                     vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
00735                 }
00736 
00737                 dCoeff.Allocate( nOverlapFaces, nP * nP );
00738 
00739                 for( int j = 0; j < nOverlapFaces; j++ )
00740                 {
00741                     for( int p = 0; p < nP; p++ )
00742                     {
00743                         for( int q = 0; q < nP; q++ )
00744                         {
00745                             dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
00746                         }
00747                     }
00748                 }
00749 
00750                 // Target volumes only partially cover source elements
00751             }
00752             else if( m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea > areaTolerance )
00753             {
00754                 double dExtraneousArea = m_meshInputCov->vecFaceArea[ixFirst] - dTargetArea;
00755 
00756                 vecTargetArea.Allocate( nOverlapFaces + 1 );
00757                 for( int j = 0; j < nOverlapFaces; j++ )
00758                 {
00759                     vecTargetArea[j] = m_meshOverlap->vecFaceArea[ixOverlap + j];
00760                 }
00761                 vecTargetArea[nOverlapFaces] = dExtraneousArea;
00762 
00763 #ifdef VERBOSE
00764                 Announce( "Partial volume: %i (%1.10e / %1.10e)", ixFirst, dTargetArea,
00765                           m_meshInputCov->vecFaceArea[ixFirst] );
00766 #endif
00767                 if( dTargetArea > m_meshInputCov->vecFaceArea[ixFirst] )
00768                 {
00769                     _EXCEPTIONT( "Partial element area exceeds total element area" );
00770                 }
00771 
00772                 dCoeff.Allocate( nOverlapFaces + 1, nP * nP );
00773 
00774                 for( int j = 0; j < nOverlapFaces; j++ )
00775                 {
00776                     for( int p = 0; p < nP; p++ )
00777                     {
00778                         for( int q = 0; q < nP; q++ )
00779                         {
00780                             dCoeff[j][p * nP + q] = dRemapCoeff[p][q][j];
00781                         }
00782                     }
00783                 }
00784                 for( int p = 0; p < nP; p++ )
00785                 {
00786                     for( int q = 0; q < nP; q++ )
00787                     {
00788                         dCoeff[nOverlapFaces][p * nP + q] = dataGLLJacobian[p][q][ixFirst];
00789                     }
00790                 }
00791                 for( int j = 0; j < nOverlapFaces; j++ )
00792                 {
00793                     for( int p = 0; p < nP; p++ )
00794                     {
00795                         for( int q = 0; q < nP; q++ )
00796                         {
00797                             dCoeff[nOverlapFaces][p * nP + q] -=
00798                                 dRemapCoeff[p][q][j] * m_meshOverlap->vecFaceArea[ixOverlap + j];
00799                         }
00800                     }
00801                 }
00802                 for( int p = 0; p < nP; p++ )
00803                 {
00804                     for( int q = 0; q < nP; q++ )
00805                     {
00806                         dCoeff[nOverlapFaces][p * nP + q] /= dExtraneousArea;
00807                     }
00808                 }
00809 
00810                 // Source elements only partially cover target volumes
00811             }
00812             else
00813             {
00814                 Announce( "Coverage area: %1.10e, and target element area: %1.10e)", ixFirst,
00815                           m_meshInputCov->vecFaceArea[ixFirst], dTargetArea );
00816                 _EXCEPTIONT( "Target grid must be a subset of source grid" );
00817             }
00818 
00819             ForceConsistencyConservation3( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType > 0 )
00820                                            /*, m_remapper->lid_to_gid_covsrc[ixFirst]*/ );
00821 
00822             for( int j = 0; j < nOverlapFaces; j++ )
00823             {
00824                 for( int p = 0; p < nP; p++ )
00825                 {
00826                     for( int q = 0; q < nP; q++ )
00827                     {
00828                         dRemapCoeff[p][q][j] = dCoeff[j][p * nP + q];
00829                     }
00830                 }
00831             }
00832         }
00833 
00834 #ifdef VERBOSE
00835         // output_file << "[" << m_remapper->lid_to_gid_covsrc[ixFirst] << "] \t";
00836         // for ( int j = 0; j < nOverlapFaces; j++ )
00837         // {
00838         //     for ( int p = 0; p < nP; p++ )
00839         //     {
00840         //         for ( int q = 0; q < nP; q++ )
00841         //         {
00842         //             output_file << dRemapCoeff[p][q][j] << " ";
00843         //         }
00844         //     }
00845         // }
00846         // output_file << std::endl;
00847 #endif
00848 
00849         // Put these remap coefficients into the SparseMatrix map
00850         for( int j = 0; j < nOverlapFaces; j++ )
00851         {
00852             int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
00853 
00854             // signal to not participate, because it is a ghost target
00855             if( ixSecondFace < 0 ) continue;  // do not do anything
00856 
00857             for( int p = 0; p < nP; p++ )
00858             {
00859                 for( int q = 0; q < nP; q++ )
00860                 {
00861                     if( fContinuousIn )
00862                     {
00863                         int ixFirstNode = dataGLLNodes[p][q][ixFirst] - 1;
00864 
00865                         smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
00866                                                                 m_meshOverlap->vecFaceArea[ixOverlap + j] /
00867                                                                 m_meshOutput->vecFaceArea[ixSecondFace];
00868                     }
00869                     else
00870                     {
00871                         int ixFirstNode = ixFirst * nP * nP + p * nP + q;
00872 
00873                         smatMap( ixSecondFace, ixFirstNode ) += dRemapCoeff[p][q][j] *
00874                                                                 m_meshOverlap->vecFaceArea[ixOverlap + j] /
00875                                                                 m_meshOutput->vecFaceArea[ixSecondFace];
00876                     }
00877                 }
00878             }
00879         }
00880         // Increment the current overlap index
00881         ixOverlap += nOverlapFaces;
00882     }
00883 #ifdef VERBOSE
00884     output_file.flush();  // required here
00885     output_file.close();
00886 #endif
00887 
00888     return;
00889 }
00890 
00891 ///////////////////////////////////////////////////////////////////////////////
00892 
00893 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_MOAB( const DataArray3D< int >& dataGLLNodesIn,
00894                                                         const DataArray3D< double >& dataGLLJacobianIn,
00895                                                         const DataArray3D< int >& dataGLLNodesOut,
00896                                                         const DataArray3D< double >& dataGLLJacobianOut,
00897                                                         const DataArray1D< double >& dataNodalAreaOut,
00898                                                         int nPin,
00899                                                         int nPout,
00900                                                         int nMonotoneType,
00901                                                         bool fContinuousIn,
00902                                                         bool fContinuousOut,
00903                                                         bool fNoConservation )
00904 {
00905     // Triangular quadrature rule
00906     TriangularQuadratureRule triquadrule( 8 );
00907 
00908     const DataArray2D< double >& dG = triquadrule.GetG();
00909     const DataArray1D< double >& dW = triquadrule.GetW();
00910 
00911     // Get SparseMatrix represntation of the OfflineMap
00912     SparseMatrix< double >& smatMap = this->GetSparseMatrix();
00913 
00914     // Sample coefficients
00915     DataArray2D< double > dSampleCoeffIn( nPin, nPin );
00916     DataArray2D< double > dSampleCoeffOut( nPout, nPout );
00917 
00918     // Announcemnets
00919     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
00920     dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_MOAB]: " );
00921     if( is_root )
00922     {
00923         dbgprint.printf( 0, "Finite Element to Finite Element Projection\n" );
00924         dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
00925         dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
00926     }
00927 
00928     // Build the integration array for each element on m_meshOverlap
00929     DataArray3D< double > dGlobalIntArray( nPin * nPin, m_meshOverlap->faces.size(), nPout * nPout );
00930 
00931     // Number of overlap Faces per source Face
00932     DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
00933 
00934     int ixOverlap = 0;
00935     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
00936     {
00937         // Determine how many overlap Faces and triangles are present
00938         int nOverlapFaces    = 0;
00939         size_t ixOverlapTemp = ixOverlap;
00940         for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
00941         {
00942             // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
00943             if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 )
00944             {
00945                 break;
00946             }
00947 
00948             nOverlapFaces++;
00949         }
00950 
00951         nAllOverlapFaces[ixFirst] = nOverlapFaces;
00952 
00953         // Increment the current overlap index
00954         ixOverlap += nAllOverlapFaces[ixFirst];
00955     }
00956 
00957     // Geometric area of each output node
00958     DataArray2D< double > dGeometricOutputArea( m_meshOutput->faces.size(), nPout * nPout );
00959 
00960     // Area of each overlap element in the output basis
00961     DataArray2D< double > dOverlapOutputArea( m_meshOverlap->faces.size(), nPout * nPout );
00962 
00963     // Loop through all faces on m_meshInput
00964     ixOverlap                      = 0;
00965     const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
00966 
00967     if( is_root ) dbgprint.printf( 0, "Building conservative distribution maps\n" );
00968 
00969     // generic triangle used for area computation, for triangles around the center of overlap face;
00970     // used for overlap faces with more than 4 edges;
00971     // nodes array will be set for each triangle;
00972     // these triangles are not part of the mesh structure, they are just temporary during
00973     //   aforementioned decomposition.
00974     Face faceTri( 3 );
00975     NodeVector nodes( 3 );
00976     faceTri.SetNode( 0, 0 );
00977     faceTri.SetNode( 1, 1 );
00978     faceTri.SetNode( 2, 2 );
00979 
00980     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
00981     {
00982         // Announce computation progress
00983         if( ixFirst % outputFrequency == 0 && is_root )
00984         {
00985             dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
00986         }
00987 
00988         // Quantities from the First Mesh
00989         const Face& faceFirst = m_meshInputCov->faces[ixFirst];
00990 
00991         const NodeVector& nodesFirst = m_meshInputCov->nodes;
00992 
00993         // Number of overlapping Faces and triangles
00994         int nOverlapFaces = nAllOverlapFaces[ixFirst];
00995 
00996         if( !nOverlapFaces ) continue;
00997 
00998         // // Calculate total element Jacobian
00999         // double dTotalJacobian = 0.0;
01000         // for (int s = 0; s < nPin; s++) {
01001         //     for (int t = 0; t < nPin; t++) {
01002         //         dTotalJacobian += dataGLLJacobianIn[s][t][ixFirst];
01003         //     }
01004         // }
01005 
01006         // Loop through all Overlap Faces
01007         for( int i = 0; i < nOverlapFaces; i++ )
01008         {
01009             // Quantities from the overlap Mesh
01010             const Face& faceOverlap = m_meshOverlap->faces[ixOverlap + i];
01011 
01012             const NodeVector& nodesOverlap = m_meshOverlap->nodes;
01013 
01014             // Quantities from the Second Mesh
01015             int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
01016 
01017             // signal to not participate, because it is a ghost target
01018             if( ixSecond < 0 ) continue;  // do not do anything
01019 
01020             const NodeVector& nodesSecond = m_meshOutput->nodes;
01021 
01022             const Face& faceSecond = m_meshOutput->faces[ixSecond];
01023 
01024             int nbEdges           = faceOverlap.edges.size();
01025             int nOverlapTriangles = 1;
01026             Node center;  // not used if nbEdges == 3
01027             if( nbEdges > 3 )
01028             {  // decompose from center in this case
01029                 nOverlapTriangles = nbEdges;
01030                 for( int k = 0; k < nbEdges; k++ )
01031                 {
01032                     const Node& node = nodesOverlap[faceOverlap[k]];
01033                     center           = center + node;
01034                 }
01035                 center = center / nbEdges;
01036                 center = center.Normalized();  // project back on sphere of radius 1
01037             }
01038 
01039             Node node0, node1, node2;
01040             double dTriArea;
01041 
01042             // Loop over all sub-triangles of this Overlap Face
01043             for( int j = 0; j < nOverlapTriangles; j++ )
01044             {
01045                 if( nbEdges == 3 )  // will come here only once, nOverlapTriangles == 1 in this case
01046                 {
01047                     node0    = nodesOverlap[faceOverlap[0]];
01048                     node1    = nodesOverlap[faceOverlap[1]];
01049                     node2    = nodesOverlap[faceOverlap[2]];
01050                     dTriArea = CalculateFaceArea( faceOverlap, nodesOverlap );
01051                 }
01052                 else  // decompose polygon in triangles around the center
01053                 {
01054                     node0    = center;
01055                     node1    = nodesOverlap[faceOverlap[j]];
01056                     int j1   = ( j + 1 ) % nbEdges;
01057                     node2    = nodesOverlap[faceOverlap[j1]];
01058                     nodes[0] = center;
01059                     nodes[1] = node1;
01060                     nodes[2] = node2;
01061                     dTriArea = CalculateFaceArea( faceTri, nodes );
01062                 }
01063 
01064                 for( int k = 0; k < triquadrule.GetPoints(); k++ )
01065                 {
01066                     // Get the nodal location of this point
01067                     double dX[3];
01068 
01069                     dX[0] = dG( k, 0 ) * node0.x + dG( k, 1 ) * node1.x + dG( k, 2 ) * node2.x;
01070                     dX[1] = dG( k, 0 ) * node0.y + dG( k, 1 ) * node1.y + dG( k, 2 ) * node2.y;
01071                     dX[2] = dG( k, 0 ) * node0.z + dG( k, 1 ) * node1.z + dG( k, 2 ) * node2.z;
01072 
01073                     double dMag = sqrt( dX[0] * dX[0] + dX[1] * dX[1] + dX[2] * dX[2] );
01074 
01075                     dX[0] /= dMag;
01076                     dX[1] /= dMag;
01077                     dX[2] /= dMag;
01078 
01079                     Node nodeQuadrature( dX[0], dX[1], dX[2] );
01080 
01081                     // Find the components of this quadrature point in the basis
01082                     // of the first Face.
01083                     double dAlphaIn;
01084                     double dBetaIn;
01085 
01086                     ApplyInverseMap( faceFirst, nodesFirst, nodeQuadrature, dAlphaIn, dBetaIn );
01087 
01088                     // Find the components of this quadrature point in the basis
01089                     // of the second Face.
01090                     double dAlphaOut;
01091                     double dBetaOut;
01092 
01093                     ApplyInverseMap( faceSecond, nodesSecond, nodeQuadrature, dAlphaOut, dBetaOut );
01094 
01095                     /*
01096                                         // Check inverse map value
01097                                         if ((dAlphaIn < 0.0) || (dAlphaIn > 1.0) ||
01098                                             (dBetaIn  < 0.0) || (dBetaIn  > 1.0)
01099                                         ) {
01100                                             _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
01101                                                 dAlphaIn, dBetaIn);
01102                                         }
01103 
01104                                         // Check inverse map value
01105                                         if ((dAlphaOut < 0.0) || (dAlphaOut > 1.0) ||
01106                                             (dBetaOut  < 0.0) || (dBetaOut  > 1.0)
01107                                         ) {
01108                                             _EXCEPTION2("Inverse Map out of range (%1.5e %1.5e)",
01109                                                 dAlphaOut, dBetaOut);
01110                                         }
01111                     */
01112                     // Sample the First finite element at this point
01113                     SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
01114 
01115                     // Sample the Second finite element at this point
01116                     SampleGLLFiniteElement( nMonotoneType, nPout, dAlphaOut, dBetaOut, dSampleCoeffOut );
01117 
01118                     // Overlap output area
01119                     for( int s = 0; s < nPout; s++ )
01120                     {
01121                         for( int t = 0; t < nPout; t++ )
01122                         {
01123                             double dNodeArea = dSampleCoeffOut[s][t] * dW[k] * dTriArea;
01124 
01125                             dOverlapOutputArea[ixOverlap + i][s * nPout + t] += dNodeArea;
01126 
01127                             dGeometricOutputArea[ixSecond][s * nPout + t] += dNodeArea;
01128                         }
01129                     }
01130 
01131                     // Compute overlap integral
01132                     int ixp = 0;
01133                     for( int p = 0; p < nPin; p++ )
01134                     {
01135                         for( int q = 0; q < nPin; q++ )
01136                         {
01137                             int ixs = 0;
01138                             for( int s = 0; s < nPout; s++ )
01139                             {
01140                                 for( int t = 0; t < nPout; t++ )
01141                                 {
01142                                     // Sample the Second finite element at this point
01143                                     dGlobalIntArray[ixp][ixOverlap + i][ixs] +=
01144                                         dSampleCoeffOut[s][t] * dSampleCoeffIn[p][q] * dW[k] * dTriArea;
01145 
01146                                     ixs++;
01147                                 }
01148                             }
01149 
01150                             ixp++;
01151                         }
01152                     }
01153                 }
01154             }
01155         }
01156 
01157         // Coefficients
01158         DataArray2D< double > dCoeff( nOverlapFaces * nPout * nPout, nPin * nPin );
01159 
01160         for( int i = 0; i < nOverlapFaces; i++ )
01161         {
01162             // int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
01163 
01164             int ixp = 0;
01165             for( int p = 0; p < nPin; p++ )
01166             {
01167                 for( int q = 0; q < nPin; q++ )
01168                 {
01169                     int ixs = 0;
01170                     for( int s = 0; s < nPout; s++ )
01171                     {
01172                         for( int t = 0; t < nPout; t++ )
01173                         {
01174                             dCoeff[i * nPout * nPout + ixs][ixp] = dGlobalIntArray[ixp][ixOverlap + i][ixs] /
01175                                                                    dOverlapOutputArea[ixOverlap + i][s * nPout + t];
01176 
01177                             ixs++;
01178                         }
01179                     }
01180 
01181                     ixp++;
01182                 }
01183             }
01184         }
01185 
01186         // Source areas
01187         DataArray1D< double > vecSourceArea( nPin * nPin );
01188 
01189         for( int p = 0; p < nPin; p++ )
01190         {
01191             for( int q = 0; q < nPin; q++ )
01192             {
01193                 vecSourceArea[p * nPin + q] = dataGLLJacobianIn[p][q][ixFirst];
01194             }
01195         }
01196 
01197         // Target areas
01198         DataArray1D< double > vecTargetArea( nOverlapFaces * nPout * nPout );
01199 
01200         for( int i = 0; i < nOverlapFaces; i++ )
01201         {
01202             // int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
01203             int ixs = 0;
01204             for( int s = 0; s < nPout; s++ )
01205             {
01206                 for( int t = 0; t < nPout; t++ )
01207                 {
01208                     vecTargetArea[i * nPout * nPout + ixs] = dOverlapOutputArea[ixOverlap + i][nPout * s + t];
01209 
01210                     ixs++;
01211                 }
01212             }
01213         }
01214 
01215         // Force consistency and conservation
01216         if( !fNoConservation )
01217         {
01218             ForceIntArrayConsistencyConservation( vecSourceArea, vecTargetArea, dCoeff, ( nMonotoneType != 0 ) );
01219         }
01220 
01221         // Update global coefficients
01222         for( int i = 0; i < nOverlapFaces; i++ )
01223         {
01224             int ixp = 0;
01225             for( int p = 0; p < nPin; p++ )
01226             {
01227                 for( int q = 0; q < nPin; q++ )
01228                 {
01229                     int ixs = 0;
01230                     for( int s = 0; s < nPout; s++ )
01231                     {
01232                         for( int t = 0; t < nPout; t++ )
01233                         {
01234                             dGlobalIntArray[ixp][ixOverlap + i][ixs] =
01235                                 dCoeff[i * nPout * nPout + ixs][ixp] * dOverlapOutputArea[ixOverlap + i][s * nPout + t];
01236 
01237                             ixs++;
01238                         }
01239                     }
01240 
01241                     ixp++;
01242                 }
01243             }
01244         }
01245 
01246 #ifdef VVERBOSE
01247         // Check column sums (conservation)
01248         for( int i = 0; i < nPin * nPin; i++ )
01249         {
01250             double dColSum = 0.0;
01251             for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
01252             {
01253                 dColSum += dCoeff[j][i] * vecTargetArea[j];
01254             }
01255             printf( "Col %i: %1.15e\n", i, dColSum / vecSourceArea[i] );
01256         }
01257 
01258         // Check row sums (consistency)
01259         for( int j = 0; j < nOverlapFaces * nPout * nPout; j++ )
01260         {
01261             double dRowSum = 0.0;
01262             for( int i = 0; i < nPin * nPin; i++ )
01263             {
01264                 dRowSum += dCoeff[j][i];
01265             }
01266             printf( "Row %i: %1.15e\n", j, dRowSum );
01267         }
01268 #endif
01269 
01270         // Increment the current overlap index
01271         ixOverlap += nOverlapFaces;
01272     }
01273 
01274     // Build redistribution map within target element
01275     if( is_root ) dbgprint.printf( 0, "Building redistribution maps on target mesh\n" );
01276     DataArray1D< double > dRedistSourceArea( nPout * nPout );
01277     DataArray1D< double > dRedistTargetArea( nPout * nPout );
01278     std::vector< DataArray2D< double > > dRedistributionMaps;
01279     dRedistributionMaps.resize( m_meshOutput->faces.size() );
01280 
01281     for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
01282     {
01283         dRedistributionMaps[ixSecond].Allocate( nPout * nPout, nPout * nPout );
01284 
01285         for( int i = 0; i < nPout * nPout; i++ )
01286         {
01287             dRedistributionMaps[ixSecond][i][i] = 1.0;
01288         }
01289 
01290         for( int s = 0; s < nPout * nPout; s++ )
01291         {
01292             dRedistSourceArea[s] = dGeometricOutputArea[ixSecond][s];
01293         }
01294 
01295         for( int s = 0; s < nPout * nPout; s++ )
01296         {
01297             dRedistTargetArea[s] = dataGLLJacobianOut[s / nPout][s % nPout][ixSecond];
01298         }
01299 
01300         if( !fNoConservation )
01301         {
01302             ForceIntArrayConsistencyConservation( dRedistSourceArea, dRedistTargetArea, dRedistributionMaps[ixSecond],
01303                                                   ( nMonotoneType != 0 ) );
01304 
01305             for( int s = 0; s < nPout * nPout; s++ )
01306             {
01307                 for( int t = 0; t < nPout * nPout; t++ )
01308                 {
01309                     dRedistributionMaps[ixSecond][s][t] *= dRedistTargetArea[s] / dRedistSourceArea[t];
01310                 }
01311             }
01312         }
01313     }
01314 
01315     // Construct the total geometric area
01316     DataArray1D< double > dTotalGeometricArea( dataNodalAreaOut.GetRows() );
01317     for( size_t ixSecond = 0; ixSecond < m_meshOutput->faces.size(); ixSecond++ )
01318     {
01319         for( int s = 0; s < nPout; s++ )
01320         {
01321             for( int t = 0; t < nPout; t++ )
01322             {
01323                 dTotalGeometricArea[dataGLLNodesOut[s][t][ixSecond] - 1] +=
01324                     dGeometricOutputArea[ixSecond][s * nPout + t];
01325             }
01326         }
01327     }
01328 
01329     // Compose the integration operator with the output map
01330     ixOverlap = 0;
01331 
01332     if( is_root ) dbgprint.printf( 0, "Assembling map\n" );
01333 
01334     // Map from source DOFs to target DOFs with redistribution applied
01335     DataArray2D< double > dRedistributedOp( nPin * nPin, nPout * nPout );
01336 
01337     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
01338     {
01339         // Announce computation progress
01340         if( ixFirst % outputFrequency == 0 && is_root )
01341         {
01342             dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
01343         }
01344 
01345         // Number of overlapping Faces and triangles
01346         int nOverlapFaces = nAllOverlapFaces[ixFirst];
01347 
01348         if( !nOverlapFaces ) continue;
01349 
01350         // Put composed array into map
01351         for( int j = 0; j < nOverlapFaces; j++ )
01352         {
01353             int ixSecondFace = m_meshOverlap->vecTargetFaceIx[ixOverlap + j];
01354 
01355             // signal to not participate, because it is a ghost target
01356             if( ixSecondFace < 0 ) continue;  // do not do anything
01357 
01358             dRedistributedOp.Zero();
01359             for( int p = 0; p < nPin * nPin; p++ )
01360             {
01361                 for( int s = 0; s < nPout * nPout; s++ )
01362                 {
01363                     for( int t = 0; t < nPout * nPout; t++ )
01364                     {
01365                         dRedistributedOp[p][s] +=
01366                             dRedistributionMaps[ixSecondFace][s][t] * dGlobalIntArray[p][ixOverlap + j][t];
01367                     }
01368                 }
01369             }
01370 
01371             int ixp = 0;
01372             for( int p = 0; p < nPin; p++ )
01373             {
01374                 for( int q = 0; q < nPin; q++ )
01375                 {
01376                     int ixFirstNode;
01377                     if( fContinuousIn )
01378                     {
01379                         ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
01380                     }
01381                     else
01382                     {
01383                         ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
01384                     }
01385 
01386                     int ixs = 0;
01387                     for( int s = 0; s < nPout; s++ )
01388                     {
01389                         for( int t = 0; t < nPout; t++ )
01390                         {
01391                             int ixSecondNode;
01392                             if( fContinuousOut )
01393                             {
01394                                 ixSecondNode = dataGLLNodesOut[s][t][ixSecondFace] - 1;
01395 
01396                                 if( !fNoConservation )
01397                                 {
01398                                     smatMap( ixSecondNode, ixFirstNode ) +=
01399                                         dRedistributedOp[ixp][ixs] / dataNodalAreaOut[ixSecondNode];
01400                                 }
01401                                 else
01402                                 {
01403                                     smatMap( ixSecondNode, ixFirstNode ) +=
01404                                         dRedistributedOp[ixp][ixs] / dTotalGeometricArea[ixSecondNode];
01405                                 }
01406                             }
01407                             else
01408                             {
01409                                 ixSecondNode = ixSecondFace * nPout * nPout + s * nPout + t;
01410 
01411                                 if( !fNoConservation )
01412                                 {
01413                                     smatMap( ixSecondNode, ixFirstNode ) +=
01414                                         dRedistributedOp[ixp][ixs] / dataGLLJacobianOut[s][t][ixSecondFace];
01415                                 }
01416                                 else
01417                                 {
01418                                     smatMap( ixSecondNode, ixFirstNode ) +=
01419                                         dRedistributedOp[ixp][ixs] / dGeometricOutputArea[ixSecondFace][s * nPout + t];
01420                                 }
01421                             }
01422 
01423                             ixs++;
01424                         }
01425                     }
01426 
01427                     ixp++;
01428                 }
01429             }
01430         }
01431 
01432         // Increment the current overlap index
01433         ixOverlap += nOverlapFaces;
01434     }
01435 
01436     return;
01437 }
01438 
01439 ///////////////////////////////////////////////////////////////////////////////
01440 
01441 void moab::TempestOnlineMap::LinearRemapGLLtoGLL2_Pointwise_MOAB( const DataArray3D< int >& dataGLLNodesIn,
01442                                                                   const DataArray3D< double >& /*dataGLLJacobianIn*/,
01443                                                                   const DataArray3D< int >& dataGLLNodesOut,
01444                                                                   const DataArray3D< double >& /*dataGLLJacobianOut*/,
01445                                                                   const DataArray1D< double >& dataNodalAreaOut,
01446                                                                   int nPin,
01447                                                                   int nPout,
01448                                                                   int nMonotoneType,
01449                                                                   bool fContinuousIn,
01450                                                                   bool fContinuousOut )
01451 {
01452     // Gauss-Lobatto quadrature within Faces
01453     DataArray1D< double > dGL;
01454     DataArray1D< double > dWL;
01455 
01456     GaussLobattoQuadrature::GetPoints( nPout, 0.0, 1.0, dGL, dWL );
01457 
01458     // Get SparseMatrix represntation of the OfflineMap
01459     SparseMatrix< double >& smatMap = this->GetSparseMatrix();
01460 
01461     // Sample coefficients
01462     DataArray2D< double > dSampleCoeffIn( nPin, nPin );
01463 
01464     // Announcemnets
01465     moab::DebugOutput dbgprint( std::cout, this->rank, 0 );
01466     dbgprint.set_prefix( "[LinearRemapGLLtoGLL2_Pointwise_MOAB]: " );
01467     if( is_root )
01468     {
01469         dbgprint.printf( 0, "Finite Element to Finite Element (Pointwise) Projection\n" );
01470         dbgprint.printf( 0, "Order of the input FE polynomial interpolant: %i\n", nPin );
01471         dbgprint.printf( 0, "Order of the output FE polynomial interpolant: %i\n", nPout );
01472     }
01473 
01474     // Number of overlap Faces per source Face
01475     DataArray1D< int > nAllOverlapFaces( m_meshInputCov->faces.size() );
01476 
01477     int ixOverlap = 0;
01478 
01479     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
01480     {
01481         size_t ixOverlapTemp = ixOverlap;
01482         for( ; ixOverlapTemp < m_meshOverlap->faces.size(); ixOverlapTemp++ )
01483         {
01484             // const Face & faceOverlap = m_meshOverlap->faces[ixOverlapTemp];
01485 
01486             if( ixFirst - m_meshOverlap->vecSourceFaceIx[ixOverlapTemp] != 0 ) break;
01487 
01488             nAllOverlapFaces[ixFirst]++;
01489         }
01490 
01491         // Increment the current overlap index
01492         ixOverlap += nAllOverlapFaces[ixFirst];
01493     }
01494 
01495     // Number of times this point was found
01496     DataArray1D< bool > fSecondNodeFound( dataNodalAreaOut.GetRows() );
01497 
01498     ixOverlap                      = 0;
01499     const unsigned outputFrequency = ( m_meshInputCov->faces.size() / 10 ) + 1;
01500 
01501     // Loop through all faces on m_meshInputCov
01502     for( size_t ixFirst = 0; ixFirst < m_meshInputCov->faces.size(); ixFirst++ )
01503     {
01504         // Announce computation progress
01505         if( ixFirst % outputFrequency == 0 && is_root )
01506         {
01507             dbgprint.printf( 0, "Element %zu/%lu\n", ixFirst, m_meshInputCov->faces.size() );
01508         }
01509 
01510         // Quantities from the First Mesh
01511         const Face& faceFirst = m_meshInputCov->faces[ixFirst];
01512 
01513         const NodeVector& nodesFirst = m_meshInputCov->nodes;
01514 
01515         // Number of overlapping Faces and triangles
01516         int nOverlapFaces = nAllOverlapFaces[ixFirst];
01517 
01518         // Loop through all Overlap Faces
01519         for( int i = 0; i < nOverlapFaces; i++ )
01520         {
01521 
01522             // Quantities from the Second Mesh
01523             int ixSecond = m_meshOverlap->vecTargetFaceIx[ixOverlap + i];
01524 
01525             // signal to not participate, because it is a ghost target
01526             if( ixSecond < 0 ) continue;  // do not do anything
01527 
01528             const NodeVector& nodesSecond = m_meshOutput->nodes;
01529 
01530             const Face& faceSecond = m_meshOutput->faces[ixSecond];
01531 
01532             // Loop through all nodes on the second face
01533             for( int s = 0; s < nPout; s++ )
01534             {
01535                 for( int t = 0; t < nPout; t++ )
01536                 {
01537                     size_t ixSecondNode;
01538                     if( fContinuousOut )
01539                     {
01540                         ixSecondNode = dataGLLNodesOut[s][t][ixSecond] - 1;
01541                     }
01542                     else
01543                     {
01544                         ixSecondNode = ixSecond * nPout * nPout + s * nPout + t;
01545                     }
01546 
01547                     if( ixSecondNode >= fSecondNodeFound.GetRows() )
01548                     {
01549                         _EXCEPTIONT( "Logic error" );
01550                     }
01551 
01552                     // Check if this node has been found already
01553                     if( fSecondNodeFound[ixSecondNode] )
01554                     {
01555                         continue;
01556                     }
01557 
01558                     // Check this node
01559                     Node node;
01560                     Node dDx1G;
01561                     Node dDx2G;
01562 
01563                     ApplyLocalMap( faceSecond, nodesSecond, dGL[t], dGL[s], node, dDx1G, dDx2G );
01564 
01565                     // Find the components of this quadrature point in the basis
01566                     // of the first Face.
01567                     double dAlphaIn;
01568                     double dBetaIn;
01569 
01570                     ApplyInverseMap( faceFirst, nodesFirst, node, dAlphaIn, dBetaIn );
01571 
01572                     // Check if this node is within the first Face
01573                     if( ( dAlphaIn < -1.0e-10 ) || ( dAlphaIn > 1.0 + 1.0e-10 ) || ( dBetaIn < -1.0e-10 ) ||
01574                         ( dBetaIn > 1.0 + 1.0e-10 ) )
01575                     {
01576                         continue;
01577                     }
01578 
01579                     // Node is within the overlap region, mark as found
01580                     fSecondNodeFound[ixSecondNode] = true;
01581 
01582                     // Sample the First finite element at this point
01583                     SampleGLLFiniteElement( nMonotoneType, nPin, dAlphaIn, dBetaIn, dSampleCoeffIn );
01584 
01585                     // Add to map
01586                     for( int p = 0; p < nPin; p++ )
01587                     {
01588                         for( int q = 0; q < nPin; q++ )
01589                         {
01590                             int ixFirstNode;
01591                             if( fContinuousIn )
01592                             {
01593                                 ixFirstNode = dataGLLNodesIn[p][q][ixFirst] - 1;
01594                             }
01595                             else
01596                             {
01597                                 ixFirstNode = ixFirst * nPin * nPin + p * nPin + q;
01598                             }
01599 
01600                             smatMap( ixSecondNode, ixFirstNode ) += dSampleCoeffIn[p][q];
01601                         }
01602                     }
01603                 }
01604             }
01605         }
01606 
01607         // Increment the current overlap index
01608         ixOverlap += nOverlapFaces;
01609     }
01610 
01611     // Check for missing samples
01612     for( size_t i = 0; i < fSecondNodeFound.GetRows(); i++ )
01613     {
01614         if( !fSecondNodeFound[i] )
01615         {
01616             _EXCEPTION1( "Can't sample point %i", i );
01617         }
01618     }
01619 
01620     return;
01621 }
01622 
01623 ///////////////////////////////////////////////////////////////////////////////
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines