Commit dcb0b1d6 authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

Merge branch 'P2-2D-prolongation' into 'master'

P2 2D prolongation

See merge request terraneo/tinyhhg!210
parents 45233aea c53f0433
Pipeline #13608 passed with stages
in 42 minutes and 28 seconds
......@@ -102,6 +102,7 @@ public:
ValueType getMaxMagnitude( uint_t level, DoFType flag = All, bool mpiReduce = true ) const;
inline BoundaryCondition getBoundaryCondition() const { return boundaryCondition_; }
inline DoFType getBoundaryTypeToSkipDuringAdditiveCommunication() const { return boundaryTypeToSkipDuringAdditiveCommunication_; }
template< typename SenderType, typename ReceiverType >
inline void startCommunication( const uint_t & level ) const
......
#include "tinyhhg_core/gridtransferoperators/P2toP2QuadraticProlongation.hpp"
#include "tinyhhg_core/gridtransferoperators/generatedKernels/GeneratedKernelsP2MacroFace2D.hpp"
namespace hhg {
void P2toP2QuadraticProlongation::prolongateAdditively( const P2Function< real_t >& function,
const uint_t& sourceLevel,
const DoFType& flag ) const
{
WALBERLA_CHECK_EQUAL( function.getVertexDoFFunction().getBoundaryTypeToSkipDuringAdditiveCommunication(),
function.getEdgeDoFFunction().getBoundaryTypeToSkipDuringAdditiveCommunication() );
WALBERLA_CHECK_EQUAL( flag, function.getVertexDoFFunction().getBoundaryTypeToSkipDuringAdditiveCommunication() ^ All );
WALBERLA_CHECK_EQUAL( flag, function.getEdgeDoFFunction().getBoundaryTypeToSkipDuringAdditiveCommunication() ^ All );
const auto storage = function.getStorage();
const uint_t fineLevel = sourceLevel + 1;
const uint_t coarseLevel = sourceLevel;
function.communicate< Vertex, Edge >( coarseLevel );
function.communicate< Edge, Face >( coarseLevel );
function.interpolate( real_c( 0 ), fineLevel, flag );
for ( const auto& faceIt : function.getStorage()->getFaces() )
{
const auto face = faceIt.second;
auto vertexFineData = face->getData( function.getVertexDoFFunction().getFaceDataID() )->getPointer( fineLevel );
auto edgeFineData = face->getData( function.getEdgeDoFFunction().getFaceDataID() )->getPointer( fineLevel );
const auto vertexCoarseData = face->getData( function.getVertexDoFFunction().getFaceDataID() )->getPointer( coarseLevel );
const auto edgeCoarseData = face->getData( function.getEdgeDoFFunction().getFaceDataID() )->getPointer( coarseLevel );
const double numNeighborFacesEdge0 =
static_cast< double >( storage->getEdge( face->neighborEdges().at( 0 ) )->getNumNeighborFaces() );
const double numNeighborFacesEdge1 =
static_cast< double >( storage->getEdge( face->neighborEdges().at( 1 ) )->getNumNeighborFaces() );
const double numNeighborFacesEdge2 =
static_cast< double >( storage->getEdge( face->neighborEdges().at( 2 ) )->getNumNeighborFaces() );
const double numNeighborFacesVertex0 =
static_cast< double >( storage->getVertex( face->neighborVertices().at( 0 ) )->getNumNeighborFaces() );
const double numNeighborFacesVertex1 =
static_cast< double >( storage->getVertex( face->neighborVertices().at( 1 ) )->getNumNeighborFaces() );
const double numNeighborFacesVertex2 =
static_cast< double >( storage->getVertex( face->neighborVertices().at( 2 ) )->getNumNeighborFaces() );
P2::macroface::generated::prolongate_2D_macroface_P2_push_from_vertexdofs( edgeFineData,
vertexCoarseData,
vertexFineData,
static_cast< int64_t >( coarseLevel ),
numNeighborFacesEdge0,
numNeighborFacesEdge1,
numNeighborFacesEdge2,
numNeighborFacesVertex0,
numNeighborFacesVertex1,
numNeighborFacesVertex2 );
P2::macroface::generated::prolongate_2D_macroface_P2_push_from_edgedofs( edgeCoarseData,
edgeFineData,
vertexFineData,
static_cast< int64_t >( coarseLevel ),
numNeighborFacesEdge0,
numNeighborFacesEdge1,
numNeighborFacesEdge2 );
}
function.getVertexDoFFunction().communicateAdditively< Face, Edge >( fineLevel );
function.getVertexDoFFunction().communicateAdditively< Face, Vertex >( fineLevel );
function.getEdgeDoFFunction().communicateAdditively< Face, Edge >( fineLevel );
}
void P2toP2QuadraticProlongation::prolongateStandard( const P2Function< real_t >& function,
const uint_t& sourceLevel,
const DoFType& flag ) const
{
const auto& vertexDoFFunction = function.getVertexDoFFunction();
const auto& edgeDoFFunction = function.getEdgeDoFFunction();
edgeDoFFunction.template communicate< Vertex, Edge >( sourceLevel );
edgeDoFFunction.template communicate< Edge, Face >( sourceLevel );
vertexDoFFunction.template communicate< Vertex, Edge >( sourceLevel );
vertexDoFFunction.template communicate< Edge, Face >( sourceLevel );
for ( const auto& it : function.getStorage()->getFaces() )
{
const Face& face = *it.second;
const DoFType faceBC = function.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
if ( testFlag( faceBC, flag ) )
{
P2::macroface::prolongate< real_t >(
sourceLevel, face, vertexDoFFunction.getFaceDataID(), edgeDoFFunction.getFaceDataID() );
}
}
for ( const auto& it : function.getStorage()->getEdges() )
{
const Edge& edge = *it.second;
const DoFType edgeBC = function.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if ( testFlag( edgeBC, flag ) )
{
P2::macroedge::prolongate< real_t >(
sourceLevel, edge, vertexDoFFunction.getEdgeDataID(), edgeDoFFunction.getEdgeDataID() );
}
}
for ( const auto& it : function.getStorage()->getVertices() )
{
const Vertex& vertex = *it.second;
const DoFType vertexBC = function.getBoundaryCondition().getBoundaryType( vertex.getMeshBoundaryFlag() );
if ( testFlag( vertexBC, flag ) )
{
P2::macrovertex::prolongate< real_t >(
sourceLevel, vertex, vertexDoFFunction.getVertexDataID(), edgeDoFFunction.getVertexDataID() );
}
}
}
} // namespace hhg
\ No newline at end of file
......@@ -12,52 +12,13 @@ public:
inline void prolongate( const P2Function <real_t> & function, const uint_t & sourceLevel, const DoFType & flag ) const override
{
const auto& vertexDoFFunction = function.getVertexDoFFunction();
const auto& edgeDoFFunction = function.getEdgeDoFFunction();
edgeDoFFunction.template communicate< Vertex, Edge >( sourceLevel );
edgeDoFFunction.template communicate< Edge, Face >( sourceLevel );
vertexDoFFunction.template communicate< Vertex, Edge >( sourceLevel );
vertexDoFFunction.template communicate< Edge, Face >( sourceLevel );
for( const auto& it : function.getStorage()->getFaces() )
{
const Face& face = *it.second;
const DoFType faceBC = function.getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
if( testFlag( faceBC, flag ) )
{
P2::macroface::prolongate< real_t >(
sourceLevel, face, vertexDoFFunction.getFaceDataID(), edgeDoFFunction.getFaceDataID() );
}
}
for( const auto& it : function.getStorage()->getEdges() )
{
const Edge& edge = *it.second;
const DoFType edgeBC = function.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if( testFlag( edgeBC, flag ) )
{
P2::macroedge::prolongate< real_t >(
sourceLevel, edge, vertexDoFFunction.getEdgeDataID(), edgeDoFFunction.getEdgeDataID() );
}
}
for( const auto& it : function.getStorage()->getVertices() )
{
const Vertex& vertex = *it.second;
const DoFType vertexBC = function.getBoundaryCondition().getBoundaryType( vertex.getMeshBoundaryFlag() );
if( testFlag( vertexBC, flag ) )
{
P2::macrovertex::prolongate< real_t >(
sourceLevel, vertex, vertexDoFFunction.getVertexDataID(), edgeDoFFunction.getVertexDataID() );
}
}
prolongateAdditively( function, sourceLevel, flag );
}
private:
void prolongateAdditively( const P2Function <real_t> & function, const uint_t & sourceLevel, const DoFType & flag ) const;
void prolongateStandard( const P2Function <real_t> & function, const uint_t & sourceLevel, const DoFType & flag ) const;
};
}
......@@ -16,6 +16,10 @@ void restrict_2D_macroface_P2_update_vertexdofs(double * _data_edgeFineSrc, doub
void restrict_2D_macroface_P2_update_edgedofs(double * _data_edgeCoarseDst, double * _data_edgeFineSrc, double * _data_vertexFineSrc, int64_t coarse_level, double num_neighbor_faces_edge0, double num_neighbor_faces_edge1, double num_neighbor_faces_edge2);
void prolongate_2D_macroface_P2_push_from_vertexdofs(double * _data_edgeFineDst, double * _data_vertexCoarseSrc, double * _data_vertexFineDst, int64_t coarse_level, double num_neighbor_faces_edge0, double num_neighbor_faces_edge1, double num_neighbor_faces_edge2, double num_neighbor_faces_vertex0, double num_neighbor_faces_vertex1, double num_neighbor_faces_vertex2);
void prolongate_2D_macroface_P2_push_from_edgedofs(double * _data_edgeCoarseSrc, double * _data_edgeFineDst, double * _data_vertexFineDst, int64_t coarse_level, double num_neighbor_faces_edge0, double num_neighbor_faces_edge1, double num_neighbor_faces_edge2);
} // namespace generated
} // namespace macroface
} // namespace P2
......
This source diff could not be displayed because it is too large. You can view the blob instead.
......@@ -119,6 +119,7 @@ class VertexDoFFunction : public Function< VertexDoFFunction< ValueType > >
ValueType getMaxMagnitude( uint_t level, DoFType flag = All, bool mpiReduce = true ) const;
BoundaryCondition getBoundaryCondition() const;
inline DoFType getBoundaryTypeToSkipDuringAdditiveCommunication() const { return boundaryTypeToSkipDuringAdditiveCommunication_; }
template < typename SenderType, typename ReceiverType >
inline void startCommunication( const uint_t& level ) const
......
......@@ -362,56 +362,6 @@ class P2Function : public Function< P2Function< ValueType > >
edgeDoFFunction_.interpolateExtended( expr, edgeDoFFunctions, level, flag );
}
inline void prolongateQuadratic( uint_t sourceLevel, DoFType flag = All )
{
WALBERLA_ABORT( "P2Function - Prolongate (quadratic) not implemented!" );
}
inline void prolongate( uint_t sourceLevel, DoFType flag = All ) const
{
edgeDoFFunction_.template communicate< Vertex, Edge >( sourceLevel );
edgeDoFFunction_.template communicate< Edge, Face >( sourceLevel );
vertexDoFFunction_.template communicate< Vertex, Edge >( sourceLevel );
vertexDoFFunction_.template communicate< Edge, Face >( sourceLevel );
for( const auto& it : this->getStorage()->getFaces() )
{
const Face& face = *it.second;
const DoFType faceBC = this->getBoundaryCondition().getBoundaryType( face.getMeshBoundaryFlag() );
if( testFlag( faceBC, flag ) )
{
P2::macroface::prolongate< ValueType >(
sourceLevel, face, vertexDoFFunction_.getFaceDataID(), edgeDoFFunction_.getFaceDataID() );
}
}
for( const auto& it : this->getStorage()->getEdges() )
{
const Edge& edge = *it.second;
const DoFType edgeBC = this->getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if( testFlag( edgeBC, flag ) )
{
P2::macroedge::prolongate< ValueType >(
sourceLevel, edge, vertexDoFFunction_.getEdgeDataID(), edgeDoFFunction_.getEdgeDataID() );
}
}
for( const auto& it : this->getStorage()->getVertices() )
{
const Vertex& vertex = *it.second;
const DoFType vertexBC = this->getBoundaryCondition().getBoundaryType( vertex.getMeshBoundaryFlag() );
if( testFlag( vertexBC, flag ) )
{
P2::macrovertex::prolongate< ValueType >(
sourceLevel, vertex, vertexDoFFunction_.getVertexDataID(), edgeDoFFunction_.getVertexDataID() );
}
}
}
inline real_t getMaxValue( uint_t level, DoFType flag = All ) const
{
auto localMax = -std::numeric_limits< ValueType >::max();
......
......@@ -5,6 +5,8 @@
#include "tinyhhg_core/p2functionspace/P2Function.hpp"
#include "tinyhhg_core/primitivestorage/SetupPrimitiveStorage.hpp"
#include "tinyhhg_core/gridtransferoperators/P2toP2QuadraticProlongation.hpp"
#include "tinyhhg_core/boundary/BoundaryConditions.hpp"
namespace hhg {
......@@ -17,11 +19,13 @@ static void testP2Prolongate() {
std::shared_ptr<SetupPrimitiveStorage> setupStorage =
std::make_shared<SetupPrimitiveStorage>(mesh, uint_c(walberla::mpi::MPIManager::instance()->numProcesses()));
std::shared_ptr<PrimitiveStorage> storage = std::make_shared<PrimitiveStorage>(*setupStorage);
auto x = std::make_shared<P2Function<real_t> >("x", storage, sourceLevel, sourceLevel + 2);
auto x = std::make_shared<P2Function<real_t> >("x", storage, sourceLevel, sourceLevel + 2, BoundaryCondition::create012BC(), None);
typedef stencilDirection sD;
std::function<real_t(const hhg::Point3D &)> values = [](const hhg::Point3D &) { return 13; };
x->interpolate(values, sourceLevel, hhg::All);
x->prolongate(sourceLevel, hhg::All);
P2toP2QuadraticProlongation prolongationOperator;
prolongationOperator.prolongate( *x, sourceLevel, All );
real_t* edgeDoFFineData = storage->getFace(PrimitiveID(6))->getData(x->getEdgeDoFFunction().getFaceDataID())->getPointer(sourceLevel + 1);
real_t* vertexDoFFineData = storage->getFace(PrimitiveID(6))->getData(x->getVertexDoFFunction().getFaceDataID())->getPointer(sourceLevel + 1);
......@@ -77,8 +81,7 @@ static void testP2Prolongate() {
}
}
x->prolongate(sourceLevel + 1, hhg::All);
prolongationOperator.prolongate( *x, sourceLevel + 1, All );
edgeDoFFineData = storage->getFace(PrimitiveID(6))->getData(x->getEdgeDoFFunction().getFaceDataID())->getPointer(sourceLevel + 2);
vertexDoFFineData = storage->getFace(PrimitiveID(6))->getData(x->getVertexDoFFunction().getFaceDataID())->getPointer(sourceLevel + 2);
......@@ -144,9 +147,11 @@ static void testP2Prolongate2() {
std::shared_ptr<SetupPrimitiveStorage> setupStorage =
std::make_shared<SetupPrimitiveStorage>(mesh, uint_c(walberla::mpi::MPIManager::instance()->numProcesses()));
std::shared_ptr<PrimitiveStorage> storage = std::make_shared<PrimitiveStorage>(*setupStorage);
auto x = std::make_shared<P2Function<real_t> >("x", storage, sourceLevel, sourceLevel + 1);
auto x = std::make_shared<P2Function<real_t> >("x", storage, sourceLevel, sourceLevel + 1, BoundaryCondition::create012BC(), None);
typedef stencilDirection sD;
P2toP2QuadraticProlongation prolongationOperator;
/// this should not be necessary but just to be save
std::function<real_t(const hhg::Point3D &)> zeros = [](const hhg::Point3D &) { return 0; };
x->interpolate(zeros, sourceLevel);
......@@ -171,7 +176,7 @@ static void testP2Prolongate2() {
storage->getFace(PrimitiveID(6))->getData(x->getEdgeDoFFunction().getFaceDataID())->getPointer(
sourceLevel)[hhg::edgedof::macroface::indexFromVertex(sourceLevel,p.first, p.second, stencilDirection::EDGE_VE_N)] = 16.0;
x->prolongate(sourceLevel, hhg::All);
prolongationOperator.prolongate( *x, sourceLevel, All );
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2, p.second *2 + 1, sD::EDGE_VE_N)], 12.);
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2, p.second *2 + 1, sD::EDGE_VE_S)], 12.);
......@@ -199,7 +204,7 @@ static void testP2Prolongate2() {
storage->getFace(PrimitiveID(6))->getData(x->getEdgeDoFFunction().getFaceDataID())->getPointer(
sourceLevel)[hhg::edgedof::macroface::indexFromVertex(sourceLevel,p.first, p.second, stencilDirection::EDGE_HO_E)] = 16.0;
x->prolongate(sourceLevel, hhg::All);
prolongationOperator.prolongate( *x, sourceLevel, All );
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2 + 1, p.second * 2, sD::EDGE_VE_N)], 8.);
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2 + 1, p.second * 2, sD::EDGE_VE_S)], 8.);
......@@ -227,7 +232,7 @@ static void testP2Prolongate2() {
storage->getFace(PrimitiveID(6))->getData(x->getEdgeDoFFunction().getFaceDataID())->getPointer(
sourceLevel)[hhg::edgedof::macroface::indexFromVertex(sourceLevel,p.first, p.second, stencilDirection::EDGE_DI_NE)] = 16.0;
x->prolongate(sourceLevel, hhg::All);
prolongationOperator.prolongate( *x, sourceLevel, All );
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2 + 1, p.second * 2 + 1, sD::EDGE_VE_N)], 8.);
WALBERLA_CHECK_FLOAT_EQUAL(edgeDoFFineData[hhg::edgedof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2 + 1, p.second * 2 + 1, sD::EDGE_VE_S)], 8.);
......@@ -254,8 +259,7 @@ static void testP2Prolongate2() {
storage->getFace(PrimitiveID(6))->getData(x->getVertexDoFFunction().getFaceDataID())->getPointer(
sourceLevel)[hhg::vertexdof::macroface::indexFromVertex(sourceLevel,p.first, p.second, stencilDirection::VERTEX_C)] = 16.0;
x->prolongate(sourceLevel, hhg::All);
prolongationOperator.prolongate( *x, sourceLevel, All );
WALBERLA_CHECK_FLOAT_EQUAL(vertexDoFFineData[hhg::vertexdof::macroface::indexFromVertex(sourceLevel + 1,p.first * 2 , p.second * 2 , sD::VERTEX_C)],
16.,
......@@ -299,17 +303,19 @@ static void testP2InterpolateAndProlongate() {
std::shared_ptr<SetupPrimitiveStorage> setupStorage =
std::make_shared<SetupPrimitiveStorage>(mesh, uint_c(walberla::mpi::MPIManager::instance()->numProcesses()));
std::shared_ptr<PrimitiveStorage> storage = std::make_shared<PrimitiveStorage>(*setupStorage);
auto x = P2Function<real_t> ("x", storage, sourceLevel, targetLevel);
auto y = P2Function<real_t> ("y", storage, sourceLevel, targetLevel);
auto error = P2Function<real_t> ("x", storage, sourceLevel, targetLevel);
auto x = P2Function<real_t> ("x", storage, sourceLevel, targetLevel, BoundaryCondition::create012BC(), None);
auto y = P2Function<real_t> ("y", storage, sourceLevel, targetLevel, BoundaryCondition::create012BC(), None);
auto error = P2Function<real_t> ("x", storage, sourceLevel, targetLevel, BoundaryCondition::create012BC(), None);
std::function<real_t(const hhg::Point3D&)> exact = [](const hhg::Point3D& xx) { return 2. * xx[0] * xx[0] + 3. * xx[0] + 13. + 4. * xx[1] + 5. * xx[1] * xx[1]; };
x.interpolate(exact,sourceLevel ,hhg::All);
y.interpolate(exact,targetLevel,hhg::All);
x.prolongate(sourceLevel ,hhg::All);
x.prolongate(sourceLevel + 1,hhg::All);
x.prolongate(sourceLevel + 2,hhg::All);
P2toP2QuadraticProlongation prolongationOperator;
prolongationOperator.prolongate( x, sourceLevel , All );
prolongationOperator.prolongate( x, sourceLevel + 1, All );
prolongationOperator.prolongate( x, sourceLevel + 2, All );
error.assign({1.0, -1.0}, {x, y}, targetLevel, hhg::All);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment