Skip to content
Snippets Groups Projects
Commit 5bb16dc4 authored by Porrmann, Maik's avatar Porrmann, Maik
Browse files

LocalTest Hermite, Bugfixes Hermite<3>

parent f010a65a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/functions/functionspacebases/hermitebasis.hh
+ 20
− 13
View file @ 5bb16dc4
......@@ -132,14 +132,14 @@ namespace Dune
static auto fillMatrix(Geometry const &geometry, LocalCoordinate const &xi, Matrix &m)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
std::size_t NumberOfVertices = dim + 1;
for (std::size_t i = 0; i < NumberOfVertices; ++i) // dim + 1 vertices
for (std::size_t i = 0; i < numberOfVertices; ++i) // dim + 1 vertices
{
m[NumberOfVertices * i][NumberOfVertices * i] = 1.;
m[numberOfVertices * i][numberOfVertices * i] = 1.;
for (std::size_t j = 0; j < dim; ++j)
for (std::size_t k = 0; k < dim; ++k)
{
m[NumberOfVertices * i + 1 + j][NumberOfVertices * i + 1 + k] = jacobianTransposed[k][j];
m[numberOfVertices * i + 1 + j][numberOfVertices * i + 1 + k]
= jacobianTransposed[k][j];
}
}
for (std::size_t i = 0; i < (dim - 1) * (dim - 1); ++i) // inner dofs
......@@ -331,8 +331,6 @@ namespace Dune
const Element *element_;
};
// template<typename GV, class MI, typename R=double>
// class HermitePreBasis;
/**
* \brief A pre-basis for a Hermitebasis
*
......@@ -345,8 +343,6 @@ namespace Dune
template <typename GV, typename R, bool useSpecialization>
class HermitePreBasis
{
static const int dim = GV::dimension;
public:
//! The grid view that the FE basis is defined on
using GridView = GV;
......@@ -361,6 +357,11 @@ namespace Dune
static constexpr size_type minMultiIndexSize = 1;
static constexpr size_type multiIndexBufferSize = 1;
private:
static const size_type dim = GV::dimension;
static constexpr size_type innerDofCodim = (dim == 2) ? 0 : 1;
public:
//! Constructor for a given grid view object
HermitePreBasis(const GV &gv)
: gridView_(gv)
......@@ -370,9 +371,7 @@ namespace Dune
}
//! Initialize the global indices
void initializeIndices()
{ //TODO check if we need to do something
}
void initializeIndices() {}
//! Obtain the grid view that the basis is defined on
const GridView &gridView() const
......@@ -398,7 +397,7 @@ namespace Dune
size_type size() const
{
if (dim != 0 && dim <= 3)
return (dim + 1) * gridView_.size(dim) + (dim - 1) * (dim - 1) * gridView_.size(0);
return (dim + 1) * gridView_.size(dim) + ((dim == 1) ? 0 : gridView_.size(innerDofCodim));
else
DUNE_THROW(Dune::NotImplemented, "No size method for " << dim << "d grids available yet!");
}
......@@ -443,7 +442,15 @@ namespace Dune
if constexpr (dim <= 3)
{
*it = {{(size_type)((localKey.codim() == dim) ? ((dim + 1) * gridIndexSet.subIndex(element, localKey.subEntity(), dim) + localKey.index()) : (dim + 1) * gridIndexSet.size(dim) + (dim - 1) * (dim - 1) * gridIndexSet.subIndex(element, localKey.subEntity(), 0) + localKey.index())}};
*it = {
{(size_type)((localKey.codim() == dim)
? ((dim + 1)
* gridIndexSet.subIndex(element, localKey.subEntity(), dim)
+ localKey.index())
: (dim + 1) * gridIndexSet.size(dim)
+ gridIndexSet.subIndex(element, localKey.subEntity(),
innerDofCodim)
+ localKey.index())}};
}
else
DUNE_THROW(Dune::NotImplemented, "indices() for Hermite Elements not implemented for dim > 3");
......
dune/functions/functionspacebases/test/hermitetest.cc
+ 4
− 4
View file @ 5bb16dc4
......@@ -81,13 +81,13 @@ int main(int argc, char *argv[])
{ // 3d
std::cout << "Hermite test in 3d" << std::endl;
auto grid
= StructuredGridFactory<UGGrid<3>>::createSimplexGrid({0., 0., 0.}, {1., 1., 1.}, {{3, 3, 3}});
auto grid = StructuredGridFactory<UGGrid<3>>::createSimplexGrid({0., 0., 0.}, {1., 1., 1.},
{{3, 3, 3}});
auto gridView = grid->leafGridView();
std::cout << "Grid has " << gridView.size(0) << " elementes and " << gridView.size(1)
<< " facettes and " << gridView.size(2) << " vertices" << std::endl;
// using GridView = decltype(gridView);
<< " facettes and " << gridView.size(2) << " edges and " << gridView.size(3)
<< " vertices " << std::endl;
{
using namespace Dune::Functions::BasisFactory;
auto basis = makeBasis(gridView, hermite());
......
dune/localfunctions/hermite.hh
+ 66
− 45
View file @ 5bb16dc4
......@@ -39,7 +39,7 @@ namespace Dune
{
public:
using Eval = StandardEvaluator<VirtualMonomialBasis<dim, D>>;
using Base = PolynomialBasisWithMatrix<Eval>;
using Base = PolynomialBasisWithMatrix<Eval, SparseCoeffMatrix<D, 1>, D, R>;
using Traits = LocalBasisTraits<D, dim, FieldVector<D, dim>, R, 1, FieldVector<R, 1>, FieldMatrix<R, 1, dim>>;
private:
......@@ -61,7 +61,45 @@ namespace Dune
DUNE_THROW(Dune::NotImplemented, "only implemented for dim <= 3");
}
HermiteLocalBasis(HermiteLocalBasis const &) : HermiteLocalBasis() {}
static unsigned int size() { return coeffSize; }
void evaluateFunction(const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
// First we evaluate the basis with default edge orientation, then we adapt the result
// accordingly This works because the edge orientation only changes the "wrong" derivative
// nodes which are multiplied by -1, and the very same happens to the nodal basis
out.resize(size());
Base::evaluateFunction(in, out);
}
/** \brief Evaluate Jacobians of all shape functions at a given point
*
* \param[in] in The evaluation point
* \param[out] out Jacobians of all shape functions at that point
*/
void evaluateJacobian(const typename Traits::DomainType &in,
std::vector<typename Traits::JacobianType> &out) const
{
out.resize(size());
Base::evaluateJacobian(in, out);
}
/** \brief Evaluate partial derivatives of all shape functions at a given point
*
* \param[in] order The partial derivative to be computed, as a multi-index
* \param[in] in The evaluation point
* \param[out] out Jacobians of all shape functions at that point
*/
void partial(const std::array<unsigned int, dim> &order,
const typename Traits::DomainType &in,
std::vector<typename Traits::RangeType> &out) const
{
out.resize(size());
Base::partial(order, in, out);
}
};
/** \brief Associations of the Hermite degrees of freedom to subentities of the
......@@ -73,6 +111,7 @@ namespace Dune
class HermiteLocalCoefficients
{
using size_type = std::size_t;
static constexpr size_type innerDofCodim = (dim == 2) ? 0 : 1;
public:
HermiteLocalCoefficients() : localKeys_(size())
......@@ -80,14 +119,15 @@ namespace Dune
if constexpr (dim > 0 && dim <= 3)
{
std::size_t numberOfVertices = dim + 1;
std::size_t numberOfInnerDofs = (dim - 1) * (dim - 1);
std::size_t numberOfInnerDofs = (dim - 1) * (dim - 1); // probably incorrect for dim > 3
for (size_type i = 0; i < numberOfVertices; ++i) // subentities: vertices
{
for (size_type k = 0; k < numberOfVertices; ++k) // dim + 1 dofs per subentity
localKeys_[numberOfVertices * i + k] = LocalKey(i, dim, k);
}
for (size_type i = 0; i < numberOfInnerDofs; ++i) // subentities: element
localKeys_[numberOfVertices * numberOfVertices + i] = LocalKey(0, 0, i); // inner dofs
localKeys_[numberOfVertices * numberOfVertices + i]
= LocalKey(i, innerDofCodim, 0); // inner dofs
}
else
DUNE_THROW(NotImplemented,
......@@ -120,6 +160,12 @@ namespace Dune
{
using size_type = std::size_t;
using LocalCoordinate = typename LocalBasis::Traits::DomainType;
static constexpr size_type dim = LocalBasis::Traits::dimDomain;
static constexpr size_type numberOfVertices = dim + 1;
static constexpr size_type innerDofCodim = (dim == 2) ? 0 : 1;
static constexpr size_type numberOfInnerDofs
= (dim - 1) * (dim - 1); // probably wrong for dim > 3
public:
/** \brief Evaluate a given function at the Lagrange nodes
......@@ -132,57 +178,30 @@ namespace Dune
template <typename F, typename C>
void interpolate(const F &ff, std::vector<C> &out) const
{
constexpr auto dim = LocalBasis::Traits::dimDomain;
// constexpr auto k = LocalBasis::order();
using D = typename LocalBasis::Traits::DomainFieldType;
auto &&f = Impl::makeFunctionWithCallOperator<LocalCoordinate>(ff);
auto &&df = derivative(ff); // makeDerivative at this position might return
// reference gradient for localfunctions!
auto &&df = makeDerivative<LocalCoordinate>(
ff); // makeDerivative at this position might return
// reference gradient for localfunctions!
// This is necessary for the test routine to compile, TODO change that
out.resize(LocalBasis::size());
if constexpr (dim == 1)
auto refElement = Dune::ReferenceElements<double, dim>::simplex();
if (dim > 0 && dim <= 3)
{
typename LocalBasis::Traits::DomainType x;
for (size_type i = 0; i < 2; i++)
for (size_type i = 0; i < numberOfVertices; ++i)
{
x[0] = (D)i;
out[2 * i] = f(x);
out[2 * i + 1] = df(x);
}
return;
}
else if constexpr (dim == 2)
{
LocalCoordinate const vertices[3]{{0, 0}, {1, 0}, {0, 1}};
LocalCoordinate x = refElement.position(i, dim);
// std::cout << "Vertex " << i << ": " << x << std::endl;
out[numberOfVertices * i] = f(x);
for (size_type i = 0; i < 3; ++i) // three vertices
{
out[3 * i] = f(vertices[i]);
auto const &grad = df(vertices[i]);
out[3 * i + 1] = getPartialDerivative(grad, 0);
out[3 * i + 2] = getPartialDerivative(grad, 1);
auto const &grad = df(x);
for (size_type d = 0; d < dim; ++d)
out[numberOfVertices * i + d + 1] = getPartialDerivative(grad, d);
}
// inner dof
out[9] = f({1. / 3., 1. / 3.});
}
else if constexpr (dim == 3)
{
LocalCoordinate const vertices[4]{{0, 0, 1}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
for (size_type i = 0; i < 4; ++i) // four vertices
for (size_type i = 0; i < numberOfInnerDofs; ++i)
{
out[3 * i] = f(vertices[i]);
auto const &grad = df(vertices[i]);
out[3 * i + 1] = getPartialDerivative(grad, 0);
out[3 * i + 2] = getPartialDerivative(grad, 1);
out[3 * i + 3] = getPartialDerivative(grad, 2);
out[numberOfVertices * numberOfVertices + i] = f(refElement.position(i, innerDofCodim));
}
D oneThird = 1. / 3.;
out[16] = f(LocalCoordinate{oneThird, oneThird, oneThird});
out[17] = f(LocalCoordinate{0., oneThird, oneThird});
out[18] = f(LocalCoordinate{oneThird, 0., oneThird});
out[19] = f(LocalCoordinate{oneThird, oneThird, 0.});
}
else
DUNE_THROW(Dune::NotImplemented,
......@@ -190,6 +209,8 @@ namespace Dune
}
private:
// TODO make this viable for general matrices
// TODO maybe make this part of a baseclass that also can be used by Argyris
template <class DerivativeType, class FieldType>
FieldType getPartialDerivative(DerivativeType const &df,
std::size_t i) const
......
dune/localfunctions/morley.hh
+ 1
− 12
View file @ 5bb16dc4
......@@ -257,24 +257,13 @@ namespace Dune
return f.normalDerivative(i);
}
// TODO there got to be a nicer way to capture Vectors than with the returntype sfinae
template <class F, decltype((derivative(std::declval<F>()), true)) = true>
auto normalDerivativeImpl(F const &f, size_type i, PriorityTag<3>) const
// -> decltype(derivative(std::declval<F>())(std::declval<LocalCoordinate>())
// .dot(std::declval<LocalCoordinate>()))
{
return matrixToVector(derivative(f)(m_[i])).dot(normal_[i]);
}
// TODO This is supposed to capture Matrix<1,N>, but it also gets called for Matrix<N,M> with
// N != 0 template <class F, decltype((derivative(std::declval<F>()), true)) = true> auto
// normalDerivativeImpl(F const &f, size_type i, PriorityTag<2>) const
// -> decltype(derivative(std::declval<F>())(std::declval<LocalCoordinate>())[0].dot(
// std::declval<LocalCoordinate>()))
// {
// return derivative(f)(m_[i])[0].dot(normal_[i]);
// }
// TODO this turns compile Errors into runtime errors. Not sure if this is clever
template <class F>
auto normalDerivativeImpl(F const &f, size_type i, PriorityTag<1>) const
{
......
dune/localfunctions/polynomialbasiscoefficients.hh
+ 12
− 11
View file @ 5bb16dc4
......@@ -86,12 +86,12 @@ namespace Dune
{1, 0, 0, 0, -3, -13, -3, -13, -13, -3, // deg 0 to 2
2, 13, 13, 2, 13, 33, 13, 13, 13, 2}, // deg 3
{0, 1, 0, 0, -2, -3, 0, -3, 0, 0,
1, 3, 2, 0, 3, 4, 0, 3, 0, 0},
1, 3, 2, 0, 3, 4, 0, 2, 0, 0},
{0, 0, 1, 0, 0, -3, -2, 0, -3, 0,
0, 2, 3, 1, 0, 4, 3, 0, 2, 0},
{0, 0, 0, 1, 0, 0, 0, -3, -3, -2,
0, 0, 0, 0, 2, 4, 2, 3, 3, 1},
{0, 0, 0, 0, 3, 7, 0, -7, 0, 0, // l_4
{0, 0, 0, 0, 3, -7, 0, -7, 0, 0, // l_4
-2, 7, 7, 0, 7, 7, 0, 7, 0, 0},
{0, 0, 0, 0, -1, 2, 0, 2, 0, 0,
1, -2, -2, 0, -2, -2, 0, -2, 0, 0},
......@@ -106,7 +106,7 @@ namespace Dune
{0, 0, 0, 0, 0, 2, -1, 0, 2, 0,
0, -2, -2, 1, 0, -2, -2, 0, -2, 0},
{0, 0, 0, 0, 0, 0, 0, 0, -1, 0,
0, 0, 0, 0, 0, 0, 3, 0, 2, 0},
0, 0, 0, 0, 0, 0, 2, 0, 1, 0},
{0, 0, 0, 0, 0, 0, 0, -7, -7, 3, // l_12
0, 0, 0, 0, 7, 7, 7, 7, 7, -2},
{0, 0, 0, 0, 0, 0, 0, -1, 0, 0,
......@@ -114,15 +114,16 @@ namespace Dune
{0, 0, 0, 0, 0, 0, 0, 0, -1, 0,
0, 0, 0, 0, 0, 0, 1, 0, 2, 0},
{0, 0, 0, 0, 0, 0, 0, 2, 2, -1,
0, 0, 0, 0, -2, -2, -2, -2, -2, -1},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, // l_16, from here on inner dofs
0, 0, 0, 0, 0, 27, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 27, 0,
0, 0, 0, 0, 0, -27, -27, 0, -27, 0},
{0, 0, 0, 0, 0, 0, 0, 27, 0, 0,
0, 0, 0, 0, -2, -2, -2, -2, -2, 1},
// l_16, from here on inner dofs
{0, 0, 0, 0, 0, 27, 0, 0, 0, 0,//bottom
0, -27, -27, 0, 0, -27, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 0, 27, 0, 0,//front
0, 0, 0, 0, -27, -27, 0, -27, 0, 0},
{0, 0, 0, 0, 0, 27, 0, 0, 0, 0,
0, -27, -27, 0, 0, -27, 0, 0, 0, 0}});
{0, 0, 0, 0, 0, 0, 0, 0, 27, 0,//left
0, 0, 0, 0, 0, -27, -27, 0, -27, 0},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //right
0, 0, 0, 0, 0, 27, 0, 0, 0, 0}});
}
else
{
......
dune/localfunctions/test/localhermitetest.cc 0 → 100644
+ 24
− 0
View file @ 5bb16dc4
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include "config.h"
#include <dune/localfunctions/hermite.hh>
#include <dune/localfunctions/test/test-localfe.hh>
using namespace Dune;
int main(int argc, char **argv)
{
bool success = true; // captured by Macros
HermiteLocalFiniteElement<double, double, 1> LFE_1d;
// FE, DisableSubTests, max order for differentiabilitytest (<=2)
TEST_FE3(LFE_1d, DisableVirtualInterface, 1); // TODO fix error for order 2
HermiteLocalFiniteElement<double, double, 2> LFE_2d;
TEST_FE3(LFE_2d, DisableVirtualInterface, 1); // TODO fix error for order 2
HermiteLocalFiniteElement<double, double, 3> LFE_3d;
TEST_FE3(LFE_3d, DisableVirtualInterface, 1); // TODO fix error for order 2
return success ? 0 : 1;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment