dune-grid-glue  2.7.0
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 <memory>
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 
24 namespace 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 
40  typedef typename GridGlue::IndexType IndexType;
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 DUNE_DEPRECATED_MSG("please use GridLocalGeometry<0> instead") = GridLocalGeometry<0>;
59  using Grid1LocalGeometry DUNE_DEPRECATED_MSG("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 DUNE_DEPRECATED_MSG("please use GridGeometry<0> instead") = GridGeometry<0>;
66  using Grid1Geometry DUNE_DEPRECATED_MSG("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 DUNE_DEPRECATED_MSG("please use GridIndexType<0> instead") = GridIndexType<0>;
72  using Grid1IndexType DUNE_DEPRECATED_MSG("please use GridIndexType<1> instead") = GridIndexType<1>;
73 
75  IntersectionData(const GridGlue& glue, unsigned int mergeindex, unsigned int offset, bool grid0local, bool grid1local);
76 
79  {
80  std::get<0>(sideData_).gridlocal = false;
81  std::get<1>(sideData_).gridlocal = false;
82  }
83 
84  /* Accessor Functions */
85 
86  template<int side>
87  const GridLocalGeometry<side>& localGeometry(unsigned int parentId = 0) const
88  { return *std::get<side>(sideData_).gridlocalgeom[parentId]; }
89 
90  template<int side>
92  { return *std::get<side>(sideData_).gridgeom; }
93 
94  template<int side>
95  bool local() const
96  { return std::get<side>(sideData_).gridlocal; }
97 
98  template<int side>
99  IndexType index(unsigned int parentId = 0) const
100  { return std::get<side>(sideData_).gridindices[parentId]; }
101 
102  template<int side>
104  { return std::get<side>(sideData_).gridindices.size(); }
105 
106  private:
107  template<int side>
108  void initializeGeometry(const GridGlue& glue, unsigned mergeindex);
109 
110  /* M E M B E R V A R I A B L E S */
111 
112  public:
115 
116  private:
117  template<int side>
118  struct SideData {
120  bool gridlocal;
121 
123  std::vector< GridIndexType<side> > gridindices;
124 
126  /* TODO [C++17 or DUNE-2.6]: use std::optional */
127  std::vector< std::unique_ptr< GridLocalGeometry<side> > > gridlocalgeom;
128 
136  /* TODO [C++17 or DUNE-2.6]: use std::optional */
137  std::unique_ptr< GridGeometry<side> > gridgeom;
138  };
139 
140  std::tuple< SideData<0>, SideData<1> > sideData_;
141  };
142 
143  template<typename P0, typename P1>
144  template<int side>
145  void IntersectionData<P0, P1>::initializeGeometry(const GridGlue& glue, unsigned mergeindex)
146  {
147  auto& data = std::get<side>(sideData_);
148 
149  const unsigned n_parents = glue.merger_->template parents<side>(mergeindex);
150 
151  // init containers
152  data.gridindices.resize(n_parents);
153  data.gridlocalgeom.resize(n_parents);
154 
155  // default values
156  data.gridindices[0] = 0;
157 
158  static constexpr int nSimplexCorners = mydim + 1;
159  using ctype = typename GridGlue::ctype;
160 
161  // initialize the local and the global geometries of grid `side`
162 
163  // compute the coordinates of the subface's corners in codim 0 entity local coordinates
164  static constexpr int elementdim = GridGlue::template GridView<side>::template Codim<0>::Geometry::mydimension;
165 
166  // coordinates within the subentity that contains the remote intersection
167  std::array<Dune::FieldVector< ctype, dim<side>() >, nSimplexCorners> corners_subEntity_local;
168 
169  for (unsigned int par = 0; par < n_parents; ++par) {
170  for (int i = 0; i < nSimplexCorners; ++i)
171  corners_subEntity_local[i] = glue.merger_->template parentLocal<side>(mergeindex, i, par);
172 
173  // Coordinates of the remote intersection corners wrt the element coordinate system
174  std::array<Dune::FieldVector<ctype, elementdim>, nSimplexCorners> corners_element_local;
175 
176  if (data.gridlocal) {
177  data.gridindices[par] = glue.merger_->template parent<side>(mergeindex,par);
178 
179  typename GridGlue::template GridPatch<side>::LocalGeometry
180  gridLocalGeometry = glue.template patch<side>().geometryLocal(data.gridindices[par]);
181  for (std::size_t i=0; i<corners_subEntity_local.size(); i++) {
182  corners_element_local[i] = gridLocalGeometry.global(corners_subEntity_local[i]);
183  }
184 
185  // set the corners of the local geometry
186 #ifdef ONLY_SIMPLEX_INTERSECTIONS
187 # if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
188  const Dune::GeometryType type = Dune::GeometryTypes::simplex(mydim);
189 # else
190  const Dune::GeometryType type(Dune::GeometryType::simplex, mydim);
191 # endif
192 #else
193 #error Not Implemented
194 #endif
195  data.gridlocalgeom[par] = std::make_unique< GridLocalGeometry<side> >(type, corners_element_local);
196 
197  // Add world geometry only for 0th parent
198  if (par == 0) {
199  typename GridGlue::template GridPatch<side>::Geometry
200  gridWorldGeometry = glue.template patch<side>().geometry(data.gridindices[par]);
201 
202  // world coordinates of the remote intersection corners
203  std::array<Dune::FieldVector<ctype, GridGlue::template GridView<side>::dimensionworld>, nSimplexCorners> corners_global;
204 
205  for (std::size_t i=0; i<corners_subEntity_local.size(); i++) {
206  corners_global[i] = gridWorldGeometry.global(corners_subEntity_local[i]);
207  }
208 
209  data.gridgeom = std::make_unique< GridGeometry<side> >(type, corners_global);
210  }
211  }
212  }
213  }
214 
216  template<typename P0, typename P1>
217  IntersectionData<P0, P1>::IntersectionData(const GridGlue& glue, unsigned int mergeindex, unsigned int offset,
218  bool grid0local, bool grid1local)
219  : index_(mergeindex+offset)
220  {
221  // if an invalid index is given do not proceed!
222  // (happens when the parent GridGlue initializes the "end"-Intersection)
223  assert (0 <= mergeindex || mergeindex < glue.index__sz);
224 
225  std::get<0>(sideData_).gridlocal = grid0local;
226  std::get<1>(sideData_).gridlocal = grid1local;
227 
228  initializeGeometry<0>(glue, mergeindex);
229  initializeGeometry<1>(glue, mergeindex);
230  }
231 
236  template<typename P0, typename P1, int inside, int outside>
238  {
241 
242  using InsideGridView = typename GridGlue::template GridView<inside>;
243  using OutsideGridView = typename GridGlue::template GridView<outside>;
244 
245  using InsideLocalGeometry = typename IntersectionData::template GridLocalGeometry<inside>;
246  using OutsideLocalGeometry = typename IntersectionData::template GridLocalGeometry<outside>;
247 
248  using Geometry = typename IntersectionData::template GridGeometry<inside>;
249  using OutsideGeometry = typename IntersectionData::template GridGeometry<outside>;
250 
251  static constexpr auto coorddim = IntersectionData::coorddim;
252  static constexpr auto mydim = IntersectionData::mydim;
253  static constexpr int insidePatch = inside;
254  static constexpr int outsidePatch = outside;
255 
256  using ctype = typename GridGlue::ctype;
257  using LocalCoordinate = Dune::FieldVector<ctype, mydim>;
258  using GlobalCoordinate = Dune::FieldVector<ctype, coorddim>;
259  };
260 
263  template<typename P0, typename P1, int I, int O>
265  {
266 
267  public:
268 
270 
271  typedef typename Traits::GridGlue GridGlue;
273 
274 
277 
281 
282  typedef typename Traits::Geometry Geometry;
283  typedef typename Traits::ctype ctype;
284 
285  typedef typename InsideGridView::Traits::template Codim<0>::Entity InsideEntity;
286  typedef typename OutsideGridView::Traits::template Codim<0>::Entity OutsideEntity;
287 
290 
292  static constexpr auto coorddim = Traits::coorddim;
293 
295  static constexpr auto mydim = Traits::mydim;
296 
298  static constexpr int insidePatch = Traits::insidePatch;
299 
301  static constexpr int outsidePatch = Traits::outsidePatch;
302 
303  // typedef unsigned int IndexType;
304 
305  private:
309  static constexpr int codimensionWorld = coorddim - mydim;
310 
311  public:
312  /* C O N S T R U C T O R S */
313 
315  Intersection(const GridGlue* glue, const IntersectionData* i) :
316  glue_(glue), i_(i) {}
317 
318  /* F U N C T I O N A L I T Y */
319 
323  inside(unsigned int parentId = 0) const
324  {
325  assert(self());
326  return glue_->template patch<I>().element(i_->template index<I>(parentId));
327  }
328 
332  outside(unsigned int parentId = 0) const
333  {
334  assert(neighbor());
335  return glue_->template patch<O>().element(i_->template index<O>(parentId));
336  }
337 
339  bool conforming() const
340  {
341  throw Dune::NotImplemented();
342  }
343 
346  const InsideLocalGeometry& geometryInInside(unsigned int parentId = 0) const
347  {
348  return i_->template localGeometry<I>(parentId);
349  }
350 
353  const OutsideLocalGeometry& geometryInOutside(unsigned int parentId = 0) const
354  {
355  return i_->template localGeometry<O>(parentId);
356  }
357 
364  const Geometry& geometry() const
365  {
366  return i_->template geometry<I>();
367  }
368 
375  const OutsideGeometry& geometryOutside() const // DUNE_DEPRECATED
376  {
377  return i_->template geometry<O>();
378  }
379 
381  Dune::GeometryType type() const
382  {
383  #ifdef ONLY_SIMPLEX_INTERSECTIONS
384  # if DUNE_VERSION_NEWER(DUNE_GEOMETRY, 2, 6)
385  return Dune::GeometryTypes::simplex(mydim);
386  # else
387  static const Dune::GeometryType type(Dune::GeometryType::simplex, mydim);
388  return type;
389  # endif
390  #else
391  #error Not Implemented
392  #endif
393  }
394 
395 
397  bool self() const
398  {
399  return i_->template local<I>();
400  }
401 
403  size_t neighbor(unsigned int g = 0) const
404  {
405  if (g == 0 && i_->template local<O>()) {
406  return i_->template parents<O>();
407  } else if (g == 1 && i_->template local<I>()) {
408  return i_->template parents<I>();
409  }
410  return 0;
411  }
412 
414  int indexInInside(unsigned int parentId = 0) const
415  {
416  assert(self());
417  return glue_->template patch<I>().indexInInside(i_->template index<I>(parentId));
418  }
419 
421  int indexInOutside(unsigned int parentId = 0) const
422  {
423  assert(neighbor());
424  return glue_->template patch<O>().indexInInside(i_->template index<O>(parentId));
425  }
426 
432  {
433  GlobalCoordinate normal;
434 
435  if (codimensionWorld == 0)
436  DUNE_THROW(Dune::Exception, "There is no normal vector to a full-dimensional intersection");
437  else if (codimensionWorld == 1) {
438  /* TODO: Implement the general n-ary cross product here */
439  const auto jacobianTransposed = geometry().jacobianTransposed(local);
440  if (mydim==1) {
441  normal[0] = - jacobianTransposed[0][1];
442  normal[1] = jacobianTransposed[0][0];
443  } else if (mydim==2) {
444  normal[0] = (jacobianTransposed[0][1] * jacobianTransposed[1][2] - jacobianTransposed[0][2] * jacobianTransposed[1][1]);
445  normal[1] = - (jacobianTransposed[0][0] * jacobianTransposed[1][2] - jacobianTransposed[0][2] * jacobianTransposed[1][0]);
446  normal[2] = (jacobianTransposed[0][0] * jacobianTransposed[1][1] - jacobianTransposed[0][1] * jacobianTransposed[1][0]);
447  } else
448  DUNE_THROW(Dune::NotImplemented, "Remote intersections don't implement the 'outerNormal' method for " << mydim << "-dimensional intersections yet");
449  } else
450  DUNE_THROW(Dune::NotImplemented, "Remote intersections don't implement the 'outerNormal' method for intersections with codim >= 2 yet");
451 
452  return normal;
453  }
454 
460  {
461  GlobalCoordinate normal = outerNormal(local);
462  normal /= normal.two_norm();
463  return normal;
464  }
465 
471  {
472  return (unitOuterNormal(local) *= geometry().integrationElement(local));
473  }
474 
480  {
481  return unitOuterNormal(ReferenceElements<ctype,mydim>::general(type()).position(0,0));
482  }
483 
488  {
489  return Intersection<P0,P1,O,I>(glue_,i_);
490  }
491 
492 #ifdef QUICKHACK_INDEX
493  typedef typename GridGlue::IndexType IndexType;
494 
495  IndexType index() const
496  {
497  return i_->index_;
498  }
499 
500 #endif
501 
502  private:
503 
504  friend class IntersectionIndexSet<P0,P1>;
505 
506  /* M E M B E R V A R I A B L E S */
507 
509  const GridGlue* glue_;
510 
512  const IntersectionData* i_;
513  };
514 
515 
516  } // end namespace GridGlue
517 } // end namespace Dune
518 
519 #endif // DUNE_GRIDGLUE_ADAPTER_INTERSECTION_HH
Central component of the module implementing the coupling of two grids.
Definition: gridglue.hh:35
bool inside(const Coordinate &x, const Field &epsilon)
Definition: projection_impl.hh:109
sequential adapter to couple two grids at specified close together boundaries
Definition: gridglue.hh:52
unsigned int IndexType
Definition: gridglue.hh:145
static constexpr int dimworld
export the world dimension This is the maximum of the extractors' world dimensions.
Definition: gridglue.hh:164
PromotionTraits< typename GridView< 0 >::ctype, typename GridView< 1 >::ctype >::PromotedType ctype
The type used for coordinates.
Definition: gridglue.hh:169
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
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:99
IndexType parents() const
Definition: intersection.hh:103
GridLocalGeometry< 1 > Grid1LocalGeometry
Definition: intersection.hh:59
const GridGeometry< side > & geometry() const
Definition: intersection.hh:91
const GridLocalGeometry< side > & localGeometry(unsigned int parentId=0) const
Definition: intersection.hh:87
bool local() const
Definition: intersection.hh:95
::Dune::GridGlue::GridGlue< P0, P1 > GridGlue
Definition: intersection.hh:38
static constexpr int mydim
Dimension of the intersection.
Definition: intersection.hh:52
IntersectionData()
Default Constructor.
Definition: intersection.hh:78
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:114
GridIndexType< 1 > Grid1IndexType
Definition: intersection.hh:72
GridGeometry< 1 > Grid1Geometry
Definition: intersection.hh:66
GridIndexType< 0 > Grid0IndexType
Definition: intersection.hh:71
The intersection of two entities of the two patches of a GridGlue.
Definition: intersection.hh:265
Traits::OutsideLocalGeometry OutsideLocalGeometry
Definition: intersection.hh:279
bool conforming() const
Return true if intersection is conforming.
Definition: intersection.hh:339
InsideGridView::Traits::template Codim< 0 >::Entity InsideEntity
Definition: intersection.hh:285
Intersection< P0, P1, O, I > flip() const
Return a copy of the intersection with inside and outside switched.
Definition: intersection.hh:487
Traits::GridGlue GridGlue
Definition: intersection.hh:271
int indexInOutside(unsigned int parentId=0) const
Local number of codim 1 entity in outside() Entity where intersection is contained in.
Definition: intersection.hh:421
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:414
const OutsideGeometry & geometryOutside() const
Geometric information about this intersection as part of the outside grid.
Definition: intersection.hh:375
static constexpr auto mydim
dimension of the intersection
Definition: intersection.hh:295
static constexpr int outsidePatch
outside patch
Definition: intersection.hh:301
const InsideLocalGeometry & geometryInInside(unsigned int parentId=0) const
Geometric information about this intersection in local coordinates of the inside() element.
Definition: intersection.hh:346
Dune::GeometryType type() const
Type of reference element for this intersection.
Definition: intersection.hh:381
InsideEntity inside(unsigned int parentId=0) const
Return element on the inside of this intersection.
Definition: intersection.hh:323
Traits::LocalCoordinate LocalCoordinate
Definition: intersection.hh:288
GlobalCoordinate unitOuterNormal(const LocalCoordinate &local) const
Return a unit outer normal.
Definition: intersection.hh:459
const OutsideLocalGeometry & geometryInOutside(unsigned int parentId=0) const
Geometric information about this intersection in local coordinates of the outside() element.
Definition: intersection.hh:353
Traits::ctype ctype
Definition: intersection.hh:283
IntersectionTraits< P0, P1, I, O > Traits
Definition: intersection.hh:269
Traits::IntersectionData IntersectionData
Definition: intersection.hh:272
Traits::GlobalCoordinate GlobalCoordinate
Definition: intersection.hh:289
Traits::OutsideGeometry OutsideGeometry
Definition: intersection.hh:280
size_t neighbor(unsigned int g=0) const
Return number of embeddings into local grid0 (grid1) entities.
Definition: intersection.hh:403
Intersection(const GridGlue *glue, const IntersectionData *i)
Constructor for a given Dataset.
Definition: intersection.hh:315
static constexpr int insidePatch
inside patch
Definition: intersection.hh:298
const Geometry & geometry() const
Geometric information about this intersection as part of the inside grid.
Definition: intersection.hh:364
Traits::Geometry Geometry
Definition: intersection.hh:282
OutsideGridView::Traits::template Codim< 0 >::Entity OutsideEntity
Definition: intersection.hh:286
OutsideEntity outside(unsigned int parentId=0) const
Return element on the outside of this intersection.
Definition: intersection.hh:332
Traits::InsideLocalGeometry InsideLocalGeometry
Definition: intersection.hh:276
GlobalCoordinate outerNormal(const LocalCoordinate &local) const
Return an outer normal (length not necessarily 1)
Definition: intersection.hh:431
Traits::OutsideGridView OutsideGridView
Definition: intersection.hh:278
GlobalCoordinate integrationOuterNormal(const LocalCoordinate &local) const
Return an outer normal with the length of the integration element.
Definition: intersection.hh:470
Traits::InsideGridView InsideGridView
Definition: intersection.hh:275
GlobalCoordinate centerUnitOuterNormal() const
Unit outer normal at the center of the intersection.
Definition: intersection.hh:479
static constexpr auto coorddim
dimension of the world space of the intersection
Definition: intersection.hh:292
Definition: intersection.hh:238
static constexpr int insidePatch
Definition: intersection.hh:253
Dune::FieldVector< ctype, mydim > LocalCoordinate
Definition: intersection.hh:257
static constexpr auto coorddim
Definition: intersection.hh:251
typename IntersectionData::template GridGeometry< outside > OutsideGeometry
Definition: intersection.hh:249
typename GridGlue::ctype ctype
Definition: intersection.hh:256
typename GridGlue::template GridView< outside > OutsideGridView
Definition: intersection.hh:243
Dune::FieldVector< ctype, coorddim > GlobalCoordinate
Definition: intersection.hh:258
typename IntersectionData::template GridLocalGeometry< outside > OutsideLocalGeometry
Definition: intersection.hh:246
typename GridGlue::template GridView< inside > InsideGridView
Definition: intersection.hh:242
typename IntersectionData::template GridGeometry< inside > Geometry
Definition: intersection.hh:248
typename IntersectionData::template GridLocalGeometry< inside > InsideLocalGeometry
Definition: intersection.hh:245
static constexpr int outsidePatch
Definition: intersection.hh:254
static constexpr auto mydim
Definition: intersection.hh:252