4 #ifndef DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
5 #define DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH
10 #include <dune/common/fmatrix.hh>
11 #include <dune/common/fvector.hh>
33 template<
class Geometry,
class RF>
35 private EdgeS0_5Common<Geometry::mydimension, typename Geometry::ctype>
48 typedef FieldVector<RangeField, dimRange>
Range;
50 typedef FieldMatrix<RangeField, dimRange, dimDomainGlobal>
Jacobian;
61 static const P1LocalBasis& p1LocalBasis;
69 std::vector<typename P1Basis::Traits::Jacobian> p1j;
71 std::vector<typename Traits::DomainField> edgel;
81 template<
typename VertexOrder>
92 for(std::size_t i = 0; i <
s; ++i) {
93 edgel[i] = (geo.corner(
refelem.subEntity(i,dim-1,0,dim))-
94 geo.corner(
refelem.subEntity(i,dim-1,1,dim))
96 const typename VertexOrder::iterator& edgeVertexOrder =
97 vertexOrder.begin(dim-1, i);
98 if(edgeVertexOrder[0] > edgeVertexOrder[1])
104 std::size_t
size ()
const {
return s; }
108 std::vector<typename Traits::Range>& out)
const
114 std::vector<typename P1LocalBasis::Traits::RangeType> p1v;
115 p1LocalBasis.evaluateFunction(xl, p1v);
117 for(std::size_t i = 0; i <
s; i++) {
118 const std::size_t i0 =
refelem.subEntity(i,dim-1,0,dim);
119 const std::size_t i1 =
refelem.subEntity(i,dim-1,1,dim);
120 out[i].axpy( p1v[i0], p1j[i1][0]);
121 out[i].axpy(-p1v[i1], p1j[i0][0]);
128 std::vector<typename Traits::Jacobian>& out)
const
132 for(std::size_t i = 0; i <
s; i++) {
133 const std::size_t i0 =
refelem.subEntity(i,dim-1,0,dim);
134 const std::size_t i1 =
refelem.subEntity(i,dim-1,1,dim);
135 for(std::size_t j = 0; j < dim; j++)
136 for(std::size_t k = 0; k < dim; k++)
137 out[i][j][k] = edgel[i] *
138 (p1j[i0][0][k]*p1j[i1][0][j]-p1j[i1][0][k]*p1j[i0][0][j]);
145 std::vector<typename Traits::Range>& out)
const
147 auto totalOrder = std::accumulate(
order.begin(),
order.end(), 0);
148 if (totalOrder == 0) {
150 }
else if (totalOrder==1) {
151 auto const k = std::distance(
order.begin(), std::find(
order.begin(),
order.end(), 1));
154 for (std::size_t i = 0; i <
s; i++)
156 const std::size_t i0 =
refelem.subEntity(i,dim-1,0,dim);
157 const std::size_t i1 =
refelem.subEntity(i,dim-1,1,dim);
158 for(std::size_t j = 0; j < dim; j++)
159 out[i][j] = edgel[i] *
160 (p1j[i0][0][k]*p1j[i1][0][j] - p1j[i1][0][k]*p1j[i0][0][j]);
163 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
168 std::size_t
order ()
const {
return 1; }
171 template<
class Geometry,
class RF>
172 const typename EdgeS0_5Basis<Geometry, RF>::P1LocalBasis&
173 EdgeS0_5Basis<Geometry, RF>::p1LocalBasis = P1LocalBasis();
177 #endif // DUNE_LOCALFUNCTIONS_WHITNEY_EDGES0_5_BASIS_HH