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.37 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

Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
41
42
43
  using Grid = Dune::YaspGrid<AMDIS_DIM>;
  auto gridPtr = MeshCreator<Grid>::create("mesh");
  auto basis = makeBasis(gridPtr->leafGridView(), lagrange<1>());

  auto prob = makeProblemStat("ellipt", *gridPtr, basis);
Praetorius, Simon's avatar
Praetorius, Simon committed
44
  prob.initialize(INIT_ALL);
45

46
  auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
Praetorius, Simon's avatar
Praetorius, Simon committed
47
  prob.addMatrixOperator(opL);
48

49
  auto opForce = makeOperator(tag::test{}, f, 6);
50
  prob.addVectorOperator(opForce, _0);
51

Praetorius, Simon's avatar
Praetorius, Simon committed
52
  // set boundary condition
53
54
55
56
  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");
57

58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
  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) {
    prob.getGrid()->globalRefine(1);
    auto gridView = prob.getGlobalBasis()->gridView();

    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
    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
    vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
    vtkWriter.write("u_" + std::to_string(i));
#endif
  }
84

85
86
87
88
  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");
89

90
91
92
93
  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
94
  }
95

96
97
98
  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]);
99

Praetorius, Simon's avatar
Praetorius, Simon committed
100
101
  AMDiS::finalize();
  return 0;
102
}