 |
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
|
Mesh Oriented datABase
(version 5.4.1)
Array-based unstructured mesh datastructure
|
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
00031 #endif
00032
00033 #include
00034 #include
00035 #include
00036 #include
00037
00038 // #define VERBOSE
00039
00040 ///
00041 /// Face index and distance metric pair.
00042 ///
00043 typedef std::pair< int, int > FaceDistancePair;
00044
00045 ///
00046 /// Vector storing adjacent Faces.
00047 ///
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 ///////////////////////////////////////////////////////////////////////////////
1.7.6.1.
Dark theme by Tilen Majerle. All rights reserved.