dune-grid-glue  2.4-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 
23 #include <dune/grid/common/mcmgmapper.hh>
24 
25 #include "extractor.hh"
26 #include "extractorpredicate.hh"
27 
28 namespace Dune {
29 
30  namespace GridGlue {
31 
32 template<typename GV>
33 class Codim0Extractor : public Extractor<GV,0>
34 {
35 
36 public:
37 
38  /* 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  typedef typename Extractor<GV,0>::ctype ctype;
44 
45  typedef typename GV::Traits::template Codim<dim>::EntityPointer VertexPtr;
46  typedef typename GV::Traits::template Codim<dim>::Entity Vertex;
47  typedef typename GV::Traits::template Codim<0>::EntityPointer ElementPtr;
48  typedef typename GV::Traits::template Codim<0>::Entity Element;
49 
50  static const Dune::PartitionIteratorType PType = Dune::Interior_Partition;
51  typedef typename GV::Traits::template Codim<0>::template Partition<PType>::Iterator ElementIter;
52 
53  // import typedefs from base class
59 
65  Codim0Extractor(const GV& gv, const ExtractorPredicate<GV,0>& descr)
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(descr);
71  }
72 
74  const bool & positiveNormalDirection() const { return positiveNormalDirection_; }
75 
76 protected:
78 private:
79  void update(const ExtractorPredicate<GV,0>& descr);
80 };
81 
82 
83 template<typename GV>
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 (ElementIter elit = this->gv_.template begin<0, PType>();
103  elit != this->gv_.template end<0, PType>(); ++elit)
104  {
105  const auto& elmt = *elit;
106  const auto geometry = elmt.geometry();
107 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
108  IndexType eindex = this->cellMapper_.index(elmt);
109 #else
110  IndexType eindex = this->cellMapper_.map(elmt);
111 #endif
112 
113  // only do sth. if this element is "interesting"
114  // implicit cast is done automatically
115  if (descr.contains(elmt,0))
116  {
117  // add an entry to the element info map, the index will be set properly later
118  this->elmtInfo_[eindex] = new ElementInfo(element_index, elmt, 1);
119 
120 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
121  unsigned int numCorners = elmt.subEntities(dim);
122 #else
123  unsigned int numCorners = elmt.template count<dim>();
124 #endif
125  unsigned int vertex_indices[numCorners]; // index in global vector
126  unsigned int vertex_numbers[numCorners]; // index in parent entity
127 
128  // try for each of the faces vertices whether it is already inserted or not
129  for (unsigned int i = 0; i < numCorners; ++i)
130  {
131  vertex_numbers[i] = i;
132 
133  // get the vertex pointer and the index from the index set
134 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
135  const Vertex vertex = elit->template subEntity<dim>(vertex_numbers[i]);
136 #else
137  const VertexPtr vertexPtr = elit->template subEntity<dim>(vertex_numbers[i]);
138  const Vertex &vertex = *vertexPtr;
139 #endif
140  IndexType vindex = this->gv_.indexSet().template index<dim>(vertex);
141 
142  // if the vertex is not yet inserted in the vertex info map
143  // it is a new one -> it will be inserted now!
144  typename VertexInfoMap::iterator vimit = this->vtxInfo_.find(vindex);
145  if (vimit == this->vtxInfo_.end())
146  {
147  // insert into the map
148  this->vtxInfo_[vindex] = new VertexInfo(vertex_index, vertex);
149  // remember this vertex' index
150  vertex_indices[i] = vertex_index;
151  // increase the current index
152  vertex_index++;
153  }
154  else
155  {
156  // only remember the vertex' index
157  vertex_indices[i] = vimit->second->idx;
158  }
159  }
160 
161  // flip cell if necessary
162  {
163  switch (int(dim))
164  {
165  case 0 :
166  break;
167  case 1 :
168  {
169  // The following test only works if the zero-th coordinate is the
170  // one that defines the orientation. A sufficient condition for
171  // this is dimworld == 1
172  /* assert(dimworld==1); */
173  bool elementNormalDirection =
174  (geometry.corner(1)[0] < geometry.corner(0)[0]);
175  if ( positiveNormalDirection_ != elementNormalDirection )
176  {
177  std::swap(vertex_indices[0], vertex_indices[1]);
178  std::swap(vertex_numbers[0], vertex_numbers[1]);
179  }
180  break;
181  }
182  case 2 :
183  {
184  Dune::FieldVector<ctype, dimworld>
185  v0 = geometry.corner(1),
186  v1 = geometry.corner(2);
187  v0 -= geometry.corner(0);
188  v1 -= geometry.corner(0);
189  ctype normal_sign = v0[0]*v1[1] - v0[1]*v1[0];
190  bool elementNormalDirection = (normal_sign < 0);
191  if ( positiveNormalDirection_ != elementNormalDirection )
192  {
193  std::cout << "swap\n";
194  if (elmt.type().isCube())
195  {
196  for (int i = 0; i < (1<<dim); i+=2)
197  {
198  // swap i and i+1
199  std::swap(vertex_indices[i], vertex_indices[i+1]);
200  std::swap(vertex_numbers[i], vertex_numbers[i+1]);
201  }
202  } else if (elmt.type().isSimplex()) {
203  std::swap(vertex_indices[0], vertex_indices[1]);
204  std::swap(vertex_numbers[0], vertex_numbers[1]);
205  } else {
206  DUNE_THROW(Dune::Exception, "Unexpected Geometrytype");
207  }
208  }
209  break;
210  }
211  }
212  }
213 
214  // add a new face to the temporary collection
215  temp_faces.push_back(SubEntityInfo(eindex,0,elmt.type()));
216  element_index++;
217  for (int i=0; i<numCorners; i++) {
218  temp_faces.back().corners[i].idx = vertex_indices[i];
219  // remember the vertices' numbers in parent element's vertices
220  temp_faces.back().corners[i].num = vertex_numbers[i];
221  }
222 
223  }
224  } // end loop over elements
225 
226  // allocate the array for the face specific information...
227  this->subEntities_.resize(element_index);
228  // ...and fill in the data from the temporary containers
229  copy(temp_faces.begin(), temp_faces.end(), this->subEntities_.begin());
230 
231  // now first write the array with the coordinates...
232  this->coords_.resize(this->vtxInfo_.size());
233  typename VertexInfoMap::const_iterator it1 = this->vtxInfo_.begin();
234  for (; it1 != this->vtxInfo_.end(); ++it1)
235  {
236  // get a pointer to the associated info object
237  CoordinateInfo* current = &this->coords_[it1->second->idx];
238  // store this coordinates index // NEEDED?
239  current->index = it1->second->idx;
240  // store the vertex' index for the index2vertex mapping
241  current->vtxindex = it1->first;
242  // store the vertex' coordinates under the associated index
243  // in the coordinates array
244 #if DUNE_VERSION_NEWER(DUNE_GRID, 2, 4)
245  const auto vtx = this->grid().entity(it1->second->p);
246 #else
247  const auto vtxPtr = this->grid().entityPointer(it1->second->p);
248  const auto& vtx = *vtxPtr;
249 #endif
250  current->coord = vtx.geometry().corner(0);
251  }
252 
253 }
254 
255 } // namespace GridGlue
256 
257 } // namespace Dune
258 
259 #endif // DUNE_GRIDGLUE_EXTRACTORS_CODIM0EXTRACTOR_HH
GV::Traits::template Codim< dim >::EntityPointer VertexPtr
Definition: codim0extractor.hh:45
static const Dune::PartitionIteratorType PType
Definition: codim0extractor.hh:50
Codim0Extractor(const GV &gv, const ExtractorPredicate< GV, 0 > &descr)
Constructor.
Definition: codim0extractor.hh:65
Extractor< GV, 0 >::VertexInfo VertexInfo
Definition: codim0extractor.hh:56
Definition: codim0extractor.hh:33
const bool & positiveNormalDirection() const
Definition: codim0extractor.hh:74
bool & positiveNormalDirection()
Definition: codim0extractor.hh:73
Extractor< GV, 0 >::CoordinateInfo CoordinateInfo
Definition: codim0extractor.hh:57
Extractor< GV, 0 >::ctype ctype
Definition: codim0extractor.hh:40
virtual bool contains(const Element &element, unsigned int subentity) const =0
Return true if a subentity should be extracted.
Extractor< GV, 0 >::SubEntityInfo SubEntityInfo
Definition: codim0extractor.hh:54
Extractor< GV, 0 >::IndexType IndexType
Definition: codim0extractor.hh:43
Base class for subentity-selecting predicates.
Definition: extractorpredicate.hh:30
Extractor< GV, 0 >::ElementInfo ElementInfo
Definition: codim0extractor.hh:55
GV::Traits::template Codim< 0 >::EntityPointer ElementPtr
Definition: codim0extractor.hh:47
Definition: gridglue.hh:33
GV::Traits::template Codim< 0 >::template Partition< PType >::Iterator ElementIter
Definition: codim0extractor.hh:51
Base class for predicates selecting the part of a grid to be extracted.
GV::Traits::template Codim< 0 >::Entity Element
Definition: codim0extractor.hh:48
Extractor< GV, 0 >::VertexInfoMap VertexInfoMap
Definition: codim0extractor.hh:58
GV::Traits::template Codim< dim >::Entity Vertex
Definition: codim0extractor.hh:46
Provides codimension-independent methods for grid extraction.
Definition: extractor.hh:47
extractor base class
bool positiveNormalDirection_
Definition: codim0extractor.hh:77