dune-grid-glue  2.6-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 
23 #include <array>
24 #include <deque>
25 #include <functional>
26 
27 #include <dune/common/deprecated.hh>
28 #include <dune/common/version.hh>
30 
31 namespace Dune {
32 
33  namespace GridGlue {
34 
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 
50 
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
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 
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 
103  void update(const Predicate& predicate);
104 
105 };
106 
107 
108 template<typename GV>
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
static constexpr auto dim
Definition: extractor.hh:49
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim1extractor.hh:60
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
extractor base class
static constexpr int simplex_corners
compile time number of corners of surface simplices
Definition: codim1extractor.hh:52
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim1extractor.hh:61
static constexpr auto codim
Definition: extractor.hh:50
static constexpr int cube_corners
Definition: extractor.hh:52
Definition: codim1extractor.hh:39
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:65
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:64
Vertex vertex(unsigned int index) const
gets the vertex for a given coordinate index throws an exception if index not valid ...
Definition: extractor.hh:394
IteratorRange<... > intersections(const GridGlue<... > &glue, const Reverse<... > &reverse=!reversed)
Iterate over all intersections of a GridGlue.
const GridView gv_
the grid object to extract the surface from
Definition: extractor.hh:199
static constexpr auto dimworld
Definition: extractor.hh:48
const Grid & grid() const
Definition: extractor.hh:366
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:67
VertexInfoMap vtxInfo_
a map enabling faster access to vertices and coordinates
Definition: extractor.hh:214
GV::Grid::ctype ctype
Definition: codim1extractor.hh:56
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim1extractor.hh:79
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:66
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:57
ElementInfoMap elmtInfo_
a map enabling faster access to elements and faces
Definition: extractor.hh:221
GV GridView
Definition: codim1extractor.hh:54
std::vector< SubEntityInfo > subEntities_
all information about the extracted subEntities
Definition: extractor.hh:207
CellMapper cellMapper_
Definition: extractor.hh:223
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim1extractor.hh:59
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:68
Definition: gridglue.hh:35
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:49
std::vector< CoordinateInfo > coords_
all information about the corner vertices of the extracted
Definition: extractor.hh:204
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:43
void clear()
delete everything build up so far and free the memory
Definition: extractor.hh:247