dune-localfunctions  2.2.1
whitney/edges0.5/basis.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set ts=8 sw=2 et sts=2:
3 
4 #ifndef DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
5 #define DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
6 
7 #include <cstddef>
8 #include <vector>
9 
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/fvector.hh>
12 
16 
17 namespace Dune {
18 
20  //
21  // Basis
22  //
23 
25 
33  template<class Geometry, class RF>
34  class EdgeS0_5Basis :
35  private EdgeS0_5Common<Geometry::mydimension, typename Geometry::ctype>
36  {
37  public:
39  struct Traits {
40  typedef typename Geometry::ctype DomainField;
41  static const std::size_t dimDomainLocal = Geometry::mydimension;
42  static const std::size_t dimDomainGlobal = Geometry::coorddimension;
43  typedef FieldVector<DomainField, dimDomainLocal> DomainLocal;
44  typedef FieldVector<DomainField, dimDomainGlobal> DomainGlobal;
45 
46  typedef RF RangeField;
47  static const std::size_t dimRange = dimDomainLocal;
48  typedef FieldVector<RangeField, dimRange> Range;
49 
50  typedef FieldMatrix<RangeField, dimRange, dimDomainGlobal> Jacobian;
51 
52  static const std::size_t diffOrder = 1;
53  };
54 
55  private:
56  typedef Dune::P1LocalBasis<typename Traits::DomainField,
57  typename Traits::RangeField,
59  > P1LocalBasis;
61 
62  static const P1LocalBasis& p1LocalBasis;
63  static const std::size_t dim = Traits::dimDomainLocal;
64 
66  using Base::refelem;
67  using Base::s;
68 
69  // global values of the Jacobians (gradients) of the p1 basis
70  std::vector<typename P1Basis::Traits::Jacobian> p1j;
71  // edge sizes and orientations
72  std::vector<typename Traits::DomainField> edgel;
73 
74  public:
76 
82  template<typename VertexOrder>
83  EdgeS0_5Basis(const Geometry& geo, const VertexOrder& vertexOrder) :
84  p1j(s, typename P1Basis::Traits::Jacobian(0)), edgel(s)
85  {
86  // use some arbitrary position to evaluate jacobians, they are constant
87  static const typename Traits::DomainLocal xl(0);
88 
89  // precompute Jacobian (gradients) of the p1 element
90  P1Basis(p1LocalBasis, geo).evaluateJacobian(xl, p1j);
91 
92  // calculate edge sizes and orientations
93  for(std::size_t i = 0; i < s; ++i) {
94  edgel[i] = (geo.corner(refelem.subEntity(i,dim-1,0,dim))-
95  geo.corner(refelem.subEntity(i,dim-1,1,dim))
96  ).two_norm();
97  const typename VertexOrder::iterator& edgeVertexOrder =
98  vertexOrder.begin(dim-1, i);
99  if(edgeVertexOrder[0] > edgeVertexOrder[1])
100  edgel[i] *= -1;
101  }
102  }
103 
105  std::size_t size () const { return s; }
106 
108  void evaluateFunction(const typename Traits::DomainLocal& xl,
109  std::vector<typename Traits::Range>& out) const
110  {
111  out.assign(s, typename Traits::Range(0));
112 
113  // compute p1 values -- use the local basis directly for that, local and
114  // global values are identical for scalars
115  std::vector<typename P1LocalBasis::Traits::RangeType> p1v;
116  p1LocalBasis.evaluateFunction(xl, p1v);
117 
118  for(std::size_t i = 0; i < s; i++) {
119  const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
120  const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
121  out[i].axpy( p1v[i0], p1j[i1][0]);
122  out[i].axpy(-p1v[i1], p1j[i0][0]);
123  out[i] *= edgel[i];
124  }
125  }
126 
128  void evaluateJacobian(const typename Traits::DomainLocal&,
129  std::vector<typename Traits::Jacobian>& out) const
130  {
131  out.resize(s);
132 
133  for(std::size_t i = 0; i < s; i++) {
134  const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
135  const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
136  for(std::size_t j = 0; j < dim; j++)
137  for(std::size_t k = 0; k < dim; k++)
138  out[i][j][k] = edgel[i] *
139  (p1j[i0][0][k]*p1j[i1][0][j]-p1j[i1][0][k]*p1j[i0][0][j]);
140  }
141  }
142 
144  std::size_t order () const { return 1; }
145  };
146 
147  template<class Geometry, class RF>
148  const typename EdgeS0_5Basis<Geometry, RF>::P1LocalBasis&
149  EdgeS0_5Basis<Geometry, RF>::p1LocalBasis = P1LocalBasis();
150 
151 } // namespace Dune
152 
153 #endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH