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