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