From b07e84cd1d604b10f75f06ae61f61fa5698b6cfa Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Thu, 1 Feb 2018 11:18:27 +0100
Subject: [PATCH] Correctly pass grid function types into the local operator
---
python/dune/perftool/pdelab/driver/solve.py | 6 +-
python/dune/perftool/pdelab/localoperator.py | 28 +++++++-
test/coeffeval/coeffeval_poisson.cc | 75 ++++++++++++++------
3 files changed, 84 insertions(+), 25 deletions(-)
diff --git a/python/dune/perftool/pdelab/driver/solve.py b/python/dune/perftool/pdelab/driver/solve.py
index 2e9b3ee1..29b0e268 100644
--- a/python/dune/perftool/pdelab/driver/solve.py
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -9,6 +9,7 @@ from dune.perftool.pdelab.driver import (get_form_ident,
name_initree,
)
from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+ type_domainfield,
type_trial_gfs,
)
from dune.perftool.pdelab.driver.constraints import (type_constraintscontainer,
@@ -82,8 +83,9 @@ def name_vector(form_ident):
@preamble
def typedef_vector(name, form_ident):
- go_type = type_gridoperator(form_ident)
- return "using {} = {}::Traits::Domain;".format(name, go_type)
+ gfs = type_trial_gfs()
+ df = type_domainfield()
+ return "using {} = Dune::PDELab::Backend::Vector<{},{}>;".format(name, gfs, df)
def type_vector(form_ident):
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index ea0dad53..336062fc 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -157,6 +157,26 @@ def localoperator_basename(form_ident):
return get_form_option("classname", form_ident)
+def name_gridfunction_member(coeff):
+ name = "gridfunction_{}".format(coeff.count())
+ define_gridfunction_member(name)
+ return name
+
+
+@class_member(classtag="operator")
+def define_gridfunction_member(name):
+ _type = type_gridfunction_template_parameter(name)
+ param = "{}_".format(name)
+ constructor_parameter("const {}&".format(_type), param, classtag="operator")
+ initializer_list(name, [param], classtag="operator")
+ return "const {}& {};".format(_type, name)
+
+
+@template_parameter(classtag="operator")
+def type_gridfunction_template_parameter(name):
+ return name.upper()
+
+
def class_type_from_cache(classtag):
from dune.perftool.generation import retrieve_cache_items
@@ -702,10 +722,14 @@ def generate_localoperator_kernels(operator):
lop_template_test_gfs()
lop_template_range_field()
- # Make sure there is always the same constructor arguments (even if parameter class is empty)
- from dune.perftool.pdelab.localoperator import name_initree_member
+ # Make sure there is always the same constructor arguments, even if some of them are
+ # not strictly needed. Also ensure the order.
name_initree_member()
+ # Iterate over the needed grid functions in correct order
+ for c in sorted(filter(lambda c: c.count() > 2, form.coefficients()), key=lambda c: c.count()):
+ name_gridfunction_member(c)
+
# Set some options!
from dune.perftool.pdelab.driver import isQuadrilateral
if isQuadrilateral(form.coefficients()[0].ufl_element().cell()):
diff --git a/test/coeffeval/coeffeval_poisson.cc b/test/coeffeval/coeffeval_poisson.cc
index 5fdf803b..cd80eb74 100644
--- a/test/coeffeval/coeffeval_poisson.cc
+++ b/test/coeffeval/coeffeval_poisson.cc
@@ -18,64 +18,95 @@
#include "poisson_localoperator.hh"
#include "dune/pdelab/stationary/linearproblem.hh"
#include "dune/grid/io/file/vtk/subsamplingvtkwriter.hh"
+#include "dune/pdelab/function/discretegridviewfunction.hh"
int main(int argc, char** argv)
{
+ // MPI helper stuff
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
- using VectorBackendP1 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>;
+
+ // Parse the ini file
+ Dune::ParameterTree initree;
+ Dune::ParameterTreeParser::readINITree(argv[1], initree);
+
+ // Build a grid
using Grid = Dune::ALUGrid<2, 2, Dune::simplex, Dune::conforming>;
using GV = Grid::LeafGridView;
+ IniGridFactory<Grid> factory(initree);
+ std::shared_ptr<Grid> grid = factory.getGrid();
+ GV gv = grid->leafGridView();
+
+ // General types and stuff
using DF = Grid::ctype;
using RangeType = double;
+
+ // Finite Element Maps
using P1_FEM = Dune::PDELab::PkLocalFiniteElementMap<GV, DF, RangeType, 1>;
+ P1_FEM p1_fem(gv);
+
+ // Grid Function Spaces
+ using VectorBackendP1 = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>;
using DirichletConstraintsAssember = Dune::PDELab::ConformingDirichletConstraints;
using P1_dirichlet_GFS = Dune::PDELab::GridFunctionSpace<GV, P1_FEM, DirichletConstraintsAssember, VectorBackendP1>;
- using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType>;
- using P1_dirichlet_GFS_CC = P1_dirichlet_GFS::ConstraintsContainer<RangeType>::Type;
- using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
- using GO_r = Dune::PDELab::GridOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, P1_dirichlet_GFS_CC, P1_dirichlet_GFS_CC>;
- using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
- using V_R = GO_r::Traits::Domain;
- using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>;
- Dune::ParameterTree initree;
- Dune::ParameterTreeParser::readINITree(argv[1], initree);
- IniGridFactory<Grid> factory(initree);
- std::shared_ptr<Grid> grid = factory.getGrid();
- GV gv = grid->leafGridView();
- P1_FEM p1_fem(gv);
P1_dirichlet_GFS p1_dirichlet_gfs_(gv, p1_fem);
p1_dirichlet_gfs_.name("p1_dirichlet_gfs_");
+ p1_dirichlet_gfs_.update();
+ std::cout << "gfs with " << p1_dirichlet_gfs_.size() << " dofs generated "<< std::endl;
+
+ // Solution vectors / Grid Functions
+ using V_R = Dune::PDELab::Backend::Vector<P1_dirichlet_GFS,DF>;
+ V_R x_r(p1_dirichlet_gfs_);
+ V_R sol(p1_dirichlet_gfs_);
+
+ using GF = Dune::PDELab::DiscreteGridViewFunction<P1_dirichlet_GFS, V_R>;
+ GF sol_gf(p1_dirichlet_gfs_, sol);
+
+ // Local Operator
+ using LOP_R = PoissonLocalOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, RangeType, GF>;
+ LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree, sol_gf);
+
+ // Constraints stuff
+ using P1_dirichlet_GFS_CC = P1_dirichlet_GFS::ConstraintsContainer<RangeType>::Type;
P1_dirichlet_GFS_CC p1_dirichlet_gfs__cc;
p1_dirichlet_gfs__cc.clear();
auto p1_bctype_lambda = [&](const auto& x){ return 1.0; };
auto p1_bctype = Dune::PDELab::makeBoundaryConditionFromCallable(gv, p1_bctype_lambda);
Dune::PDELab::constraints(p1_bctype, p1_dirichlet_gfs_, p1_dirichlet_gfs__cc);
- LOP_R lop_r(p1_dirichlet_gfs_, p1_dirichlet_gfs_, initree);
- p1_dirichlet_gfs_.update();
+ std::cout << "cc with " << p1_dirichlet_gfs__cc.size() << " dofs generated "<< std::endl;
+
+ // Matrix Backend
+ using MatrixBackend = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
int generic_dof_estimate = 6 * p1_dirichlet_gfs_.maxLocalSize();
int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
MatrixBackend mb(dofestimate);
+
+ // Grid Operator
+ using GO_r = Dune::PDELab::GridOperator<P1_dirichlet_GFS, P1_dirichlet_GFS, LOP_R, MatrixBackend, DF, RangeType, RangeType, P1_dirichlet_GFS_CC, P1_dirichlet_GFS_CC>;
GO_r go_r(p1_dirichlet_gfs_, p1_dirichlet_gfs__cc, p1_dirichlet_gfs_, p1_dirichlet_gfs__cc, lop_r, mb);
- std::cout << "gfs with " << p1_dirichlet_gfs_.size() << " dofs generated "<< std::endl;
- std::cout << "cc with " << p1_dirichlet_gfs__cc.size() << " dofs generated "<< std::endl;
+
+ // Solver
+ using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_SuperLU;
LinearSolver ls(false);
- V_R x_r(p1_dirichlet_gfs_);
- V_R sol(p1_dirichlet_gfs_);
+ using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO_r, LinearSolver, V_R>;
+
+ // Interpolation
x_r = 0.0;
std::size_t seed = 0;
auto rng = std::mt19937_64(seed);
auto dist = std::uniform_real_distribution<>(-1., 1.);
for (auto& v : sol)
v = dist(rng);
-
auto lambda_0000 = [&](const auto& x){ return (double)exp((-1.0) * ((0.5 - x[1]) * (0.5 - x[1]) + (0.5 - x[0]) * (0.5 - x[0]))); };
auto func_0000 = Dune::PDELab::makeGridFunctionFromCallable(gv, lambda_0000);
Dune::PDELab::interpolate(func_0000, p1_dirichlet_gfs_, sol);
+ // Solving
double reduction = initree.get<double>("reduction", 1e-12);
SLP slp(go_r, ls, x_r, reduction);
slp.apply();
+
+ // VTK visualization
using VTKWriter = Dune::SubsamplingVTKWriter<GV>;
Dune::RefinementIntervals subint(initree.get<int>("vtk.subsamplinglevel", 1));
VTKWriter vtkwriter(gv, subint);
@@ -83,6 +114,8 @@ int main(int argc, char** argv)
CuttingPredicate predicate;
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, p1_dirichlet_gfs_, x_r, Dune::PDELab::vtk::defaultNameScheme(), predicate);
vtkwriter.write(vtkfile, Dune::VTK::ascii);
+
+ // Error calculation
using P1_DIRICHLET_GFS__DGF = Dune::PDELab::DiscreteGridFunction<decltype(p1_dirichlet_gfs_),decltype(x_r)>;
P1_DIRICHLET_GFS__DGF p1_dirichlet_gfs__dgf(p1_dirichlet_gfs_,x_r);
RangeType l2error(0.0);
--
GitLab