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>::EntityPointer
VertexPtr;
59 typedef typename GV::Traits::template Codim<dim>::Entity
Vertex;
60 typedef typename GV::Traits::template Codim<0>::EntityPointer
ElementPtr;
61 typedef typename GV::Traits::template Codim<0>::Entity
Element;
62 typedef std::function<bool(const Element&, unsigned int subentity)>
Predicate;
64 static const Dune::PartitionIteratorType
PType = Dune::Interior_Partition;
65 typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator
ElementIter;
67 typedef typename GV::IntersectionIterator
IsIter;
85 DUNE_DEPRECATED_MSG(
"Please use a std::function<bool(const Element&, unsigned int)> in favor of the ExtractorPredicate.")
89 std::cout <<
"This is Codim1Extractor on a <" <<
dim 92 const auto predicate = [&](
const Element&
element,
unsigned int subentity) ->
bool {
93 return descr.contains(element, subentity);
106 std::cout <<
"This is Codim1Extractor on a <" <<
dim 127 void update(
const Predicate& predicate);
132 template<
typename GV>
144 int simplex_index = 0;
145 int vertex_index = 0;
153 std::deque<SubEntityInfo> temp_faces;
157 elit != this->gv_.template end<0, PType>(); ++elit)
159 const auto& elmt = *elit;
160 Dune::GeometryType gt = elmt.type();
163 if (elit->hasBoundaryIntersections())
167 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) 176 for (
IsIter is = this->
gv_.ibegin(elmt); is != this->
gv_.iend(elmt); ++is)
179 if (!is->boundary() or !predicate(elmt, is->indexInInside()))
182 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3) 183 const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
185 const Dune::GenericReferenceElement<ctype, dim>& refElement = Dune::GenericReferenceElements<ctype, dim>::general(gt);
188 const int face_corners = refElement.size(is->indexInInside(), 1,
dim);
192 switch (face_corners)
203 temp_faces.push_back(
SubEntityInfo(eindex, is->indexInInside(),
204 Dune::GeometryType(Dune::GeometryType::simplex,
dim-
codim)));
206 std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
209 for (
int i = 0; i < face_corners; ++i)
212 int vertex_number = refElement.subEntity(is->indexInInside(), 1, i,
dim);
215 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) 216 const Vertex vertex = elit->template subEntity<dim>(vertex_number);
218 const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_number);
219 const Vertex& vertex = *vertexPtr;
221 cornerCoords[i] = vertex.geometry().corner(0);
226 temp_faces.back().corners[i].num = vertex_number;
230 typename VertexInfoMap::iterator vimit = this->
vtxInfo_.find(vindex);
236 temp_faces.back().corners[i].idx = vertex_index;
243 temp_faces.back().corners[i].idx = vimit->second->idx;
251 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
254 FieldVector<ctype,dimworld> reconstructedNormal;
257 reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
258 reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
260 FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
261 FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
264 reconstructedNormal /= reconstructedNormal.two_norm();
266 if (realNormal * reconstructedNormal < 0.0)
267 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
277 unsigned int vertex_indices[4];
278 unsigned int vertex_numbers[4];
283 std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
290 vertex_numbers[i] = refElement.subEntity(is->indexInInside(), 1, i,
dim);
293 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) 294 const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
296 const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_numbers[i]);
297 const Vertex &vertex = *vertexPtr;
299 cornerCoords[i] = vertex.geometry().corner(0);
305 typename VertexInfoMap::iterator vimit = this->
vtxInfo_.find(vindex);
311 vertex_indices[i] = vertex_index;
318 vertex_indices[i] = vimit->second->idx;
327 temp_faces.push_back(
SubEntityInfo(eindex, is->indexInInside(),
328 Dune::GeometryType(Dune::GeometryType::simplex,
dim-
codim)));
329 temp_faces.back().corners[0].idx = vertex_indices[0];
330 temp_faces.back().corners[1].idx = vertex_indices[1];
331 temp_faces.back().corners[2].idx = vertex_indices[2];
333 temp_faces.back().corners[0].num = vertex_numbers[0];
334 temp_faces.back().corners[1].num = vertex_numbers[1];
335 temp_faces.back().corners[2].num = vertex_numbers[2];
341 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
344 FieldVector<ctype,dimworld> reconstructedNormal =
crossProduct(cornerCoords[1] - cornerCoords[0],
345 cornerCoords[2] - cornerCoords[0]);
346 reconstructedNormal /= reconstructedNormal.two_norm();
348 if (realNormal * reconstructedNormal < 0.0)
349 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
353 temp_faces.push_back(
SubEntityInfo(eindex, is->indexInInside(),
354 Dune::GeometryType(Dune::GeometryType::simplex,
dim-
codim)));
355 temp_faces.back().corners[0].idx = vertex_indices[3];
356 temp_faces.back().corners[1].idx = vertex_indices[2];
357 temp_faces.back().corners[2].idx = vertex_indices[1];
359 temp_faces.back().corners[0].num = vertex_numbers[3];
360 temp_faces.back().corners[1].num = vertex_numbers[2];
361 temp_faces.back().corners[2].num = vertex_numbers[1];
368 reconstructedNormal =
crossProduct(cornerCoords[2] - cornerCoords[3],
369 cornerCoords[1] - cornerCoords[3]);
370 reconstructedNormal /= reconstructedNormal.two_norm();
372 if (realNormal * reconstructedNormal < 0.0)
373 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
379 DUNE_THROW(Dune::NotImplemented,
"the extractor does only work for triangle and quadrilateral faces (" << face_corners <<
" corners)");
386 std::cout <<
"added " << simplex_index <<
" subfaces\n";
391 copy(temp_faces.begin(), temp_faces.end(), this->
subEntities_.begin());
397 typename VertexInfoMap::const_iterator it1 = this->
vtxInfo_.begin();
398 for (; it1 != this->
vtxInfo_.end(); ++it1)
403 current->index = it1->second->idx;
405 current->vtxindex = it1->first;
408 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4) 409 const auto vtx = this->
grid().entity(it1->second->p);
411 const auto vtxPtr = this->
grid().entityPointer(it1->second->p);
412 const auto& vtx = *vtxPtr;
414 current->coord = vtx.geometry().corner(0);
423 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH static const Dune::PartitionIteratorType PType
Definition: codim1extractor.hh:64
GV::Grid::ctype ctype
Definition: codim1extractor.hh:55
const GridView gv_
the grid object to extract the surface from
Definition: extractor.hh:206
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim1extractor.hh:62
GV::IntersectionIterator IsIter
Definition: codim1extractor.hh:67
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:71
void clear()
delete everything build up so far and free the memory
Definition: extractor.hh:252
std::vector< CoordinateInfo > coords_
all information about the corner vertices of the extracted
Definition: extractor.hh:211
std::vector< SubEntityInfo > subEntities_
all information about the extracted subEntities
Definition: extractor.hh:214
Definition: gridglue.hh:33
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:72
Definition: extractor.hh:54
Definition: codim1extractor.hh:50
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:396
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:56
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim1extractor.hh:103
const Grid & grid() const
Definition: extractor.hh:380
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:70
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:74
Definition: extractor.hh:53
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim1extractor.hh:59
GV::Traits::template Codim< dim >::EntityPointer VertexPtr
Definition: codim1extractor.hh:58
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:45
GV GridView
Definition: codim1extractor.hh:53
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:30
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:73
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim1extractor.hh:65
VertexInfoMap vtxInfo_
a map enabling faster access to vertices and coordinates
Definition: extractor.hh:221
Definition: extractor.hh:52
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim1extractor.hh:61
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim1extractor.hh:60
Base class for predicates selecting the part of a grid to be extracted.
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
CellMapper cellMapper_
Definition: extractor.hh:230
Vertex vertex(unsigned int index) const
gets the vertex for a given coordinate index throws an exception if index not valid ...
Definition: extractor.hh:420
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:47
Definition: extractor.hh:58
Definition: codim1extractor.hh:35
ElementInfoMap elmtInfo_
a map enabling faster access to elements and faces
Definition: extractor.hh:228