dune-grid-glue 2.5-git
codim1extractor.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3/*
4 * Filename: codim1extractor.hh
5 * Version: 1.0
6 * Created on: Jun 23, 2009
7 * Author: Oliver Sander, Christian Engwer
8 * ---------------------------------
9 * Project: dune-grid-glue
10 * Description: class for grid extractors extracting surface grids
11 *
12 */
18#ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
19#define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
20
21#include "extractor.hh"
22#include "extractorpredicate.hh"
23
24#include <deque>
25#include <functional>
26
27#include <dune/common/deprecated.hh>
29
30namespace Dune {
31
32 namespace GridGlue {
33
34template<typename GV>
35class Codim1Extractor : public Extractor<GV,1>
36{
37public:
38
39 /* E X P O R T E D T Y P E S A N D C O N S T A N T S */
40
41 using Extractor<GV,1>::dimworld;
42 using Extractor<GV,1>::dim;
43 using Extractor<GV,1>::codim;
44 using Extractor<GV,1>::cube_corners;
46
48 enum
49 {
51 };
52
53 typedef GV GridView;
54
55 typedef typename GV::Grid::ctype ctype;
56 typedef Dune::FieldVector<ctype, dimworld> Coords;
57
58 typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
59 typedef typename GV::Traits::template Codim<0>::Entity Element;
60 typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
61
62 static const Dune::PartitionIteratorType PType = Dune::Interior_Partition;
63 typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator ElementIter;
64
65 typedef typename GV::IntersectionIterator IsIter;
66
67 // import typedefs from base class
73
74public:
75
76 /* C O N S T R U C T O R S A N D D E S T R U C T O R S */
77
83 DUNE_DEPRECATED_MSG("Please use a std::function<bool(const Element&, unsigned int)> in favor of the ExtractorPredicate.")
84 Codim1Extractor(const GV& gv, const ExtractorPredicate<GV,1>& descr)
85 : Extractor<GV,1>(gv)
86 {
87 std::cout << "This is Codim1Extractor on a <" << dim
88 << "," << dimworld << "> grid!"
89 << std::endl;
90 const auto predicate = [&](const Element& element, unsigned int subentity) -> bool {
91 return descr.contains(element, subentity);
92 };
93 update(predicate);
94 }
95
101 Codim1Extractor(const GV& gv, const Predicate& predicate)
102 : Extractor<GV,1>(gv)
103 {
104 std::cout << "This is Codim1Extractor on a <" << dim
105 << "," << dimworld << "> grid!"
106 << std::endl;
107 update(predicate);
108 }
109
110private:
111
125 void update(const Predicate& predicate);
126
127};
128
129
130template<typename GV>
131void Codim1Extractor<GV>::update(const Predicate& predicate)
132{
133 // free everything there is in this object
134 this->clear();
135
136 // In this first pass iterate over all entities of codim 0.
137 // For each codim 1 intersection check if it is part of the boundary and if so,
138 // get its corner vertices, find resp. store them together with their associated index,
139 // and remember the indices of the boundary faces' corners.
140 {
141 // several counter for consecutive indexing are needed
142 int simplex_index = 0;
143 int vertex_index = 0;
144 IndexType eindex = 0; // supress warning
145
146 // needed later for insertion into a std::set which only
147 // works with const references
148
149 // a temporary container where newly acquired face
150 // information can be stored at first
151 std::deque<SubEntityInfo> temp_faces;
152
153 // iterate over interior codim 0 elements on the grid
154 for (ElementIter elit = this->gv_.template begin<0, PType>();
155 elit != this->gv_.template end<0, PType>(); ++elit)
156 {
157 const auto& elmt = *elit;
158 Dune::GeometryType gt = elmt.type();
159
160 // if some face is part of the surface add it!
161 if (elit->hasBoundaryIntersections())
162 {
163 // add an entry to the element info map, the index will be set properly later,
164 // whereas the number of faces is already known
165 eindex = this->cellMapper_.index(elmt);
166 this->elmtInfo_[eindex] = new ElementInfo(simplex_index, elmt, 0);
167
168 // now add the faces in ascending order of their indices
169 // (we are only talking about 1-4 faces here, so O(n^2) is ok!)
170 for (IsIter is = this->gv_.ibegin(elmt); is != this->gv_.iend(elmt); ++is)
171 {
172 // Stop only at selected boundary faces
173 if (!is->boundary() or !predicate(elmt, is->indexInInside()))
174 continue;
175
176 const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
177 // get the corner count of this face
178 const int face_corners = refElement.size(is->indexInInside(), 1, dim);
179
180 // now we only have to care about the 3D case, i.e. a triangle face can be
181 // inserted directly whereas a quadrilateral face has to be divided into two triangles
182 switch (face_corners)
183 {
184 case 2 :
185 case 3:
186 {
187 // we have a simplex here
188
189 // register the additional face(s)
190 this->elmtInfo_[eindex]->faces++;
191
192 // add a new face to the temporary collection
193 temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
194 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
195
196 std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
197
198 // try for each of the faces vertices whether it is already inserted or not
199 for (int i = 0; i < face_corners; ++i)
200 {
201 // get the number of the vertex in the parent element
202 int vertex_number = refElement.subEntity(is->indexInInside(), 1, i, dim);
203
204 // get the vertex pointer and the index from the index set
205 const Vertex vertex = elit->template subEntity<dim>(vertex_number);
206 cornerCoords[i] = vertex.geometry().corner(0);
207
208 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
209
210 // remember the vertex' number in parent element's vertices
211 temp_faces.back().corners[i].num = vertex_number;
212
213 // if the vertex is not yet inserted in the vertex info map
214 // it is a new one -> it will be inserted now!
215 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
216 if (vimit == this->vtxInfo_.end())
217 {
218 // insert into the map
219 this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
220 // remember the vertex as a corner of the current face in temp_faces
221 temp_faces.back().corners[i].idx = vertex_index;
222 // increase the current index
223 vertex_index++;
224 }
225 else
226 {
227 // only insert the index into the simplices array
228 temp_faces.back().corners[i].idx = vimit->second->idx;
229 }
230 }
231
232 // Now we have the correct vertices in the last entries of temp_faces, but they may
233 // have the wrong orientation. We want them to be oriented such that all boundary edges
234 // point in the counterclockwise direction. Therefore, we check the orientation of the
235 // new face and possibly switch the two vertices.
236 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
237
238 // Compute segment normal
239 FieldVector<ctype,dimworld> reconstructedNormal;
240 if (dim==2) // boundary face is a line segment
241 {
242 reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
243 reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
244 } else { // boundary face is a triangle
245 FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
246 FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
247 reconstructedNormal = crossProduct(segment1, segment2);
248 }
249 reconstructedNormal /= reconstructedNormal.two_norm();
250
251 if (realNormal * reconstructedNormal < 0.0)
252 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
253
254 // now increase the current face index
255 simplex_index++;
256 break;
257 }
258 case 4 :
259 {
260 assert(dim == 3);
261 // we have a quadrilateral here
262 unsigned int vertex_indices[4];
263 unsigned int vertex_numbers[4];
264
265 // register the additional face(s) (2 simplices)
266 this->elmtInfo_[eindex]->faces += 2;
267
268 std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
269
270 // get the vertex pointers for the quadrilateral's corner vertices
271 // and try for each of them whether it is already inserted or not
272 for (int i = 0; i < cube_corners; ++i)
273 {
274 // get the number of the vertex in the parent element
275 vertex_numbers[i] = refElement.subEntity(is->indexInInside(), 1, i, dim);
276
277 // get the vertex pointer and the index from the index set
278 const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
279 cornerCoords[i] = vertex.geometry().corner(0);
280
281 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
282
283 // if the vertex is not yet inserted in the vertex info map
284 // it is a new one -> it will be inserted now!
285 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
286 if (vimit == this->vtxInfo_.end())
287 {
288 // insert into the map
289 this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
290 // remember this vertex' index
291 vertex_indices[i] = vertex_index;
292 // increase the current index
293 vertex_index++;
294 }
295 else
296 {
297 // only remember the vertex' index
298 vertex_indices[i] = vimit->second->idx;
299 }
300 }
301
302 // now introduce the two triangles subdividing the quadrilateral
303 // ATTENTION: the order of vertices given by "orientedSubface" corresponds to the order
304 // of a Dune quadrilateral, i.e. the triangles are given by 0 1 2 and 3 2 1
305
306 // add a new face to the temporary collection for the first tri
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];
312 // remember the vertices' numbers in parent element's vertices
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];
316
317 // Now we have the correct vertices in the last entries of temp_faces, but they may
318 // have the wrong orientation. We want the triangle vertices on counterclockwise order,
319 // when viewed from the outside of the grid. Therefore, we check the orientation of the
320 // new face and possibly switch two vertices.
321 FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
322
323 // Compute segment normal
324 FieldVector<ctype,dimworld> reconstructedNormal = crossProduct(cornerCoords[1] - cornerCoords[0],
325 cornerCoords[2] - cornerCoords[0]);
326 reconstructedNormal /= reconstructedNormal.two_norm();
327
328 if (realNormal * reconstructedNormal < 0.0)
329 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
330
331
332 // add a new face to the temporary collection for the second tri
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];
338 // remember the vertices' numbers in parent element's vertices
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];
342
343 // Now we have the correct vertices in the last entries of temp_faces, but they may
344 // have the wrong orientation. We want the triangle vertices on counterclockwise order,
345 // when viewed from the outside of the grid. Therefore, we check the orientation of the
346 // new face and possibly switch two vertices.
347 // Compute segment normal
348 reconstructedNormal = crossProduct(cornerCoords[2] - cornerCoords[3],
349 cornerCoords[1] - cornerCoords[3]);
350 reconstructedNormal /= reconstructedNormal.two_norm();
351
352 if (realNormal * reconstructedNormal < 0.0)
353 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
354
355 simplex_index+=2;
356 break;
357 }
358 default :
359 DUNE_THROW(Dune::NotImplemented, "the extractor does only work for triangle and quadrilateral faces (" << face_corners << " corners)");
360 break;
361 }
362 } // end loop over found surface parts
363 }
364 } // end loop over elements
365
366 std::cout << "added " << simplex_index << " subfaces\n";
367
368 // allocate the array for the face specific information...
369 this->subEntities_.resize(simplex_index);
370 // ...and fill in the data from the temporary containers
371 copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
372 }
373
374
375 // now first write the array with the coordinates...
376 this->coords_.resize(this->vtxInfo_.size());
377 typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
378 for (; it1 != this->vtxInfo_.end(); ++it1)
379 {
380 // get a pointer to the associated info object
381 CoordinateInfo* current = &this->coords_[it1->second->idx];
382 // store this coordinates index // NEEDED?
383 current->index = it1->second->idx;
384 // store the vertex' index for the index2vertex mapping
385 current->vtxindex = it1->first;
386 // store the vertex' coordinates under the associated index
387 // in the coordinates array
388 const auto vtx = this->grid().entity(it1->second->p);
389 current->coord = vtx.geometry().corner(0);
390 }
391
392}
393
394} // namespace GridGlue
395
396} // namespace Dune
397
398#endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
extractor base class
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
int IndexType
Definition extractor.hh:74
@ dimworld
Definition extractor.hh:47
std::map< IndexType, VertexInfo * > VertexInfoMap
Definition extractor.hh:184
@ codim
Definition extractor.hh:49
@ cube_corners
Definition extractor.hh:53
simple struct holding a vertex pointer and an index
Definition extractor.hh:117
simple struct holding an element seed and an index
Definition extractor.hh:129
Holds some information about an element's subentity involved in a coupling.
Definition extractor.hh:148
Base class for subentity-selecting predicates.
Definition extractorpredicate.hh:31