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