Commit cc9f10bc authored by Benjamin Mann's avatar Benjamin Mann
Browse files

Merge branch 'mogli/refactoring_P1Operators' into 'master'

Replace usage of P1PolynomialBlendingOperator

See merge request !454
parents 7f615a6f 36d5c898
Pipeline #34726 passed with stages
in 190 minutes and 2 seconds
This diff is collapsed.
......@@ -32,7 +32,6 @@
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1VariableOperator.hpp"
#include "hyteg/p1functionspace/P1PolynomialBlendingOperator.hpp"
#include "hyteg/p2functionspace/P2Function.hpp"
#include "hyteg/p2functionspace/P2VariableOperator.hpp"
......
......@@ -29,7 +29,7 @@
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/P1SurrogateOperator.hpp"
#include "hyteg/p1functionspace/P1VariableOperator_new.hpp"
#include "hyteg/p1functionspace/P1VariableOperator.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/primitivestorage/Visualization.hpp"
......@@ -249,8 +249,8 @@ int main( int argc, char* argv[] )
// ===================
logSectionHeader( "Standard FEM Part" );
logMessage( "Preparing P1VariableOperator_new" );
P1AffineDivkGradOperator_new varOp( storage, minLevel, maxLevel, form );
logMessage( "Preparing P1VariableOperator" );
P1AffineDivkGradOperator varOp( storage, minLevel, maxLevel, form );
logMessage( "Preparing Stencil Generation" );
const PrimitiveDataID< StencilMemory< real_t >, Face >& stencilID = varOp.getFaceStencilID();
......
......@@ -27,10 +27,10 @@
#include "hyteg/functions/FunctionProperties.hpp"
#include "hyteg/geometry/AnnulusMap.hpp"
#include "hyteg/geometry/IcosahedralShellMap.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator_new.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/P1SurrogateOperator.hpp"
#include "hyteg/p1functionspace/P1VariableOperator_new.hpp"
#include "hyteg/p1functionspace/P1VariableOperator.hpp"
#include "hyteg/primitivestorage/PrimitiveStorage.hpp"
#include "hyteg/primitivestorage/SetupPrimitiveStorage.hpp"
#include "hyteg/primitivestorage/loadbalancing/SimpleBalancer.hpp"
......@@ -49,6 +49,28 @@ enum KerType
SURROGATE_FAST,
};
// enable access to protected member functions
template < class Op >
class OperatorWrapper : public Op
{
public:
using Op::Op;
using Op::apply_cell;
using Op::apply_cell_generated;
using Op::apply_face;
using Op::apply_face_generated;
using Op::smooth_sor_cell;
using Op::smooth_sor_cell_generated;
using Op::smooth_sor_face;
using Op::smooth_sor_face_generated;
};
using LaplaceConst = OperatorWrapper< P1ConstantLaplaceOperator >;
using LaplaceBlending = OperatorWrapper< hyteg::P1VariableOperator< hyteg::forms::p1_diffusion_blending_q1 > >;
using LaplaceAffine = OperatorWrapper< hyteg::P1VariableOperator< hyteg::forms::p1_diffusion_affine_q1 > >;
using LaplaceSurrogate = OperatorWrapper< P1SurrogateOperator< hyteg::forms::p1_diffusion_blending_q1, false > >;
using LaplaceSurrogateFast = OperatorWrapper< P1SurrogateOperator< hyteg::forms::p1_diffusion_blending_q1, true > >;
const std::map< std::string, KerType > strToKerType = { { "variable", VARIABLE },
{ "constant", CONSTANT },
{ "generated", GENERATED },
......@@ -191,11 +213,11 @@ void benchmark( KerType kertype, uint_t dim, bool blending, uint_t q, uint_t lev
WALBERLA_LOG_INFO_ON_ROOT( "Intitialize operators" );
P1ConstantLaplaceOperator_new L_const( storage, level, level );
hyteg::P1VariableOperator_new< hyteg::forms::p1_diffusion_blending_q1 > L_blend( storage, level, level );
hyteg::P1VariableOperator_new< hyteg::forms::p1_diffusion_affine_q1 > L_aff( storage, level, level );
P1SurrogateOperator< hyteg::forms::p1_diffusion_blending_q1, false > L_q( storage, level, level );
P1SurrogateOperator< hyteg::forms::p1_diffusion_blending_q1, true > L_q_fast( storage, level, level );
LaplaceConst L_const( storage, level, level );
LaplaceBlending L_blend( storage, level, level );
LaplaceAffine L_aff( storage, level, level );
LaplaceSurrogate L_q( storage, level, level );
LaplaceSurrogateFast L_q_fast( storage, level, level );
L_q.interpolateStencils( q, sample );
L_q_fast.interpolateStencils( q, sample );
......
......@@ -31,7 +31,7 @@
#include "hyteg/p1functionspace/P1Function.hpp"
#include "hyteg/p1functionspace/P1ConstantOperator.hpp"
#include "hyteg/p1functionspace/P1VariableOperator.hpp"
#include "hyteg/p1functionspace/P1PolynomialBlendingOperator.hpp"
#include "hyteg/p1functionspace/P1SurrogateOperator.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearRestriction.hpp"
#include "hyteg/gridtransferoperators/P1toP1LinearProlongation.hpp"
#include "hyteg/gridtransferoperators/P1toP1QuadraticProlongation.hpp"
......@@ -140,7 +140,7 @@ int main(int argc, char* argv[])
typedef hyteg::P1BlendingMassOperator MassOperator;
typedef hyteg::P1BlendingLaplaceOperator SolveOperatorNodal;
typedef hyteg::P1PolynomialBlendingLaplaceOperator SolveOperatorPoly;
typedef hyteg::P1SurrogateLaplaceOperator SolveOperatorPoly;
std::function<real_t(const hyteg::Point3D&)> exact = [](const hyteg::Point3D& x) { return sin(x[0])*sinh(x[1]); };
// std::function<real_t(const hyteg::Point3D&)> rhs = [](const hyteg::Point3D& x) { return -2*(x[0] + 1)*cos(x[0])*sinh(x[1]) - 3*sin(x[0])*cosh(x[1]); };
......@@ -166,17 +166,15 @@ int main(int argc, char* argv[])
MassOperator M(storage, minLevel, maxMemoryLevel);
std::shared_ptr<SolveOperatorPoly> Lpoly;
std::shared_ptr<SolveOperatorNodal> L;
uint_t useDegree;
uint_t useDegree = maxPolyDegree;
auto start = walberla::timing::getWcTime();
L = std::make_shared<SolveOperatorNodal>(storage, minLevel, maxMemoryLevel);
if (polynomialOperator) {
Lpoly = std::make_shared<SolveOperatorPoly>(storage, minLevel, maxLevel, interpolationLevel);
Lpoly = std::make_shared<SolveOperatorPoly>(storage, minLevel, maxLevel);
Lpoly->interpolateStencils(maxPolyDegree);
useDegree = maxPolyDegree;
Lpoly->useDegree(useDegree);
Lpoly->interpolateStencils(useDegree, interpolationLevel);
// real_t polyError34 = Lpoly->lInfinityError(3, 4, 5);
// real_t polyError45 = Lpoly->lInfinityError(4, 5, 5);
......@@ -310,8 +308,7 @@ int main(int argc, char* argv[])
WALBERLA_LOG_INFO_ON_ROOT("Increasing polynomial degree to " << useDegree + 1);
++useDegree;
Lpoly->interpolateStencils(maxLevel, useDegree);
Lpoly->useDegree(useDegree);
Lpoly->interpolateStencils(useDegree, interpolationLevel);
updatedDegree = true;
} else {
updatedDegree = false;
......
......@@ -7,12 +7,10 @@ Parameters
elementType 1;
//operatorType:
// 0: constant;
// 0: old constant;
// 1: variable;
// 2: polynomial;
// 3: new_constant;
// 4: new_variable;
// 5: new_polynomial;
// 3: new constant;
operatorType 4;
// parameters for surrogate polynomials;
......
......@@ -81,6 +81,16 @@ class P1LinearCombinationForm : public P1Form
}
}
void integrateRow0( const std::array< Point3D, 3 >& coords, Matrixr< 1, 3 >& elMat ) const override
{
Point3D row;
integrate(coords, row);
for (int i = 0; i < 3; ++i)
{
elMat(0,i) = row[i];
}
}
void integrateAll( const std::array< Point3D, 3 >& coords, Matrix3r& elMat ) const
{
elMat.setAll( 0 );
......@@ -118,6 +128,16 @@ class P1LinearCombinationForm : public P1Form
}
}
void integrateRow0( const std::array< Point3D, 4 >& coords, Matrixr< 1, 4 >& elMat ) const override
{
Point4D row;
integrate(coords, row);
for (int i = 0; i < 4; ++i)
{
elMat(0,i) = row[i];
}
}
void integrateAll( const std::array< Point3D, 4 >& coords, Matrix4r& elMat ) const
{
elMat.setAll( 0 );
......
......@@ -113,6 +113,27 @@ class P1RowSumForm : public P1Form
bool assembly3DDefined() const override { return true; }
private:
void integrateRow0( const std::array< Point3D, 3 >& coords, Matrixr< 1, 3 >& elMat ) const override
{
Point3D row;
integrate( coords, row );
for ( int i = 0; i < 3; ++i )
{
elMat( 0, i ) = row[i];
}
}
void integrateRow0( const std::array< Point3D, 4 >& coords, Matrixr< 1, 4 >& elMat ) const override
{
Point4D row;
integrate( coords, row );
for ( int i = 0; i < 4; ++i )
{
elMat( 0, i ) = row[i];
}
}
std::shared_ptr< P1Form > form_;
};
......
/*
* Copyright (c) 2021 Benjamin Mann
*
* This file is part of HyTeG
* (see https://i10git.cs.fau.de/hyteg/hyteg).
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#pragma once
namespace hyteg {
// wrapper to enable combining P1 operators with higher order forms
template < class WrappedForm >
class P1WrapperForm : public WrappedForm
{
public:
// use Ctor from WrappedForm
using WrappedForm::WrappedForm;
// conversion WrappedForm -> P1Wrapper
P1WrapperForm( const WrappedForm& wf )
: WrappedForm( wf )
{}
// conversion WrappedForm -> P1Wrapper
P1WrapperForm( WrappedForm&& wf )
: WrappedForm( wf )
{}
void integrateRow( const uint_t& row, const std::array< Point3D, 3 >& coords, Matrixr< 1, 3 >& elMat ) const
{
switch ( row )
{
case 0:
integrateRow0( coords, elMat );
break;
default:
WALBERLA_ABORT( "P1Form::integrateRow() not implemented for row " << row );
}
}
void integrateRow( const uint_t& row, const std::array< Point3D, 4 >& coords, Matrixr< 1, 4 >& elMat ) const
{
switch ( row )
{
case 0:
integrateRow0( coords, elMat );
break;
default:
WALBERLA_ABORT( "P1Form::integrateRow() not implemented for row " << row );
}
}
private:
// extract vertex to vertex DoF part of WrappedForm
void integrateRow0( const std::array< Point3D, 3 >& coords, Matrixr< 1, 3 >& elMat ) const
{
Point3D row;
this->integrate( coords, row );
for ( int i = 0; i < 3; ++i )
{
elMat( 0, i ) = row[i];
}
}
// extract vertex to vertex DoF part of WrappedForm
void integrateRow0( const std::array< Point3D, 4 >& coords, Matrixr< 1, 4 >& elMat ) const
{
Point4D row;
this->integrate( coords, row );
for ( int i = 0; i < 4; ++i )
{
elMat( 0, i ) = row[i];
}
}
};
} // namespace hyteg
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_0_affine_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_0_blending_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_1_affine_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_1_blending_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_2_affine_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_div_2_blending_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_divt_0_affine_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_divt_0_blending_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_divt_1_affine_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
......@@ -20,7 +20,7 @@
/*
* The entire file was generated with the HyTeG form generator.
*
*
* Software:
*
* - quadpy version: 0.16.5
......@@ -39,7 +39,7 @@ namespace forms {
/// Implementation of the integration of a weak form over an element.
///
/// - name: p1_divt_1_blending_q1
/// - description:
/// - description:
/// - trial space: Lagrange, degree: 1
/// - test space: Lagrange, degree: 1
///
......
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