dune-localfunctions  2.6-git
whitney/edges0.5/basis.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 
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 
53  private:
54  typedef Dune::P1LocalBasis<typename Traits::DomainField,
55  typename Traits::RangeField,
57  > P1LocalBasis;
59 
60  static const P1LocalBasis& p1LocalBasis;
61  static const std::size_t dim = Traits::dimDomainLocal;
62 
64  using Base::refelem;
65  using Base::s;
66 
67  // global values of the Jacobians (gradients) of the p1 basis
68  std::vector<typename P1Basis::Traits::Jacobian> p1j;
69  // edge sizes and orientations
70  std::vector<typename Traits::DomainField> edgel;
71 
72  public:
74 
80  template<typename VertexOrder>
81  EdgeS0_5Basis(const Geometry& geo, const VertexOrder& vertexOrder) :
82  p1j(s, typename P1Basis::Traits::Jacobian(0)), edgel(s)
83  {
84  // use some arbitrary position to evaluate jacobians, they are constant
85  static const typename Traits::DomainLocal xl(0);
86 
87  // precompute Jacobian (gradients) of the p1 element
88  P1Basis(p1LocalBasis, geo).evaluateJacobian(xl, p1j);
89 
90  // calculate edge sizes and orientations
91  for(std::size_t i = 0; i < s; ++i) {
92  edgel[i] = (geo.corner(refelem.subEntity(i,dim-1,0,dim))-
93  geo.corner(refelem.subEntity(i,dim-1,1,dim))
94  ).two_norm();
95  const typename VertexOrder::iterator& edgeVertexOrder =
96  vertexOrder.begin(dim-1, i);
97  if(edgeVertexOrder[0] > edgeVertexOrder[1])
98  edgel[i] *= -1;
99  }
100  }
101 
103  std::size_t size () const { return s; }
104 
106  void evaluateFunction(const typename Traits::DomainLocal& xl,
107  std::vector<typename Traits::Range>& out) const
108  {
109  out.assign(s, typename Traits::Range(0));
110 
111  // compute p1 values -- use the local basis directly for that, local and
112  // global values are identical for scalars
113  std::vector<typename P1LocalBasis::Traits::RangeType> p1v;
114  p1LocalBasis.evaluateFunction(xl, p1v);
115 
116  for(std::size_t i = 0; i < s; i++) {
117  const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
118  const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
119  out[i].axpy( p1v[i0], p1j[i1][0]);
120  out[i].axpy(-p1v[i1], p1j[i0][0]);
121  out[i] *= edgel[i];
122  }
123  }
124 
126  void evaluateJacobian(const typename Traits::DomainLocal&,
127  std::vector<typename Traits::Jacobian>& out) const
128  {
129  out.resize(s);
130 
131  for(std::size_t i = 0; i < s; i++) {
132  const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
133  const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
134  for(std::size_t j = 0; j < dim; j++)
135  for(std::size_t k = 0; k < dim; k++)
136  out[i][j][k] = edgel[i] *
137  (p1j[i0][0][k]*p1j[i1][0][j]-p1j[i1][0][k]*p1j[i0][0][j]);
138  }
139  }
140 
142  void partial (const std::array<unsigned int, dim>& order,
143  const typename Traits::DomainLocal& in, // position
144  std::vector<typename Traits::Range>& out) const // return value
145  {
146  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
147  if (totalOrder == 0) {
148  evaluateFunction(in, out);
149  } else if (totalOrder==1) {
150  auto const k = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
151  out.resize(size());
152 
153  for (std::size_t i = 0; i < s; i++)
154  {
155  const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
156  const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
157  for(std::size_t j = 0; j < dim; j++)
158  out[i][j] = edgel[i] *
159  (p1j[i0][0][k]*p1j[i1][0][j] - p1j[i1][0][k]*p1j[i0][0][j]);
160  }
161  } else {
162  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
163  }
164  }
165 
167  std::size_t order () const { return 1; }
168  };
169 
170  template<class Geometry, class RF>
173 
174 } // namespace Dune
175 
176 #endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
Geometry::ctype DomainField
Definition: whitney/edges0.5/basis.hh:40
static const std::size_t dimRange
Definition: whitney/edges0.5/basis.hh:47
RF RangeField
Definition: whitney/edges0.5/basis.hh:46
Common base class for edge elements.
Definition: common.hh:17
FieldVector< DomainField, dimDomainGlobal > DomainGlobal
Definition: whitney/edges0.5/basis.hh:44
Basis for order 0.5 (lowest order) edge elements on simplices.
Definition: whitney/edges0.5/basis.hh:34
Linear Lagrange shape functions on the simplex.
Definition: p1localbasis.hh:27
std::size_t size() const
number of shape functions
Definition: whitney/edges0.5/basis.hh:103
std::size_t s
The number of base functions.
Definition: common.hh:32
void evaluateJacobian(const typename Traits::DomainLocal &, std::vector< typename Traits::Jacobian > &out) const
Evaluate all Jacobians.
Definition: whitney/edges0.5/basis.hh:126
EdgeS0_5Basis(const Geometry &geo, const VertexOrder &vertexOrder)
Construct an EdgeS0_5Basis.
Definition: whitney/edges0.5/basis.hh:81
static const std::size_t dimDomainLocal
Definition: whitney/edges0.5/basis.hh:41
RefElem refelem
The reference element for this edge element.
Definition: common.hh:24
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:15
static const std::size_t dimDomainGlobal
Definition: whitney/edges0.5/basis.hh:42
FieldMatrix< RangeField, dimRange, dimDomainGlobal > Jacobian
Definition: whitney/edges0.5/basis.hh:50
FieldVector< RangeField, dimRange > Range
Definition: whitney/edges0.5/basis.hh:48
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: p1localbasis.hh:41
FieldVector< DomainField, dimDomainLocal > DomainLocal
Definition: whitney/edges0.5/basis.hh:43
Convert a simple scalar local basis into a global basis.
Definition: localtoglobaladaptors.hh:63
export type traits for function signature
Definition: whitney/edges0.5/basis.hh:39
void evaluateFunction(const typename Traits::DomainLocal &xl, std::vector< typename Traits::Range > &out) const
Evaluate all shape functions.
Definition: whitney/edges0.5/basis.hh:106
void partial(const std::array< unsigned int, dim > &order, const typename Traits::DomainLocal &in, std::vector< typename Traits::Range > &out) const
Evaluate partial derivatives of all shape functions.
Definition: whitney/edges0.5/basis.hh:142
std::size_t order() const
Polynomial order of the shape functions.
Definition: whitney/edges0.5/basis.hh:167