From 98f285b174cd9c4b2663cc0bf9ba646de07fe74b Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 3 May 2016 17:16:00 +0200
Subject: [PATCH] Poisson DG test cases + reference solutions
---
python/dune/perftool/pdelab/driver.py | 13 +-
python/dune/perftool/pdelab/parameter.py | 5 +-
test/laplace/CMakeLists.txt | 3 +-
test/laplace/laplace_dg_numdiff.mini | 1 +
test/laplace/laplace_dg_symdiff.mini | 2 +
test/laplace/laplace_numdiff.mini | 1 +
test/laplace/laplace_symdiff.mini | 1 +
test/laplace/reference_driver.hh | 144 ++++++++++++++++++++++
test/laplace/reference_main.cc | 33 ++++++
test/poisson/CMakeLists.txt | 45 ++++++-
test/poisson/poisson_dg.ufl | 22 ++--
test/poisson/poisson_dg_numdiff.mini | 11 ++
test/poisson/poisson_dg_ref.vtu | 3 +
test/poisson/reference_driver.hh | 145 +++++++++++++++++++++++
test/poisson/reference_main.cc | 33 ++++++
15 files changed, 440 insertions(+), 22 deletions(-)
create mode 100644 test/laplace/reference_driver.hh
create mode 100644 test/laplace/reference_main.cc
create mode 100644 test/poisson/poisson_dg_numdiff.mini
create mode 100644 test/poisson/poisson_dg_ref.vtu
create mode 100644 test/poisson/reference_driver.hh
create mode 100644 test/poisson/reference_main.cc
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index 91d0e76d..46d54a8d 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -774,14 +774,17 @@ def name_vtkfile():
@preamble
def print_matrix():
+ ini = name_initree()
t_go = type_gridoperator()
n_go = name_gridoperator()
v = name_vector()
- return ["typedef typename {}::Traits::Jacobian M;".format(t_go),
- "M m({});".format(n_go),
- "{}.jacobian({},m);".format(n_go, v),
- "using Dune::PDELab::Backend::native;",
- "Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);"]
+ return ["if ({}.get<bool>(\"printmatrix\", false)) {{".format(ini),
+ " typedef typename {}::Traits::Jacobian M;".format(t_go),
+ " M m({});".format(n_go),
+ " {}.jacobian({},m);".format(n_go, v),
+ " using Dune::PDELab::Backend::native;",
+ " Dune::printmatrix(std::cout, native(m),\"global stiffness matrix\",\"row\",9,1);",
+ "}"]
@preamble
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index abffb46d..d8ac5fd7 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -1,6 +1,7 @@
""" Generators for parameter functions """
-from dune.perftool.generation import (class_basename,
+from dune.perftool.generation import (cached,
+ class_basename,
class_member,
constructor_parameter,
generator_factory,
@@ -105,12 +106,14 @@ def define_parameter_temporary(t):
return _define
+@cached
def cell_parameter_function(name, expr, restriction, t='double'):
temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
define_parameter_function_class_member(name, expr, t, True)
evaluate_cell_parameter_function(name, restriction)
+@cached
def intersection_parameter_function(name, expr, t='double'):
temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
define_parameter_function_class_member(name, expr, t, False)
diff --git a/test/laplace/CMakeLists.txt b/test/laplace/CMakeLists.txt
index d430d221..d58b9c54 100644
--- a/test/laplace/CMakeLists.txt
+++ b/test/laplace/CMakeLists.txt
@@ -30,4 +30,5 @@ add_generated_executable(UFLFILE laplace_dg.ufl
dune_add_system_test(TARGET laplace_dg_symdiff
INIFILE laplace_dg_symdiff.mini)
-add_executable(dgref reference_main.cc)
+add_executable(laplace_dg_ref reference_main.cc)
+set_target_properties(laplace_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/laplace/laplace_dg_numdiff.mini b/test/laplace/laplace_dg_numdiff.mini
index 269235ec..ba1e0014 100644
--- a/test/laplace/laplace_dg_numdiff.mini
+++ b/test/laplace/laplace_dg_numdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 2 2
elementType = simplical
+printmatrix = true
diff --git a/test/laplace/laplace_dg_symdiff.mini b/test/laplace/laplace_dg_symdiff.mini
index 6544b4d1..b03b5ada 100644
--- a/test/laplace/laplace_dg_symdiff.mini
+++ b/test/laplace/laplace_dg_symdiff.mini
@@ -4,3 +4,5 @@ lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 2 2
elementType = simplical
+printmatrix = true
+
diff --git a/test/laplace/laplace_numdiff.mini b/test/laplace/laplace_numdiff.mini
index 60d37cd0..8ee48e24 100644
--- a/test/laplace/laplace_numdiff.mini
+++ b/test/laplace/laplace_numdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 4 4
elementType = simplical
+printmatrix = true
diff --git a/test/laplace/laplace_symdiff.mini b/test/laplace/laplace_symdiff.mini
index ec27540f..c0ee0dd3 100644
--- a/test/laplace/laplace_symdiff.mini
+++ b/test/laplace/laplace_symdiff.mini
@@ -4,3 +4,4 @@ lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 4 4
elementType = simplical
+printmatrix = true
diff --git a/test/laplace/reference_driver.hh b/test/laplace/reference_driver.hh
new file mode 100644
index 00000000..db1b5594
--- /dev/null
+++ b/test/laplace/reference_driver.hh
@@ -0,0 +1,144 @@
+#ifndef _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+#define _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <string>
+#include <dune/common/parametertree.hh>
+#include <dune/pdelab/finiteelementmap/opbfem.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/pdelab/stationary/linearproblem.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/testtools/gridconstruction.hh>
+#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
+
+template<typename GV, typename RF>
+class CDProb
+{
+ typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+ public:
+ typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+ //! tensor diffusion coefficient
+ typename Traits::PermTensorType
+ A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ typename Traits::PermTensorType I;
+ for (std::size_t i=0; i<Traits::dimDomain; i++)
+ for (std::size_t j=0; j<Traits::dimDomain; j++)
+ I[i][j] = (i==j) ? 1 : 0;
+ return I;
+ }
+
+ //! velocity field
+ typename Traits::RangeType
+ b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ typename Traits::RangeType v(0.0);
+ return v;
+ }
+
+ //! sink term
+ typename Traits::RangeFieldType
+ c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! source term
+ typename Traits::RangeFieldType
+ f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! boundary condition type function
+ BCType
+ bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+ }
+
+ //! Dirichlet boundary condition value
+ typename Traits::RangeFieldType
+ g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! Neumann boundary condition
+ typename Traits::RangeFieldType
+ j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! outflow boundary condition
+ typename Traits::RangeFieldType
+ o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return 0.0;
+ }
+};
+
+
+void driver(int argc, char** argv){ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1> VectorBackend;
+ static const int dim = 2;
+ typedef Dune::UGGrid<dim> Grid;
+ typedef Grid::LeafGridView GV;
+ typedef Grid::ctype DF;
+ typedef double R;
+ typedef Dune::PDELab::OPBLocalFiniteElementMap<DF, R, 1, dim, Dune::GeometryType::simplex> DG1_FEM;
+ typedef Dune::PDELab::NoConstraints NoConstraintsAssembler;
+ typedef Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackend> DG1_DIRICHLET_GFS;
+ Dune::ParameterTree initree;
+ Dune::ParameterTreeParser::readINITree(argv[1], initree);
+ IniGridFactory<Grid> factory(initree);
+ std::shared_ptr<Grid> grid = factory.getGrid();
+ GV gv = grid->leafGridView();
+ DG1_FEM dg1_fem;
+ DG1_DIRICHLET_GFS dg1_dirichlet_gfs(gv, dg1_fem);
+ dg1_dirichlet_gfs.name("bla");
+ typedef Dune::SubsamplingVTKWriter<GV> VTKWriter;
+ int sublevel = initree.get<int>("vtk.subsamplinglevel", 0);
+ VTKWriter vtkwriter(gv, sublevel);
+ using LOP = Dune::PDELab::ConvectionDiffusionDG<CDProb<GV, R>, DG1_FEM>;
+ typedef DG1_DIRICHLET_GFS::ConstraintsContainer<R>::Type DG1_CC;
+ typedef Dune::PDELab::istl::BCRSMatrixBackend<> MatrixBackend;
+ typedef Dune::PDELab::GridOperator<DG1_DIRICHLET_GFS, DG1_DIRICHLET_GFS, LOP, MatrixBackend, DF, R, R, DG1_CC, DG1_CC> GO;
+ typedef GO::Traits::Domain V;
+ V x(dg1_dirichlet_gfs);
+ x = 0.0;
+ std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
+ typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LinearSolver;
+ typedef Dune::PDELab::StationaryLinearProblemSolver<GO, LinearSolver, V> SLP;
+ DG1_CC dg1_cc;
+ dg1_cc.clear();
+ CDProb<GV, R> params;
+ LOP lop(params, Dune::PDELab::ConvectionDiffusionDGMethod::SIPG);
+ int generic_dof_estimate = 6 * dg1_dirichlet_gfs.maxLocalSize();
+ int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
+ MatrixBackend mb(dofestimate);
+ GO go(dg1_dirichlet_gfs, dg1_cc, dg1_dirichlet_gfs, dg1_cc, lop, mb);
+ std::cout << "gfs with " << dg1_dirichlet_gfs.size() << " dofs generated "<< std::endl;
+ std::cout << "cc with " << dg1_cc.size() << " dofs generated "<< std::endl;
+ LinearSolver ls(false);
+ double reduction = initree.get<double>("reduction", 1e-12);
+ SLP slp(go, ls, x, reduction);
+ slp.apply();
+ typedef typename GO::Traits::Jacobian M;
+ M m(go);
+ go.jacobian(x,m);
+ using Dune::PDELab::Backend::native;
+ Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
+ Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_dirichlet_gfs, x);
+ vtkwriter.write(vtkfile, Dune::VTK::ascii);
+}
+
+#endif //_HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
diff --git a/test/laplace/reference_main.cc b/test/laplace/reference_main.cc
new file mode 100644
index 00000000..ff00d5c0
--- /dev/null
+++ b/test/laplace/reference_main.cc
@@ -0,0 +1,33 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+
+#include"/home/dominic/dune/dune-perftool/test/laplace/reference_driver.hh"
+
+int main(int argc, char** argv)
+{
+ try{
+ //Maybe initialize Mpi
+ Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+ if(Dune::MPIHelper::isFake)
+ std::cout<< "This is a sequential program." << std::endl;
+ else
+ std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+ <<" processes!"<<std::endl;
+
+ driver(argc, argv);
+
+ return 0;
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 1;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ return 1;
+ }
+}
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index d18236b9..e00a32f4 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -1,3 +1,9 @@
+# Symlink any reference solutions
+dune_symlink_to_source_files(FILES poisson_ref.vtu
+ poisson_dg_ref.vtu
+ )
+
+
# 1. Poisson Test Case: source f, bc g + numerical differentiation
add_generated_executable(UFLFILE poisson.ufl
TARGET poisson_numdiff
@@ -9,9 +15,6 @@ dune_add_system_test(TARGET poisson_numdiff
SCRIPT dune_vtkcompare.py
)
-# Symlink any reference solutions
-dune_symlink_to_source_files(FILES poisson_ref.vtu)
-
# 2. Poisson Test Case: source f, bc g + symbolic differentiation
add_generated_executable(UFLFILE poisson.ufl
@@ -42,3 +45,39 @@ add_generated_executable(UFLFILE poisson_neumann.ufl
dune_add_system_test(TARGET poisson_neumann_symdiff
INIFILE poisson_neumann_symdiff.mini
)
+
+# 5. Poisson Test Case: DG, f + pure dirichlet, numeric differentiation
+add_generated_executable(UFLFILE poisson_dg.ufl
+ TARGET poisson_dg_numdiff
+ FORM_COMPILER_ARGS --numerical-jacobian
+ )
+
+dune_add_system_test(TARGET poisson_dg_numdiff
+ INIFILE poisson_dg_numdiff.mini
+ SCRIPT dune_vtkcompare.py
+ )
+
+# 6. Poisson Test Case: DG, f + pure dirichlet, symbolic differentiation
+add_generated_executable(UFLFILE poisson_dg.ufl
+ TARGET poisson_dg_symdiff
+ )
+
+dune_add_system_test(TARGET poisson_dg_symdiff
+ INIFILE poisson_dg_symdiff.mini
+ SCRIPT dune_vtkcompare.py
+ )
+
+# 7. Poisson Test Case: DG, mixed bc, numeric differentiation
+add_generated_executable(UFLFILE poisson_dg_neumann.ufl
+ TARGET poisson_dg_neumann_numdiff
+ FORM_COMPILER_ARGS --numerical-jacobian
+ )
+
+dune_add_system_test(TARGET poisson_dg_neumann_numdiff
+ INIFILE poisson_dg_neumann_numdiff.mini
+ )
+
+# Add the reference code with the PDELab localoperator that produced
+# the reference vtk file
+add_executable(poisson_dg_ref reference_main.cc)
+set_target_properties(poisson_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
diff --git a/test/poisson/poisson_dg.ufl b/test/poisson/poisson_dg.ufl
index 416ccf8f..0ee5a857 100644
--- a/test/poisson/poisson_dg.ufl
+++ b/test/poisson/poisson_dg.ufl
@@ -1,7 +1,5 @@
f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
-g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
-j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 1; else return 0;", on_intersection=True)
+g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());", on_intersection=True)
cell = triangle
V = FiniteElement("DG", cell, 1)
@@ -9,20 +7,20 @@ V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
-n = FacetNormal(cell)
+n = FacetNormal(cell)('+')
gamma = 1.0
theta = 1.0
-ds = ds(subdomain_data=bctype)
-
r = inner(grad(u), grad(v))*dx \
- - inner(n, avg(grad(u)))*jump(v)*(dS+ds(0)) \
- + gamma*jump(u)*jump(v)*(dS+ds(0)) \
- - theta*jump(u)*inner(grad(v), n)*(dS+ds(0)) \
+ + inner(n, avg(grad(u)))*jump(v)*dS \
+ + gamma*jump(u)*jump(v)*dS \
+ + theta*jump(u)*inner(avg(grad(v)), n)*dS \
+ - inner(n, grad(u))*v*ds \
+ + gamma*u*v*ds \
+ - theta*u*inner(grad(v), n)*ds \
- f*v*dx \
- + theta*g*inner(grad(v), n)*ds(0) \
- - gamma*g*v*ds(0) \
- + j*v*ds(1)
+ + theta*g*inner(grad(v), n)*ds \
+ - gamma*g*v*ds
forms = [r]
diff --git a/test/poisson/poisson_dg_numdiff.mini b/test/poisson/poisson_dg_numdiff.mini
new file mode 100644
index 00000000..146955e8
--- /dev/null
+++ b/test/poisson/poisson_dg_numdiff.mini
@@ -0,0 +1,11 @@
+__name = poisson_dg_numdiff
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 32 32
+elementType = simplical
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = poisson_dg_ref
+extension = vtu
diff --git a/test/poisson/poisson_dg_ref.vtu b/test/poisson/poisson_dg_ref.vtu
new file mode 100644
index 00000000..06396d8a
--- /dev/null
+++ b/test/poisson/poisson_dg_ref.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:34ba9a3065699851e3bdc68c2ef14f525c0bdfb42fc46c07c4218109e2d30c97
+size 226320
diff --git a/test/poisson/reference_driver.hh b/test/poisson/reference_driver.hh
new file mode 100644
index 00000000..6c3e45f2
--- /dev/null
+++ b/test/poisson/reference_driver.hh
@@ -0,0 +1,145 @@
+#ifndef _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+#define _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
+
+#include <dune/pdelab/gridoperator/gridoperator.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <dune/pdelab/gridfunctionspace/vtk.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/pdelab/backend/istl.hh>
+#include <string>
+#include <dune/common/parametertree.hh>
+#include <dune/pdelab/finiteelementmap/opbfem.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/pdelab/stationary/linearproblem.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/testtools/gridconstruction.hh>
+#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
+#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
+
+template<typename GV, typename RF>
+class CDProb
+{
+ typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
+
+ public:
+ typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
+
+ //! tensor diffusion coefficient
+ typename Traits::PermTensorType
+ A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ typename Traits::PermTensorType I;
+ for (std::size_t i=0; i<Traits::dimDomain; i++)
+ for (std::size_t j=0; j<Traits::dimDomain; j++)
+ I[i][j] = (i==j) ? 1 : 0;
+ return I;
+ }
+
+ //! velocity field
+ typename Traits::RangeType
+ b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ typename Traits::RangeType v(0.0);
+ return v;
+ }
+
+ //! sink term
+ typename Traits::RangeFieldType
+ c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! source term
+ typename Traits::RangeFieldType
+ f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ auto global = e.geometry().global(x);
+ Dune::FieldVector<double,2> c(0.5);
+ c -= global;
+ return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());
+ }
+
+ //! boundary condition type function
+ BCType
+ bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
+ }
+
+ //! Dirichlet boundary condition value
+ typename Traits::RangeFieldType
+ g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
+ {
+ auto global = e.geometry().global(x);
+ Dune::FieldVector<double,2> c(0.5);
+ c-= global;
+ return std::exp(-1.*c.two_norm2());
+ }
+
+ //! Neumann boundary condition
+ typename Traits::RangeFieldType
+ j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return 0.0;
+ }
+
+ //! outflow boundary condition
+ typename Traits::RangeFieldType
+ o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
+ {
+ return 0.0;
+ }
+};
+
+
+void driver(int argc, char** argv){ typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none, 1> VectorBackend;
+ static const int dim = 2;
+ typedef Dune::UGGrid<dim> Grid;
+ typedef Grid::LeafGridView GV;
+ typedef Grid::ctype DF;
+ typedef double R;
+ typedef Dune::PDELab::OPBLocalFiniteElementMap<DF, R, 1, dim, Dune::GeometryType::simplex> DG1_FEM;
+ typedef Dune::PDELab::NoConstraints NoConstraintsAssembler;
+ typedef Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackend> DG1_DIRICHLET_GFS;
+ Dune::ParameterTree initree;
+ Dune::ParameterTreeParser::readINITree(argv[1], initree);
+ IniGridFactory<Grid> factory(initree);
+ std::shared_ptr<Grid> grid = factory.getGrid();
+ GV gv = grid->leafGridView();
+ DG1_FEM dg1_fem;
+ DG1_DIRICHLET_GFS dg1_dirichlet_gfs(gv, dg1_fem);
+ dg1_dirichlet_gfs.name("bla");
+ typedef Dune::SubsamplingVTKWriter<GV> VTKWriter;
+ int sublevel = initree.get<int>("vtk.subsamplinglevel", 0);
+ VTKWriter vtkwriter(gv, sublevel);
+ using LOP = Dune::PDELab::ConvectionDiffusionDG<CDProb<GV, R>, DG1_FEM>;
+ typedef DG1_DIRICHLET_GFS::ConstraintsContainer<R>::Type DG1_CC;
+ typedef Dune::PDELab::istl::BCRSMatrixBackend<> MatrixBackend;
+ typedef Dune::PDELab::GridOperator<DG1_DIRICHLET_GFS, DG1_DIRICHLET_GFS, LOP, MatrixBackend, DF, R, R, DG1_CC, DG1_CC> GO;
+ typedef GO::Traits::Domain V;
+ V x(dg1_dirichlet_gfs);
+ x = 0.0;
+ std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
+ typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LinearSolver;
+ typedef Dune::PDELab::StationaryLinearProblemSolver<GO, LinearSolver, V> SLP;
+ DG1_CC dg1_cc;
+ dg1_cc.clear();
+ CDProb<GV, R> params;
+ LOP lop(params, Dune::PDELab::ConvectionDiffusionDGMethod::SIPG);
+ int generic_dof_estimate = 6 * dg1_dirichlet_gfs.maxLocalSize();
+ int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
+ MatrixBackend mb(dofestimate);
+ GO go(dg1_dirichlet_gfs, dg1_cc, dg1_dirichlet_gfs, dg1_cc, lop, mb);
+ std::cout << "gfs with " << dg1_dirichlet_gfs.size() << " dofs generated "<< std::endl;
+ std::cout << "cc with " << dg1_cc.size() << " dofs generated "<< std::endl;
+ LinearSolver ls(false);
+ double reduction = initree.get<double>("reduction", 1e-12);
+ SLP slp(go, ls, x, reduction);
+ slp.apply();
+ Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_dirichlet_gfs, x);
+ vtkwriter.write(vtkfile, Dune::VTK::ascii);
+}
+
+#endif //_HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
diff --git a/test/poisson/reference_main.cc b/test/poisson/reference_main.cc
new file mode 100644
index 00000000..4304d9b0
--- /dev/null
+++ b/test/poisson/reference_main.cc
@@ -0,0 +1,33 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/exceptions.hh>
+
+#include "reference_driver.hh"
+
+int main(int argc, char** argv)
+{
+ try{
+ //Maybe initialize Mpi
+ Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
+ if(Dune::MPIHelper::isFake)
+ std::cout<< "This is a sequential program." << std::endl;
+ else
+ std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
+ <<" processes!"<<std::endl;
+
+ driver(argc, argv);
+
+ return 0;
+ }
+ catch (Dune::Exception &e){
+ std::cerr << "Dune reported error: " << e << std::endl;
+ return 1;
+ }
+ catch (...){
+ std::cerr << "Unknown exception thrown!" << std::endl;
+ return 1;
+ }
+}
--
GitLab