dune-grid-glue  2.8-git
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 
29 namespace Dune {
30 
31  namespace GridGlue {
32 
36 template<typename GV>
37 class Codim0Extractor : public Extractor<GV,0>
38 {
39 
40 public:
41 
42  /* 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 */
44  typedef typename Extractor<GV,0>::ctype ctype;
48 
49  typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
50  typedef typename GV::Traits::template Codim<0>::Entity Element;
51  typedef std::function<bool(const Element&, unsigned int subentity)> Predicate;
52 
53  // import typedefs from base class
59 
65  Codim0Extractor(const GV& gv, const Predicate& predicate)
66  : Extractor<GV,0>(gv), positiveNormalDirection_(false)
67  {
68  std::cout << "This is Codim0Extractor on a <"
69  << GV::dimension << "," << GV::dimensionworld << "> grid!" << std::endl;
70  update(predicate);
71  }
72 
74  const bool & positiveNormalDirection() const { return positiveNormalDirection_; }
75 
76 protected:
78 private:
79  void update(const Predicate& predicate);
80 };
81 
82 
83 template<typename GV>
84 void Codim0Extractor<GV>::update(const Predicate& predicate)
85 {
86  // In this first pass iterate over all entities of codim 0.
87  // Get its corner vertices, find resp. store them together with their associated index,
88  // and remember the indices of the corners.
89 
90  // free everything there is in this object
91  this->clear();
92 
93  // several counter for consecutive indexing are needed
94  size_t element_index = 0;
95  size_t vertex_index = 0;
96 
97  // a temporary container where newly acquired face
98  // information can be stored at first
99  std::deque<SubEntityInfo> temp_faces;
100 
101  // iterate over all codim 0 elements on the grid
102  for (const auto& elmt : elements(this->gv_, Partitions::interior))
103  {
104  const auto geometry = elmt.geometry();
105  IndexType eindex = this->cellMapper_.index(elmt);
106 
107  // only do sth. if this element is "interesting"
108  // implicit cast is done automatically
109  if (predicate(elmt,0))
110  {
111  // add an entry to the element info map, the index will be set properly later
112  this->elmtInfo_.emplace(eindex, ElementInfo(element_index, elmt, 1));
113 
114  unsigned int numCorners = elmt.subEntities(dim);
115  unsigned int vertex_indices[numCorners]; // index in global vector
116  unsigned int vertex_numbers[numCorners]; // index in parent entity
117 
118  // try for each of the faces vertices whether it is already inserted or not
119  for (unsigned int i = 0; i < numCorners; ++i)
120  {
121  vertex_numbers[i] = i;
122 
123  // get the vertex pointer and the index from the index set
124  const Vertex vertex = elmt.template subEntity<dim>(vertex_numbers[i]);
125  IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
126 
127  // if the vertex is not yet inserted in the vertex info map
128  // it is a new one -> it will be inserted now!
129  typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
130  if (vimit == this->vtxInfo_.end())
131  {
132  // insert into the map
133  this->vtxInfo_.emplace(vindex, VertexInfo(vertex_index, vertex));
134  // remember this vertex' index
135  vertex_indices[i] = vertex_index;
136  // increase the current index
137  vertex_index++;
138  }
139  else
140  {
141  // only remember the vertex' index
142  vertex_indices[i] = vimit->second.idx;
143  }
144  }
145 
146  // flip cell if necessary
147  {
148  switch (int(dim))
149  {
150  case 0 :
151  break;
152  case 1 :
153  {
154  // The following test only works if the zero-th coordinate is the
155  // one that defines the orientation. A sufficient condition for
156  // this is dimworld == 1
157  /* assert(dimworld==1); */
158  bool elementNormalDirection =
159  (geometry.corner(1)[0] < geometry.corner(0)[0]);
160  if ( positiveNormalDirection_ != elementNormalDirection )
161  {
162  std::swap(vertex_indices[0], vertex_indices[1]);
163  std::swap(vertex_numbers[0], vertex_numbers[1]);
164  }
165  break;
166  }
167  case 2 :
168  {
169  Dune::FieldVector<ctype, dimworld>
170  v0 = geometry.corner(1),
171  v1 = geometry.corner(2);
172  v0 -= geometry.corner(0);
173  v1 -= geometry.corner(0);
174  ctype normal_sign = v0[0]*v1[1] - v0[1]*v1[0];
175  bool elementNormalDirection = (normal_sign < 0);
176  if ( positiveNormalDirection_ != elementNormalDirection )
177  {
178  std::cout << "swap\n";
179  if (elmt.type().isCube())
180  {
181  for (int i = 0; i < (1<<dim); i+=2)
182  {
183  // swap i and i+1
184  std::swap(vertex_indices[i], vertex_indices[i+1]);
185  std::swap(vertex_numbers[i], vertex_numbers[i+1]);
186  }
187  } else if (elmt.type().isSimplex()) {
188  std::swap(vertex_indices[0], vertex_indices[1]);
189  std::swap(vertex_numbers[0], vertex_numbers[1]);
190  } else {
191  DUNE_THROW(Dune::Exception, "Unexpected Geometrytype");
192  }
193  }
194  break;
195  }
196  }
197  }
198 
199  // add a new face to the temporary collection
200  temp_faces.emplace_back(eindex, 0, elmt.type());
201  element_index++;
202  for (unsigned int i=0; i<numCorners; i++) {
203  temp_faces.back().corners[i].idx = vertex_indices[i];
204  // remember the vertices' numbers in parent element's vertices
205  temp_faces.back().corners[i].num = vertex_numbers[i];
206  }
207 
208  }
209  } // end loop over elements
210 
211  // allocate the array for the face specific information...
212  this->subEntities_.resize(element_index);
213  // ...and fill in the data from the temporary containers
214  copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
215 
216  // now first write the array with the coordinates...
217  this->coords_.resize(this->vtxInfo_.size());
218  for (const auto& vinfo : this->vtxInfo_)
219  {
220  // get a pointer to the associated info object
221  CoordinateInfo* current = &this->coords_[vinfo.second.idx];
222  // store this coordinates index // NEEDED?
223  current->index = vinfo.second.idx;
224  // store the vertex' index for the index2vertex mapping
225  current->vtxindex = vinfo.first;
226  // store the vertex' coordinates under the associated index
227  // in the coordinates array
228  const auto vtx = this->grid().entity(vinfo.second.p);
229  current->coord = vtx.geometry().corner(0);
230  }
231 
232 }
233 
234 } // namespace GridGlue
235 
236 } // namespace Dune
237 
238 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM0EXTRACTOR_HH
extractor base class
Definition: gridglue.hh:35
Definition: codim0extractor.hh:38
Extractor< GV, 0 >::IndexType IndexType
Definition: codim0extractor.hh:47
bool & positiveNormalDirection()
Definition: codim0extractor.hh:73
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim0extractor.hh:49
Extractor< GV, 0 >::CoordinateInfo CoordinateInfo
Definition: codim0extractor.hh:57
Extractor< GV, 0 >::VertexInfo VertexInfo
Definition: codim0extractor.hh:56
Extractor< GV, 0 >::ctype ctype
Definition: codim0extractor.hh:44
bool positiveNormalDirection_
Definition: codim0extractor.hh:77
Extractor< GV, 0 >::VertexInfoMap VertexInfoMap
Definition: codim0extractor.hh:58
std::function< bool(const Element &, unsigned int subentity)> Predicate
Definition: codim0extractor.hh:51
const bool & positiveNormalDirection() const
Definition: codim0extractor.hh:74
Codim0Extractor(const GV &gv, const Predicate &predicate)
Constructor.
Definition: codim0extractor.hh:65
Extractor< GV, 0 >::ElementInfo ElementInfo
Definition: codim0extractor.hh:55
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim0extractor.hh:50
Extractor< GV, 0 >::SubEntityInfo SubEntityInfo
Definition: codim0extractor.hh:54
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:44