stokes3.cc 1.76 KB
Newer Older
1
2
3
4
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

5
6
7
8
#include <iostream>
#include <ctime>
#include <cmath>

9
10
11
12
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/assembler/StokesOperator.hpp>
13
14
15

using namespace AMDiS;

16
17
18
19
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

20
21
22
23
int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

24
  StokesProblem prob("stokes");
25
26
27
28
29
  prob.initialize(INIT_ALL);

  double viscosity = 1.0;
  Parameters::get("stokes->viscosity", viscosity);

30
  // tree-paths for components
31
32
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
33

34
  auto opStokes = makeOperator(tag::stokes{}, viscosity);
35
  prob.addMatrixOperator(opStokes);
36

37
38
39
  auto opZero = makeOperator(tag::test_trial{}, 0.0);
  prob.addMatrixOperator(opZero, _p, _p);

40
41
42
43
44
  // define boundary regions
  auto left = [](auto const& x) { return x[0] < 1.e-8; };
  auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };

  // define boundary values
45
  auto parabolic_y = [](auto const& x) -> FieldVector<double,2>
46
47
48
49
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };

50
  auto zero = [](auto const& x) -> FieldVector<double,2>
51
52
53
54
55
  {
    return {0.0, 0.0};
  };

  // set boundary conditions for velocity
56
57
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
58
59

  // set point constraint for pressure
60
  prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);
61

62
  AdaptInfo adaptInfo("adapt");
63
64

  // assemble and solve system
Praetorius, Simon's avatar
Praetorius, Simon committed
65
  prob.assemble(adaptInfo);
66
67
68
  prob.solve(adaptInfo);

  // output solution
69
  prob.writeFiles(adaptInfo);
70
71
72
73

  AMDiS::finalize();
  return 0;
}