Skip to content

GitLab

Commit 7d9dcbca authored by Maxi Dechant's avatar Maxi Dechant
Browse files

made RobinLaplace work

parent c76d9e2c
Pipeline #36153 failed with stages
in 10 minutes and 34 seconds
Hide whitespace changes
Inline Side-by-side
apps/RobinLaplace/RobinLaplace.cpp
View file @ 7d9dcbca
......@@ -24,10 +24,6 @@
#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 "hyteg/p1functionspace/p1boundary/P1OperatorWithBC.hpp"
#include "hyteg/p1functionspace/p1boundary/P1OperatorWithBC.cpp"
#include "hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.hpp"
......@@ -37,14 +33,6 @@ 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,
......@@ -63,8 +51,7 @@ 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( customBoundaryCheckOfEdge( edge.getCoordinates()[0], edge.getCoordinates()[1] ) )
if ( testFlag( edgeBC, flag ) )
{
PrimitiveDataID< FunctionMemory< real_t >, Edge > rhsEdgeDoFIdx = rhsFunction.getEdgeDataID();
auto rhsData = edge.getData( rhsEdgeDoFIdx )->getPointer( level );
......@@ -88,12 +75,14 @@ void generateRightRobinSide( const std::shared_ptr< PrimitiveStorage >&
integralMicro0 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * quadratureNodes_[j];
integralMicro1 += (quadratureWeights_[j] * microEdgeLength.norm() ) * boundaryValue * (1 - quadratureNodes_[j]);
}
rhsData[ microEdgeIndices[0].col() ] += integralMicro0;
rhsData[ microEdgeIndices[1].col() ] += integralMicro1;
rhsData[ microEdgeIndices[0].col() ] = integralMicro0;
rhsData[ microEdgeIndices[1].col() ] = integralMicro1;
// std::cout << "rhs[ " << microEdgeIndices[0].col() << " ] = " << rhsData[ microEdgeIndices[0].col() ] << " , "
// << "rhs[ " << microEdgeIndices[1].col() << " ] = " << rhsData[ microEdgeIndices[1].col() ] << std::endl;
}
}
}
rhsFunction.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage, false );
rhsFunction.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage, true );
rhsFunction.communicate< Vertex, Edge >( level );
rhsFunction.communicate< Edge, Face >( level );
}
......@@ -115,9 +104,9 @@ int main( int argc, char** argv )
const real_t robinFactor = parameters.getParameter< real_t >( "robinFactor" );
const bool useVTK = parameters.getParameter< bool >( "vtkOutput" );
hyteg::MeshInfo meshInfo = hyteg::MeshInfo::meshRectangle( hyteg::Point2D( {0.0, 0.0} ), hyteg::Point2D( {1.0, 1.0} ), hyteg::MeshInfo::CRISS, 2, 2 );
hyteg::MeshInfo meshInfo = hyteg::MeshInfo::meshRectangle( hyteg::Point2D( {0.0, 0.0} ), hyteg::Point2D( {1.0, 1.0} ), hyteg::MeshInfo::CRISS, 10, 10 );
hyteg::SetupPrimitiveStorage setupStorage( meshInfo, numProcesses );
setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, false );
setupStorage.setMeshBoundaryFlagsOnBoundary( 2, 0, true );
hyteg::loadbalancing::roundRobin( setupStorage );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage );
......@@ -134,25 +123,14 @@ int main( int argc, char** argv )
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;
real_t gradY = -x[0] - 2 * x[1] + 1.0;
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< 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 boundaryTOL = 0.0005;
real_t boundaryEvaluate = 0.0;
std::array< real_t, 2 > boundaryNormal = { 0.0, 0.0 };
if( x[0] >= 1-boundaryTOL )
......@@ -177,45 +155,22 @@ int main( int argc, char** argv )
}
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;
if( x[0] >= 1-boundaryTOL )
{
boundaryEvaluate = 2.0 - x[1] + robinFactor * (2.0 - x[1]*x[1]);
}
else if( x[1] >= 1-boundaryTOL )
{
boundaryEvaluate = -x[0] - 1.0 + robinFactor * (x[0]*x[0] - x[0] + 1.0);
}
else if( x[0] <= boundaryTOL )
{
boundaryEvaluate = x[1] + robinFactor * (-x[1]*x[1] + x[1] + 1.0);
}
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 };
hyteg::P1ConstantLaplaceOperator innerOperator( storage, minLevel, maxLevel );
P1CustomBoundaryOperator boundaryOperator( storage, minLevel, maxLevel, robinFactor, boundaryFunction, boundaryQuadNodes, boundaryQuadWeights );
P1CustomBoundaryOperator boundaryOperator( storage, minLevel, maxLevel, robinFactor, boundaryFunction, false );
P1OperatorWithBC< P1ConstantLaplaceOperator, P1CustomBoundaryOperator > CustomLaplaceOperator( storage, minLevel, maxLevel, innerOperator, boundaryOperator );
P1OperatorWithBC< P1ConstantLaplaceOperator, P1CustomBoundaryOperator > CustomLaplaceOperator( storage, minLevel, maxLevel, innerOperator, boundaryOperator, hyteg::NeumannBoundary );
hyteg::P1Function< real_t > f( "f", storage, minLevel, maxLevel );
f.interpolate( 0.0, maxLevel, hyteg::NeumannBoundary );
generateRightRobinSide( storage, f, boundaryFunction, maxLevel, hyteg::NeumannBoundary );
// f.interpolate( 42.0, maxLevel, hyteg::NeumannBoundary ); // <- PROBLEM
std::vector< PrimitiveID > vertexIDs = storage->getVertexIDs();
/* std::vector< PrimitiveID > edgeIDs = storage->getEdgeIDs();
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage->getVertex( vertexIDs[uint_c( i )] );
......@@ -224,7 +179,31 @@ int main( int argc, char** argv )
real_t fAtVertex = fData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " : f( vertex.getCoordinates() ) = " << fAtVertex << std::endl;
}
} */
generateRightRobinSide( storage, f, boundaryFunction, maxLevel, hyteg::NeumannBoundary );
/* for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
{
Vertex& vertex = *storage->getVertex( vertexIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > fVertexDoFIdx = f.getVertexDataID();
auto fData = vertex.getData( fVertexDoFIdx )->getPointer( maxLevel );
real_t fAtVertex = fData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " : f( vertex.getCoordinates() ) = " << fAtVertex << std::endl;
}
for ( int i = 0; i < int_c( edgeIDs.size() ); i++ )
{
Edge& edge = *storage->getEdge( edgeIDs[uint_c( i )] );
PrimitiveDataID< FunctionMemory< real_t >, Edge > fEdgeDoFIdx = f.getEdgeDataID();
auto fData = edge.getData( fEdgeDoFIdx )->getPointer( maxLevel );
for( uint_t j = 0; j < levelinfo::num_microvertices_per_edge( maxLevel ); j++ )
{
std::cout << "fData[ " << j << " ] = " << fData[j] << " | ";
}
std::cout << std::endl;
} */
hyteg::P1Function< real_t > u( "u", storage, minLevel, maxLevel );
u.interpolate( 0.5, maxLevel, hyteg::DirichletBoundary );
......@@ -268,7 +247,6 @@ int main( int argc, char** argv )
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( vertexIDs[uint_c( i )] );
......@@ -298,19 +276,6 @@ int main( int argc, char** argv )
" : 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() <<
" : exactOperator( vertex.getCoordinates() ) = " << exactOpAtVertex << std::endl;
} */
}
apps/RobinLaplace/RobinLaplace.prm
View file @ 7d9dcbca
......@@ -2,7 +2,7 @@ Parameters
{
minLevel 0;
maxLevel 0;
maxKrylowDim 5;
maxKrylowDim 105;
robinFactor -1;
vtkOutput true;
}
\ No newline at end of file
src/hyteg/p1functionspace/VertexDoFMacroVertex.hpp
View file @ 7d9dcbca
......@@ -167,7 +167,7 @@ inline void apply( Vertex&
auto opr_data = vertex.getData( operatorId )->getPointer( level );
auto src = vertex.getData( srcId )->getPointer( level );
auto dst = vertex.getData( dstId )->getPointer( level );
if ( update == Replace )
{
dst[0] = opr_data[0] * src[0];
......
src/hyteg/p1functionspace/p1boundary/MicroEdgeIterator.hpp
View file @ 7d9dcbca
......@@ -29,8 +29,11 @@ public:
, step_( 0 )
{
WALBERLA_ASSERT_GREATER( width, 0, "Size of edge must be larger than zero!" );
WALBERLA_ASSERT_LESS( offsetToCenter - 1, width, "Offset to center is beyond edge width!" );
if( offsetToCenter != 0 )
{
WALBERLA_ASSERT_LESS( offsetToCenter - 1, width, "Offset to center is beyond edge width!" );
}
coordinates_[0].dep() = 0;
coordinates_[0].row() = 0;
coordinates_[1].dep() = 0;
......
src/hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.cpp
View file @ 7d9dcbca
......@@ -8,10 +8,12 @@ P1CustomBoundaryOperator::P1CustomBoundaryOperator( const std::shared_ptr< Pri
size_t minLevel,
size_t maxLevel,
real_t robinFactor,
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction )
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
bool replaceOnVertices )
: Operator( storage, minLevel, maxLevel )
, robinFactor_( robinFactor )
, boundaryFunction_( boundaryFunction )
, replaceOnVertices_( replaceOnVertices )
{
quadratureNodes_.push_back( 0.5 );
quadratureWeights_.push_back( 1.0 );
......@@ -22,23 +24,17 @@ P1CustomBoundaryOperator::P1CustomBoundaryOperator( const std::shared_ptr< Pri
size_t maxLevel,
real_t robinFactor,
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
bool replaceOnVertices,
std::vector< real_t > quadratureNodes,
std::vector< real_t > quadratureWeights )
: Operator( storage, minLevel, maxLevel )
, robinFactor_( robinFactor )
, boundaryFunction_( boundaryFunction )
, replaceOnVertices_( replaceOnVertices )
, quadratureNodes_( quadratureNodes )
, 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,
......@@ -51,7 +47,6 @@ void P1CustomBoundaryOperator::apply( const P1Function< real_t >& src,
}
else
{
std::cout << "\n\nBoundary operator was called..." << std::endl;
communication::syncFunctionBetweenPrimitives( src, level );
std::vector< hyteg::PrimitiveID > edgeIDs = this->getStorage()->getEdgeIDs();
......@@ -61,7 +56,6 @@ 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 ) ))
{
......@@ -69,9 +63,7 @@ void P1CustomBoundaryOperator::apply( const P1Function< real_t >& src,
}
}
}
dst.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage_, false );
dst.communicate< Vertex, Edge >( level );
dst.communicate< Edge, Face >( level );
dst.communicateAdditively< Edge, Vertex >( level, hyteg::All ^ hyteg::NeumannBoundary, *storage_, replaceOnVertices_ );
}
}
......@@ -82,13 +74,13 @@ void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t
const P1Function< real_t >& dst,
UpdateType update ) const
{
std::cout << "\nquadratureOnBoundaryEdge() was called..." << std::endl;
// std::cout << "\nquadratureOnBoundaryEdge() was called..." << std::endl;
real_t edgeNum = real_t( levelinfo::num_microedges_per_edge( level ) );
Point3D microEdgeLength = (edge.getCoordinates()[1] - edge.getCoordinates()[0]) / edgeNum;
Point3D coord0 = edge.getCoordinates()[0] + real_t( microedgeIdx[0].col() ) * microEdgeLength;
Point3D coord1 = edge.getCoordinates()[0] + real_t( microedgeIdx[1].col() ) * microEdgeLength;
/* Point3D coord0 = edge.getCoordinates()[0] + real_t( microedgeIdx[0].col() ) * microEdgeLength;
Point3D coord1 = edge.getCoordinates()[0] + real_t( microedgeIdx[1].col() ) * microEdgeLength; */
std::cout << "Quadrature on boundary with parameters..." << std::endl;
/* std::cout << "Quadrature on boundary with parameters..." << std::endl;
std::cout << " edge.getCoordinates()[0] = " << edge.getCoordinates()[0] << std::endl;
std::cout << " edge.getCoordinates()[1] = " << edge.getCoordinates()[1] << std::endl;
std::cout << " microedgeIdx[0].col() = " << microedgeIdx[0].col() << std::endl;
......@@ -96,44 +88,38 @@ void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t
std::cout << " edgeNum = " << edgeNum << std::endl;
std::cout << " coord0 = " << coord0 << std::endl;
std::cout << " coord1 = " << coord1 << std::endl;
std::cout << " microEdgeLength = " << microEdgeLength << std::endl;
std::cout << " microEdgeLength = " << microEdgeLength << std::endl; */
real_t integralUU = 0;
real_t integralUV = 0;
real_t integralVU = 0;
real_t integralVV = 0;
real_t integralMicro0 = 0;
real_t integralMicro1 = 0;
for(uint_t i = 0; i < quadratureNodes_.size(); i++)
{
Point3D currentNodeCoord({ coord0[0] * quadratureNodes_[i] + coord1[0] * (1 - quadratureNodes_[i]),
/* Point3D currentNodeCoord({ coord0[0] * quadratureNodes_[i] + coord1[0] * (1 - quadratureNodes_[i]),
coord0[1] * quadratureNodes_[i] + coord1[1] * (1 - quadratureNodes_[i])});
real_t boundaryValue = boundaryFunction_( currentNodeCoord );
real_t boundaryValue = boundaryFunction_( currentNodeCoord ); */
integralUU += quadratureWeights_[i] * microEdgeLength.norm() * robinFactor_ * quadratureNodes_[i] * quadratureNodes_[i];
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]);
}
std::cout << "Boundary operator came up with the following integrals:" << std::endl;
/* std::cout << "Boundary operator came up with the following integrals:" << std::endl;
std::cout << " integralUU = " << integralUU << std::endl;
std::cout << " integralUV = " << integralUV << std::endl;
std::cout << " integralVU = " << integralVU << std::endl;
std::cout << " integralVV = " << integralVV << std::endl;
std::cout << " integralMicro0 = " << integralMicro0 << std::endl;
std::cout << " integralMicro1 = " << integralMicro1 << std::endl;
std::cout << " integralVV = " << integralVV << std::endl; */
PrimitiveDataID< FunctionMemory< real_t >, Edge > dstEdgeDoFIdx = dst.getEdgeDataID();
PrimitiveDataID< FunctionMemory< real_t >, Edge > srcEdgeDoFIdx = src.getEdgeDataID();
auto dstData = edge.getData( dstEdgeDoFIdx )->getPointer( level );
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;
// 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)
{
......@@ -150,9 +136,7 @@ void P1CustomBoundaryOperator::quadratureOnBoundaryEdge( uint_t
dstData[ microedgeIdx[0].col() ] += integralUU * srcData[ microedgeIdx[0].col() ] + integralUV * srcData[ microedgeIdx[1].col() ];
dstData[ microedgeIdx[1].col() ] += integralVU * srcData[ microedgeIdx[0].col() ] + integralVV * srcData[ microedgeIdx[1].col() ];
}
dstData[ microedgeIdx[0].col() ] += integralMicro0;
dstData[ microedgeIdx[1].col() ] += integralMicro1;
std::cout << "NEW: dstData[" << microedgeIdx[0].col() << "] = " << dstData[microedgeIdx[0].col()] << ", dstData[" << microedgeIdx[1].col() << "] = " << dstData[microedgeIdx[1].col()] << std::endl;
// std::cout << "NEW: dstData[" << microedgeIdx[0].col() << "] = " << dstData[microedgeIdx[0].col()] << ", dstData[" << microedgeIdx[1].col() << "] = " << dstData[microedgeIdx[1].col()] << std::endl;
}
......
src/hyteg/p1functionspace/p1boundary/P1CustomBoundaryOperator.hpp
View file @ 7d9dcbca
......@@ -21,27 +21,29 @@ public:
size_t minLevel,
size_t maxLevel,
real_t robinFactor,
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction );
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction,
bool replaceOnVertices );
P1CustomBoundaryOperator( const std::shared_ptr< hyteg::PrimitiveStorage >& storage,
size_t minLevel,
size_t maxLevel,
real_t robinFactor,
std::function< real_t( const hyteg::Point3D& )>& boundaryFunction,
bool replaceOnVertices,
std::vector< real_t > quadratureNodes,
std::vector< real_t > quadratureWeights );
virtual ~P1CustomBoundaryOperator() = default;
bool customBoundaryCheckOfEdge( const hyteg::Point3D& endpoint0, const hyteg::Point3D& endpoint1 ) const;
// 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;
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,
......@@ -51,8 +53,9 @@ public:
private:
real_t robinFactor_;
std::function< real_t( const hyteg::Point3D& ) >& boundaryFunction_;
std::vector< real_t > quadratureNodes_;
std::vector< real_t > quadratureWeights_;
bool replaceOnVertices_;
std::vector< real_t > quadratureNodes_;
std::vector< real_t > quadratureWeights_;
};
}
\ No newline at end of file
src/hyteg/p1functionspace/p1boundary/P1OperatorWithBC.cpp
View file @ 7d9dcbca
......@@ -10,10 +10,14 @@ P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
size_t minLevel,
size_t maxLevel,
InnerOperatorType& innerOperator,
BoundaryOperatorType& boundaryOperator )
BoundaryOperatorType& boundaryOperator,
hyteg::DoFType boundaryDoFs )
: Operator( storage, minLevel, maxLevel )
, dstInner_( "dstInner", storage, minLevel, maxLevel )
, dstBoundary_( "dstBoundary", storage, minLevel, maxLevel )
, innerOperator_( innerOperator )
, boundaryOperator_( boundaryOperator )
, boundaryDoFs_( boundaryDoFs )
{}
template< class InnerOperatorType, class BoundaryOperatorType >
......@@ -35,7 +39,7 @@ void P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
{
communication::syncFunctionBetweenPrimitives( src, level );
std::cout << "\n Functions before application of P1OperatorWithBC: \n" << std::endl;
/* 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++ )
{
......@@ -47,11 +51,12 @@ void P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
real_t srcAtVertex = srcData[0];
real_t dstAtVertex = dstData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , src = " << srcAtVertex << " , dst = " << dstAtVertex << std::endl;
}
} */
innerOperator_.apply( src, dst, level, hyteg::All, updateType );
dstInner_.interpolate( 0.0, maxLevel_, hyteg::All );
innerOperator_.apply( src, dstInner_, level, flag ); // DEBUG
std::cout << "\n Functions 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 )] );
......@@ -63,9 +68,28 @@ void P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
real_t srcAtVertex = srcData[0];
std::cout << "vertex.getCoordinates() = " << vertex.getCoordinates() << " , src = " << srcAtVertex << " , dst = " << dstAtVertex << std::endl;
} */
dstBoundary_.interpolate( 0.0, maxLevel_, hyteg::All );
boundaryOperator_.apply( src, dstBoundary_, level, boundaryDoFs_, Replace ); // DEBUG
if( updateType == Replace ) {
dst.assign( {1.0}, {dstInner_}, level, hyteg::All );
} else {
dst.assign( {1.0, 1.0}, {dst, dstInner_}, level, hyteg::All );
}
dst.assign( {1.0, 1.0}, {dst, dstBoundary_}, level, boundaryDoFs_ );
boundaryOperator_.apply( src, dst, level, flag, Add );
/* std::cout << "\n inner- and boundary-dst after application of operator: \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 > innerDstVertexDoFIdx = dstInner_.getVertexDataID();
auto innerDstData = vertex.getData( innerDstVertexDoFIdx )->getPointer( level );
PrimitiveDataID< FunctionMemory< real_t >, Vertex > boundaryDstVertexDoFIdx = dstBoundary_.getVertexDataID();
auto boundaryDstData = vertex.getData( boundaryDstVertexDoFIdx )->getPointer( level );
std::cout << "vertex.getID() = " << vertex.getID() << " , innerDst = " << innerDstData[0] << " , boundaryDst = " << boundaryDstData[0] << std::endl;
}
std::cout << "\n Functions after application of Boundaryoperator: \n" << std::endl;
for ( int i = 0; i < int_c( vertexIDs.size() ); i++ )
......@@ -79,7 +103,7 @@ void P1OperatorWithBC< InnerOperatorType, BoundaryOperatorType >::
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++ ) {
......
src/hyteg/p1functionspace/p1boundary/P1OperatorWithBC.hpp
View file @ 7d9dcbca
......@@ -16,7 +16,8 @@ public:
size_t minLevel,
size_t maxLevel,
InnerOperatorType& innerOperator,
BoundaryOperatorType& boundaryOperator );
BoundaryOperatorType& boundaryOperator,
hyteg::DoFType boundaryDoFs_ = hyteg::RobinBoundary );
virtual ~P1OperatorWithBC() = default;
void apply( const P1Function< real_t >& src,
......@@ -26,10 +27,11 @@ public:
UpdateType updateType = Replace ) const override final;
private:
// P1Function< real_t > solutionInner_;
// P1Function< real_t > solutionBoundary_;
InnerOperatorType innerOperator_;
BoundaryOperatorType boundaryOperator_;
P1Function< real_t > dstInner_;
P1Function< real_t > dstBoundary_;
InnerOperatorType innerOperator_;
BoundaryOperatorType boundaryOperator_;
hyteg::DoFType boundaryDoFs_;
};
}
......
src/hyteg/solvers/preconditioners/IdentityPreconditioner.hpp
View file @ 7d9dcbca
......@@ -33,7 +33,7 @@ class IdentityPreconditioner : public Solver< OperatorType >
public:
IdentityPreconditioner()
: updateType_( Replace )
, flag_( hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary )
, flag_( hyteg::Inner | hyteg::NeumannBoundary | hyteg::FreeslipBoundary | hyteg::RobinBoundary )
{}
void solve( const OperatorType&,
......
tests/hyteg/operators/P1BoundaryOperatorTest.cpp
View file @ 7d9dcbca
......@@ -21,7 +21,7 @@ using namespace hyteg;
void checkLinearBoundaryIntegral( std::shared_ptr< hyteg::PrimitiveStorage > storage )
{
real_t expectedResult = 6.0;
real_t expectedResult = 8.0;
uint_t minLevel = 0;
uint_t maxLevel = 2;
......@@ -39,28 +39,15 @@ void checkLinearBoundaryIntegral( std::shared_ptr< hyteg::PrimitiveStorage > sto
hyteg::P1Function< real_t > testFunctionSRC( "testFunctionSRC", storage, minLevel, maxLevel );
testFunctionSRC.interpolate( templateTestFunction, maxLevel, hyteg::All );
hyteg::P1Function< real_t > testFunctionDST( "testFunctionDST", storage, minLevel, maxLevel);
testFunctionDST.interpolate( 69.0, maxLevel, hyteg::All );
P1CustomBoundaryOperator boundaryOp( storage, minLevel, maxLevel, robinFactor, dummyFunction );
P1CustomBoundaryOperator boundaryOp( storage, minLevel, maxLevel, robinFactor, dummyFunction, true );
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 )[