dune-pdelab  2.0.0
intersectionindexset.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 #ifndef DUNE_PDELAB_INTERSECTIONINDEXSET_HH
4 #define DUNE_PDELAB_INTERSECTIONINDEXSET_HH
5 
6 //===============================================================
7 // Index set for intersections
8 // implementation only for
9 // - hanging node grids, but conforming macro grid
10 // - no multiple element types !
11 //===============================================================
12 
13 #include<vector>
14 #include <iostream>
15 
16 #include<dune/common/exceptions.hh>
17 
18 namespace Dune {
19  namespace PDELab {
20 
21  template<typename GV>
23  {
24  public:
25  typedef typename GV::IndexSet IndexSet;
26  typedef typename IndexSet::IndexType IndexType;
27  typedef typename GV::Intersection Intersection;
28  typedef typename GV::Traits::template Codim<0>::Entity Element;
29 
30  IntersectionIndexSet (const GV& gv_)
31  : gv(gv_), is(gv.indexSet())
32  {
33  typedef typename GV::Traits::template Codim<0>::Iterator ElementIterator;
34 
35  pre();
36  for (ElementIterator it = gv.template begin<0>();
37  it!=gv.template end<0>(); ++it)
38  {
39  visit(*it);
40  }
41  post();
42  }
43 
44  // number of intersections in index set
45  // (intersections do not have a geometry type)
46  IndexType size () const
47  {
48  return intersection_counter;
49  }
50 
51  // number of intersections associated with given element
52  IndexType size (const Element& element) const
53  {
54  return number_of_intersections[is.index(element)];
55  }
56 
57  // get index assigned to intersection
58  IndexType index (const Intersection& intersection) const
59  {
60  // on the boundary, take own codim 1 subentity
61  if (intersection.boundary() && (!intersection.neighbor()))
62  return codim1index_to_intersectionindex[is.subIndex(*(intersection.inside()),intersection.indexInInside(),1)];
63 
64  // if we have a neighbor, take higher level
65  if (intersection.neighbor())
66  {
67  if (intersection.inside()->level()>=intersection.outside()->level())
68  return codim1index_to_intersectionindex[is.subIndex(*(intersection.inside()),intersection.indexInInside(),1)];
69  else
70  return codim1index_to_intersectionindex[is.subIndex(*(intersection.outside()),intersection.indexInOutside(),1)];
71  }
72 
73  // we are at a processor boundary
74  DUNE_THROW(Dune::Exception,"intersection index at processor boundary requested");
75  }
76 
77  // get index of i'th intersection of element
78  // (in order they are visited by intersection iterator)
79  IndexType subIndex (const Element& element, int i) const
80  {
81  return element_intersection_subindex[entry[is.index(element)]+i];
82  }
83 
84  private:
85 
86  // prepare loop over elements
87  void pre ()
88  {
89  codim1index_to_intersectionindex.resize(gv.size(1));
90  invalidIndex = gv.size(1)+1;
91  for (size_t i=0; i<codim1index_to_intersectionindex.size(); ++i)
92  codim1index_to_intersectionindex[i] = invalidIndex;
93  number_of_intersections.resize(gv.size(0));
94  entry.resize(gv.size(0));
95  element_intersection_subindex.resize(2*gv.size(1));
96  intersection_counter = 0;
97  oriented_intersection_counter = 0;
98  std::cout << "number of codim 1 entities is " << gv.size(1) << std::endl;
99  }
100 
101  // process given element
102  void visit (const Element& element)
103  {
104  typedef typename GV::IntersectionIterator IntersectionIterator;
105  size_t count = 0;
106  entry[is.index(element)] = oriented_intersection_counter;
107  IntersectionIterator endit = gv.iend(element);
108  for (IntersectionIterator iit = gv.ibegin(element); iit!=endit; ++iit)
109  {
110  if (iit->neighbor())
111  {
112  IndexType c1index;
113  if (iit->inside()->level()>=iit->outside()->level())
114  c1index = is.subIndex(*(iit->inside()),iit->indexInInside(),1);
115  else
116  c1index = is.subIndex(*(iit->outside()),iit->indexInOutside(),1);
117  if (codim1index_to_intersectionindex[c1index]==invalidIndex)
118  codim1index_to_intersectionindex[c1index]=intersection_counter++;
119  element_intersection_subindex[oriented_intersection_counter] = codim1index_to_intersectionindex[c1index];
120  }
121  else if (iit->boundary())
122  {
123  IndexType c1index = is.subIndex(*(iit->inside()),iit->indexInInside(),1);
124  if (codim1index_to_intersectionindex[c1index]==invalidIndex)
125  codim1index_to_intersectionindex[c1index]=intersection_counter++;
126  element_intersection_subindex[oriented_intersection_counter] = codim1index_to_intersectionindex[c1index];
127  }
128  count++;
129  oriented_intersection_counter++;
130  }
131  number_of_intersections[is.index(element)] = static_cast<unsigned char>(count);
132  }
133 
134  // finalize computation of index set
135  void post ()
136  {
137  std::cout << "number of oriented intersections " << oriented_intersection_counter << std::endl;
138  std::cout << "number of intersections " << intersection_counter << std::endl;
139  }
140 
141  GV gv; // our grid view
142  const IndexSet& is; // index set of the grid view
143 
144  // we assume that the mesh is conforming in the
145  // sense that there is a codim 1 entity for each intersection
146  // the following vector assigns intersection to the codim 1 entity on the "higher" side
147  std::vector<IndexType> codim1index_to_intersectionindex;
148 
149  // number of intersections of an element
150  std::vector<unsigned char> number_of_intersections;
151 
152  // map (element index,interection number) to index
153  std::vector<IndexType> element_intersection_subindex;
154 
155  // entry point of element into element_intersection_subindex
156  std::vector<size_t> entry;
157 
158  size_t intersection_counter;
159  size_t oriented_intersection_counter;
160  IndexType invalidIndex;
161  };
162  }
163 }
164 
165 #endif
IndexType subIndex(const Element &element, int i) const
Definition: intersectionindexset.hh:79
IndexSet::IndexType IndexType
Definition: intersectionindexset.hh:26
IndexType size(const Element &element) const
Definition: intersectionindexset.hh:52
IndexType size() const
Definition: intersectionindexset.hh:46
GV::Traits::template Codim< 0 >::Entity Element
Definition: intersectionindexset.hh:28
Definition: intersectionindexset.hh:22
IndexType index(const Intersection &intersection) const
Definition: intersectionindexset.hh:58
IntersectionIndexSet(const GV &gv_)
Definition: intersectionindexset.hh:30
GV::IndexSet IndexSet
Definition: intersectionindexset.hh:25
GV::Intersection Intersection
Definition: intersectionindexset.hh:27