Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

ellipt.cc 2.98 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
11
using namespace AMDiS;

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

17
  int numLevels = GRIDDIM == 2 ? 6 : 4;
18
19
20
21
22
23
24
25
26
  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)); };
27
  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,GRIDDIM> {
28
29
    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
  };
30

31
32
  using Param = LagrangeBasis<Grid, 2>;
  ProblemStat<Param> prob("ellipt");
Praetorius, Simon's avatar
Praetorius, Simon committed
33
  prob.initialize(INIT_ALL);
34

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

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

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

  AdaptInfo adaptInfo("adapt");
45

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

    double h = 0;
53
54
55
56
57
58
59
60
61
    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());
      }
    }
62
63
    widths.push_back(h);

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

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

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

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

81
82
83
84
  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
85
  }
86

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

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