Commit 5838616d authored by oliver.sander_at_tu-dresden.de's avatar oliver.sander_at_tu-dresden.de
Browse files

Make examples build and run

parent 140820bc
......@@ -1353,7 +1353,7 @@ on the icon in the margin.
\marginpar{\attachfile[author=Carsten Gr\"aser and Oliver Sander,
color = 1 0 0,
mimetype=text/plain,
description=Complete source code of the Poisson example]{../../src/poisson.cc}}
description=Complete source code of the Poisson example]{../../examples/poisson.cc}}
\bigskip
......@@ -1365,7 +1365,7 @@ for this first example we discuss all five of them in turn.
The grid setup part looks like this:
%
\lstinputlisting[linerange={poisson_grid_setup_begin-poisson_grid_setup_end},
numbers=left]{../../src/poisson.cc}
numbers=left]{../../examples/poisson.cc}
%
The actual grid creation takes place in Lines~\ref{li:poisson_grid_creation_begin}--\ref{li:poisson_grid_creation_end}.
We use the \cpp{YaspGrid} grid implementation, which is a structured cube grid. We create a grid with $2 \times 2$
......@@ -1389,7 +1389,7 @@ The next part of the code creates the stiffness matrix $A$ and the load vector $
module.
%
\lstinputlisting[linerange={poisson_assembler_begin-poisson_assembler_end},
numbers=left]{../../src/poisson.cc}
numbers=left]{../../examples/poisson.cc}
%
Lines~\ref{li:poisson_fe_space_setup_begin}--\ref{li:poisson_fe_space_setup_end} set up a first-order Lagrangian function space
......@@ -1407,13 +1407,13 @@ Before we begin, we construct the initial iterate:
%
\lstinputlisting[linerange={initial_iterate_begin-initial_iterate_end},
numbers=left]
{../../src/poisson.cc}
{../../examples/poisson.cc}
%
The next block sets up the functional $J$. As it is simply quadratic, we can use the prebuilt
\cpp{QuadraticFunctional} class.
%
\lstinputlisting[linerange={poisson_create_functional_begin-poisson_create_functional_end},numbers=left]{../../src/poisson.cc}
\lstinputlisting[linerange={poisson_create_functional_begin-poisson_create_functional_end},numbers=left]{../../examples/poisson.cc}
%
First, we set up a coarse correction step. We choose a linear geometric $V(3,3)$ multigrid step.
......@@ -1428,11 +1428,13 @@ and independent from \modulename{dune-tnnmg}.
%
\lstinputlisting[linerange={poisson_setup_mg_begin-poisson_setup_mg_end},
numbers=left]
{../../src/poisson.cc}
{../../examples/poisson.cc}
%
The following code block then sets up an object representing one iteration of the TNNMG method.
%
\lstinputlisting[linerange={poisson_create_tnnmgstep_begin-poisson_create_tnnmgstep_end},numbers=left]{../../src/poisson.cc}
\lstinputlisting[linerange={poisson_create_tnnmgstep_begin-poisson_create_tnnmgstep_end},
numbers=left]
{../../examples/poisson.cc}
%
\todo[inline]{Erkläre den vorangegangenen Code!}
......@@ -1447,7 +1449,9 @@ method of the \cpp{TNNMGStep} object repeatedly until a prescribed error criteri
This work is conveniently done for us by the \cpp{LoopSolver} class from \modulename{dune-solvers}.
The energy norm object is used by the loop solver to measure the error.
%
\lstinputlisting[linerange={poisson_create_loopsolver_begin-poisson_create_loopsolver_end},numbers=left]{../../src/poisson.cc}
\lstinputlisting[linerange={poisson_create_loopsolver_begin-poisson_create_loopsolver_end},
numbers=left]
{../../examples/poisson.cc}
%
......@@ -1455,12 +1459,15 @@ The energy norm object is used by the loop solver to measure the error.
Once this is done it is easy to actually call the solver
%
\lstinputlisting[linerange={poisson_solve_begin-poisson_solve_end},
numbers=left]{../../src/poisson.cc}
numbers=left]
{../../examples/poisson.cc}
%
The solution is then contained in the variable \cpp{x}, of type \cpp{BlockVector<FieldVector<double,1> >}.
It is written into a vtu file with the last code block:
%
\lstinputlisting[linerange={poisson_write_begin-poisson_write_end},numbers=left]{../../src/poisson.cc}
\lstinputlisting[linerange={poisson_write_begin-poisson_write_end},
numbers=left]
{../../examples/poisson.cc}
%
When running the program, it produces the screen output:
......@@ -1607,13 +1614,16 @@ $\overline{u} \equiv 0.3$. The complete source code can be accessed through the
\marginpar{\attachfile[author=Carsten Gr\"aser and Oliver Sander,
color = 1 0 0,
mimetype=text/plain,
description=Complete source code of the Poisson example]{../../src/obstacle.cc}}
description=Complete source code of the Poisson example]
{../../examples/obstacle.cc}}
\bigskip
Remember the code snippet from the Poisson example that created the functional $J$
%
\lstinputlisting[linerange={poisson_create_functional_begin-poisson_create_functional_end},numbers=left]{../../src/poisson.cc}
\lstinputlisting[linerange={poisson_create_functional_begin-poisson_create_functional_end},
numbers=left]
{../../examples/poisson.cc}
%
Note how the class \lstinline!ZeroNonlinearity! is used to implement the zero functions $\phi_i \equiv 0$.
All we need to do is to replace this class by another that implements the obstacle functional $\chi$ for us.
......@@ -1630,7 +1640,9 @@ and implements the functional
In other words it implements an upper and a lower obstacle. If you do not want the lower obstacle just set all the entries to $-\infty$.
This is what the following code does.
%
\lstinputlisting[linerange={obstacle_create_functional_begin-obstacle_create_functional_end},numbers=left]{../../src/obstacle.cc}
\lstinputlisting[linerange={obstacle_create_functional_begin-obstacle_create_functional_end},
numbers=left]
{../../examples/obstacle.cc}
%
......@@ -1843,7 +1855,8 @@ The complete source code of this example can again be accessed through the icon
\marginpar{\attachfile[author=Carsten Gr\"aser and Oliver Sander,
color = 1 0 0,
mimetype=text/plain,
description=Complete source code of the one-body contact example]{../../src/contact.cc}}
description=Complete source code of the one-body contact example]
{../../examples/contact.cc}}
%
In our code, we again use a uniform grid with quadrilateral elements to discretize the domain $\Omega$.
This grid is implemented as an object of type \lstinline!YaspGrid!, and grid creation is identical
......@@ -1856,7 +1869,8 @@ we want to prescribe the contact conditions. In our example this is $\Gamma_C =
code snippet extracts both, and stores the information in two \lstinline!BitSetVector!s \lstinline!dirichletNodes! and
\lstinline!contactNodes!, respectively.
%
\lstinputlisting[linerange={contact_contactnodes_begin-contact_contactnodes_end}]{../../src/contact.cc}
\lstinputlisting[linerange={contact_contactnodes_begin-contact_contactnodes_end}]
{../../examples/contact.cc}
%
Note how the Dirichlet nodes get stored in a \lstinline!BitSetVector! with block size \lstinline!d!, whereas the
\lstinline!BitSetVector! for the contact nodes has only block size 1. Since we have a vector-valued problem, the property
......@@ -1875,7 +1889,9 @@ differs from the previous examples because we are discretizing a different diffe
but there are also some technicalities related to the fact that we work with a vector-valued
space this time. The first step is again to create a function space object.
%
\lstinputlisting[linerange={contact_fe_space_setup_begin-contact_fe_space_setup_end},numbers=left]{../../src/contact.cc}
\lstinputlisting[linerange={contact_fe_space_setup_begin-contact_fe_space_setup_end},
numbers=left]
{../../examples/contact.cc}
%
The \lstinline!FiniteElementMap! object created in Line~\ref{li:contact_create_fe_map} specifies
that we want to work with first-order Lagrangian finite elements. The \lstinline!VectorGridFunctionSpace!
......@@ -1893,13 +1909,16 @@ Next we construct the differential operator, i.e., the object that assembles the
For this we need to combine the FE space we just constructed with an assembler the computes the linear
elasticity stiffness matrix on individual elements. This is done by the following code:
%
\lstinputlisting[linerange={contact_operator_setup_begin-contact_operator_setup_end},numbers=left]{../../src/contact.cc}
\lstinputlisting[linerange={contact_operator_setup_begin-contact_operator_setup_end},
numbers=left]
{../../examples/contact.cc}
%
Note how the material parameters, represented in form of the two Lam\'e constants $\lambda$ and $\mu$,
are handed to a parameter class \lstinline!StVenantKirchhoffParameters!. This class is defined
separately at the top of the main source file. It looks like this:
%
\lstinputlisting[linerange={contact_material_model_begin-contact_material_model_end}]{../../src/contact.cc}
\lstinputlisting[linerange={contact_material_model_begin-contact_material_model_end}]
{../../examples/contact.cc}
%
In the simple
form shown here, it appears as an overly complicated way to store two floating point numbers. However, in more general cases,
......@@ -1907,7 +1926,8 @@ where $\lambda$ and $\mu$ can depend on $x$, this construction becomes useful.
Finally, we can assemble the stiffness matrix. This happens as in the previous examples
%
\lstinputlisting[linerange={contact_assembly_begin-contact_assembly_end}]{../../src/contact.cc}
\lstinputlisting[linerange={contact_assembly_begin-contact_assembly_end}]
{../../examples/contact.cc}
%
Note how we do not have to assemble an algebraic load vector \lstinline!rhs!: as we have neither a
volume source term nor a Neumann boundary term in this example, the linear term \lstinline!rhs! is zero
......@@ -1918,7 +1938,8 @@ anyway.
Next we set up the initial iterate. As previously, we start from zero; however, as we now have non-zero
Dirichlet values we need to add them at the correct positions.
%
\lstinputlisting[linerange={contact_initial_iterate_begin-contact_initial_iterate_end}]{../../src/contact.cc}
\lstinputlisting[linerange={contact_initial_iterate_begin-contact_initial_iterate_end}]
{../../examples/contact.cc}
%
\subsubsection{Creating the TNNMG solver step}
......@@ -1931,7 +1952,8 @@ differs from the previous obstacle problem example code. Indeed, the algebraic
set~\eqref{eq:contact_algebraic_admissible_set} is again an axis-aligned hypercube. Hence we can
again implement its indicator functional by the \lstinline!CubeIndicatorFunctional! class.
%
\lstinputlisting[linerange={contact_create_tnnmg_begin-contact_create_tnnmg_end}]{../../src/contact.cc}
\lstinputlisting[linerange={contact_create_tnnmg_begin-contact_create_tnnmg_end}]
{../../examples/contact.cc}
\subsection{One-Body Contact with Tresca Friction}
\label{sec:tresca_friction}
......@@ -2100,7 +2122,7 @@ the source code of the complete example by clicking on the icon in the margin.
color = 1 0 0,
mimetype=text/plain,
description=Complete source code of the Tresca friction example]
{../../src/trescafriction.cc}}
{../../examples/trescafriction.cc}}
The Tresca friction example consists of more code than previous examples, because we
can use less ready-made code from the \modulename{dune-tnnmg} module. On the other hand,
......@@ -2128,7 +2150,7 @@ This is done by the following piece of code:
%
\lstinputlisting[linerange={trescafriction_integrate_over_shapefunctions_begin-trescafriction_integrate_over_shapefunctions_end},
numbers=left]
{../../src/trescafriction.cc}
{../../examples/trescafriction.cc}
%
\todo[inline]{Achtung: aktuell wird hier noch über $\Omega$ integriert, und nicht über $\Gamma_C$.
Das muss noch korrigiert werden.}
......@@ -2137,7 +2159,7 @@ Setting up the functional is done by the following code.
\lstinputlisting[linerange={trescafriction_create_functional_begin-trescafriction_create_functional_end},
numbers=left]
{../../src/trescafriction.cc}
{../../examples/trescafriction.cc}
Note how the fact that the non-quadratic part of the functional is a sum $\chi_K + j_T$ is reflected in the code.
Lines~\ref{li:trescafriction_indicator_functional_begin}--\ref{li:trescafriction_indicator_functional_end} set up a \lstinline!CubeIndicatorFunctional! for the contact constraint
......@@ -2160,7 +2182,7 @@ Here, we change that to something taylor-made, which we call \lstinline!TrescaFr
We discuss the \lstinline!TrescaFrictionTNNMGStepImp! in more detail in Section~\ref{}.
\lstinputlisting[linerange={trescafriction_create_tnnmgstep_begin-trescafriction_create_tnnmgstep_end}]
{../../src/trescafriction.cc}
{../../examples/trescafriction.cc}
\subsubsection{The nonlinear functional}
\label{sec:tresca_functional}
......@@ -2172,7 +2194,7 @@ in the margin
color = 1 0 0,
mimetype=text/plain,
description=Complete source code of the TrescaFrictionFunctional class]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}}
%
This file is not part of the \modulename{dune-tnnmg} module, but rather belongs to the examples
distributed with this article.
......@@ -2187,7 +2209,7 @@ The class has a single template parameter \lstinline!dim!, which is the dimensio
It inherits from the abstract base class \lstinline!NonsmoothConvexFunctional!
%
\lstinputlisting[linerange={trescafunctional_header_begin-trescafunctional_header_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
The first method evaluates the functional $j_T(v)$ for a given $v \in (\R^d)^n$.
It is not actually used anywhere in the solver, but can be helpful for debugging.
......@@ -2195,14 +2217,14 @@ It is not actually used anywhere in the solver, but can be helpful for debugging
\todo[inline]{Im Code fehlt das Integral über die Formfunktionen}
%
\lstinputlisting[linerange={trescafunctional_energy_begin-trescafunctional_energy_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
The code first defines a small helper function \lstinline!tangentialPart!, which computes
$v_T \in \R^{d-1}$ from $v \in \R^d$.
%
\lstinputlisting[linerange={trescafunctional_tangentialpart_begin-trescafunctional_tangentialpart_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
The next two methods compute the gradient and the Hesse matrix of the Tresca friction functional.
......@@ -2237,11 +2259,11 @@ the content of the \lstinline!gradient! and \lstinline!hessian! variables, respe
the calling algorithm to save one copy operation.
%
\lstinputlisting[linerange={trescafunctional_addgradient_begin-trescafunctional_addgradient_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
%
\lstinputlisting[linerange={trescafunctional_addhessian_begin-trescafunctional_addhessian_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
\todo[inline]{Was macht die leere Methode addHessianIndices?}
......@@ -2257,7 +2279,7 @@ the direction $u + tv$. For the Tresca friction functional $j_T$ this subdiffer
This general method is used by the line search algorithm.
%
\lstinputlisting[linerange={trescafunctional_directional_subdiff_begin-trescafunctional_directional_subdiff_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
In addition, the Gau\ss--Seidel-type algorithm we use to solve the local block problems
......@@ -2269,7 +2291,7 @@ if $v= e_{i,j}$. The corresponding code is:
%
\lstinputlisting[linerange={trescafunctional_subdiff_begin-trescafunctional_subdiff_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
Note that there is no argument $u$.
\todo[inline]{Discuss internal state variable}
......@@ -2282,7 +2304,7 @@ In the case of the Tresca friction functional, which is defined everywhere on $(
the method simply returns $(-inf, inf)$ regardless of the input arguments.
%
\lstinputlisting[linerange={trescafunctional_domain_begin-trescafunctional_domain_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
The second method is used to truncate away components where the functional has very large
......@@ -2291,21 +2313,19 @@ where $j_T$ is not differentiable.
\todo[inline]{Begründung?}
%
\lstinputlisting[linerange={trescafunctional_regularity_begin-trescafunctional_regularity_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
\todo[inline]{Methods for the internal state}
%
\lstinputlisting[linerange={trescafunctional_state_begin-trescafunctional_state_end}]
{../../dune/tnnmg-examples/trescafrictionfunctional.hh}
{../../dune/tnnmg/functionals/trescafrictionfunctional.hh}
%
\subsubsection{The local block solver}
\label{sec:tresca_local_block_solver}
%\lstinputlisting{../../dune/tnnmg-examples/trescafrictiontnnmgstepimp.hh}
\bibliography{dune-tnnmg}
......
......@@ -8,3 +8,6 @@ Version: 2.8-git
Maintainer: graeser@mi.fu-berlin.de
#depending on
Depends: dune-common (>=2.7) dune-istl (>=2.7) dune-solvers (>=2.7)
# Optional dependencies. These are needed for the examples.
Suggests: dune-pdelab
#include <config.h>
#include <array>
#include <dune/common/parametertree.hh>
#include <dune/grid/yaspgrid.hh>
......@@ -97,7 +99,7 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
FieldVector<double,dim> h(1); // bounding box
array<int,dim> n;
std::array<int,dim> n;
std::fill(n.begin(), n.end(), 10);
GridType grid(h,n);
......@@ -132,8 +134,8 @@ int main (int argc, char *argv[]) try
typedef PDELab::VectorGridFunctionSpace<GridView, /*@\label{li:contact_create_gfs_begin}@*/
FEM,
dim, /*@\label{li:contact_value_space_dimension}@*/
PDELab::istl::VectorBackend<PDELab::istl::Blocking::fixed>,
PDELab::istl::VectorBackend<>,
PDELab::ISTL::VectorBackend<PDELab::ISTL::Blocking::fixed>,
PDELab::ISTL::VectorBackend<>,
PDELab::ConformingDirichletConstraints,
PDELab::EntityBlockedOrderingTag,
PDELab::DefaultLeafOrderingTag
......@@ -153,7 +155,7 @@ int main (int argc, char *argv[]) try
typedef PDELab::GridOperator<GFS,
GFS,
PDELab::LinearElasticity<StVenantKirchhoffParameters<GridView> >,
PDELab::istl::BCRSMatrixBackend<>,
PDELab::ISTL::BCRSMatrixBackend<>,
double,double,double> GridOperator;
GridOperator gridOperator(gfs, gfs, lop, {9});
......@@ -222,7 +224,7 @@ int main (int argc, char *argv[]) try
auto linearMultigridStep = std::make_shared<MultigridStep<MatrixType, VectorType> >();
linearMultigridStep->setMGType(1, nu1, nu2);
linearMultigridStep->basesolver_ = linearBaseSolver.get();
linearMultigridStep->setBaseSolver(linearBaseSolver);
linearMultigridStep->setSmoother(linearSmoother);
using TransferOperator = TruncatedDenseMGTransfer<VectorType>;
......@@ -235,7 +237,7 @@ int main (int argc, char *argv[]) try
// create transfer operator from level i to i+1
transferOperators[i] = std::make_shared<TransferOperator>();
transferOperators[i]->setup(grid, i, i+1);
linearMultigridStep->mgTransfer_[i] = transferOperators[i].get();
linearMultigridStep->mgTransfer_[i] = transferOperators[i];
}
// { contact_setup_mg_end }
......@@ -266,17 +268,17 @@ int main (int argc, char *argv[]) try
DefectProjection,
LineSearchSolver>;
auto step = Step(J, x, nonlinearSmoother, linearMultigridStep, 1, DefectProjection(), LineSearchSolver());
auto step = std::make_shared<Step>(J, x, nonlinearSmoother, linearMultigridStep, 1, DefectProjection(), LineSearchSolver());
auto norm = EnergyNorm<MatrixType, VectorType>(stiffnessMatrix);
::LoopSolver<VectorType> solver(&step, numIt, tolerance, &norm, Solver::FULL);
auto norm = std::make_shared<EnergyNorm<MatrixType, VectorType> >(stiffnessMatrix);
::LoopSolver<VectorType> solver(step, numIt, tolerance, norm, Solver::FULL);
step.setIgnore(dirichletNodes);
step.setPreSmoothingSteps(3);
step->setIgnore(dirichletNodes);
step->setPreSmoothingSteps(3);
solver.addCriterion(
[&](){
return Dune::formatString(" % 12d", step.linearization().truncated().count());
return Dune::formatString(" % 12d", step->linearization().truncated().count());
},
" truncated ");
......@@ -293,7 +295,7 @@ int main (int argc, char *argv[]) try
solver.addCriterion(
[&](){
return Dune::formatString(" % 12.5e", step.lastDampingFactor());
return Dune::formatString(" % 12.5e", step->lastDampingFactor());
},
" step size ");
......@@ -322,6 +324,6 @@ int main (int argc, char *argv[]) try
Dune::PDELab::addSolutionToVTKWriter(vtkWriter,gfs,xV);
vtkWriter.write("contact_result");
} catch (Exception e) {
std::cout << e << std::endl;
}
} catch (Exception& e) {
std::cout << e.what() << std::endl;
}
......@@ -118,14 +118,19 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////////////
// Assemble the stiffness matrix and load vector
////////////////////////////////////////////////////////////////
using Basis = Functions::PQkNodalBasis<GridView,1>;
using Basis = Functions::LagrangeBasis<GridView,1>;
auto basis = std::make_shared<Basis>(gridView);
using VectorType = BlockVector<FieldVector<double,1> >;
typedef PDELab::Experimental::GridFunctionSpace<Basis,
VectorType,
PDELab::NoConstraints> GridFunctionSpace;
GridFunctionSpace gfs(basis);
using FiniteElementMap = PDELab::QkLocalFiniteElementMap<GridView,double,double,1>;
FiniteElementMap fem(gridView);
typedef PDELab::GridFunctionSpace<GridView,
FiniteElementMap,
PDELab::ConformingDirichletConstraints,
PDELab::ISTL::VectorBackend<>
> GridFunctionSpace;
GridFunctionSpace gfs(grid.leafGridView(), fem);
gfs.name("value");
// Container for the Dirichlet boundary conditions
......@@ -144,7 +149,7 @@ int main (int argc, char *argv[]) try
typedef PDELab::GridOperator<GridFunctionSpace,
GridFunctionSpace,
LOP,
PDELab::istl::BCRSMatrixBackend<>,
PDELab::ISTL::BCRSMatrixBackend<>,
double,double,double,C,C> GO;
GO go(gfs,cg,gfs,cg,lop, {9});
......@@ -182,19 +187,19 @@ int main (int argc, char *argv[]) try
&baseEnergyNorm,
Solver::QUIET);
auto smoother = TruncatedBlockGSStep<MatrixType, VectorType, BitVector>();
auto smoother = std::make_shared<TruncatedBlockGSStep<MatrixType, VectorType, BitVector> >();
auto linearMultigridStep = std::make_shared<Solvers::MultigridStep<MatrixType, VectorType> >();
linearMultigridStep->setMGType(1, nu1, nu2);
linearMultigridStep->basesolver_ = linearBaseSolver.get();
linearMultigridStep->setSmoother(&smoother);
linearMultigridStep->setBaseSolver(linearBaseSolver);
linearMultigridStep->setSmoother(smoother);
linearMultigridStep->mgTransfer_.resize(grid.maxLevel());
for (size_t i=0; i<linearMultigridStep->mgTransfer_.size(); i++)
{
// create transfer operator from level i to i+1
auto transferOperator = new TruncatedDenseMGTransfer<VectorType>;
auto transferOperator = std::make_shared<TruncatedDenseMGTransfer<VectorType> >();
transferOperator->setup(grid, i, i+1);
linearMultigridStep->mgTransfer_[i] = transferOperator;
......@@ -224,17 +229,17 @@ int main (int argc, char *argv[]) try
DefectProjection,
LineSearchSolver>;
auto step = Step(J, x, nonlinearSmoother, linearMultigridStep, mu, DefectProjection(), LineSearchSolver());
auto step = std::make_shared<Step>(J, x, nonlinearSmoother, linearMultigridStep, mu, DefectProjection(), LineSearchSolver());
auto norm = EnergyNorm<MatrixType, VectorType>(stiffnessMatrix);
::LoopSolver<VectorType> solver(&step, numIt, tolerance, &norm, Solver::FULL);
auto norm = std::make_shared<EnergyNorm<MatrixType, VectorType> >(stiffnessMatrix);
::LoopSolver<VectorType> solver(step, numIt, tolerance, norm, Solver::FULL);
step.setIgnore(dirichletNodes);
step.setPreSmoothingSteps(3);
step->setIgnore(dirichletNodes);
step->setPreSmoothingSteps(3);
solver.addCriterion(
[&](){
return Dune::formatString(" % 12d", step.linearization().truncated().count());
return Dune::formatString(" % 12d", step->linearization().truncated().count());
},
" truncated ");
......@@ -251,7 +256,7 @@ int main (int argc, char *argv[]) try
solver.addCriterion(
[&](){
return Dune::formatString(" % 12.5e", step.lastDampingFactor());
return Dune::formatString(" % 12.5e", step->lastDampingFactor());
},
" step size ");
......@@ -264,7 +269,7 @@ int main (int argc, char *argv[]) try
std::vector<double> correctionNorms;
solver.addCriterion(Dune::Solvers::correctionNormCriterion(step, norm, 1e-10, correctionNorms));
solver.addCriterion(Dune::Solvers::correctionNormCriterion(*step, *norm, 1e-10, correctionNorms));
// { obstacle_create_functional_end }
/////////////////////////////
......@@ -285,6 +290,6 @@ int main (int argc, char *argv[]) try
vtkWriter.addVertexData(xFunction, VTK::FieldInfo("x", VTK::FieldInfo::Type::scalar, 1));
vtkWriter.write("obstacle_result");
} catch (Exception e) {
std::cout << e << std::endl;
}
} catch (Exception& e) {
std::cout << e.what() << std::endl;
}
......@@ -20,7 +20,7 @@
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include <dune/solvers/iterationsteps/blockgsstep.hh>
#include <dune/solvers/iterationsteps/blockgssteps.hh>
#include <dune/solvers/transferoperators/truncateddensemgtransfer.hh>
#include <dune/tnnmg/functionals/quadraticfunctional.hh>
......@@ -96,22 +96,25 @@ int main (int argc, char *argv[]) try
// { poisson_assembler_begin }
// Construct first-order Lagrangian finite element space for a cube grid
using Basis = Functions::PQkNodalBasis<GridView,1>;
using Basis = Functions::LagrangeBasis<GridView,1>;
auto basis = std::make_shared<Basis>(gridView);
using VectorType = BlockVector<FieldVector<double,1> >;
typedef PDELab::Experimental::GridFunctionSpace<Basis,
VectorType,
PDELab::NoConstraints> GridFunctionSpace;
GridFunctionSpace gfs(basis);
using FiniteElementMap = PDELab::QkLocalFiniteElementMap<GridView,double,double,1>;
FiniteElementMap fem(gridView);
typedef PDELab::GridFunctionSpace<GridView,
FiniteElementMap,
PDELab::ConformingDirichletConstraints,
PDELab::ISTL::VectorBackend<>
> GridFunctionSpace;
GridFunctionSpace gfs(grid.leafGridView(), fem);
gfs.name("value"); /*@\label{li:poisson_fe_space_setup_end}@*/
// make grid operator
PoissonProblem<GridView,double> poissonProblem; /*@\label{li:poisson_assembler_setup_begin}@*/
typedef PDELab::ConvectionDiffusionFEM<decltype(poissonProblem),
typename GridFunctionSpace::Traits::FiniteElementMap> LocalOperator;
FiniteElementMap> LocalOperator;
LocalOperator localOperator(poissonProblem);
typedef PDELab::GridOperator<GridFunctionSpace,
......@@ -122,7 +125,9 @@ int main (int argc, char *argv[]) try
GlobalLaplaceOperator globalLaplaceOperator(gfs,gfs,localOperator, {9}); /*@\label{li:poisson_assembler_setup_end}@*/
// make coefficent vector and initialize it from a function
typedef GlobalLaplaceOperator::Traits::Domain VectorContainer; /*@\label{li:poisson_assembly_begin}@*/
// typedef GlobalLaplaceOperator::Traits::Domain VectorContainer;
using VectorContainer = PDELab::Backend::Vector<GridFunctionSpace,double>; /*@\label{li:poisson_assembly_begin}@*/
VectorContainer x0(gfs);
x0 = 0.0;
......@@ -163,7 +168,12 @@ int main (int argc, char *argv[]) try
// { poisson_create_functional_end }
// { poisson_setup_mg_begin }
BlockGSStep<MatrixType,VectorType> linearBaseSolverStep;
auto localLinearSolver = Solvers::BlockGS::LocalSolvers::direct();
using BlockGSStep = Dune::Solvers::BlockGSStep<std::decay_t<decltype(localLinearSolver)>,
MatrixType,
VectorType,
BitVector>;
BlockGSStep linearBaseSolverStep;
EnergyNorm<MatrixType,VectorType> baseEnergyNorm(linearBaseSolverStep);
auto linearBaseSolver = std::make_shared<Solvers::LoopSolver<VectorType> >(&linearBaseSolverStep,
......@@ -172,19 +182,19 @@ int main (int argc, char *argv[]) try
&baseEnergyNorm,
Solver::QUIET);
auto smoother = BlockGSStep<MatrixType, VectorType, BitVector>();
auto smoother = std::make_shared<BlockGSStep>();
auto linearMultigridStep = std::make_shared<Solvers::MultigridStep<MatrixType, VectorType> >();
linearMultigridStep->setMGType(1, nu1, nu2);
linearMultigridStep->basesolver_ = linearBaseSolver.get();
linearMultigridStep->setSmoother(&smoother);
linearMultigridStep->setBaseSolver(linearBaseSolver);
linearMultigridStep->setSmoother(smoother);
linearMultigridStep->mgTransfer_.resize(grid.maxLevel());
for (size_t i=0; i<linearMultigridStep->mgTransfer_.size(); i++)
{
// create transfer operator from level i to i+1
auto transferOperator = new TruncatedDenseMGTransfer<VectorType>;
auto transferOperator = std::make_shared<TruncatedDenseMGTransfer<VectorType> >();
transferOperator->setup(grid, i, i+1);
linearMultigridStep->mgTransfer_[i] = transferOperator;
......@@ -206,20 +216,20 @@ int main (int argc, char *argv[]) try