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