convergence.cc 8.7 KB
Newer Older
1
#ifdef HAVE_CONFIG_H
2
#include "config.h"
3
#endif
4

5
#include <iostream>
6

7
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
8
#include <dune/common/test/testsuite.hh>
9
#include <dune/curvedsurfacegrid/curvedsurfacegrid.hh>
10
11
12
#include <dune/curvedsurfacegrid/geometries/ellipsoid.hh>
#include <dune/curvedsurfacegrid/geometries/sphere.hh>
#include <dune/curvedsurfacegrid/geometries/torus.hh>
13
#include <dune/foamgrid/foamgrid.hh>
14
15
16
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/localfunctions/lagrange/pk.hh>
17
18
19
#include <dune/vtk/vtkwriter.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>

20
21
22
23
24
25
// 1 .. Sphere
// 2 .. Ellipsoid
// 3 .. Torus
#ifndef GEOMETRY_TYPE
#define GEOMETRY_TYPE 1
#endif
26

27
28
const int order = 3;
const int quad_order = order+6;
29
30
31
32

#if GEOMETRY_TYPE < 3
const int num_levels = 4;
#else
33
const int num_levels = 3;
34
#endif
35
36

// Test Hausdorf-distance
37
template <class Grid>
38
typename Grid::ctype inf_error (const Grid& grid)
39
{
40
  std::cout << "  inf_error..." << std::endl;
41
42
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;
43
  auto projection = grid.gridFunction().f();
44
  for (const auto& element : elements(grid.leafGridView()))
45
46
47
  {
    auto geometry = element.geometry();

48
49
    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
50
51
52
53
54
55
56
57
58
59
      auto x = geometry.global(qp.position());
      auto y = projection(x);
      dist = std::max(dist, (x - y).two_norm());
    }
  }

  return dist;
}

// Test integrated Hausdorf-distance
60
template <class Grid>
61
typename Grid::ctype L2_error (const Grid& grid)
62
{
63
  std::cout << "  L2_error..." << std::endl;
64
65
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;
66
  auto projection = grid.gridFunction().f();
67
  for (const auto& element : elements(grid.leafGridView()))
68
69
70
  {
    auto geometry = element.geometry();

71
72
    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
73
74
75
76
77
78
79
80
81
82
      auto x = geometry.global(qp.position());
      auto y = projection(x);
      dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
    }
  }

  return std::sqrt(dist);
}

// Test integrated error in normal vectors
83
template <class Grid>
84
typename Grid::ctype normal_error (const Grid& grid)
85
{
86
  std::cout << "  normal_error..." << std::endl;
87
88
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
  typename Grid::ctype dist = 0;
89
90

  auto projection = grid.gridFunction().f();
91
  for (const auto& element : elements(grid.leafGridView()))
92
93
  {
    auto geometry = element.geometry();
94
95
    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
96
      auto x = geometry.impl().normal(qp.position());
97
      auto y = projection.normal(geometry.global(qp.position()));
98
99
100
101
102
103
104
105
      dist += (x - y).two_norm2() * qp.weight() * geometry.integrationElement(qp.position());
    }
  }

  return std::sqrt(dist);
}

// Test integrated error in mean curvature
106
template <class Grid>
107
typename Grid::ctype curvature_error (const Grid& grid)
108
{
109
  std::cout << "  curvature_error..." << std::endl;
110
111
  using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
  using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
112
  using FiniteElementType = Dune::PkLocalFiniteElement<typename Grid::ctype, typename Grid::ctype, 2, order+4>;
113
114
115
116
117
118
119
  FiniteElementType fe;

  using LBTraits = typename FiniteElementType::Traits::LocalBasisType::Traits;
  std::vector<typename LBTraits::JacobianType> shapeGradients;
  std::vector<GlobalCoordinate> gradients;
  std::vector<GlobalCoordinate> normals;

120
121
  auto projection = grid.gridFunction().f();

122
  typename Grid::ctype dist = 0;
123
  for (const auto& element : elements(grid.leafGridView()))
124
125
  {
    auto geometry = element.geometry();
126
127
128
129
    const auto& localBasis = fe.localBasis();

    { // get coefficients of normal vectors
      const auto& localInterpolation = fe.localInterpolation();
130

131
132
133
134
135
136
137
      auto n = [&](const auto& local) -> GlobalCoordinate
      {
        auto x = geometry.impl().normal(local);
        return x;
      };
      localInterpolation.interpolate(n, normals);
    }
138
139
140
141

    shapeGradients.resize(localBasis.size());
    gradients.resize(localBasis.size());

142
143
    const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
    for (const auto& qp : quadRule) {
144
145
146
147
148
149
150
151
152
153
      auto jInvT = geometry.jacobianInverseTransposed(qp.position());
      localBasis.evaluateJacobian(qp.position(), shapeGradients);

      for (std::size_t i = 0; i < gradients.size(); ++i)
        jInvT.mv(shapeGradients[i][0], gradients[i]);

      typename Grid::ctype H = 0;
      for (std::size_t i = 0; i < gradients.size(); ++i)
        H += normals[i].dot(gradients[i]);

154
155
      auto H_exact = projection.mean_curvature(geometry.global(qp.position()));
      dist += std::abs(H/2 - H_exact) * qp.weight() * geometry.integrationElement(qp.position());
156
157
158
159
160
161
162
163
    }
  }

  return dist;
}

