dune-grid-glue  2.8-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 
8 namespace Dune {
9 namespace GridGlue {
10 
11 namespace AreaWriterImplementation {
12 
13 template<int dimgrid>
15 {
16  bool contains(Dune::GeometryType gt) const
17  {
18  return gt.dim() == dimgrid - 1;
19  }
20 };
21 
22 template<typename GridView>
23 void 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 
81 template<int side, typename Glue>
82 void 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 
119 template<int side, typename Glue>
120 void 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 
126 template<typename Glue>
127 void 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:35
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