stokes3.cc 1.73 KB
Newer Older
1
2
3
4
#include <iostream>
#include <ctime>
#include <cmath>

5
6
7
8
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/assembler/StokesOperator.hpp>
9
10
11

using namespace AMDiS;

12
13
14
15
using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

16
17
18
19
int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);

20
  StokesProblem prob("stokes");
21
22
23
24
25
  prob.initialize(INIT_ALL);

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

26
27
28
29
  // tree-paths for components
  auto _v = 0_c;
  auto _p = 1_c;

30
31
  auto opStokes = makeOperator(tag::stokes{}, viscosity);
  prob.addMatrixOperator(opStokes, treepath(), treepath());
32

33
34
35
  auto opZero = makeOperator(tag::test_trial{}, 0.0);
  prob.addMatrixOperator(opZero, _p, _p);

36
37
38
39
40
  // 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
41
  auto parabolic_y = [](auto const& x) -> FieldVector<double,2>
42
43
44
45
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };

46
  auto zero = [](auto const& x) -> FieldVector<double,2>
47
48
49
50
51
  {
    return {0.0, 0.0};
  };

  // set boundary conditions for velocity
52
53
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
54
55

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

58
  AdaptInfo adaptInfo("adapt");
59
60
61

  // assemble and solve system
  prob.buildAfterCoarsen(adaptInfo, Flag(0));
62

63
64
65
  prob.solve(adaptInfo);

  // output solution
66
  prob.writeFiles(adaptInfo);
67
68
69
70

  AMDiS::finalize();
  return 0;
}