DUNE PDELab (git)

monomialset.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
5// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
6
7#ifndef DUNE_FUNCTIONS_ANALYTICFUNCTIONS_MONOMIALSET_HH
8#define DUNE_FUNCTIONS_ANALYTICFUNCTIONS_MONOMIALSET_HH
9
10#include <array>
11
14#include <dune/common/math.hh>
16
17
18
19namespace Dune::Functions {
20
21 namespace Impl {
22
23 // Computes `y = [1, x, x^2, x^3, ..., x^maxDegree]` with `x` a single coordinate
24 template<int maxDegree, class DomainFieldType, class RangeType>
25 void computePowers(const DomainFieldType& x, RangeType& y)
26 {
27 if constexpr(maxDegree >= 0)
28 {
29 y[0] = 1;
30 for(auto k : Dune::range(maxDegree))
31 y[k+1] = y[k]*x;
32 }
33 }
34
35 } // namespace Impl
36
37
62 template<class RangeFieldType, int dimension, int maxDegree>
64 {
65 static constexpr int dim = dimension;
66 static constexpr int size = Dune::binomial(maxDegree + dim, dim);
67 static_assert(1 <= dim && dim <= 3, "MonomialSet is specialized for dimensions 1,2,3 only.");
68
80 static constexpr std::array<std::array<std::size_t,dim>,size> exponents();
81
91 template<class DomainFieldType>
93
99 {
109 template<class DomainFieldType>
111
116 struct Hessian
117 {
127 template<class DomainFieldType>
128 constexpr std::array<FieldMatrix<RangeFieldType,dim, dim>, size> operator()(const Dune::FieldVector<DomainFieldType,dim>& x) const;
129
130 };
131
135 constexpr friend auto derivative(const Derivative & d)
136 {
137 return Hessian{};
138 }
139
140 };
141
145 constexpr friend auto derivative(const MonomialSet& m)
146 {
147 return Derivative{};
148 }
149
150 };
151
152
153#ifndef DOXYGEN
154 // Specialization for dim = 1
155 template<class RangeFieldType, int maxDegree>
156 struct MonomialSet<RangeFieldType, 1, maxDegree>
157 {
158 static constexpr int dim = 1;
159 static constexpr int size = maxDegree+1;
160
161 static constexpr auto exponents()
162 {
163 auto p = std::array<std::array<std::size_t,1>,size>{};
164 for(auto k : Dune::range(size))
165 p[k][0] = k;
166 return p;
167 }
168
169 template<class DomainFieldType>
170 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
171 {
173 Impl::computePowers<maxDegree>(x[0], y);
174 return y;
175 }
176
177 struct Derivative
178 {
179 template<class DomainFieldType>
180 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
181 {
183 Impl::computePowers<maxDegree-1>(x[0], xPowers);
184
186 for(auto degree : Dune::range(1, maxDegree+1))
187 y[degree][0] = degree*xPowers[degree-1];
188 return y;
189 }
190
191 struct Hessian
192 {
193 template<class DomainFieldType>
194 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,1>& x) const
195 {
196 auto xPowers = std::array<RangeFieldType,maxDegree+1>{};
197 Impl::computePowers<maxDegree-2>(x[0],xPowers);
198
199 auto y = std::array<Dune::FieldMatrix<RangeFieldType,1,1>,size>{};
200 for(auto degree : Dune::range(maxDegree+1))
201 if (degree-1 > 0)
202 y[degree][0][0] = degree*(degree-1)*xPowers[degree-2];
203 return y;
204 }
205
206 };
207
208 constexpr friend auto derivative(const Derivative& d)
209 {
210 return Hessian{};
211 }
212
213 };
214
215 constexpr friend auto derivative(const MonomialSet& m)
216 {
217 return Derivative{};
218 }
219
220 };
221
222
223 // Specialization for dim = 2
224 template<class RangeFieldType, int maxDegree>
225 struct MonomialSet<RangeFieldType, 2, maxDegree>
226 {
227 static constexpr int dim = 2;
228 static constexpr int size = (maxDegree+1)*(maxDegree+2)/2;
229
230 static constexpr auto exponents()
231 {
232 auto p = std::array<std::array<std::size_t,2>,size>{};
233 std::size_t index = 0;
234 for(auto degree : Dune::range(maxDegree+1))
235 {
236 for(auto k : Dune::range(degree+1))
237 {
238 p[index][0] = degree-k;
239 p[index][1] = k;
240 ++index;
241 }
242 }
243 return p;
244 }
245
246 template<class DomainFieldType>
247 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
248 {
249 auto xPowers = std::array<Dune::FieldVector<RangeFieldType,size>,dim>{};
250 for(auto j : Dune::range(dim))
251 Impl::computePowers<maxDegree>(x[j], xPowers[j]);
252
254 std::size_t index = 0;
255 for(auto degree : Dune::range(maxDegree+1))
256 {
257 for(auto k : Dune::range(degree+1))
258 {
259 y[index] = xPowers[0][degree-k]*xPowers[1][k];
260 ++index;
261 }
262 }
263 return y;
264 }
265
266 struct Derivative
267 {
268 template<class DomainFieldType>
269 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
270 {
271 auto xPowers = std::array<Dune::FieldVector<RangeFieldType,size>,dim>{};
272 for(auto j : Dune::range(dim))
273 Impl::computePowers<maxDegree-1>(x[j], xPowers[j]);
274
276 std::size_t index = 0;
277 for(auto degree : Dune::range(maxDegree+1))
278 {
279 for(auto k : Dune::range(degree+1))
280 {
281 if (degree-k > 0)
282 y[index][0] = (degree-k)*xPowers[0][degree-k-1]*xPowers[1][k];
283 if (k > 0)
284 y[index][1] = k*xPowers[0][degree-k]*xPowers[1][k-1];
285 ++index;
286 }
287 }
288 return y;
289 }
290
291 struct Hessian
292 {
293 template<class DomainFieldType>
294 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,2>& x) const
295 {
296 auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
297 for(auto j : Dune::range(dim))
298 Impl::computePowers<maxDegree-2>(x[j], xPowers[j]);
299
300 auto y = std::array<Dune::FieldMatrix<RangeFieldType,2,2>,size>{};
301 std::size_t index = 0;
302 for(auto degree : Dune::range(maxDegree+1))
303 {
304 for(auto k : Dune::range(degree+1))
305 {
306 if (degree-k > 1)
307 y[index][0][0] = (degree-k-1)*(degree-k)*xPowers[0][degree-k-2]*xPowers[1][k];
308 if (k > 0 and degree-k > 0){
309 auto mixed = k*(degree-k)*xPowers[0][degree-k-1]*xPowers[1][k-1];
310 y[index][0][1]= mixed;
311 y[index][1][0]= mixed;
312 }
313 if (k > 1)
314 y[index][1][1] = k*(k-1)*xPowers[0][degree-k]*xPowers[1][k-2];
315
316 ++index;
317 }
318 }
319 return y;
320 }
321
322 };
323
324 constexpr friend auto derivative(const Derivative & d)
325 {
326 return Hessian{};
327 }
328
329 };
330
331 constexpr friend auto derivative(const MonomialSet& m)
332 {
333 return Derivative{};
334 }
335
336 };
337
338
339 // Specialization for dim = 3
340 template<class RangeFieldType, int maxDegree>
341 struct MonomialSet<RangeFieldType, 3, maxDegree>
342 {
343 static constexpr int dim = 3;
344 static constexpr int size = Dune::binomial(maxDegree + dim, dim);
345
346 static constexpr auto exponents()
347 {
348 auto p = std::array<std::array<std::size_t,3>,size>{};
349 std::size_t index = 0;
350 for(auto degree : Dune::range(maxDegree+1))
351 {
352 for(auto k : Dune::range(degree+1))
353 {
354 for (auto l : Dune::range(degree-k+1))
355 {
356 p[index][0] = degree-k-l;
357 p[index][1] = l;
358 p[index][2] = k;
359 ++index;
360 }
361
362 }
363 }
364 return p;
365 }
366
367 template<class DomainFieldType>
368 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
369 {
370 auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
371 for(auto j : Dune::range(dim))
372 Impl::computePowers<maxDegree>(x[j], xPowers[j]);
373
375 std::size_t index = 0;
376 for(auto degree : Dune::range(maxDegree+1))
377 {
378 for(auto k : Dune::range(degree+1))
379 {
380 for (auto l : Dune::range(degree-k+1))
381 {
382 y[index] = xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k];
383 ++index;
384 }
385 }
386 }
387 return y;
388 }
389
390 struct Derivative
391 {
392 template<class DomainFieldType>
393 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
394 {
395 auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
396 for(auto j : Dune::range(dim))
397 {
398 xPowers[j][0] = 1.0;
399 for(auto k : Dune::range(maxDegree))
400 xPowers[j][k+1] = xPowers[j][k]*x[j];
401 }
402
404 std::size_t index = 0;
405 for(auto degree : Dune::range(maxDegree+1))
406 {
407 for(auto k : Dune::range(degree+1))
408 {
409 for (auto l : Dune::range(degree-k+1))
410 {
411 if (degree-k-l > 0)
412 y[index][0] = (degree-k-l)*xPowers[0][degree-k-l-1]*xPowers[1][l]*xPowers[2][k];
413 if (l > 0)
414 y[index][1] = l*xPowers[0][degree-k-l]*xPowers[1][l-1]*xPowers[2][k];
415 if (k > 0)
416 y[index][2] = k*xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k-1];
417 ++index;
418 }
419 }
420 }
421 return y;
422 }
423
424 struct Hessian
425 {
426 template<class DomainFieldType>
427 constexpr auto operator()(const Dune::FieldVector<DomainFieldType,3>& x) const
428 {
429 auto xPowers = std::array<std::array<RangeFieldType,maxDegree+1>,dim>{};
430 for(auto j : Dune::range(dim))
431 Impl::computePowers<maxDegree>(x[j], xPowers[j]);
432
433 auto y = std::array<Dune::FieldMatrix<RangeFieldType,3,3>,size>{};
434 std::size_t index = 0;
435 for(auto degree : Dune::range(maxDegree+1))
436 {
437 for(auto k : Dune::range(degree+1))
438 {
439 for (auto l : Dune::range(degree-k+1))
440 {
441 // xx
442 if (degree-k-l-1 > 0)
443 y[index][0][0] = (degree-k-l)*(degree-k-l-1)*xPowers[0][degree-k-l-2]*xPowers[1][l]*xPowers[2][k];
444 // xy and yx
445 if (degree-k-l > 0 and l > 0){
446 y[index][0][1] = (degree-k-l)*l*xPowers[0][degree-k-l-1]*xPowers[1][l-1]*xPowers[2][k];
447 y[index][1][0] = y[index][0][1];
448 }
449 // yy
450 if (l-1 > 0)
451 y[index][1][1] = l*(l-1)*xPowers[0][degree-k-l]*xPowers[1][l-2]*xPowers[2][k];
452 // xz and zx
453 if (k > 0 and degree-k-l > 0)
454 {
455 y[index][0][2] = (degree-k-l)*k*xPowers[0][degree-k-l-1]*xPowers[1][l]*xPowers[2][k-1];
456 y[index][2][0] = y[index][0][2];
457 }
458 // yz
459 if (l > 0 and k > 0)
460 {
461 y[index][1][2] = l*k*xPowers[0][degree-k-l]*xPowers[1][l-1]*xPowers[2][k-1];
462 y[index][2][1] = y[index][1][2];
463 }
464 // zz
465 if (k-1 > 0)
466 y[index][2][2] = (k-1)*k*xPowers[0][degree-k-l]*xPowers[1][l]*xPowers[2][k-2];
467 ++index;
468 }
469 }
470 }
471 return y;
472 }
473
474 };
475
476 constexpr friend auto derivative(const Derivative & d)
477 {
478 return Hessian{};
479 }
480
481 };
482
483 constexpr friend auto derivative(const MonomialSet& m)
484 {
485 return Derivative{};
486 }
487
488 };
489 #endif
490} // namespace Dune::Functions
491
492#endif // DUNE_FUNCTIONS_ANALYTICFUNCTIONS_MONOMIALSET_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
std::size_t degree(const Node &node)
Returns the degree of node as run time information.
Definition: nodeinterface.hh:79
Some useful basic math stuff.
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:131
Utilities for reduction like operations on ranges.
Set of all second order derivatives of monomials up to degree maxDegree as vector of matrix valued fu...
Definition: monomialset.hh:117
constexpr std::array< FieldMatrix< RangeFieldType, dim, dim >, size > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of Matrices of monomial derivatives.
Set of all first order derivatives of monomials up to degree maxDegree as vector of vector valued fun...
Definition: monomialset.hh:99
constexpr Dune::FieldMatrix< RangeFieldType, size, dim > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of arrays of monomial derivatives.
constexpr friend auto derivative(const Derivative &d)
Construct the Hessian object from a Derivative.
Definition: monomialset.hh:135
Function, which evaluates all monomials up to degree maxDegree in a given coordinate.
Definition: monomialset.hh:64
constexpr Dune::FieldVector< RangeFieldType, size > operator()(const Dune::FieldVector< DomainFieldType, dim > &x) const
Return array of monomial evaluations.
static constexpr std::array< std::array< std::size_t, dim >, size > exponents()
Return array of monomial exponents with shape size x dim
constexpr friend auto derivative(const MonomialSet &m)
Construct the Derivative object from a MonomialSet.
Definition: monomialset.hh:145
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)