Skip to content

GitLab

Commit c76d9e2c authored by Maxi Dechant's avatar Maxi Dechant
Browse files

added BoundaryOperatorTest, improved RobinLaplaceApp

parent e7cfe321
Pipeline #36044 failed with stages
in 7 minutes and 23 seconds
Hide whitespace changes
Inline Side-by-side
apps/RobinLaplace/RobinLaplace.cpp
View file @ c76d9e2c
......@@ -24,15 +24,27 @@
#include "hyteg/solvers/GaussSeidelSmoother.hpp"
#include "hyteg/solvers/GeometricMultigridSolver.hpp"
#include "helpers/P1OperatorWithBC.hpp"
#include "helpers/P1OperatorWithBC.cpp"
#include "helpers/P1CustomBoundaryOperator.hpp"
#include "helpers/P1CustomBoundaryOperator.cpp"
//#include "helpers/P1OperatorWithBC.hpp"
//#include "helpers/P1OperatorWithBC.cpp"
//#include "helpers/P1CustomBoundaryOperator.hpp"
//#include "helpers/P1CustomBoundaryOperator.cpp"
#include "hyteg/p1functionspace/p1boundary/P1OperatorWithBC.hpp"
#include "hyteg/p1functionspace/p1boundary/P1OperatorWithBC.cpp"
#include "hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.hpp"
#include "hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.cpp"
using walberla::real_c;
using walberla::real_t;
using namespace hyteg;
bool customBoundaryCheckOfEdge( const hyteg::Point3D& endpoint0, const hyteg::Point3D& endpoint1 )
{
real_t midPointX = (endpoint0[0] + ( endpoint1[0] - endpoint0[0] ) / 2);
real_t midPointY = (endpoint0[1] + ( endpoint1[1] - endpoint0[1] ) / 2);
real_t tol = 0.05;
return (midPointX < tol) || (midPointX > 1 - tol) || (midPointY < tol) || (midPointY > 1 - tol);
}
void generateRightRobinSide( const std::shared_ptr< PrimitiveStorage >& storage,
hyteg::P1Function< real_t >& rhsFunction,
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
......@@ -51,8 +63,12 @@ void generateRightRobinSide( const std::shared_ptr< PrimitiveStorage >&
hyteg::Edge& edge = *storage->getEdge( edgeIDs[uint_c( i )] );
const DoFType edgeBC = rhsFunction.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if ( testFlag( edgeBC, flag ) )
// if ( testFlag( edgeBC, flag ) )
if( customBoundaryCheckOfEdge( edge.getCoordinates()[0], edge.getCoordinates()[1] ) )
{
PrimitiveDataID< FunctionMemory< real_t >, Edge > rhsEdgeDoFIdx = rhsFunction.getEdgeDataID();
auto rhsData = edge.getData( rhsEdgeDoFIdx )->getPointer( level );
for( const auto& microEdgeIndices: hyteg::indexing::MicroEdgeIterator( levelinfo::num_microedges_per_edge( level ) ))
{
real_t edgeNum = real_t( levelinfo::num_microedges_per_edge( level ) );
......@@ -70,18 +86,16 @@ void generateRightRobinSide( const std::shared_ptr< PrimitiveStorage >&
real_t boundaryValue = boundaryFunction( currentNodeCoord );
integralMicro0 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[j];
integralMicro1 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[j]);
integralMicro1 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[j]);
}
PrimitiveDataID< FunctionMemory< real_t >, Edge > rhsEdgeDoFIdx = rhsFunction.getEdgeDataID();
auto rhsData = edge.getData( rhsEdgeDoFIdx )->getPointer( level );
rhsData[ microEdgeIndices[0].col() ] += integralMicro0;
rhsData[ microEdgeIndices[1].col() ] += integralMicro1;
}
rhsFunction.communicateAdditively< Edge, Vertex >( level, hyteg::NeumannBoundary, *storage, false );
rhsFunction.communicateAdditively< Edge, Face >( level, hyteg::NeumannBoundary, *storage, false );
}
}
}
rhsFunction.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage, false );
rhsFunction.communicate< Vertex, Edge >( level );
rhsFunction.communicate< Edge, Face >( level );
}
int main( int argc, char** argv )
......@@ -115,9 +129,59 @@ int main( int argc, char** argv )
std::function< real_t( const hyteg::Point3D& ) > exactSolutionInterpolate = [] ( const hyteg::Point3D& x )
{
return std::pow( x[0], 2 ) - x[0]*x[1] - std::pow( x[1], 2 ) + x[1] + 1.0;
};
std::function< std::array< real_t, 2 >( const hyteg::Point3D& ) > exactSolutionGradient = [] ( const hyteg::Point3D& x )
{
real_t gradX = 2 * x[0] - x[1];
real_t gradY = -x[0] - 2 * x[1] + 1;
std::array< real_t, 2 > response = { gradX, gradY };
return response;
};
/* std::function< real_t( const hyteg::Point3D& ) > exactSolutionInterpolate = [] ( const hyteg::Point3D& x )
{
return 2 * x[0] + x[1] + 1;
};
std::function< real_t( const hyteg::Point3D& ) > boundaryFunction = [ robinFactor ]( const hyteg::Point3D& x )
std::function< std::array< real_t, 2 >( const hyteg::Point3D& ) > exactSolutionGradient = [] ( const hyteg::Point3D& x )
{
std::array< real_t, 2 > response = { 2.0, 1.0 };
return response;
}; */
std::function< real_t( const hyteg::Point3D& ) > boundaryFunction = [robinFactor, exactSolutionInterpolate, exactSolutionGradient ] ( const hyteg::Point3D& x )
{
real_t boundaryTOL = 0.05;
real_t boundaryEvaluate = 0.0;
std::array< real_t, 2 > boundaryNormal = { 0.0, 0.0 };
if( x[0] >= 1-boundaryTOL )
{
boundaryNormal[0] = 1;
boundaryNormal[1] = 0;
}
else if( x[1] >= 1-boundaryTOL )
{
boundaryNormal[0] = 0;
boundaryNormal[1] = 1;
}
else if( x[0] <= boundaryTOL )
{
boundaryNormal[0] = -1;
boundaryNormal[1] = 0;
}
else if( x[1] <= boundaryTOL )
{
boundaryNormal[0] = 0;
boundaryNormal[1] = -1;
}
std::array< real_t, 2 > gradient = exactSolutionGradient( x );
boundaryEvaluate = boundaryNormal[0] * gradient[0] + boundaryNormal[1] * gradient[1] + robinFactor * exactSolutionInterpolate( x );
std::cout << "boundaryFunction( " << x << " ) = " << boundaryEvaluate << std::endl;
return boundaryEvaluate;
};
/* std::function< real_t( const hyteg::Point3D& ) > boundaryFunction = [ robinFactor ]( const hyteg::Point3D& x )
{
real_t boundaryTOL = 0.05;
real_t boundaryEvaluate = 0.0;
......@@ -136,9 +200,9 @@ int main( int argc, char** argv )
else if( x[1] <= boundaryTOL )
{
boundaryEvaluate = x[0] - 1.0 + robinFactor * (x[0]*x[0] + 1.0);
}
}
return boundaryEvaluate;
};
}; */
std::vector< real_t > boundaryQuadNodes = { 0.5 };
std::vector< real_t > boundaryQuadWeights = { 1.0 };
......@@ -149,6 +213,7 @@ int main( int argc, char** argv )
P1OperatorWithBC< P1ConstantLaplaceOperator, P1CustomBoundaryOperator > CustomLaplaceOperator( storage, minLevel, maxLevel, innerOperator, boundaryOperator );
hyteg::P1Function< real_t > f( "f", storage, minLevel, maxLevel );
f.interpolate( 0.0, maxLevel, hyteg::NeumannBoundary );
generateRightRobinSide( storage, f, boundaryFunction, maxLevel, hyteg::NeumannBoundary );
std::vector< PrimitiveID > vertexIDs = storage->getVertexIDs();
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
......@@ -175,6 +240,8 @@ int main( int argc, char** argv )
hyteg::MinResSolver minresSolver = hyteg::MinResSolver< P1OperatorWithBC< hyteg::P1ConstantLaplaceOperator, P1CustomBoundaryOperator > >( storage, minLevel, maxLevel, maxKrylowDim );
minresSolver.solve( CustomLaplaceOperator, u, f, maxLevel );
std::cout << "\n\nSolver finished" << std::endl;
if( useVTK )
{
vtkOutput.write( maxLevel, 1 );
......@@ -182,6 +249,10 @@ int main( int argc, char** argv )
hyteg::P1Function< real_t > exactSolution( "exactSolution", storage, minLevel, maxLevel );
exactSolution.interpolate( exactSolutionInterpolate, maxLevel, hyteg::All );
hyteg::P1Function< real_t > computedRHS( "computedRHS", storage, minLevel, maxLevel );
CustomLaplaceOperator.apply( u, computedRHS, maxLevel, hyteg::All );
std::cout << "\n\nFinished!" << std::endl;
std::cout << "sup( exactSolution ) = " << exactSolution.getMaxMagnitude( maxLevel, hyteg::All, false ) << std::endl;
std::cout << "sup( u|complete ) = " << u.getMaxMagnitude( maxLevel, hyteg::All, false ) << std::endl;
std::cout << "sup( u|boundary ) = " << u.getMaxMagnitude( maxLevel, hyteg::NeumannBoundary, false ) << std::endl;
......@@ -192,10 +263,15 @@ int main( int argc, char** argv )
real_t numDoFs = real_t( hyteg::numberOfGlobalDoFs<P1FunctionTag>( *storage, maxLevel));
std::cout << " p1-norm per DoF of residual (complete) : " << (residual.dotGlobal( residual, maxLevel, hyteg::All ) / numDoFs) << std::endl;
std::vector< PrimitiveID > vertexIDsDebug = storage->getVertexIDs();
for ( int i = 0; i < int_c( vertexIDsDebug.size() ); i++ )
hyteg::P1Function< real_t > rhsResidual( "rhsResidual", storage, minLevel, maxLevel );
rhsResidual.assign( {1.0, -1.0}, {f, computedRHS}, maxLevel, hyteg::All );
std::cout << "\nsup norm of rhsResidual : " << (rhsResidual.dotGlobal( rhsResidual, maxLevel, hyteg::All )) << std::endl;
std::cout << "\nSolutions on vertices:" << std::endl;
// std::vector< PrimitiveID > vertexIDsDebug = storage->getVertexIDs();
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage->getVertex( vertexIDsDebug[uint_c( i )] );
Vertex& vertex = *storage->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > exactVertexDoFIdx = exactSolution.getVertexDataID();
PrimitiveDataID< FunctionMemory< real_t >, Vertex > uVertexDoFIdx = u.getVertexDataID();
auto exactData = vertex.getData( exactVertexDoFIdx )->getPointer( maxLevel );
......@@ -203,9 +279,37 @@ int main( int argc, char** argv )
real_t exactAtVertex = exactData[0];
real_t uAtVertex = uData[0];
std::cout << " vertexID = " << vertexIDs[uint_c( i )] << " , vertex.getCoordinates() = " << vertex.getCoordinates() <<
" : exactSolution = " << exactAtVertex << " , u = " << uAtVertex << std::endl;
}
std::cout << "\nResiduals on vertices:" << std::endl;
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > resVertexDoFIdx = residual.getVertexDataID();
PrimitiveDataID< FunctionMemory< real_t >, Vertex > rhsResVertexDoFIdx = rhsResidual.getVertexDataID();
auto resData = vertex.getData( resVertexDoFIdx )->getPointer( maxLevel );
auto rhsResData = vertex.getData( rhsResVertexDoFIdx )->getPointer( maxLevel );
real_t resAtVertex = resData[0];
real_t rhsResAtVertex = rhsResData[0];
std::cout << " vertexID = " << vertexIDs[uint_c( i )] << " , vertex.getCoordinates() = " << vertex.getCoordinates() <<
" : residual = " << resAtVertex << " , rhsResidual = " << rhsResAtVertex << std::endl;
}
/* hyteg::P1Function< real_t > operatorExact( "operatorExact", storage, minLevel, maxLevel );
CustomLaplaceOperator.apply( exactSolution, operatorExact, maxLevel, hyteg::All);
for ( int i = 0; i < int_c( vertexIDsDebug.size() ); i++ )
{
Vertex& vertex = *storage->getVertex( vertexIDsDebug[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > exactOpVertexDoFIdx = operatorExact.getVertexDataID();
auto exactOpData = vertex.getData( exactOpVertexDoFIdx )->getPointer( maxLevel );
real_t exactOpAtVertex = exactOpData[0];
std::cout << "vertexID = " << vertexIDsDebug[uint_c( i )] << " , vertex.getCoordinates() = " << vertex.getCoordinates() <<
" : exactSolution( vertex.getCoordinates() ) = " << exactAtVertex << " , u( vertex.getCoordinates() ) = " << uAtVertex << std::endl;
}
" : exactOperator( vertex.getCoordinates() ) = " << exactOpAtVertex << std::endl;
} */
}
......
apps/RobinLaplace/helpers/boundaryForms/P1BoundaryForm.hpp deleted 100644 → 0
View file @ e7cfe321
#include "hyteg/forms/Form.hpp"
namespace hyteg {
class P1Form : public Form
{
public:
virtual ~P1Form() {}
// 2D P1
virtual void integrateMicroEdge( const std::array< Point3D, 2 >& mircoEdgecoords,
Matrixr< 2, 2 >& multIntegrals,
std::array< real_t, 2 >& addIntegrals ) const
{ WALBERLA_ABORT( "Not implemented." ); }
// 3D P1
// TO IMPLEMENT
};
} // namespace hyteg
apps/RobinLaplace/helpers/boundaryForms/P1RobinLaplaceBoundaryForm.cpp deleted 100644 → 0
View file @ e7cfe321
namespace hyteg
{
class P1RobinLaplaceBoundaryForm : public P1BoundaryForm
{
public:
virtual ~P1RobinLaplaceBoundaryForm() {}
void integrateMicroEdge( const std::array< Point3D, 2 >& mircoEdgecoords,
Matrixr< 2, 2 >& multIntegrals,
std::array< real_t, 2 >& addIntegrals ) const
{
}
// private:
}
} // namespace hyteg
apps/RobinLaplace/helpers/P1CustomBoundaryOperator.cpp → src/hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.cpp
View file @ c76d9e2c
......@@ -31,6 +31,14 @@ P1CustomBoundaryOperator::P1CustomBoundaryOperator( const std::shared_ptr< Pri
, quadratureWeights_( quadratureWeights )
{}
bool P1CustomBoundaryOperator::customBoundaryCheckOfEdge( const hyteg::Point3D& endpoint0, const hyteg::Point3D& endpoint1 ) const //Temporary solution for debugging!!!
{
real_t midPointX = (endpoint0[0] + ( endpoint1[0] - endpoint0[0] ) / 2);
real_t midPointY = (endpoint0[1] + ( endpoint1[1] - endpoint0[1] ) / 2);
real_t tol = 0.05;
return (midPointX < tol) || (midPointX > 1 - tol) || (midPointY < tol) || (midPointY > 1 - tol);
}
void P1CustomBoundaryOperator::apply( const P1Function< real_t >& src,
const P1Function< real_t >& dst,
size_t level,
......@@ -53,6 +61,7 @@ void P1CustomBoundaryOperator::apply( const P1Function< real_t >& src,
const DoFType edgeBC = dst.getBoundaryCondition().getBoundaryType( edge.getMeshBoundaryFlag() );
if ( testFlag( edgeBC, flag ) )
// if( customBoundaryCheckOfEdge( edge.getCoordinates()[0], edge.getCoordinates()[1] ))
{
for( const auto& microEdgeIndices: hyteg::indexing::MicroEdgeIterator( levelinfo::num_microedges_per_edge( level ) ))
{
......@@ -106,8 +115,8 @@ void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t
integralUV += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * quadratureNodes_[i] * (1 - quadratureNodes_[i]);
integralVU += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * (1 - quadratureNodes_[i]) * quadratureNodes_[i];
integralVV += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * (1 - quadratureNodes_[i]) * (1 - quadratureNodes_[i]);
integralMicro0 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[i];
integralMicro1 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[i]);
// integralMicro0 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[i];
// integralMicro1 += - (quadratureWeights_[i] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[i]);
}
std::cout << "Boundary operator came up with the following integrals:" << std::endl;
......@@ -124,6 +133,7 @@ void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t
auto srcData = edge.getData( srcEdgeDoFIdx )->getPointer( level );
std::cout << "OLD: dstData[" << microedgeIdx[0].col() << "] = " << dstData[microedgeIdx[0].col()] << ", dstData[" << microedgeIdx[1].col() << "] = " << dstData[microedgeIdx[1].col()] << std::endl;
std::cout << "OLD: srcData[" << microedgeIdx[0].col() << "] = " << srcData[microedgeIdx[0].col()] << ", srcData[" << microedgeIdx[1].col() << "] = " << srcData[microedgeIdx[1].col()] << std::endl;
if (update == Replace)
{
......
apps/RobinLaplace/helpers/P1CustomBoundaryOperator.hpp → src/hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.hpp
View file @ c76d9e2c
......@@ -31,15 +31,17 @@ public:
std::vector< real_t > quadratureNodes,
std::vector< real_t > quadratureWeights );
virtual ~P1CustomBoundaryOperator() = default;
virtual ~P1CustomBoundaryOperator() = default;
void apply( const P1Function< real_t >& src,
bool customBoundaryCheckOfEdge( const hyteg::Point3D& endpoint0, const hyteg::Point3D& endpoint1 ) const;
void apply( const P1Function< real_t >& src,
const P1Function< real_t >& dst,
size_t level,
DoFType flag,
UpdateType updateType = Replace ) const override final;
void quadratureOnBoundaryEdge( uint_t level,
void quadratureOnBoundaryEdge( uint_t level,
const Edge& edge,
const std::array< indexing::Index, 2 > microedgeIdx,
const P1Function< real_t >& src,
......
apps/RobinLaplace/helpers/P1OperatorWithBC.cpp → src/hyteg/p1functionspace/p1boundary/P1OperatorWithBC.cpp
View file @ c76d9e2c
......@@ -12,8 +12,6 @@ P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
InnerOperatorType& innerOperator,
BoundaryOperatorType& boundaryOperator )
: Operator( storage, minLevel, maxLevel )
, solutionInner_( "solutionInner", storage, minLevel, maxLevel )
, solutionBoundary_( "solutionBoundary", storage, minLevel, maxLevel )
, innerOperator_( innerOperator )
, boundaryOperator_( boundaryOperator )
{}
......@@ -35,52 +33,72 @@ void P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
WALBERLA_ABORT( "Operator with custom boundary conditions is not available for global cells yet" );
} else
{
std::cout << "\nsrc-Function before application of P1OperatorWithBC: \n" << std::endl;
communication::syncFunctionBetweenPrimitives( src, level );
std::cout << "\n Functions before application of P1OperatorWithBC: \n" << std::endl;
std::vector< PrimitiveID > vertexIDs = storage_->getVertexIDs();
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage_->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > srcVertexDoFIdx = src.getVertexDataID();
auto srcData = vertex.getData( srcVertexDoFIdx )->getPointer( level );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > dstVertexDoFIdx = dst.getVertexDataID();
auto dstData = vertex.getData( dstVertexDoFIdx )->getPointer( level );
real_t srcAtVertex = srcData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , dst( vertex.getCoordinates() ) = " << srcAtVertex << std::endl;
real_t dstAtVertex = dstData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , src = " << srcAtVertex << " , dst = " << dstAtVertex << std::endl;
}
solutionInner_.interpolate( 0.0, maxLevel_, hyteg::All );
solutionBoundary_.interpolate( 0.0, maxLevel_, hyteg::All );
innerOperator_.apply( src, dst, level, hyteg::All, updateType );
std::cout << "\n dst-Function after application of Inneroperator: \n" << std::endl;
std::cout << "\n Functions after application of Inneroperator: \n" << std::endl;
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage_->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > dstVertexDoFIdx = dst.getVertexDataID();
auto dstData = vertex.getData( dstVertexDoFIdx )->getPointer( level );
real_t dstAtVertex = dstData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , dst at VERTEX = " << dstAtVertex << std::endl;
PrimitiveDataID< FunctionMemory< real_t >, Vertex > srcVertexDoFIdx = src.getVertexDataID();
auto srcData = vertex.getData( srcVertexDoFIdx )->getPointer( level );
real_t srcAtVertex = srcData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , src = " << srcAtVertex << " , dst = " << dstAtVertex << std::endl;
}
boundaryOperator_.apply( src, dst, level, flag, Add );
std::cout << "\n dst-Function after application of Boundaryoperator: \n" << std::endl;
std::cout << "\n Functions after application of Boundaryoperator: \n" << std::endl;
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage_->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > dstVertexDoFIdx = dst.getVertexDataID();
auto dstData = vertex.getData( dstVertexDoFIdx )->getPointer( level );
real_t dstAtVertex = dstData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , dst at VERTEX = " << dstAtVertex << std::endl;
}
std::cout << "\n" << std::endl;
PrimitiveDataID< FunctionMemory< real_t >, Vertex > srcVertexDoFIdx = src.getVertexDataID();
auto srcData = vertex.getData( srcVertexDoFIdx )->getPointer( level );
real_t srcAtVertex = srcData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , src = " << srcAtVertex << " , dst = " << dstAtVertex << std::endl;
}
/* std::cout << "\n" << std::endl;
std::vector< PrimitiveID > edgeIDs = storage_->getEdgeIDs();
for( int i = 0; i < int_c( edgeIDs.size() ); i++ ) {
Edge& edge = *storage_->getEdge( edgeIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Edge > dstEdgeDoFIdx = dst.getEdgeDataID();
auto dstData = edge.getData( dstEdgeDoFIdx )->getPointer( level );
real_t dstAtEdge0 = dstData[0];
real_t dstAtEdge1 = dstData[1];
std::cout << "edge.coord0 = " << edge.getCoordinates()[0] << " , edge.coord1 = " << edge.getCoordinates()[1] << " : dstAtEdge0 = " << dstAtEdge0 << " , dstAtEdge1 = " << dstAtEdge1 << std::endl;
}
PrimitiveDataID< FunctionMemory< real_t >, Edge > srcEdgeDoFIdx = src.getEdgeDataID();
auto srcData = edge.getData( srcEdgeDoFIdx )->getPointer( level );
// real_t dstAtEdge0 = dstData[0];
// real_t dstAtEdge1 = dstData[1];
// std::cout << "edge.coord0 = " << edge.getCoordinates()[0] << " , edge.coord1 = " << edge.getCoordinates()[1] << " : dstAtEdge0 = " << dstAtEdge0 << " , dstAtEdge1 = " << dstAtEdge1 << std::endl;
std::cout << "Functions on edge from " << edge.getCoordinates()[0] << " to " << edge.getCoordinates()[1] << " : " << std::endl;
for( int j = 0; j < levelinfo::num_microvertices_per_edge( level ); j++ )
{
std::cout << " j = " << j << " : dst = " << dstData[vertexdof::macroedge::indexFromVertex( level, j, stencilDirection::VERTEX_C )] << ", src = "
<< srcData[vertexdof::macroedge::indexFromVertex( level, j, stencilDirection::VERTEX_C )] << std::endl;
}
} */
}
this->stopTiming( "apply" );
......
apps/RobinLaplace/helpers/P1OperatorWithBC.hpp → src/hyteg/p1functionspace/p1boundary/P1OperatorWithBC.hpp
View file @ c76d9e2c
......@@ -2,6 +2,7 @@
#include "hyteg/p1functionspace/P1Operator.hpp"
#include "hyteg/operators/Operator.hpp"
#include "hyteg/communication/Syncing.hpp"
namespace hyteg {
......@@ -25,8 +26,8 @@ public:
UpdateType updateType = Replace ) const override final;
private:
P1Function< real_t > solutionInner_;
P1Function< real_t > solutionBoundary_;
// P1Function< real_t > solutionInner_;
// P1Function< real_t > solutionBoundary_;
InnerOperatorType innerOperator_;
BoundaryOperatorType boundaryOperator_;
};
......
tests/hyteg/CMakeLists.txt
View file @ c76d9e2c
......@@ -770,5 +770,8 @@ endif()
## BoundaryCommunication ##
waLBerla_compile_test(FILES vertexdofspace/P1BoundaryCommunication2DTest.cpp DEPENDS hyteg core)
waLBerla_execute_test(NAME P1BoundaryCommunication2D)
waLBerla_compile_test(FILES vertexdofspace/P1BoundaryCommunication2DTest.cpp DEPENDS hyteg core)
waLBerla_execute_test(NAME P1BoundaryCommunication2DTest)
waLBerla_compile_test(FILES operators/P1BoundaryOperatorTest.cpp DEPENDS hyteg core)
waLBerla_execute_test(NAME P1BoundaryOperatorTest)
tests/hyteg/operators/P1BoundaryOperatorTest.cpp 0 → 100644
View file @ c76d9e2c
//Add header
#include "core/mpi/all.h"
#include "core/debug/all.h"
#include "core/Environment.h"
#include "core/debug/CheckFunctions.h"
#include "core/debug/TestSubsystem.h"
#include "core/timing/all.h"
#include "hyteg/mesh/MeshInfo.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.hpp"
namespace hyteg {
using walberla::real_t;
using namespace hyteg;
void checkLinearBoundaryIntegral( std::shared_ptr< hyteg::PrimitiveStorage > storage )
{
real_t expectedResult = 6.0;
uint_t minLevel = 0;
uint_t maxLevel = 2;
real_t robinFactor = 1;
std::function< real_t( const hyteg::Point3D& ) > dummyFunction = [] ( const hyteg::Point3D& x )
{
return 0.0;
};
std::function< real_t( const hyteg::Point3D& ) > templateTestFunction = [] ( const hyteg::Point3D& x )
{
return ( 2 * x[0] + 1 );
};
hyteg::P1Function< real_t > testFunctionSRC( "testFunctionSRC", storage, minLevel, maxLevel );
testFunctionSRC.interpolate( templateTestFunction, maxLevel, hyteg::All );
hyteg::P1Function< real_t > testFunctionDST( "testFunctionDST", storage, minLevel, maxLevel);
P1CustomBoundaryOperator boundaryOp( storage, minLevel, maxLevel, robinFactor, dummyFunction );
boundaryOp.apply( testFunctionSRC, testFunctionDST, maxLevel, hyteg::NeumannBoundary );
hyteg::P1Function< real_t > biLinearFormFunction( "biLinFunction", storage, minLevel, maxLevel );
biLinearFormFunction.interpolate( 1.0, maxLevel, hyteg::All );
// real_t result = testFunctionDST.dotGlobal( biLinearFormFunction, maxLevel );
real_t result = 0.0;
std::vector< PrimitiveID > edgeIDs = storage->getEdgeIDs();
for( int i = 0; i < int_c( edgeIDs.size() ); i++ ) {
Edge& edge = *storage->getEdge( edgeIDs[uint_c( i )] );
walberla::math::KahanAccumulator< real_t > scalarProduct;
size_t rowsize = levelinfo::num_microvertices_per_edge( maxLevel );
for (size_t j = 1; j < rowsize - 1; ++j) {
scalarProduct += edge.getData(testFunctionDST.getEdgeDataID())->getPointer( maxLevel )[vertexdof::macroedge::indexFromVertex( maxLevel, j, stencilDirection::VERTEX_C )];
}
result += scalarProduct.get();
}
// WALBERLA_LOG_INFO_ON_ROOT( " [P1BoundaryOperatorTest] Computed integral = " << result );
WALBERLA_CHECK_FLOAT_EQUAL( result, expectedResult );
}
} // namespace hyteg
int main( int argc, char** argv )
{
walberla::debug::enterTestMode();
walberla::Environment env( argc, argv );
walberla::mpi::MPIManager::instance()->useWorldComm();
uint_t numProcesses = walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() );
hyteg::MeshInfo meshInfo = hyteg::MeshInfo::meshRectangle( hyteg::Point2D( {0.0, 0.0} ), hyteg::Point2D( {1.0, 1.0} ), hyteg::MeshInfo::CRISS, 2, 2 );
hyteg::SetupPrimitiveStorage setupStorage( meshInfo, numProcesses );
setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, false );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage );
checkLinearBoundaryIntegral( storage );
WALBERLA_LOG_INFO_ON_ROOT(" [P1BoundaryOperatorTest] Success! ");
return EXIT_SUCCESS;
}
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