DUNE-FEM (unstable)

localfiniteelementgeometry.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_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
6#define DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
7
8#include <cassert>
9#include <functional>
10#include <limits>
11#include <type_traits>
12#include <vector>
13
16#include <dune/common/math.hh>
18#include <dune/common/std/type_traits.hh>
19
20#include <dune/geometry/affinegeometry.hh> // for FieldMatrixHelper
22#include <dune/geometry/referenceelements.hh>
23#include <dune/geometry/type.hh>
24#include <dune/geometry/utility/algorithms.hh>
25#include <dune/geometry/utility/convergence.hh>
26
27namespace Dune {
28
38template <class LFE, int cdim>
40{
41 using LocalFiniteElement = LFE;
42 using LocalBasis = typename LFE::Traits::LocalBasisType;
43 using LocalBasisTraits = typename LocalBasis::Traits;
44
45public:
48
51
53 static const int coorddimension = cdim;
54
57
60
62 using Volume = decltype(power(std::declval<ctype>(),mydimension));
63
66
69
72
75
76public:
79 using ReferenceElement = typename ReferenceElements::ReferenceElement;
80
81protected:
82 using MatrixHelper = Impl::FieldMatrixHelper<ctype>;
83
84public:
87
103 LocalFiniteElementGeometry (const ReferenceElement& refElement,
104 const LocalFiniteElement& localFE,
105 std::vector<GlobalCoordinate> vertices)
106 : refElement_(refElement)
107 , localFE_(localFE)
108 , vertices_(std::move(vertices))
109 {
110 assert(localFE_.size() == vertices_.size());
111 }
112
126 template <class Param,
127 std::enable_if_t<std::is_invocable_r_v<GlobalCoordinate,Param,LocalCoordinate>, int> = 0>
128 LocalFiniteElementGeometry (const ReferenceElement& refElement,
129 const LocalFiniteElement& localFE,
130 Param&& parametrization)
131 : refElement_(refElement)
132 , localFE_(localFE)
133 {
134 localFE_.localInterpolation().interpolate(parametrization, vertices_);
135 }
136
143 template <class... Args>
144 explicit LocalFiniteElementGeometry (GeometryType gt, Args&&... args)
145 : LocalFiniteElementGeometry(ReferenceElements::general(gt), std::forward<Args>(args)...)
146 {}
147
149 int order () const
150 {
151 return localBasis().order();
152 }
153
159 bool affine () const
160 {
161 if (!affine_)
162 affine_.emplace(affineImpl());
163 return *affine_;
164 }
165
168 {
169 return refElement_.type();
170 }
171
173 int corners () const
174 {
175 return refElement_.size(mydimension);
176 }
177
180 {
181 assert( (i >= 0) && (i < corners()) );
182 return global(refElement_.position(i, mydimension));
183 }
184
187 {
188 return global(refElement_.position(0, 0));
189 }
190
203 {
204 thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
205 localBasis().evaluateFunction(local, shapeValues);
206 assert(shapeValues.size() == vertices_.size());
207
208 GlobalCoordinate out(0);
209 for (std::size_t i = 0; i < shapeValues.size(); ++i)
210 out.axpy(shapeValues[i], vertices_[i]);
211
212 return out;
213 }
214
232 LocalCoordinate local (const GlobalCoordinate& y, Impl::GaussNewtonOptions<ctype> opts = {}) const
233 {
234 LocalCoordinate x = refElement_.position(0,0);
235 Impl::GaussNewtonErrorCode err = Impl::gaussNewton(
236 [&](const LocalCoordinate& local) { return this->global(local); },
237 [&](const LocalCoordinate& local) { return this->jacobianTransposed(local); },
238 y, x, opts
239 );
240
241 if (err != Impl::GaussNewtonErrorCode::OK)
243 "Local coordinate can not be recovered from global coordinate, error code = " << int(err) << "\n"
244 << " (global(x) - y).two_norm() = " << (global(x) - y).two_norm()
245 << " > tol = " << opts.absTol);
246
247 return x;
248 }
249
261 {
262 return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
263 }
264
276 Volume volume (Impl::ConvergenceOptions<ctype> opts = {}) const
277 {
279 if (affine())
280 return vol0;
281
282 using std::abs;
283 for (int p = 2; p < opts.maxIt; ++p) {
285 if (abs(vol1 - vol0) < opts.absTol)
286 return vol1;
287
288 vol0 = vol1;
289 }
290 return vol0;
291 }
292
295 {
296 Volume vol(0);
297 for (const auto& qp : quadRule)
298 vol += integrationElement(qp.position()) * qp.weight();
299 return vol;
300 }
301
307 Jacobian jacobian (const LocalCoordinate& local) const
308 {
309 thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
310 localBasis().evaluateJacobian(local, shapeJacobians);
311 assert(shapeJacobians.size() == vertices_.size());
312
313 Jacobian out(0);
314 for (std::size_t i = 0; i < shapeJacobians.size(); ++i) {
315 for (int j = 0; j < Jacobian::rows; ++j) {
316 shapeJacobians[i].umtv(vertices_[i][j], out[j]);
317 }
318 }
319 return out;
320 }
321
328 {
329 return jacobian(local).transposed();
330 }
331
340 {
341 JacobianInverse out;
342 MatrixHelper::leftInvA(jacobian(local), out);
343 return out;
344 }
345
354 {
355 return jacobianInverse(local).transposed();
356 }
357
359 friend ReferenceElement referenceElement (const LocalFiniteElementGeometry& geometry)
360 {
361 return geometry.refElement_;
362 }
363
365 const LocalFiniteElement& finiteElement () const
366 {
367 return localFE_;
368 }
369
371 const std::vector<GlobalCoordinate>& coefficients () const
372 {
373 return vertices_;
374 }
375
377 const LocalBasis& localBasis () const
378 {
379 return localFE_.localBasis();
380 }
381
382private:
383
384 bool affineImpl () const
385 {
386 if constexpr(mydimension == 0)
387 // point geometries are always affine mappings
388 return true;
389 else {
390 if (order() > 1)
391 // higher-order parametrizations are by definition not affine
392 return false;
393 if constexpr(mydimension == 1)
394 // linear line geometries are affine mappings
395 return true;
396 else {
397 if (type().isSimplex())
398 // linear simplex geometries are affine mappings
399 return true;
400
401 // multi-linear mappings on non-simplex geometries might be affine
402 // as well. This is tested explicitly for all vertices by constructing
403 // an affine mapping from dim+1 affine-independent corners and evaluating
404 // at the other corners.
405 auto refSimplex = referenceElement<ctype,mydimension>(GeometryTypes::simplex(mydimension));
406 auto simplexIndices = Dune::range(refSimplex.size(mydimension));
407 auto simplexCorners = Dune::transformedRangeView(simplexIndices,
408 [&](int i) { return this->global(refSimplex.position(i,mydimension)); });
409 AffineGeometry<ctype,mydimension,coorddimension> affineGeo(refSimplex,simplexCorners);
410 using std::sqrt;
411 const ctype tol = sqrt(std::numeric_limits<ctype>::epsilon());
412 for (int i = 0; i < corners(); ++i) {
413 const auto corner = refElement_.position(i,mydimension);
414 if ((affineGeo.global(corner) - global(corner)).two_norm() > tol)
415 return false;
416 }
417 return true;
418 }
419 }
420 }
421
422private:
424 ReferenceElement refElement_{};
425
427 LocalFiniteElement localFE_{};
428
430 std::vector<GlobalCoordinate> vertices_{};
431
432 mutable std::optional<bool> affine_ = std::nullopt;
433};
434
435namespace Impl {
436
437// extract the LocalCoordinate type from a LocalFiniteElement
438template <class LFE>
439using LocalCoordinate_t
440 = FieldVector<typename LFE::Traits::LocalBasisType::Traits::DomainFieldType,
441 LFE::Traits::LocalBasisType::Traits::dimDomain>;
442
443} // end namespace Impl
444
445
446// deduction guides
447template <class I, class LFE, class GlobalCoordinate>
448LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, std::vector<GlobalCoordinate>)
449 -> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
450
451template <class I, class LFE, class F,
452 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
453LocalFiniteElementGeometry (Geo::ReferenceElement<I>, const LFE&, const F&)
454 -> LocalFiniteElementGeometry<LFE, Range::dimension>;
455
456template <class LFE, class GlobalCoordinate>
457LocalFiniteElementGeometry (GeometryType, const LFE& localFE, std::vector<GlobalCoordinate>)
458 -> LocalFiniteElementGeometry<LFE, GlobalCoordinate::dimension>;
459
460template <class LFE, class F,
461 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
462LocalFiniteElementGeometry (GeometryType, const LFE&, const F&)
463 -> LocalFiniteElementGeometry<LFE, Range::dimension>;
464
465} // namespace Dune
466
467#endif // DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
An implementation of the Geometry interface for affine geometries.
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
static constexpr int rows
The number of rows.
Definition: fmatrix.hh:123
FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition: fmatrix.hh:171
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition: type.hh:319
Geometry implementation based on local-basis function parametrization.
Definition: localfiniteelementgeometry.hh:40
GeometryType type() const
Obtain the name of the reference element.
Definition: localfiniteelementgeometry.hh:167
decltype(power(std::declval< ctype >(), mydimension)) Volume
type of volume
Definition: localfiniteelementgeometry.hh:62
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, std::vector< GlobalCoordinate > vertices)
Constructor from a vector of coefficients of the LocalBasis parametrizing the geometry.
Definition: localfiniteelementgeometry.hh:103
JacobianTransposed jacobianTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian.
Definition: localfiniteelementgeometry.hh:327
const LocalBasis & localBasis() const
The local basis of the stored local finite-element.
Definition: localfiniteelementgeometry.hh:377
typename LocalBasisTraits::DomainFieldType ctype
coordinate type
Definition: localfiniteelementgeometry.hh:47
JacobianInverseTransposed jacobianInverseTransposed(const LocalCoordinate &local) const
Obtain the transposed of the Jacobian's inverse.
Definition: localfiniteelementgeometry.hh:353
friend ReferenceElement referenceElement(const LocalFiniteElementGeometry &geometry)
Obtain the reference-element related to this geometry.
Definition: localfiniteelementgeometry.hh:359
LocalFiniteElementGeometry()=default
Default constructed geometry results in an empty/invalid representation.
bool affine() const
Is this mapping affine? Geometries of order 1 might be affine, but it needs to be checked on non-simp...
Definition: localfiniteelementgeometry.hh:159
Jacobian jacobian(const LocalCoordinate &local) const
Obtain the Jacobian.
Definition: localfiniteelementgeometry.hh:307
Volume volume(const QuadratureRule< ctype, mydimension > &quadRule) const
Obtain the volume of the mapping's image by given quadrature rules.
Definition: localfiniteelementgeometry.hh:294
const std::vector< GlobalCoordinate > & coefficients() const
Obtain the coefficients of the parametrization.
Definition: localfiniteelementgeometry.hh:371
LocalFiniteElementGeometry(GeometryType gt, Args &&... args)
Constructor, forwarding to the other constructors that take a reference-element.
Definition: localfiniteelementgeometry.hh:144
static const int coorddimension
coordinate dimension
Definition: localfiniteelementgeometry.hh:53
int corners() const
Obtain number of corners of the corresponding reference element.
Definition: localfiniteelementgeometry.hh:173
GlobalCoordinate center() const
Obtain the centroid of the mapping's image.
Definition: localfiniteelementgeometry.hh:186
const LocalFiniteElement & finiteElement() const
Obtain the local finite-element.
Definition: localfiniteelementgeometry.hh:365
GlobalCoordinate global(const LocalCoordinate &local) const
Evaluate the coordinate mapping.
Definition: localfiniteelementgeometry.hh:202
GlobalCoordinate corner(int i) const
Obtain coordinates of the i-th corner.
Definition: localfiniteelementgeometry.hh:179
Volume volume(Impl::ConvergenceOptions< ctype > opts={}) const
Obtain the volume of the mapping's image.
Definition: localfiniteelementgeometry.hh:276
static const int mydimension
geometry dimension
Definition: localfiniteelementgeometry.hh:50
JacobianInverse jacobianInverse(const LocalCoordinate &local) const
Obtain the Jacobian's inverse.
Definition: localfiniteelementgeometry.hh:339
ctype integrationElement(const LocalCoordinate &local) const
Obtain the integration element.
Definition: localfiniteelementgeometry.hh:260
int order() const
Obtain the (highest) polynomial order of the parametrization.
Definition: localfiniteelementgeometry.hh:149
FieldVector< ctype, mydimension > LocalCoordinate
type of local coordinates
Definition: localfiniteelementgeometry.hh:56
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, Param &&parametrization)
Constructor from a local parametrization function, mapping local to (curved) global coordinates.
Definition: localfiniteelementgeometry.hh:128
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
A container for all quadrature rules of dimension dim
Definition: quadraturerules.hh:260
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
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
auto transformedRangeView(R &&range, F &&f)
Create a TransformedRangeView.
Definition: rangeutilities.hh:670
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
STL namespace.
Class providing access to the singletons of the reference elements.
Definition: referenceelements.hh:128
typename Container::ReferenceElement ReferenceElement
The reference element type.
Definition: referenceelements.hh:146
DF DomainFieldType
Export type for domain field.
Definition: localbasis.hh:37
static constexpr int dimDomain
dimension of the domain
Definition: localbasis.hh:40
A unique label for each type of element that can occur in a grid.
Traits for type conversions and type information.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)