FileWriter.hpp 6.66 KB
Newer Older
1
2
3
#pragma once

#include <string>
4
5
//#include <vector>
//#include <array>
6
7
#include <memory>

8
#include <dune/functions/functionspacebases/lagrangebasis.hh>
9
#include <dune/grid/io/file/vtk/vtkwriter.hh>
10
#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
11
12
//#include <dune/geometry/referenceelements.hh>
#include <dune/typetree/childextraction.hh>
13

14
#include <dune/amdis/Initfile.hpp>
15
#include <dune/amdis/io/FileWriterInterface.hpp>
16
#include <dune/amdis/utility/Filesystem.hpp>
17
18
19

namespace AMDiS
{
20
  template <class GlobalBasis, class TreePath>
21
  class FileWriter
22
      : public FileWriterInterface
23
24
  {
  private: // typedefs and static constants
25

26
    using GridView = typename GlobalBasis::GridView;
27

28
    /// Dimension of the mesh
29
    static constexpr int dim = GridView::dimension;
30

31
    /// Dimension of the world
32
    static constexpr int dow = GridView::dimensionworld;
33

34
    using Vector = DOFVector<GlobalBasis>;
35
36
37



38
  public:
39
    /// Constructor.
40
41
42
43
44
45
46
47
    FileWriter(std::string baseName,
               std::shared_ptr<GlobalBasis> const& basis,
               TreePath const& tp,
               std::shared_ptr<Vector> const& vector)
      : FileWriterInterface(baseName)
      , basis_(basis)
      , treePath_(tp)
      , vector_(vector)
48
    {
49
      //int subSampling = Parameters::get<int>(baseName + "->subsampling").value_or(0);
50

51
52
53
54
      //if (subSampling > 0)
      //  vtkWriter_ = std::make_shared<Dune::SubsamplingVTKWriter<GridView>>(basis_->gridView(), subSampling);
      //else
      vtkWriter_ = std::make_shared<Dune::VTKWriter<GridView>>(basis_->gridView());
55

56
57
58
      vtkSeqWriter_ = std::make_shared<Dune::VTKSequenceWriter<GridView>>(vtkWriter_, filename_, dir_, "");

      mode_ = Parameters::get<int>(baseName + "->ParaView mode").value_or(0);
59
60
    }

61

62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    /// Implements \ref FileWriterInterface::writeFiles
    virtual void writeFiles(AdaptInfo& adaptInfo, bool force) override
    {
      using Tree = typename GlobalBasis::LocalView::Tree;
      using Node = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;

      writeVertexData(typename Node::NodeTag{}, index_<Node::CHILDREN>, [&,this]()
      {
        if (mode_ == 0)
          vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::ascii);
        else if (mode_ == 1)
          vtkSeqWriter_->write(adaptInfo.getTime(), Dune::VTK::appendedraw);
      });
    }

  protected:
    template <class W>
    void writeVertexData(Dune::TypeTree::LeafNodeTag, index_t<0>, W write)
    {
      using namespace Dune::Functions::BasisBuilder;
      auto fct = makeDiscreteFunction(basis_,treePath_,vector_);

      auto p1basis = makeBasis(basis_->gridView(), lagrange<1>());
      auto data = makeDOFVector(p1basis, name_);
      interpolate(p1basis, data, *fct);

      auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<double>(p1basis,data);

      vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::scalar, 1));
      write();
      vtkWriter_->clear();
    }

    template <std::size_t C, class W>
    void writeVertexData(Dune::TypeTree::PowerNodeTag, index_t<C>, W write)
    {
      using namespace Dune::Functions::BasisBuilder;
      assert( C == dow );
      auto fct = makeDiscreteFunction(basis_,treePath_,vector_);

      auto p1basis = makeBasis(basis_->gridView(), power<C>(lagrange<1>(), flatLexicographic()));
      auto data = makeDOFVector(p1basis, name_);
      interpolate(p1basis, data, *fct);

      using Range = Dune::FieldVector<double,C>;
      auto dataFct = Dune::Functions::makeDiscreteGlobalBasisFunction<Range>(p1basis,data);

      vtkWriter_->addVertexData(dataFct, Dune::VTK::FieldInfo(name_, Dune::VTK::FieldInfo::Type::vector, C));
      write();
      vtkWriter_->clear();
    }

    template <class NodeTag, std::size_t C, class W>
    void writeVertexData(NodeTag, index_t<C>, W) {}

#if 0
118
    /// default write method for time-depended data
119
    template <class SystemVectorType>
120
    void write(double time, SystemVectorType const& solutions)
121
    {
122
123
      vtkWriter->clear();
      // copy dofvector to vertex data
124
      forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
125
      {
126
        this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
127
128
129
130
        vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
      });
      vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
    }
131

132
133
134

    /// default write method for stationary data
    template <class SystemVectorType>
135
    void write(SystemVectorType const& solutions)
136
137
138
    {
      vtkWriter->clear();
      // copy dofvector to vertex data
139
      forEach(range_<0, nComponents>, [this, &solutions](const auto _i)
140
      {
141
        this->dofVector2vertexVector(solutions, std::get<_i>(data_vectors));
142
143
144
        vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
      });
      vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
145
    }
146

147
    template <class DOFVector, class Vector>
148
    void dofVector2vertexVector(DOFVector const& dofvector, Vector& data)
149
    {
150
151
      using Geometry = typename MeshView::template Codim<0>::Geometry;
      using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
152

153
      data.resize(meshView.size(dim));
154

155
156
      auto const& indexSet = meshView.indexSet();
      auto const& feSpace = dofvector.getFeSpace();
157

158
159
      auto localView = feSpace.localView();
      auto localIndexSet = feSpace.localIndexSet();
160

161
162
163
164
165
      // copy data to P1-vector
      for (auto const& element : elements(meshView)) {
        localView.bind(element);
        localIndexSet.bind(localView);
        auto const& localBasis = localView.tree().finiteElement().localBasis();
166

167
        auto const& refElement = RefElements::general(element.type());
168

169
        std::vector<Dune::FieldVector<double,1> > shapeValues;
170

171
172
        std::size_t nVertices = element.subEntities(dim);
        for (std::size_t i = 0; i < nVertices; ++i) {
173
174
          auto const& v = element.template subEntity<dim>(i);
          auto pos = refElement.position(i, dim);
175

176
          localBasis.evaluateFunction(pos, shapeValues);
177

178
          std::size_t idx = indexSet.index(v);
179
          data[idx] = 0.0;
180
          for (std::size_t j = 0; j < shapeValues.size(); ++j) {
181
182
183
184
185
            const auto global_idx = localIndexSet.index(j);
            data[idx] += dofvector[global_idx] * shapeValues[j];
          }
        }
      }
186
    }
187
#endif
188

189
  private:
190
191
192
193
194
195
    std::shared_ptr<GlobalBasis> basis_;
    TreePath treePath_;
    std::shared_ptr<Vector> vector_;

    std::shared_ptr<Dune::VTKWriter<GridView>> vtkWriter_;
    std::shared_ptr<Dune::VTKSequenceWriter<GridView>> vtkSeqWriter_;
196

197
    // std::vector<double> data_vector;
198

199
200
    // represents VTK::OutputType: 0...ascii, 1...appendedraw
    int mode_;
201
  };
202
203

} // end namespace AMDiS