dune-grid-glue  2.9
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 
33 namespace Dune {
34 
35  namespace GridGlue {
36 
40 template<typename GV>
41 class Codim1Extractor : public Extractor<GV,1>
42 {
43 public:
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 
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 
72 public:
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 
90 private:
91 
105  void update(const Predicate& predicate);
106 
107 };
108 
109 
110 template<typename GV>
111 void 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
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
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
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
static constexpr auto dim
Definition: extractor.hh:51