dune-localfunctions 3.0-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
17namespace Dune {
18
20 //
21 // Basis
22 //
23
25
33 template<class Geometry, class RF>
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:
57 typename Traits::RangeField,
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
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 void partial (const std::array<unsigned int, dim>& order,
145 const typename Traits::DomainLocal& in, // position
146 std::vector<typename Traits::Range>& out) const // return value
147 {
148 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
149 if (totalOrder == 0) {
150 evaluateFunction(in, out);
151 } else if (totalOrder==1) {
152 auto const k = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
153 out.resize(size());
154
155 for (std::size_t i = 0; i < s; i++)
156 {
157 const std::size_t i0 = refelem.subEntity(i,dim-1,0,dim);
158 const std::size_t i1 = refelem.subEntity(i,dim-1,1,dim);
159 for(std::size_t j = 0; j < dim; j++)
160 out[i][j] = edgel[i] *
161 (p1j[i0][0][k]*p1j[i1][0][j] - p1j[i1][0][k]*p1j[i0][0][j]);
162 }
163 } else {
164 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
165 }
166 }
167
169 std::size_t order () const { return 1; }
170 };
171
172 template<class Geometry, class RF>
173 const typename EdgeS0_5Basis<Geometry, RF>::P1LocalBasis&
174 EdgeS0_5Basis<Geometry, RF>::p1LocalBasis = P1LocalBasis();
175
176} // namespace Dune
177
178#endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
Definition brezzidouglasmarini1cube2d.hh:14
Convert a simple scalar local basis into a global basis.
Definition localtoglobaladaptors.hh:65
void evaluateJacobian(const typename Traits::DomainLocal &in, std::vector< typename Traits::Jacobian > &out) const
Definition localtoglobaladaptors.hh:125
Linear Lagrange shape functions on the simplex.
Definition p1localbasis.hh:28
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition p1localbasis.hh:41
Basis for order 0.5 (lowest order) edge elements on simplices.
Definition whitney/edges0.5/basis.hh:36
EdgeS0_5Basis(const Geometry &geo, const VertexOrder &vertexOrder)
Construct an EdgeS0_5Basis.
Definition whitney/edges0.5/basis.hh:83
void evaluateJacobian(const typename Traits::DomainLocal &, std::vector< typename Traits::Jacobian > &out) const
Evaluate all Jacobians.
Definition whitney/edges0.5/basis.hh:128
void evaluateFunction(const typename Traits::DomainLocal &xl, std::vector< typename Traits::Range > &out) const
Evaluate all shape functions.
Definition whitney/edges0.5/basis.hh:108
std::size_t size() const
number of shape functions
Definition whitney/edges0.5/basis.hh:105
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:144
std::size_t order() const
Polynomial order of the shape functions.
Definition whitney/edges0.5/basis.hh:169
export type traits for function signature
Definition whitney/edges0.5/basis.hh:39
FieldVector< DomainField, dimDomainGlobal > DomainGlobal
Definition whitney/edges0.5/basis.hh:44
static const std::size_t dimRange
Definition whitney/edges0.5/basis.hh:47
static const std::size_t dimDomainLocal
Definition whitney/edges0.5/basis.hh:41
FieldMatrix< RangeField, dimRange, dimDomainGlobal > Jacobian
Definition whitney/edges0.5/basis.hh:50
static const std::size_t diffOrder
Definition whitney/edges0.5/basis.hh:52
static const std::size_t dimDomainGlobal
Definition whitney/edges0.5/basis.hh:42
RF RangeField
Definition whitney/edges0.5/basis.hh:46
FieldVector< DomainField, dimDomainLocal > DomainLocal
Definition whitney/edges0.5/basis.hh:43
FieldVector< RangeField, dimRange > Range
Definition whitney/edges0.5/basis.hh:48
Geometry::ctype DomainField
Definition whitney/edges0.5/basis.hh:40
Common base class for edge elements.
Definition common.hh:15
static const ReferenceElement< DF, dim > & refelem
The reference element for this edge element.
Definition common.hh:17
static const std::size_t s
The number of base functions.
Definition common.hh:23