Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

stokes0.cc 2.43 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
#include <amdis/AMDiS.hpp>
10
#include <amdis/LocalOperators.hpp>
11
#include <amdis/ProblemStat.hpp>
12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

#ifdef DOW
#undef DOW
#endif
#define DOW AMDIS_DOW

using namespace AMDiS;

// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1

using Grid = Dune::YaspGrid<AMDIS_DIM>;
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

int main(int argc, char** argv)
{
Praetorius, Simon's avatar
Praetorius, Simon committed
28
  AMDiS::init(argc, argv);
29

Praetorius, Simon's avatar
Praetorius, Simon committed
30 31
  StokesProblem prob("stokes");
  prob.initialize(INIT_ALL);
32

Praetorius, Simon's avatar
Praetorius, Simon committed
33 34
  double viscosity = 1.0;
  Parameters::get("stokes->viscosity", viscosity);
35

Praetorius, Simon's avatar
Praetorius, Simon committed
36 37 38
  // tree-paths for components
  auto _v = Dune::Indices::_0;
  auto _p = Dune::Indices::_1;
39

Praetorius, Simon's avatar
Praetorius, Simon committed
40 41 42 43
  // <viscosity*grad(u_i), grad(v_i)>
  for (std::size_t i = 0; i < DOW; ++i) {
    auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
    prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
44

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

Praetorius, Simon's avatar
Praetorius, Simon committed
49 50 51 52
    // <q, d_i(u_i)>
    auto opDiv = makeOperator(tag::test_partialtrial{i}, 1.0);
    prob.addMatrixOperator(opDiv, _p, treepath(_v,i));
  }
53

Praetorius, Simon's avatar
Praetorius, Simon committed
54 55
  auto opZero = makeOperator(tag::test_trial{}, 0.0);
  prob.addMatrixOperator(opZero, _p, _p);
56

Praetorius, Simon's avatar
Praetorius, Simon committed
57 58 59
  // 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; };
60

Praetorius, Simon's avatar
Praetorius, Simon committed
61 62 63 64 65
  // define boundary values
  auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,DOW>
  {
    return {0.0, x[1]*(1.0 - x[1])};
  };
66

Praetorius, Simon's avatar
Praetorius, Simon committed
67 68 69 70
  auto zero = [](auto const& x) -> Dune::FieldVector<double,DOW>
  {
    return {0.0, 0.0};
  };
71

Praetorius, Simon's avatar
Praetorius, Simon committed
72 73 74
  // set boundary conditions for velocity
  prob.addDirichletBC(left, _v, _v, parabolic_y);
  prob.addDirichletBC(not_left, _v, _v, zero);
75

Praetorius, Simon's avatar
Praetorius, Simon committed
76
  prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);
77

Praetorius, Simon's avatar
Praetorius, Simon committed
78
  AdaptInfo adaptInfo("adapt");
79

Praetorius, Simon's avatar
Praetorius, Simon committed
80 81
  // assemble and solve system
  prob.assemble(adaptInfo);
82

83
#ifdef DEBUG_MTL
Praetorius, Simon's avatar
Praetorius, Simon committed
84 85 86 87
  // write matrix to file
  mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
  out << prob.systemMatrix().matrix();
  std::cout << prob.systemMatrix().matrix() << '\n';
88
#endif
89

Praetorius, Simon's avatar
Praetorius, Simon committed
90
  prob.solve(adaptInfo);
91

Praetorius, Simon's avatar
Praetorius, Simon committed
92 93
  // output solution
  prob.writeFiles(adaptInfo);
94

Praetorius, Simon's avatar
Praetorius, Simon committed
95 96
  AMDiS::finalize();
  return 0;
97
}