Skip to content

GitLab

Commit ab4ed4c0 authored by Marcus Mohr's avatar Marcus Mohr
Browse files

Merge branch 'kohl/toMatrix' into 'master' 

Adds a PetscMatrixAssemblyTest

See merge request !450
parents 51a991b9 1d192b38
Pipeline #34637 canceled with stages
in 27 seconds
Hide whitespace changes
Inline Side-by-side
apps/2020-moc/MantleConvection.cpp
View file @ ab4ed4c0
......@@ -32,6 +32,7 @@
#include "hyteg/dataexport/SQL.hpp"
#include "hyteg/dataexport/TimingOutput.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/elementwiseoperators/P2P1ElementwiseBlendingStokesOperator.hpp"
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/geometry/AnnulusMap.hpp"
#include "hyteg/geometry/IcosahedralShellMap.hpp"
......
apps/2021-tokamak/Tokamak.cpp
View file @ ab4ed4c0
......@@ -551,7 +551,7 @@ void tokamak( TokamakDomain tokamakDomain,
{
const auto relativeResidualToleranceCoarseGrid = 1e-30;
const auto absoluteResidualToleranceCoarseGrid = 1e-12;
const auto maxIterationsCoarseGrid = static_cast< PetscInt >( solverSettings.maxCoarseGridSolverIterations );
const auto maxIterationsCoarseGrid = static_cast< idx_t >( solverSettings.maxCoarseGridSolverIterations );
auto actualCoarseGridSolver =
std::make_shared< PETScCGSolver< LaplaceOperator_T > >( storage,
minLevel,
......
apps/StokesFlowSolverComparison.cpp
View file @ ab4ed4c0
......@@ -100,7 +100,7 @@ public:
storage_( storage ), velocityUBC_( velocityUBC ), velocityVBC_( velocityVBC )
{
tmpRHS_ = std::make_shared< Function_T< real_t > >( "tmpRHS", storage, minLevel, maxLevel );
numerator_ = std::make_shared< Function_T< PetscInt > >( "numerator", storage, minLevel, maxLevel );
numerator_ = std::make_shared< Function_T< idx_t > >( "numerator", storage, minLevel, maxLevel );
}
void solve( const Operator_T & A,
......@@ -119,11 +119,11 @@ public:
}
private:
std::shared_ptr< Function_T< PetscInt > > numerator_;
std::shared_ptr< Function_T< real_t > > tmpRHS_;
std::shared_ptr< PrimitiveStorage > storage_;
std::function< real_t ( const hyteg::Point3D & ) > velocityUBC_;
std::function< real_t ( const hyteg::Point3D & ) > velocityVBC_;
std::shared_ptr< Function_T< idx_t > > numerator_;
std::shared_ptr< Function_T< real_t > > tmpRHS_;
std::shared_ptr< PrimitiveStorage > storage_;
std::function< real_t( const hyteg::Point3D& ) > velocityUBC_;
std::function< real_t( const hyteg::Point3D& ) > velocityVBC_;
#else
public:
PetscSolver( const std::shared_ptr< hyteg::PrimitiveStorage > &,
......
apps/annulus.cpp
View file @ ab4ed4c0
......@@ -23,7 +23,6 @@
#include <core/timing/Timer.h>
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
#include "hyteg/geometry/AnnulusMap.hpp"
......@@ -84,20 +83,19 @@ void solveProblem( std::shared_ptr< hyteg::PrimitiveStorage >& storage, uint_t l
opType lapOp( storage, level, level );
// determine indices and dimensions
funcType< PetscInt > enumerator( "enumerator", storage, level, level );
funcType< idx_t > enumerator( "enumerator", storage, level, level );
enumerator.enumerate( level );
typedef typename FunctionTrait< funcType< PetscInt > >::Tag enumTag;
typedef typename FunctionTrait< funcType< idx_t > >::Tag enumTag;
uint_t globalDoFs = numberOfGlobalDoFs< enumTag >( *storage, level );
uint_t localDoFs = numberOfLocalDoFs< enumTag >( *storage, level );
// assemble matrices
PETScSparseMatrix< opType, funcType > lapMat( localDoFs, globalDoFs );
PETScSparseMatrix< opType > lapMat( localDoFs, globalDoFs );
switch ( verbosity )
{
case 2:
{
case 2: {
lapMat.createMatrixFromOperator( lapOp, level, enumerator, All );
lapMat.applyDirichletBC( enumerator, level );
......@@ -117,17 +115,16 @@ void solveProblem( std::shared_ptr< hyteg::PrimitiveStorage >& storage, uint_t l
// << globalDoFCount );
// break;
}
[[fallthrough]];
case 1:
{
[[fallthrough]];
case 1: {
WALBERLA_LOG_INFO_ON_ROOT( "* no. of global DoFs (HyTeG) ............. " << globalDoFs );
WALBERLA_LOG_INFO_ON_ROOT( "* no. of local DoFs (HyTeG) .............. " << localDoFs );
uint_t globalDoFCount = numberOfGlobalInnerDoFs< enumTag >( *storage, level );
WALBERLA_LOG_INFO_ON_ROOT( "* no. of global inner DoFs (HyTeG) ....... " << globalDoFCount );
break;
}
default:
{}
default: {
}
}
WALBERLA_LOG_INFO_ON_ROOT( " Entering PETSc solve " );
......@@ -138,7 +135,7 @@ void solveProblem( std::shared_ptr< hyteg::PrimitiveStorage >& storage, uint_t l
// check size of error
funcType< real_t > error( "error", storage, level, level );
error.assign( {1.0, -1.0}, {u_exact, u}, level, All );
error.assign( { 1.0, -1.0 }, { u_exact, u }, level, All );
VTKOutput vtkOutput( "../output", fileName, storage );
// VTKOutput vtkOutput( "../output", "annulus", storage );
......
apps/benchmarks/ElementwiseOps/PolarLaplacianBenchmark.cpp
View file @ ab4ed4c0
......@@ -22,7 +22,6 @@
#include <core/timing/Timer.h>
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
#include "hyteg/forms/form_fenics_base/P1FenicsForm.hpp"
......@@ -180,7 +179,7 @@ caseResult analyseCase( std::shared_ptr< PrimitiveStorage > storage,
funcType& u,
funcType& error )
{
typedef typename FunctionTrait< typename funcType::template FunctionType< PetscInt > >::Tag funcTag;
typedef typename FunctionTrait< typename funcType::template FunctionType< idx_t > >::Tag funcTag;
// embed numeric solution in next finer space
embedInRefinedSpace( u, level );
......@@ -379,7 +378,7 @@ void solve_using_pimped_form( uint_t minLevel, uint_t maxLevel, bool outputVTK )
{
// perform some template magic
typedef typename opType::srcType funcType;
// typedef typename FunctionTrait< typename opType::srcType::template FunctionType<PetscInt> >::Tag funcTag;
// typedef typename FunctionTrait< typename opType::srcType::template FunctionType<idx_t> >::Tag funcTag;
// generate annulus mesh in polar coordinates
real_t rmin = 1.0;
......
apps/benchmarks/ElementwiseOps/StencilScalingBenchmark.cpp
View file @ ab4ed4c0
......@@ -29,7 +29,6 @@
#include <core/timing/Timer.h>
#include "hyteg/dataexport/VTKOutput.hpp"
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
#include "hyteg/gridtransferoperators/P2toP2QuadraticProlongation.hpp"
#include "hyteg/petsc/PETScLUSolver.hpp"
......@@ -63,7 +62,7 @@ caseResult analyseCase( std::shared_ptr< PrimitiveStorage > storage,
const funcType& u,
funcType& error )
{
typedef typename FunctionTrait< typename funcType::template FunctionType< PetscInt > >::Tag funcTag;
typedef typename FunctionTrait< typename funcType::template FunctionType< idx_t > >::Tag funcTag;
// embed numeric solution in next finer space
P2toP2QuadraticProlongation embeddor;
......
apps/benchmarks/PetscCompare-2D-P2-Apply/PetscCompare-2D-P2-Apply.cpp
View file @ ab4ed4c0
......@@ -111,7 +111,7 @@ int main( int argc, char* argv[] )
hyteg::P2Function< double > y( "y", storage, level, level );
hyteg::P2Function< double > z( "z", storage, level, level );
hyteg::P2Function< double > diff( "diff", storage, level, level );
hyteg::P2Function< PetscInt > numerator( "numerator", storage, level, level );
hyteg::P2Function< idx_t > numerator( "numerator", storage, level, level );
wcTimingTreeApp.stop( "Function allocation" );
const uint_t totalDoFs = numberOfGlobalDoFs< hyteg::P2FunctionTag >( *storage, level );
......@@ -133,7 +133,7 @@ int main( int argc, char* argv[] )
LIKWID_MARKER_START( "PETSc-setup" );
wcTimingTreeApp.start( "Petsc setup" );
hyteg::PETScSparseMatrix< hyteg::P2ConstantLaplaceOperator, hyteg::P2Function > matPetsc( localDoFs, globalDoFs );
hyteg::PETScSparseMatrix< hyteg::P2ConstantLaplaceOperator > matPetsc( localDoFs, globalDoFs );
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P2Function > vecPetsc( localDoFs );
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
......
apps/benchmarks/PetscCompare-3D-P1-Apply/PetscCompare-3D-P1-Apply.cpp
View file @ ab4ed4c0
......@@ -21,9 +21,9 @@
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/config/Config.h"
#include "core/math/Constants.h"
#include "core/mpi/MPIManager.h"
#include "core/timing/TimingJSON.h"
#include "core/math/Constants.h"
#include "hyteg/LikwidWrapper.hpp"
#include "hyteg/dataexport/VTKOutput.hpp"
......@@ -55,12 +55,13 @@ int main( int argc, char* argv[] )
PETScManager petscManager( &argc, &argv );
auto cfg = std::make_shared< walberla::config::Config >();
if( env.config() == nullptr )
if ( env.config() == nullptr )
{
auto defaultFile = "./PetscCompare-3D-P1-Apply.prm";
WALBERLA_LOG_PROGRESS_ON_ROOT( "No Parameter file given loading default parameter file: " << defaultFile );
cfg->readParameterFile( defaultFile );
} else
}
else
{
cfg = env.config();
}
......@@ -70,14 +71,15 @@ int main( int argc, char* argv[] )
wcTimingTreeApp.start( "Mesh setup + load balancing" );
std::shared_ptr< hyteg::MeshInfo > meshInfo;
if( mainConf.getParameter< bool >( "useMeshFile" ) )
if ( mainConf.getParameter< bool >( "useMeshFile" ) )
{
std::string meshFileName = mainConf.getParameter< std::string >( "mesh" );
meshInfo = std::make_shared< hyteg::MeshInfo >( hyteg::MeshInfo::fromGmshFile( meshFileName ) );
} else
}
else
{
uint_t numberOfFaces = mainConf.getParameter< uint_t >( "numberOfFaces" );
if( mainConf.getParameter< bool >( "facesTimesProcs" ) )
if ( mainConf.getParameter< bool >( "facesTimesProcs" ) )
{
meshInfo = std::make_shared< hyteg::MeshInfo >(
hyteg::MeshInfo::meshFaceChain( numberOfFaces * uint_c( walberla::MPIManager::instance()->numProcesses() ) ) );
......@@ -85,13 +87,13 @@ int main( int argc, char* argv[] )
}
hyteg::SetupPrimitiveStorage setupStorage( *meshInfo,
walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
walberla::uint_c( walberla::mpi::MPIManager::instance()->numProcesses() ) );
uint_t numberOfFaces = setupStorage.getNumberOfFaces();
hyteg::loadbalancing::roundRobin( setupStorage );
std::shared_ptr< walberla::WcTimingTree > timingTree( new walberla::WcTimingTree() );
std::shared_ptr< walberla::WcTimingTree > timingTree( new walberla::WcTimingTree() );
std::shared_ptr< hyteg::PrimitiveStorage > storage = std::make_shared< hyteg::PrimitiveStorage >( setupStorage, timingTree );
wcTimingTreeApp.stop( "Mesh setup + load balancing" );
......@@ -108,7 +110,7 @@ int main( int argc, char* argv[] )
hyteg::P1Function< double > y( "y", storage, level, level );
hyteg::P1Function< double > z( "z", storage, level, level );
hyteg::P1Function< double > diff( "diff", storage, level, level );
hyteg::P1Function< PetscInt > numerator( "numerator", storage, level, level );
hyteg::P1Function< idx_t > numerator( "numerator", storage, level, level );
wcTimingTreeApp.stop( "Function allocation" );
const uint_t totalDoFs = numberOfGlobalDoFs< hyteg::P1FunctionTag >( *storage, level );
......@@ -130,7 +132,7 @@ int main( int argc, char* argv[] )
LIKWID_MARKER_START( "PETSc-setup" );
wcTimingTreeApp.start( "Petsc setup" );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator, hyteg::P1Function > matPetsc( localDoFs, globalDoFs );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc( localDoFs, globalDoFs );
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc( localDoFs );
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
......@@ -162,7 +164,7 @@ int main( int argc, char* argv[] )
//dstvecPetsc.print("../output/vector1.vec");
if( mainConf.getParameter< bool >( "printTiming" ) )
if ( mainConf.getParameter< bool >( "printTiming" ) )
{
auto wcTPReduced = wcTimingTreeApp.getReduced();
WALBERLA_LOG_INFO_ON_ROOT( wcTPReduced );
......@@ -179,12 +181,12 @@ int main( int argc, char* argv[] )
jsonOutput.close();
}
diff.assign( {1.0, -1.0}, {z, y}, level, hyteg::All );
diff.assign( { 1.0, -1.0 }, { z, y }, level, hyteg::All );
if( mainConf.getParameter< bool >( "VTKOutput" ) )
if ( mainConf.getParameter< bool >( "VTKOutput" ) )
{
WALBERLA_LOG_INFO_ON_ROOT( "writing VTK output" );
hyteg::VTKOutput vtkOutput("./output", "PetscCompare-2D-P2-Apply", storage);
hyteg::VTKOutput vtkOutput( "./output", "PetscCompare-2D-P2-Apply", storage );
vtkOutput.add( x );
vtkOutput.add( z );
vtkOutput.add( y );
......
apps/benchmarks/PetscCompare/PetscCompare.cpp
View file @ ab4ed4c0
......@@ -114,7 +114,7 @@ int main( int argc, char* argv[] )
hyteg::P1Function< double > diff( "diff", storage, level, level );
x.interpolate( exact, level, hyteg::Inner );
//hyteg::communication::syncFunctionBetweenPrimitives(x,level);
hyteg::P1Function< PetscInt > numerator( "numerator", storage, level, level );
hyteg::P1Function< idx_t > numerator( "numerator", storage, level, level );
hyteg::P1ConstantLaplaceOperator mass( storage, level, level );
// for (const auto & faceIT : storage->getFaces()) {
......@@ -155,14 +155,14 @@ int main( int argc, char* argv[] )
// WALBERLA_CRITICAL_SECTION_START
// for (auto &edgeIT : storage->getEdges()) {
// auto edge = edgeIT.second;
// hyteg::vertexdof::macroedge::printFunctionMemory< PetscInt >(level, *edge, numerator.getEdgeDataID());
// hyteg::vertexdof::macroedge::printFunctionMemory< idx_t >(level, *edge, numerator.getEdgeDataID());
// }
// WALBERLA_CRITICAL_SECTION_END
numerator.enumerate( level );
LIKWID_MARKER_START( "PETSc-setup" );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator, hyteg::P1Function > matPetsc( localDoFs, totalDoFs );
hyteg::PETScSparseMatrix< hyteg::P1ConstantLaplaceOperator > matPetsc( localDoFs, totalDoFs );
matPetsc.createMatrixFromOperator( mass, level, numerator, hyteg::Inner );
hyteg::PETScVector< real_t, hyteg::P1Function > vecPetsc( localDoFs );
vecPetsc.createVectorFromFunction( x, numerator, level, hyteg::Inner );
......
apps/stokesSphere/StokesSphere.cpp
View file @ ab4ed4c0
......@@ -281,7 +281,7 @@ int main( int argc, char* argv[] )
}
#if 0
auto numerator = std::make_shared< hyteg::P1StokesFunction< PetscInt > >( "numerator", storage, level, level );
auto numerator = std::make_shared< hyteg::P1StokesFunction< idx_t > >( "numerator", storage, level, level );
uint_t globalSize = 0;
const uint_t localSize = numerator->enumerate(level, globalSize);
PETScManager petscManager( &argc, &argv );
......
src/hyteg/boundary/BoundaryConditions.hpp
View file @ ab4ed4c0
......@@ -21,7 +21,8 @@
#include "core/DataTypes.h"
#include "core/uid/all.h"
#include "hyteg/types/flags.hpp"
#include "hyteg/types/types.hpp"
namespace hyteg {
......
src/hyteg/communication/DoFSpacePackInfo.cpp
View file @ ab4ed4c0
......@@ -58,6 +58,7 @@ DoFSpacePackInfo< ValueType >::DoFSpacePackInfo( uint_t
template class DoFSpacePackInfo< double >;
template class DoFSpacePackInfo< int >;
template class DoFSpacePackInfo< long >;
template class DoFSpacePackInfo< long long >;
template class DoFSpacePackInfo< uint_t >;
} // namespace communication
......
src/hyteg/composites/ConcatenatedOperator.hpp
View file @ ab4ed4c0
......@@ -109,7 +109,7 @@ class ConcatenatedOperator : public Operator< typename OpType1::srcType, typenam
op1_->computeDiagonalOperatorValues();
op2_->computeDiagonalOperatorValues();
for ( uint_t level = op1_->getMinLevel(); level <= op1_->getMaxLevel(); level += 1 )
diagonalValues_->multElementwise( { *op1_->getDiagonalValues(), *op2_->getDiagonalValues() }, level, All );
diagonalValues_->multElementwise( {*op1_->getDiagonalValues(), *op2_->getDiagonalValues()}, level, All );
}
/// Trigger (re)computation of inverse diagonal matrix entries (central operator weights)
......@@ -123,10 +123,10 @@ class ConcatenatedOperator : public Operator< typename OpType1::srcType, typenam
op2_->computeDiagonalOperatorValues();
for ( uint_t level = op1_->getMinLevel(); level <= op1_->getMaxLevel(); level += 1 )
{
inverseDiagonalValues_->multElementwise( { *op1_->getDiagonalValues(), *op2_->getDiagonalValues() }, level, All );
inverseDiagonalValues_->multElementwise( {*op1_->getDiagonalValues(), *op2_->getDiagonalValues()}, level, All );
// invert:
inverseDiagonalValues_->interpolate(
[]( auto, auto val ) { return 1. / val[0]; }, { *inverseDiagonalValues_, *inverseDiagonalValues_ }, level, All );
[]( auto, auto val ) { return 1. / val[0]; }, {*inverseDiagonalValues_, *inverseDiagonalValues_}, level, All );
}
}
......
src/hyteg/composites/P1BlendingStokesOperator.hpp
View file @ ab4ed4c0
......@@ -19,11 +19,11 @@
*/
#pragma once
#include "hyteg/composites/P1StokesBlockPreconditioner.hpp"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/StokesOperatorTraits.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1VariableOperator.hpp"
#include "hyteg/composites/P1StokesBlockPreconditioner.hpp"
namespace hyteg {
......@@ -50,6 +50,8 @@ class P1BlendingStokesOperator : public Operator< P1StokesFunction< real_t >, P1
const uint_t level,
const DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1BlendingStokesOperator not implemented for 3D." );
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A_uu.apply( src.uvw[0], dst.uvw[0], level, flag, Replace );
......@@ -65,6 +67,28 @@ class P1BlendingStokesOperator : public Operator< P1StokesFunction< real_t >, P1
pspg.apply( src.p, dst.p, level, flag | DirichletBoundary, Add );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< idx_t >& src,
const P1StokesFunction< idx_t >& dst,
size_t level,
DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1BlendingStokesOperator not implemented for 3D." );
A_uu.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
A_uv.toMatrix( mat, src.uvw[1], dst.uvw[0], level, flag );
divT_x.toMatrix( mat, src.p, dst.uvw[0], level, flag );
A_vu.toMatrix( mat, src.uvw[0], dst.uvw[1], level, flag );
A_vv.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
divT_y.toMatrix( mat, src.p, dst.uvw[1], level, flag );
div_x.toMatrix( mat, src.uvw[0], dst.p, level, flag | DirichletBoundary );
div_y.toMatrix( mat, src.uvw[1], dst.p, level, flag | DirichletBoundary );
pspg.toMatrix( mat, src.p, dst.p, level, flag | DirichletBoundary );
}
P1BlendingEpsilonOperator_11 A_uu;
P1BlendingEpsilonOperator_12 A_uv;
P1BlendingEpsilonOperator_21 A_vu;
......
src/hyteg/composites/P1EpsilonStokesOperator.hpp
View file @ ab4ed4c0
......@@ -43,6 +43,7 @@ class P1EpsilonStokesOperator : public Operator< P1StokesFunction< real_t >, P1S
void apply( const P1StokesFunction< real_t >& src, const P1StokesFunction< real_t >& dst, size_t level, DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1EpsilonStokesOperator not implemented for 3D." );
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A_uu.apply( src.uvw[0], dst.uvw[0], level, flag, Replace );
......@@ -58,6 +59,28 @@ class P1EpsilonStokesOperator : public Operator< P1StokesFunction< real_t >, P1S
pspg.apply( src.p, dst.p, level, flag | DirichletBoundary, Add );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< idx_t >& src,
const P1StokesFunction< idx_t >& dst,
size_t level,
DoFType flag ) const
{
WALBERLA_CHECK( !src.uvw[0].getStorage()->hasGlobalCells(), "P1EpsilonStokesOperator not implemented for 3D." );
A_uu.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
A_uv.toMatrix( mat, src.uvw[1], dst.uvw[0], level, flag );
divT_x.toMatrix( mat, src.p, dst.uvw[0], level, flag );
A_vu.toMatrix( mat, src.uvw[0], dst.uvw[1], level, flag );
A_vv.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
divT_y.toMatrix( mat, src.p, dst.uvw[1], level, flag );
div_x.toMatrix( mat, src.uvw[0], dst.p, level, flag | DirichletBoundary );
div_y.toMatrix( mat, src.uvw[1], dst.p, level, flag | DirichletBoundary );
pspg.toMatrix( mat, src.p, dst.p, level, flag | DirichletBoundary );
}
P1ConstantEpsilonOperator_11 A_uu;
P1ConstantEpsilonOperator_12 A_uv;
P1ConstantEpsilonOperator_21 A_vu;
......
src/hyteg/composites/P1StokesBlockLaplaceOperator.hpp
View file @ ab4ed4c0
......@@ -44,6 +44,21 @@ class P1StokesBlockLaplaceOperator : public Operator< P1StokesFunction< real_t >
}
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< idx_t >& src,
const P1StokesFunction< idx_t >& dst,
size_t level,
DoFType flag ) const
{
A.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
A.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
A.toMatrix( mat, src.uvw[2], dst.uvw[2], level, flag );
}
}
P1ConstantLaplaceOperator A;
};
......
src/hyteg/composites/P1StokesBlockPreconditioner.hpp
View file @ ab4ed4c0
#pragma once
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
namespace hyteg {
......@@ -16,6 +17,23 @@ class P1StokesBlockPreconditioner : public Operator< P1StokesFunction< real_t >,
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< idx_t >& src,
const P1StokesFunction< idx_t >& dst,
size_t level,
DoFType flag ) const
{
A.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
A.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
A.toMatrix( mat, src.uvw[2], dst.uvw[2], level, flag );
}
P.toMatrix( mat, src.p, dst.p, level, flag );
}
P1ConstantLaplaceOperator A;
P1LumpedMassOperator P;
......
src/hyteg/composites/P1StokesOperator.hpp
View file @ ab4ed4c0
......@@ -22,6 +22,7 @@
#include "hyteg/composites/P1StokesBlockPreconditioner.hpp"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/composites/StokesOperatorTraits.hpp"
#include "hyteg/operators/VectorLaplaceOperator.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1ScalarToP1VectorOperator.hpp"
#include "hyteg/p1functionspace/P1VectorToP1ScalarOperator.hpp"
......@@ -31,9 +32,10 @@ namespace hyteg {
class P1StokesOperator : public Operator< P1StokesFunction< real_t >, P1StokesFunction< real_t > >
{
public:
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
typedef P1ConstantVectorLaplaceOperator VelocityBlockOperator_T;
typedef P1ConstantLaplaceOperator VelocityOperator_T;
typedef P1ConstantLaplaceOperator PressureOperator_T;
typedef P1StokesBlockPreconditioner BlockPreconditioner_T;
P1StokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, size_t minLevel, size_t maxLevel )
: Operator( storage, minLevel, maxLevel )
......@@ -81,6 +83,34 @@ class P1StokesOperator : public Operator< P1StokesFunction< real_t >, P1StokesFu
pspg.apply( src.p, dst.p, level, flag, Add );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< idx_t >& src,
const P1StokesFunction< idx_t >& dst,
size_t level,
DoFType flag ) const
{
A.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
divT_x.toMatrix( mat, src.p, dst.uvw[0], level, flag );
A.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
divT_y.toMatrix( mat, src.p, dst.uvw[1], level, flag );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
A.toMatrix( mat, src.uvw[2], dst.uvw[2], level, flag );
divT_z.toMatrix( mat, src.p, dst.uvw[2], level, flag );
}
div_x.toMatrix( mat, src.uvw[0], dst.p, level, flag | DirichletBoundary );
div_y.toMatrix( mat, src.uvw[1], dst.p, level, flag | DirichletBoundary );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
div_z.toMatrix( mat, src.uvw[2], dst.p, level, flag | DirichletBoundary );
}
pspg.toMatrix( mat, src.p, dst.p, level, flag | DirichletBoundary );
}
P1ConstantLaplaceOperator A;
P1DivxOperator div_x;
P1DivyOperator div_y;
......
src/hyteg/composites/P2P1BlendingTaylorHoodStokesOperator.hpp
View file @ ab4ed4c0
......@@ -26,68 +26,66 @@
#include "hyteg/mixedoperators/P2ToP1VariableOperator.hpp"
#include "hyteg/p2functionspace/P2VariableOperator.hpp"
namespace hyteg {
class P2P1BlendingTaylorHoodStokesOperator
: public Operator<P2P1TaylorHoodFunction<real_t>, P2P1TaylorHoodFunction<real_t>>
class P2P1BlendingTaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< real_t >, P2P1TaylorHoodFunction< real_t > >
{
public:
typedef P2BlendingLaplaceOperator VelocityOperator_T;
P2P1BlendingTaylorHoodStokesOperator(const std::shared_ptr< PrimitiveStorage >& storage,
uint_t minLevel, uint_t maxLevel)
: Operator(storage, minLevel, maxLevel)
, A(storage, minLevel, maxLevel)
, div_x(storage, minLevel, maxLevel)
, div_y(storage, minLevel, maxLevel)
, div_z(storage, minLevel, maxLevel)
, divT_x(storage, minLevel, maxLevel)
, divT_y(storage, minLevel, maxLevel)
, divT_z(storage, minLevel, maxLevel)
, pspg_inv_diag_(storage, minLevel, maxLevel)
, hasGlobalCells_(storage->hasGlobalCells())
P2P1BlendingTaylorHoodStokesOperator( const std::shared_ptr< PrimitiveStorage >& storage, uint_t minLevel, uint_t maxLevel )