ellipt.cc 3.05 KB
Newer Older
1
#include <amdis/AMDiS.hpp>
2
#include <amdis/Integrate.hpp>
3
#include <amdis/LocalOperators.hpp>
4
5
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
6

7
8
9
#include <dune/grid/yaspgrid.hh>
using Grid = Dune::YaspGrid<GRIDDIM>;

10
using namespace AMDiS;
11
using namespace Dune::Functions::BasisFactory;
12

Praetorius, Simon's avatar
Praetorius, Simon committed
13
14
int main(int argc, char** argv)
{
15
  Environment env(argc, argv);
16
  Dune::Timer t;
17

18
  int numLevels = GRIDDIM == 2 ? 6 : 4;
19
20
21
22
23
24
25
26
27
  if (argc > 2)
    numLevels = std::atoi(argv[2]);

  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 g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
28
  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,GRIDDIM> {
29
30
    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
  };
31

32
33
34
  auto grid = MeshCreator<Grid>{"elliptMesh"}.create();

  ProblemStat prob("ellipt", *grid, lagrange<2>());
Praetorius, Simon's avatar
Praetorius, Simon committed
35
  prob.initialize(INIT_ALL);
36

37
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
38
  prob.addMatrixOperator(opL);
39

40
  auto opForce = makeOperator(tag::test{}, f, 6);
41
  prob.addVectorOperator(opForce);
42

Praetorius, Simon's avatar
Praetorius, Simon committed
43
  // set boundary condition
44
  prob.addDirichletBC(1, g);
45
46

  AdaptInfo adaptInfo("adapt");
47

48
49
50
  std::vector<double> errL2; errL2.reserve(numLevels);
  std::vector<double> errH1; errH1.reserve(numLevels);
  std::vector<double> widths; widths.reserve(numLevels);
51
52
  for (int l = 0; l < numLevels; ++l) {
    prob.globalRefine(1);
53
54

    double h = 0;
55
56
57
58
59
60
61
62
63
    for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) {
      auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type());
      auto geo = e.geometry();
      for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges
        auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM));
        auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM));
        h = std::max(h, (v0 - v1).two_norm());
      }
    }
64
65
    widths.push_back(h);

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

69
    double errorL2 = integrate(sqr(g - prob.solution()), prob.gridView(), 6);
70
    errL2.push_back(std::sqrt(errorL2));
71
    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution())), prob.gridView(), 6);
72
73
    errH1.push_back(std::sqrt(errorH1));
  }
74

Praetorius, Simon's avatar
Praetorius, Simon committed
75
76
77
  // write last solution to file
  prob.writeFiles(adaptInfo);

78
  msg("");
79
  msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
80
      "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
81
  msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
82

83
84
85
86
  std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
  for (int i = 1; i < numLevels; ++i) {
    eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
    eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
Praetorius, Simon's avatar
Praetorius, Simon committed
87
  }
88

89
  for (int i = 0; i < numLevels; ++i)
90
    msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
91
        i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
92

93
  msg("elapsed time: {} seconds", t.elapsed());
Praetorius, Simon's avatar
Praetorius, Simon committed
94
  return 0;
95
}