stokes0.cc 1.54 KB
Newer Older
1
#include <amdis/AMDiS.hpp>
2
#include <amdis/LocalOperators.hpp>
3
#include <amdis/ProblemStat.hpp>
4
5
6

using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
7
using Grid = Dune::YaspGrid<2>;
8
using StokesParam = TaylorHoodBasis<Grid>;
9
10
11
12
using StokesProblem = ProblemStat<StokesParam>;

int main(int argc, char** argv)
{
13
  Environment env(argc, argv);
14

Praetorius, Simon's avatar
Praetorius, Simon committed
15
16
  StokesProblem prob("stokes");
  prob.initialize(INIT_ALL);
17

Praetorius, Simon's avatar
Praetorius, Simon committed
18
19
  double viscosity = 1.0;
  Parameters::get("stokes->viscosity", viscosity);
20

Praetorius, Simon's avatar
Praetorius, Simon committed
21
22
23
  // tree-paths for components
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
24

Praetorius, Simon's avatar
Praetorius, Simon committed
25
26
  for (std::size_t i = 0; i < WORLDDIM; ++i) {
    // <viscosity*grad(u_i), grad(v_i)>
Praetorius, Simon's avatar
Praetorius, Simon committed
27
28
    auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
    prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
29

Praetorius, Simon's avatar
Praetorius, Simon committed
30
31
32
    // <d_i(v_i), p>
    auto opP = makeOperator(tag::partialtest_trial{i}, 1.0);
    prob.addMatrixOperator(opP, treepath(_v,i), _p);
33

Praetorius, Simon's avatar
Praetorius, Simon committed
34
35
36
37
    // <q, d_i(u_i)>
    auto opDiv = makeOperator(tag::test_partialtrial{i}, 1.0);
    prob.addMatrixOperator(opDiv, _p, treepath(_v,i));
  }
38

Praetorius, Simon's avatar
Praetorius, Simon committed
39
  // define boundary values
Praetorius, Simon's avatar
Praetorius, Simon committed
40
  auto parabolic_y = [](auto const& x)
Praetorius, Simon's avatar
Praetorius, Simon committed
41
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
42
    return FieldVector<double,2>{0.0, x[1]*(1.0 - x[1])};
Praetorius, Simon's avatar
Praetorius, Simon committed
43
  };
44

Praetorius, Simon's avatar
Praetorius, Simon committed
45
  FieldVector<double,2> zero(0);
46

Praetorius, Simon's avatar
Praetorius, Simon committed
47
  // set boundary conditions for velocity
Praetorius, Simon's avatar
Praetorius, Simon committed
48
49
50
  prob.boundaryManager()->setBoxBoundary({1,2,2,2});
  prob.addDirichletBC(1, _v, _v, parabolic_y);
  prob.addDirichletBC(2, _v, _v, zero);
51

Praetorius, Simon's avatar
Praetorius, Simon committed
52
  AdaptInfo adaptInfo("adapt");
53

Praetorius, Simon's avatar
Praetorius, Simon committed
54
55
56
  // assemble and solve system
  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
57

Praetorius, Simon's avatar
Praetorius, Simon committed
58
59
  // output solution
  prob.writeFiles(adaptInfo);
60

Praetorius, Simon's avatar
Praetorius, Simon committed
61
  return 0;
62
}