Loading...
Searching...
No Matches
argyrisbasis.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
#ifndef DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ARGYRISBASIS_HH
#define DUNE_FUNCTIONS_FUNCTIONSPACEBASES_ARGYRISBASIS_HH
#include <algorithm>
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localkey.hh>
#include <dune/functions/analyticfunctions/monomialset.hh>
#include <dune/functions/common/mapperutilities.hh>
#include <dune/functions/functionspacebases/cubichermitebasis.hh>
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/functionaldescriptor.hh>
#include <dune/functions/functionspacebases/leafprebasismappermixin.hh>
#include <dune/functions/functionspacebases/nodes.hh>
#include <dune/functions/functionspacebases/transformedfiniteelementmixin.hh>
namespace Dune::Functions
{
template<class GV, class R>
class ArgyrisPreBasis;
template<class GV, class R = double>
using ArgyrisBasis = DefaultGlobalBasis<ArgyrisPreBasis<GV, R> >;
namespace Impl
{
class ArgyrisLocalCoefficients
{
public:
using size_type = std::size_t;
static constexpr int dim = 2;
ArgyrisLocalCoefficients()
{
for (unsigned int i = 0; i < 3; i++) // subentities: three vertices
for (unsigned int k = 0; k < 6; k++) // 6 basis functions per vertex
localKeys_[6 * i + k]
= LocalKey(i, dim, k); //(subentity, codim, number of dof for this subentity)
for (unsigned int i = 0; i < 3; ++i) // subentities: three edges
localKeys_[18 + i] = LocalKey(i, dim - 1, 0); // one node per edge
}
static constexpr size_type size()
{
return 21;
}
LocalKey const &localKey(size_type i) const
{
return localKeys_[i];
}
private:
std::array<Dune::LocalKey, 21> localKeys_;
};
template<class D, class R>
class ArgyrisReferenceLocalBasis
{
static constexpr int dim = 2;
public:
using Traits = H2LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1, FieldVector<R, 1>,
FieldMatrix<R, 1, dim>, FieldMatrix<R, dim, dim> >;
private:
static constexpr auto getArgyrisCoefficients()
{
// Define std::sqrt(2.) manually in double precision,
// because std::sqrt is not constexpr before C++26.
D sqrt2 = -8. * 1.414213562373095;
return Dune::FieldMatrix<D, 21,21>{
// vertex functionals
// l_0
{1, /*0th order*/
0, 0, /*1th order*/
0, 0, 0, /*2th order*/
-10, 0, 0, -10, /*3th order*/
15, 0, -30, 0, 15, /*4th order*/
-6, 0, 30, 30, 0, -6}, /*5th order*/
// l_1
{0,
1, 0,
0, 0, 0,
-6, 0, -11, 0,
8, 0, 10, 18, 0,
-3, 0, 1, -10, -8, 0},
// l_2, l_1 mirrored
{
0,
0, 1,
0, 0, 0,
0, -11, 0, -6,
0, 18, 10, 0, 8,
0, -8, -10, 1, 0, -3},
// l_3
{
0,
0, 0,
0.5, 0, 0,
-1.5, 0, 0, 0,
1.5, 0, -1.5, 0, 0,
-0.5, 0, 1.5, 1, 0, 0},
// l_4
{
0,
0, 0,
0, 1, 0,
0, -4, -4, 0,
0, 5, 10, 5, 0,
0, -2, -6, -6, -2, 0},
// l_5, l_3 mirrored
{
0,
0, 0,
0, 0, 0.5,
0, 0, 0, -1.5,
0, 0, -1.5, 0, 1.5,
0, 0, 1, 1.5, 0, -0.5},
// l_6
{0,
0, 0,
0, 0, 0,
10, 0, 0, 0,
-15, 0, 15, 0, 0,
6, 0, -15, -15, 0, 0},
// l_7
{0,
0, 0,
0, 0, 0,
-4, 0, 0, 0,
7, 0, -3.5, 0, 0,
-3, 0, 3.5, 3.5, 0, 0},
// l_8
{0,
0, 0,
0, 0, 0,
0, -5, 0, 0,
0, 14, 18.5, 0, 0,
0, -8, -18.5, -13.5, 0, 0},
// l_9
{
0,
0, 0,
0, 0, 0,
0.5, 0, 0, 0,
-1, 0, 0.25, 0, 0,
0.5, 0, -0.25, -0.25, 0, 0},
// l_10
{
0,
0, 0,
0, 0, 0,
0, 1, 0, 0,
0, -3, -3.5, 0, 0,
0, 2, 3.5, 2.5, 0, 0},
// l_11
{
0,
0, 0,
0, 0, 0,
0, 0, 0, 0,
0, 0, 1.25, 0, 0,
0, 0, -0.75, -1.25, 0, 0},
// l_12 mirrors l_6
{
0,
0, 0,
0, 0, 0,
0, 0, 0, 10,
0, 0, 15, 0, -15,
0, 0, -15, -15, 0, 6},
// l_13 mirrors l_8
{
0,
0, 0,
0, 0, 0,
0, 0, -5, 0,
0, 0, 18.5, 14, 0,
0, 0, -13.5, -18.5, -8, 0},
// l_14 mirrors l_7
{
0,
0, 0,
0, 0, 0,
0, 0, 0, -4,
0, 0, -3.5, 0, 7,
0, 0, 3.5, 3.5, 0, -3},
// l_15 mirrors l_11
{
0,
0, 0,
0, 0, 0,
0, 0, 0, 0,
0, 0, 1.25, 0, 0,
0, 0, -1.25, -0.75, 0, 0},
// l_16 mirrors l_10
{
0,
0, 0,
0, 0, 0,
0, 0, 1, 0,
0, 0, -3.5, -3, 0,
0, 0, 2.5, 3.5, 2, 0},
// l_17 mirrors l_9
{
0,
0, 0,
0, 0, 0,
0, 0, 0, 0.5,
0, 0, 0.25, 0, -1,
0, 0, -0.25, -0.25, 0, 0.5},
// edge functionals
// l_18
{
0,
0, 0,
0, 0, 0,
0, 16, 0, 0,
0, -32, -32, 0, 0,
0, 16, 32, 16, 0, 0},
// l_19
{
0,
0, 0,
0, 0, 0,
0, 0, -16, 0,
0, 0, 32, 32, 0,
0, 0, -16, -32, -16, 0},
// l_20
{
0,
0, 0,
0, 0, 0,
0, 0, 0, 0,
0, 0, -1. * sqrt2, 0, 0,
0, 0, sqrt2, sqrt2, 0, 0}};
}
static constexpr auto referenceBasisCoefficients = getArgyrisCoefficients();
static constexpr MonomialSet<typename Traits::RangeFieldType, dim, 5> monomials = {};
public:
static constexpr unsigned int size()
{
return ArgyrisLocalCoefficients::size();
}
static constexpr unsigned int order()
{
return 5;
}
void evaluateFunction(const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
out.resize(size());
auto monomialValues = monomials(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void evaluateJacobian(const typename Traits::DomainType &in,
std::vector<typename Traits::JacobianType> &out) const
{
out.resize(size());
auto monomialValues = derivative(monomials)(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void evaluateHessian(const typename Traits::DomainType &in,
std::vector<typename Traits::HessianType> &out) const
{
out.resize(size());
auto monomialValues = derivative(derivative(monomials))(in);
multiplyWithCoefficentMatrix(referenceBasisCoefficients, monomialValues, out);
}
void partial(std::array<unsigned int, dim> order, const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
out.resize(size());
auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
if (totalOrder == 0)
evaluateFunction(in, out);
else if (totalOrder == 1)
{
evaluateJacobian(in,jacobiansBuffer_);
std::size_t which = std::max_element(order.begin(), order.end()) - order.begin();
for (auto i : Dune::range(size()))
out[i] = jacobiansBuffer_[i][0][which];
}
else if (totalOrder == 2)
{
evaluateHessian(in, hessianBuffer_);
std::size_t first, second;
first = std::max_element(order.begin(), order.end()) - order.begin();
if (order[first] == 2)
second = first;
else
{
order[first] = 0;
second = std::max_element(order.begin(), order.end()) - order.begin();
}
for (auto i : Dune::range(size()))
out[i] = hessianBuffer_[i][first][second];
}
else
DUNE_THROW(RangeError, "partial() not implemented for given order");
}
private:
mutable std::vector<typename Traits::JacobianType> jacobiansBuffer_;
mutable std::vector<typename Traits::HessianType> hessianBuffer_;
};
template<class D>
class ArgyrisLocalInterpolation
{
using size_type = std::size_t;
static constexpr int dim = 2;
static constexpr unsigned int size()
{
return ArgyrisLocalCoefficients::size();
}
using FunctionalDescriptor = Dune::Functions::Impl::FunctionalDescriptor<dim>;
public:
ArgyrisLocalInterpolation()
{
// first vertex
// function evaluation
descriptors_[0] = FunctionalDescriptor();
// gradient evaluations
descriptors_[1] = FunctionalDescriptor({1,0});
descriptors_[2] = FunctionalDescriptor({0,1});
// hessian evaluations
descriptors_[3] = FunctionalDescriptor({2,0});
descriptors_[4] = FunctionalDescriptor({1,1});
descriptors_[5] = FunctionalDescriptor({0,2});
// second vertex
descriptors_[6] = FunctionalDescriptor();
descriptors_[7] = FunctionalDescriptor({1,0});
descriptors_[8] = FunctionalDescriptor({0,1});
descriptors_[9] = FunctionalDescriptor({2,0});
descriptors_[10] = FunctionalDescriptor({1,1});
descriptors_[11] = FunctionalDescriptor({0,2});
// third vertex
descriptors_[12] = FunctionalDescriptor();
descriptors_[13] = FunctionalDescriptor({1,0});
descriptors_[14] = FunctionalDescriptor({0,1});
descriptors_[15] = FunctionalDescriptor({2,0});
descriptors_[16] = FunctionalDescriptor({1,1});
descriptors_[17] = FunctionalDescriptor({0,2});
// normal derivatives at edge midpoints
descriptors_[18] = FunctionalDescriptor(1);
descriptors_[19] = FunctionalDescriptor(1);
descriptors_[20] = FunctionalDescriptor(1);
}
template<class Element>
void bind( Element const& element, std::array<D, 3>const& averageVertexMeshSize, std::bitset<3>const& edgeOrientation)
{
averageVertexMeshSize_ = &averageVertexMeshSize;
auto geometry = element.geometry();
// get global Normals and midpoints
auto refElement = Dune::referenceElement<D, 2>(geometry.type());
for (std::size_t i = 0; i < 3; ++i)
{
localVertices_[i] = refElement.position(i, 2);
localMidpoints_[i] = refElement.position(i, 1);
std::size_t lower = (i == 2) ? 1 : 0;
std::size_t upper = (i == 0) ? 1 : 2;
auto edge = geometry.global(refElement.position(upper, 2))
- geometry.global(refElement.position(lower, 2));
// normalize and orient
edge /= edge.two_norm() * (edgeOrientation[i] ? -1. : 1.);
// Rotation by pi/2. Note that Kirby rotates by 3*pi/2
globalNormals_[i] = {-edge[1], edge[0]};
}
}
template<class F, class C>
void interpolate(F const& f, std::vector<C>& out) const
{
out.resize(size());
auto&& df = derivative(f);
auto&& Hf = derivative(df);
int offset = 0;
// iterate over vertices
for (int i = 0; i < 3; i++, offset += 6)
{
auto&& x = localVertices_[i];
auto derivativeValue = squeezeTensor(df(x));
auto hessianValue = squeezeTensor(Hf(x));
auto && h = (*averageVertexMeshSize_)[i];
out[offset ] = f(x);
out[offset + 1] = derivativeValue[0] * h;
out[offset + 2] = derivativeValue[1] * h;
out[offset + 3] = hessianValue[0][0] * h*h;
out[offset + 4] = hessianValue[0][1] * h*h;
out[offset + 5] = hessianValue[1][1] * h*h;
}
// iterate over edges
for (int i = 0; i < 3; i++)
out[18 + i] = squeezeTensor(df(localMidpoints_[i])).dot(globalNormals_[i]);
}
const FunctionalDescriptor& functionalDescriptor(size_type i) const
{
return descriptors_[i];
}
protected:
std::array<Dune::FieldVector<D, 2>, 3> globalNormals_;
std::array<Dune::FieldVector<D, 2>, 3> localMidpoints_;
std::array<Dune::FieldVector<D, 2>, 3> localVertices_;
std::array<D, 3> const* averageVertexMeshSize_;
std::array<FunctionalDescriptor, size()> descriptors_;
};
// TODO: There is a copy of this further up
template<class D, class R>
struct ArgyrisLocalBasisTraits
: public H2LocalBasisTraits<D, 2, Dune::FieldVector<D,2>, R, 1,
Dune::FieldVector<R,1>, Dune::FieldMatrix<R,1,2>, Dune::FieldMatrix<R,2,2> >
{};
template<class D, class R>
class ArgyrisLocalFiniteElement
: public Impl::TransformedFiniteElementMixin<ArgyrisLocalFiniteElement<D,R>, ArgyrisLocalBasisTraits<D, R> >
{
using Base = Impl::TransformedFiniteElementMixin< ArgyrisLocalFiniteElement<D,R>, ArgyrisLocalBasisTraits<D, R> >;
friend class Impl::TransformedLocalBasis<ArgyrisLocalFiniteElement<D,R>, ArgyrisLocalBasisTraits<D, R> >;
static constexpr int dim = 2;
public:
using size_type = std::size_t;
using Traits = LocalFiniteElementTraits<
Impl::TransformedLocalBasis<ArgyrisLocalFiniteElement<D,R>, ArgyrisLocalBasisTraits<D, R> >,
Impl::ArgyrisLocalCoefficients,
Impl::ArgyrisLocalInterpolation<D> >;
const typename Traits::LocalCoefficientsType &localCoefficients() const
{
return coefficients_;
}
const typename Traits::LocalInterpolationType &localInterpolation() const
{
return interpolation_;
}
static constexpr GeometryType type()
{
return GeometryTypes::simplex(dim);
}
static constexpr size_type size()
{
return Impl::ArgyrisLocalCoefficients::size();
}
template<class VertexMapper,class ElementMapper, class Element>
void bind(VertexMapper const& vertexMapper, std::vector<D> const& globalAverageVertexMeshSize,
ElementMapper const& elementMapper, std::vector<std::bitset<3> >const& edgeOrientations, Element const &e)
{
// Cache average mesh size for each vertex
for (auto i : range(dim+1))
averageVertexMeshSize_[i] = globalAverageVertexMeshSize[vertexMapper.subIndex(e, i, dim)];
// Cache orientation for each mesh
edgeOrientation_ = edgeOrientations[elementMapper.index(e)];
// Bind LocalInterpolation to updated local state
interpolation_.bind(e, averageVertexMeshSize_, edgeOrientation_);
// Compute local transformation matrices for each vertex
const auto& geometry = e.geometry();
const auto& refElement = ReferenceElements<typename Element::Geometry::ctype, dim>::simplex();
for (auto i : range(dim+1))
{
vertexJacobians_[i] = geometry.jacobian(refElement.position(i, dim));
}
// get geometrical information
std::array<FieldVector<R, 2>, 3> referenceTangents;
// get local and global Tangents
for (std::size_t i = 0; i < 3; ++i)
{
std::size_t lower = (i == 2) ? 1 : 0;
std::size_t upper = (i == 0) ? 1 : 2;
auto edge = refElement.position(upper, 2) - refElement.position(lower, 2);
// store normalized reference Tangent vectors
referenceTangents[i] = edge / edge.two_norm();
auto globalEdge = geometry.global(refElement.position(upper, 2))
- geometry.global(refElement.position(lower, 2));
// store length of global tangents and normalized global tangent vectors
l[i] = globalEdge.two_norm();
globalTangents[i] = globalEdge / l[i];
tau[i] = FieldVector<R, 3>{
globalTangents[i][0] * globalTangents[i][0],
2. * globalTangents[i][0] * globalTangents[i][1],
globalTangents[i][1] * globalTangents[i][1]
};
// The following variables are named as in Kirby's paper
// two g matrices
auto&& referenceG = FieldMatrix<R, 2, 2>{
{-referenceTangents[i][1], referenceTangents[i][0]},
{referenceTangents[i][0], referenceTangents[i][1]}
};
auto&& globalG = FieldMatrix<R, 2, 2>{
{-globalTangents[i][1], globalTangents[i][0]},
{globalTangents[i][0], globalTangents[i][1]}
};
//
b[i] = globalG * geometry.jacobian(refElement.position(i, 2)) * referenceG;
theta[i][0][0] = vertexJacobians_[i][0][0] * vertexJacobians_[i][0][0];
theta[i][0][1] = vertexJacobians_[i][0][0] * vertexJacobians_[i][0][1];
theta[i][0][2] = vertexJacobians_[i][0][1] * vertexJacobians_[i][0][1];
theta[i][1][0] = 2. * vertexJacobians_[i][0][0] * vertexJacobians_[i][1][0];
theta[i][1][1] = vertexJacobians_[i][0][0] * vertexJacobians_[i][1][1] + vertexJacobians_[i][0][1] * vertexJacobians_[i][1][0];
theta[i][1][2] = 2. * vertexJacobians_[i][0][1] * vertexJacobians_[i][1][1];
theta[i][2][0] = vertexJacobians_[i][1][0] * vertexJacobians_[i][1][0];
theta[i][2][1] = vertexJacobians_[i][1][0] * vertexJacobians_[i][1][1];
theta[i][2][2] = vertexJacobians_[i][1][1] * vertexJacobians_[i][1][1];
}
}
protected:
Impl::ArgyrisReferenceLocalBasis<D, R> const& referenceLocalBasis() const
{
return basis_;
}
template<class InputValues, class OutputValues>
void transform(InputValues const &inValues, OutputValues &outValues) const
{
using std::pow;
assert(inValues.size() == size());
assert(outValues.size() == inValues.size());
// compatibility with sympy code below
auto &[b_0, b_1, b_2] = b;
auto &[J_0, J_1, J_2] = vertexJacobians_;
auto &[theta_0, theta_1, theta_2] = theta;
auto & h = averageVertexMeshSize_;
auto & o = edgeOrientation_;
std::array<Dune::FieldVector<R, 2>, 3> const &t = globalTangents;
// This code is generated with sympy.
// It a matrix free implementation of the matrix V in Kirbys paper.
outValues[0] = -15.0/8.0*b_0[1][0]*inValues[18]/l[0] - 15.0/8.0*b_1[1][0]*inValues[19]/l[1] + inValues[0];
outValues[1] = (J_0[0][0]*inValues[1] + J_0[0][1]*inValues[2] - 0.4375*b_0[1][0]*inValues[18]*t[0][0] - 0.4375*b_1[1][0]*inValues[19]*t[1][0])/h[0];
outValues[2] = (J_0[1][0]*inValues[1] + J_0[1][1]*inValues[2] - 0.4375*b_0[1][0]*inValues[18]*t[0][1] - 0.4375*b_1[1][0]*inValues[19]*t[1][1])/h[0];
outValues[3] = (-1.0/32.0*b_0[1][0]*inValues[18]*l[0]*tau[0][0] - 1.0/32.0*b_1[1][0]*inValues[19]*l[1]*tau[1][0] + inValues[3]*theta_0[0][0] + inValues[4]*theta_0[0][1] + inValues[5]*theta_0[0][2])/h[0] /h[0];
outValues[4] = (-1.0/32.0*b_0[1][0]*inValues[18]*l[0]*tau[0][1] - 1.0/32.0*b_1[1][0]*inValues[19]*l[1]*tau[1][1] + inValues[3]*theta_0[1][0] + inValues[4]*theta_0[1][1] + inValues[5]*theta_0[1][2])/h[0] /h[0];
outValues[5] = (-1.0/32.0*b_0[1][0]*inValues[18]*l[0]*tau[0][2] - 1.0/32.0*b_1[1][0]*inValues[19]*l[1]*tau[1][2] + inValues[3]*theta_0[2][0] + inValues[4]*theta_0[2][1] + inValues[5]*theta_0[2][2])/h[0] /h[0];
outValues[6] = (15.0/8.0)*b_0[1][0]*inValues[18]/l[0] - 15.0/8.0*b_2[1][0]*inValues[20]/l[2] + inValues[6];
outValues[7] = (J_1[0][0]*inValues[7] + J_1[0][1]*inValues[8] - 0.4375*b_0[1][0]*inValues[18]*t[0][0] - 0.4375*b_2[1][0]*inValues[20]*t[2][0])/h[1];
outValues[8] = (J_1[1][0]*inValues[7] + J_1[1][1]*inValues[8] - 0.4375*b_0[1][0]*inValues[18]*t[0][1] - 0.4375*b_2[1][0]*inValues[20]*t[2][1])/h[1];
outValues[9] = ((1.0/32.0)*b_0[1][0]*inValues[18]*l[0]*tau[0][0] - 1.0/32.0*b_2[1][0]*inValues[20]*l[2]*tau[2][0] + inValues[9]*theta_1[0][0] + inValues[10]*theta_1[0][1] + inValues[11]*theta_1[0][2])/h[1] /h[1];
outValues[10] = ((1.0/32.0)*b_0[1][0]*inValues[18]*l[0]*tau[0][1] - 1.0/32.0*b_2[1][0]*inValues[20]*l[2]*tau[2][1] + inValues[9]*theta_1[1][0] + inValues[10]*theta_1[1][1] + inValues[11]*theta_1[1][2])/h[1] /h[1];
outValues[11] = ((1.0/32.0)*b_0[1][0]*inValues[18]*l[0]*tau[0][2] - 1.0/32.0*b_2[1][0]*inValues[20]*l[2]*tau[2][2] + inValues[9]*theta_1[2][0] + inValues[10]*theta_1[2][1] + inValues[11]*theta_1[2][2])/h[1] /h[1];
outValues[12] = (15.0/8.0)*b_1[1][0]*inValues[19]/l[1] + (15.0/8.0)*b_2[1][0]*inValues[20]/l[2] + inValues[12];
outValues[13] = (J_2[0][0]*inValues[13] + J_2[0][1]*inValues[14] - 0.4375*b_1[1][0]*inValues[19]*t[1][0] - 0.4375*b_2[1][0]*inValues[20]*t[2][0])/h[2];
outValues[14] = (J_2[1][0]*inValues[13] + J_2[1][1]*inValues[14] - 0.4375*b_1[1][0]*inValues[19]*t[1][1] - 0.4375*b_2[1][0]*inValues[20]*t[2][1])/h[2];
outValues[15] = ((1.0/32.0)*b_1[1][0]*inValues[19]*l[1]*tau[1][0] + (1.0/32.0)*b_2[1][0]*inValues[20]*l[2]*tau[2][0] + inValues[15]*theta_2[0][0] + inValues[16]*theta_2[0][1] + inValues[17]*theta_2[0][2])/h[2] /h[2];
outValues[16] = ((1.0/32.0)*b_1[1][0]*inValues[19]*l[1]*tau[1][1] + (1.0/32.0)*b_2[1][0]*inValues[20]*l[2]*tau[2][1] + inValues[15]*theta_2[1][0] + inValues[16]*theta_2[1][1] + inValues[17]*theta_2[1][2])/h[2] /h[2];
outValues[17] = ((1.0/32.0)*b_1[1][0]*inValues[19]*l[1]*tau[1][2] + (1.0/32.0)*b_2[1][0]*inValues[20]*l[2]*tau[2][2] + inValues[15]*theta_2[2][0] + inValues[16]*theta_2[2][1] + inValues[17]*theta_2[2][2])/h[2] /h[2];
outValues[18] = b_0[0][0]*inValues[18]*(o[0] ? -1 : 1);
outValues[19] = b_1[0][0]*inValues[19]*(o[1] ? -1 : 1);
outValues[20] = b_2[0][0]*inValues[20]*(o[2] ? -1 : 1);
}
private:
typename Impl::ArgyrisReferenceLocalBasis<D, R> basis_;
typename Traits::LocalCoefficientsType coefficients_;
typename Traits::LocalInterpolationType interpolation_;
// the transformation to correct the lack of affine equivalence boils down to
// matrix without blockstructure, because the normal derivative dofs interact nontrivially
// with the others, for details see Kirby's paper.
// The following geometric information is needed
// the jacobians per vertex
std::array<Dune::FieldMatrix<R, dim, dim>, dim+1> vertexJacobians_;
// a three by three matrix per vertex formalizing the rotation of the hessian in voigt notation
std::array<FieldMatrix<R, 3, 3>, 3> theta;
// the edge lengths
FieldVector<R, 3> l;
// the global tangent vectors (normalized)
std::array<FieldVector<R, 2>, 3> globalTangents;
// Geometric quantities without trivial meaning
std::array<FieldVector<R, 3>, 3> tau;
std::array<FieldMatrix<R, 2, 2>, 3> b;
// Additionally, we collect some global information
// the local state, i.e. a collection of global information restricted to this element
std::array<D, dim+1> averageVertexMeshSize_;
std::bitset<3> edgeOrientation_;
};
} // end namespace Impl
// *****************************************************************************
// This is the reusable part of the basis. It contains
//
// ArgyrisPreBasis
// ArgyrisNode
//
// The pre-basis allows to create the others and is the owner of possible shared
// state. These components do _not_ depend on the global basis and local view
// and can be used without a global basis.
// *****************************************************************************
template<class GV, class R>
class ArgyrisNode
: public LeafBasisNode
{
using Mapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
public:
using size_type = std::size_t;
using Element = typename GV::template Codim<0>::Entity;
using FiniteElement = typename Impl::ArgyrisLocalFiniteElement<typename GV::ctype, R>;
ArgyrisNode(Mapper const& vertexMapper, Mapper const& elementMapper,
std::vector<typename GV::ctype> const& averageVertexMeshSize,
std::vector<std::bitset<3> > const& edgeOrientation)
: element_(nullptr)
, vertexMapper_(&vertexMapper)
, elementMapper_(&elementMapper)
, averageVertexMeshSize_(&averageVertexMeshSize)
, edgeOrientation_(&edgeOrientation)
{}
Element const &element() const
{
return *element_;
}
FiniteElement const &finiteElement() const
{
return finiteElement_;
}
void bind(Element const &e)
{
element_ = &e;
finiteElement_.bind(*vertexMapper_, *averageVertexMeshSize_,
*elementMapper_,* edgeOrientation_,
*element_ );
this->setSize(finiteElement_.size());
}
unsigned int order() const { return finiteElement_.localBasis().order(); }
protected:
FiniteElement finiteElement_;
Element const* element_;
Mapper const* vertexMapper_;
Mapper const* elementMapper_;
std::vector<typename GV::ctype> const* averageVertexMeshSize_;
std::vector<std::bitset<3> > const* edgeOrientation_;
};
template<class GV, class R>
class ArgyrisPreBasis
: public LeafPreBasisMapperMixin<GV>
{
using Base = LeafPreBasisMapperMixin<GV>;
using SubEntityMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
using Element = typename GV::template Codim<0>::Entity;
using D = typename GV::ctype;
static const std::size_t dim = GV::dimension;
// helper methods to assign each subentity the number of dofs. Used by the LeafPreBasisMapperMixin.
static constexpr auto ArgyrisMapperLayout(Dune::GeometryType type, int gridDim)
{
assert(gridDim == 2);
if (type.isVertex())
return 6; // one evaluation dof and two derivative dofs and three hessian dofs per vertex
else if (type.isLine())
return 1;
else if (type.isTriangle())
return 0;
else
DUNE_THROW(Dune::Exception, "Invalid Geometry type for Argyris Element!");
}
public:
using GridView = GV;
using size_type = std::size_t;
using Node = ArgyrisNode<GridView, R>;
public:
ArgyrisPreBasis(const GV &gv)
: Base(gv, ArgyrisMapperLayout)
, vertexMapper_({gv, mcmgVertexLayout()})
, elementMapper_({gv, mcmgElementLayout()})
{
averageVertexMeshSize_ = Impl::computeAverageSubEntityMeshSize<D>(vertexMapper_);
edgeOrientations_ = Impl::computeEdgeOrientations(elementMapper_);
}
void update(GridView const &gv)
{
Base::update(gv);
vertexMapper_.update(this->gridView());
elementMapper_.update(this->gridView());
averageVertexMeshSize_ = Impl::computeAverageSubEntityMeshSize<D>(vertexMapper_);
edgeOrientations_ = Impl::computeEdgeOrientations(elementMapper_);
}
Node makeNode() const
{
return Node{vertexMapper_, elementMapper_, averageVertexMeshSize_, edgeOrientations_};
}
protected:
SubEntityMapper vertexMapper_;
std::vector<D> averageVertexMeshSize_;
SubEntityMapper elementMapper_;
std::vector<std::bitset<3> > edgeOrientations_;
}; // class ArgyrisPreBasis
namespace BasisFactory
{
template<class R = double>
auto argyris()
{
return [=](auto const &gridView) {
return ArgyrisPreBasis<std::decay_t<decltype(gridView)>, R>(gridView);
};
}
} // end namespace BasisFactory
} // end namespace Dune::Functions
#endif
nodes.hh
leafprebasismappermixin.hh
functionaldescriptor.hh
cubichermitebasis.hh
This file provides an implementation of the cubic Hermite finite element in 1 to 3 dimensions.
defaultglobalbasis.hh
transformedfiniteelementmixin.hh
monomialset.hh
mapperutilities.hh
Dune::Functions::derivative
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition trigonometricfunction.hh:43
Dune::Functions::ArgyrisBasis
DefaultGlobalBasis< ArgyrisPreBasis< GV, R > > ArgyrisBasis
Nodal basis of a scalar cubic Argyris finite element space.
Definition argyrisbasis.hh:69
Dune::Functions
Definition monomialset.hh:19
Dune::Functions::interpolate
void interpolate(const B &basis, C &&coeff, const F &f, const BV &bv, const NTRE &nodeToRangeEntry)
Interpolate given function in discrete function space.
Definition interpolate.hh:205
Dune::Functions::BasisFactory
Definition argyrisbasis.hh:926
Dune::Functions::MonomialSet
Function, which evaluates all monomials up to degree maxDegree in a given coordinate.
Definition monomialset.hh:64
Dune::Functions::ArgyrisPreBasis
A pre-basis for a Argyrisbasis.
Definition argyrisbasis.hh:856
Dune::Functions::ArgyrisPreBasis::vertexMapper_
SubEntityMapper vertexMapper_
Definition argyrisbasis.hh:918
Dune::Functions::ArgyrisPreBasis::elementMapper_
SubEntityMapper elementMapper_
Definition argyrisbasis.hh:920
Dune::Functions::ArgyrisPreBasis::edgeOrientations_
std::vector< std::bitset< 3 > > edgeOrientations_
Definition argyrisbasis.hh:921
Dune::Functions::ArgyrisPreBasis::makeNode
Node makeNode() const
Create tree node.
Definition argyrisbasis.hh:911
Dune::Functions::ArgyrisPreBasis::GridView
GV GridView
The grid view that the FE basis is defined on.
Definition argyrisbasis.hh:880
Dune::Functions::ArgyrisPreBasis::averageVertexMeshSize_
std::vector< D > averageVertexMeshSize_
Definition argyrisbasis.hh:919
Dune::Functions::ArgyrisPreBasis::Node
ArgyrisNode< GridView, R > Node
Template mapping root tree path to type of created tree node.
Definition argyrisbasis.hh:886
Dune::Functions::ArgyrisPreBasis::update
void update(GridView const &gv)
Update the stored grid view, to be called if the grid has changed.
Definition argyrisbasis.hh:901
Dune::Functions::ArgyrisPreBasis::ArgyrisPreBasis
ArgyrisPreBasis(const GV &gv)
Constructor for a given grid view object.
Definition argyrisbasis.hh:891
Dune::Functions::ArgyrisPreBasis::size_type
std::size_t size_type
Type used for indices and size information.
Definition argyrisbasis.hh:883
Dune::Functions::ArgyrisNode
Definition argyrisbasis.hh:786
Dune::Functions::ArgyrisNode< GridView, R >::averageVertexMeshSize_
std::vector< typename GridView::ctype > const * averageVertexMeshSize_
Definition argyrisbasis.hh:839
Dune::Functions::ArgyrisNode::bind
void bind(Element const &e)
Bind to element.
Definition argyrisbasis.hh:821
Dune::Functions::ArgyrisNode< GridView, R >::edgeOrientation_
std::vector< std::bitset< 3 > > const * edgeOrientation_
Definition argyrisbasis.hh:840
Dune::Functions::ArgyrisNode::FiniteElement
typename Impl::ArgyrisLocalFiniteElement< typename GV::ctype, R > FiniteElement
Definition argyrisbasis.hh:792
Dune::Functions::ArgyrisNode::element
Element const & element() const
Return current element, throw if unbound.
Definition argyrisbasis.hh:805
Dune::Functions::ArgyrisNode::size_type
std::size_t size_type
Definition argyrisbasis.hh:790
Dune::Functions::ArgyrisNode< GridView, R >::finiteElement_
FiniteElement finiteElement_
Definition argyrisbasis.hh:835
Dune::Functions::ArgyrisNode::ArgyrisNode
ArgyrisNode(Mapper const &vertexMapper, Mapper const &elementMapper, std::vector< typename GV::ctype > const &averageVertexMeshSize, std::vector< std::bitset< 3 > > const &edgeOrientation)
Definition argyrisbasis.hh:794
Dune::Functions::ArgyrisNode< GridView, R >::elementMapper_
Mapper const * elementMapper_
Definition argyrisbasis.hh:838
Dune::Functions::ArgyrisNode::finiteElement
FiniteElement const & finiteElement() const
Return the LocalFiniteElement for the element we are bound to.
Definition argyrisbasis.hh:815
Dune::Functions::ArgyrisNode< GridView, R >::vertexMapper_
Mapper const * vertexMapper_
Definition argyrisbasis.hh:837
Dune::Functions::ArgyrisNode::order
unsigned int order() const
The order of the local basis.
Definition argyrisbasis.hh:832
Dune::Functions::ArgyrisNode< GridView, R >::element_
Element const * element_
Definition argyrisbasis.hh:836
Dune::Functions::ArgyrisNode::Element
typename GV::template Codim< 0 >::Entity Element
Definition argyrisbasis.hh:791
Dune::Functions::H2LocalBasisTraits
Definition cubichermitebasis.hh:84
Dune::Functions::DefaultGlobalBasis
Global basis for given pre-basis.
Definition defaultglobalbasis.hh:53
Dune::Functions::LeafPreBasisMapperMixin::gridView
const GridView & gridView() const
Export the stored GridView.
Definition leafprebasismappermixin.hh:95
Dune::Functions::LeafPreBasisMapperMixin::update
void update(const GridView &gv)
Update the stored GridView.
Definition leafprebasismappermixin.hh:101
Dune::Functions::LeafPreBasisMapperMixin::LeafPreBasisMapperMixin
LeafPreBasisMapperMixin(const GridView &gv, Dune::MCMGLayout layout)
Construct from GridView and local DOF layout.
Definition leafprebasismappermixin.hh:74
Dune::Functions::BasisNodeMixin::setSize
void setSize(const size_type size)
Definition nodes.hh:197
Dune::Functions::LeafBasisNode
Definition nodes.hh:218