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