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