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

globalInterpolation and rotation derivative dofs for hermite

parent 596e4f40
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
+ 269
− 30
View file @ b97a2780
......@@ -48,37 +48,33 @@ namespace Dune
class ElementInformation
{
using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
using ctype = typename Element::Geometry::ctype;
static_assert(std::is_same<GlobalCoordinate, FieldVector<ctype, dim>>::
value); // not sure my code works for other types of global coordinates
public:
ElementInformation()
{
for (std::size_t i = 0; i < dim + 1; ++i)
{
// default directions are global coordinates
if constexpr (dim == 1)
derivativeDirections[i][0] = GlobalCoordinate{1.};
else if constexpr (dim == 2)
{
derivativeDirections[i][0] = GlobalCoordinate{1., 0.};
derivativeDirections[i][1] = GlobalCoordinate{0., 1.};
}
else if constexpr (dim == 3)
{
derivativeDirections[i][0] = GlobalCoordinate{1., 0., 0.};
derivativeDirections[i][1] = GlobalCoordinate{0., 1., 0.};
derivativeDirections[i][2] = GlobalCoordinate{0., 0., 1.};
}
derivativeDirections[i] = 0;
for (std::size_t j = 0; j < dim; ++j)
derivativeDirections[i][j][j] = 1.;
}
}
ElementInformation(std::array<FieldMatrix<ctype, dim, dim>, dim + 1> directions)
: derivativeDirections(directions)
{
}
std::array<std::array<GlobalCoordinate, dim>, dim + 1> getDerivativeDirections() const
std::array<FieldMatrix<ctype, dim, dim>, dim + 1> const &getDerivativeDirections() const
{
return derivativeDirections;
}
private:
std::array<std::array<GlobalCoordinate, dim>, dim + 1> derivativeDirections;
std::array<FieldMatrix<ctype, dim, dim>, dim + 1> derivativeDirections;
};
private:
......@@ -95,8 +91,7 @@ namespace Dune
template <class Element>
void bind(Element const &e, ElementInformation<Element> const &elementInformation)
{
fillMatrix(Dune::ReferenceElements<double, dim>::simplex().position(0, 0),
e.geometry()); // barycenter, because we need some value.
fillMatrix(e.geometry(), elementInformation);
}
template <typename Values, typename LocalCoordinate, typename Geometry>
......@@ -147,18 +142,25 @@ namespace Dune
* |0 1| } repeat (dim-1)^2 times
*
*/
template <class Geometry, typename LocalCoordinate>
void fillMatrix(LocalCoordinate const &xi, Geometry const &geometry)
template <class Geometry, typename Element>
void fillMatrix(Geometry const &geometry, ElementInformation<Element> const &elementInfo)
{
auto jacobianTransposed = geometry.jacobianTransposed(xi);
auto const &refElement
= Dune::ReferenceElements<typename Geometry::ctype, dim>::simplex();
auto const &directions = elementInfo.getDerivativeDirections();
// auto jacobianTransposed = geometry.jacobianTransposed(xi);
for (std::size_t i = 0; i < numberOfVertices; ++i) // dim + 1 vertices
{
FieldMatrix<R, dim, dim> dir = directions[i];
dir.invert();
auto directionalJacobianTransposed
= geometry.jacobianTransposed(refElement.position(i, dim)) * dir;
mat_[numberOfVertices * i][numberOfVertices * i] = 1.;
for (std::size_t j = 0; j < dim; ++j)
for (std::size_t k = 0; k < dim; ++k)
{
mat_[numberOfVertices * i + 1 + j][numberOfVertices * i + 1 + k]
= jacobianTransposed[k][j];
= directionalJacobianTransposed[k][j]; // jacobianTransposed[k][j];
}
}
for (std::size_t i = 0; i < (dim - 1) * (dim - 1); ++i) // inner dofs
......@@ -170,8 +172,95 @@ namespace Dune
BCRSMatrix<R> mat_;
public:
static constexpr bool globalInterpolation = false;
static constexpr bool globalInterpolation = true;
template <class LocalBasis, class Element>
class GlobalInterpolation
{
using size_type = std::size_t;
using LocalCoordinate = typename LocalBasis::Traits::DomainType;
using ctype = typename Element::Geometry::ctype;
static constexpr size_type numberOfVertices = dim + 1;
static constexpr size_type innerDofCodim
= (dim == 2) ? 0 : 1; // probably wrong for dim > 3
static constexpr size_type numberOfInnerDofs
= (dim - 1) * (dim - 1); // probably wrong for dim > 3
public:
GlobalInterpolation() {}
template <class LocalInterpolation>
void bind(const LocalInterpolation &dummyForCompability, const Element &element,
const ElementInformation<Element> &elementInfo)
{
element_ = &element;
elementInfo_ = elementInfo;
}
/** \brief Evaluate a given function at the Lagrange nodes
*
* \tparam F Type of function to evaluate
* \tparam C Type used for the values of the function
* \param[in] ff Function to evaluate
* \param[out] out Array of function values
*/
template <typename F, typename C>
void interpolate(const F &ff, std::vector<C> &out) const
{
auto &&f
= Dune::Impl::makeFunctionWithCallOperator<typename LocalBasis::Traits::DomainType>(
ff);
auto df = derivative(ff);
out.resize(LocalBasis::size());
auto const &refElement = Dune::ReferenceElements<ctype, dim>::simplex();
// Iterate over vertices, dim dofs per vertex
for (size_type i = 0; i < dim + 1; ++i)
{
LocalCoordinate x = refElement.position(i, dim);
// matrix storing directions as columns (!)
auto const &directionMatrix = elementInfo_.getDerivativeDirections()[i];
auto const &derivativeValue = vectorToMatrix(df(x)) * transpose(directionMatrix);
out[i * numberOfVertices] = f(x);
for (size_type d = 0; d < dim; ++d)
out[i * numberOfVertices + d + 1] = matrixToVector(derivativeValue)[d];
}
for (size_type i = 0; i < numberOfInnerDofs; ++i)
{
out[numberOfVertices * numberOfVertices + i]
= f(refElement.position(i, innerDofCodim));
}
}
private:
template <class V, int N>
FieldVector<V, N> const &matrixToVector(FieldVector<V, N> const &v) const
{
return v;
}
template <class V, int N>
FieldVector<V, N> const &matrixToVector(FieldMatrix<V, 1, N> const &m) const
{
return m[0];
}
template <class V, int N>
FieldMatrix<V, 1, N> &&vectorToMatrix(FieldVector<V, N> const &v) const
{
return FieldMatrix<V, 1, N>({v});
}
template <class V, int N>
FieldMatrix<V, 1, N> const &vectorToMatrix(FieldMatrix<V, 1, N> const &m) const
{
return m;
}
protected:
const Element *element_; // TODO remove
ElementInformation<Element> elementInfo_;
};
template <class Function, class LocalCoordinate, class Element>
class LocalValuedFunction
{
......@@ -340,6 +429,8 @@ namespace Dune
static_assert(GV::dimension == dim);
using ElementInformation =
typename HermiteTransformator<dim, false, R>::template ElementInformation<Element>;
using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
using IndexType = typename GV::IndexSet::IndexType;
public:
HermiteElementInformationMap(GV const &gv)
......@@ -362,20 +453,168 @@ namespace Dune
}
private:
void fill(const GV &gv)
void fill(const GV &gv, std::enable_if_t<dim == 2, int> = 0)
{
// TODO calculate tangetials for boundary elements
for (auto const &element : elements(gv))
elementInformation_[elementMapper_.index(element)] = ElementInformation();
// compute orientation for all elements
const auto &indexSet = gv.indexSet();
// vector with directions and bools to indicate whether this dir was already set
FieldMatrix<D, dim, dim> eye = 0;
for (std::size_t i = 0; i < dim; ++i)
eye[i][i] = 1.;
std::tuple<FieldMatrix<D, dim, dim>, std::array<bool, dim>> defaultInfo;
std::get<0>(defaultInfo) = eye;
std::get<1>(defaultInfo).fill(false);
std::vector<std::tuple<FieldMatrix<D, dim, dim>, std::array<bool, dim>>>
directionPerVertex(indexSet.size(dim), defaultInfo);
// interate over intersections
for (const auto &element : elements(gv))
{
for (const auto &intersection : intersections(gv, element))
{
// fill vertex Data with linear independent tangentials
if (intersection.boundary())
{
auto facet
= intersection.inside().template subEntity<1>(intersection.indexInInside());
auto tangential
= intersection.geometry().corner(1) - intersection.geometry().corner(0);
auto startIndex = indexSet.subIndex(facet, 0, dim);
auto endIndex = indexSet.subIndex(facet, 1, dim);
setVertexData(directionPerVertex[startIndex], tangential);
setVertexData(directionPerVertex[endIndex], tangential);
}
// inner dofs, keep global axes
}
}
// interate over vector and set the remaining direction to normals
for (auto &[dir, isSet] : directionPerVertex)
if (isSet[0])
{
if (isSet[1])
continue; // both are tangetials, i.e. vertex is corner
// else, set second direction normal to first
dir[1] = Dune::FieldVector<D, dim>{-dir[0][1], dir[0][0]};
}
// do nothing if non was set
// booleans now also encode, whether the dof is tangential and thus to be used for
// dirichlet interpolation
// TODO hand this to the element information
for (const auto &element : elements(gv))
{
auto elementIndex = elementMapper_.index(element);
std::array<Dune::FieldMatrix<D, dim, dim>, dim + 1> directions;
for (std::size_t i = 0; i < element.subEntities(dim); ++i)
{
auto vertexIndex = indexSet.subIndex(element, i, dim);
directions[i] = std::get<0>(directionPerVertex[vertexIndex]);
}
elementInformation_[elementIndex] = ElementInformation(directions);
}
}
Dune::MultipleCodimMultipleGeomTypeMapper<GV> elementMapper_;
std::vector<ElementInformation> elementInformation_;
bool linearIndependent(Dune::FieldVector<D, 2> a, Dune::FieldVector<D, 2> b)
{
return std::abs(a.dot(Dune::FieldVector<D, dim>{-b[1], b[0]})) > 1e-14;
}
void setVertexData(
std::tuple<FieldMatrix<D, dim, dim>, std::array<bool, dim>> &directionPerVertex,
GlobalCoordinate const &tangential)
{
auto &[dir, isSet] = directionPerVertex;
if (!isSet[0])
{
dir[0] = tangential;
isSet[0] = true;
}
else if (linearIndependent(dir[0], tangential))
{
if (isSet[1])
if (linearIndependent(dir[1], tangential))
DUNE_THROW(Dune::NotImplemented,
"More than two linear independent tangentials at vertex");
else
; // linear dependent to second, do nothing
else // second not set
{
dir[1] = tangential;
isSet[1] = true;
}
}
// else; // linear dependent to first, do nothing!
}
};
} // namespace Impl
template <typename GV, typename R, bool useSpecialization>
class HermiteNode;
// spezialization for dim == 1, no need to do anything here
template <typename GV, typename R>
class HermiteElementInformationMap<GV, 1, R>
{
using D = typename GV::ctype;
using Element = typename GV::template Codim<0>::Entity;
static_assert(GV::dimension == 1);
using ElementInformation =
typename HermiteTransformator<1, false, R>::template ElementInformation<Element>;
using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
using IndexType = typename GV::IndexSet::IndexType;
public:
HermiteElementInformationMap(GV const &gv)
{
defaultElementInfo_[0] = Dune::FieldMatrix<D, 1, 1>({{1.}});
defaultElementInfo_[1] = Dune::FieldMatrix<D, 1, 1>({{1.}});
}
void update(GV const &gv) {}
const auto &find(const Element &element) const { return defaultElementInfo_; }
private:
std::array<Dune::FieldMatrix<D, 1, 1>, 2> defaultElementInfo_;
};
// spezialization for dim == 3, no default solution here, since one vertex on boundary has
// arbitrarily many boundary tangentials
template <typename GV, typename R>
class HermiteElementInformationMap<GV, 3, R>
{
using D = typename GV::ctype;
using Element = typename GV::template Codim<0>::Entity;
static_assert(GV::dimension == 3);
using ElementInformation =
typename HermiteTransformator<3, false, R>::template ElementInformation<Element>;
using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
using IndexType = typename GV::IndexSet::IndexType;
public:
HermiteElementInformationMap(GV const &gv)
{
for (std::size_t vertex = 0; vertex < 4; ++vertex)
{
defaultElementInfo_[vertex] = 0;
for (std::size_t i = 0; i < 3; ++i)
defaultElementInfo_[vertex][i][i] = 1.;
}
}
void update(GV const &gv) {}
const auto &find(const Element &element) const { return defaultElementInfo_; }
private:
std::array<FieldMatrix<D, 3, 3>, 4> defaultElementInfo_;
};
} // namespace Impl
// primary template
template <typename GV, typename R, bool useSpezialization>
class HermiteNode
{
};
// specialization for globalvaluedFiniteElement
template <typename GV, typename R>
class HermiteNode<GV, R, true>: public LeafBasisNode
......
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