dune-grid-glue 2.5-git
areawriter_impl.hh
Go to the documentation of this file.
1#include <fstream>
2#include <vector>
3
4#include <dune/common/fvector.hh>
5#include <dune/geometry/type.hh>
6#include <dune/grid/common/mcmgmapper.hh>
7
8namespace Dune {
9namespace GridGlue {
10
11namespace AreaWriterImplementation {
12
13template<int dimgrid>
15{
16 bool contains(Dune::GeometryType gt) const
17 {
18 return gt.dim() == dimgrid - 1;
19 }
20};
21
22template<typename GridView>
23void write_facet_geometry(const GridView& gv, std::ostream& out)
24{
25 using Coordinate = Dune::FieldVector<double, 3>;
26
27 std::vector<Coordinate> corners;
28 for (const auto& facet : facets(gv)) {
29 const auto geometry = facet.geometry();
30 for (int i = 0; i < geometry.corners(); ++i) {
31 /* VTK always needs 3-dim coordinates... */
32 const auto c0 = geometry.corner(i);
33 Coordinate c1;
34 for (int d = 0; d < GridView::dimensionworld; ++d)
35 c1[d] = c0[d];
36 for (int d = GridView::dimensionworld; d < Coordinate::dimension; ++d)
37 c1[d] = double(0);
38 corners.push_back(c1);
39 }
40 }
41
42 {
43 out << "DATASET UNSTRUCTURED_GRID\n"
44 << "POINTS " << corners.size() << " double\n";
45 for (const auto& c : corners)
46 out << c << "\n";
47 }
48 {
49 out << "CELLS " << gv.size(1) << " " << (gv.size(1) + corners.size()) << "\n";
50 std::size_t c = 0;
51 for (const auto& facet : facets(gv)) {
52 const auto geometry = facet.geometry();
53 out << geometry.corners();
54 for (int i = 0; i < geometry.corners(); ++i, ++c)
55 out << " " << c;
56 out << "\n";
57 }
58 }
59 {
60 out << "CELL_TYPES " << gv.size(1) << "\n";
61 for (const auto& facet : facets(gv)) {
62 const auto type = facet.type();
63 if (type.isVertex())
64 out << "1\n";
65 else if (type.isLine())
66 out << "2\n";
67 else if (type.isTriangle())
68 out << "5\n";
69 else if (type.isQuadrilateral())
70 out << "9\n";
71 else if (type.isTetrahedron())
72 out << "10\n";
73 else
74 DUNE_THROW(Dune::Exception, "Unhandled geometry type");
75 }
76 }
77}
78
79} /* namespace AreaWriterImplementation */
80
81template<int side, typename Glue>
82void write_glue_area_vtk(const Glue& glue, std::ostream& out)
83{
84 using GridView = typename std::decay< decltype(glue.template gridView<side>()) >::type;
85 using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView, AreaWriterImplementation::FacetLayout>;
86 using ctype = typename GridView::ctype;
87
88 const GridView gv = glue.template gridView<side>();
89 Mapper mapper(gv);
90 std::vector<ctype> coveredArea(mapper.size(), ctype(0));
91 std::vector<ctype> totalArea(mapper.size(), ctype(1));
92
93 for (const auto& in : intersections(glue, Reverse<side == 1>())) {
94 const auto element = in.inside();
95 const auto index = mapper.subIndex(element, in.indexInInside(), 1);
96 coveredArea[index] += in.geometryInInside().volume();
97
98 const auto& refElement = Dune::ReferenceElements<ctype, GridView::dimension>::general(element.type());
99 const auto& subGeometry = refElement.template geometry<1>(in.indexInInside());
100 totalArea[index] = subGeometry.volume();
101 }
102
103 for (std::size_t i = 0; i < coveredArea.size(); ++i)
104 coveredArea[i] /= totalArea[i];
105
106 out << "# vtk DataFile Version 2.0\n"
107 << "Filename: Glue Area\n"
108 << "ASCII\n";
109
111
112 out << "CELL_DATA " << coveredArea.size() << "\n"
113 << "SCALARS CoveredArea double 1\n"
114 << "LOOKUP_TABLE default\n";
115 for (const auto& value : coveredArea)
116 out << value << "\n";
117}
118
119template<int side, typename Glue>
120void write_glue_area_vtk(const Glue& glue, const std::string& filename)
121{
122 std::ofstream out(filename.c_str());
123 write_glue_area_vtk<side>(glue, out);
124}
125
126template<typename Glue>
127void write_glue_areas_vtk(const Glue& glue, const std::string& base)
128{
129 {
130 std::string filename = base;
131 filename += "-inside.vtk";
132 write_glue_area_vtk<0>(glue, filename);
133 }
134 {
135 std::string filename = base;
136 filename += "-outside.vtk";
137 write_glue_area_vtk<1>(glue, filename);
138 }
139}
140
141} /* namespace GridGlue */
142} /* namespace Dune */
Definition gridglue.hh:30
void write_glue_area_vtk(const Glue &glue, std::ostream &out)
Definition areawriter_impl.hh:82
void write_glue_areas_vtk(const Glue &glue, const std::string &base)
Definition areawriter_impl.hh:127
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
void write_facet_geometry(const GridView &gv, std::ostream &out)
Definition areawriter_impl.hh:23
Definition rangegenerators.hh:15
bool contains(Dune::GeometryType gt) const
Definition areawriter_impl.hh:16