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