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.55 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
10
11
12
// -*- 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
13
#include <dune/common/version.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
14
15
16
17
18
19

#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>
Praetorius, Simon's avatar
Praetorius, Simon committed
20
#if HAVE_DUNE_UGGRID
Praetorius, Simon's avatar
Praetorius, Simon committed
21
#include <dune/grid/uggrid.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
22
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
23
24
#include <dune/grid/yaspgrid.hh>

25
#include <dune/vtk/vtkwriter.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
26
27
28
29

using namespace Dune;
using namespace Dune::Functions;

30
31
32
33
34
35
36
37
38
39
40
41
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
42
{
43
#if DUNE_VERSION_LT(DUNE_FUNCTIONS, 2, 6)
44
45
  using namespace BasisBuilder;
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
46
  using namespace BasisFactory;
47
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
48
49
  auto basis = makeBasis(gridView, lagrange<1>());

Praetorius, Simon's avatar
Praetorius, Simon committed
50
51
52
53
  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
54

Praetorius, Simon's avatar
Praetorius, Simon committed
55
  assert(basis.dimension() > 0);
56
57
  std::vector<double> vec(basis.dimension());
  interpolate(basis, vec, [&c](auto const& x) { return c.dot(x); });
Praetorius, Simon's avatar
Praetorius, Simon committed
58
59

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

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

  for (auto const& test_case : test_cases) {
66
67
68
69
70
    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
71
    vtkWriter.write(prefix + "_" + std::to_string(GridView::dimensionworld) + "d_" + std::get<0>(test_case) + ".vtu");
72
  }
Praetorius, Simon's avatar
Praetorius, Simon committed
73
}
74
75
76
77
78
79
80
81

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
82
#if HAVE_DUNE_UGGRID
83
84
85
86
87
88
89
  // 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
90
      auto numElements = filledArray<dim.value,unsigned int>(8);
91
      auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
92
      gridPtr->loadBalance();
Praetorius, Simon's avatar
Praetorius, Simon committed
93
      write("vtkwriter_ug", gridPtr->leafGridView());
94
95
96
97
98
99
100
101
102
    }
  });
#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
103
104
    auto numElements = filledArray<dim.value,int>(8);
    GridType grid(upperRight, numElements, 0, 0);
Praetorius, Simon's avatar
Praetorius, Simon committed
105
    write("vtkwriter_yasp", grid.leafGridView());
106
  });
107
}