Skip to content
Snippets Groups Projects
Commit 3e63e4da authored by Simon Praetorius's avatar Simon Praetorius
Browse files

updated the test

parent 315ab02a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/curvedgeometry/curvedgeometry.hh
+ 92
− 88
View file @ 3e63e4da
......@@ -10,6 +10,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/std/type_traits.hh>
#include <dune/geometry/affinegeometry.hh>
#include <dune/geometry/quadraturerules.hh>
......@@ -22,33 +23,6 @@ namespace Dune
{
namespace Impl
{
template <class T, int n, int m>
FieldMatrix<T,n,m> outerProduct(const FieldVector<T,n>& a, const FieldVector<T,m>& b)
{
FieldMatrix<T,n,m> res;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
res[i][j] = a[i] * b[j];
return res;
}
template <class T, int n, int m>
FieldMatrix<T,n,m> outerProduct(const FieldMatrix<T,n,1>& a, const FieldVector<T,m>& b)
{
FieldMatrix<T,n,m> res;
for (int i = 0; i < n; ++i)
for (int j = 0; j < m; ++j)
res[i][j] = a[i][0] * b[j];
return res;
}
template <class T, int n, int m,
std::enable_if_t<(n > 1), int> = 0>
FieldMatrix<T,n,m> outerProduct(const FieldMatrix<T,1,n>& a, const FieldVector<T,m>& b)
{
return outerProduct(a[0], b);
}
template <class T, int n, int m>
void outerProductAccumulate(const FieldVector<T,n>& a, const FieldVector<T,m>& b, FieldMatrix<T,n,m>& res)
{
......@@ -85,8 +59,9 @@ namespace Dune
* This structure provides the default values.
*
* \tparam ct coordinate type
* \tparam LFECache A LocalFiniteElementVariantCache implementation, e.g. LagrangeLocalFiniteElementCache
*/
template <class ct>
template <class ct, class LFECache>
struct CurvedGeometryTraits
{
/// \brief helper structure containing some matrix routines. See affinegeometry.hh
......@@ -98,26 +73,7 @@ namespace Dune
/// \brief maximal number of Newton iteration in `geometry.local(global)`
static int maxIteration () { return 100; }
/// \brief template specifying the storage for the vertices
/**
* Internally, the CurvedGeometry needs to store the lagrange vertices of the
* geometry.
*
* \tparam coorddim coordinate dimension
*/
template <int coorddim>
struct VertexStorage
{
using type = std::vector<FieldVector<ct, coorddim>>;
};
/// \brief will there be only one geometry type for a dimension?
template <int dim>
struct hasSingleGeometryType
{
static const bool v = false;
static const unsigned int topologyId = ~0u;
};
using LocalFiniteElementCache = LFECache;
};
......@@ -131,28 +87,23 @@ namespace Dune
* \tparam ct coordinate type
* \tparam mydim geometry dimension
* \tparam cdim coordinate dimension
* \tparam polynomial_order
* \tparam Traits Parametrization of the geometry, see \ref CurvedGeometryTraits
*
* The requirements on the traits are documented along with their default,
* CurvedGeometryTraits.
*/
template <class ct, int mydim, int cdim, int polynomial_order>
template <class ct, int mydim, int cdim, class Traits>
class CurvedGeometry
{
public:
/// coordinate type
using ctype = ct;
using Traits = CurvedGeometryTraits<ctype>;
/// geometry dimension
static const int mydimension = mydim;
/// coordinate dimension
static const int coorddimension = cdim;
/// Polynomial order of geometry
static const int order = polynomial_order;
/// type of local coordinates
using LocalCoordinate = FieldVector<ctype, mydimension>;
/// type of global coordinates
......@@ -176,43 +127,91 @@ namespace Dune
protected:
using MatrixHelper = typename Traits::MatrixHelper;
using LocalFECache = LagrangeLocalFiniteElementCache<ctype, ctype, mydimension, order>;
using LocalFECache = typename Traits::LocalFiniteElementCache;
using LocalFiniteElement = typename LocalFECache::FiniteElementType;
using LocalBasis = typename LocalFiniteElement::Traits::LocalBasisType;
using LocalBasisTraits = typename LocalBasis::Traits;
protected:
CurvedGeometry (const ReferenceElement& refElement)
: refElement_(refElement)
, localFECache_()
, localFE_(localFECache_.get(refElement.type()))
{}
public:
/// \brief constructor
/**
* \param[in] refElement reference element for the geometry
* \param[in] corners corners to store internally
* \param[in] vertices vertices to store internally
*
* \note The type of vertices is actually a template argument.
* It is only required that the internal vertex storage can be
* constructed from this object.
*/
template <class Vertices>
CurvedGeometry (const ReferenceElement& refElement, const Vertices& vertices)
: refElement_(refElement)
, vertices_(vertices)
, localFECache_()
, localBasis_(localFECache_.get(refElement.type()).localBasis())
{}
CurvedGeometry (const ReferenceElement& refElement, std::vector<GlobalCoordinate> vertices)
: CurvedGeometry(refElement)
{
vertices_ = std::move(vertices);
assert(localFE_.size() == vertices_.size());
}
template <class Parametrization,
std::enable_if_t<Std::is_callable<Parametrization(LocalCoordinate), GlobalCoordinate>::value, bool> = true>
CurvedGeometry (const ReferenceElement& refElement, Parametrization&& param)
: CurvedGeometry(refElement)
{
auto const& localInterpolation = localFE_.localInterpolation();
localInterpolation.interpolate(param, vertices_);
}
/// \brief constructor
/**
* \param[in] gt geometry type
* \param[in] corners corners to store internally
* \param[in] gt geometry type
* \param[in] param either a vector of vertices, or a functor that can be used to construct the vertices
*/
template <class Vertices>
CurvedGeometry (Dune::GeometryType gt, const Vertices& vertices)
: CurvedGeometry(ReferenceElements::general(gt), vertices)
template <class Parametrization>
CurvedGeometry (Dune::GeometryType gt, Parametrization&& param)
: CurvedGeometry(ReferenceElements::general(gt), std::forward<Parametrization>(param))
{}
/// \brief Copy constructor
CurvedGeometry (const CurvedGeometry& that)
: CurvedGeometry(that.refElement_, that.vertices_)
{}
/// \brief Move constructor
CurvedGeometry (CurvedGeometry&& that)
: CurvedGeometry(that.refElement_, std::move(that.vertices_))
{}
/// \brief Copy assignment operator
CurvedGeometry& operator=(const CurvedGeometry& that)
{
assert(refElement_ == that.refElement_);
vertices_ = that.vertices_;
return *this;
}
/// \brief Move assignment operator
CurvedGeometry& operator=(CurvedGeometry&& that)
{
assert(refElement_ == that.refElement_);
vertices_ = std::move(that.vertices_);
return *this;
}
/// \brief obtain the polynomial order of the parametrization
int order () const
{
return localBasis().order();
}
/// \brief is this mapping affine?
bool affine () const
{
return refElement_.type().isSimplex() && int(vertices_.size()) == corners();
return refElement_.type().isSimplex() && order() == 1;
}
/// \brief obtain the name of the reference element
......@@ -221,7 +220,7 @@ namespace Dune
return refElement_.type();
}
/// \brief obtain number of corners of the corresponding reference element */
/// \brief obtain number of corners of the corresponding reference element
int corners () const
{
return refElement_.size(mydimension);
......@@ -257,11 +256,8 @@ namespace Dune
*/
GlobalCoordinate global (const LocalCoordinate& local) const
{
if (mydimension == 0)
return vertices_[0];
thread_local std::vector<typename LocalBasisTraits::RangeType> shapeValues;
localBasis_.evaluateFunction(local, shapeValues);
localBasis().evaluateFunction(local, shapeValues);
assert(shapeValues.size() == vertices_.size());
GlobalCoordinate out(0);
......@@ -293,12 +289,11 @@ namespace Dune
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 = (*this).global( x ) - globalCoord;
const bool invertible =
MatrixHelper::template xTRightInvA<mydimension, coorddimension>(jacobianTransposed(x), dglobal, dx);
const GlobalCoordinate dglobal = global(x) - globalCoord;
const bool invertible = MatrixHelper::xTRightInvA(jacobianTransposed(x), dglobal, dx);
if (!invertible)
return LocalCoordinate(std::numeric_limits<ctype>::max());
break;
// update x with correction
x -= dx;
......@@ -306,6 +301,7 @@ namespace Dune
// for affine mappings only one iteration is needed
if (affineMapping)
break;
if (dx.two_norm2() < tolerance)
break;
}
......@@ -343,7 +339,7 @@ namespace Dune
*/
ctype integrationElement (const LocalCoordinate& local) const
{
return MatrixHelper::template sqrtDetAAT<mydimension, coorddimension>(jacobianTransposed(local));
return MatrixHelper::sqrtDetAAT(jacobianTransposed(local));
}
/// \brief Obtain the volume of the mapping's image
......@@ -352,8 +348,8 @@ namespace Dune
*/
Volume volume () const
{
int p = 2*order; // TODO: needs to be checked
auto const& quadRule = QuadratureRules<ctype, mydimension>::rule(type(), p);
const int p = 2*localBasis().order(); // TODO: needs to be checked
const auto& quadRule = QuadratureRules<ctype, mydimension>::rule(type(), p);
Volume vol(0);
for (auto const& qp : quadRule)
vol += integrationElement(qp.position()) * qp.weight();
......@@ -373,7 +369,7 @@ namespace Dune
JacobianTransposed jacobianTransposed (const LocalCoordinate& local) const
{
std::vector<typename LocalBasisTraits::JacobianType> shapeJacobians;
localBasis_.evaluateJacobian(local, shapeJacobians);
localBasis().evaluateJacobian(local, shapeJacobians);
assert(shapeJacobians.size() == vertices_.size());
JacobianTransposed out(0);
......@@ -401,7 +397,7 @@ namespace Dune
return geometry.refElement();
}
auto const& vertices () const
const std::vector<GlobalCoordinate>& vertices () const
{
return vertices_;
}
......@@ -412,6 +408,11 @@ namespace Dune
return refElement_;
}
const LocalBasis& localBasis() const
{
return localFE_.localBasis();
}
GlobalCoordinate normalEdge (const LocalCoordinate& local, const JacobianTransposed& J) const
{
GlobalCoordinate res{
......@@ -433,13 +434,16 @@ namespace Dune
private:
ReferenceElement refElement_;
std::vector<GlobalCoordinate> vertices_;
//typename Traits::template VertexStorage<coorddimension>::type vertices_;
LocalFECache localFECache_;
LocalBasis const& localBasis_;
LocalFiniteElement const& localFE_;
std::vector<GlobalCoordinate> vertices_;
};
template <class ctype, int mydim, int cdim, int order>
using LagrangeCurvedGeometry = CurvedGeometry<ctype,mydim,cdim,
CurvedGeometryTraits<ctype, LagrangeLocalFiniteElementCache<ctype, ctype, mydim, order>> >;
} // namespace Dune
......
src/dune-curvedgeometry.cc
+ 15
− 8
View file @ 3e63e4da
......@@ -9,6 +9,7 @@
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/curvedgeometry/curvedgeometry.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/yaspgrid.hh>
......@@ -21,13 +22,14 @@ int main(int argc, char** argv)
using namespace Dune;
MPIHelper& helper = MPIHelper::instance(argc, argv);
YaspGrid<2> grid({1.0,1.0}, {4,4});
YaspGrid<2> grid({8.0,8.0}, {32,32});
using LocalCoordinate = FieldVector<double,2>;
using GlobalCoordinate = FieldVector<double,2>;
using WorldCoordinate = FieldVector<double,3>;
using LocalFECache = LagrangeLocalFiniteElementCache<double, double, 2, order>;
LocalFECache localFeCache;
using Geometry = LagrangeCurvedGeometry<double, 2, 3, order>;
FieldMatrix<double,2,2> I{{1,0},{0,1}};
// coordinate projection
auto project = [](GlobalCoordinate const& global) -> WorldCoordinate
......@@ -37,22 +39,27 @@ int main(int argc, char** argv)
std::sin(global[0])*std::cos(global[1]) };
};
TestSuite test("curved geometry");
for (auto const& e : elements(grid.leafGridView()))
{
// projection from local coordinates
auto X = [project,geo=e.geometry()](LocalCoordinate const& local) -> WorldCoordinate { return project(geo.global(local)); };
// construct lagrange vertices
std::vector<WorldCoordinate> vertices;
localFeCache.get(e.type()).localInterpolation().interpolate(X, vertices);
// create a curved geometry
CurvedGeometry<double,2,3,order> geometry(e.type(), vertices);
Geometry geometry(e.type(), X);
auto const& quadRule = QuadratureRules<double,2>::rule(e.type(), 4);
for (auto const& qp : quadRule) {
auto Jt = geometry.jacobianTransposed(qp.position());
auto Jtinv = geometry.jacobianInverseTransposed(qp.position());
FieldMatrix<double, 2, 2> res;
FMatrixHelp::multMatrix(Jt, Jtinv, res);
res -= I;
test.check(res.frobenius_norm() < std::sqrt(std::numeric_limits<double>::epsilon()), "J^-1 * J == I");
}
}
return test.report() ? 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