dune-foamgrid  2.7-git
foamgridintersections.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_INTERSECTIONS_HH
4 #define DUNE_FOAMGRID_INTERSECTIONS_HH
5 
9 #include <memory>
10 #include <algorithm>
11 
12 #include <dune/grid/common/intersection.hh>
13 
18 
19 namespace Dune {
20 
21 template <class GridImp>
22 class FoamGridLevelIntersectionIterator;
23 
24 template <class GridImp>
25 class FoamGridLeafIntersectionIterator;
26 
27 template <class GridImp>
28 class FoamGridLevelIntersection;
29 
30 template <class GridImp>
31 class FoamGridLeafIntersection;
32 
33 
37 template<class GridImp>
39 {
40  enum {dimgrid = GridImp::dimension};
41  enum {dimworld = GridImp::dimensionworld};
42 
43  // The type used to store coordinates
44  typedef typename GridImp::ctype ctype;
45  typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
46 
47  // The iterators need access to all members
48  friend class FoamGridLevelIntersectionIterator<GridImp>;
49  friend class FoamGridLeafIntersectionIterator<GridImp>;
50 
51  // As weĺl as the implementations
52  friend class FoamGridLevelIntersection<GridImp>;
53  friend class FoamGridLeafIntersection<GridImp>;
54 
55  template<typename, typename>
56  friend class Dune::Intersection;
57 
59  : center_(nullptr)
60  , facetIndex_(-1)
61  , neighbor_(FoamGridNullIteratorFactory<dimgrid, dimworld, ctype>::null())
62  {}
63 
72  int facet)
73  : center_(center), facetIndex_(facet), neighbor_(FoamGridNullIteratorFactory<dimgrid, dimworld, ctype>::null())
74  {}
75 
76 public:
77 
78  typedef typename GridImp::template Codim<0>::Entity Entity;
79 
82  Entity inside() const {
84  }
85 
86 
89  Entity outside(/*std::size_t neighborIndex = 0*/) const {
90  // Return the 'other' element on the current facet
92  }
93 
95  bool equals(const FoamGridIntersection<GridImp>& i) const {
96  return center_==i.center_ && neighbor_ == i.neighbor_ && facetIndex_ == i.facetIndex_;
97  }
98 
101  bool boundary () const {
102  return center_->facet_[facetIndex_]->elements_.size()==1;
103  }
104 
106  int boundarySegmentIndex () const
107  {
108  return center_->facet_[facetIndex_]->boundarySegmentIndex();
109  }
110 
112  GeometryType type () const
113  {
114  return Dune::GeometryTypes::simplex(dimgrid-1);
115  }
116 
118  int indexInInside () const
119  {
120  assert(facetIndex_ < center_->corners());
121  return facetIndex_;
122  }
123 
124  virtual int indexInOutside(std::size_t neighborIndex = 0) const=0;
125 
127  FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dimgrid-1>& local) const
128  {
129  // The intersection normal is a vector that is orthogonal to the element normal
130  // and to the intersection itself.
131 
132  // only works for triangles and lines
133  assert(center_->type().isTriangle() || center_->type().isLine());
134 
135  // Compute vertices
136  const auto refElement = ReferenceElements<ctype, dimgrid>::general(center_->type());
137 
138  // facet vertices (vertex in 1d, edge in 2d)
139  int v0 = refElement.subEntity(facetIndex_, 1, 0, dimgrid);
140 
141  if(dimgrid == 2) {
142  //second vertex only exists for an edge
143  int v1 = refElement.subEntity(facetIndex_, 1, 1, dimgrid);
144  // opposite vertex
145  int v2 = (v1+1)%3;
146  if (v2==v0)
147  v2 = (v0+1)%3;
148  assert(v2!=v0 and v2!=v1);
149 
150  FieldVector<ctype, dimworld> facet = center_->vertex_[v0]->pos_ - center_->vertex_[v1]->pos_;
151  FieldVector<ctype, dimworld> otherEdge = center_->vertex_[v2]->pos_ - center_->vertex_[v1]->pos_;
152 
153  //Cross product of facet and otherEdge is a scaled element normal
154  FieldVector<ctype, dimworld> scaledElementNormal;
155 
156  if(dimworld == 3) //dimworld==3
157  {
158  scaledElementNormal[0] = facet[1]*otherEdge[2] - facet[2]*otherEdge[1];
159  scaledElementNormal[1] = facet[2]*otherEdge[0] - facet[0]*otherEdge[2];
160  scaledElementNormal[2] = facet[0]*otherEdge[1] - facet[1]*otherEdge[0];
161  outerNormal_[0] = facet[1]*scaledElementNormal[2] - facet[2]*scaledElementNormal[1];
162  outerNormal_[1] = facet[2]*scaledElementNormal[0] - facet[0]*scaledElementNormal[2];
163  outerNormal_[2] = facet[0]*scaledElementNormal[1] - facet[1]*scaledElementNormal[0];
164  }
165  else //dimworld==2
166  {
167  outerNormal_[0] = facet[1];
168  outerNormal_[1] = -facet[0];
169  }
170 
171  //Check if scaled EdgeNormal is inner normal, if yes flip
172  otherEdge = center_->vertex_[v0]->pos_ - center_->vertex_[v2]->pos_;
173  if(otherEdge*outerNormal_ < 0)
174  outerNormal_ *= -1.0;
175 
176  return outerNormal_;
177  }
178  else //dimgrid ==1
179  {
180  int v1 = 1-v0;
181  assert(v0!=v1);
182  //automatically has right orientation
183  outerNormal_ = center_->vertex_[v0]->pos_ - center_->vertex_[v1]->pos_;
184 
185  return outerNormal_;
186  }
187  DUNE_THROW(GridError, "Non-existing grid dimension requested! Has to be 1 or 2!");
188  }
189 
191  FieldVector<ctype, dimworld> integrationOuterNormal (const FieldVector<ctype, dimgrid-1>& local) const
192  {
193 
194  if(dimgrid == 2)
195  {
196  const auto refElement = ReferenceElements<ctype, dimgrid>::general(center_->type());
197 
198  // facet vertices
199  int v0 = refElement.subEntity(facetIndex_, 1, 0, dimgrid);
200  int v1 = refElement.subEntity(facetIndex_, 1, 1, dimgrid);
201 
202  //facet length
203  ctype facetLength = (center_->vertex_[v0]->pos_- center_->vertex_[v1]->pos_).two_norm();
204 
205  FieldVector<ctype, dimworld> integrationOuterNormal_ = unitOuterNormal(local);
206  integrationOuterNormal_ *= facetLength;
207  return integrationOuterNormal_;
208  }
209  else //dimgrid == 1
210  {
211  integrationOuterNormal_ = this->unitOuterNormal(local);
212  return integrationOuterNormal_;
213  }
214  DUNE_THROW(GridError, "Non-existing grid dimension requested! Has to be 1 or 2!");
215  }
216 
218  FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dimgrid-1>& local) const
219  {
220  unitOuterNormal_ = this->outerNormal(local);
221  unitOuterNormal_ /= unitOuterNormal_.two_norm();
222  return unitOuterNormal_;
223  }
224 
226  FieldVector<ctype, dimworld> centerUnitOuterNormal () const
227  {
228  return unitOuterNormal(FieldVector<ctype,dimgrid-1>(0.5));
229  }
230  private:
231 
233  mutable FieldVector<ctype, dimworld> outerNormal_;
234  mutable FieldVector<ctype, dimworld> unitOuterNormal_;
235  mutable FieldVector<ctype, dimworld> integrationOuterNormal_;
236 
237  protected:
238 
240 
243 
245  typename std::vector<const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*>::const_iterator neighbor_;
246 
247 };
248 
249 template<class GridImp>
251  : public FoamGridIntersection<GridImp>
252 {
253  friend class FoamGridLevelIntersectionIterator<GridImp>;
254 
255  template<typename, typename>
256  friend class Dune::Intersection;
257 
259 
260  enum {dimgrid = GridImp::dimension};
261  enum{ dimworld = GridImp::dimensionworld };
262 
263  // Geometry is a CachedMultiLinearGeometry
264  typedef typename GridImp::template Codim<1>::Geometry Geometry;
265  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
266 
267  typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
268  typedef typename GridImp::Traits::template Codim<1>::LocalGeometryImpl LocalGeometryImpl;
269 
270  typedef typename GridImp::ctype ctype;
271 
272  FoamGridLevelIntersection(const FoamGridEntityImp<dimgrid, dimgrid ,dimworld, ctype>* center, std::size_t facet)
273  : FoamGridIntersection<GridImp>(center, facet)
274  {}
275 
276 public:
277 
279  bool conforming () const {
280  // FoamGrid level intersections are always conforming
281  return true;
282  }
283 
285  int indexInOutside (std::size_t neighborIndex = 0) const {
286  return std::find((*this->neighbor_)->facet_.begin(), (*this->neighbor_)->facet_.end(),
287  this->center_->facet_[this->facetIndex_])
288  - (*this->neighbor_)->facet_.begin();
289  }
290 
292  bool neighbor () const {
293  return this->neighbor_!=neighborEnd_;
294  }
295 
300  LocalGeometry geometryInInside () const
301  {
302  std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
303 
304  // Get reference element
305  const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
306 
307  for (int idx = 0; idx < dimgrid; ++idx)
308  coordinates[idx] = refElement.position(refElement.subEntity(this->facetIndex_, 1, idx, dimgrid), dimgrid);
309 
310  geometryInInside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
311 
312  return LocalGeometry(*geometryInInside_);
313  }
314 
318  LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const {
319 
320  // Get two vertices of the intersection
321  const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
322 
323  std::array<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*, dimgrid> vtx;
324 
325  for (int idx = 0; idx < dimgrid; ++idx)
326  vtx[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)];
327 
328  std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
329 
330  // Find the intersection vertices in local numbering of the outside element
331  // That way we get the local orientation correctly.
332  const auto refElementOther = ReferenceElements<ctype, dimgrid>::general((*this->neighbor_)->type());
333 
334  for (std::size_t j=0; j<dimgrid; j++)
335  for (int i=0; i<refElementOther.size(dimgrid); i++)
336  if (vtx[j] == (*this->neighbor_)->vertex_[refElementOther.subEntity(0, 0, i, dimgrid)])
337  coordinates[j] = refElement.position(refElement.subEntity(0, 0, i, dimgrid), dimgrid);
338 
339  geometryInOutside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
340 
341  return LocalGeometry(*geometryInOutside_);
342  }
343 
346  Geometry geometry () const {
347 
348  std::vector<FieldVector<ctype, dimworld> > coordinates(dimgrid);
349 
350  // Get two vertices of the intersection
351  const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
352 
353  for (std::size_t idx = 0; idx < dimgrid; ++idx)
354  coordinates[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)]->pos_;
355 
356  geometry_ = std::make_shared<GeometryImpl>(this->type(), coordinates);
357 
358  return Geometry(*geometry_);
359  }
360 
361 
362  private:
364  typename std::vector<const FoamGridEntityImp<dimgrid, dimgrid ,dimworld, ctype>*>::const_iterator neighborEnd_;
366  mutable std::shared_ptr<GeometryImpl> geometry_;
367  mutable std::shared_ptr<LocalGeometryImpl> geometryInInside_;
368  mutable std::shared_ptr<LocalGeometryImpl> geometryInOutside_;
369 
370 };
371 
372 
373 
374 
383 template<class GridImp>
385  : public FoamGridIntersection<GridImp>
386 {
387 
388  friend class FoamGridLeafIntersectionIterator<GridImp>;
389 
390  template<typename, typename>
391  friend class Dune::Intersection;
392 
394 
395  enum {dimworld = GridImp::dimensionworld};
396  enum {dimgrid = GridImp::dimension};
397 
398  typedef typename GridImp::ctype ctype;
399 
400  typedef typename GridImp::template Codim<1>::Geometry Geometry;
401  typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
402 
403  typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
404  typedef typename GridImp::Traits::template Codim<1>::LocalGeometryImpl LocalGeometryImpl;
405 
406  FoamGridLeafIntersection(const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* center,
407  std::size_t facet)
408  : FoamGridIntersection<GridImp>(center, facet)
409  {}
410 
411 public:
412 
414  bool conforming () const {
415  // FoamGrid leaf? intersections are always? conforming
416  return true;
417  }
418 
420  int indexInOutside (std::size_t neighborIndex = 0) const
421  {
422  // get our facet object
423  auto* facet = this->center_->facet_[this->facetIndex_];
424 
425  // if the neighbor level is greater than the facet level
426  // walk until our facet's father that has the same level as the neighbor
427  // and look for the father facet in the neighbor to get the index
428  while(facet->level() > (*this->neighbor_)->level())
429  {
430  assert(facet->father_ != nullptr);
431  facet = facet->father_;
432  }
433 
434  if (facet->level() == (*this->neighbor_)->level())
435  {
436  return std::distance((*this->neighbor_)->facet_.begin(),
437  std::find((*this->neighbor_)->facet_.begin(),
438  (*this->neighbor_)->facet_.end(),
439  facet)
440  );
441 
442  }
443 
444  // the neighbor level is smaller than the facet level
445  // find the neighbor's facet that has the center element in it's element list
446  else
447  {
448  return std::distance((*this->neighbor_)->facet_.begin(),
449  std::find_if((*this->neighbor_)->facet_.begin(),
450  (*this->neighbor_)->facet_.end(),
451  [this](const auto& facet) -> bool {
452  for (const auto& e : facet->elements_)
453  if (this->center_ == e)
454  return true;
455  return false;
456  })
457  );
458  }
459  }
460 
465  LocalGeometry geometryInInside () const
466  {
467  std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
468 
469  // Get reference element
470  const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
471 
472  for (std::size_t idx = 0; idx < dimgrid; ++idx)
473  coordinates[idx] = refElement.position(refElement.subEntity(this->facetIndex_, 1, idx, dimgrid), dimgrid);
474 
475  geometryInInside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
476 
477  return LocalGeometry(*geometryInInside_);
478  }
479 
484  LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const
485  {
486  // Get reference element
487  const auto refElement = Dune::ReferenceElements<ctype, dimgrid>::general(this->center_->type());
488 
489  // Map the global position of the intersection vertices to the local position in the neighbor
490  std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
491  for (std::size_t vIdx = 0; vIdx < dimgrid; ++vIdx)
492  {
493  const auto vIdxInElement = refElement.subEntity(this->facetIndex_, 1, vIdx, dimgrid);
494  coordinates[vIdx] = (*this->neighbor_)->globalToLocal(this->center_->vertex_[vIdxInElement]->pos_);
495  }
496 
497  geometryInOutside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
498 
499  return LocalGeometry(*geometryInOutside_);
500  }
501 
504  Geometry geometry () const {
505 
506  std::vector<FieldVector<ctype, dimworld> > coordinates(dimgrid);
507 
508  // Get two vertices of the intersection
509  const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
510 
511  for (std::size_t idx = 0; idx < dimgrid; ++idx)
512  coordinates[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)]->pos_;
513 
514  geometry_ = std::make_shared<GeometryImpl>(this->type(), coordinates);
515 
516  return Geometry(*geometry_);
517  }
518 
520  bool neighbor () const {
522  }
523 
524 private:
525 
527  mutable std::shared_ptr<GeometryImpl> geometry_;
528  mutable std::shared_ptr<LocalGeometryImpl> geometryInInside_;
529  mutable std::shared_ptr<LocalGeometryImpl> geometryInOutside_;
530 };
531 
532 
533 } // namespace Dune
534 
535 #endif
The FoamGridEntity class.
The FoamGridGeometry class.
The null iterator factory for intersections.
Definition: dgffoam.cc:6
Element specialization of FoamGridEntityImp. Element is a grid entity of topological codimension 0 an...
Definition: foamgridelements.hh:18
The implementation of entities in a FoamGrid.
Definition: foamgridentity.hh:54
Definition: foamgridintersectioniterators.hh:239
Iterator over all element neighborsMesh entities of codimension 0 ("elements") allow to visit all nei...
Definition: foamgridintersectioniterators.hh:28
Definition: foamgridintersections.hh:252
bool conforming() const
Return true if this is a conforming intersection.
Definition: foamgridintersections.hh:279
int indexInOutside(std::size_t neighborIndex=0) const
local number of codim 1 entity in neighbor where intersection is contained
Definition: foamgridintersections.hh:285
LocalGeometry geometryInOutside(std::size_t neighborIndex=0) const
Definition: foamgridintersections.hh:318
bool neighbor() const
return true if across the facet a neighbor on this level exists
Definition: foamgridintersections.hh:292
Geometry geometry() const
Definition: foamgridintersections.hh:346
LocalGeometry geometryInInside() const
Definition: foamgridintersections.hh:300
Iterator over all element neighborsMesh entities of codimension 0 ("elements") allow to visit all nei...
Definition: foamgridintersections.hh:386
int indexInOutside(std::size_t neighborIndex=0) const
local number of codim 1 entity in neighbor where intersection is contained
Definition: foamgridintersections.hh:420
LocalGeometry geometryInInside() const
Definition: foamgridintersections.hh:465
LocalGeometry geometryInOutside(std::size_t neighborIndex=0) const
Definition: foamgridintersections.hh:484
bool conforming() const
Return true if this is a conforming intersection.
Definition: foamgridintersections.hh:414
bool neighbor() const
return true if across the facet a neighbor on this level exists
Definition: foamgridintersections.hh:520
Geometry geometry() const
Definition: foamgridintersections.hh:504
Base class of all intersections within FoamGrid.
Definition: foamgridintersections.hh:39
virtual int indexInOutside(std::size_t neighborIndex=0) const =0
GridImp::template Codim< 0 >::Entity Entity
Definition: foamgridintersections.hh:78
FieldVector< ctype, dimworld > outerNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return outer normal
Definition: foamgridintersections.hh:127
int facetIndex_
Count on which facet we are lookin' at.
Definition: foamgridintersections.hh:242
Entity inside() const
Definition: foamgridintersections.hh:82
std::vector< const FoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * >::const_iterator neighbor_
Iterator to the other neighbor of the intersection.
Definition: foamgridintersections.hh:245
const FoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * center_
Definition: foamgridintersections.hh:239
FieldVector< ctype, dimworld > integrationOuterNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return outer normal multiplied by the integration element
Definition: foamgridintersections.hh:191
bool boundary() const
return true if intersection is with boundary.
Definition: foamgridintersections.hh:101
int boundarySegmentIndex() const
return information about the Boundary
Definition: foamgridintersections.hh:106
int indexInInside() const
local number of codim 1 entity in self where intersection is contained in
Definition: foamgridintersections.hh:118
Entity outside() const
Definition: foamgridintersections.hh:89
GeometryType type() const
Geometry type of an intersection.
Definition: foamgridintersections.hh:112
FieldVector< ctype, dimworld > centerUnitOuterNormal() const
return unit outer normal at the intersection center
Definition: foamgridintersections.hh:226
bool equals(const FoamGridIntersection< GridImp > &i) const
equality
Definition: foamgridintersections.hh:95
FieldVector< ctype, dimworld > unitOuterNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return unit outer normal
Definition: foamgridintersections.hh:218
Definition: foamgridnulliteratorfactory.hh:16
static std::vector< const FoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * >::const_iterator null()
Definition: foamgridnulliteratorfactory.hh:18