boundary.cc 3.02 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
#include <config.h>

3
4
5
6
7
8
9
10
11
#include <iostream>

#if HAVE_DUNE_SPGRID
#include <dune/grid/spgrid.hh>
#endif

#include <amdis/AMDiS.hpp>
#include <amdis/Integrate.hpp>
#include <amdis/ProblemStat.hpp>
12
#include <amdis/LocalOperators.hpp>
13
14
15
16

using namespace AMDiS;

// 1 component with polynomial degree 2
17
using Param   = YaspGridBasis<GRIDDIM, 2>;
18
19
20
21
22
23
24
using ElliptProblem = ProblemStat<Param>;

template <class SetBoundary>
void run(SetBoundary setBoundary)
{
  ElliptProblem prob("ellipt");
  prob.initialize(INIT_ALL);
25
  setBoundary(*prob.boundaryManager());
26
27

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
28
  prob.addMatrixOperator(opL);
29
30
31
32
33
34
35

  auto f = [](auto const& x) {
    double r2 = dot(x,x);
    double ux = std::exp(-10.0 * r2);
    return -(400.0 * r2 - 20.0 * x.size()) * ux;
  };
  auto opForce = makeOperator(tag::test{}, f, 6);
36
  prob.addVectorOperator(opForce);
37
38
39

  // set boundary condition
  auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
40
  prob.addDirichletBC(BoundaryType{1}, g);
41
42
43
44
45
46

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);

47
  double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
  msg("error_L2 = {}", errorL2);
}

void run_periodic()
{
#if HAVE_DUNE_SPGRID
  Dune::SPCube<double,2> cube({0.0,0.0},{1.0,1.0});
  Dune::SPDomain<double,2> domain({cube}, Dune::SPTopology<2>(0b01));
  Dune::SPGrid<double,2> grid(domain, Dune::SPMultiIndex<2>({2,2}));
  using Grid = Dune::SPGrid<double,2>;
#else
  Dune::YaspGrid<2> grid({1.0,1.0},{2,2},std::bitset<2>("10"),0);
  using Grid = Dune::YaspGrid<2>;
#endif

63
  using Traits = LagrangeBasis<Grid, 2>;
64
65
  ProblemStat<Traits> prob("ellipt", grid);
  prob.initialize(INIT_ALL);
66
  prob.boundaryManager()->setBoxBoundary({-1,-1,1,1});
67
68

  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
69
  prob.addMatrixOperator(opL);
70
71

  auto opForce = makeOperator(tag::test{}, 1.0, 6);
72
  prob.addVectorOperator(opForce);
73
74
75
76
77
78

  // set boundary condition
  auto g = [](auto const& x) {
    return std::sin(0.5*M_PI*x[1])*std::sin(2*M_PI * (0.25 + x[0]))
         + std::cos(0.5*M_PI*x[1])*std::sin(2*M_PI * (-0.25 + x[0]));
  };
79
  prob.addDirichletBC(BoundaryType{1}, g);
80
81
82
83
84
85
86
87
88

  typename ElliptProblem::WorldMatrix A{{1.0,0.0}, {0.0,1.0}};
  typename ElliptProblem::WorldVector b{1.0, 0.0};
  prob.addPeriodicBC(BoundaryType{-1}, A, b);

  AdaptInfo adaptInfo("adapt");

  prob.assemble(adaptInfo);
  prob.solve(adaptInfo);
Praetorius, Simon's avatar
Praetorius, Simon committed
89
  prob.writeFiles(adaptInfo);
90
91
92
93
94
}


int main(int argc, char** argv)
{
95
  Environment env(argc, argv);
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

  auto b = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };

  auto setBoundary1 = [](auto& boundaryManager)
  {
    boundaryManager.setBoxBoundary({1,1,1,1});
  };
  run(setBoundary1);

  auto setBoundary2 = [b](auto& boundaryManager)
  {
    boundaryManager.setIndicator([b](auto const& x) { return b(x) ? 1 : 0; });
  };
  run(setBoundary2);

  auto setBoundary3 = [b](auto& boundaryManager)
  {
    boundaryManager.setPredicate(b, 1);
  };
  run(setBoundary3);

  run_periodic();

  return 0;
}