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 3.21 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::Indices;
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
  using Param = LagrangeBasis<Grid, 2>;
  ProblemStat<Param> prob("ellipt");
Praetorius, Simon's avatar
Praetorius, Simon committed
34
  prob.initialize(INIT_ALL);
35

36
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
Praetorius, Simon's avatar
Praetorius, Simon committed
37
  prob.addMatrixOperator(opL, 0, 0);
38

39
  auto opForce = makeOperator(tag::test{}, f, 6);
Praetorius, Simon's avatar
Praetorius, Simon committed
40
  prob.addVectorOperator(opForce, 0);
41

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

  AdaptInfo adaptInfo("adapt");
46

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

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

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

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

#if WRITE_FILES
74
    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
Praetorius, Simon's avatar
Praetorius, Simon committed
75
    vtkWriter.addVertexData(prob.solution(0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
76
    vtkWriter.write("u_" + std::to_string(l));
77
78
#endif
  }
79

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

85
86
87
88
  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
89
  }
90

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

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