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 */
13 /**
14 * @file
15 * @brief Grid extractor class for codim 1 subgrids
16 */
17
18 #ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
19 #define DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
20
21 #include "extractor.hh"
22
23 #include <array>
24 #include <deque>
25 #include <functional>
26
27 #include <dune/common/deprecated.hh>
28 #include <dune/common/version.hh>
29 #include <dune/grid-glue/common/crossproduct.hh>
30
31 namespace Dune {
32
33 namespace GridGlue {
34
35 /**
36 * extractor for codim-1 entities (facets)
37 */
38 template<typename GV>
39 class Codim1Extractor : public Extractor<GV,1>
40 {
41 public:
42
43 /* 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 */
44
45 using Extractor<GV,1>::dimworld;
46 using Extractor<GV,1>::dim;
47 using Extractor<GV,1>::codim;
48 using Extractor<GV,1>::cube_corners;
49 typedef typename Extractor<GV,1>::IndexType IndexType;
50
51 /// @brief compile time number of corners of surface simplices
52 static constexpr int simplex_corners = dim;
53
54 typedef GV GridView;
55
56 typedef typename GV::Grid::ctype ctype;
57 typedef Dune::FieldVector<ctype, dimworld> Coords;
58
59 typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
60 typedef typename GV::Traits::template Codim<0>::Entity Element;
61 typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
62
63 // import typedefs from base class
64 typedef typename Extractor<GV,1>::SubEntityInfo SubEntityInfo;
65 typedef typename Extractor<GV,1>::ElementInfo ElementInfo;
66 typedef typename Extractor<GV,1>::VertexInfo VertexInfo;
67 typedef typename Extractor<GV,1>::CoordinateInfo CoordinateInfo;
68 typedef typename Extractor<GV,1>::VertexInfoMap VertexInfoMap;
69
70 public:
71
72 /* 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 */
73
74 /**
75 * @brief Constructor
76 * @param gv the grid view object to work with
77 * @param predicate a predicate to mark entities for extraction (unary functor returning bool)
78 */
Codim1Extractor(const GV & gv,const Predicate & predicate)79 Codim1Extractor(const GV& gv, const Predicate& predicate)
80 : Extractor<GV,1>(gv)
81 {
82 std::cout << "This is Codim1Extractor on a <" << dim
83 << "," << dimworld << "> grid!"
84 << std::endl;
85 update(predicate);
86 }
87
88 private:
89
90 /**
91 * Extracts a codimension 1 surface from the grid @c g and builds up two arrays
92 * with the topology of the surface written to them. The description of the
93 * surface part that is to be extracted is given in form of a predicate class.
94 *
95 * Assumed that we are in 2D the coords array will have the structure
96 * x0 y0 x1 y1 ... x(n-1) y(n-1)
97 * Values in the @c _indices array then refer to the indices of the coordinates, e.g.
98 * index 1 is associated with the position x1. If the surface consists of triangles
99 * we have always groups of 3 indices describing one triangle.
100 *
101 * @param predicate a predicate that "selects" the faces to add to the surface
102 */
103 void update(const Predicate& predicate);
104
105 };
106
107
108 template<typename GV>
update(const Predicate & predicate)109 void Codim1Extractor<GV>::update(const Predicate& predicate)
110 {
111 // free everything there is in this object
112 this->clear();
113
114 // In this first pass iterate over all entities of codim 0.
115 // For each codim 1 intersection check if it is part of the boundary and if so,
116 // get its corner vertices, find resp. store them together with their associated index,
117 // and remember the indices of the boundary faces' corners.
118 {
119 // several counter for consecutive indexing are needed
120 int simplex_index = 0;
121 int vertex_index = 0;
122 IndexType eindex = 0; // supress warning
123
124 // needed later for insertion into a std::set which only
125 // works with const references
126
127 // a temporary container where newly acquired face
128 // information can be stored at first
129 std::deque<SubEntityInfo> temp_faces;
130
131 // iterate over interior codim 0 elements on the grid
132 for (const auto& elmt : elements(this->gv_, Partitions::interior))
133 {
134 Dune::GeometryType gt = elmt.type();
135
136 // if some face is part of the surface add it!
137 if (elmt.hasBoundaryIntersections())
138 {
139 // add an entry to the element info map, the index will be set properly later,
140 // whereas the number of faces is already known
141 eindex = this->cellMapper_.index(elmt);
142 this->elmtInfo_.emplace(eindex, ElementInfo(simplex_index, elmt, 0));
143
144 // now add the faces in ascending order of their indices
145 // (we are only talking about 1-4 faces here, so O(n^2) is ok!)
146 for (const auto& in : intersections(this->gv_, elmt))
147 {
148 // Stop only at selected boundary faces
149 if (!in.boundary() or !predicate(elmt, in.indexInInside()))
150 continue;
151
152 const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
153 // get the corner count of this face
154 const int face_corners = refElement.size(in.indexInInside(), 1, dim);
155
156 // now we only have to care about the 3D case, i.e. a triangle face can be
157 // inserted directly whereas a quadrilateral face has to be divided into two triangles
158 switch (face_corners)
159 {
160 case 2 :
161 case 3:
162 {
163 // we have a simplex here
164
165 // register the additional face(s)
166 this->elmtInfo_.at(eindex).faces++;
167
168 // add a new face to the temporary collection
169 temp_faces.emplace_back(eindex, in.indexInInside(),
170 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
171 Dune::GeometryTypes::simplex(dim-codim)
172 #else
173 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
174 #endif
175 );
176
177 std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
178
179 // try for each of the faces vertices whether it is already inserted or not
180 for (int i = 0; i < face_corners; ++i)
181 {
182 // get the number of the vertex in the parent element
183 int vertex_number = refElement.subEntity(in.indexInInside(), 1, i, dim);
184
185 // get the vertex pointer and the index from the index set
186 const Vertex vertex = elmt.template subEntity<dim>(vertex_number);
187 cornerCoords[i] = vertex.geometry().corner(0);
188
189 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
190
191 // remember the vertex' number in parent element's vertices
192 temp_faces.back().corners[i].num = vertex_number;
193
194 // if the vertex is not yet inserted in the vertex info map
195 // it is a new one -> it will be inserted now!
196 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
197 if (vimit == this->vtxInfo_.end())
198 {
199 // insert into the map
200 this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
201 // remember the vertex as a corner of the current face in temp_faces
202 temp_faces.back().corners[i].idx = vertex_index;
203 // increase the current index
204 vertex_index++;
205 }
206 else
207 {
208 // only insert the index into the simplices array
209 temp_faces.back().corners[i].idx = vimit->second.idx;
210 }
211 }
212
213 // Now we have the correct vertices in the last entries of temp_faces, but they may
214 // have the wrong orientation. We want them to be oriented such that all boundary edges
215 // point in the counterclockwise direction. Therefore, we check the orientation of the
216 // new face and possibly switch the two vertices.
217 FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
218
219 // Compute segment normal
220 FieldVector<ctype,dimworld> reconstructedNormal;
221 if (dim==2) // boundary face is a line segment
222 {
223 reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
224 reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
225 } else { // boundary face is a triangle
226 FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
227 FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
228 reconstructedNormal = crossProduct(segment1, segment2);
229 }
230 reconstructedNormal /= reconstructedNormal.two_norm();
231
232 if (realNormal * reconstructedNormal < 0.0)
233 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
234
235 // now increase the current face index
236 simplex_index++;
237 break;
238 }
239 case 4 :
240 {
241 assert(dim == 3 && cube_corners == 4);
242 // we have a quadrilateral here
243 std::array<unsigned int, 4> vertex_indices;
244 std::array<unsigned int, 4> vertex_numbers;
245
246 // register the additional face(s) (2 simplices)
247 this->elmtInfo_.at(eindex).faces += 2;
248
249 std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
250
251 // get the vertex pointers for the quadrilateral's corner vertices
252 // and try for each of them whether it is already inserted or not
253 for (int i = 0; i < cube_corners; ++i)
254 {
255 // get the number of the vertex in the parent element
256 vertex_numbers[i] = refElement.subEntity(in.indexInInside(), 1, i, dim);
257
258 // get the vertex pointer and the index from the index set
259 const Vertex vertex = elmt.template subEntity<dim>(vertex_numbers[i]);
260 cornerCoords[i] = vertex.geometry().corner(0);
261
262 IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
263
264 // if the vertex is not yet inserted in the vertex info map
265 // it is a new one -> it will be inserted now!
266 typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
267 if (vimit == this->vtxInfo_.end())
268 {
269 // insert into the map
270 this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
271 // remember this vertex' index
272 vertex_indices[i] = vertex_index;
273 // increase the current index
274 vertex_index++;
275 }
276 else
277 {
278 // only remember the vertex' index
279 vertex_indices[i] = vimit->second.idx;
280 }
281 }
282
283 // now introduce the two triangles subdividing the quadrilateral
284 // ATTENTION: the order of vertices given by "orientedSubface" corresponds to the order
285 // of a Dune quadrilateral, i.e. the triangles are given by 0 1 2 and 3 2 1
286
287 // add a new face to the temporary collection for the first tri
288 temp_faces.emplace_back(eindex, in.indexInInside(),
289 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
290 Dune::GeometryTypes::simplex(dim-codim)
291 #else
292 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
293 #endif
294 );
295 temp_faces.back().corners[0].idx = vertex_indices[0];
296 temp_faces.back().corners[1].idx = vertex_indices[1];
297 temp_faces.back().corners[2].idx = vertex_indices[2];
298 // remember the vertices' numbers in parent element's vertices
299 temp_faces.back().corners[0].num = vertex_numbers[0];
300 temp_faces.back().corners[1].num = vertex_numbers[1];
301 temp_faces.back().corners[2].num = vertex_numbers[2];
302
303 // Now we have the correct vertices in the last entries of temp_faces, but they may
304 // have the wrong orientation. We want the triangle vertices on counterclockwise order,
305 // when viewed from the outside of the grid. Therefore, we check the orientation of the
306 // new face and possibly switch two vertices.
307 FieldVector<ctype,dimworld> realNormal = in.centerUnitOuterNormal();
308
309 // Compute segment normal
310 FieldVector<ctype,dimworld> reconstructedNormal = crossProduct(cornerCoords[1] - cornerCoords[0],
311 cornerCoords[2] - cornerCoords[0]);
312 reconstructedNormal /= reconstructedNormal.two_norm();
313
314 if (realNormal * reconstructedNormal < 0.0)
315 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
316
317
318 // add a new face to the temporary collection for the second tri
319 temp_faces.emplace_back(eindex, in.indexInInside(),
320 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
321 Dune::GeometryTypes::simplex(dim-codim)
322 #else
323 Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)
324 #endif
325 );
326 temp_faces.back().corners[0].idx = vertex_indices[3];
327 temp_faces.back().corners[1].idx = vertex_indices[2];
328 temp_faces.back().corners[2].idx = vertex_indices[1];
329 // remember the vertices' numbers in parent element's vertices
330 temp_faces.back().corners[0].num = vertex_numbers[3];
331 temp_faces.back().corners[1].num = vertex_numbers[2];
332 temp_faces.back().corners[2].num = vertex_numbers[1];
333
334 // Now we have the correct vertices in the last entries of temp_faces, but they may
335 // have the wrong orientation. We want the triangle vertices on counterclockwise order,
336 // when viewed from the outside of the grid. Therefore, we check the orientation of the
337 // new face and possibly switch two vertices.
338 // Compute segment normal
339 reconstructedNormal = crossProduct(cornerCoords[2] - cornerCoords[3],
340 cornerCoords[1] - cornerCoords[3]);
341 reconstructedNormal /= reconstructedNormal.two_norm();
342
343 if (realNormal * reconstructedNormal < 0.0)
344 std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
345
346 simplex_index+=2;
347 break;
348 }
349 default :
350 DUNE_THROW(Dune::NotImplemented, "the extractor does only work for triangle and quadrilateral faces (" << face_corners << " corners)");
351 break;
352 }
353 } // end loop over found surface parts
354 }
355 } // end loop over elements
356
357 std::cout << "added " << simplex_index << " subfaces\n";
358
359 // allocate the array for the face specific information...
360 this->subEntities_.resize(simplex_index);
361 // ...and fill in the data from the temporary containers
362 copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
363 }
364
365
366 // now first write the array with the coordinates...
367 this->coords_.resize(this->vtxInfo_.size());
368 for (const auto& vinfo : this->vtxInfo_)
369 {
370 // get a pointer to the associated info object
371 CoordinateInfo* current = &this->coords_[vinfo.second.idx];
372 // store this coordinates index // NEEDED?
373 current->index = vinfo.second.idx;
374 // store the vertex' index for the index2vertex mapping
375 current->vtxindex = vinfo.first;
376 // store the vertex' coordinates under the associated index
377 // in the coordinates array
378 const auto vtx = this->grid().entity(vinfo.second.p);
379 current->coord = vtx.geometry().corner(0);
380 }
381
382 }
383
384 } // namespace GridGlue
385
386 } // namespace Dune
387
388 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
389