dune-grid-glue  2.4.0
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 
30 namespace Dune {
31 
32  namespace GridGlue {
33 
34 template<typename GV>
35 class Codim1Extractor : public Extractor<GV,1>
36 {
37 public:
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 
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>::EntityPointer VertexPtr;
59  typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
60  typedef typename GV::Traits::template Codim<0>::EntityPointer ElementPtr;
61  typedef typename GV::Traits::template Codim<0>::Entity Element;
62  typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
63 
64  static const Dune::PartitionIteratorType PType = Dune::Interior_Partition;
65  typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator ElementIter;
66 
67  typedef typename GV::IntersectionIterator IsIter;
68 
69  // import typedefs from base class
75 
76 public:
77 
78  /* 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 */
79 
85  DUNE_DEPRECATED_MSG("Please use a std::function<bool(const Element&, unsigned int)> in favor of the ExtractorPredicate.")
86  Codim1Extractor(const GV& gv, const ExtractorPredicate<GV,1>& descr)
87  : Extractor<GV,1>(gv)
88  {
89  std::cout << "This is Codim1Extractor on a <" << dim
90  << "," << dimworld << "> grid!"
91  << std::endl;
92  const auto predicate = [&](const Element& element, unsigned int subentity) -> bool {
93  return descr.contains(element, subentity);
94  };
95  update(predicate);
96  }
97 
103  Codim1Extractor(const GV& gv, const Predicate& predicate)
104  : Extractor<GV,1>(gv)
105  {
106  std::cout << "This is Codim1Extractor on a <" << dim
107  << "," << dimworld << "> grid!"
108  << std::endl;
109  update(predicate);
110  }
111 
112 private:
113 
127  void update(const Predicate& predicate);
128 
129 };
130 
131 
132 template<typename GV>
133 void Codim1Extractor<GV>::update(const Predicate& predicate)
134 {
135  // free everything there is in this object
136  this->clear();
137 
138  // In this first pass iterate over all entities of codim 0.
139  // For each codim 1 intersection check if it is part of the boundary and if so,
140  // get its corner vertices, find resp. store them together with their associated index,
141  // and remember the indices of the boundary faces' corners.
142  {
143  // several counter for consecutive indexing are needed
144  int simplex_index = 0;
145  int vertex_index = 0;
146  IndexType eindex = 0; // supress warning
147 
148  // needed later for insertion into a std::set which only
149  // works with const references
150 
151  // a temporary container where newly acquired face
152  // information can be stored at first
153  std::deque<SubEntityInfo> temp_faces;
154 
155  // iterate over interior codim 0 elements on the grid
156  for (ElementIter elit = this->gv_.template begin<0, PType>();
157  elit != this->gv_.template end<0, PType>(); ++elit)
158  {
159  const auto& elmt = *elit;
160  Dune::GeometryType gt = elmt.type();
161 
162  // if some face is part of the surface add it!
163  if (elit->hasBoundaryIntersections())
164  {
165  // add an entry to the element info map, the index will be set properly later,
166  // whereas the number of faces is already known
167 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
168  eindex = this->cellMapper_.index(elmt);
169 #else
170  eindex = this->cellMapper_.map(elmt);
171 #endif
172  this->elmtInfo_[eindex] = new ElementInfo(simplex_index, elmt, 0);
173 
174  // now add the faces in ascending order of their indices
175  // (we are only talking about 1-4 faces here, so O(n^2) is ok!)
176  for (IsIter is = this->gv_.ibegin(elmt); is != this->gv_.iend(elmt); ++is)
177  {
178  // Stop only at selected boundary faces
179  if (!is->boundary() or !predicate(elmt, is->indexInInside()))
180  continue;
181 
182 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
183  const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
184 #else
185  const Dune::GenericReferenceElement<ctype, dim>& refElement = Dune::GenericReferenceElements<ctype, dim>::general(gt);
186 #endif
187  // get the corner count of this face
188  const int face_corners = refElement.size(is->indexInInside(), 1, dim);
189 
190  // now we only have to care about the 3D case, i.e. a triangle face can be
191  // inserted directly whereas a quadrilateral face has to be divided into two triangles
192  switch (face_corners)
193  {
194  case 2 :
195  case 3:
196  {
197  // we have a simplex here
198 
199  // register the additional face(s)
200  this->elmtInfo_[eindex]->faces++;
201 
202  // add a new face to the temporary collection
203  temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
204  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
205 
206  std::vector<FieldVector<ctype,dimworld> > cornerCoords(face_corners);
207 
208  // try for each of the faces vertices whether it is already inserted or not
209  for (int i = 0; i < face_corners; ++i)
210  {
211  // get the number of the vertex in the parent element
212  int vertex_number = refElement.subEntity(is->indexInInside(), 1, i, dim);
213 
214  // get the vertex pointer and the index from the index set
215 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
216  const Vertex vertex = elit->template subEntity<dim>(vertex_number);
217 #else
218  const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_number);
219  const Vertex& vertex = *vertexPtr;
220 #endif
221  cornerCoords[i] = vertex.geometry().corner(0);
222 
223  IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
224 
225  // remember the vertex' number in parent element's vertices
226  temp_faces.back().corners[i].num = vertex_number;
227 
228  // if the vertex is not yet inserted in the vertex info map
229  // it is a new one -> it will be inserted now!
230  typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
231  if (vimit == this->vtxInfo_.end())
232  {
233  // insert into the map
234  this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
235  // remember the vertex as a corner of the current face in temp_faces
236  temp_faces.back().corners[i].idx = vertex_index;
237  // increase the current index
238  vertex_index++;
239  }
240  else
241  {
242  // only insert the index into the simplices array
243  temp_faces.back().corners[i].idx = vimit->second->idx;
244  }
245  }
246 
247  // Now we have the correct vertices in the last entries of temp_faces, but they may
248  // have the wrong orientation. We want them to be oriented such that all boundary edges
249  // point in the counterclockwise direction. Therefore, we check the orientation of the
250  // new face and possibly switch the two vertices.
251  FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
252 
253  // Compute segment normal
254  FieldVector<ctype,dimworld> reconstructedNormal;
255  if (dim==2) // boundary face is a line segment
256  {
257  reconstructedNormal[0] = cornerCoords[1][1] - cornerCoords[0][1];
258  reconstructedNormal[1] = cornerCoords[0][0] - cornerCoords[1][0];
259  } else { // boundary face is a triangle
260  FieldVector<ctype,dimworld> segment1 = cornerCoords[1] - cornerCoords[0];
261  FieldVector<ctype,dimworld> segment2 = cornerCoords[2] - cornerCoords[0];
262  reconstructedNormal = crossProduct(segment1, segment2);
263  }
264  reconstructedNormal /= reconstructedNormal.two_norm();
265 
266  if (realNormal * reconstructedNormal < 0.0)
267  std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
268 
269  // now increase the current face index
270  simplex_index++;
271  break;
272  }
273  case 4 :
274  {
275  assert(dim == 3);
276  // we have a quadrilateral here
277  unsigned int vertex_indices[4];
278  unsigned int vertex_numbers[4];
279 
280  // register the additional face(s) (2 simplices)
281  this->elmtInfo_[eindex]->faces += 2;
282 
283  std::array<FieldVector<ctype,dimworld>, 4> cornerCoords;
284 
285  // get the vertex pointers for the quadrilateral's corner vertices
286  // and try for each of them whether it is already inserted or not
287  for (int i = 0; i < cube_corners; ++i)
288  {
289  // get the number of the vertex in the parent element
290  vertex_numbers[i] = refElement.subEntity(is->indexInInside(), 1, i, dim);
291 
292  // get the vertex pointer and the index from the index set
293 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
294  const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
295 #else
296  const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_numbers[i]);
297  const Vertex &vertex = *vertexPtr;
298 #endif
299  cornerCoords[i] = vertex.geometry().corner(0);
300 
301  IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
302 
303  // if the vertex is not yet inserted in the vertex info map
304  // it is a new one -> it will be inserted now!
305  typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
306  if (vimit == this->vtxInfo_.end())
307  {
308  // insert into the map
309  this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
310  // remember this vertex' index
311  vertex_indices[i] = vertex_index;
312  // increase the current index
313  vertex_index++;
314  }
315  else
316  {
317  // only remember the vertex' index
318  vertex_indices[i] = vimit->second->idx;
319  }
320  }
321 
322  // now introduce the two triangles subdividing the quadrilateral
323  // ATTENTION: the order of vertices given by "orientedSubface" corresponds to the order
324  // of a Dune quadrilateral, i.e. the triangles are given by 0 1 2 and 3 2 1
325 
326  // add a new face to the temporary collection for the first tri
327  temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
328  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
329  temp_faces.back().corners[0].idx = vertex_indices[0];
330  temp_faces.back().corners[1].idx = vertex_indices[1];
331  temp_faces.back().corners[2].idx = vertex_indices[2];
332  // remember the vertices' numbers in parent element's vertices
333  temp_faces.back().corners[0].num = vertex_numbers[0];
334  temp_faces.back().corners[1].num = vertex_numbers[1];
335  temp_faces.back().corners[2].num = vertex_numbers[2];
336 
337  // Now we have the correct vertices in the last entries of temp_faces, but they may
338  // have the wrong orientation. We want the triangle vertices on counterclockwise order,
339  // when viewed from the outside of the grid. Therefore, we check the orientation of the
340  // new face and possibly switch two vertices.
341  FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
342 
343  // Compute segment normal
344  FieldVector<ctype,dimworld> reconstructedNormal = crossProduct(cornerCoords[1] - cornerCoords[0],
345  cornerCoords[2] - cornerCoords[0]);
346  reconstructedNormal /= reconstructedNormal.two_norm();
347 
348  if (realNormal * reconstructedNormal < 0.0)
349  std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
350 
351 
352  // add a new face to the temporary collection for the second tri
353  temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
354  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
355  temp_faces.back().corners[0].idx = vertex_indices[3];
356  temp_faces.back().corners[1].idx = vertex_indices[2];
357  temp_faces.back().corners[2].idx = vertex_indices[1];
358  // remember the vertices' numbers in parent element's vertices
359  temp_faces.back().corners[0].num = vertex_numbers[3];
360  temp_faces.back().corners[1].num = vertex_numbers[2];
361  temp_faces.back().corners[2].num = vertex_numbers[1];
362 
363  // Now we have the correct vertices in the last entries of temp_faces, but they may
364  // have the wrong orientation. We want the triangle vertices on counterclockwise order,
365  // when viewed from the outside of the grid. Therefore, we check the orientation of the
366  // new face and possibly switch two vertices.
367  // Compute segment normal
368  reconstructedNormal = crossProduct(cornerCoords[2] - cornerCoords[3],
369  cornerCoords[1] - cornerCoords[3]);
370  reconstructedNormal /= reconstructedNormal.two_norm();
371 
372  if (realNormal * reconstructedNormal < 0.0)
373  std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
374 
375  simplex_index+=2;
376  break;
377  }
378  default :
379  DUNE_THROW(Dune::NotImplemented, "the extractor does only work for triangle and quadrilateral faces (" << face_corners << " corners)");
380  break;
381  }
382  } // end loop over found surface parts
383  }
384  } // end loop over elements
385 
386  std::cout << "added " << simplex_index << " subfaces\n";
387 
388  // allocate the array for the face specific information...
389  this->subEntities_.resize(simplex_index);
390  // ...and fill in the data from the temporary containers
391  copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
392  }
393 
394 
395  // now first write the array with the coordinates...
396  this->coords_.resize(this->vtxInfo_.size());
397  typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
398  for (; it1 != this->vtxInfo_.end(); ++it1)
399  {
400  // get a pointer to the associated info object
401  CoordinateInfo* current = &this->coords_[it1->second->idx];
402  // store this coordinates index // NEEDED?
403  current->index = it1->second->idx;
404  // store the vertex' index for the index2vertex mapping
405  current->vtxindex = it1->first;
406  // store the vertex' coordinates under the associated index
407  // in the coordinates array
408 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
409  const auto vtx = this->grid().entity(it1->second->p);
410 #else
411  const auto vtxPtr = this->grid().entityPointer(it1->second->p);
412  const auto& vtx = *vtxPtr;
413 #endif
414  current->coord = vtx.geometry().corner(0);
415  }
416 
417 }
418 
419 } // namespace GridGlue
420 
421 } // namespace Dune
422 
423 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
static const Dune::PartitionIteratorType PType
Definition: codim1extractor.hh:64
GV::Grid::ctype ctype
Definition: codim1extractor.hh:55
const GridView gv_
the grid object to extract the surface from
Definition: extractor.hh:206
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim1extractor.hh:62
GV::IntersectionIterator IsIter
Definition: codim1extractor.hh:67
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:71
void clear()
delete everything build up so far and free the memory
Definition: extractor.hh:252
std::vector< CoordinateInfo > coords_
all information about the corner vertices of the extracted
Definition: extractor.hh:211
std::vector< SubEntityInfo > subEntities_
all information about the extracted subEntities
Definition: extractor.hh:214
Definition: gridglue.hh:33
Extractor< GV, 1 >::VertexInfo VertexInfo
Definition: codim1extractor.hh:72
extractor base class
Definition: codim1extractor.hh:50
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:396
Dune::FieldVector< ctype, dimworld > Coords
Definition: codim1extractor.hh:56
Codim1Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim1extractor.hh:103
const Grid & grid() const
Definition: extractor.hh:380
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:70
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:74
Definition: extractor.hh:53
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim1extractor.hh:59
GV::Traits::template Codim< dim >::EntityPointer VertexPtr
Definition: codim1extractor.hh:58
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:45
GV GridView
Definition: codim1extractor.hh:53
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:30
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:73
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim1extractor.hh:65
VertexInfoMap vtxInfo_
a map enabling faster access to vertices and coordinates
Definition: extractor.hh:221
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim1extractor.hh:61
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim1extractor.hh:60
Base class for predicates selecting the part of a grid to be extracted.
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
CellMapper cellMapper_
Definition: extractor.hh:230
Vertex vertex(unsigned int index) const
gets the vertex for a given coordinate index throws an exception if index not valid ...
Definition: extractor.hh:420
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:47
Definition: codim1extractor.hh:35
ElementInfoMap elmtInfo_
a map enabling faster access to elements and faces
Definition: extractor.hh:228