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.07 KB
Newer Older
1
2
3
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
4
5
#include <iostream>

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

12
using namespace AMDiS;
13
using namespace Dune::Indices;
14

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

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

23
24
25
26
27
28
29
30
31
32
33
34
35
  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};
  };
36

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

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
46
  // set boundary condition
47
  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; };
Praetorius, Simon's avatar
Praetorius, Simon committed
48
  prob.addDirichletBC(boundary, 0, 0, g);
49
50

  AdaptInfo adaptInfo("adapt");
51

52
53
54
55
  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) {
56
    prob.grid().globalRefine(1);
57
    auto gridView = prob.gridView();
58
59
60
61
62
63

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

64
    prob.globalBasis().update(gridView);
65
    prob.solutionVector().compress();
Praetorius, Simon's avatar
Praetorius, Simon committed
66
    prob.assemble(adaptInfo);
67
68
    prob.solve(adaptInfo);

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

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

81
  std::cout << "\n";
82
  msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}",
83
      "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
84
  msg_("{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
  for (int i = 0; i < numLevels; ++i)
93
    msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}",
94
        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
}