dune-grid-glue  2.3.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 
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<0>::EntityPointer ElementPtr;
58 
59  static const Dune::PartitionIteratorType PType = Dune::Interior_Partition;
60  typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator ElementIter;
61 
62  typedef typename GV::IntersectionIterator IsIter;
63 
64  // import typedefs from base class
70 
71 public:
72 
73  /* 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 */
74 
80  Codim1Extractor(const GV& gv, const ExtractorPredicate<GV,1>& descr)
81  : Extractor<GV,1>(gv)
82  {
83  std::cout << "This is Codim1Extractor on a <" << dim
84  << "," << dimworld << "> grid!"
85  << std::endl;
86  update(descr);
87  }
88 
89 private:
90 
104  void update(const ExtractorPredicate<GV,1>& descr);
105 
106 };
107 
108 
109 template<typename GV>
110 void Codim1Extractor<GV>::update(const ExtractorPredicate<GV,1>& descr)
111 {
112  // free everything there is in this object
113  this->clear();
114 
115  // In this first pass iterate over all entities of codim 0.
116  // For each codim 1 intersection check if it is part of the boundary and if so,
117  // get its corner vertices, find resp. store them together with their associated index,
118  // and remember the indices of the boundary faces' corners.
119  {
120  // several counter for consecutive indexing are needed
121  int simplex_index = 0;
122  int vertex_index = 0;
123  IndexType eindex = 0; // supress warning
124 
125  // needed later for insertion into a std::set which only
126  // works with const references
127 
128  // a temporary container where newly acquired face
129  // information can be stored at first
130  std::deque<SubEntityInfo> temp_faces;
131 
132  // iterate over interior codim 0 elements on the grid
133  for (ElementIter elit = this->gv_.template begin<0, PType>();
134  elit != this->gv_.template end<0, PType>(); ++elit)
135  {
136  ElementPtr eptr(elit);
137  Dune::GeometryType gt = elit->type();
138 
139  // if some face is part of the surface add it!
140  if (elit->hasBoundaryIntersections())
141  {
142  // add an entry to the element info map, the index will be set properly later,
143  // whereas the number of faces is already known
144  eindex = this->cellMapper_.map(*elit);
145  this->elmtInfo_[eindex] = new ElementInfo(simplex_index, eptr, 0);
146 
147  // now add the faces in ascending order of their indices
148  // (we are only talking about 1-4 faces here, so O(n^2) is ok!)
149  for (IsIter is = this->gv_.ibegin(*elit); is != this->gv_.iend(*elit); ++is)
150  {
151  // Stop only at selected boundary faces
152  if (!is->boundary() or !descr.contains(eptr, is->indexInInside()))
153  continue;
154 
155 #if DUNE_VERSION_NEWER(DUNE_GEOMETRY,2,3)
156  const Dune::ReferenceElement<ctype, dim>& refElement = Dune::ReferenceElements<ctype, dim>::general(gt);
157 #else
158  const Dune::GenericReferenceElement<ctype, dim>& refElement = Dune::GenericReferenceElements<ctype, dim>::general(gt);
159 #endif
160  // get the corner count of this face
161  const int face_corners = refElement.size(is->indexInInside(), 1, dim);
162 
163  // now we only have to care about the 3D case, i.e. a triangle face can be
164  // inserted directly whereas a quadrilateral face has to be divided into two triangles
165  switch (face_corners)
166  {
167  case 2 :
168  case 3:
169  {
170  // we have a simplex here
171 
172  // register the additional face(s)
173  this->elmtInfo_[eindex]->faces++;
174 
175  // add a new face to the temporary collection
176  temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
177  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
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(is->indexInInside(), 1, i, dim);
186 
187  // get the vertex pointer and the index from the index set
188  VertexPtr vptr(elit->template subEntity<dim>(vertex_number));
189  cornerCoords[i] = vptr->geometry().corner(0);
190 
191  IndexType vindex = this->gv_.indexSet().template index<dim>(*vptr);
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_[vindex] = new VertexInfo(vertex_index, vptr);
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 = is->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 = Projection::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);
244  // we have a quadrilateral here
245  unsigned int vertex_indices[4];
246  unsigned int vertex_numbers[4];
247 
248  // register the additional face(s) (2 simplices)
249  this->elmtInfo_[eindex]->faces += 2;
250 
251  Dune::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(is->indexInInside(), 1, i, dim);
259 
260  // get the vertex pointer and the index from the index set
261  VertexPtr vptr(elit->template subEntity<dim>(vertex_numbers[i]));
262  cornerCoords[i] = vptr->geometry().corner(0);
263 
264  IndexType vindex = this->gv_.indexSet().template index<dim>(*vptr);
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_[vindex] = new VertexInfo(vertex_index, vptr);
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.push_back(SubEntityInfo(eindex, is->indexInInside(),
291  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
292  temp_faces.back().corners[0].idx = vertex_indices[0];
293  temp_faces.back().corners[1].idx = vertex_indices[1];
294  temp_faces.back().corners[2].idx = vertex_indices[2];
295  // remember the vertices' numbers in parent element's vertices
296  temp_faces.back().corners[0].num = vertex_numbers[0];
297  temp_faces.back().corners[1].num = vertex_numbers[1];
298  temp_faces.back().corners[2].num = vertex_numbers[2];
299 
300  // Now we have the correct vertices in the last entries of temp_faces, but they may
301  // have the wrong orientation. We want the triangle vertices on counterclockwise order,
302  // when viewed from the outside of the grid. Therefore, we check the orientation of the
303  // new face and possibly switch two vertices.
304  FieldVector<ctype,dimworld> realNormal = is->centerUnitOuterNormal();
305 
306  // Compute segment normal
307  FieldVector<ctype,dimworld> reconstructedNormal = Projection::crossProduct(cornerCoords[1] - cornerCoords[0],
308  cornerCoords[2] - cornerCoords[0]);
309  reconstructedNormal /= reconstructedNormal.two_norm();
310 
311  if (realNormal * reconstructedNormal < 0.0)
312  std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
313 
314 
315  // add a new face to the temporary collection for the second tri
316  temp_faces.push_back(SubEntityInfo(eindex, is->indexInInside(),
317  Dune::GeometryType(Dune::GeometryType::simplex,dim-codim)));
318  temp_faces.back().corners[0].idx = vertex_indices[3];
319  temp_faces.back().corners[1].idx = vertex_indices[2];
320  temp_faces.back().corners[2].idx = vertex_indices[1];
321  // remember the vertices' numbers in parent element's vertices
322  temp_faces.back().corners[0].num = vertex_numbers[3];
323  temp_faces.back().corners[1].num = vertex_numbers[2];
324  temp_faces.back().corners[2].num = vertex_numbers[1];
325 
326  // Now we have the correct vertices in the last entries of temp_faces, but they may
327  // have the wrong orientation. We want the triangle vertices on counterclockwise order,
328  // when viewed from the outside of the grid. Therefore, we check the orientation of the
329  // new face and possibly switch two vertices.
330  // Compute segment normal
331  reconstructedNormal = Projection::crossProduct(cornerCoords[2] - cornerCoords[3],
332  cornerCoords[1] - cornerCoords[3]);
333  reconstructedNormal /= reconstructedNormal.two_norm();
334 
335  if (realNormal * reconstructedNormal < 0.0)
336  std::swap(temp_faces.back().corners[0], temp_faces.back().corners[1]);
337 
338  simplex_index+=2;
339  break;
340  }
341  default :
342  DUNE_THROW(Dune::NotImplemented, "the extractor does only work for triangle and quadrilateral faces (" << face_corners << " corners)");
343  break;
344  }
345  } // end loop over found surface parts
346  }
347  } // end loop over elements
348 
349  std::cout << "added " << simplex_index << " subfaces\n";
350 
351  // allocate the array for the face specific information...
352  this->subEntities_.resize(simplex_index);
353  // ...and fill in the data from the temporary containers
354  copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
355  }
356 
357 
358  // now first write the array with the coordinates...
359  this->coords_.resize(this->vtxInfo_.size());
360  typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
361  for (; it1 != this->vtxInfo_.end(); ++it1)
362  {
363  // get a pointer to the associated info object
364  CoordinateInfo* current = &this->coords_[it1->second->idx];
365  // store this coordinates index // NEEDED?
366  current->index = it1->second->idx;
367  // store the vertex' index for the index2vertex mapping
368  current->vtxindex = it1->first;
369  // store the vertex' coordinates under the associated index
370  // in the coordinates array
371  current->coord = it1->second->p->geometry().corner(0);
372  }
373 
374 }
375 
376 } // namespace GridGlue
377 
378 } // namespace Dune
379 
380 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM1EXTRACTOR_HH
Extractor< GV, 1 >::VertexInfoMap VertexInfoMap
Definition: codim1extractor.hh:69
GV::IntersectionIterator IsIter
Definition: codim1extractor.hh:62
Definition: gridglue.hh:34
Extractor< GV, 1 >::ElementInfo ElementInfo
Definition: codim1extractor.hh:66
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:54
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:67
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim1extractor.hh:57
Codim1Extractor(const GV &gv, const ExtractorPredicate< GV, 1 > &descr)
Constructor.
Definition: codim1extractor.hh:80
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:19
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim1extractor.hh:60
extractor base class
static const Dune::PartitionIteratorType PType
Definition: codim1extractor.hh:59
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:48
Definition: codim1extractor.hh:33
Extractor< GV, 1 >::SubEntityInfo SubEntityInfo
Definition: codim1extractor.hh:65
Extractor< GV, 1 >::IndexType IndexType
Definition: codim1extractor.hh:43
Base class for predicates selecting the part of a grid to be extracted.
Extractor< GV, 1 >::CoordinateInfo CoordinateInfo
Definition: codim1extractor.hh:68
GV GridView
Definition: codim1extractor.hh:51