│ │ │ │
│ │ │ │
│ │ │ │
│ │ │ │ -
│ │ │ │ -
6#ifndef DUNE_GEOMETRY_QUADRATURERULES_HH
│ │ │ │ -
7#define DUNE_GEOMETRY_QUADRATURERULES_HH
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
16#include <dune/common/fvector.hh>
│ │ │ │ -
17#include <dune/common/exceptions.hh>
│ │ │ │ -
18#include <dune/common/stdstreams.hh>
│ │ │ │ -
19#include <dune/common/stdthread.hh>
│ │ │ │ -
20#include <dune/common/visibility.hh>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
33 template<
typename ct,
int dim>
│ │ │ │ -
34 class QuadraturePoint;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
39 template<
typename ct,
int dim>
│ │ │ │ -
40 struct tuple_size<
Dune::QuadraturePoint<ct,dim>> :
public std::integral_constant<std::size_t,2> {};
│ │ │ │ -
│ │ │ │ -
42 template<
typename ct,
int dim>
│ │ │ │ -
43 struct tuple_element<0,
Dune::QuadraturePoint<ct,dim>> {
using type = Dune::FieldVector<ct, dim>; };
│ │ │ │ +
5#ifndef DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
│ │ │ │ +
6#define DUNE_GEOMETRY_LOCALFINITEELEMENTGEOMETRY_HH
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
14#include <dune/common/fmatrix.hh>
│ │ │ │ +
15#include <dune/common/fvector.hh>
│ │ │ │ +
16#include <dune/common/math.hh>
│ │ │ │ +
17#include <dune/common/typetraits.hh>
│ │ │ │ +
18#include <dune/common/std/type_traits.hh>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
38template <
class LFE,
int cdim>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
41 using LocalFiniteElement = LFE;
│ │ │ │ +
42 using LocalBasis =
typename LFE::Traits::LocalBasisType;
│ │ │ │ +
43 using LocalBasisTraits =
typename LocalBasis::Traits;
│ │ │ │
│ │ │ │ -
45 template<
typename ct,
int dim>
│ │ │ │ -
46 struct tuple_element<1,
Dune::QuadraturePoint<ct,dim>> {
using type = ct; };
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
47 using ctype =
typename LocalBasisTraits::DomainFieldType;
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
65 template<
typename ct,
int dim>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
75 typedef Dune::FieldVector<ct,dim>
Vector;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
65 using Jacobian = FieldMatrix<ctype, coorddimension, mydimension>;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
111 template<std::
size_t index, std::enable_if_t<(index<=1),
int> = 0>
│ │ │ │ -
│ │ │ │ -
112 std::tuple_element_t<index, QuadraturePo
int<ct, dim>> get() const
│ │ │ │ -
│ │ │ │ -
114 if constexpr (index == 0) {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
123 FieldVector<ct, dim> local;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
130 namespace QuadratureType {
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
104 const LocalFiniteElement& localFE,
│ │ │ │ +
105 std::vector<GlobalCoordinate> vertices)
│ │ │ │ +
106 : refElement_(refElement)
│ │ │ │ +
│ │ │ │ +
108 , vertices_(
std::move(vertices))
│ │ │ │ +
│ │ │ │ +
110 assert(localFE_.size() == vertices_.size());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
126 template <
class Param,
│ │ │ │ +
127 std::enable_if_t<std::is_invocable_r_v<GlobalCoordinate,Param,LocalCoordinate>,
int> = 0>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
129 const LocalFiniteElement& localFE,
│ │ │ │ +
130 Param&& parametrization)
│ │ │ │ +
131 : refElement_(refElement)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
134 localFE_.localInterpolation().interpolate(parametrization, vertices_);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
143 template <
class... Args>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
162 affine_.emplace(affineImpl());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
169 return refElement_.type();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
212 template<
typename ct,
int dim>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
231 constexpr static int d = dim;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
237 virtual int order ()
const {
return delivered_order; }
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
245 typedef typename std::vector<QuadraturePoint<ct,dim> >::const_iterator
iterator;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
181 assert( (i >= 0) && (i <
corners()) );
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
188 return global(refElement_.position(0, 0));
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
204 thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
│ │ │ │ +
│ │ │ │ +
206 assert(shapeValues.size() == vertices_.size());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
209 for (std::size_t i = 0; i < shapeValues.size(); ++i)
│ │ │ │ +
210 out.axpy(shapeValues[i], vertices_[i]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
235 Impl::GaussNewtonErrorCode err = Impl::gaussNewton(
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
241 if (err != Impl::GaussNewtonErrorCode::OK)
│ │ │ │ +
242 DUNE_THROW(Dune::Exception,
│ │ │ │ +
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);
│ │ │ │
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
259 template<
typename ctype,
int dim>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │
│ │ │ │ -
│ │ │ │ -
266 using QuadratureOrderVector = std::vector<std::pair<std::once_flag, QuadratureRule> >;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
269 using GeometryTypeVector = std::vector<std::pair<std::once_flag, QuadratureOrderVector> >;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
272 using QuadratureCacheVector = std::vector<std::pair<std::once_flag, GeometryTypeVector> >;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
277 assert(t.
dim()==dim);
│ │ │ │ -
│ │ │ │ -
279 DUNE_ASSERT_CALL_ONCE();
│ │ │ │ -
│ │ │ │ -
281 static QuadratureCacheVector quadratureCache(QuadratureType::size);
│ │ │ │ -
│ │ │ │ -
283 auto& [ onceFlagQuadratureType, geometryTypes ] = quadratureCache[qt];
│ │ │ │ -
│ │ │ │ -
285 std::call_once(onceFlagQuadratureType, [&types = geometryTypes]{
│ │ │ │ -
286 types = GeometryTypeVector(LocalGeometryTypeIndex::size(dim));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
289 auto& [ onceFlagGeometryType, quadratureOrders ] = geometryTypes[LocalGeometryTypeIndex::index(t)];
│ │ │ │ -
│ │ │ │ -
291 std::call_once(onceFlagGeometryType, [&, &orders = quadratureOrders]{
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
294 orders = QuadratureOrderVector(numRules);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
298 auto& [ onceFlagQuadratureOrder, quadratureRule ] = quadratureOrders[dim == 0 ? 0 : p];
│ │ │ │ -
│ │ │ │ -
300 std::call_once(onceFlagQuadratureOrder, [&, &rule = quadratureRule]{
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
304 return quadratureRule;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
328 return instance()._rule(t,p,qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
335 return instance()._rule(gt,p,qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
341#define DUNE_INCLUDING_IMPLEMENTATION
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
344#include "quadraturerules/pointquadrature.hh"
│ │ │ │ -
│ │ │ │ -
346#include "quadraturerules/gausslobattoquadrature.hh"
│ │ │ │ -
347#include "quadraturerules/gaussquadrature.hh"
│ │ │ │ -
348#include "quadraturerules/gaussradauleftquadrature.hh"
│ │ │ │ -
349#include "quadraturerules/gaussradaurightquadrature.hh"
│ │ │ │ -
350#include "quadraturerules/jacobi1quadrature.hh"
│ │ │ │ -
351#include "quadraturerules/jacobi2quadrature.hh"
│ │ │ │ -
352#include "quadraturerules/jacobiNquadrature.hh"
│ │ │ │ -
│ │ │ │ -
354#include "quadraturerules/prismquadrature.hh"
│ │ │ │ -
│ │ │ │ -
356#include "quadraturerules/simplexquadrature.hh"
│ │ │ │ -
357#include "quadraturerules/tensorproductquadrature.hh"
│ │ │ │ -
│ │ │ │ -
359#undef DUNE_INCLUDING_IMPLEMENTATION
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
369 template<
typename ctype,
int dim>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
375 return TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
379 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
383 template<
typename ctype>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
386 constexpr static int dim = 0;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
392 return std::numeric_limits<int>::max();
│ │ │ │ -
│ │ │ │ -
394 DUNE_THROW(Exception,
"Unknown GeometryType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
400 return PointQuadratureRule<ctype>();
│ │ │ │ -
│ │ │ │ -
402 DUNE_THROW(Exception,
"Unknown GeometryType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
406 template<
typename ctype>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
409 constexpr static int dim = 1;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
417 return GaussQuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
419 return Jacobi1QuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
421 return Jacobi2QuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
423 return GaussLobattoQuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
425 return JacobiNQuadratureRule1D<ctype>::maxOrder();
│ │ │ │ -
│ │ │ │ -
427 return GaussRadauLeftQuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
429 return GaussRadauRightQuadratureRule1D<ctype>::highest_order;
│ │ │ │ -
│ │ │ │ -
431 DUNE_THROW(Exception,
"Unknown QuadratureType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
434 DUNE_THROW(Exception,
"Unknown GeometryType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
442 return GaussQuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
444 return Jacobi1QuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
446 return Jacobi2QuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
448 return GaussLobattoQuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
450 return JacobiNQuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
452 return GaussRadauLeftQuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
454 return GaussRadauRightQuadratureRule1D<ctype>(p);
│ │ │ │ -
│ │ │ │ -
456 DUNE_THROW(Exception,
"Unknown QuadratureType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
459 DUNE_THROW(Exception,
"Unknown GeometryType");
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
463 template<
typename ctype>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
466 constexpr static int dim = 2;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
471 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
474 (order,
static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
481 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
│ │ │ │ -
│ │ │ │ -
483 return SimplexQuadratureRule<ctype,dim>(p);
│ │ │ │ -
│ │ │ │ -
485 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
489 template<
typename ctype>
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
492 constexpr static int dim = 3;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
497 TensorProductQuadratureRule<ctype,dim>::maxOrder(t.
id(), qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
500 (order,
static_cast<unsigned>(SimplexQuadratureRule<ctype,dim>::highest_order));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
503 (order,
static_cast<unsigned>(PrismQuadratureRule<ctype,dim>::highest_order));
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
511 && p <= SimplexQuadratureRule<ctype,dim>::highest_order)
│ │ │ │ -
│ │ │ │ -
513 return SimplexQuadratureRule<ctype,dim>(p);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
517 && p <= PrismQuadratureRule<ctype,dim>::highest_order)
│ │ │ │ -
│ │ │ │ -
519 return PrismQuadratureRule<ctype,dim>(p);
│ │ │ │ -
│ │ │ │ -
521 return TensorProductQuadratureRule<ctype,dim>(t.
id(), p, qt);
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
525#ifndef DUNE_NO_EXTERN_QUADRATURERULES
│ │ │ │ -
526 extern template class GaussLobattoQuadratureRule<double, 1>;
│ │ │ │ -
527 extern template class GaussQuadratureRule<double, 1>;
│ │ │ │ -
528 extern template class GaussRadauLeftQuadratureRule<double, 1>;
│ │ │ │ -
529 extern template class GaussRadauRightQuadratureRule<double, 1>;
│ │ │ │ -
530 extern template class Jacobi1QuadratureRule<double, 1>;
│ │ │ │ -
531 extern template class Jacobi2QuadratureRule<double, 1>;
│ │ │ │ -
532 extern template class JacobiNQuadratureRule<double, 1>;
│ │ │ │ -
533 extern template class PrismQuadratureRule<double, 3>;
│ │ │ │ -
534 extern template class SimplexQuadratureRule<double, 2>;
│ │ │ │ -
535 extern template class SimplexQuadratureRule<double, 3>;
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
│ │ │ │ -
A unique label for each type of element that can occur in a grid.
│ │ │ │ -
Helper classes to provide indices for geometrytypes for use in a vector.
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
283 for (
int p = 2; p < opts.maxIt; ++p) {
│ │ │ │ +
│ │ │ │ +
285 if (abs(vol1 - vol0) < opts.absTol)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
297 for (
const auto& qp : quadRule)
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
309 thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
│ │ │ │ +
│ │ │ │ +
311 assert(shapeJacobians.size() == vertices_.size());
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
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]);
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
361 return geometry.refElement_;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
379 return localFE_.localBasis();
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
384 bool affineImpl ()
const
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
405 auto refSimplex = referenceElement<ctype,mydimension>(GeometryTypes::simplex(
mydimension));
│ │ │ │ +
406 auto simplexIndices = Dune::range(refSimplex.size(
mydimension));
│ │ │ │ +
407 auto simplexCorners = Dune::transformedRangeView(simplexIndices,
│ │ │ │ +
│ │ │ │ +
409 AffineGeometry<ctype,mydimension,coorddimension> affineGeo(refSimplex,simplexCorners);
│ │ │ │ +
│ │ │ │ +
411 const ctype tol = sqrt(std::numeric_limits<ctype>::epsilon());
│ │ │ │ +
412 for (
int i = 0; i <
corners(); ++i) {
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
427 LocalFiniteElement localFE_{};
│ │ │ │ +
│ │ │ │ +
430 std::vector<GlobalCoordinate> vertices_{};
│ │ │ │ +
│ │ │ │ +
432 mutable std::optional<bool> affine_ = std::nullopt;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
439using LocalCoordinate_t
│ │ │ │ +
440 = FieldVector<
typename LFE::Traits::LocalBasisType::Traits::DomainFieldType,
│ │ │ │ +
441 LFE::Traits::LocalBasisType::Traits::dimDomain>;
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
447template <
class I,
class LFE,
class GlobalCoordinate>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
451template <
class I,
class LFE,
class F,
│ │ │ │ +
452 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
456template <
class LFE,
class GlobalCoordinate>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
460template <
class LFE,
class F,
│ │ │ │ +
461 class Range = std::invoke_result_t<F,Impl::LocalCoordinate_t<LFE>>>
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
An implementation of the Geometry interface for affine geometries.
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
│ │ │ │ +
A unique label for each type of element that can occur in a grid.
│ │ │ │
│ │ │ │
Definition affinegeometry.hh:21
│ │ │ │ -
Enum
Definition quadraturerules.hh:131
│ │ │ │ -
@ GaussJacobi_n_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:168
│ │ │ │ -
@ GaussJacobi_2_0
Gauss-Legendre rules with .
Definition quadraturerules.hh:155
│ │ │ │ -
@ GaussRadauRight
Gauss-Radau rules including the right endpoint.
Definition quadraturerules.hh:193
│ │ │ │ -
@ GaussJacobi_1_0
Gauss-Jacobi rules with .
Definition quadraturerules.hh:148
│ │ │ │ -
@ GaussLobatto
Gauss-Lobatto rules.
Definition quadraturerules.hh:176
│ │ │ │ -
@ GaussRadauLeft
Gauss-Radau rules including the left endpoint.
Definition quadraturerules.hh:184
│ │ │ │ -
@ GaussLegendre
Gauss-Legendre rules (default)
Definition quadraturerules.hh:141
│ │ │ │ -
Single evaluation point in a quadrature rule.
Definition quadraturerules.hh:66
│ │ │ │ -
const Vector & position() const
return local coordinates of integration point i
Definition quadraturerules.hh:82
│ │ │ │ -
Dune::FieldVector< ct, dim > Vector
Type used for the position of a quadrature point.
Definition quadraturerules.hh:75
│ │ │ │ -
ct Field
Number type used for coordinates and quadrature weights.
Definition quadraturerules.hh:72
│ │ │ │ -
const ct & weight() const
return weight associated with integration point i
Definition quadraturerules.hh:88
│ │ │ │ -
ct weight_
Definition quadraturerules.hh:124
│ │ │ │ -
static constexpr int dimension
Dimension of the integration domain.
Definition quadraturerules.hh:69
│ │ │ │ -
QuadraturePoint(const Vector &x, ct w)
set up quadrature of given order in d dimensions
Definition quadraturerules.hh:78
│ │ │ │ -
FieldVector< ct, dim > local
Definition quadraturerules.hh:123
│ │ │ │ -
Dune::FieldVector< ct, dim > type
Definition quadraturerules.hh:43
│ │ │ │ -
ct type
Definition quadraturerules.hh:46
│ │ │ │ -
Exception thrown if a desired QuadratureRule is not available, because the requested order is to high...
Definition quadraturerules.hh:55
│ │ │ │ +
This class provides access to geometric and topological properties of a reference element.
Definition referenceelement.hh:52
│ │ │ │ +
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
│ │ │ │ +
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
│ │ │ │ +
LocalCoordinate local(const GlobalCoordinate &y, Impl::GaussNewtonOptions< ctype > opts={}) const
Evaluate the inverse coordinate mapping.
Definition localfiniteelementgeometry.hh:232
│ │ │ │ +
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
│ │ │ │ +
Impl::FieldMatrixHelper< ctype > MatrixHelper
Definition localfiniteelementgeometry.hh:82
│ │ │ │ +
typename ReferenceElements::ReferenceElement ReferenceElement
Definition localfiniteelementgeometry.hh:79
│ │ │ │ +
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
│ │ │ │ +
FieldVector< ctype, coorddimension > GlobalCoordinate
type of global coordinates
Definition localfiniteelementgeometry.hh:59
│ │ │ │ +
FieldMatrix< ctype, coorddimension, mydimension > Jacobian
type of jacobian
Definition localfiniteelementgeometry.hh:65
│ │ │ │ +
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
type of jacobian inverse transposed
Definition localfiniteelementgeometry.hh:74
│ │ │ │ +
FieldMatrix< ctype, mydimension, coorddimension > JacobianInverse
type of jacobian inverse
Definition localfiniteelementgeometry.hh:71
│ │ │ │ +
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
│ │ │ │ +
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
type of jacobian transposed
Definition localfiniteelementgeometry.hh:68
│ │ │ │ +
LocalFiniteElementGeometry(const ReferenceElement &refElement, const LocalFiniteElement &localFE, Param &¶metrization)
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
│ │ │ │ -
virtual ~QuadratureRule()
Definition quadraturerules.hh:241
│ │ │ │ -
virtual GeometryType type() const
return type of element
Definition quadraturerules.hh:240
│ │ │ │ -
int delivered_order
Definition quadraturerules.hh:249
│ │ │ │ -
QuadratureRule(GeometryType t, int order)
Constructor for a given geometry type and a given quadrature order.
Definition quadraturerules.hh:228
│ │ │ │ -
GeometryType geometry_type
Definition quadraturerules.hh:248
│ │ │ │ -
ct CoordType
The type used for coordinates.
Definition quadraturerules.hh:234
│ │ │ │ -
QuadratureRule()
Default constructor.
Definition quadraturerules.hh:221
│ │ │ │ -
virtual int order() const
return order
Definition quadraturerules.hh:237
│ │ │ │ -
QuadratureRule(GeometryType t)
Constructor for a given geometry type. Leaves the quadrature order invalid
Definition quadraturerules.hh:225
│ │ │ │ -
std::vector< QuadraturePoint< ct, dim > >::const_iterator iterator
Definition quadraturerules.hh:245
│ │ │ │ -
Factory class for creation of quadrature rules, depending on GeometryType, order and QuadratureType.
Definition quadraturerules.hh:370
│ │ │ │
A container for all quadrature rules of dimension dim
Definition quadraturerules.hh:260
│ │ │ │ -
static unsigned maxOrder(const GeometryType &t, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
maximum quadrature order for given geometry type and quadrature type
Definition quadraturerules.hh:319
│ │ │ │ -
static const QuadratureRule & rule(const GeometryType::BasicType t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition quadraturerules.hh:332
│ │ │ │ -
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition quadraturerules.hh:326
│ │ │ │
Unique label for each type of entities that can occur in DUNE grids.
Definition type.hh:114
│ │ │ │ -
constexpr bool isPrism() const
Return true if entity is a prism.
Definition type.hh:309
│ │ │ │ -
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition type.hh:279
│ │ │ │ -
constexpr unsigned int dim() const
Return dimension of the type.
Definition type.hh:360
│ │ │ │ -
BasicType
Each entity can be tagged by one of these basic types plus its space dimension.
Definition type.hh:120
│ │ │ │ -
constexpr bool isLine() const
Return true if entity is a line segment.
Definition type.hh:284
│ │ │ │ -
constexpr unsigned int id() const
Return the topology id of the type.
Definition type.hh:365
│ │ │ │
constexpr bool isSimplex() const
Return true if entity is a simplex of any dimension.
Definition type.hh:319
│ │ │ │