dune-foamgrid  2.7-git
foamgridfactory.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set ts=8 sw=4 et sts=4:
3 #ifndef DUNE_FOAMGRID_FACTORY_HH
4 #define DUNE_FOAMGRID_FACTORY_HH
5 
11 #include <vector>
12 #include <memory>
13 #include <unordered_map>
14 
15 #include <dune/common/version.hh>
16 #include <dune/common/function.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/hash.hh>
19 #include <dune/common/to_unique_ptr.hh>
20 
21 #include <dune/grid/common/gridfactory.hh>
23 
24 namespace Dune {
25 
28 template <int dimgrid, int dimworld, class ct>
30  : public GridFactoryInterface<FoamGrid<dimgrid, dimworld, ct> >
31  {
33  using ctype = ct;
34 
35  public:
36 
39  : factoryOwnsGrid_(true)
40  {
42  grid_->entityImps_.resize(1);
43  }
44 
55  : factoryOwnsGrid_(false)
56  {
57  grid_ = grid;
58  grid_->entityImps_.resize(1);
59  }
60 
62  ~GridFactoryBase() override {
63  if (grid_ && factoryOwnsGrid_)
64  delete grid_;
65  }
66 
68  void insertVertex(const FieldVector<ctype,dimworld>& pos) override {
69  std::get<0>(grid_->entityImps_[0]).emplace_back(
70  0, // level
71  pos, // position
72  grid_->getNextFreeId()
73  );
74  vertexArray_.push_back(&*std::get<0>(grid_->entityImps_[0]).rbegin());
75  }
76 
79  unsigned int
80  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<0>::Entity &entity) const override
81  {
82  return entity.impl().target_->leafIndex_;
83  }
84 
87  unsigned int
88  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<dimgrid>::Entity &vertex) const override
89  {
90  return vertex.impl().target_->leafIndex_;
91  }
92 
95  unsigned int
96  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection& intersection) const override
97  {
98  return intersection.boundarySegmentIndex();
99  }
100 
101  protected:
102 
103  // Pointer to the grid being built
105 
106  // True if the factory allocated the grid itself, false if the
107  // grid was handed over from the outside
109 
111  std::vector<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*> vertexArray_;
112 
114  unsigned int boundarySegmentCounter_ = 0;
115  };
116 
117 template <int dimgrid, int dimworld, class ct>
118  class GridFactory<FoamGrid<dimgrid, dimworld, ct> >
119  : public GridFactoryBase<dimgrid, dimworld, ct>
120  {};
121 
124 template <int dimworld, class ct>
125  class GridFactory<FoamGrid<1, dimworld, ct> >
126  : public GridFactoryBase<1, dimworld, ct>
127  {
129  enum {dimgrid = 1};
130 
132  using ctype = ct;
133 
135  template<int mydim>
137 
138  public:
139 
141 
143  GridFactoryBase<1,dimworld,ctype>(grid)
144  {}
145 
149  void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
150  {
151  boundarySegmentIndices_[ vertices[0] ] = this->boundarySegmentCounter_++;
152  }
153 
157  void insertBoundarySegment(const std::vector<unsigned int>& vertices,
158  const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
159  {
160  DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
161  }
162 
165  bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
166  {
167  if ( !intersection.boundary() || intersection.inside().level() != 0 )
168  return false;
169 
170  const auto& vertex = intersection.inside().template subEntity<1>(intersection.indexInInside());
171  const auto& it = boundarySegmentIndices_.find( this->insertionIndex(vertex) );
172  return (it != boundarySegmentIndices_.end());
173  }
174 
179  void insertElement(const GeometryType& type,
180  const std::vector<unsigned int>& vertices) override
181  {
182  assert(type.isLine());
183  // insert new element
184  std::get<1>(this->grid_->entityImps_[0]).emplace_back(
185  this->vertexArray_[vertices[0]],
186  this->vertexArray_[vertices[1]],
187  0, // level
188  this->grid_->getNextFreeId()
189  );
190  }
191 
197  void insertElement(const GeometryType& type,
198  const std::vector<unsigned int>& vertices,
199  const std::shared_ptr<VirtualFunction<FieldVector<ctype,dimgrid>,FieldVector<ctype,dimworld> > >& elementParametrization) override
200  {
201  assert(type.isLine());
202  EntityImp<1> newElement(this->vertexArray_[vertices[0]],
203  this->vertexArray_[vertices[1]],
204  0, // level
205  this->grid_->getNextFreeId());
206  // save the pointer to the element parametrization
207  newElement.elementParametrization_ = elementParametrization;
208 
209  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
210  }
211 
216 #if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
217  FoamGrid<1, dimworld, ctype>* createGrid() override {
218 #else
219  ToUniquePtr<FoamGrid<1, dimworld, ctype>> createGrid() override {
220 #endif
221  // Prevent a crash when this method is called twice in a row
222  // You never know who may do this...
223  if (this->grid_==nullptr)
224  return nullptr;
225 
226  // make facets (vertices) know about the element
227  for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
228  {
229  element.vertex_[0]->elements_.push_back(&element);
230  element.vertex_[1]->elements_.push_back(&element);
231  }
232 
233  // Create the index sets
234  this->grid_->setIndices();
235 
236  // ////////////////////////////////////////////////
237  // Set the boundary ids
238  // ////////////////////////////////////////////////
239 
240  for (auto& facet : std::get<0>(this->grid_->entityImps_[0]))
241  {
242  if (facet.elements_.size() == 1) // if boundary facet
243  {
244  const auto it = boundarySegmentIndices_.find( facet.vertex_[0]->leafIndex_ );
245  if (it != boundarySegmentIndices_.end())
246  facet.boundarySegmentIndex_ = it->second;
247  else
248  facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
249  }
250  }
251 
252  // ////////////////////////////////////////////////
253  // Hand over the new grid
254  // ////////////////////////////////////////////////
255 
256  Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
257  tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
258  this->grid_ = nullptr;
259  return tmp;
260  }
261 
262  private:
263  std::unordered_map<unsigned int, unsigned int> boundarySegmentIndices_;
264  };
265 
268 template <int dimworld, class ct>
269  class GridFactory<FoamGrid<2, dimworld, ct> >
270  : public GridFactoryBase<2, dimworld, ct>
271  {
273  enum {dimgrid = 2};
274 
276  using ctype = ct;
277 
279  template<int mydim>
281 
282  public:
283 
285 
287  GridFactoryBase<2, dimworld, ctype>(grid)
288  {}
289 
293  void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
294  {
295  std::array<unsigned int, 2> vertexIndices {{ vertices[0], vertices[1] }};
296 
297  // sort the indices
298  if ( vertexIndices[0] > vertexIndices[1] )
299  std::swap( vertexIndices[0], vertexIndices[1] );
300 
301  boundarySegmentIndices_[ vertexIndices ] = this->boundarySegmentCounter_++;
302  }
303 
307  void insertBoundarySegment(const std::vector<unsigned int>& vertices,
308  const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
309  {
310  DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
311  }
312 
315  bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
316  {
317  if ( !intersection.boundary() || intersection.inside().level() != 0 )
318  return false;
319 
320  // Get the vertices of the intersection
321  const auto refElement = ReferenceElements<ctype, dimgrid>::general(intersection.inside().type());
322 
323  const int subIdx0 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/0, dimgrid);
324  const auto vertex0 = intersection.inside().template subEntity<2>( subIdx0 );
325  const int subIdx1 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/1, dimgrid);
326  const auto vertex1 = intersection.inside().template subEntity<2>( subIdx1 );
327 
328  std::array<unsigned int, 2> vertexIndices {{
329  this->insertionIndex( vertex0 ),
330  this->insertionIndex( vertex1 )
331  }};
332 
333  // sort the indices
334  if ( vertexIndices[0] > vertexIndices[1] )
335  std::swap( vertexIndices[0], vertexIndices[1] );
336 
337  const auto& it = boundarySegmentIndices_.find( vertexIndices );
338  return (it != boundarySegmentIndices_.end());
339  }
340 
345  void insertElement(const GeometryType& type,
346  const std::vector<unsigned int>& vertices) override {
347 
348  assert(type.isTriangle());
349 
350  EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
351  newElement.vertex_[0] = this->vertexArray_[vertices[0]];
352  newElement.vertex_[1] = this->vertexArray_[vertices[1]];
353  newElement.vertex_[2] = this->vertexArray_[vertices[2]];
354 
355  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
356  }
357 
363  void insertElement(const GeometryType& type,
364  const std::vector<unsigned int>& vertices,
365  const std::shared_ptr<VirtualFunction<FieldVector<ctype,dimgrid>,FieldVector<ctype,dimworld> > >& elementParametrization) override
366  {
367  assert(type.isTriangle());
368  EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
369  newElement.vertex_[0] = this->vertexArray_[vertices[0]];
370  newElement.vertex_[1] = this->vertexArray_[vertices[1]];
371  newElement.vertex_[2] = this->vertexArray_[vertices[2]];
372  // save the pointer to the element parametrization
373  newElement.elementParametrization_ = elementParametrization;
374 
375  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
376  }
377 
381 #if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
382  FoamGrid<dimgrid, dimworld, ctype>* createGrid() override {
383 #else
384  ToUniquePtr<FoamGrid<dimgrid, dimworld, ctype>> createGrid() override {
385 #endif
386  // Prevent a crash when this method is called twice in a row
387  // You never know who may do this...
388  if (this->grid_==nullptr)
389  return nullptr;
390 
391  // ////////////////////////////////////////////////////
392  // Create the edges
393  // ////////////////////////////////////////////////////
394 
395  // For fast retrieval: a map from pairs of vertices to the edge that connects them
396  std::unordered_map<std::pair<const EntityImp<0>*, const EntityImp<0>*>, EntityImp<1>*,
397  HashPair<const EntityImp<0>*>> edgeMap;
398  for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
399  {
400  const auto refElement = ReferenceElements<ctype, dimgrid>::general(element.type());
401 
402  // Loop over all edges of this element
403  for (std::size_t i=0; i<element.facet_.size(); ++i)
404  {
405  // get pointer to the edge (insert edge if it's not in the entity list yet)
406  auto edge = [&](){
407  // Get two vertices of the potential edge
408  auto v0 = element.vertex_[refElement.subEntity(i, 1, 0, 2)];
409  auto v1 = element.vertex_[refElement.subEntity(i, 1, 1, 2)];
410  // sort pointers
411  if ( v0 > v1 )
412  std::swap( v0, v1 );
413 
414  // see if the edge was already inserted
415  // the pointer pair hash is symmetric so we don't have to check for pair(v1,v0)
416  auto e = edgeMap.find(std::make_pair(v0, v1));
417  if (e != edgeMap.end())
418  return e->second;
419 
420  // edge was not inserted yet: insert the edge now
421  std::get<1>(this->grid_->entityImps_[0]).emplace_back(v0, v1, /*level=*/0, this->grid_->getNextFreeId());
422  // insert it into the map of inserted edges
423  auto newEdge = &*std::get<1>(this->grid_->entityImps_[0]).rbegin();
424  edgeMap.insert(std::make_pair(std::make_pair(v0, v1), newEdge));
425  return newEdge;
426  }();
427 
428  // make element know about the edge
429  element.facet_[i] = edge;
430 
431  // make edge know about the element
432  edge->elements_.push_back(&element);
433  }
434  }
435 
436  // Create the index sets
437  this->grid_->setIndices();
438 
439 
440  // ////////////////////////////////////////////////
441  // Set the boundary ids
442  // ////////////////////////////////////////////////
443 
444  for (auto& facet : std::get<1>(this->grid_->entityImps_[0]))
445  {
446  if (facet.elements_.size() == 1) // if boundary facet
447  {
448  std::array<unsigned int, 2> vertexIndices {{ facet.vertex_[0]->leafIndex_, facet.vertex_[1]->leafIndex_ }};
449  // sort the indices
450  if ( vertexIndices[0] > vertexIndices[1] )
451  std::swap( vertexIndices[0], vertexIndices[1] );
452 
453  const auto it = boundarySegmentIndices_.find( vertexIndices );
454  if (it != boundarySegmentIndices_.end())
455  facet.boundarySegmentIndex_ = it->second;
456  else
457  facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
458  }
459  }
460 
461  // ////////////////////////////////////////////////
462  // Hand over the new grid
463  // ////////////////////////////////////////////////
464 
465  Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
466  tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
467  this->grid_ = nullptr;
468  return tmp;
469  }
470 
471  private:
472  template<class T, class U = T>
473  struct HashPair {
474  std::size_t operator() (const std::pair<T, U>& a) const {
475  std::size_t seed = 0;
476  hash_combine(seed, a.first);
477  hash_combine(seed, a.second);
478  return seed;
479  }
480  };
481 
482  struct HashUIntArray {
483  std::size_t operator() (const std::array<unsigned int, 2>& a) const {
484  return hash_range(a.begin(), a.end());
485  }
486  };
487 
488  std::unordered_map<std::array<unsigned int, 2>, unsigned int, HashUIntArray> boundarySegmentIndices_;
489  };
490 
491 } // end namespace Dune
492 
493 #endif
The FoamGrid class.
Definition: dgffoam.cc:6
Specialization of the generic GridFactory for FoamGrid<dimgrid, dimworld>
Definition: foamgridfactory.hh:31
~GridFactoryBase() override
Destructor.
Definition: foamgridfactory.hh:62
void insertVertex(const FieldVector< ctype, dimworld > &pos) override
Insert a vertex into the coarse grid.
Definition: foamgridfactory.hh:68
unsigned int insertionIndex(const typename FoamGrid< dimgrid, dimworld, ctype >::LeafIntersection &intersection) const override
Obtain a boundary's insertion index.
Definition: foamgridfactory.hh:96
GridFactoryBase(FoamGrid< dimgrid, dimworld, ctype > *grid)
Constructor for a given grid object.
Definition: foamgridfactory.hh:54
std::vector< FoamGridEntityImp< 0, dimgrid, dimworld, ctype > * > vertexArray_
Array containing all vertices.
Definition: foamgridfactory.hh:111
unsigned int boundarySegmentCounter_
Counter that creates the boundary segment indices.
Definition: foamgridfactory.hh:114
bool factoryOwnsGrid_
Definition: foamgridfactory.hh:108
FoamGrid< dimgrid, dimworld, ctype > * grid_
Definition: foamgridfactory.hh:104
unsigned int insertionIndex(const typename FoamGrid< dimgrid, dimworld, ctype >::Traits::template Codim< dimgrid >::Entity &vertex) const override
Obtain a vertex' insertion index.
Definition: foamgridfactory.hh:88
GridFactoryBase()
Default constructor.
Definition: foamgridfactory.hh:38
unsigned int insertionIndex(const typename FoamGrid< dimgrid, dimworld, ctype >::Traits::template Codim< 0 >::Entity &entity) const override
Obtain an element's insertion index.
Definition: foamgridfactory.hh:80
ToUniquePtr< FoamGrid< 1, dimworld, ctype > > createGrid() override
Finalize grid creation and hand over the grid.
Definition: foamgridfactory.hh:219
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices, const std::shared_ptr< VirtualFunction< FieldVector< ctype, dimgrid >, FieldVector< ctype, dimworld > > > &elementParametrization) override
Insert a parametrized element into the coarse grid.
Definition: foamgridfactory.hh:197
void insertBoundarySegment(const std::vector< unsigned int > &vertices, const std::shared_ptr< BoundarySegment< dimgrid, dimworld > > &boundarySegment) override
Insert a boundary segment and the boundary segment geometry This influences the ordering of the bound...
Definition: foamgridfactory.hh:157
bool wasInserted(const typename FoamGrid< dimgrid, dimworld, ctype >::LeafIntersection &intersection) const override
Return true if leaf intersection was inserted as boundary segment.
Definition: foamgridfactory.hh:165
void insertBoundarySegment(const std::vector< unsigned int > &vertices) override
Insert a boundary segment. This is only needed if you want to control the numbering of the boundary s...
Definition: foamgridfactory.hh:149
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices) override
Insert an element into the coarse grid.
Definition: foamgridfactory.hh:179
GridFactory(FoamGrid< 1, dimworld, ctype > *grid)
Definition: foamgridfactory.hh:142
GridFactory()
Definition: foamgridfactory.hh:140
GridFactory()
Definition: foamgridfactory.hh:284
void insertBoundarySegment(const std::vector< unsigned int > &vertices) override
Insert a boundary segment. This is only needed if you want to control the numbering of the boundary s...
Definition: foamgridfactory.hh:293
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices, const std::shared_ptr< VirtualFunction< FieldVector< ctype, dimgrid >, FieldVector< ctype, dimworld > > > &elementParametrization) override
Insert a parametrized element into the coarse grid.
Definition: foamgridfactory.hh:363
ToUniquePtr< FoamGrid< dimgrid, dimworld, ctype > > createGrid() override
Finalize grid creation and hand over the grid The receiver takes responsibility of the memory allocat...
Definition: foamgridfactory.hh:384
void insertBoundarySegment(const std::vector< unsigned int > &vertices, const std::shared_ptr< BoundarySegment< dimgrid, dimworld > > &boundarySegment) override
Insert a boundary segment (== a line) and the boundary segment geometry This influences the ordering ...
Definition: foamgridfactory.hh:307
bool wasInserted(const typename FoamGrid< dimgrid, dimworld, ctype >::LeafIntersection &intersection) const override
Return true if leaf intersection was inserted as boundary segment.
Definition: foamgridfactory.hh:315
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices) override
Insert an element into the coarse grid.
Definition: foamgridfactory.hh:345
GridFactory(FoamGrid< 2, dimworld, ctype > *grid)
Definition: foamgridfactory.hh:286
The actual entity implementation.
Definition: foamgridvertex.hh:47