dune-grid-glue 2.8.0
Loading...
Searching...
No Matches
intersection.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:
9#ifndef DUNE_GRIDGLUE_ADAPTER_INTERSECTION_HH
10#define DUNE_GRIDGLUE_ADAPTER_INTERSECTION_HH
11
12#include <algorithm>
13#include <optional>
14#include <tuple>
15
16#include <dune/common/deprecated.hh>
17#include <dune/common/version.hh>
18#include <dune/geometry/affinegeometry.hh>
19#include <dune/geometry/referenceelements.hh>
21
22#define ONLY_SIMPLEX_INTERSECTIONS
23
24namespace Dune {
25 namespace GridGlue {
26
27 // forward declaration
28 template<typename P0, typename P1>
29 class IntersectionIndexSet;
30
34 template<typename P0, typename P1>
36 {
37 public:
38 typedef ::Dune::GridGlue::GridGlue<P0, P1> GridGlue;
39
41
43 static constexpr int coorddim = GridGlue::dimworld;
44
45 private:
46 // intermediate quantities
47 template<int side>
48 static constexpr int dim() { return GridGlue::template GridView<side>::Grid::dimension - GridGlue::template GridPatch<side>::codim; }
49
50 public:
52 static constexpr int mydim = dim<0>() < dim<1>() ? dim<0>() : dim<1>();
53
54 template<int side>
55 using GridLocalGeometry = AffineGeometry<
56 typename GridGlue::template GridView<side>::ctype, mydim, GridGlue::template GridView<side>::dimension>;
57
58 using Grid0LocalGeometry [[deprecated("please use GridLocalGeometry<0> instead")]] = GridLocalGeometry<0>;
59 using Grid1LocalGeometry [[deprecated("please use GridLocalGeometry<1> instead")]] = GridLocalGeometry<1>;
60
61 template<int side>
62 using GridGeometry = AffineGeometry<
63 typename GridGlue::template GridView<side>::ctype, mydim, GridGlue::template GridView<side>::dimensionworld>;
64
65 using Grid0Geometry [[deprecated("please use GridGeometry<0> instead")]] = GridGeometry<0>;
66 using Grid1Geometry [[deprecated("please use GridGeometry<1> instead")]] = GridGeometry<1>;
67
68 template<int side>
69 using GridIndexType = typename GridGlue::template GridView<side>::IndexSet::IndexType;
70
71 using Grid0IndexType [[deprecated("please use GridIndexType<0> instead")]] = GridIndexType<0>;
72 using Grid1IndexType [[deprecated("please use GridIndexType<1> instead")]] = GridIndexType<1>;
73
75 IntersectionData(const GridGlue& glue, unsigned int mergeindex, unsigned int offset, bool grid0local, bool grid1local);
76
78 IntersectionData() = default;
79
80 /* Accessor Functions */
81
82 template<int side>
83 const GridLocalGeometry<side>& localGeometry(unsigned int parentId = 0) const
84 { return *std::get<side>(sideData_).gridlocalgeom[parentId]; }
85
86 template<int side>
88 { return *std::get<side>(sideData_).gridgeom; }
89
90 template<int side>
91 bool local() const
92 { return std::get<side>(sideData_).gridlocal; }
93
94 template<int side>
95 IndexType index(unsigned int parentId = 0) const
96 { return std::get<side>(sideData_).gridindices[parentId]; }
97
98 template<int side>
100 { return std::get<side>(sideData_).gridindices.size(); }
101
102 private:
103 template<int side>
104 void initializeGeometry(const GridGlue& glue, unsigned mergeindex);
105
106 /* M E M B E R V A R I A B L E S */
107
108 public:
111
112 private:
113 template<int side>
114 struct SideData {
116 bool gridlocal = false;
117
119 std::vector< GridIndexType<side> > gridindices;
120
122 std::vector< std::optional< GridLocalGeometry<side> > > gridlocalgeom;
123
131 std::optional< GridGeometry<side> > gridgeom;
132 };
133
134 std::tuple< SideData<0>, SideData<1> > sideData_;
135 };
136
137 template<typename P0, typename P1>
138 template<int side>
139 void IntersectionData<P0, P1>::initializeGeometry(const GridGlue& glue, unsigned mergeindex)
140 {
141 auto& data = std::get<side>(sideData_);
142
143 const unsigned n_parents = glue.merger_->template parents<side>(mergeindex);
144
145 // init containers
146 data.gridindices.resize(n_parents);
147 data.gridlocalgeom.resize(n_parents);
148
149 // default values
150 data.gridindices[0] = 0;
151
152 static constexpr int nSimplexCorners = mydim + 1;
153 using ctype = typename GridGlue::ctype;
154
155 // initialize the local and the global geometries of grid `side`
156
157 // compute the coordinates of the subface's corners in codim 0 entity local coordinates
158 static constexpr int elementdim = GridGlue::template GridView<side>::template Codim<0>::Geometry::mydimension;
159
160 // coordinates within the subentity that contains the remote intersection
161 std::array<Dune::FieldVector< ctype, dim<side>() >, nSimplexCorners> corners_subEntity_local;
162
163 for (unsigned int par = 0; par < n_parents; ++par) {
164 for (int i = 0; i < nSimplexCorners; ++i)
165 corners_subEntity_local[i] = glue.merger_->template parentLocal<side>(mergeindex, i, par);
166
167 // Coordinates of the remote intersection corners wrt the element coordinate system
168 std::array<Dune::FieldVector<ctype, elementdim>, nSimplexCorners> corners_element_local;
169
170 if (data.gridlocal) {
171 data.gridindices[par] = glue.merger_->template parent<side>(mergeindex,par);
172
173 typename GridGlue::template GridPatch<side>::LocalGeometry
174 gridLocalGeometry = glue.template patch<side>().geometryLocal(data.gridindices[par]);
175 for (std::size_t i=0; i<corners_subEntity_local.size(); i++) {
176 corners_element_local[i] = gridLocalGeometry.global(corners_subEntity_local[i]);
177 }
178
179 // set the corners of the local geometry
180#ifdef ONLY_SIMPLEX_INTERSECTIONS
181# if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
182 const Dune::GeometryType type = Dune::GeometryTypes::simplex(mydim);
183# else
184 const Dune::GeometryType type(Dune::GeometryType::simplex, mydim);
185# endif
186#else
187#error Not Implemented
188#endif
189 data.gridlocalgeom[par].emplace(type, corners_element_local);
190
191 // Add world geometry only for 0th parent
192 if (par == 0) {
193 typename GridGlue::template GridPatch<side>::Geometry
194 gridWorldGeometry = glue.template patch<side>().geometry(data.gridindices[par]);
195
196 // world coordinates of the remote intersection corners
197 std::array<Dune::FieldVector<ctype, GridGlue::template GridView<side>::dimensionworld>, nSimplexCorners> corners_global;
198
199 for (std::size_t i=0; i<corners_subEntity_local.size(); i++) {
200 corners_global[i] = gridWorldGeometry.global(corners_subEntity_local[i]);
201 }
202
203 data.gridgeom.emplace(type, corners_global);
204 }
205 }
206 }
207 }
208
210 template<typename P0, typename P1>
211 IntersectionData<P0, P1>::IntersectionData(const GridGlue& glue, unsigned int mergeindex, unsigned int offset,
212 bool grid0local, bool grid1local)
213 : index_(mergeindex+offset)
214 {
215 // if an invalid index is given do not proceed!
216 // (happens when the parent GridGlue initializes the "end"-Intersection)
217 assert (0 <= mergeindex || mergeindex < glue.index__sz);
218
219 std::get<0>(sideData_).gridlocal = grid0local;
220 std::get<1>(sideData_).gridlocal = grid1local;
221
222 initializeGeometry<0>(glue, mergeindex);
223 initializeGeometry<1>(glue, mergeindex);
224 }
225
230 template<typename P0, typename P1, int inside, int outside>
232 {
235
236 using InsideGridView = typename GridGlue::template GridView<inside>;
237 using OutsideGridView = typename GridGlue::template GridView<outside>;
238
239 using InsideLocalGeometry = typename IntersectionData::template GridLocalGeometry<inside>;
240 using OutsideLocalGeometry = typename IntersectionData::template GridLocalGeometry<outside>;
241
242 using Geometry = typename IntersectionData::template GridGeometry<inside>;
243 using OutsideGeometry = typename IntersectionData::template GridGeometry<outside>;
244
245 static constexpr auto coorddim = IntersectionData::coorddim;
246 static constexpr auto mydim = IntersectionData::mydim;
247 static constexpr int insidePatch = inside;
248 static constexpr int outsidePatch = outside;
249
250 using ctype = typename GridGlue::ctype;
251 using LocalCoordinate = Dune::FieldVector<ctype, mydim>;
252 using GlobalCoordinate = Dune::FieldVector<ctype, coorddim>;
253 };
254
257 template<typename P0, typename P1, int I, int O>
259 {
260
261 public:
262
264
265 typedef typename Traits::GridGlue GridGlue;
267
268
271
275
276 typedef typename Traits::Geometry Geometry;
277 typedef typename Traits::ctype ctype;
278
279 typedef typename InsideGridView::Traits::template Codim<0>::Entity InsideEntity;
280 typedef typename OutsideGridView::Traits::template Codim<0>::Entity OutsideEntity;
281
284
286 static constexpr auto coorddim = Traits::coorddim;
287
289 static constexpr auto mydim = Traits::mydim;
290
292 static constexpr int insidePatch = Traits::insidePatch;
293
295 static constexpr int outsidePatch = Traits::outsidePatch;
296
297 // typedef unsigned int IndexType;
298
299 private:
303 static constexpr int codimensionWorld = coorddim - mydim;
304
305 public:
306 /* C O N S T R U C T O R S */
307
309 Intersection(const GridGlue* glue, const IntersectionData* i) :
310 glue_(glue), i_(i) {}
311
312 /* F U N C T I O N A L I T Y */
313
317 inside(unsigned int parentId = 0) const
318 {
319 assert(self());
320 return glue_->template patch<I>().element(i_->template index<I>(parentId));
321 }
322
326 outside(unsigned int parentId = 0) const
327 {
328 assert(neighbor());
329 return glue_->template patch<O>().element(i_->template index<O>(parentId));
330 }
331
333 bool conforming() const
334 {
335 throw Dune::NotImplemented();
336 }
337
340 const InsideLocalGeometry& geometryInInside(unsigned int parentId = 0) const
341 {
342 return i_->template localGeometry<I>(parentId);
343 }
344
347 const OutsideLocalGeometry& geometryInOutside(unsigned int parentId = 0) const
348 {
349 return i_->template localGeometry<O>(parentId);
350 }
351
358 const Geometry& geometry() const
359 {
360 return i_->template geometry<I>();
361 }
362
369 const OutsideGeometry& geometryOutside() const // DUNE_DEPRECATED
370 {
371 return i_->template geometry<O>();
372 }
373
375 Dune::GeometryType type() const
376 {
377 #ifdef ONLY_SIMPLEX_INTERSECTIONS
378 # if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
379 return Dune::GeometryTypes::simplex(mydim);
380 # else
381 static const Dune::GeometryType type(Dune::GeometryType::simplex, mydim);
382 return type;
383 # endif
384 #else
385 #error Not Implemented
386 #endif
387 }
388
389
391 bool self() const
392 {
393 return i_->template local<I>();
394 }
395
397 size_t neighbor(unsigned int g = 0) const
398 {
399 if (g == 0 && i_->template local<O>()) {
400 return i_->template parents<O>();
401 } else if (g == 1 && i_->template local<I>()) {
402 return i_->template parents<I>();
403 }
404 return 0;
405 }
406
408 int indexInInside(unsigned int parentId = 0) const
409 {
410 assert(self());
411 return glue_->template patch<I>().indexInInside(i_->template index<I>(parentId));
412 }
413
415 int indexInOutside(unsigned int parentId = 0) const
416 {
417 assert(neighbor());
418 return glue_->template patch<O>().indexInInside(i_->template index<O>(parentId));
419 }
420
426 {
427 GlobalCoordinate normal;
428
429 if (codimensionWorld == 0)
430 DUNE_THROW(Dune::Exception, "There is no normal vector to a full-dimensional intersection");
431 else if (codimensionWorld == 1) {
432 /* TODO: Implement the general n-ary cross product here */
433 const auto jacobianTransposed = geometry().jacobianTransposed(local);
434 if (mydim==1) {
435 normal[0] = - jacobianTransposed[0][1];
436 normal[1] = jacobianTransposed[0][0];
437 } else if (mydim==2) {
438 normal[0] = (jacobianTransposed[0][1] * jacobianTransposed[1][2] - jacobianTransposed[0][2] * jacobianTransposed[1][1]);
439 normal[1] = - (jacobianTransposed[0][0] * jacobianTransposed[1][2] - jacobianTransposed[0][2] * jacobianTransposed[1][0]);
440 normal[2] = (jacobianTransposed[0][0] * jacobianTransposed[1][1] - jacobianTransposed[0][1] * jacobianTransposed[1][0]);
441 } else
442 DUNE_THROW(Dune::NotImplemented, "Remote intersections don't implement the 'outerNormal' method for " << mydim << "-dimensional intersections yet");
443 } else
444 DUNE_THROW(Dune::NotImplemented, "Remote intersections don't implement the 'outerNormal' method for intersections with codim >= 2 yet");
445
446 return normal;
447 }
448
454 {
455 GlobalCoordinate normal = outerNormal(local);
456 normal /= normal.two_norm();
457 return normal;
458 }
459
465 {
466 return (unitOuterNormal(local) *= geometry().integrationElement(local));
467 }
468
474 {
475 return unitOuterNormal(ReferenceElements<ctype,mydim>::general(type()).position(0,0));
476 }
477
482 {
483 return Intersection<P0,P1,O,I>(glue_,i_);
484 }
485
486#ifdef QUICKHACK_INDEX
487 typedef typename GridGlue::IndexType IndexType;
488
489 IndexType index() const
490 {
491 return i_->index_;
492 }
493
494#endif
495
496 private:
497
498 friend class IntersectionIndexSet<P0,P1>;
499
500 /* M E M B E R V A R I A B L E S */
501
503 const GridGlue* glue_;
504
506 const IntersectionData* i_;
507 };
508
509
510 } // end namespace GridGlue
511} // end namespace Dune
512
513#endif // DUNE_GRIDGLUE_ADAPTER_INTERSECTION_HH
Central component of the module implementing the coupling of two grids.
Definition: gridglue.hh:35
sequential adapter to couple two grids at specified close together boundaries
Definition: gridglue.hh:65
unsigned int IndexType
Definition: gridglue.hh:145
PromotionTraits< typenameGridView< 0 >::ctype, typenameGridView< 1 >::ctype >::PromotedType ctype
The type used for coordinates.
Definition: gridglue.hh:169
static constexpr int dimworld
export the world dimension This is the maximum of the extractors' world dimensions.
Definition: gridglue.hh:164
storage class for Dune::GridGlue::Intersection related data
Definition: intersection.hh:36
GridGlue::IndexType IndexType
Definition: intersection.hh:40
AffineGeometry< typename GridGlue::template GridView< side >::ctype, mydim, GridGlue::template GridView< side >::dimensionworld > GridGeometry
Definition: intersection.hh:63
const GridLocalGeometry< side > & localGeometry(unsigned int parentId=0) const
Definition: intersection.hh:83
static constexpr int coorddim
Dimension of the world space of the intersection.
Definition: intersection.hh:43
GridGeometry< 0 > Grid0Geometry
Definition: intersection.hh:65
typename GridGlue::template GridView< side >::IndexSet::IndexType GridIndexType
Definition: intersection.hh:69
IndexType index(unsigned int parentId=0) const
Definition: intersection.hh:95
IndexType parents() const
Definition: intersection.hh:99
GridLocalGeometry< 1 > Grid1LocalGeometry
Definition: intersection.hh:59
const GridGeometry< side > & geometry() const
Definition: intersection.hh:87
bool local() const
Definition: intersection.hh:91
::Dune::GridGlue::GridGlue< P0, P1 > GridGlue
Definition: intersection.hh:38
static constexpr int mydim
Dimension of the intersection.
Definition: intersection.hh:52
AffineGeometry< typename GridGlue::template GridView< side >::ctype, mydim, GridGlue::template GridView< side >::dimension > GridLocalGeometry
Definition: intersection.hh:56
GridLocalGeometry< 0 > Grid0LocalGeometry
Definition: intersection.hh:58
IndexType index_
index of this intersection after GridGlue interface
Definition: intersection.hh:110
GridIndexType< 1 > Grid1IndexType
Definition: intersection.hh:72
GridGeometry< 1 > Grid1Geometry
Definition: intersection.hh:66
GridIndexType< 0 > Grid0IndexType
Definition: intersection.hh:71
IntersectionData()=default
Default Constructor.
The intersection of two entities of the two patches of a GridGlue.
Definition: intersection.hh:259
Traits::OutsideLocalGeometry OutsideLocalGeometry
Definition: intersection.hh:273
Intersection< P0, P1, O, I > flip() const
Return a copy of the intersection with inside and outside switched.
Definition: intersection.hh:481
bool conforming() const
Return true if intersection is conforming.
Definition: intersection.hh:333
InsideGridView::Traits::template Codim< 0 >::Entity InsideEntity
Definition: intersection.hh:279
Traits::GridGlue GridGlue
Definition: intersection.hh:265
int indexInOutside(unsigned int parentId=0) const
Local number of codim 1 entity in outside() Entity where intersection is contained in.
Definition: intersection.hh:415
int indexInInside(unsigned int parentId=0) const
Local number of codim 1 entity in the inside() Entity where intersection is contained in.
Definition: intersection.hh:408
bool self() const
For parallel computations: Return true if inside() entity exists locally.
Definition: intersection.hh:391
static constexpr auto mydim
dimension of the intersection
Definition: intersection.hh:289
static constexpr int outsidePatch
outside patch
Definition: intersection.hh:295
Dune::GeometryType type() const
Type of reference element for this intersection.
Definition: intersection.hh:375
InsideEntity inside(unsigned int parentId=0) const
Return element on the inside of this intersection.
Definition: intersection.hh:317
Traits::LocalCoordinate LocalCoordinate
Definition: intersection.hh:282
const OutsideLocalGeometry & geometryInOutside(unsigned int parentId=0) const
Geometric information about this intersection in local coordinates of the outside() element.
Definition: intersection.hh:347
GlobalCoordinate unitOuterNormal(const LocalCoordinate &local) const
Return a unit outer normal.
Definition: intersection.hh:453
Traits::ctype ctype
Definition: intersection.hh:277
IntersectionTraits< P0, P1, I, O > Traits
Definition: intersection.hh:263
Traits::IntersectionData IntersectionData
Definition: intersection.hh:266
Traits::GlobalCoordinate GlobalCoordinate
Definition: intersection.hh:283
Traits::OutsideGeometry OutsideGeometry
Definition: intersection.hh:274
size_t neighbor(unsigned int g=0) const
Return number of embeddings into local grid0 (grid1) entities.
Definition: intersection.hh:397
Intersection(const GridGlue *glue, const IntersectionData *i)
Constructor for a given Dataset.
Definition: intersection.hh:309
static constexpr int insidePatch
inside patch
Definition: intersection.hh:292
Traits::Geometry Geometry
Definition: intersection.hh:276
OutsideGridView::Traits::template Codim< 0 >::Entity OutsideEntity
Definition: intersection.hh:280
OutsideEntity outside(unsigned int parentId=0) const
Return element on the outside of this intersection.
Definition: intersection.hh:326
const Geometry & geometry() const
Geometric information about this intersection as part of the inside grid.
Definition: intersection.hh:358
Traits::InsideLocalGeometry InsideLocalGeometry
Definition: intersection.hh:270
GlobalCoordinate outerNormal(const LocalCoordinate &local) const
Return an outer normal (length not necessarily 1)
Definition: intersection.hh:425
const OutsideGeometry & geometryOutside() const
Geometric information about this intersection as part of the outside grid.
Definition: intersection.hh:369
Traits::OutsideGridView OutsideGridView
Definition: intersection.hh:272
GlobalCoordinate integrationOuterNormal(const LocalCoordinate &local) const
Return an outer normal with the length of the integration element.
Definition: intersection.hh:464
Traits::InsideGridView InsideGridView
Definition: intersection.hh:269
const InsideLocalGeometry & geometryInInside(unsigned int parentId=0) const
Geometric information about this intersection in local coordinates of the inside() element.
Definition: intersection.hh:340
GlobalCoordinate centerUnitOuterNormal() const
Unit outer normal at the center of the intersection.
Definition: intersection.hh:473
static constexpr auto coorddim
dimension of the world space of the intersection
Definition: intersection.hh:286
Definition: intersection.hh:232
static constexpr int insidePatch
Definition: intersection.hh:247
Dune::FieldVector< ctype, mydim > LocalCoordinate
Definition: intersection.hh:251
static constexpr auto coorddim
Definition: intersection.hh:245
typename IntersectionData::template GridGeometry< outside > OutsideGeometry
Definition: intersection.hh:243
typename GridGlue::ctype ctype
Definition: intersection.hh:250
typename GridGlue::template GridView< outside > OutsideGridView
Definition: intersection.hh:237
Dune::FieldVector< ctype, coorddim > GlobalCoordinate
Definition: intersection.hh:252
typename IntersectionData::template GridLocalGeometry< outside > OutsideLocalGeometry
Definition: intersection.hh:240
typename GridGlue::template GridView< inside > InsideGridView
Definition: intersection.hh:236
typename IntersectionData::template GridGeometry< inside > Geometry
Definition: intersection.hh:242
typename IntersectionData::template GridLocalGeometry< inside > InsideLocalGeometry
Definition: intersection.hh:239
static constexpr int outsidePatch
Definition: intersection.hh:248
static constexpr auto mydim
Definition: intersection.hh:246