DUNE PDELab (git)

cubichermitebasis.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// 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#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
7#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_CUBICHERMITEBASIS_HH
8
9#include <algorithm>
10#include <array>
11#include <type_traits>
12#include <vector>
13
18
20
21#include <dune/localfunctions/common/localbasis.hh>
22#include <dune/localfunctions/common/localfiniteelementtraits.hh>
23#include <dune/localfunctions/common/localkey.hh>
24
25#include <dune/functions/common/mapperutilities.hh>
26#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
27#include <dune/functions/functionspacebases/functionaldescriptor.hh>
28#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
29#include <dune/functions/functionspacebases/nodes.hh>
30#include <dune/functions/functionspacebases/transformedfiniteelementmixin.hh>
31
32
33
34#include <dune/functions/analyticfunctions/monomialset.hh>
35
51namespace Dune::Functions
52{
53
54 template<class GV, class R, bool reduced>
55 struct CubicHermitePreBasis;
56
76 template<class GV,bool reduced = false, class R = double>
78
79 template<class DF, int n, class D, class RF, int m, class R, class J, class H>
80 struct H2LocalBasisTraits
81 : public LocalBasisTraits<DF, n, D, RF, m, R, J>
82 {
88 using HessianType = H;
89 };
90
91 namespace Impl
92 {
93
94 // *****************************************************************************
95 // * Some helper functions for building polynomial bases from monomials
96 // *****************************************************************************
97
101 template<class K>
102 constexpr K const& squeeze(Dune::FieldVector<K,1> const& v){ return v[0]; }
103
106 template<class K>
107 constexpr K& squeeze(Dune::FieldVector<K,1>& v){ return v[0]; }
108
111 template<class K, int N>
112 constexpr Dune::FieldVector<K, N> const& squeeze(Dune::FieldMatrix<K,1,N> const& m){ return m[0]; }
113
116 template<class K, int N>
117 constexpr Dune::FieldVector<K, N>& squeeze(Dune::FieldMatrix<K,1,N>& m){ return m[0]; }
118
121 template<class Object>
122 constexpr Object& squeeze(Object& o){ return o; }
123
126 template<class Object>
127 constexpr Object const& squeeze(Object const& o){ return o; }
128
138 template<class KCoeff, int sizePolynom, int sizeMonom, class In, class Out>
139 void multiplyWithCoefficentMatrix(Dune::FieldMatrix<KCoeff, sizePolynom, sizeMonom> const& coefficients,
140 In const& monomialValues,
141 Out& polynomialValues)
142 {
143 for (int i = 0; i < sizePolynom; ++i)
144 {
145 squeeze(polynomialValues[i]) = 0;
146 for (int j = 0; j < sizeMonom; ++j)
147 squeeze(polynomialValues[i]) += coefficients[i][j]*monomialValues[j];
148 }
149 }
150
156 template<int dim, bool reduced>
157 class CubicHermiteLocalCoefficients
158 {
159 public:
160 using size_type = std::size_t;
161
162 CubicHermiteLocalCoefficients()
163 : localKeys_(size())
164 {
165 static_assert((dim > 0) and (dim <= 3), "CubicHermiteLocalCoefficients only implemented for dim=1,2,3");
166 static_assert((not reduced) or (dim == 2), "Reduced version of CubicHermiteLocalCoefficients only implemented for dim=2");
167 for (size_type i = 0; i < (dim +1); ++i)
168 {
169 // dim derivatives + 1 evaluation dofs per vertex
170 for (size_type k = 0; k < (dim +1); ++k)
171 localKeys_[(dim +1) * i + k] = LocalKey(i, dim, k);
172 }
173 if constexpr (not reduced)
174 {
175 // 1 evaluation per element (2d) / facets (3d)
176 for (size_type i = 0; i < (dim - 1) * (dim - 1); ++i)
177 localKeys_[(dim +1) * (dim +1) + i] = LocalKey(i, (dim == 2) ? 0 : 1, 0); // inner dofs
178 }
179 }
180
183 static constexpr size_type size()
184 {
185 if constexpr (dim==1)
186 return 4;
187 if constexpr ((dim==2) and (reduced))
188 return 9;
189 if constexpr ((dim==2) and (not reduced))
190 return 10;
191 if constexpr (dim==3)
192 return 20;
193 return 0;
194 }
195
198 LocalKey const &localKey(size_type i) const
199 {
200 return localKeys_[i];
201 }
202
203 private:
204 std::vector<LocalKey> localKeys_;
205 };
206
207
213 template<class D, class R, int dim, bool reduced>
214 class CubicHermiteReferenceLocalBasis
215 {
216 public:
217 using Traits = H2LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1, FieldVector<R, 1>,
218 FieldMatrix<R, 1, dim>, FieldMatrix<R, dim, dim>>;
219
220 private:
221
234 static constexpr auto getCubicHermiteCoefficients()
235 {
236 if constexpr (dim == 1)
237 return Dune::FieldMatrix<D, 4, 4>({{1, 0, -3, 2}, {0, 1, -2, 1}, {0, 0, 3, -2}, {0, 0, -1, 1}});
238 else if constexpr (dim == 2)
239 {
240 if constexpr (reduced) {
241 auto w = std::array<D, 9>{1. / 3, 1. / 18, 1. / 18, 1. / 3, -1. / 9,
242 1. / 18, 1. / 3, 1. / 18, -1. / 9};
244 {1, 0, 0, -3, -13 + w[0] * 27, -3, 2, 13 - w[0] * 27, 13 - w[0] * 27, 2},
245 {0, 1, 0, -2, -3 + w[1] * 27, 0, 1, 3 - w[1] * 27, 2 - w[1] * 27, 0},
246 {0, 0, 1, 0, -3 + w[2] * 27, -2, 0, 2 - w[2] * 27, 3 - w[2] * 27, 1},
247 {0, 0, 0, 3, -7 + w[3] * 27, 0, -2, 7 - w[3] * 27, 7 - w[3] * 27, 0},
248 {0, 0, 0, -1, 2 + w[4] * 27, 0, 1, -2 - w[4] * 27, -2 - w[4] * 27, 0},
249 {0, 0, 0, 0, -1 + w[5] * 27, 0, 0, 2 - w[5] * 27, 1 - w[5] * 27, 0},
250 {0, 0, 0, 0, -7 + w[6] * 27, 3, 0, 7 - w[6] * 27, 7 - w[6] * 27, -2},
251 {0, 0, 0, 0, -1 + w[7] * 27, 0, 0, 1 - w[7] * 27, 2 - w[7] * 27, 0},
252 {0, 0, 0, 0, 2 + w[8] * 27, -1, 0, -2 - w[8] * 27, -2 - w[8] * 27, 1},
253 });
254 }
255 else
257 {1, 0, 0, -3, -13, -3, 2, 13, 13, 2},
258 {0, 1, 0, -2, -3, 0, 1, 3, 2, 0},
259 {0, 0, 1, 0, -3, -2, 0, 2, 3, 1}, // l_2
260 {0, 0, 0, 3, -7, 0, -2, 7, 7, 0},
261 {0, 0, 0, -1, 2, 0, 1, -2, -2, 0},
262 {0, 0, 0, 0, -1, 0, 0, 2, 1, 0},
263 {0, 0, 0, 0, -7, 3, 0, 7, 7, -2}, // l_6
264 {0, 0, 0, 0, -1, 0, 0, 1, 2, 0},
265 {0, 0, 0, 0, 2, -1, 0, -2, -2, 1},
266 {0, 0, 0, 0, 27, 0, 0, -27, -27, 0}}); // l_9, inner dof
267 }
268 else if constexpr (dim == 3)
269 {
270 return Dune::FieldMatrix<D, 20,20>({{1, 0, 0, 0, -3, -13, -3, -13, -13, -3, // deg 0 to 2
271 2, 13, 13, 2, 13, 33, 13, 13, 13, 2}, // deg 3
272 {0, 1, 0, 0,/*xx*/ -2, /*xy*/-3,/*yy*/ 0,/*xz*/ -3,/*yz*/ 0,/*zz*/ 0, 1, 3, 2, 0, 3, 4, 0, 2, 0, 0},
273 {0, 0, 1, 0, 0, -3, -2, 0, -3, 0, 0, 2, 3, 1, 0, 4, 3, 0, 2, 0},
274 {0, 0, 0, 1, 0, 0, 0, -3, -3, -2, 0, 0, 0, 0, 2, 4, 2, 3, 3, 1},
275 {0, 0, 0, 0, 3, -7, 0, -7, 0, 0, // l_4
276 -2, 7, 7, 0, 7, 7, 0, 7, 0, 0},
277 {0, 0, 0, 0, -1, 2, 0, 2, 0, 0, 1, -2, -2, 0, -2, -2, 0, -2, 0, 0},
278 {0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 2, 1, 0, 0, 0, 0, 0, 0, 0},
279 {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0},
280 {0, 0, 0, 0, 0, -7, 3, 0, -7, 0, // l_8
281 0, 7, 7, -2, 0, 7, 7, 0, 7, 0},
282 {0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0},
283 {0, 0, 0, 0, 0, 2, -1, 0, 2, 0, 0, -2, -2, 1, 0, -2, -2, 0, -2, 0},
284 {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 1, 0},
285 {0, 0, 0, 0, 0, 0, 0, -7, -7, 3, // l_12
286 0, 0, 0, 0, 7, 7, 7, 7, 7, -2},
287 {0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0},
288 {0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 2, 0},
289 {0, 0, 0, 0, 0, 0, 0, 2, 2, -1, 0, 0, 0, 0, -2, -2, -2, -2, -2, 1},
290 // l_16, from here on inner dofs
291 {0, 0, 0, 0, 0, 27, 0, 0, 0, 0, // bottom
292 0, -27, -27, 0, 0, -27, 0, 0, 0, 0},
293 {0, 0, 0, 0, 0, 0, 0, 27, 0, 0, // front
294 0, 0, 0, 0, -27, -27, 0, -27, 0, 0},
295 {0, 0, 0, 0, 0, 0, 0, 0, 27, 0, // left
296 0, 0, 0, 0, 0, -27, -27, 0, -27, 0},
297 {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // right
298 0, 0, 0, 0, 0, 27, 0, 0, 0, 0}});
299 }
300 }
301
302 static constexpr auto referenceBasisCoefficients = getCubicHermiteCoefficients();
303 static constexpr MonomialSet<typename Traits::RangeFieldType, dim, 3> monomials = {};
304
305 public:
306
307 CubicHermiteReferenceLocalBasis()
308 {
309 static_assert((dim > 0) and (dim <= 3), "CubicHermiteReferenceLocalBasis only implemented for dim=1,2,3");
310 static_assert((not reduced) or (dim == 2), "Reduced version of CubicHermiteReferenceLocalBasis only implemented for dim=2");
311 }
312
315 static constexpr unsigned int size()
316 {
317 return CubicHermiteLocalCoefficients<dim,reduced>::size();
318 }
319
322 unsigned int order() const
323 {
324 return 3;
325 }
326
332 void evaluateFunction(const typename Traits::DomainType &in,
333 std::vector<typename Traits::RangeType> &out) const
334 {
335 out.resize(size());
336 auto monomialValues = monomials(in);
337 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
338 }
339
345 void evaluateJacobian(const typename Traits::DomainType &in,
346 std::vector<typename Traits::JacobianType> &out) const
347 {
348 out.resize(size());
349 auto monomialValues = derivative(monomials)(in);
350 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
351 }
352
358 void evaluateHessian(const typename Traits::DomainType &in,
359 std::vector<typename Traits::HessianType> &out) const
360 {
361 out.resize(size());
362 auto monomialValues = derivative(derivative(monomials))(in);
363 multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
364 }
365
372 void partial(std::array<unsigned int, dim> order, const typename Traits::DomainType &in,
373 std::vector<typename Traits::RangeType> &out) const
374 {
375 out.resize(size());
376 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
377 if (totalOrder == 0)
378 evaluateFunction(in, out);
379 else if (totalOrder == 1)
380 {
381 evaluateJacobian(in,jacobiansBuffer_);
382 std::size_t which = std::max_element(order.begin(), order.end()) - order.begin();
383 for (auto i : Dune::range(size()))
384 out[i] = jacobiansBuffer_[i][0][which];
385 }
386 else if (totalOrder == 2)
387 {
388 evaluateHessian(in, hessianBuffer_);
389 std::size_t first, second;
390 first = std::max_element(order.begin(), order.end()) - order.begin();
391 if (order[first] == 2)
392 second = first;
393 else
394 {
395 order[first] = 0;
396 second = std::max_element(order.begin(), order.end()) - order.begin();
397 }
398 for (auto i : Dune::range(size()))
399 out[i] = hessianBuffer_[i][first][second];
400 }
401 else
402 DUNE_THROW(RangeError, "partial() not implemented for given order");
403 }
404
405 private:
406 mutable std::vector<typename Traits::JacobianType> jacobiansBuffer_;
407 mutable std::vector<typename Traits::HessianType> hessianBuffer_;
408 };
409
410
416 template<class D, int dim, bool reduced = false>
417 class CubicHermiteLocalInterpolation
418 {
419 using size_type = std::size_t;
420
421 static constexpr unsigned int size()
422 {
423 return CubicHermiteLocalCoefficients<dim,reduced>::size();
424 }
425
426 using FunctionalDescriptor = Dune::Functions::Impl::FunctionalDescriptor<dim>;
427
428 public:
429
430 CubicHermiteLocalInterpolation()
431 {
432 static_assert((dim > 0) and (dim <= 3), "CubicHermiteLocalInterpolation only implemented for dim=1,2,3");
433 static_assert((not reduced) or (dim == 2), "Reduced version of CubicHermiteLocalInterpolation only implemented for dim=2");
434 if constexpr (dim==1)
435 {
436 descriptors_[0] = FunctionalDescriptor();
437 descriptors_[1] = FunctionalDescriptor({1});
438 descriptors_[2] = FunctionalDescriptor();
439 descriptors_[3] = FunctionalDescriptor({1});
440 }
441 if constexpr (dim==2)
442 {
443 descriptors_[0] = FunctionalDescriptor();
444 descriptors_[1] = FunctionalDescriptor({1,0});
445 descriptors_[2] = FunctionalDescriptor({0,1});
446 descriptors_[3] = FunctionalDescriptor();
447 descriptors_[4] = FunctionalDescriptor({1,0});
448 descriptors_[5] = FunctionalDescriptor({0,1});
449 descriptors_[6] = FunctionalDescriptor();
450 descriptors_[7] = FunctionalDescriptor({1,0});
451 descriptors_[8] = FunctionalDescriptor({0,1});
452 if (not reduced)
453 descriptors_[9] = FunctionalDescriptor();
454 }
455 if constexpr (dim==3)
456 {
457 descriptors_[0] = FunctionalDescriptor();
458 descriptors_[1] = FunctionalDescriptor({1,0,0});
459 descriptors_[2] = FunctionalDescriptor({0,1,0});
460 descriptors_[3] = FunctionalDescriptor({0,0,1});
461 descriptors_[4] = FunctionalDescriptor();
462 descriptors_[5] = FunctionalDescriptor({1,0,0});
463 descriptors_[6] = FunctionalDescriptor({0,1,0});
464 descriptors_[7] = FunctionalDescriptor({0,0,1});
465 descriptors_[8] = FunctionalDescriptor();
466 descriptors_[9] = FunctionalDescriptor({1,0,0});
467 descriptors_[10] = FunctionalDescriptor({0,1,0});
468 descriptors_[11] = FunctionalDescriptor({0,0,1});
469 descriptors_[12] = FunctionalDescriptor();
470 descriptors_[13] = FunctionalDescriptor({1,0,0});
471 descriptors_[14] = FunctionalDescriptor({0,1,0});
472 descriptors_[15] = FunctionalDescriptor({0,0,1});
473 descriptors_[16] = FunctionalDescriptor();
474 descriptors_[17] = FunctionalDescriptor();
475 descriptors_[18] = FunctionalDescriptor();
476 descriptors_[19] = FunctionalDescriptor();
477 }
478 }
479
482 template<class Element>
483 void bind( Element const &element, std::array<D, dim+1>const& averageVertexMeshSize)
484 {
485 averageVertexMeshSize_ = &averageVertexMeshSize;
486 }
487
495 template<class F, class C>
496 void interpolate(const F &f, std::vector<C> &out) const
497 {
498 out.resize(size());
499 auto df = derivative(f);
500 auto const &refElement = Dune::ReferenceElements<D, dim>::simplex();
501
502 // Iterate over vertices, dim derivative +1 evaluation dofs per vertex
503 for (int i = 0; i < (dim+1); ++i)
504 {
505 auto x = refElement.position(i, dim);
506 auto&& derivativeValue = df(x);
507 out[i * (dim +1)] = f(x);
508 for (int d = 0; d < dim; ++d)
509 out[i * (dim+1) + d + 1] = squeeze(derivativeValue)[d] * (*averageVertexMeshSize_)[i];
510 }
511
512 if constexpr (not reduced)
513 {
514 for (size_type i = 0; i < (dim - 1) * (dim - 1); ++i)
515 out[(dim +1) * (dim +1) + i] = f(refElement.position(i, (dim == 2) ? 0 : 1));
516 }
517 }
518
522 const FunctionalDescriptor& functionalDescriptor(size_type i) const
523 {
524 return descriptors_[i];
525 }
526
527 protected:
528 std::array<D, dim+1> const* averageVertexMeshSize_;
529 std::array<FunctionalDescriptor, size()> descriptors_;
530 };
531
532 template<class D, class R, int dim , bool reduced>
533 struct CubicHermiteLocalBasisTraits
534 : public H2LocalBasisTraits<D, dim, Dune::FieldVector<D,dim>, R, 1,
535 Dune::FieldVector<R,1>, Dune::FieldMatrix<R,1,dim>, Dune::FieldMatrix<R,dim,dim>>
536 {};
537
546 template<class D, class R, int dim, bool reduced = false>
547 class CubicHermiteLocalFiniteElement
548 : public Impl::TransformedFiniteElementMixin<CubicHermiteLocalFiniteElement<D,R,dim,reduced>, CubicHermiteLocalBasisTraits<D, R, dim, reduced>>
549 {
550 using Base = Impl::TransformedFiniteElementMixin< CubicHermiteLocalFiniteElement<D,R,dim,reduced>, CubicHermiteLocalBasisTraits<D, R, dim, reduced>>;
551 friend class Impl::TransformedLocalBasis<CubicHermiteLocalFiniteElement<D,R,dim,reduced>, CubicHermiteLocalBasisTraits<D, R, dim, reduced>>;
552
553 public:
554
555 CubicHermiteLocalFiniteElement()
556 : Base()
557 {
558 static_assert((dim > 0) and (dim <= 3), "CubicHermiteLocalFiniteElement only implemented for dim=1,2,3");
559 static_assert((not reduced) or (dim == 2), "Reduced version of CubicHermiteLocalFiniteElement only implemented for dim=2");
560 }
561
564 using size_type = std::size_t;
565 using Traits = LocalFiniteElementTraits<
566 Impl::TransformedLocalBasis<CubicHermiteLocalFiniteElement<D,R,dim,reduced>, CubicHermiteLocalBasisTraits<D, R, dim, reduced>>,
567 Impl::CubicHermiteLocalCoefficients<dim, reduced>,
568 Impl::CubicHermiteLocalInterpolation<D, dim, reduced>>;
569
573 const typename Traits::LocalCoefficientsType &localCoefficients() const
574 {
575 return coefficients_;
576 }
577
580 const typename Traits::LocalInterpolationType &localInterpolation() const
581 {
582 return interpolation_;
583 }
584
587 static constexpr GeometryType type()
588 {
589 return GeometryTypes::simplex(dim);
590 }
591
594 static constexpr size_type size()
595 {
596 return Impl::CubicHermiteLocalCoefficients<dim,reduced>::size();
597 }
598
601 template<class Mapper, class Element>
602 void bind(Mapper const& vertexMapper, std::vector<D> const& globalAverageVertexMeshSize, Element const &e)
603 {
604 // Cache average mesh size for each vertex
605 for (auto i : range(dim+1))
606 averageVertexMeshSize_[i] = globalAverageVertexMeshSize[vertexMapper.subIndex(e, i, dim)];
607
608 // Bind LocalInterpolation to updated local state
609 interpolation_.bind(e, averageVertexMeshSize_);
610
611 // Compute local transformation matrices for each vertex
612 const auto& geometry = e.geometry();
614 for (auto i : range(dim+1))
615 {
616 scaledVertexJacobians_[i] = geometry.jacobian(refElement.position(i, dim));
617 scaledVertexJacobians_[i] /= averageVertexMeshSize_[i];
618 }
619 }
620
621 protected:
622
625 Impl::CubicHermiteReferenceLocalBasis<D, R, dim, reduced> const& referenceLocalBasis() const
626 {
627 return basis_;
628 }
629
634 template<class InputValues, class OutputValues>
635 void transform(InputValues const &inValues, OutputValues &outValues) const
636 {
637 assert(inValues.size() == size());
638 assert(outValues.size() == inValues.size());
639 auto inIt = inValues.begin();
640 auto outIt = outValues.begin();
641
642 for (auto vertex : Dune::range((dim +1)))
643 {
644 *outIt = *inIt; // value dof is not transformed
645 outIt++, inIt++;
646 // transform the gradient dofs together
647 for (auto &&[row_i, i] : sparseRange(scaledVertexJacobians_[vertex]))
648 {
649 outIt[i] = 0.;
650 for (auto &&[val_i_j, j] : sparseRange(row_i))
651 outIt[i] += val_i_j * inIt[j];
652 }
653 // increase pointer by size of gradient = dim
654 outIt += dim, inIt += dim;
655 }
656
657 // For the non-reduced case: Copy all remaining inner dofs
658 if constexpr (dim > 1 and (not reduced))
659 std::copy(inIt, inValues.end(), outIt);
660 }
661
662 private:
663
664 typename Impl::CubicHermiteReferenceLocalBasis<D, R, dim, reduced> basis_;
665 typename Traits::LocalCoefficientsType coefficients_;
666 typename Traits::LocalInterpolationType interpolation_;
667 // the transformation to correct the lack of affine equivalence boils down to
668 // one transformation matrix per vertex
669 std::array<Dune::FieldMatrix<R, dim, dim>, dim+1> scaledVertexJacobians_;
670 // the local state, i.e. a collection of global information restricted to this element
671 std::array<D, dim+1> averageVertexMeshSize_;
672
673 };
674
675 } // end namespace Impl
676
677
678
679 // *****************************************************************************
680 // This is the reusable part of the basis. It contains
681 //
682 // CubicHermitePreBasis
683 // CubicHermiteNode
684 //
685 // The pre-basis allows to create the others and is the owner of possible shared
686 // state. These components do _not_ depend on the global basis and local view
687 // and can be used without a global basis.
688 // *****************************************************************************
689
690 template<class GV, class R, bool reduced>
691 class CubicHermiteNode
692 : public LeafBasisNode
693 {
695
696 public:
697 using size_type = std::size_t;
698 using Element = typename GV::template Codim<0>::Entity;
699 using FiniteElement = typename Impl::CubicHermiteLocalFiniteElement<typename GV::ctype, R, GV::dimension, reduced>;
700
701 CubicHermiteNode(Mapper const& m, std::vector<typename GV::ctype> const& averageVertexMeshSize)
702 : element_(nullptr)
703 , vertexMapper_(&m)
704 , averageVertexMeshSize_(&averageVertexMeshSize)
705 {
706 this->setSize(finiteElement_.size());
707 }
708
710 Element const &element() const
711 {
712 return *element_;
713 }
714
720 FiniteElement const &finiteElement() const
721 {
722 return finiteElement_;
723 }
724
726 void bind(Element const &e)
727 {
728 element_ = &e;
729 finiteElement_.bind(*vertexMapper_, *averageVertexMeshSize_, *element_);
730 }
731
733 unsigned int order() const { return finiteElement_.localBasis().order(); }
734
735 protected:
736 FiniteElement finiteElement_;
737 Element const* element_;
738 Mapper const* vertexMapper_;
739 std::vector<typename GV::ctype> const* averageVertexMeshSize_;
740 };
741
742
752 template<class GV, class R, bool reduced = false>
753 class CubicHermitePreBasis
754 : public LeafPreBasisMapperMixin<GV>
755 {
756 using Base = LeafPreBasisMapperMixin<GV>;
757 using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
758 using Element = typename GV::template Codim<0>::Entity;
759 using D = typename GV::ctype;
760 static const std::size_t dim = GV::dimension;
761
762 // helper methods to assign each subentity the number of dofs. Used by the LeafPreBasisMapperMixin.
763 static constexpr auto cubicHermiteMapperLayout(Dune::GeometryType type, int gridDim)
764 {
765 if (type.isVertex())
766 return 1 + gridDim; // one evaluation dof and gridDim derivative dofs per vertex
767 if (gridDim == 1) // in 1d there are no other dofs
768 return 0;
769 // in 2d we have one inner dof (i.e. on the triangle) or non for the reduced case
770 // and in 3d we have one dof on each face (i.e. on each triangle)
771 if ((type.isTriangle()) and (not reduced))
772 return 1;
773 else
774 return 0; // this case is only entered for the interior of the 3d element. There are no dofs.
775 }
776
777 public:
779 using GridView = GV;
780
782 using size_type = std::size_t;
783
785 using Node = CubicHermiteNode<GridView, R,reduced>;
786
787 public:
788
791 : Base(gv, cubicHermiteMapperLayout)
792 , vertexMapper_({gv, mcmgVertexLayout()})
793 {
794 static_assert((dim > 0) and (dim <= 3), "CubicHermitePreBasis only implemented for dim=1,2,3");
795 static_assert((not reduced) or (dim == 2), "Reduced version of CubicHermitePreBasis only implemented for dim=2");
796 averageVertexMeshSize_ = Impl::computeAverageSubEntityMeshSize<D>(vertexMapper_);
797 }
798
800 void update(GridView const &gv)
801 {
802 Base::update(gv);
803 vertexMapper_.update(this->gridView());
804 averageVertexMeshSize_ = Impl::computeAverageSubEntityMeshSize<D>(vertexMapper_);
805 }
806
811 {
812 return Node{vertexMapper_, averageVertexMeshSize_};
813 }
814
815 protected:
816
817 SubEntityMapper vertexMapper_;
818 std::vector<D> averageVertexMeshSize_;
819
820 }; // class CubicHermitePreBasis
821
822 namespace BasisFactory
823 {
824
832 template<class R = double>
833 auto cubicHermite()
834 {
835 return [=](auto const &gridView) {
836 return CubicHermitePreBasis<std::decay_t<decltype(gridView)>, R>(gridView);
837 };
838 }
839
847 template<class R = double>
849 {
850 return [=](auto const &gridView) {
851 return CubicHermitePreBasis<std::decay_t<decltype(gridView)>, R, true>(gridView);
852 };
853 }
854
855 } // end namespace BasisFactory
856
857
858} // end namespace Dune::Functions
859
860#endif
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Global basis for given pre-basis.
Definition: defaultglobalbasis.hh:50
A generic MixIn class for PreBasis with flat indices computed from a mapper.
Definition: leafprebasismappermixin.hh:62
const GridView & gridView() const
Export the stored GridView.
Definition: leafprebasismappermixin.hh:95
void update(const GridView &gv)
Update the stored GridView.
Definition: leafprebasismappermixin.hh:101
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:279
constexpr bool isTriangle() const
Return true if entity is a triangle.
Definition: type.hh:289
Implementation class for a multiple codim and multiple geometry type mapper.
Definition: mcmgmapper.hh:129
void update(const GV &gridView)
Recalculates indices after grid adaptation.
Definition: mcmgmapper.hh:308
A few common exception classes.
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
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.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:43
auto cubicHermite()
Create a pre-basis factory that can create a CubicHermite pre-basis.
Definition: cubichermitebasis.hh:511
auto reducedCubicHermite()
Create a pre-basis factory that can create a CubicHermite pre-basis.
Definition: cubichermitebasis.hh:524
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
MCMGLayout mcmgVertexLayout()
layout for vertices (dim-0 entities)
Definition: mcmgmapper.hh:107
auto sparseRange(Range &&range)
Allow structured-binding for-loops for sparse iterators.
Definition: rangeutilities.hh:722
Mapper for multiple codim and multiple geometry types.
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Utilities for reduction like operations on ranges.
CubicHermiteNode< GridView, R, reduced > Node
Template mapping root tree path to type of created tree node.
Definition: cubichermitebasis.hh:785
std::size_t size_type
Type used for indices and size information.
Definition: cubichermitebasis.hh:782
Node makeNode() const
Create tree node.
Definition: cubichermitebasis.hh:810
GV GridView
The grid view that the FE basis is defined on.
Definition: cubichermitebasis.hh:779
void update(GridView const &gv)
Update the stored grid view, to be called if the grid has changed.
Definition: cubichermitebasis.hh:800
CubicHermitePreBasis(const GV &gv)
Constructor for a given grid view object.
Definition: cubichermitebasis.hh:790
static const ReferenceElement & simplex()
get simplex reference elements
Definition: referenceelements.hh:162
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)