18 #ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
19 #define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
53 typedef typename GV::Grid::ctype
ctype;
54 typedef Dune::FieldVector<ctype, dimworld>
Coords;
56 typedef typename GV::Traits::template Codim<dim>::EntityPointer
VertexPtr;
57 typedef typename GV::Traits::template Codim<0>::EntityPointer
ElementPtr;
59 static const Dune::PartitionIteratorType
PType = Dune::Interior_Partition;
60 typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator
ElementIter;
62 typedef typename GV::IntersectionIterator
IsIter;
83 std::cout <<
"This is Codim1Extractor on a <" <<
dim
109 template<
typename GV>
110 void Codim1Extractor<GV>::update(
const ExtractorPredicate<GV,1>& descr)
121 int simplex_index = 0;
122 int vertex_index = 0;
123 IndexType eindex = 0;
130 std::deque<SubEntityInfo> temp_faces;
133 for (ElementIter elit = this->gv_.template begin<0, PType>();
134 elit != this->gv_.template end<0, PType>(); ++elit)
136 ElementPtr eptr(elit);
137 Dune::GeometryType gt = elit->type();
140 if (elit->hasBoundaryIntersections())
144 eindex = this->cellMapper_.map(*elit);
145 this->elmtInfo_[eindex] =
new ElementInfo(simplex_index, eptr, 0);
149 for (IsIter is = this->gv_.ibegin(*elit); is != this->gv_.iend(*elit); ++is)
152 if (!is->boundary() or !descr.contains(eptr, is->indexInInside()))
155 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
156 const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
158 const Dune::GenericReferenceElement<ctype, dim>& refElement = Dune::GenericReferenceElements<ctype, dim>::general(gt);
161 const int face_corners = refElement.size(is->indexInInside(), 1, dim);
165 switch (face_corners)
173 this->elmtInfo_[eindex]->faces++;
176 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
177 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(is->indexInInside(), 1, i, dim);
188 VertexPtr vptr(elit->template subEntity<dim>(vertex_number));
189 cornerCoords[i] = vptr->geometry().corner(0);
191 IndexType vindex = this->gv_.indexSet().template index<dim>(*vptr);
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_[vindex] =
new VertexInfo(vertex_index, vptr);
204 temp_faces.back().corners[i].idx = vertex_index;
211 temp_faces.back().corners[i].idx = vimit->second->idx;
219 FieldVector<ctype,dimworld> realNormal = is->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]);
245 unsigned int vertex_indices[4];
246 unsigned int vertex_numbers[4];
249 this->elmtInfo_[eindex]->faces += 2;
251 Dune::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
255 for (
int i = 0; i < cube_corners; ++i)
258 vertex_numbers[i] = refElement.subEntity(is->indexInInside(), 1, i, dim);
261 VertexPtr vptr(elit->template subEntity<dim>(vertex_numbers[i]));
262 cornerCoords[i] = vptr->geometry().corner(0);
264 IndexType vindex = this->gv_.indexSet().template index<dim>(*vptr);
268 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
269 if (vimit == this->vtxInfo_.end())
272 this->vtxInfo_[vindex] =
new VertexInfo(vertex_index, vptr);
274 vertex_indices[i] = vertex_index;
281 vertex_indices[i] = vimit->second->idx;
290 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
291 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
292 temp_faces.back().corners[0].idx = vertex_indices[0];
293 temp_faces.back().corners[1].idx = vertex_indices[1];
294 temp_faces.back().corners[2].idx = vertex_indices[2];
296 temp_faces.back().corners[0].num = vertex_numbers[0];
297 temp_faces.back().corners[1].num = vertex_numbers[1];
298 temp_faces.back().corners[2].num = vertex_numbers[2];
304 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
308 cornerCoords[2] - cornerCoords[0]);
309 reconstructedNormal /= reconstructedNormal.two_norm();
311 if (realNormal * reconstructedNormal < 0.0)
312 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
316 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
317 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
318 temp_faces.back().corners[0].idx = vertex_indices[3];
319 temp_faces.back().corners[1].idx = vertex_indices[2];
320 temp_faces.back().corners[2].idx = vertex_indices[1];
322 temp_faces.back().corners[0].num = vertex_numbers[3];
323 temp_faces.back().corners[1].num = vertex_numbers[2];
324 temp_faces.back().corners[2].num = vertex_numbers[1];
332 cornerCoords[1] - cornerCoords[3]);
333 reconstructedNormal /= reconstructedNormal.two_norm();
335 if (realNormal * reconstructedNormal < 0.0)
336 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
342 DUNE_THROW(Dune::NotImplemented,
"the extractor does only work for triangle and quadrilateral faces (" << face_corners <<
" corners)");
349 std::cout <<
"added " << simplex_index <<
" subfaces\n";
352 this->subEntities_.resize(simplex_index);
354 copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
359 this->coords_.resize(this->vtxInfo_.size());
360 typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
361 for (; it1 != this->vtxInfo_.end(); ++it1)
364 CoordinateInfo* current = &this->coords_[it1->second->idx];
366 current->index = it1->second->idx;
368 current->vtxindex = it1->first;
371 current->coord = it1->second->p->geometry().corner(0);
380 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:69
GV::IntersectionIterator IsIter
Definition: codim1extractor.hh:62
Definition: gridglue.hh:34
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:66
Definition: codim1extractor.hh:48
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:54
GV::Traits::template Codim< dim >::EntityPointer VertexPtr
Definition: codim1extractor.hh:56
Definition: extractor.hh:54
GV::Grid::ctype ctype
Definition: codim1extractor.hh:53
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:30
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:67
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim1extractor.hh:57
Codim1Extractor(const GV &gv, const ExtractorPredicate< GV, 1 > &descr)
Constructor.
Definition: codim1extractor.hh:80
static Dune::FieldVector< T, dim > crossProduct(const Dune::FieldVector< T, dim > &a, const Dune::FieldVector< T, dim > &b)
Compute the cross product of two vectors.
Definition: projectionhelper.hh:19
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim1extractor.hh:60
Definition: extractor.hh:53
static const Dune::PartitionIteratorType PType
Definition: codim1extractor.hh:59
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:48
Definition: codim1extractor.hh:33
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:65
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:43
Base class for predicates selecting the part of a grid to be extracted.
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:68
GV GridView
Definition: codim1extractor.hh:51