Commit 0768f4a0 authored by Nils Kohl's avatar Nils Kohl 🌝
Browse files

[WIP] Started refactoring of sparse matrix creation API throughout the entire framework.

Implementing toMatrix() methods for many operators. Also got rid of some free functions
but this needs some more work. Almost all cases in PetscMatrixAssemblyTest are compiling
but many still abort. Also refactored PETScSparseMatrix a little so that only one template
argument is necessary now (only the operator).

Added support for non-square sparse matrices.

Some things that are still open:
- get PetscMatrixAssemblyTest green, enabling all test cases
- remove all former createMatrix() free function calls and definitions
- replace matIdx_t with the new index type before/when merging
- check that no PetscInts remain
- remove all the Petsc ifdefs that are not required anymore
- remove the petsc namespace where they technically only wrap sparse matrix stuff
- discuss what to do with the P1*Operator_new versions
parent ea620af6
Pipeline #34461 failed with stages
in 2 minutes and 6 seconds
......@@ -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"
......
......@@ -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"
......
......@@ -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"
......
......@@ -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< matIdx_t >& src,
const P1StokesFunction< matIdx_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;
......
......@@ -44,6 +44,21 @@ class P1StokesBlockLaplaceOperator : public Operator< P1StokesFunction< real_t >
}
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< matIdx_t >& src,
const P1StokesFunction< matIdx_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;
};
......
#pragma once
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/composites/P1StokesFunction.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
namespace hyteg {
......@@ -17,6 +17,23 @@ class P1StokesBlockPreconditioner : public Operator< P1StokesFunction< real_t >,
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P1StokesFunction< matIdx_t >& src,
const P1StokesFunction< matIdx_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;
......
......@@ -81,6 +81,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< matIdx_t >& src,
const P1StokesFunction< matIdx_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;
......
......@@ -19,6 +19,23 @@ class P2P1TaylorHoodStokesBlockPreconditioner
, hasGlobalCells_( storage->hasGlobalCells() )
{}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2P1TaylorHoodFunction< matIdx_t >& src,
const P2P1TaylorHoodFunction< matIdx_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 );
}
P2ConstantLaplaceOperator A;
P1LumpedMassOperator P;
......
......@@ -63,18 +63,52 @@ class P2P1TaylorHoodStokesOperator : public Operator< P2P1TaylorHoodFunction< re
div.apply( src.uvw, dst.p, level, flag, Replace );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2P1TaylorHoodFunction< matIdx_t >& src,
const P2P1TaylorHoodFunction< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
A.toMatrix( mat, src.uvw[0], dst.uvw[0], level, flag );
divT_x.getVertexToVertexOpr().toMatrix( mat, src.p, dst.uvw[0].getVertexDoFFunction(), level, flag );
divT_x.getVertexToEdgeOpr().toMatrix( mat, src.p, dst.uvw[0].getEdgeDoFFunction(), level, flag );
A.toMatrix( mat, src.uvw[1], dst.uvw[1], level, flag );
divT_y.getVertexToVertexOpr().toMatrix( mat, src.p, dst.uvw[1].getVertexDoFFunction(), level, flag );
divT_y.getVertexToEdgeOpr().toMatrix( mat, src.p, dst.uvw[1].getEdgeDoFFunction(), level, flag );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
A.toMatrix( mat, src.uvw[2], dst.uvw[2], level, flag );
divT_z.getVertexToVertexOpr().toMatrix( mat, src.p, dst.uvw[2].getVertexDoFFunction(), level, flag );
divT_z.getVertexToEdgeOpr().toMatrix( mat, src.p, dst.uvw[2].getEdgeDoFFunction(), level, flag );
}
div_x.getVertexToVertexOpr().toMatrix( mat, src.uvw[0].getVertexDoFFunction(), dst.p, level, flag | DirichletBoundary );
div_x.getEdgeToVertexOpr().toMatrix( mat, src.uvw[0].getEdgeDoFFunction(), dst.p, level, flag | DirichletBoundary );
div_y.getVertexToVertexOpr().toMatrix( mat, src.uvw[1].getVertexDoFFunction(), dst.p, level, flag | DirichletBoundary );
div_y.getEdgeToVertexOpr().toMatrix( mat, src.uvw[1].getEdgeDoFFunction(), dst.p, level, flag | DirichletBoundary );
if ( src.uvw[0].getStorage()->hasGlobalCells() )
{
div_z.getVertexToVertexOpr().toMatrix( mat, src.uvw[2].getVertexDoFFunction(), dst.p, level, flag | DirichletBoundary );
div_z.getEdgeToVertexOpr().toMatrix( mat, src.uvw[2].getEdgeDoFFunction(), dst.p, level, flag | DirichletBoundary );
}
}
P2ConstantVectorLaplaceOperator Lapl;
P2ToP1ConstantDivOperator div;
P1ToP2ConstantDivTOperator divT;
P2ToP1ConstantDivOperator div;
P1ToP2ConstantDivTOperator divT;
// currently need these for being able to call createMatrix()
P2ConstantLaplaceOperator A;
P2ToP1ConstantDivxOperator div_x;
P2ToP1ConstantDivyOperator div_y;
P2ToP1ConstantDivzOperator div_z;
P1ToP2ConstantDivTxOperator divT_x;
P1ToP2ConstantDivTyOperator divT_y;
P1ToP2ConstantDivTzOperator divT_z;
P2ConstantLaplaceOperator A;
P2ToP1ConstantDivxOperator div_x;
P2ToP1ConstantDivyOperator div_y;
P2ToP1ConstantDivzOperator div_z;
P1ToP2ConstantDivTxOperator divT_x;
P1ToP2ConstantDivTyOperator divT_y;
P1ToP2ConstantDivTzOperator divT_z;
/// this operator is need in the uzawa smoother
P1PSPGOperator pspg_;
......
/*
* Copyright (c) 2017-2019 Dominik Thoennes, Nils Kohl.
* Copyright (c) 2017-2021 Dominik Thoennes, Nils Kohl.
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
......@@ -72,6 +72,34 @@ class P2P2StabilizedStokesOperator : public Operator< P2P2StokesFunction< real_t
pspg.apply( src.p, dst.p, level, flag, Add );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2P2StokesFunction< matIdx_t >& src,
const P2P2StokesFunction< matIdx_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 );
}
P2ConstantLaplaceOperator A;
P2ConstantDivxOperator div_x;
P2ConstantDivyOperator div_y;
......
......@@ -48,7 +48,6 @@ class P2P2UnstableStokesOperator : public Operator< P2P2StokesFunction< real_t >
const size_t level,
DoFType flag ) const
{
WALBERLA_ASSERT_NOT_IDENTICAL( std::addressof( src ), std::addressof( dst ) );
A.apply( src.uvw[0], dst.uvw[0], level, flag, Replace );
......@@ -72,6 +71,32 @@ class P2P2UnstableStokesOperator : public Operator< P2P2StokesFunction< real_t >
}
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const P2P2StokesFunction< matIdx_t >& src,
const P2P2StokesFunction< matIdx_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 );
}
}
P2ConstantLaplaceOperator A;
P2ConstantDivxOperator div_x;
P2ConstantDivyOperator div_y;
......
......@@ -22,18 +22,17 @@
#include <type_traits>
#include <vector>
#include "hyteg/operators/Operator.hpp"
#include "hyteg/composites/P1StokesOperator.hpp"
#include "hyteg/composites/P1BlendingStokesOperator.hpp"
#include "hyteg/composites/P1StokesOperator.hpp"
#include "hyteg/composites/P2P1TaylorHoodStokesOperator.hpp"
#include "hyteg/composites/petsc/P1StokesPetsc.hpp"
#include "hyteg/composites/petsc/P2P1TaylorHoodPetsc.hpp"
#include "hyteg/petsc/PETScWrapper.hpp"
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
#include "hyteg/operators/Operator.hpp"
#include "hyteg/p1functionspace/P1ProjectNormalOperator.hpp"
#include "hyteg/p2functionspace/P2ProjectNormalOperator.hpp"
#include "hyteg/petsc/PETScWrapper.hpp"
#include "hyteg/sparseassembly/SparseMatrixProxy.hpp"
namespace hyteg {
......@@ -58,7 +57,6 @@ template < typename OpType, typename ProjOpType, bool PreProjection = false >
class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typename OpType::dstType >
{
public:
typedef typename OpType::BlockPreconditioner_T BlockPreconditioner_T;
StrongFreeSlipWrapper( std::shared_ptr< OpType > op, std::shared_ptr< ProjOpType > projOp, DoFType projFlag )
......@@ -81,7 +79,7 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
if ( PreProjection && projFlag_ != None )
{
tmp_->assign( {1}, {src}, level, All );
tmp_->assign( { 1 }, { src }, level, All );
projOp_->project( *tmp_, level, projFlag_ );
op_->apply( *tmp_, dst, level, flag );
}
......@@ -94,7 +92,6 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
projOp_->project( dst, level, projFlag_ );
}
#ifdef HYTEG_BUILD_WITH_PETSC
/// Assemble operator as sparse matrix
///
/// \param mat a sparse matrix proxy
......@@ -103,18 +100,18 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
/// \param level level in mesh hierarchy for which local operator is to be assembled
/// \param flag determines on which primitives this operator is assembled
///
void assembleLocalMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const typename OpType::srcType::template FunctionType< PetscInt >& numeratorSrc,
const typename OpType::dstType::template FunctionType< PetscInt >& numeratorDst,
uint_t level,
DoFType flag ) const
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const typename OpType::srcType::template FunctionType< PetscInt >& numeratorSrc,
const typename OpType::dstType::template FunctionType< PetscInt >& numeratorDst,
uint_t level,
DoFType flag ) const
{
auto matProxyOp = mat->createCopy();
hyteg::petsc::createMatrix< OpType >( *op_, numeratorSrc, numeratorDst, matProxyOp, level, flag );
op_->toMatrix( matProxyOp, numeratorSrc, numeratorDst, level, flag );
auto matProxyProjectionPost = mat->createCopy();
// projOp_->assembleLocalMatrix( matProxyProjectionPost, numeratorSrc.uvw[0], numeratorSrc.uvw[1], numeratorSrc.uvw[2], level, projFlag_ );
projOp_->assembleLocalMatrix( matProxyProjectionPost, numeratorSrc.uvw, level, projFlag_ );
projOp_->toMatrix( matProxyProjectionPost, numeratorSrc.uvw, numeratorDst.uvw, level, projFlag_ );
// we need the Id also in the pressure block
petsc::saveIdentityOperator( numeratorDst.p, matProxyProjectionPost, level, All );
......@@ -126,8 +123,8 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
if ( PreProjection )
{
auto matProxyProjectionPre = mat->createCopy();
// projOp_->assembleLocalMatrix( matProxyProjectionPre, numeratorSrc.uvw[0], numeratorSrc.uvw[1], numeratorSrc.uvw[2], level, projFlag_ );
projOp_->assembleLocalMatrix( matProxyProjectionPre, numeratorSrc.uvw, level, projFlag_ );
projOp_->toMatrix( matProxyProjectionPre, numeratorSrc.uvw, numeratorDst.uvw, level, projFlag_ );
petsc::saveIdentityOperator( numeratorDst.p, matProxyProjectionPre, level, All );
......@@ -136,7 +133,6 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
mat->createFromMatrixProduct( matrices );
}
#endif
private:
std::shared_ptr< OpType > op_;
......@@ -147,77 +143,4 @@ class StrongFreeSlipWrapper : public Operator< typename OpType::srcType, typenam
std::shared_ptr< typename OpType::dstType > tmp_;
};
#ifdef HYTEG_BUILD_WITH_PETSC
namespace petsc {
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P1BlendingStokesOperator, P1ProjectNormalOperator, true >& opr,
const P1StokesFunction< PetscInt >& src,
const P1StokesFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
WALBERLA_ABORT( "Not implemented." );
}
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P1BlendingStokesOperator, P1ProjectNormalOperator, false >& opr,
const P1StokesFunction< PetscInt >& src,
const P1StokesFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
WALBERLA_ABORT( "Not implemented." );
}
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P2P1TaylorHoodStokesOperator, P2ProjectNormalOperator, true >& opr,
const P2P1TaylorHoodFunction< PetscInt >& src,
const P2P1TaylorHoodFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
opr.assembleLocalMatrix( mat, src, dst, level, flag );
}
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P2P1TaylorHoodStokesOperator, P2ProjectNormalOperator, false >& opr,
const P2P1TaylorHoodFunction< PetscInt >& src,
const P2P1TaylorHoodFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
opr.assembleLocalMatrix( mat, src, dst, level, flag );
}
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P2P1ElementwiseBlendingStokesOperator, P2ProjectNormalOperator, true >& opr,
const P2P1TaylorHoodFunction< PetscInt >& src,
const P2P1TaylorHoodFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
opr.assembleLocalMatrix( mat, src, dst, level, flag );
}
template <>
inline void createMatrix( const StrongFreeSlipWrapper< P2P1ElementwiseBlendingStokesOperator, P2ProjectNormalOperator, false >& opr,
const P2P1TaylorHoodFunction< PetscInt >& src,
const P2P1TaylorHoodFunction< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
size_t level,
DoFType flag )
{
opr.assembleLocalMatrix( mat, src, dst, level, flag );
}
} // namespace petsc
#endif
} // namespace hyteg
......@@ -26,7 +26,6 @@
#include "hyteg/edgedofspace/EdgeDoFIndexing.hpp"
#include "hyteg/edgedofspace/EdgeDoFMacroEdge.hpp"
#include "hyteg/edgedofspace/EdgeDoFMacroFace.hpp"
#include "hyteg/elementwiseoperators/ElementwiseOperatorPetsc.hpp"
#include "hyteg/elementwiseoperators/P1ElementwiseOperator.hpp"
#include "hyteg/elementwiseoperators/P2ElementwiseOperator.hpp"
#include "hyteg/forms/P1LinearCombinationForm.hpp"
......@@ -205,6 +204,15 @@ class UnsteadyDiffusionOperator : public Operator< FunctionType, FunctionType >,
LinearCombinationForm_T( { 1.0, dtScaling() * dt * diffusionCoefficient_ }, { massForm_.get(), laplaceForm_.get() } ) );
}
void toMatrix( const std::shared_ptr< SparseMatrixProxy >& mat,
const typename Operator_T< LinearCombinationForm_T >::srcType::template FunctionType< matIdx_t >& src,
const typename Operator_T< LinearCombinationForm_T >::srcType::template FunctionType< matIdx_t >& dst,
size_t level,
DoFType flag ) const
{
unsteadyDiffusionOperator_->toMatrix( mat, src, dst, level, flag );
}
DiffusionTimeIntegrator getTimeIntegrator() const { return timeIntegrator_; }
const Operator_T< LinearCombinationForm_T >& getOperator() const { return *unsteadyDiffusionOperator_; }
......@@ -263,44 +271,6 @@ typedef UnsteadyDiffusionOperator< P2Function< real_t >,
P2LinearCombinationForm >
P2ElementwiseUnsteadyDiffusionOperator;
#ifdef HYTEG_BUILD_WITH_PETSC
namespace petsc {
template <>
inline void createMatrix< P1ConstantUnsteadyDiffusionOperator >( const P1ConstantUnsteadyDiffusionOperator& opr,
const P1Function< PetscInt >& src,
const P1Function< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
uint_t level,
DoFType flag )
{
createMatrix( opr.getOperator(), src, dst, mat, level, flag );
}
template <>
inline void createMatrix< P2ConstantUnsteadyDiffusionOperator >( const P2ConstantUnsteadyDiffusionOperator& opr,
const P2Function< PetscInt >& src,
const P2Function< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
uint_t level,
DoFType flag )
{
createMatrix( opr.getOperator(), src, dst, mat, level, flag );
}
template <>
inline void createMatrix< P2ElementwiseUnsteadyDiffusionOperator >( const P2ElementwiseUnsteadyDiffusionOperator& opr,
const P2Function< PetscInt >& src,
const P2Function< PetscInt >& dst,
const std::shared_ptr< SparseMatrixProxy >& mat,
uint_t level,
DoFType flag )
{
createMatrix< P2ElementwiseOperator< P2LinearCombinationForm >, P2LinearCombinationForm >(
opr.getOperator(), src, dst, mat, level, flag );
}
} // namespace petsc
#endif