FileWriter.hpp 3.21 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
56
57
58
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
#pragma once

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

#include <dune/grid/io/file/vtk/vtksequencewriter.hh>
#include <dune/geometry/genericgeometry/referenceelements.hh>

#include "Loops.hpp"
#include "Initfile.hpp"

namespace AMDiS
{

template <class Traits>
class FileWriter
{
    using Mesh     = typename Traits::Mesh;
    using MeshView = typename Mesh::LeafGridView;
    
    /// Dimension of the mesh
    static constexpr int dim = Traits::dim;
    
    /// Dimension of the world
    static constexpr int dow = Traits::dimworld;
    
    /// Number of problem components
    static constexpr int nComponents = Traits::nComponents;
    
    
public:
    FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names_)
	: meshView(meshView)
	, names(names_)
    {
	std::string filename = "solution";
	Parameters::get(base + "->filename", filename);
	std::string dir = "output";
	Parameters::get(base + "->output directory", dir);
	
	vtkWriter = std::make_shared<Dune::VTKWriter<MeshView>>(meshView);
	vtkSeqWriter = std::make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
    }

    template <class SystemVectorType>
    void write(double time, SystemVectorType&& solutions)
    {
	vtkWriter->clear();
	// copy dofvector to vertex data
	For<0, nComponents>::loop([this, &solutions](const auto _i) 
	{
	    this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
	    vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
	});
	vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
    }
    
    
protected:
    template <class DOFVector, class Vector>
    void dofVector2vertexVector(DOFVector dofvector, Vector& data)
    {
	using Geometry = typename MeshView::template Codim<0>::Geometry;
	using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
	
	data.resize(meshView.size(dim));
	
	auto const& indexSet = meshView.indexSet();
	auto const& feSpace = dofvector.getFeSpace();
	
	auto localView = feSpace.localView();
	auto localIndexSet = feSpace.localIndexSet();
	
	// copy data to P1-vector
	for (auto const& element : elements(meshView)) {
	    localView.bind(element);
	    localIndexSet.bind(localView);
	    auto const& localBasis = localView.tree().finiteElement().localBasis();
	    
	    auto const& refElement = RefElements::general(element.type());
	    
	    std::vector<Dune::FieldVector<double,1> > shapeValues;
	    
	    size_t nVertices = element.subEntities(dim);
	    for (size_t i = 0; i < nVertices; ++i) {
		auto const& v = element.template subEntity<dim>(i);
		auto pos = refElement.position(i, dim);
		
		localBasis.evaluateFunction(pos, shapeValues);
		
		size_t idx = indexSet.index(v);
		data[idx] = 0.0;
		for (size_t j = 0; j < shapeValues.size(); ++j) {
		    const auto global_idx = localIndexSet.index(j);
		    data[idx] += dofvector[global_idx] * shapeValues[j];
		}
	    }
	}
    }

private:
    MeshView const& meshView;
    std::shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
    std::shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
    
    std::vector<std::string> names;
    std::array<std::vector<double>, nComponents> data_vectors;
};

} // end namespace AMDiS