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