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