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.

FileWriter.hpp 4.22 KB
Newer Older
1
2
3
4
5
6
7
#pragma once

#include <string>
#include <vector>
#include <array>
#include <memory>

8
#include <dune/grid/io/file/vtk/vtkwriter.hh>
9
#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
10
#include <dune/geometry/referenceelements.hh>
11

12
13
#include <dune/amdis/Initfile.hpp>
#include <dune/amdis/common/Loops.hpp>
14
#include <dune/amdis/utility/Filesystem.hpp>
15
16
17
18

namespace AMDiS
{

19
20
21
22
  template <class Traits>
  class FileWriter
  {
  private: // typedefs and static constants
23

24
    using MeshView = typename Traits::GridView;
25

26
    /// Dimension of the mesh
27
    static constexpr int dim = MeshView::dimension;
28

29
    /// Dimension of the world
30
    static constexpr int dow = MeshView::dimensionworld;
31

32
    /// Number of problem components
33
    static constexpr int nComponents = 1;
34
35


36
  public:
37

38
    /// Constructor.
39
    FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names)
40
      : meshView(meshView)
41
      , names(names)
42
    {
43
44
45
46
      filename = "solution";
      Parameters::get(base + "->filename", filename);
      dir = "output";
      Parameters::get(base + "->output directory", dir);
47

48
49
50
51
52
      if (!filesystem::exists(dir))
        error_exit("Output directory '",dir,"' does not exist!");

      vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
      vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
53
54
    }

55

56
    /// default write method for time-depended data
57
    template <class SystemVectorType>
58
    void write(double time, SystemVectorType const& solutions)
59
    {
60
61
      vtkWriter->clear();
      // copy dofvector to vertex data
62
      forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
63
      {
64
        this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
65
66
67
68
        vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
      });
      vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
    }
69

70
71
72

    /// default write method for stationary data
    template <class SystemVectorType>
73
    void write(SystemVectorType const& solutions)
74
75
76
    {
      vtkWriter->clear();
      // copy dofvector to vertex data
77
      forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
78
      {
79
        this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
80
81
82
        vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
      });
      vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
83
    }
84
85


86
  protected:
87

88
    template <class DOFVector, class Vector>
89
    void dofVector2vertexVector(DOFVector const& dofvector, Vector& data)
90
    {
91
92
      using Geometry = typename MeshView::template Codim<0>::Geometry;
      using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
93

94
      data.resize(meshView.size(dim));
95

96
97
      auto const& indexSet = meshView.indexSet();
      auto const& feSpace = dofvector.getFeSpace();
98

99
100
      auto localView = feSpace.localView();
      auto localIndexSet = feSpace.localIndexSet();
101
#if 0
102
103
104
105
106
      // copy data to P1-vector
      for (auto const& element : elements(meshView)) {
        localView.bind(element);
        localIndexSet.bind(localView);
        auto const& localBasis = localView.tree().finiteElement().localBasis();
107

108
        auto const& refElement = RefElements::general(element.type());
109

110
        std::vector<Dune::FieldVector<double,1> > shapeValues;
111

112
113
        std::size_t nVertices = element.subEntities(dim);
        for (std::size_t i = 0; i < nVertices; ++i) {
114
115
          auto const& v = element.template subEntity<dim>(i);
          auto pos = refElement.position(i, dim);
116

117
          localBasis.evaluateFunction(pos, shapeValues);
118

119
          std::size_t idx = indexSet.index(v);
120
          data[idx] = 0.0;
121
          for (std::size_t j = 0; j < shapeValues.size(); ++j) {
122
123
124
125
126
            const auto global_idx = localIndexSet.index(j);
            data[idx] += dofvector[global_idx] * shapeValues[j];
          }
        }
      }
127
#endif
128
129
    }

130
  private:
131
    MeshView meshView;
132
133
    shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
    shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
134

135
136
    std::vector<std::string> names;
    std::array<std::vector<double>, nComponents> data_vectors;
137

138
139
140
    std::string filename;
    std::string dir;
  };
141
142

} // end namespace AMDiS