18 #ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
19 #define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
27 #include <dune/common/deprecated.hh>
55 typedef typename GV::Grid::ctype
ctype;
56 typedef Dune::FieldVector<ctype, dimworld>
Coords;
58 typedef typename GV::Traits::template Codim<dim>::Entity
Vertex;
59 typedef typename GV::Traits::template Codim<0>::Entity
Element;
62 static const Dune::PartitionIteratorType
PType = Dune::Interior_Partition;
63 typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator
ElementIter;
65 typedef typename GV::IntersectionIterator
IsIter;
83 DUNE_DEPRECATED_MSG(
"Please use a std::function<bool(const Element&, unsigned int)> in favor of the ExtractorPredicate.")
87 std::cout <<
"This is Codim1Extractor on a <" <<
dim
90 const auto predicate = [&](
const Element&
element,
unsigned int subentity) ->
bool {
91 return descr.contains(
element, subentity);
104 std::cout <<
"This is Codim1Extractor on a <" <<
dim
130 template<
typename GV>
131 void Codim1Extractor<GV>::update(
const Predicate& predicate)
142 int simplex_index = 0;
143 int vertex_index = 0;
144 IndexType eindex = 0;
151 std::deque<SubEntityInfo> temp_faces;
154 for (ElementIter elit = this->gv_.template begin<0, PType>();
155 elit != this->gv_.template end<0, PType>(); ++elit)
157 const auto& elmt = *elit;
158 Dune::GeometryType gt = elmt.type();
161 if (elit->hasBoundaryIntersections())
165 eindex = this->cellMapper_.index(elmt);
166 this->elmtInfo_[eindex] =
new ElementInfo(simplex_index, elmt, 0);
170 for (IsIter is = this->gv_.ibegin(elmt); is != this->gv_.iend(elmt); ++is)
173 if (!is->boundary() or !predicate(elmt, is->indexInInside()))
176 const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
178 const int face_corners = refElement.size(is->indexInInside(), 1, dim);
182 switch (face_corners)
190 this->elmtInfo_[eindex]->faces++;
193 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
194 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
196 std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
199 for (
int i = 0; i < face_corners; ++i)
202 int vertex_number = refElement.subEntity(is->indexInInside(), 1, i, dim);
205 const Vertex vertex = elit->template subEntity<dim>(vertex_number);
206 cornerCoords[i] = vertex.geometry().corner(0);
208 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
211 temp_faces.back().corners[i].num = vertex_number;
215 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
216 if (vimit == this->vtxInfo_.end())
219 this->vtxInfo_[vindex] =
new VertexInfo(vertex_index, vertex);
221 temp_faces.back().corners[i].idx = vertex_index;
228 temp_faces.back().corners[i].idx = vimit->second->idx;
236 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
239 FieldVector<ctype,dimworld> reconstructedNormal;
242 reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
243 reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
245 FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
246 FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
249 reconstructedNormal /= reconstructedNormal.two_norm();
251 if (realNormal * reconstructedNormal < 0.0)
252 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
262 unsigned int vertex_indices[4];
263 unsigned int vertex_numbers[4];
266 this->elmtInfo_[eindex]->faces += 2;
268 std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
272 for (
int i = 0; i < cube_corners; ++i)
275 vertex_numbers[i] = refElement.subEntity(is->indexInInside(), 1, i, dim);
278 const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
279 cornerCoords[i] = vertex.geometry().corner(0);
281 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
285 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
286 if (vimit == this->vtxInfo_.end())
289 this->vtxInfo_[vindex] =
new VertexInfo(vertex_index, vertex);
291 vertex_indices[i] = vertex_index;
298 vertex_indices[i] = vimit->second->idx;
307 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
308 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
309 temp_faces.back().corners[0].idx = vertex_indices[0];
310 temp_faces.back().corners[1].idx = vertex_indices[1];
311 temp_faces.back().corners[2].idx = vertex_indices[2];
313 temp_faces.back().corners[0].num = vertex_numbers[0];
314 temp_faces.back().corners[1].num = vertex_numbers[1];
315 temp_faces.back().corners[2].num = vertex_numbers[2];
321 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
324 FieldVector<ctype,dimworld> reconstructedNormal =
crossProduct(cornerCoords[1] - cornerCoords[0],
325 cornerCoords[2] - cornerCoords[0]);
326 reconstructedNormal /= reconstructedNormal.two_norm();
328 if (realNormal * reconstructedNormal < 0.0)
329 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
333 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
334 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
335 temp_faces.back().corners[0].idx = vertex_indices[3];
336 temp_faces.back().corners[1].idx = vertex_indices[2];
337 temp_faces.back().corners[2].idx = vertex_indices[1];
339 temp_faces.back().corners[0].num = vertex_numbers[3];
340 temp_faces.back().corners[1].num = vertex_numbers[2];
341 temp_faces.back().corners[2].num = vertex_numbers[1];
348 reconstructedNormal =
crossProduct(cornerCoords[2] - cornerCoords[3],
349 cornerCoords[1] - cornerCoords[3]);
350 reconstructedNormal /= reconstructedNormal.two_norm();
352 if (realNormal * reconstructedNormal < 0.0)
353 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
359 DUNE_THROW(Dune::NotImplemented,
"the extractor does only work for triangle and quadrilateral faces (" << face_corners <<
" corners)");
366 std::cout <<
"added " << simplex_index <<
" subfaces\n";
369 this->subEntities_.resize(simplex_index);
371 copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
376 this->coords_.resize(this->vtxInfo_.size());
377 typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
378 for (; it1 != this->vtxInfo_.end(); ++it1)
381 CoordinateInfo* current = &this->coords_[it1->second->idx];
383 current->index = it1->second->idx;
385 current->vtxindex = it1->first;
388 const auto vtx = this->grid().entity(it1->second->p);
389 current->coord = vtx.geometry().corner(0);
Base class for predicates selecting the part of a grid to be extracted.
Definition: gridglue.hh:30
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
compute cross product
Definition: crossproduct.hh:13
Definition: codim1extractor.hh:36
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:45
GV GridView
Definition: codim1extractor.hh:53
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim1extractor.hh:59
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:56
GV::Grid::ctype ctype
Definition: codim1extractor.hh:55
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim1extractor.hh:101
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:70
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:71
@ simplex_corners
Definition: codim1extractor.hh:50
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:69
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim1extractor.hh:63
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:68
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:72
static const Dune::PartitionIteratorType PType
Definition: codim1extractor.hh:62
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim1extractor.hh:60
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim1extractor.hh:58
GV::IntersectionIterator IsIter
Definition: codim1extractor.hh:65
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:43
Element element(unsigned int index) const
gets the parent element for a given face index, throws an exception if index not valid
Definition: extractor.hh:375
@ dim
Definition: extractor.hh:48
@ dimworld
Definition: extractor.hh:47
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:31