20 #ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
21 #define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
29 #include <dune/common/deprecated.hh>
30 #include <dune/common/version.hh>
58 typedef typename GV::Grid::ctype
ctype;
59 typedef Dune::FieldVector<ctype, dimworld>
Coords;
61 typedef typename GV::Traits::template Codim<dim>::Entity
Vertex;
62 typedef typename GV::Traits::template Codim<0>::Entity
Element;
84 std::cout <<
"This is Codim1Extractor on a <" <<
dim
110 template<
typename GV>
111 void Codim1Extractor<GV>::update(
const Predicate& predicate)
122 int simplex_index = 0;
123 int vertex_index = 0;
124 IndexType eindex = 0;
131 std::deque<SubEntityInfo> temp_faces;
134 for (
const auto& elmt : elements(this->gv_, Partitions::interior))
136 Dune::GeometryType gt = elmt.type();
139 if (elmt.hasBoundaryIntersections())
143 eindex = this->cellMapper_.index(elmt);
144 this->elmtInfo_.emplace(eindex, ElementInfo(simplex_index, elmt, 0));
151 if (!in.boundary() or !predicate(elmt, in.indexInInside()))
154 const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
156 const int face_corners = refElement.size(in.indexInInside(), 1, dim);
160 switch (face_corners)
168 this->elmtInfo_.at(eindex).faces++;
171 temp_faces.emplace_back(eindex, in.indexInInside(),
172 #
if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
173 Dune::GeometryTypes::simplex(dim-codim)
175 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
179 std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
182 for (
int i = 0; i < face_corners; ++i)
185 int vertex_number = refElement.subEntity(in.indexInInside(), 1, i, dim);
188 const Vertex vertex = elmt.template subEntity<dim>(vertex_number);
189 cornerCoords[i] = vertex.geometry().corner(0);
191 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
194 temp_faces.back().corners[i].num = vertex_number;
198 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
199 if (vimit == this->vtxInfo_.end())
202 this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
204 temp_faces.back().corners[i].idx = vertex_index;
211 temp_faces.back().corners[i].idx = vimit->second.idx;
219 FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
222 FieldVector<ctype,dimworld> reconstructedNormal;
225 reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
226 reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
228 FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
229 FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
232 reconstructedNormal /= reconstructedNormal.two_norm();
234 if (realNormal * reconstructedNormal < 0.0)
235 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
243 assert(dim == 3 && cube_corners == 4);
245 std::array<unsigned int, 4> vertex_indices;
246 std::array<unsigned int, 4> vertex_numbers;
249 this->elmtInfo_.at(eindex).faces += 2;
251 std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
255 for (
int i = 0; i < cube_corners; ++i)
258 vertex_numbers[i] = refElement.subEntity(in.indexInInside(), 1, i, dim);
261 const Vertex vertex = elmt.template subEntity<dim>(vertex_numbers[i]);
262 cornerCoords[i] = vertex.geometry().corner(0);
264 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
268 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
269 if (vimit == this->vtxInfo_.end())
272 this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
274 vertex_indices[i] = vertex_index;
281 vertex_indices[i] = vimit->second.idx;
290 temp_faces.emplace_back(eindex, in.indexInInside(),
291 #
if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
292 Dune::GeometryTypes::simplex(dim-codim)
294 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
297 temp_faces.back().corners[0].idx = vertex_indices[0];
298 temp_faces.back().corners[1].idx = vertex_indices[1];
299 temp_faces.back().corners[2].idx = vertex_indices[2];
301 temp_faces.back().corners[0].num = vertex_numbers[0];
302 temp_faces.back().corners[1].num = vertex_numbers[1];
303 temp_faces.back().corners[2].num = vertex_numbers[2];
309 FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
312 FieldVector<ctype,dimworld> reconstructedNormal =
crossProduct(cornerCoords[1] - cornerCoords[0],
313 cornerCoords[2] - cornerCoords[0]);
314 reconstructedNormal /= reconstructedNormal.two_norm();
316 if (realNormal * reconstructedNormal < 0.0)
317 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
321 temp_faces.emplace_back(eindex, in.indexInInside(),
322 #
if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
323 Dune::GeometryTypes::simplex(dim-codim)
325 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
328 temp_faces.back().corners[0].idx = vertex_indices[3];
329 temp_faces.back().corners[1].idx = vertex_indices[2];
330 temp_faces.back().corners[2].idx = vertex_indices[1];
332 temp_faces.back().corners[0].num = vertex_numbers[3];
333 temp_faces.back().corners[1].num = vertex_numbers[2];
334 temp_faces.back().corners[2].num = vertex_numbers[1];
341 reconstructedNormal =
crossProduct(cornerCoords[2] - cornerCoords[3],
342 cornerCoords[1] - cornerCoords[3]);
343 reconstructedNormal /= reconstructedNormal.two_norm();
345 if (realNormal * reconstructedNormal < 0.0)
346 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
352 DUNE_THROW(Dune::NotImplemented,
"the extractor does only work for triangle and quadrilateral faces (" << face_corners <<
" corners)");
359 std::cout <<
"added " << simplex_index <<
" subfaces\n";
362 this->subEntities_.resize(simplex_index);
364 copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
369 this->coords_.resize(this->vtxInfo_.size());
370 for (
const auto& vinfo : this->vtxInfo_)
373 CoordinateInfo* current = &this->coords_[vinfo.second.idx];
375 current->index = vinfo.second.idx;
377 current->vtxindex = vinfo.first;
380 const auto vtx = this->grid().entity(vinfo.second.p);
381 current->coord = vtx.geometry().corner(0);
Definition: gridglue.hh:37
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
compute cross product
Definition: crossproduct.hh:15
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
Definition: codim1extractor.hh:42
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:51
GV GridView
Definition: codim1extractor.hh:56
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim1extractor.hh:62
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:59
GV::Grid::ctype ctype
Definition: codim1extractor.hh:58
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim1extractor.hh:81
static constexpr int simplex_corners
compile time number of corners of surface simplices
Definition: codim1extractor.hh:54
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:68
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:69
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:67
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:66
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:70
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim1extractor.hh:63
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim1extractor.hh:61
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:46
static constexpr auto dimworld
Definition: extractor.hh:50
static constexpr auto dim
Definition: extractor.hh:51