dune-foamgrid  2.8-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 #include <functional>
15 
16 #include <dune/common/deprecated.hh>
17 #include <dune/common/version.hh>
18 // TODO remove header and macro after release Dune 2.8
19 #define DUNE_FUNCTION_HH_SILENCE_DEPRECATION // silence deprecation warning from <dune/common/function.hh>
20 #include <dune/common/function.hh>
21 #include <dune/common/fvector.hh>
22 #include <dune/common/hash.hh>
23 
24 #if DUNE_VERSION_EQUAL(DUNE_COMMON, 2, 7)
25 #include <dune/common/to_unique_ptr.hh>
26 #endif
27 
28 #include <dune/grid/common/gridfactory.hh>
30 
31 namespace Dune {
32 
35 template <int dimgrid, int dimworld, class ct>
37  : public GridFactoryInterface<FoamGrid<dimgrid, dimworld, ct> >
38  {
40  using ctype = ct;
41 
42  public:
43 
44 #if DUNE_VERSION_LT(DUNE_COMMON, 2, 7)
46 #elif DUNE_VERSION_LT(DUNE_COMMON, 2, 8)
47  using GridPtrType = ToUniquePtr<FoamGrid<dimgrid, dimworld, ctype>>;
48 #else
49  using GridPtrType = std::unique_ptr<FoamGrid<dimgrid, dimworld,ctype>>;
50 #endif
51 
52  public:
53 
56  : factoryOwnsGrid_(true)
57  {
59  grid_->entityImps_.resize(1);
60  }
61 
72  : factoryOwnsGrid_(false)
73  {
74  grid_ = grid;
75  grid_->entityImps_.resize(1);
76  }
77 
79  ~GridFactoryBase() override {
80  if (grid_ && factoryOwnsGrid_)
81  delete grid_;
82  }
83 
85  void insertVertex(const FieldVector<ctype,dimworld>& pos) override {
86  std::get<0>(grid_->entityImps_[0]).emplace_back(
87  0, // level
88  pos, // position
89  grid_->getNextFreeId()
90  );
91  vertexArray_.push_back(&*std::get<0>(grid_->entityImps_[0]).rbegin());
92  }
93 
96  unsigned int
97  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<0>::Entity &entity) const override
98  {
99  return entity.impl().target_->leafIndex_;
100  }
101 
104  unsigned int
105  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<dimgrid>::Entity &vertex) const override
106  {
107  return vertex.impl().target_->leafIndex_;
108  }
109 
112  unsigned int
113  insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection& intersection) const override
114  {
115  return intersection.boundarySegmentIndex();
116  }
117 
118  protected:
119 
120  // Pointer to the grid being built
122 
123  // True if the factory allocated the grid itself, false if the
124  // grid was handed over from the outside
126 
128  std::vector<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*> vertexArray_;
129 
131  unsigned int boundarySegmentCounter_ = 0;
132  };
133 
134 template <int dimgrid, int dimworld, class ct>
135  class GridFactory<FoamGrid<dimgrid, dimworld, ct> >
136  : public GridFactoryBase<dimgrid, dimworld, ct>
137  {};
138 
141 template <int dimworld, class ct>
142  class GridFactory<FoamGrid<1, dimworld, ct> >
143  : public GridFactoryBase<1, dimworld, ct>
144  {
146  enum {dimgrid = 1};
147 
149  using ctype = ct;
150 
152  template<int mydim>
154 
155  using GridPtrType = typename GridFactoryBase<1, dimworld, ct>::GridPtrType;
156 
157  public:
158 
160 
162  GridFactoryBase<1,dimworld,ctype>(grid)
163  {}
164 
168  void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
169  {
170  boundarySegmentIndices_[ vertices[0] ] = this->boundarySegmentCounter_++;
171  }
172 
176  void insertBoundarySegment(const std::vector<unsigned int>& vertices,
177  const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
178  {
179  DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
180  }
181 
184  bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
185  {
186  if ( !intersection.boundary() || intersection.inside().level() != 0 )
187  return false;
188 
189  const auto& vertex = intersection.inside().template subEntity<1>(intersection.indexInInside());
190  const auto& it = boundarySegmentIndices_.find( this->insertionIndex(vertex) );
191  return (it != boundarySegmentIndices_.end());
192  }
193 
198  void insertElement(const GeometryType& type,
199  const std::vector<unsigned int>& vertices) override
200  {
201  assert(type.isLine());
202  // insert new element
203  std::get<1>(this->grid_->entityImps_[0]).emplace_back(
204  this->vertexArray_[vertices[0]],
205  this->vertexArray_[vertices[1]],
206  0, // level
207  this->grid_->getNextFreeId()
208  );
209  }
210 
211 // disable the deprecated interface when using Dune > 2.8
212 #if DUNE_VERSION_LTE(DUNE_COMMON, 2, 8)
213 DUNE_NO_DEPRECATED_BEGIN
219  [[deprecated("Signature with VirtualFunction is deprecated and will be removed after Dune 2.8. Use signature with std::function.")]]
220  void insertElement(const GeometryType& type,
221  const std::vector<unsigned int>& vertices,
222  const std::shared_ptr<VirtualFunction<FieldVector<ctype,dimgrid>,FieldVector<ctype,dimworld> > >& elementParametrization) override
223  {
224  assert(type.isLine());
225  EntityImp<1> newElement(this->vertexArray_[vertices[0]],
226  this->vertexArray_[vertices[1]],
227  0, // level
228  this->grid_->getNextFreeId());
229  // save the pointer to the element parametrization
230  newElement.elementParametrization_ =
231  [elementParametrization](const FieldVector<ctype,dimgrid>& x){
232  FieldVector<ctype,dimworld> y;
233  elementParametrization->evaluate(x, y);
234  return y;
235  };
236 
237  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
238  }
239 DUNE_NO_DEPRECATED_END
240 #endif
241 
247  void insertElement(const GeometryType& type,
248  const std::vector<unsigned int>& vertices,
249  std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dimgrid>)> elementParametrization)
250 #if DUNE_VERSION_GT(DUNE_COMMON, 2, 7)
251  override
252 #endif
253  {
254  assert(type.isLine());
255  EntityImp<1> newElement(this->vertexArray_[vertices[0]],
256  this->vertexArray_[vertices[1]],
257  0, // level
258  this->grid_->getNextFreeId());
259 
260  newElement.elementParametrization_ = elementParametrization;
261  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
262  }
263 
268  GridPtrType createGrid() override
269  {
270  // Prevent a crash when this method is called twice in a row
271  // You never know who may do this...
272  if (this->grid_==nullptr)
273  return nullptr;
274 
275  // make facets (vertices) know about the element
276  for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
277  {
278  element.vertex_[0]->elements_.push_back(&element);
279  element.vertex_[1]->elements_.push_back(&element);
280  }
281 
282  // Create the index sets
283  this->grid_->setIndices();
284 
285  // ////////////////////////////////////////////////
286  // Set the boundary ids
287  // ////////////////////////////////////////////////
288 
289  for (auto& facet : std::get<0>(this->grid_->entityImps_[0]))
290  {
291  if (facet.elements_.size() == 1) // if boundary facet
292  {
293  const auto it = boundarySegmentIndices_.find( facet.vertex_[0]->leafIndex_ );
294  if (it != boundarySegmentIndices_.end())
295  facet.boundarySegmentIndex_ = it->second;
296  else
297  facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
298  }
299  }
300 
301  // ////////////////////////////////////////////////
302  // Hand over the new grid
303  // ////////////////////////////////////////////////
304 
305  Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
306  tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
307  this->grid_ = nullptr;
308  return GridPtrType(tmp);
309  }
310 
311  private:
312  std::unordered_map<unsigned int, unsigned int> boundarySegmentIndices_;
313  };
314 
317 template <int dimworld, class ct>
318  class GridFactory<FoamGrid<2, dimworld, ct> >
319  : public GridFactoryBase<2, dimworld, ct>
320  {
322  enum {dimgrid = 2};
323 
325  using ctype = ct;
326 
328  template<int mydim>
330 
331  using GridPtrType = typename GridFactoryBase<2, dimworld, ct>::GridPtrType;
332 
333  public:
334 
336 
338  GridFactoryBase<2, dimworld, ctype>(grid)
339  {}
340 
344  void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
345  {
346  std::array<unsigned int, 2> vertexIndices {{ vertices[0], vertices[1] }};
347 
348  // sort the indices
349  if ( vertexIndices[0] > vertexIndices[1] )
350  std::swap( vertexIndices[0], vertexIndices[1] );
351 
352  boundarySegmentIndices_[ vertexIndices ] = this->boundarySegmentCounter_++;
353  }
354 
358  void insertBoundarySegment(const std::vector<unsigned int>& vertices,
359  const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
360  {
361  DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
362  }
363 
366  bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
367  {
368  if ( !intersection.boundary() || intersection.inside().level() != 0 )
369  return false;
370 
371  // Get the vertices of the intersection
372  const auto refElement = ReferenceElements<ctype, dimgrid>::general(intersection.inside().type());
373 
374  const int subIdx0 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/0, dimgrid);
375  const auto vertex0 = intersection.inside().template subEntity<2>( subIdx0 );
376  const int subIdx1 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/1, dimgrid);
377  const auto vertex1 = intersection.inside().template subEntity<2>( subIdx1 );
378 
379  std::array<unsigned int, 2> vertexIndices {{
380  this->insertionIndex( vertex0 ),
381  this->insertionIndex( vertex1 )
382  }};
383 
384  // sort the indices
385  if ( vertexIndices[0] > vertexIndices[1] )
386  std::swap( vertexIndices[0], vertexIndices[1] );
387 
388  const auto& it = boundarySegmentIndices_.find( vertexIndices );
389  return (it != boundarySegmentIndices_.end());
390  }
391 
396  void insertElement(const GeometryType& type,
397  const std::vector<unsigned int>& vertices) override {
398 
399  assert(type.isTriangle());
400 
401  EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
402  newElement.vertex_[0] = this->vertexArray_[vertices[0]];
403  newElement.vertex_[1] = this->vertexArray_[vertices[1]];
404  newElement.vertex_[2] = this->vertexArray_[vertices[2]];
405 
406  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
407  }
408 
409 // disable the deprecated interface when using Dune > 2.8
410 #if DUNE_VERSION_LTE(DUNE_COMMON, 2, 8)
411 DUNE_NO_DEPRECATED_BEGIN
417  [[deprecated("Signature with VirtualFunction is deprecated and will be removed after Dune 2.8. Use signature with std::function.")]]
418  void insertElement(const GeometryType& type,
419  const std::vector<unsigned int>& vertices,
420  const std::shared_ptr<VirtualFunction<FieldVector<ctype,dimgrid>,FieldVector<ctype,dimworld> > >& elementParametrization) override
421  {
422  assert(type.isTriangle());
423  EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
424  newElement.vertex_[0] = this->vertexArray_[vertices[0]];
425  newElement.vertex_[1] = this->vertexArray_[vertices[1]];
426  newElement.vertex_[2] = this->vertexArray_[vertices[2]];
427  // save the pointer to the element parametrization
428  newElement.elementParametrization_ =
429  [elementParametrization](const FieldVector<ctype,dimgrid>& x){
430  FieldVector<ctype,dimworld> y;
431  elementParametrization->evaluate(x, y);
432  return y;
433  };
434 
435  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
436  }
437 DUNE_NO_DEPRECATED_END
438 #endif
439 
445  void insertElement(const GeometryType& type,
446  const std::vector<unsigned int>& vertices,
447  std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dimgrid>)> elementParametrization)
448 #if DUNE_VERSION_GT(DUNE_COMMON, 2, 7)
449  override
450 #endif
451  {
452  assert(type.isTriangle());
453  EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
454  newElement.vertex_[0] = this->vertexArray_[vertices[0]];
455  newElement.vertex_[1] = this->vertexArray_[vertices[1]];
456  newElement.vertex_[2] = this->vertexArray_[vertices[2]];
457  newElement.elementParametrization_ = elementParametrization;
458  std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
459  }
460 
464  GridPtrType createGrid() override
465  {
466  // Prevent a crash when this method is called twice in a row
467  // You never know who may do this...
468  if (this->grid_==nullptr)
469  return nullptr;
470 
471  // ////////////////////////////////////////////////////
472  // Create the edges
473  // ////////////////////////////////////////////////////
474 
475  // For fast retrieval: a map from pairs of vertices to the edge that connects them
476  std::unordered_map<std::pair<const EntityImp<0>*, const EntityImp<0>*>, EntityImp<1>*,
477  HashPair<const EntityImp<0>*>> edgeMap;
478  for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
479  {
480  const auto refElement = ReferenceElements<ctype, dimgrid>::general(element.type());
481 
482  // Loop over all edges of this element
483  for (std::size_t i=0; i<element.facet_.size(); ++i)
484  {
485  // get pointer to the edge (insert edge if it's not in the entity list yet)
486  auto edge = [&](){
487  // Get two vertices of the potential edge
488  auto v0 = element.vertex_[refElement.subEntity(i, 1, 0, 2)];
489  auto v1 = element.vertex_[refElement.subEntity(i, 1, 1, 2)];
490  // sort pointers
491  if ( v0 > v1 )
492  std::swap( v0, v1 );
493 
494  // see if the edge was already inserted
495  // the pointer pair hash is symmetric so we don't have to check for pair(v1,v0)
496  auto e = edgeMap.find(std::make_pair(v0, v1));
497  if (e != edgeMap.end())
498  return e->second;
499 
500  // edge was not inserted yet: insert the edge now
501  std::get<1>(this->grid_->entityImps_[0]).emplace_back(v0, v1, /*level=*/0, this->grid_->getNextFreeId());
502  // insert it into the map of inserted edges
503  auto newEdge = &*std::get<1>(this->grid_->entityImps_[0]).rbegin();
504  edgeMap.insert(std::make_pair(std::make_pair(v0, v1), newEdge));
505  return newEdge;
506  }();
507 
508  // make element know about the edge
509  element.facet_[i] = edge;
510 
511  // make edge know about the element
512  edge->elements_.push_back(&element);
513  }
514  }
515 
516  // Create the index sets
517  this->grid_->setIndices();
518 
519 
520  // ////////////////////////////////////////////////
521  // Set the boundary ids
522  // ////////////////////////////////////////////////
523 
524  for (auto& facet : std::get<1>(this->grid_->entityImps_[0]))
525  {
526  if (facet.elements_.size() == 1) // if boundary facet
527  {
528  std::array<unsigned int, 2> vertexIndices {{ facet.vertex_[0]->leafIndex_, facet.vertex_[1]->leafIndex_ }};
529  // sort the indices
530  if ( vertexIndices[0] > vertexIndices[1] )
531  std::swap( vertexIndices[0], vertexIndices[1] );
532 
533  const auto it = boundarySegmentIndices_.find( vertexIndices );
534  if (it != boundarySegmentIndices_.end())
535  facet.boundarySegmentIndex_ = it->second;
536  else
537  facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
538  }
539  }
540 
541  // ////////////////////////////////////////////////
542  // Hand over the new grid
543  // ////////////////////////////////////////////////
544 
545  Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
546  tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
547  this->grid_ = nullptr;
548  return GridPtrType(tmp);
549  }
550 
551  private:
552  template<class T, class U = T>
553  struct HashPair {
554  std::size_t operator() (const std::pair<T, U>& a) const {
555  std::size_t seed = 0;
556  hash_combine(seed, a.first);
557  hash_combine(seed, a.second);
558  return seed;
559  }
560  };
561 
562  struct HashUIntArray {
563  std::size_t operator() (const std::array<unsigned int, 2>& a) const {
564  return hash_range(a.begin(), a.end());
565  }
566  };
567 
568  std::unordered_map<std::array<unsigned int, 2>, unsigned int, HashUIntArray> boundarySegmentIndices_;
569  };
570 
571 } // end namespace Dune
572 
573 #endif
The FoamGrid class.
Definition: dgffoam.cc:6
Specialization of the generic GridFactory for FoamGrid<dimgrid, dimworld>
Definition: foamgridfactory.hh:38
~GridFactoryBase() override
Destructor.
Definition: foamgridfactory.hh:79
void insertVertex(const FieldVector< ctype, dimworld > &pos) override
Insert a vertex into the coarse grid.
Definition: foamgridfactory.hh:85
unsigned int insertionIndex(const typename FoamGrid< dimgrid, dimworld, ctype >::LeafIntersection &intersection) const override
Obtain a boundary's insertion index.
Definition: foamgridfactory.hh:113
GridFactoryBase(FoamGrid< dimgrid, dimworld, ctype > *grid)
Constructor for a given grid object.
Definition: foamgridfactory.hh:71
std::vector< FoamGridEntityImp< 0, dimgrid, dimworld, ctype > * > vertexArray_
Array containing all vertices.
Definition: foamgridfactory.hh:128
unsigned int boundarySegmentCounter_
Counter that creates the boundary segment indices.
Definition: foamgridfactory.hh:131
bool factoryOwnsGrid_
Definition: foamgridfactory.hh:125
FoamGrid< dimgrid, dimworld, ctype > * grid_
Definition: foamgridfactory.hh:121
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:105
GridFactoryBase()
Default constructor.
Definition: foamgridfactory.hh:55
std::unique_ptr< FoamGrid< dimgrid, dimworld, ctype > > GridPtrType
Definition: foamgridfactory.hh:49
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:97
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices, std::function< FieldVector< ctype, dimworld >(FieldVector< ctype, dimgrid >)> elementParametrization)
Insert a parametrized element into the coarse grid.
Definition: foamgridfactory.hh:247
GridPtrType createGrid() override
Finalize grid creation and hand over the grid.
Definition: foamgridfactory.hh:268
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:176
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:184
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:168
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices) override
Insert an element into the coarse grid.
Definition: foamgridfactory.hh:198
GridFactory(FoamGrid< 1, dimworld, ctype > *grid)
Definition: foamgridfactory.hh:161
GridFactory()
Definition: foamgridfactory.hh:159
GridFactory()
Definition: foamgridfactory.hh:335
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:344
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:358
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:366
GridPtrType createGrid() override
Finalize grid creation and hand over the grid The receiver takes responsibility of the memory allocat...
Definition: foamgridfactory.hh:464
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices, std::function< FieldVector< ctype, dimworld >(FieldVector< ctype, dimgrid >)> elementParametrization)
Insert a parametrized element into the coarse grid.
Definition: foamgridfactory.hh:445
void insertElement(const GeometryType &type, const std::vector< unsigned int > &vertices) override
Insert an element into the coarse grid.
Definition: foamgridfactory.hh:396
GridFactory(FoamGrid< 2, dimworld, ctype > *grid)
Definition: foamgridfactory.hh:337
The actual entity implementation.
Definition: foamgridvertex.hh:47