Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
/********************************************************/
// Beware of line number changes, they may corrupt docu!
//! \brief Driver function to set up and solve the problem
/********************************************************/
template<typename GV, typename FEM>
void driver (const GV& gv, const FEM& fem,
Dune::ParameterTree& ptree)
{
// dimension and important types
const int dim = GV::dimension;
typedef double RF; // type for computations
// make PDE parameter class
RF eta = ptree.get("problem.eta",(RF)1.0);
Problem<RF> problem(eta);
auto glambda = [&](const auto& e, const auto& x)
{return problem.g(e,x);};
auto g = Dune::PDELab::
makeGridFunctionFromCallable(gv,glambda);
auto blambda = [&](const auto& i, const auto& x)
{return problem.b(i,x);};
auto b = Dune::PDELab::
makeBoundaryConditionFromCallable(gv,blambda);
// Make grid function space
//== Exercise 2 {
// typedef Dune::PDELab::ConformingDirichletConstraints CON;
typedef Dune::PDELab::NoConstraints CON;
//== }
typedef Dune::PDELab::ISTL::VectorBackend<> VBE;
typedef Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE> GFS;
GFS gfs(gv,fem);
gfs.name("Vh");
// Assemble constraints
//== Exercise 2 {
// typedef typename GFS::template
// ConstraintsContainer<RF>::Type CC;
// CC cc;
// Dune::PDELab::constraints(b,gfs,cc); // assemble constraints
// std::cout << "constrained dofs=" << cc.size() << " of "
// << gfs.globalSize() << std::endl;
typedef Dune::PDELab::EmptyTransformation CC;
//== }
// A coefficient vector
using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
Z z(gfs); // initial value
// Make a grid function out of it
typedef Dune::PDELab::DiscreteGridFunction<GFS,Z> ZDGF;
ZDGF zdgf(gfs,z);
// Fill the coefficient vector
Dune::PDELab::interpolate(g,gfs,z);
// Make a local operator
//== Exercise 2 {
RF stab = ptree.get("fem.stab",(RF)1);
typedef NitscheNonlinearPoissonFEM<Problem<RF>,FEM> LOP;
LOP lop(problem,stab);
#else
typedef NonlinearPoissonFEM<Problem<RF>,FEM> LOP;
LOP lop(problem);
#endif
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
//== }
// Make a global operator
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE;
int degree = ptree.get("fem.degree",(int)1);
MBE mbe((int)pow(1+2*degree,dim));
typedef Dune::PDELab::GridOperator<
GFS,GFS, /* ansatz and test space */
LOP, /* local operator */
MBE, /* matrix backend */
RF,RF,RF, /* domain, range, jacobian field type*/
CC,CC /* constraints for ansatz and test space */
> GO;
//== Exercise 2 {
// GO go(gfs,cc,gfs,cc,lop,mbe);
GO go(gfs,gfs,lop,mbe);
//== }
// Select a linear solver backend
typedef Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GO> LS;
LS ls(100,2);
// set up nonlinear solver
Dune::PDELab::Newton<GO,LS,Z> newton(go,z,ls);
newton.setReassembleThreshold(0.0); // always reassemble J
newton.setVerbosityLevel(3); // be verbose
newton.setReduction(1e-10); // total reduction
newton.setMinLinearReduction(1e-4); // min. red. in lin. solve
newton.setMaxIterations(25); // limit number of its
newton.setLineSearchMaxIterations(10); // limit line search
// solve nonlinear problem
newton.apply();
// Write VTK output file
Z w(gfs); // Lagrange interpolation of exact solution
Dune::PDELab::interpolate(g,gfs,w);
ZDGF wdgf(gfs,w);
int subsampling = ptree.get("output.subsampling",(int)1);
Dune::SubsamplingVTKWriter<GV> vtkwriter(gv,Dune::refinementIntervals(subsampling));
typedef Dune::PDELab::VTKGridFunctionAdapter<ZDGF> VTKF;
vtkwriter.addVertexData(std::shared_ptr<VTKF>(new
VTKF(zdgf,"fesol")));
vtkwriter.addVertexData(std::shared_ptr<VTKF>(new
VTKF(wdgf,"exact")));
vtkwriter.write(ptree.get("output.filename","output"),
Dune::VTK::appendedraw);
}