6#include <dune/common/fvector.hh>
7#include <dune/geometry/type.hh>
8#include <dune/grid/common/mcmgmapper.hh>
13namespace AreaWriterImplementation {
20 return gt.dim() == dimgrid - 1;
24template<
typename Gr
idView>
27 using Coordinate = Dune::FieldVector<double, 3>;
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) {
34 const auto c0 = geometry.corner(i);
36 for (
int d = 0; d < GridView::dimensionworld; ++d)
38 for (
int d = GridView::dimensionworld; d < Coordinate::dimension; ++d)
40 corners.push_back(c1);
45 out <<
"DATASET UNSTRUCTURED_GRID\n"
46 <<
"POINTS " << corners.size() <<
" double\n";
47 for (
const auto& c : corners)
51 out <<
"CELLS " << gv.size(1) <<
" " << (gv.size(1) + corners.size()) <<
"\n";
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)
62 out <<
"CELL_TYPES " << gv.size(1) <<
"\n";
63 for (
const auto& facet : facets(gv)) {
64 const auto type = facet.type();
67 else if (type.isLine())
69 else if (type.isTriangle())
71 else if (type.isQuadrilateral())
73 else if (type.isTetrahedron())
76 DUNE_THROW(Dune::Exception,
"Unhandled geometry type");
83template<
int s
ide,
typename Glue>
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;
90 const GridView gv = glue.template gridView<side>();
92 std::vector<ctype> coveredArea(mapper.size(), ctype(0));
93 std::vector<ctype> totalArea(mapper.size(), ctype(1));
96 const auto element = in.inside();
97 const auto index = mapper.subIndex(element, in.indexInInside(), 1);
98 coveredArea[index] += in.geometryInInside().volume();
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();
105 for (std::size_t i = 0; i < coveredArea.size(); ++i)
106 coveredArea[i] /= totalArea[i];
108 out <<
"# vtk DataFile Version 2.0\n"
109 <<
"Filename: Glue Area\n"
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";
121template<
int s
ide,
typename Glue>
124 std::ofstream out(filename.c_str());
125 write_glue_area_vtk<side>(glue, out);
128template<
typename Glue>
132 std::string filename = base;
133 filename +=
"-inside.vtk";
134 write_glue_area_vtk<0>(glue, filename);
137 std::string filename = base;
138 filename +=
"-outside.vtk";
139 write_glue_area_vtk<1>(glue, filename);
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
Definition: areawriter_impl.hh:17
bool contains(Dune::GeometryType gt) const
Definition: areawriter_impl.hh:18