Commit 47a93b82 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

implement a local-element test

parent e3b79c32
......@@ -385,7 +385,7 @@ public:
/// \brief Obtain the volume of the mapping's image
/**
* Calculates the volume of the entity by numerical integration. Since the
* polynomial order of the Volume element is not known, iterativly compute
* polynomial order of the Volume element is not known, iteratively compute
* numerical integrals with increasing order of the quadrature rules, until
* tolerance is reached.
**/
......@@ -428,6 +428,12 @@ public:
for (std::size_t i = 0; i < shapeJacobians.size(); ++i)
Impl::outerProductAccumulate(shapeJacobians[i], vertices_[i], out);
JacobianTransposed out1;
out1[0] = vertices_[1] - vertices_[0];
out1[1] = vertices_[2] - vertices_[0];
assert((out - out1).frobenius_norm() < 1.e-10);
return out;
}
......
......@@ -238,18 +238,21 @@ public:
**/
std::optional<LocalCoordinate> checkedLocal (const GlobalCoordinate& globalCoord) const
{
LocalCoordinate x = flatGeometry().local(globalCoord);
LocalCoordinate x = findClosestLocal(globalCoord);
LocalCoordinate dx;
for (int i = 0; i < Traits::maxIteration(); ++i)
{
// Newton's method: DF^n dx^n = F^n, x^{n+1} -= dx^n
const GlobalCoordinate dglobal = global(x) - globalCoord;
const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);
const bool invertible
= MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);
// break if jacobian is not invertible
if (!invertible)
if (!invertible) {
std::cout << "not invertible!" << std::endl;
return std::nullopt;
}
// update x with correction
x -= dx;
......@@ -259,9 +262,30 @@ public:
return x;
}
if (dx.two_norm2() > Traits::tolerance())
if (dx.two_norm2() > Traits::tolerance()) {
std::cout << "tolerance not reached!" << std::endl;
return std::nullopt;
}
return x;
}
LocalCoordinate findClosestLocal (const GlobalCoordinate& globalCoord, int order = 4) const
{
LocalCoordinate x = refElement().position(0,0);
ctype dGlobal = (global(x) - globalCoord).two_norm();
auto const& quad = QuadratureRules<ctype, mydimension>::rule(type(), order);
for (auto qp : quad) {
ctype d = (global(qp.position()) - globalCoord).two_norm();
if (d < dGlobal) {
dGlobal = d;
x = qp.position();
if (d < Traits::tolerance())
break;
}
}
return x;
}
......@@ -299,7 +323,7 @@ public:
/// \brief Obtain the volume of the mapping's image
/**
* Calculates the volume of the entity by numerical integration. Since the polynomial order of the
* Volume element is not known, iterativly compute numerical integrals with increasing order
* Volume element is not known, iteratively compute numerical integrals with increasing order
* of the quadrature rules, until tolerance is reached.
**/
Volume volume () const
......
......@@ -370,7 +370,9 @@ public:
**/
ctype integrationElement (const LocalCoordinate& local) const
{
return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
return isFlatAffine && order() == 1
? flatGeometry().integrationElement(local)
: MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
}
/// \brief Obtain the volume of the mapping's image
......@@ -382,6 +384,9 @@ public:
**/
Volume volume () const
{
if (isFlatAffine && order() == 1)
return flatGeometry().volume();
using std::abs;
Volume vol0 = volume(QuadratureRules<ctype, mydimension>::rule(type(), 1));
for (int p = 2; p < 10; ++p) {
......@@ -411,6 +416,9 @@ public:
**/
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
if (isFlatAffine && order() == 1)
return flatGeometry().jacobianTransposed(local);
thread_local std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
localBasis().evaluateJacobian(local, shapeJacobians);
assert(shapeJacobians.size() == vertices_.size());
......@@ -430,6 +438,9 @@ public:
**/
JacobianInverseTransposed jacobianInverseTransposed (const LocalCoordinate& local) const
{
if (isFlatAffine && order() == 1)
return flatGeometry().jacobianInverseTransposed(local);
JacobianInverseTransposed out;
MatrixHelper::rightInvA(jacobianTransposed(local), out);
return out;
......
dune_add_test(NAME test-element
SOURCES test-element.cc)
dune_add_test(NAME test-curvedgeometry-double
SOURCES test-curvedgeometry.cc
COMPILE_DEFINITIONS USE_FLOAT128=0
......
#include <dune/common/float_cmp.hh>
#include <dune/common/fmatrix.hh>
namespace Dune {
namespace FloatCmp {
template<class T, int n, int m>
struct EpsilonType<FieldMatrix<T, n, m> > {
typedef typename EpsilonType<FieldVector<T,n*m>>::Type Type;
};
namespace Impl {
template<class T, int n, int m, CmpStyle cstyle>
struct eq_t_fmat
{
typedef FieldMatrix<T, n, m> M;
static bool eq(const M &first,
const M &second,
typename EpsilonType<M>::Type epsilon = DefaultEpsilon<M>::value())
{
for(int i = 0; i < n; ++i)
for(int j = 0; j < m; ++j)
if(!eq_t<T, cstyle>::eq(first[i][j], second[i][j], epsilon))
return false;
return true;
}
};
template<class T, int n, int m>
struct eq_t<FieldMatrix<T,n,m>,relativeWeak>:eq_t_fmat<T,n,m,relativeWeak> {};
template<class T, int n, int m>
struct eq_t<FieldMatrix<T,n,m>,relativeStrong>:eq_t_fmat<T,n,m,relativeStrong> {};
template<class T, int n, int m>
struct eq_t<FieldMatrix<T,n,m>,absolute>:eq_t_fmat<T,n,m,absolute> {};
template<class I, class T, int n, int m, CmpStyle cstyle, RoundingStyle rstyle>
struct round_t<FieldMatrix<I,n,m>, FieldMatrix<T,n,m>, cstyle, rstyle>
{
static FieldMatrix<I,n,m>
round(const T &val,
typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value()))
{
FieldMatrix<I,n,m> res;
for(int i = 0; i < n; ++i)
for(int j = 0; j < m; ++j)
res[i][j] = round_t<I, T, cstyle, rstyle>::round(val[i][j], epsilon);
return res;
}
};
template<class I, class T, int n, int m, CmpStyle cstyle, RoundingStyle rstyle>
struct trunc_t<FieldMatrix<I,n,m>, FieldMatrix<T,n,m>, cstyle, rstyle>
{
static FieldMatrix<I,n,m>
trunc(const FieldMatrix<T,n,m> &val,
typename EpsilonType<T>::Type epsilon = (DefaultEpsilon<T, cstyle>::value()))
{
FieldMatrix<I,n,m> res;
for(int i = 0; i < n; ++i)
for(int j = 0; j < m; ++j)
res[i][j] = trunc_t<I, T, cstyle, rstyle>::trunc(val[i][j], epsilon);
return res;
}
};
}}} // end namespace FloatCmp::Impl
#pragma once
#include <cassert>
#include <optional>
#include <type_traits>
namespace Dune {
//! LocalFunction associated to the \ref AnalyticGridFunction
template <class LocalContext, class F>
class LocalAnalyticGridFunction;
//! Generator for \ref LocalAnalyticGridFunction
template <class LocalContext, class FF>
auto localAnalyticGridFunction (FF&& ff)
{
using F = std::decay_t<FF>;
return LocalAnalyticGridFunction<LocalContext, F>{std::forward<FF>(ff)};
}
template <class LC, class Functor>
class LocalAnalyticGridFunction
{
using Self = LocalAnalyticGridFunction;
public:
using LocalContext = LC;
using Geometry = typename LocalContext::Geometry;
using Domain = typename Geometry::GlobalCoordinate;
using LocalDomain = typename Geometry::LocalCoordinate;
using Range = std::result_of_t<Functor(Domain)>;
using Signature = Range(LocalDomain);
public:
//! Constructor. Stores the functor f by value
template <class FF,
std::enable_if_t<not std::is_same_v<Self, std::decay_t<FF>>, bool> = true>
LocalAnalyticGridFunction (FF&& f)
: f_(std::forward<FF>(f))
{}
//! bind the LocalFunction to the local context
void bind (const LocalContext& localContext)
{
localContext_.emplace(localContext);
geometry_.emplace(localContext.geometry());
}
//! unbind the localContext from the localFunction
void unbind ()
{
geometry_.reset();
localContext_.reset();
}
//! evaluate the stored function in local coordinates
Range operator() (const LocalDomain& x) const
{
assert(!!geometry_);
return f_(geometry_->global(x));
}
//! return the bound localContext.
const LocalContext& localContext () const
{
assert(!!localContext_);
return *localContext_;
}
//! obtain the functor
const Functor& f () const
{
return f_;
}
private:
Functor f_;
// some caches
std::optional<LocalContext> localContext_;
std::optional<Geometry> geometry_;
};
//! Derivative of a \ref LocalAnalyticGridFunction
template <class LC, class F>
auto derivative (const LocalAnalyticGridFunction<LC,F>& t)
-> decltype(localAnalyticGridFunction<LC>(derivative(t.f())))
{
return localAnalyticGridFunction<LC>(derivative(t.f()));
}
} // end namespace Dune
#pragma once
#include <type_traits>
#include <dune/common/typeutilities.hh>
#include <dune/curvedgeometry/curvedgeometry.hh>
namespace Dune {
template <class K>
FieldVector<K,3> cross(FieldVector<K,3> const& a, FieldVector<K,3> const& b)
{
return {a[1]*b[2] - a[2]*b[1],
a[2]*b[0] - a[0]*b[2],
a[0]*b[1] - a[1]*b[0]};
}
namespace Impl {
template <class Geometry, class LocalCoordinate>
auto normal(Geometry const& geometry, LocalCoordinate const& local, PriorityTag<2>)
-> decltype(geometry.normal(local))
{
return geometry.normal(local);
}
template <class Geometry, class LocalCoordinate>
auto normal(Geometry const& geometry, LocalCoordinate const& local, PriorityTag<1>)
-> decltype(geometry.impl().normal(local))
{
return geometry.impl().normal(local);
}
template <class Geometry, class LocalCoordinate>
typename Geometry::GlobalCoordinate normal(Geometry const& geometry, LocalCoordinate const& local, PriorityTag<0>)
{
auto edge1 = geometry.corner(0) - geometry.corner(1);
auto edge2 = geometry.corner(0) - geometry.corner(2);
auto n = cross(edge1, edge2);
return n / n.two_norm();
}
} // end namespace Impl
template <class Geometry, class LocalCoordinate>
auto normal(Geometry const& geometry, LocalCoordinate const& local)
{
return Impl::normal(geometry, local, PriorityTag<3>{});
}
} // end namespace Dune
#include <config.h>
#include <dune/common/fmatrixev.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/geometry/affinegeometry.hh>
#include <dune/curvedgeometry/parametrizedgeometry.hh>
#include <dune/curvedgeometry/localfunctiongeometry.hh>
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
#include <dune/localfunctions/lagrange/lfecache.hh>
#include "floatcmp.hh"
#include "localanalyticgridfunction.hh"
#include "normals.hh"
using namespace Dune;
using LocalCoordinate = FieldVector<double,2>;
using GlobalCoordinate = FieldVector<double,3>;
/// Functor representing a sphere
struct SphereProjection
{
using Jacobian = FieldMatrix<double,3,3>;
/// project the coordinate to the sphere at origin with `radius`
GlobalCoordinate operator() (const GlobalCoordinate& x) const
{
return x / x.two_norm();
}
/// derivative of the projection
friend auto derivative (const SphereProjection& sphere)
{
return [](const GlobalCoordinate& x) -> Jacobian
{
Jacobian out;
auto nrm = x.two_norm();
for (int i = 0; i < 3; ++i)
for (int j = 0; j < 3; ++j)
out[i][j] = ((i == j ? 1 : 0) - (x[i]/nrm) * (x[j]/nrm)) / nrm;
return out;
};
}
};
template <class Geo>
struct RefElement
{
using Geometry = Geo;
RefElement(Geometry const& geometry)
: geometry_(&geometry)
{}
Geometry const& geometry() const
{
return *geometry_;
}
GeometryType type() const
{
return geometry_->type();
}
Geometry const* geometry_;
};
// convert the value into a string
template <class T>
std::string to_str(T const& value)
{
std::ostringstream oss;
oss << value;
return oss.str();
}
// macro to test for equality of two objects
#define EXPECT_EQ(A,B) \
test.check(FloatCmp::eq<std::decay_t<decltype(A)>, FloatCmp::absolute>(A,B), std::string(#A " != " #B ": ") + to_str(A) + " != " + to_str(B));
// projection matrix
template <class T, int m>
FieldMatrix<T,m,m> projection_mat(FieldVector<T,m> const& n)
{
FieldMatrix<T,m,m> P;
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
P[i][j] = (i == j ? 1 : 0) - n[i]*n[j];
return P;
}
// projection matrix
template <class T, int m>
FieldMatrix<T,m,m> curvature_mat(FieldVector<T,m> const& x)
{
FieldMatrix<T,m,m> H;
T nrm = x.two_norm();
for (int i = 0; i < m; ++i)
for (int j = 0; j < m; ++j)
H[i][j] = ((i == j ? 1 : 0) - x[i]/nrm*x[j]/nrm)/nrm;
return H;
}
using MatrixHelper = Impl::FieldMatrixHelper<double>;
template <int order, bool flat = false, class FlatGeometry>
void run_test(TestSuite& testSuite, FlatGeometry const& flatGeometry)
{
SphereProjection projection;
TestSuite test{"order" + std::to_string(order)};
using LFECache = LagrangeLFECache<double,double,2>;
LFECache lfeCache(order);
using Element = RefElement<FlatGeometry>;
Element element0(flatGeometry);
// local finite-elements
using LFE1 = LagrangeSimplexLocalFiniteElement<double,double,2,order>;
LFE1 lagrange1;
using LFE2 = typename LFECache::FiniteElementType;
LFE2 lagrange2 = lfeCache.get(GeometryTypes::triangle);
// local function defined on the flat elements
using LF0 = LocalAnalyticGridFunction<Element,SphereProjection>;
LF0 localFct0{projection};
localFct0.bind(element0);
// geometries defined on the flat element
using ExactGeometry0 = ElementLocalFunctionGeometry<LF0,double,2>;
ExactGeometry0 exactGeometry0{GeometryTypes::triangle, localFct0};
auto curvedGeometry1 = [&]{
if constexpr (order == 1 && flat) {
return flatGeometry;
} else {
return ParametrizedGeometry<LFE1,3>{GeometryTypes::triangle, lagrange1, localFct0};
}
}();
using CurvedGeometry1 = decltype(curvedGeometry1);
using CurvedGeometry2 = ParametrizedGeometry<LFE2,3>;
CurvedGeometry2 curvedGeometry2{GeometryTypes::triangle, lagrange2, localFct0};
using CurvedElement1 = RefElement<CurvedGeometry1>;
using CurvedElement2 = RefElement<CurvedGeometry2>;
CurvedElement1 curvedElement1{curvedGeometry1};
CurvedElement2 curvedElement2{curvedGeometry2};
// local functions defined on the curved elements
using LF1 = LocalAnalyticGridFunction<CurvedElement1,SphereProjection>;
LF1 localFct1{projection};
localFct1.bind(curvedElement1);
using LF2 = LocalAnalyticGridFunction<CurvedElement2,SphereProjection>;
LF2 localFct2{projection};
localFct2.bind(curvedElement2);
// geometries defined on the curved elements
using ExactGeometry1 = ElementLocalFunctionGeometry<LF1,double,2>;
using ExactGeometry2 = ElementLocalFunctionGeometry<LF2,double,2>;
ExactGeometry1 exactGeometry1{GeometryTypes::triangle, localFct1};
ExactGeometry2 exactGeometry2{GeometryTypes::triangle, localFct2};
// ----------------------------------------------------
FieldMatrix<double,2,2> I2{{1,0},{0,1}};
FieldMatrix<double,3,3> I3{{1,0,0},{0,1,0},{0,0,1}};
FieldVector<double,2> v{1.0, 1.0};
const auto& quadRule = QuadratureRules<double,2>::rule(GeometryTypes::triangle, 6);
for (const auto& qp : quadRule)
{
// 1. test whether local-to-global mapping is correct
// --------------------------------------------------
auto global0 = flatGeometry.global(qp.position());
auto global1 = curvedGeometry1.global(qp.position());
auto global2 = curvedGeometry2.global(qp.position());
EXPECT_EQ(global1, global2);
auto global = exactGeometry0.global(qp.position());
EXPECT_EQ(global, projection(global0));
EXPECT_EQ(global, projection(global));
// 2. test whether global-to-local mapping is correct
// --------------------------------------------------
auto local0 = flatGeometry.local(global0);
auto local1 = curvedGeometry1.local(global1);
EXPECT_EQ(local1, qp.position());
auto local2 = curvedGeometry2.local(global2);
EXPECT_EQ(local2, qp.position());
auto local = exactGeometry0.local(global);
EXPECT_EQ(local, qp.position());
// 3. test whether Jacobian is correct
// -----------------------------------
auto Jt0 = flatGeometry.jacobianTransposed(qp.position());
auto Jt1 = curvedGeometry1.jacobianTransposed(qp.position());
auto Jt2 = curvedGeometry2.jacobianTransposed(qp.position());
EXPECT_EQ(Jt1, Jt2);
auto Jt = exactGeometry0.jacobianTransposed(qp.position());
// 3.1 transform vectors and check whether tangential
// --------------------------------------------------
FieldVector<double,3> V1,V2,V;
Jt1.mtv(v, V1);
Jt2.mtv(v, V2);
Jt.mtv(v, V);
auto nh0 = normal(flatGeometry, qp.position());
auto nh1 = normal(curvedGeometry1, qp.position());
EXPECT_EQ(V1.dot(nh1), 0.0);
auto nh2 = curvedGeometry2.normal(qp.position());
EXPECT_EQ(V2.dot(nh2), 0.0);
auto n0 = exactGeometry0.normal(qp.position());
auto n1 = exactGeometry1.normal(qp.position());
auto n2 = exactGeometry2.normal(qp.position());
EXPECT_EQ(V.dot(n0), 0.0);
// 4. test the integration element
// -------------------------------
auto dS0 = flatGeometry.integrationElement(qp.position());
auto dS1 = curvedGeometry1.integrationElement(qp.position());
auto dS2 = curvedGeometry2.integrationElement(qp.position());
EXPECT_EQ(dS1, dS2);
auto dS = exactGeometry0.integrationElement(qp.position());
EXPECT_EQ(dS, MatrixHelper::sqrtDetAAT(Jt));
// 4.1 check dSi is product of singular values
// -----