DUNE-FEM (unstable)

monomiallocalbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
4// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5#ifndef DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
6#define DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
7
8#include <array>
9#include <cassert>
10#include <numeric>
11
13#include <dune/common/math.hh>
14
15#include "../common/localbasis.hh"
16
17namespace Dune
18{
19 namespace MonomImp
20 {
22 template <typename Traits>
23 class EvalAccess {
24 std::vector<typename Traits::RangeType> &out;
25#ifndef NDEBUG
26 unsigned int first_unused_index;
27#endif
28
29 public:
30 EvalAccess(std::vector<typename Traits::RangeType> &out_)
31 : out(out_)
32#ifndef NDEBUG
33 , first_unused_index(0)
34#endif
35 { }
36#ifndef NDEBUG
37 ~EvalAccess() {
38 assert(first_unused_index == out.size());
39 }
40#endif
41 typename Traits::RangeFieldType &operator[](unsigned int index)
42 {
43 assert(index < out.size());
44#ifndef NDEBUG
45 if(first_unused_index <= index)
46 first_unused_index = index+1;
47#endif
48 return out[index][0];
49 }
50 };
51
53 template <typename Traits>
55 std::vector<typename Traits::JacobianType> &out;
56 unsigned int row;
57#ifndef NDEBUG
58 unsigned int first_unused_index;
59#endif
60
61 public:
62 JacobianAccess(std::vector<typename Traits::JacobianType> &out_,
63 unsigned int row_)
64 : out(out_), row(row_)
65#ifndef NDEBUG
66 , first_unused_index(0)
67#endif
68 { }
69#ifndef NDEBUG
71 assert(first_unused_index == out.size());
72 }
73#endif
74 typename Traits::RangeFieldType &operator[](unsigned int index)
75 {
76 assert(index < out.size());
77#ifndef NDEBUG
78 if(first_unused_index <= index)
79 first_unused_index = index+1;
80#endif
81 return out[index][0][row];
82 }
83 };
84
97 template <typename Traits, int c>
98 struct Evaluate
99 {
101 constexpr static int d = Traits::dimDomain - c;
102
109 template <typename Access>
110 static void eval (
111 const typename Traits::DomainType &in,
114 const std::array<unsigned int, Traits::dimDomain> &derivatives,
117 typename Traits::RangeFieldType prod,
119 int bound,
121 int& index,
123 Access &access)
124 {
125 // start with the highest exponent for this dimension, then work down
126 for (int e = bound; e >= 0; --e)
127 {
128 // the rest of the available exponents, to be used by the other
129 // dimensions
130 int newbound = bound - e;
131 if(e < (int)derivatives[d])
133 eval(in, derivatives, 0, newbound, index, access);
134 else {
135 int coeff = 1;
136 for(int i = e - derivatives[d] + 1; i <= e; ++i)
137 coeff *= i;
138 // call the evaluator for the next dimension
140 eval( // pass the coordinate and the derivatives unchanged
141 in, derivatives,
142 // also pass the product accumulated so far, but also
143 // include the current dimension
144 prod * power(in[d], e-derivatives[d]) * coeff,
145 // pass the number of remaining exponents to the next
146 // dimension
147 newbound,
148 // pass the next index to fill and the output access
149 // wrapper
150 index, access);
151 }
152 }
153 }
154 };
155
160 template <typename Traits>
161 struct Evaluate<Traits, 1>
162 {
163 constexpr static int d = Traits::dimDomain-1;
165 template <typename Access>
166 static void eval (const typename Traits::DomainType &in,
167 const std::array<unsigned int, Traits::dimDomain> &derivatives,
168 typename Traits::RangeFieldType prod,
169 int bound, int& index, Access &access)
170 {
171 if(bound < (int)derivatives[d])
172 prod = 0;
173 else {
174 int coeff = 1;
175 for(int i = bound - derivatives[d] + 1; i <= bound; ++i)
176 coeff *= i;
177 prod *= power(in[d], bound-derivatives[d]) * coeff;
178 }
179 access[index] = prod;
180 ++index;
181 }
182 };
183
184 } //namespace MonomImp
185
200 template<class D, class R, unsigned int d, unsigned int p>
202 {
203 // Helper: Number of shape functions for a k-th order element in dimension dd
204 static constexpr unsigned int size (int dd, int k)
205 {
206 if (dd==0 || k==0)
207 return 1;
208 return size(dd,k-1) + size(dd-1,k);
209 }
210
211 public:
215
217 static constexpr unsigned int size ()
218 {
219 return size(d,p);
220 }
221
223 inline void evaluateFunction (const typename Traits::DomainType& in,
224 std::vector<typename Traits::RangeType>& out) const
225 {
226 out.resize(size());
227 int index = 0;
228 std::array<unsigned int, d> derivatives;
229 std::fill(derivatives.begin(), derivatives.end(), 0);
231 for (unsigned int lp = 0; lp <= p; ++lp)
232 MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
233 }
234
240 inline void partial(const std::array<unsigned int,d>& order,
241 const typename Traits::DomainType& in,
242 std::vector<typename Traits::RangeType>& out) const
243 {
244 out.resize(size());
245 int index = 0;
247 for (unsigned int lp = 0; lp <= p; ++lp)
248 MonomImp::Evaluate<Traits, d>::eval(in, order, 1, lp, index, access);
249 }
250
252 inline void
253 evaluateJacobian (const typename Traits::DomainType& in, // position
254 std::vector<typename Traits::JacobianType>& out) const // return value
255 {
256 out.resize(size());
257 std::array<unsigned int, d> derivatives;
258 for(unsigned int i = 0; i < d; ++i)
259 derivatives[i] = 0;
260 for(unsigned int i = 0; i < d; ++i)
261 {
262 derivatives[i] = 1;
263 int index = 0;
265 for(unsigned int lp = 0; lp <= p; ++lp)
266 MonomImp::Evaluate<Traits, d>::eval(in, derivatives, 1, lp, index, access);
267 derivatives[i] = 0;
268 }
269 }
270
272 unsigned int order () const
273 {
274 return p;
275 }
276 };
277
278}
279
280#endif // DUNE_LOCALFUNCTIONS_MONOMIAL_MONOMIALLOCALBASIS_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Access output vector of evaluateFunction() and evaluate()
Definition: monomiallocalbasis.hh:23
Access output vector of evaluateJacobian()
Definition: monomiallocalbasis.hh:54
Definition: monomiallocalbasis.hh:202
unsigned int order() const
Polynomial order of the shape functions.
Definition: monomiallocalbasis.hh:272
void partial(const std::array< unsigned int, d > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition: monomiallocalbasis.hh:240
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: monomiallocalbasis.hh:253
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: monomiallocalbasis.hh:223
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
export type traits for function signature
Definition: monomiallocalbasis.hh:214
static constexpr unsigned int size()
Number of shape functions.
Definition: monomiallocalbasis.hh:217
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Some useful basic math stuff.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr Base power(Base m, Exponent p)
Power method for integer exponents.
Definition: math.hh:75
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
static void eval(const typename Traits::DomainType &in, const std::array< unsigned int, Traits::dimDomain > &derivatives, typename Traits::RangeFieldType prod, int bound, int &index, Access &access)
Definition: monomiallocalbasis.hh:166
Definition: monomiallocalbasis.hh:99
static constexpr int d
The next dimension to try for factors.
Definition: monomiallocalbasis.hh:101
static void eval(const typename Traits::DomainType &in, const std::array< unsigned int, Traits::dimDomain > &derivatives, typename Traits::RangeFieldType prod, int bound, int &index, Access &access)
Definition: monomiallocalbasis.hh:110
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)