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