FileWriter.hpp 4.27 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 25
    using Mesh     = typename Traits::Mesh;
    using MeshView = typename Mesh::LeafGridView;
26

27 28
    /// Dimension of the mesh
    static constexpr int dim = Traits::dim;
29

30 31
    /// Dimension of the world
    static constexpr int dow = Traits::dimworld;
32

33 34
    /// Number of problem components
    static constexpr int nComponents = Traits::nComponents;
35 36


37
  public:
38

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

49 50 51 52 53
      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, "");
54 55
    }

56

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

71 72 73

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


87
  protected:
88

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

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

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

100 101
      auto localView = feSpace.localView();
      auto localIndexSet = feSpace.localIndexSet();
102

103 104 105 106 107
      // copy data to P1-vector
      for (auto const& element : elements(meshView)) {
        localView.bind(element);
        localIndexSet.bind(localView);
        auto const& localBasis = localView.tree().finiteElement().localBasis();
108

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

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

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

118
          localBasis.evaluateFunction(pos, shapeValues);
119

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

130
  private:
131
    MeshView const& 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