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.15 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>

6
7
#include <fmt/core.h>

8
9
10
11
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/common/Literals.hpp>
12
#include <amdis/gridfunctions/Integrate.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
13

14
using namespace AMDiS;
15
using namespace Dune::Indices;
16

Praetorius, Simon's avatar
Praetorius, Simon committed
17
// 1 component with polynomial degree 1
18
19
using Param   = YaspGridBasis<AMDIS_DIM, 1>;
using ElliptProblem = ProblemStat<Param>;
Praetorius, Simon's avatar
Praetorius, Simon committed
20
21
22
23

int main(int argc, char** argv)
{
  AMDiS::init(argc, argv);
24

25
26
27
28
29
30
31
32
33
34
35
36
37
  int numLevels = AMDIS_DIM == 2 ? 8 : 5;
  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)); };
  auto grad_g = [](auto const& x) -> FieldMatrix<double,1,AMDIS_DIM> {
    return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
  };
38

39
  ElliptProblem prob("ellipt");
Praetorius, Simon's avatar
Praetorius, Simon committed
40
  prob.initialize(INIT_ALL);
41

42
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
43
  prob.addMatrixOperator(opL, _0, _0);
44

45
  auto opForce = makeOperator(tag::test{}, f, 6);
46
  prob.addVectorOperator(opForce, _0);
47

Praetorius, Simon's avatar
Praetorius, Simon committed
48
  // set boundary condition
49
50
51
52
  auto boundary = [](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; };
  prob.addDirichletBC(boundary, _0, _0, g);

  AdaptInfo adaptInfo("adapt");
53

54
55
56
57
  std::vector<double> errL2; errL2.reserve(numLevels);
  std::vector<double> errH1; errH1.reserve(numLevels);
  std::vector<double> widths; widths.reserve(numLevels);
  for (int i = 0; i < numLevels; ++i) {
58
    prob.grid().globalRefine(1);
59
    auto gridView = prob.gridView();
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74

    double h = 0;
    for (auto const& e : edges(gridView))
      h = std::max(h, e.geometry().volume());
    widths.push_back(h);

    prob.buildAfterCoarsen(adaptInfo, Flag(0));
    prob.solve(adaptInfo);

    double errorL2 = integrate(sqr(g - prob.getSolution(_0)), gridView, 6);
    errL2.push_back(std::sqrt(errorL2));
    double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.getSolution(_0))), gridView, 6);
    errH1.push_back(std::sqrt(errorH1));

#if WRITE_FILES
75
    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
76
77
78
79
    vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
    vtkWriter.write("u_" + std::to_string(i));
#endif
  }
80

81
82
83
84
  std::cout << "\n";
  fmt::print("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n",
             "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
  fmt::print("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
85

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

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

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