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.04 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
7
8
9
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/common/Literals.hpp>
10
#include <amdis/gridfunctions/Integrate.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);
41
  prob.addMatrixOperator(opL, _0, _0);
42

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

Praetorius, Simon's avatar
Praetorius, Simon committed
46
  // set boundary condition
47
48
49
50
  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");
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.buildAfterAdapt(adaptInfo, Flag(0));
65
66
67
68
69
70
71
72
    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
73
    Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView);
74
75
76
77
    vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
    vtkWriter.write("u_" + std::to_string(i));
#endif
  }
78

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

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

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

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