dune-foamgrid 2.8.0
Loading...
Searching...
No Matches
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
19namespace Dune {
20
21template <class GridImp>
22class FoamGridLevelIntersectionIterator;
23
24template <class GridImp>
25class FoamGridLeafIntersectionIterator;
26
27template <class GridImp>
28class FoamGridLevelIntersection;
29
30template <class GridImp>
31class FoamGridLeafIntersection;
32
33
37template<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
76public:
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
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
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
249template<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
276public:
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
383template<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
411public:
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
524private:
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 null iterator factory for intersections.
The FoamGridEntity class.
The FoamGridGeometry class.
Definition: dgffoam.cc:6
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
int facetIndex_
Count on which facet we are lookin' at.
Definition: foamgridintersections.hh:242
std::vector< constFoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * >::const_iterator neighbor_
Iterator to the other neighbor of the intersection.
Definition: foamgridintersections.hh:245
Entity inside() const
Definition: foamgridintersections.hh:82
FieldVector< ctype, dimworld > integrationOuterNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return outer normal multiplied by the integration element
Definition: foamgridintersections.hh:191
const FoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * center_
Definition: foamgridintersections.hh:239
FieldVector< ctype, dimworld > outerNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return outer normal
Definition: foamgridintersections.hh:127
FieldVector< ctype, dimworld > centerUnitOuterNormal() const
return unit outer normal at the intersection center
Definition: foamgridintersections.hh:226
FieldVector< ctype, dimworld > unitOuterNormal(const FieldVector< ctype, dimgrid-1 > &local) const
return unit outer normal
Definition: foamgridintersections.hh:218
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
bool equals(const FoamGridIntersection< GridImp > &i) const
equality
Definition: foamgridintersections.hh:95
Definition: foamgridnulliteratorfactory.hh:16
static std::vector< constFoamGridEntityImp< dimgrid, dimgrid, dimworld, ctype > * >::const_iterator null()
Definition: foamgridnulliteratorfactory.hh:18
The actual entity implementation.
Definition: foamgridvertex.hh:47