datareader.cc 5.58 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
// -*- 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 <cmath>
#include <iostream>
#include <vector>

#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/common/filledarray.hh>
#include <dune/functions/gridfunctions/analyticgridviewfunction.hh>
#include <dune/grid/utility/structuredgridfactory.hh>
#include <dune/vtk/vtkreader.hh>
#include <dune/vtk/datacollectors/lagrangedatacollector.hh>
#include <dune/vtk/gridcreators/continuousgridcreator.hh>
#include <dune/vtk/gridcreators/lagrangegridcreator.hh>
#include <dune/vtk/writers/vtkunstructuredgridwriter.hh>

#if HAVE_DUNE_UGGRID
#include <dune/grid/uggrid.hh>
#endif

#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif

using namespace Dune;

template <class GridType, class GridCreator>
void run(std::string const& prefix)
{
  std::tuple f_tuple{0,
    [](auto const& x) { return std::sin(x[0]*x[0]*x[0] + 20); },
    [](auto const& x) { return std::sin(x[0]*x[0]) + std::cos(x[1] + x[0]); },
    [](auto const& x) { return std::sin(x[0] + x[1] + x[2]) + std::cos(x[0]*x[1]*x[2]) + std::sin(x[0]*x[2])*std::cos(x[1] + x[2]); }
  };

  constexpr int dim = GridType::dimension;
  auto f = std::get<dim>(f_tuple);

  using GridView = typename GridType::LeafGridView;
  {
    FieldVector<double,dim> lowerLeft; lowerLeft = 0.0;
    FieldVector<double,dim> upperRight; upperRight = 1.0;
    auto numElements = filledArray<dim,unsigned int>(4);
    auto gridPtr = StructuredGridFactory<GridType>::createSimplexGrid(lowerLeft, upperRight, numElements);
    auto& grid = *gridPtr;

    GridView gridView = grid.leafGridView();
    auto data = Functions::makeAnalyticGridViewFunction(f, gridView);

56
    VtkUnstructuredGridWriter<GridView> vtkWriter1(gridView, Vtk::FormatTypes::ASCII);
57
58
59
60
    vtkWriter1.addPointData(data, "f");
    vtkWriter1.addCellData(data, "g");
    vtkWriter1.write(prefix + "_ascii.vtu");

61
    VtkUnstructuredGridWriter<GridView> vtkWriter2(gridView, Vtk::FormatTypes::BINARY);
62
63
64
65
    vtkWriter2.addPointData(data, "f");
    vtkWriter2.addCellData(data, "g");
    vtkWriter2.write(prefix + "_binary.vtu");

Praetorius, Simon's avatar
Praetorius, Simon committed
66
    VtkUnstructuredGridWriter<GridView> vtkWriter3(gridView, Vtk::FormatTypes::COMPRESSED, Vtk::DataTypes::FLOAT64);
67
68
69
70
71
    vtkWriter3.addPointData(data, "f");
    vtkWriter3.addCellData(data, "g");
    vtkWriter3.write(prefix + "_compressed.vtu");

    if constexpr(std::is_same_v<GridCreator, Vtk::LagrangeGridCreator<GridType>>) {
72
      VtkUnstructuredGridWriter<GridView, Vtk::LagrangeDataCollector<GridView,3>> vtkWriter4(gridView, Vtk::FormatTypes::BINARY, Vtk::DataTypes::FLOAT64);
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
      vtkWriter4.addPointData(data, "f");
      vtkWriter4.addCellData(data, "g");
      vtkWriter4.write(prefix + "_order3.vtu");
    }
  }

  {
    std::cout << "read '" << prefix << "_ascii.vtu'..." << std::endl;
    VtkReader<GridType, GridCreator, float> reader;
    reader.read(prefix + "_ascii.vtu");
    auto gridPtr = reader.createGrid();
    auto& grid = *gridPtr;

    auto gf = reader.getPointData("f");
    auto gg = reader.getCellData("g");

89
    VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::FormatTypes::ASCII);
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
    vtkWriter.addPointData(gf, "f");
    vtkWriter.addCellData(gg, "g");
    vtkWriter.write(prefix + "_ascii_out.vtu");
  }

  {
    std::cout << "read '" << prefix << "_binary.vtu'..." << std::endl;
    VtkReader<GridType, GridCreator, float> reader;
    reader.read(prefix + "_binary.vtu");
    auto gridPtr = reader.createGrid();
    auto& grid = *gridPtr;

    auto gf = reader.getPointData("f");
    auto gg = reader.getCellData("g");

105
    VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::FormatTypes::ASCII);
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
    vtkWriter.addPointData(gf, "f");
    vtkWriter.addCellData(gg, "g");
    vtkWriter.write(prefix + "_binary_out.vtu");
  }

  {
    std::cout << "read '" << prefix << "_compressed.vtu'..." << std::endl;
    VtkReader<GridType, GridCreator> reader;
    reader.read(prefix + "_compressed.vtu");
    auto gridPtr = reader.createGrid();
    auto& grid = *gridPtr;

    auto gf = reader.getPointData("f");
    auto gg = reader.getCellData("g");

121
    VtkUnstructuredGridWriter<GridView> vtkWriter(grid.leafGridView(), Vtk::FormatTypes::ASCII);
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
    vtkWriter.addPointData(gf, "f");
    vtkWriter.addCellData(gg, "g");
    vtkWriter.write(prefix + "_compressed_out.vtu");
  }

  if constexpr(std::is_same_v<GridCreator, Vtk::LagrangeGridCreator<GridType>>) {
    std::cout << "read '" << prefix << "_order3.vtu'..." << std::endl;
    VtkReader<GridType, GridCreator> reader;
    reader.read(prefix + "_order3.vtu");
    auto gridPtr = reader.createGrid();
    auto& grid = *gridPtr;

    auto gf = reader.getPointData("f");
    auto gg = reader.getCellData("g");

137
    VtkUnstructuredGridWriter<GridView, Vtk::LagrangeDataCollector<GridView,3>> vtkWriter(grid.leafGridView(), Vtk::FormatTypes::ASCII);
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
    vtkWriter.addPointData(gf, "f");
    vtkWriter.addCellData(gg, "g");
    vtkWriter.write(prefix + "_order3_out.vtu");
  }
}

template <class GridType>
void run(std::string const& prefix)
{
  run<GridType, Vtk::ContinuousGridCreator<GridType>>(prefix + "_continuous");
  run<GridType, Vtk::LagrangeGridCreator<GridType>>(prefix + "_lagrange");
}

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

  run<UGGrid<2>>("uggrid_2d");
  run<UGGrid<3>>("uggrid_3d");

#if HAVE_DUNE_ALUGRID
  run<Dune::ALUGrid<2,2,Dune::simplex,Dune::conforming>>("alugrid_2d");
  run<Dune::ALUGrid<3,3,Dune::simplex,Dune::conforming>>("alugrid_3d");
#endif
}