1#ifndef __ELLIPTIC_FEMSCHEME_HH__
2#define __ELLIPTIC_FEMSCHEME_HH__
8#include <dune/fem/gridpart/adaptiveleafgridpart.hh>
11#include <dune/fem/space/lagrange.hh>
14#include <dune/fem/function/adaptivefunction.hh>
15#include <dune/fem/space/common/adaptationmanager.hh>
18#include <dune/fem/function/blockvectorfunction.hh>
21#include <dune/fem/solver/newtoninverseoperator.hh>
24#include <dune/fem/space/common/interpolate.hh>
27#include <dune/fem/misc/l2norm.hh>
28#include <dune/fem/misc/h1norm.hh>
31#include <dune/fem/io/parameter.hh>
32#include <dune/fem/solver/parameter.hh>
35#include <dune/fem/io/file/dataoutput.hh>
38#include "../common/matrixhelper.hh"
39#include "../algorithms/femschemeinterface.hh"
40#include "../algorithms/operatorselector.hh"
41#include "../models/modeltraits.hh"
42#include "../operators/ellipticoperator.hh"
43#include "../operators/functionals/functionals.hh"
44#include "../operators/constraints/dirichletconstraints.hh"
45#include "../functions/functions.hh"
47#include "../estimators/ellipticestimator.hh"
48#include "../algorithms/marking.hh"
49#include "../common/dataoutput.hh"
52#include "../common/solverselector.hh"
90 template<
class DiscreteFunction,
class Model,
class InitialGuess,
class RHSFunctional>
98 typedef RHSFunctional RHSFunctionalType;
101 typedef typename DiscreteFunctionType::GridType
GridType;
133 typedef typename SolverSelectorType::LinearOperatorType LinearOperatorType;
134 typedef typename SolverSelectorType::LinearInverseOperatorType LinearInverseOperatorType;
135 typedef typename LinearInverseOperatorType::SolverParameterType LinearSolverParameterType;
139 Fem::NewtonInverseOperator<LinearOperatorType, LinearInverseOperatorType>
154 typedef DifferentiableEllipticOperator<LinearOperatorType, DiscreteModelType, const ConstraintsOperatorType&>
OperatorType;
165 typedef std::tuple<DiscreteFunctionType *, InitialGuessType *>
IOTupleType;
189 const RHSFunctionalType& rhsFunctional,
190 const std::string&
name)
191 : grid_(solution.space().gridPart().grid()),
192 model_(
asExpression(discreteModel(model, solution.space()))),
194 gridPart_(solution.space().gridPart()),
195 discreteSpace_(solution.space()),
198 estimator_(model_, discreteSpace_),
199 markingStrategy_(gridPart_, estimator_, name_ +
".adaptation"),
203 adaptationManager_(0),
205 constraints_(discreteSpace_, model_),
207 operator_(model_, constraints_),
209 linearOperator_(
"assembled elliptic operator", discreteSpace_, discreteSpace_),
210 solverParams_(TraitsType::isAffineLinear ? name_ +
".solver." : name_ +
".linear."),
213 initialGuess_(initialGuess),
215 rhsFunctional_(rhsFunctional),
217 ioTuple_(&solution_, &initialGuess_),
221 chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
230 : grid_(other.grid_),
231 model_(other.model_),
233 gridPart_(other.gridPart_),
234 discreteSpace_(other.discreteSpace_),
235 solution_(other.solution_),
237 estimator_(model_, discreteSpace_),
238 markingStrategy_(other.markingStrategy_),
242 adaptationManager_(0),
244 constraints_(discreteSpace_, model_),
246 operator_(model_, constraints_),
248 linearOperator_(
"assembled elliptic operator", discreteSpace_, discreteSpace_),
249 solverParams_(other.solverParams_),
252 initialGuess_(other.initialGuess_),
254 rhsFunctional_(other.rhsFunctional_),
256 ioTuple_(&solution_, &initialGuess_),
260 chooseDefaultSolverMethod<SolverSelectorType>(solverParams_);
268 if (adaptationManager_) {
269 delete adaptationManager_;
271 if (restrictProlong_) {
272 delete restrictProlong_;
276 virtual std::string
name ()
const
288 interpolate(initialGuess_, solution_);
292 virtual void assemble() {
293 operator_.jacobian(solution_, linearOperator_);
297 virtual void solve(
bool forceMatrixAssembling =
true)
300 constraints_.constrain(solution_);
308 rhs += rhsFunctional_;
309 constraints_.zeroConstrain(rhs);
311 if (TraitsType::isAffineLinear) {
323 assert(!TraitsType::isAffineLinear);
325 Dune::Fem::NewtonParameter newtonParam(solverParams_, name_ +
".newton.");
327 newton.bind(operator_);
329 newton(rhs, solution_);
331 converged_ = newton.converged();
332 assert(newton.converged());
341 assert(TraitsType::isAffineLinear);
344 if (sequence_ != discreteSpace_.sequence() || forceMatrixAssembling) {
350 sequence_ = discreteSpace_.sequence();
354 LinearInverseOperatorType solver(solverParams_);
355 solver.bind(linearOperator_);
366 operator_(solution_, update);
375 converged_ = (
static_cast<int>(solver.iterations()) < solverParams_.maxIterations());
382 virtual bool mark (
const double tolerance)
384 return markingStrategy_.mark(tolerance);
390 return estimator_.
estimate(solution_);
404 if (!adaptationManager_) {
413 adaptationManager_->adapt();
414 if (adaptationManager_->loadBalance()) {
430 if (!dataOutput_->willWrite()) {
435 dataOutput_->write();
437 return dataOutput_->writeStep() - 1;
453 rhs += rhsFunctional_;
454 constraints_.zeroConstrain(rhs);
457 operator_(solution_, update);
459 return std::sqrt(rhs.scalarProductDofs(rhs));
466 typedef Fem::H1Norm< GridPartType > NormType;
474 NormType norm(gridPart_, 2*discreteSpace_.order());
475 return norm.distance(initialGuess_, solution_);
481 size_t numberOfDofs = discreteSpace_.size() - discreteSpace_.slaveDofs().size() + 1;
483 numberOfDofs = grid_.comm().sum(numberOfDofs);
489 virtual void printSystemMatrix(
bool sparse =
false, std::ostream& out = std::clog) {
491 MatrixHelper<LinearOperatorType> mh(linearOperator_);
503 const std::string name_;
510 EstimatorType estimator_;
520 LinearOperatorType linearOperator_;
522 LinearSolverParameterType solverParams_;
523 bool converged_ =
false;
524 mutable int sequence_;
527 RHSFunctionalType rhsFunctional_;
538 class DiscreteFunction,
class Model,
class InitialGuess,
class RHSFunctional,
539 std::enable_if_t<(IsWrappableByConstLocalFunction<InitialGuess>::value
540 && IsLinearFunctional<RHSFunctional>::value
544 const InitialGuess& initialGuess,
545 const RHSFunctional& rhsFunctional,
546 const std::string name =
"acfem.schemes.elliptic")
550 return ReturnType(solution, model, initialGuess, rhsFunctional, name);
554 template<
class DiscreteFunction,
class Model,
class InitialGuess,
555 std::enable_if_t<IsWrappableByConstLocalFunction<InitialGuess>::value,
int> = 0>
559 const InitialGuess& initialGuess,
560 const std::string name =
"acfem.schemes.elliptic")
569 return ReturnType(solution,
572 zeroFunctional(solution.space()),
577 template<
class DiscreteFunction,
class Model,
class RHSFunctional,
578 std::enable_if_t<IsLinearFunctional<RHSFunctional>::value,
int> = 0>
581 const RHSFunctional& rhsFunctional,
582 const std::string name =
"acfem.schemes.elliptic")
587 ZeroGridFunction<
typename DiscreteFunction::FunctionSpaceType,
588 typename DiscreteFunction::GridPartType>,
592 return ReturnType(solution,
594 GridFunction::zeros(solution.space()),
600 template<
class DiscreteFunction,
class Model>
603 const std::string name =
"acfem.schemes.elliptic")
608 ZeroGridFunction<
typename DiscreteFunction::FunctionSpaceType,
609 typename DiscreteFunction::GridPartType>,
613 return ReturnType(solution,
615 GridFunction::zeros(solution.space()),
616 zeroFunctional(solution.space()),
Abstract space adaptative FEM scheme.
Definition: femschemeinterface.hh:70
Potentially overwrite some parameters.
Definition: dataoutput.hh:37
RangeFieldType estimate(const DiscreteFunctionType &uh)
calculate estimator
Definition: ellipticestimator.hh:259
Adaptive fem-scheme for "elliptic" problems.
Definition: ellipticfemscheme.hh:93
virtual void linearSolve(DiscreteFunctionType &rhs, bool forceMatrixAssembling)
Perform only one step of the Newton scheme for the affine-linear case.
Definition: ellipticfemscheme.hh:339
SolverSelector< DiscreteFunctionType, DiscreteModelType > SolverSelectorType
choose type of discrete function, Matrix implementation and solver implementation
Definition: ellipticfemscheme.hh:132
Fem::RestrictProlongDefault< DiscreteFunctionType > RestrictionProlongationType
type of restriction/prolongation projection for adaptive simulations (use default here,...
Definition: ellipticfemscheme.hh:144
InitialGuess InitialGuessType
Initial guess/exact solution.
Definition: ellipticfemscheme.hh:162
virtual size_t size() const
return some measure about the number of DOFs in use
Definition: ellipticfemscheme.hh:478
virtual double error() const
calculate L2/H1 error
Definition: ellipticfemscheme.hh:463
virtual void solve(bool forceMatrixAssembling=true)
Solve the system.
Definition: ellipticfemscheme.hh:297
virtual bool converged() const
check whether solver has converged
Definition: ellipticfemscheme.hh:441
virtual double residual() const
calculate residual (in small l^2)
Definition: ellipticfemscheme.hh:447
virtual double estimate()
calculate error estimator
Definition: ellipticfemscheme.hh:388
virtual int output()
data I/O
Definition: ellipticfemscheme.hh:420
MarkingStrategy< GridPartType > MarkingStrategyType
type of marking strategy
Definition: ellipticfemscheme.hh:159
DirichletConstraints< DiscreteFunctionSpaceType, DiscreteModelType > ConstraintsOperatorType
type of constraints operator
Definition: ellipticfemscheme.hh:150
DiscreteFunctionType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType
type of the discrete function space
Definition: ellipticfemscheme.hh:107
EllipticFemScheme(const EllipticFemScheme &other)
Copy-constructor.
Definition: ellipticfemscheme.hh:229
DiscreteFunctionSpaceType::FunctionSpaceType FunctionSpaceType
type of function space
Definition: ellipticfemscheme.hh:129
virtual void nonLinearSolve(DiscreteFunctionType &rhs)
Run the full Newton-scheme ...
Definition: ellipticfemscheme.hh:321
DiscreteFunctionType::GridType GridType
type of hierarchic grid
Definition: ellipticfemscheme.hh:101
virtual bool mark(const double tolerance)
mark elements for adaptation
Definition: ellipticfemscheme.hh:382
DiscreteFunctionType::GridPartType GridPartType
type of the grid view
Definition: ellipticfemscheme.hh:104
Fem::AdaptationManager< GridType, RestrictionProlongationType > AdaptationManagerType
type of adaptation manager handling adapation and DoF compression
Definition: ellipticfemscheme.hh:147
virtual void adapt()
do the adaptation for a given marking
Definition: ellipticfemscheme.hh:394
virtual std::string name() const
name of the Fem Scheme
Definition: ellipticfemscheme.hh:276
DiscreteFunction DiscreteFunctionType
Type of the discrete solution function.
Definition: ellipticfemscheme.hh:96
DifferentiableEllipticOperator< LinearOperatorType, DiscreteModelType, const ConstraintsOperatorType & > OperatorType
type of error estimator define differential operator
Definition: ellipticfemscheme.hh:154
virtual void initialize()
initialize solution
Definition: ellipticfemscheme.hh:282
Fem::NewtonInverseOperator< LinearOperatorType, LinearInverseOperatorType > NonLinearInverseOperatorType
Non-linear solver.
Definition: ellipticfemscheme.hh:140
DataOutput< GridType, IOTupleType > DataOutputType
type of data writer
Definition: ellipticfemscheme.hh:168
EllipticFemScheme(DiscreteFunctionType &solution, const ModelType &model, const InitialGuessType &initialGuess, const RHSFunctionalType &rhsFunctional, const std::string &name)
Constructor for the elliptic fem-scheme.
Definition: ellipticfemscheme.hh:186
Model ModelType
type of the provided model
Definition: ellipticfemscheme.hh:110
std::tuple< DiscreteFunctionType *, InitialGuessType * > IOTupleType
type of input/output tuple
Definition: ellipticfemscheme.hh:165
The zero functional.
Definition: zero.hh:24
DiscreteModel< Model, DiscreteFunctionSpaceType > DiscreteModelType
In the DG-case the resulting ModelType will a NitscheDirichletBoundaryModel.
Definition: ellipticfemscheme.hh:121
constexpr decltype(auto) asExpression(T &&t)
Return a non-closure expression as is.
Definition: interface.hh:122
auto ellipticFemScheme(DiscreteFunction &solution, const Model &model, const InitialGuess &initialGuess, const RHSFunctional &rhsFunctional, const std::string name="acfem.schemes.elliptic")
Adaptive fem-scheme for "elliptic" problems.
Definition: ellipticfemscheme.hh:542
ModelIntrospection::Traits< Model > ModelTraits
Traits class for models.
Definition: modeltraits.hh:898
Default expression traits definition is a recursion in order to ease disambiguation.
Definition: expressiontraits.hh:54
Select one appropriate (linear) solver depending on whether the model is symmetric and/or semidefinit...
Definition: solverselector.hh:91