// longest edge in the grid
template <class Grid>
164
typename Grid::ctype edge_length (const Grid& grid)
165
166
{
  typename Grid::ctype h = 0;
167
  for (const auto& e : edges(grid.hostGrid().leafGridView()))
168
169
170
171
172
    h = std::max(h, e.geometry().volume());

  return h;
}

173

174
int main (int argc, char** argv)
175
176
{
  using namespace Dune;
177
  MPIHelper::instance(argc, argv);
178

Praetorius, Simon's avatar
Praetorius, Simon committed
179
180
  auto grid_order = std::integral_constant<int,order>{};

181
  using HostGrid = FoamGrid<2,3>;
182
183
#if GEOMETRY_TYPE == 1
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "sphere.msh");
184
  auto gridFct = sphereGridFunction<HostGrid>(1.0);
185
186
#elif GEOMETRY_TYPE == 2
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "ellipsoid3.msh");
187
  auto gridFct = ellipsoidGridFunction<HostGrid>(1.0, 0.75, 1.25);
188
189
#elif GEOMETRY_TYPE == 3
  std::unique_ptr<HostGrid> hostGrid = GmshReader<HostGrid>::read( DUNE_GRID_PATH "torus.msh");
190
  auto gridFct = torusGridFunction<HostGrid>(2.0, 1.0);
191
#endif
192

193
194
  CurvedSurfaceGrid grid(*hostGrid, gridFct, grid_order);
  using Grid = decltype(grid);
195
196
197
198
199
  std::vector<typename Grid::ctype> inf_errors, L2_errors, normal_errors, curvature_errors, edge_lengths;
  for (int i = 0; i < num_levels; ++i) {
    if (i > 0)
      grid.globalRefine(1);

200
201
202
203
    using DC = LagrangeDataCollector<typename Grid::LeafGridView, order>;
    VtkUnstructuredGridWriter<typename Grid::LeafGridView, DC> writer(grid.leafGridView());
    writer.write("level_" + std::to_string(i) + ".vtu");

204
    std::cout << "level = " << i << std::endl;
205
206
207
208
    inf_errors.push_back( inf_error(grid) );
    L2_errors.push_back( L2_error(grid) );
    normal_errors.push_back( normal_error(grid) );
    curvature_errors.push_back( curvature_error(grid) );
209
210
211
    edge_lengths.push_back( edge_length(grid) );
  }

212
213
  using std::log;
  using std::abs;
214

215
  // calculate the experimental order of convergence
216
217
  std::vector<typename Grid::ctype> eocInf(num_levels, 0), eocL2(num_levels, 0), eocNormal(num_levels, 0), eocCurvature(num_levels, 0);
  for (int i = 1; i < num_levels; ++i) {
218
219
220
221
    eocInf[i]       = log(inf_errors[i]/inf_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocL2[i]        = log(L2_errors[i]/L2_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocNormal[i]    = log(normal_errors[i]/normal_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
    eocCurvature[i] = log(curvature_errors[i]/curvature_errors[i-1]) / log(edge_lengths[i]/edge_lengths[i-1]);
222
223
  }

224
225
226
227
228
229
230
231
  TestSuite test;
  test.check(abs(eocInf.back()    - (order+1)) < 0.5, "|d|_L^inf != order+1");
  test.check(abs(eocL2.back()     - (order+1)) < 0.5, "|d|_L^2 != order+1");
  test.check(abs(eocNormal.back() - (order+0)) < 0.5, "|n - nh|_L^2 != order");
  test.check(abs(eocCurvature.back()   - (order-1)) < 0.5, "|H - Hh|_L^2 != order-1");


  auto print_line = [](auto i, const auto& data, auto accept) {
232
    std::cout << i;
233
234
235
236
237
238
239
    for (const auto& d : data) {
      std::cout << '\t' << "| ";
      if (accept(d))
        std::cout << d;
      else
        std::cout << '\t';
    }
240
241
242
243
244
245
246
247
248
249
250
    std::cout << std::endl;
  };

  auto print_break = [] {
    std::cout.width(8 + 8*16);
    std::cout.fill('-');
    std::cout << '-' << std::endl;
    std::cout.fill(' ');
  };

  std::cout.setf(std::ios::scientific);
251
  print_line("level", std::vector<std::string>{"err_inf", "eoc_inf", "err_L2", "eoc_L2", "err_norm", "eoc_norm", "err_curv", "eoc_curv"}, [](auto&&) { return true; });
252
253
254
  print_break();

  for (int i = 0; i < num_levels; ++i)
255
256
257
258
259
260
    print_line(i, std::vector<typename Grid::ctype>{inf_errors[i], eocInf[i], L2_errors[i], eocL2[i], normal_errors[i], eocNormal[i], curvature_errors[i], eocCurvature[i]},
      [](auto const& v) { return v != 0.0; });
  print_break();


  return test.exit();
261
}