dune-grid-glue  2.4.0
codim0extractor.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: codim0extractor.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: base class for grid extractors extracting surface grids
11  *
12  */
18 #ifndef DUNE_GRIDGLUE_EXTRACTORS_CODIM0EXTRACTOR_HH
19 #define DUNE_GRIDGLUE_EXTRACTORS_CODIM0EXTRACTOR_HH
20 
21 #include <deque>
22 #include <functional>
23 
24 #include <dune/common/deprecated.hh>
25 #include <dune/grid/common/mcmgmapper.hh>
26 
27 #include "extractor.hh"
28 #include "extractorpredicate.hh"
29 
30 namespace Dune {
31 
32  namespace GridGlue {
33 
34 template<typename GV>
35 class Codim0Extractor : public Extractor<GV,0>
36 {
37 
38 public:
39 
40  /* 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 */
42  typedef typename Extractor<GV,0>::ctype ctype;
46 
47  typedef typename GV::Traits::template Codim<dim>::EntityPointer VertexPtr;
48  typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
49  typedef typename GV::Traits::template Codim<0>::EntityPointer ElementPtr;
50  typedef typename GV::Traits::template Codim<0>::Entity Element;
51  typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
52 
53  static const Dune::PartitionIteratorType PType = Dune::Interior_Partition;
54  typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator ElementIter;
55 
56  // import typedefs from base class
62 
68  DUNE_DEPRECATED_MSG("Please use a std::function<bool(const Element&, unsigned int)> in favor of the ExtractorPredicate.")
69  Codim0Extractor(const GV& gv, const ExtractorPredicate<GV,0>& descr)
70  : Extractor<GV,0>(gv), positiveNormalDirection_(false)
71  {
72  std::cout << "This is Codim0Extractor on a <"
73  << GV::dimension << "," << GV::dimensionworld << "> grid!" << std::endl;
74  const auto predicate = [&](const Element& element, unsigned int subentity) -> bool {
75  return descr.contains(element, subentity);
76  };
77  update(predicate);
78  }
79 
85  Codim0Extractor(const GV& gv, const Predicate& predicate)
86  : Extractor<GV,0>(gv), positiveNormalDirection_(false)
87  {
88  std::cout << "This is Codim0Extractor on a <"
89  << GV::dimension << "," << GV::dimensionworld << "> grid!" << std::endl;
90  update(predicate);
91  }
92 
94  const bool & positiveNormalDirection() const { return positiveNormalDirection_; }
95 
96 protected:
98 private:
99  void update(const Predicate& predicate);
100 };
101 
102 
103 template<typename GV>
104 void Codim0Extractor<GV>::update(const Predicate& predicate)
105 {
106  // In this first pass iterate over all entities of codim 0.
107  // Get its corner vertices, find resp. store them together with their associated index,
108  // and remember the indices of the corners.
109 
110  // free everything there is in this object
111  this->clear();
112 
113  // several counter for consecutive indexing are needed
114  size_t element_index = 0;
115  size_t vertex_index = 0;
116 
117  // a temporary container where newly acquired face
118  // information can be stored at first
119  std::deque<SubEntityInfo> temp_faces;
120 
121  // iterate over all codim 0 elements on the grid
122  for (ElementIter elit = this->gv_.template begin<0, PType>();
123  elit != this->gv_.template end<0, PType>(); ++elit)
124  {
125  const auto& elmt = *elit;
126  const auto geometry = elmt.geometry();
127 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
128  IndexType eindex = this->cellMapper_.index(elmt);
129 #else
130  IndexType eindex = this->cellMapper_.map(elmt);
131 #endif
132 
133  // only do sth. if this element is "interesting"
134  // implicit cast is done automatically
135  if (predicate(elmt,0))
136  {
137  // add an entry to the element info map, the index will be set properly later
138  this->elmtInfo_[eindex] = new ElementInfo(element_index, elmt, 1);
139 
140 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
141  unsigned int numCorners = elmt.subEntities(dim);
142 #else
143  unsigned int numCorners = elmt.template count<dim>();
144 #endif
145  unsigned int vertex_indices[numCorners]; // index in global vector
146  unsigned int vertex_numbers[numCorners]; // index in parent entity
147 
148  // try for each of the faces vertices whether it is already inserted or not
149  for (unsigned int i = 0; i < numCorners; ++i)
150  {
151  vertex_numbers[i] = i;
152 
153  // get the vertex pointer and the index from the index set
154 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
155  const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
156 #else
157  const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_numbers[i]);
158  const Vertex &vertex = *vertexPtr;
159 #endif
160  IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
161 
162  // if the vertex is not yet inserted in the vertex info map
163  // it is a new one -> it will be inserted now!
164  typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
165  if (vimit == this->vtxInfo_.end())
166  {
167  // insert into the map
168  this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
169  // remember this vertex' index
170  vertex_indices[i] = vertex_index;
171  // increase the current index
172  vertex_index++;
173  }
174  else
175  {
176  // only remember the vertex' index
177  vertex_indices[i] = vimit->second->idx;
178  }
179  }
180 
181  // flip cell if necessary
182  {
183  switch (int(dim))
184  {
185  case 0 :
186  break;
187  case 1 :
188  {
189  // The following test only works if the zero-th coordinate is the
190  // one that defines the orientation. A sufficient condition for
191  // this is dimworld == 1
192  /* assert(dimworld==1); */
193  bool elementNormalDirection =
194  (geometry.corner(1)[0] < geometry.corner(0)[0]);
195  if ( positiveNormalDirection_ != elementNormalDirection )
196  {
197  std::swap(vertex_indices[0], vertex_indices[1]);
198  std::swap(vertex_numbers[0], vertex_numbers[1]);
199  }
200  break;
201  }
202  case 2 :
203  {
204  Dune::FieldVector<ctype, dimworld>
205  v0 = geometry.corner(1),
206  v1 = geometry.corner(2);
207  v0 -= geometry.corner(0);
208  v1 -= geometry.corner(0);
209  ctype normal_sign = v0[0]*v1[1] - v0[1]*v1[0];
210  bool elementNormalDirection = (normal_sign < 0);
211  if ( positiveNormalDirection_ != elementNormalDirection )
212  {
213  std::cout << "swap\n";
214  if (elmt.type().isCube())
215  {
216  for (int i = 0; i < (1<<dim); i+=2)
217  {
218  // swap i and i+1
219  std::swap(vertex_indices[i], vertex_indices[i+1]);
220  std::swap(vertex_numbers[i], vertex_numbers[i+1]);
221  }
222  } else if (elmt.type().isSimplex()) {
223  std::swap(vertex_indices[0], vertex_indices[1]);
224  std::swap(vertex_numbers[0], vertex_numbers[1]);
225  } else {
226  DUNE_THROW(Dune::Exception, "Unexpected Geometrytype");
227  }
228  }
229  break;
230  }
231  }
232  }
233 
234  // add a new face to the temporary collection
235  temp_faces.push_back(SubEntityInfo(eindex,0,elmt.type()));
236  element_index++;
237  for (unsigned int i=0; i<numCorners; i++) {
238  temp_faces.back().corners[i].idx = vertex_indices[i];
239  // remember the vertices' numbers in parent element's vertices
240  temp_faces.back().corners[i].num = vertex_numbers[i];
241  }
242 
243  }
244  } // end loop over elements
245 
246  // allocate the array for the face specific information...
247  this->subEntities_.resize(element_index);
248  // ...and fill in the data from the temporary containers
249  copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
250 
251  // now first write the array with the coordinates...
252  this->coords_.resize(this->vtxInfo_.size());
253  typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
254  for (; it1 != this->vtxInfo_.end(); ++it1)
255  {
256  // get a pointer to the associated info object
257  CoordinateInfo* current = &this->coords_[it1->second->idx];
258  // store this coordinates index // NEEDED?
259  current->index = it1->second->idx;
260  // store the vertex' index for the index2vertex mapping
261  current->vtxindex = it1->first;
262  // store the vertex' coordinates under the associated index
263  // in the coordinates array
264 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
265  const auto vtx = this->grid().entity(it1->second->p);
266 #else
267  const auto vtxPtr = this->grid().entityPointer(it1->second->p);
268  const auto& vtx = *vtxPtr;
269 #endif
270  current->coord = vtx.geometry().corner(0);
271  }
272 
273 }
274 
275 } // namespace GridGlue
276 
277 } // namespace Dune
278 
279 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM0EXTRACTOR_HH
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim0extractor.hh:50
const GridView gv_
the grid object to extract the surface from
Definition: extractor.hh:206
Extractor< GV, 0 >::ElementInfo ElementInfo
Definition: codim0extractor.hh:58
void clear()
delete everything build up so far and free the memory
Definition: extractor.hh:252
Definition: codim0extractor.hh:35
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 base class
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim0extractor.hh:54
const bool & positiveNormalDirection() const
Definition: codim0extractor.hh:94
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
static const Dune::PartitionIteratorType PType
Definition: codim0extractor.hh:53
const Grid & grid() const
Definition: extractor.hh:380
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim0extractor.hh:48
Extractor< GV, 0 >::SubEntityInfo SubEntityInfo
Definition: codim0extractor.hh:57
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim0extractor.hh:49
Extractor< GV, 0 >::ctype ctype
Definition: codim0extractor.hh:42
Extractor< GV, 0 >::IndexType IndexType
Definition: codim0extractor.hh:45
Definition: extractor.hh:53
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim0extractor.hh:51
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:30
Extractor< GV, 0 >::CoordinateInfo CoordinateInfo
Definition: codim0extractor.hh:60
Codim0Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim0extractor.hh:85
VertexInfoMap vtxInfo_
a map enabling faster access to vertices and coordinates
Definition: extractor.hh:221
bool positiveNormalDirection_
Definition: codim0extractor.hh:97
Geometry geometry(unsigned int index) const
Get world geometry of the extracted face.
Base class for predicates selecting the part of a grid to be extracted.
CellMapper cellMapper_
Definition: extractor.hh:230
Extractor< GV, 0 >::VertexInfo VertexInfo
Definition: codim0extractor.hh:59
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
GV::Traits::template Codim< dim >::EntityPointer VertexPtr
Definition: codim0extractor.hh:47
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:47
bool & positiveNormalDirection()
Definition: codim0extractor.hh:93
ElementInfoMap elmtInfo_
a map enabling faster access to elements and faces
Definition: extractor.hh:228
Extractor< GV, 0 >::VertexInfoMap VertexInfoMap
Definition: codim0extractor.hh:61