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.

vtkwriter.cc 3.35 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#ifdef HAVE_CONFIG_H
# include "config.h"
#endif

#include <iostream>
#include <vector>

#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions

#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>

22
#include <dune/vtk/vtkwriter.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
23
24
25
26

using namespace Dune;
using namespace Dune::Functions;

27
28
29
30
31
32
33
34
35
36
37
38
using TestCases = std::set<std::tuple<std::string,Vtk::FormatTypes,Vtk::DataTypes>>;
static TestCases test_cases = {
  {"ascii32", Vtk::ASCII, Vtk::FLOAT32},
  {"bin32", Vtk::BINARY, Vtk::FLOAT32},
  {"zlib32", Vtk::COMPRESSED, Vtk::FLOAT32},
  {"ascii64", Vtk::ASCII, Vtk::FLOAT64},
  {"bin64", Vtk::BINARY, Vtk::FLOAT64},
  {"zlib64", Vtk::COMPRESSED, Vtk::FLOAT64},
};

template <class GridView>
void write (std::string prefix, GridView const& gridView)
Praetorius, Simon's avatar
Praetorius, Simon committed
39
40
41
42
{
  using namespace BasisFactory;
  auto basis = makeBasis(gridView, lagrange<1>());

Praetorius, Simon's avatar
Praetorius, Simon committed
43
44
45
46
  FieldVector<double,GridView::dimensionworld> c;
  if (GridView::dimensionworld > 0) c[0] = 11.0;
  if (GridView::dimensionworld > 1) c[1] = 7.0;
  if (GridView::dimensionworld > 2) c[2] = 3.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
47

Praetorius, Simon's avatar
Praetorius, Simon committed
48
  assert(basis.dimension() > 0);
49
50
  std::vector<double> vec(basis.dimension());
  interpolate(basis, vec, [&c](auto const& x) { return c.dot(x); });
Praetorius, Simon's avatar
Praetorius, Simon committed
51
52

  // write discrete global-basis function
53
  auto p1Interpol = makeDiscreteGlobalBasisFunction<double>(basis, vec);
Praetorius, Simon's avatar
Praetorius, Simon committed
54
55

  // write analytic function
56
57
58
  auto p1Analytic = makeAnalyticGridViewFunction([&c](auto const& x) { return c.dot(x); }, gridView);

  for (auto const& test_case : test_cases) {
59
60
61
62
63
    VtkWriter<GridView> vtkWriter(gridView, std::get<1>(test_case), std::get<2>(test_case));
    vtkWriter.addPointData(p1Interpol, "p1");
    vtkWriter.addCellData(p1Interpol, "p0");
    vtkWriter.addPointData(p1Analytic, "q1");
    vtkWriter.addCellData(p1Analytic, "q0");
Praetorius, Simon's avatar
Praetorius, Simon committed
64
    vtkWriter.write(prefix + "_" + std::to_string(GridView::dimensionworld) + "d_" + std::get<0>(test_case) + ".vtu");
65
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
66
}
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

template <int I>
using int_ = std::integral_constant<int,I>;

int main (int argc, char** argv)
{
  Dune::MPIHelper::instance(argc, argv);

#ifdef HAVE_UG
  // Test VtkWriter for UGGrid
  Hybrid::forEach(std::make_tuple(int_<2>{}, int_<3>{}), [](auto dim)
  {
    using GridType = UGGrid<dim.value>;
    {
      FieldVector<double,dim.value> lowerLeft; lowerLeft = 0.0;
      FieldVector<double,dim.value> upperRight; upperRight = 1.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
83
      auto numElements = filledArray<dim.value,unsigned int>(8);
84
85
86
87
88
89
90
91
92
93
94
95
      auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);

      write("ug", gridPtr->leafGridView());
    }
  });
#endif

  // Test VtkWriter for YaspGrid
  Hybrid::forEach(std::make_tuple(int_<1>{}, int_<2>{}, int_<3>{}), [](auto dim)
  {
    using GridType = YaspGrid<dim.value>;
    FieldVector<double,dim.value> upperRight; upperRight = 1.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
96
97
    auto numElements = filledArray<dim.value,int>(8);
    GridType grid(upperRight, numElements, 0, 0);
98
99
100
    write("yasp", grid.leafGridView());
  });
}