4 #include <dune/common/fvector.hh> 5 #include <dune/geometry/type.hh> 6 #include <dune/grid/common/mcmgmapper.hh> 11 namespace AreaWriterImplementation {
18 return gt.dim() == dimgrid - 1;
22 template<
typename Gr
idView>
25 using Coordinate = Dune::FieldVector<double, 3>;
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) {
32 const auto c0 = geometry.corner(i);
34 for (
int d = 0; d < GridView::dimensionworld; ++d)
36 for (
int d = GridView::dimensionworld; d < Coordinate::dimension; ++d)
38 corners.push_back(c1);
43 out <<
"DATASET UNSTRUCTURED_GRID\n" 44 <<
"POINTS " << corners.size() <<
" double\n";
45 for (
const auto& c : corners)
49 out <<
"CELLS " << gv.size(1) <<
" " << (gv.size(1) + corners.size()) <<
"\n";
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)
60 out <<
"CELL_TYPES " << gv.size(1) <<
"\n";
61 for (
const auto& facet : facets(gv)) {
62 const auto type = facet.type();
65 else if (type.isLine())
67 else if (type.isTriangle())
69 else if (type.isQuadrilateral())
71 else if (type.isTetrahedron())
74 DUNE_THROW(Dune::Exception,
"Unhandled geometry type");
81 template<
int s
ide,
typename Glue>
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;
88 const GridView gv = glue.template gridView<side>();
90 std::vector<ctype> coveredArea(mapper.size(), ctype(0));
91 std::vector<ctype> totalArea(mapper.size(), ctype(1));
94 const auto element = in.inside();
95 const auto index = mapper.subIndex(element, in.indexInInside(), 1);
96 coveredArea[index] += in.geometryInInside().volume();
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();
103 for (std::size_t i = 0; i < coveredArea.size(); ++i)
104 coveredArea[i] /= totalArea[i];
106 out <<
"# vtk DataFile Version 2.0\n" 107 <<
"Filename: Glue Area\n" 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";
119 template<
int s
ide,
typename Glue>
122 std::ofstream out(filename.c_str());
123 write_glue_area_vtk<side>(glue, out);
126 template<
typename Glue>
130 std::string filename = base;
131 filename +=
"-inside.vtk";
132 write_glue_area_vtk<0>(glue, filename);
135 std::string filename = base;
136 filename +=
"-outside.vtk";
137 write_glue_area_vtk<1>(glue, filename);
bool contains(Dune::GeometryType gt) const
Definition: areawriter_impl.hh:16
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.
Definition: rangegenerators.hh:13
void write_glue_area_vtk(const Glue &glue, std::ostream &out)
Definition: areawriter_impl.hh:82
Definition: areawriter_impl.hh:14
Definition: gridglue.hh:30
void write_facet_geometry(const GridView &gv, std::ostream &out)
Definition: areawriter_impl.hh:23