Commit 38f79c23 authored by Simon Praetorius's avatar Simon Praetorius
Browse files

Merge branch 'issue/changes-in-geometry-and-localfunctions' into 'master'

Update LagrangePoints due to changes in dune-geometry and dune-localfunctions

See merge request !30
parents cb00bc20 a9022011
Pipeline #37175 passed with stage
in 13 minutes and 43 seconds
......@@ -226,7 +226,11 @@ namespace Dune
int order (GeometryType type, std::size_t nNodes) const
{
for (int o = 1; o <= int(nNodes); ++o)
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
if (numLagrangePoints(type.id(), type.dim(), o) == std::size_t(nNodes))
#else
if (numLagrangePoints(type, o) == std::size_t(nNodes))
#endif
return o;
return 1;
......
......@@ -4,6 +4,7 @@
#include <vector>
#include <dune/common/dynvector.hh>
#include <dune/common/version.hh>
#include <dune/localfunctions/lagrange.hh>
#include <dune/vtk/gridfunctions/common.hh>
#include <dune/vtk/gridfunctions/continuousgridfunction.hh>
......@@ -83,7 +84,11 @@ namespace Dune
std::int64_t shift = (insertionIndex == 0 ? 0 : (*offsets_)[insertionIndex-1]);
std::int64_t numNodes = (*offsets_)[insertionIndex] - shift;
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
[[maybe_unused]] std::int64_t maxNumNodes = numLagrangePoints(element.type().id(), element.type().dim(), 20);
#else
[[maybe_unused]] std::int64_t maxNumNodes = numLagrangePoints(element.type(), 20);
#endif
VTK_ASSERT(numNodes > 0 && numNodes < maxNumNodes);
int order = creator_->order(element.type(), numNodes);
......
......@@ -4,6 +4,7 @@
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
......@@ -47,16 +48,29 @@ namespace Dune
}
/// Fill the lagrange points for the given topology type `Topology`
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
template <class Topology>
bool build ()
{
build(GeometryType(Topology{}));
return true;
}
#else
template <GeometryType::Id geometryId>
bool build ()
{
build(GeometryType(geometryId));
return true;
}
#endif
/// Returns whether the point set support the given topology type `Topology` and can
/// generate point for the given order.
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
template <class Topology>
#else
template <GeometryType::Id geometryId>
#endif
static bool supports (std::size_t order)
{
return true;
......
......@@ -4,6 +4,7 @@
#include <array>
#include <dune/common/exceptions.hh>
#include <dune/common/version.hh>
#include <dune/geometry/type.hh>
#include <dune/localfunctions/lagrange/equidistantpoints.hh>
......@@ -89,7 +90,11 @@ template <class K>
template <class Points>
void LagrangePointSetBuilder<K,2>::operator() (GeometryType gt, int order, Points& points) const
{
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
#else
std::size_t nPoints = numLagrangePoints(gt, order);
#endif
if (gt.isTriangle())
buildTriangle(nPoints, order, points);
......@@ -200,8 +205,8 @@ void LagrangePointSetBuilder<K,2>::buildQuad(std::size_t nPoints, int order, Poi
const K r = K(m) / orders[0];
const K s = K(n) / orders[1];
Vec p = (1.0 - r) * (nodes[3] * s + nodes[0] * (1.0 - s)) +
r * (nodes[2] * s + nodes[1] * (1.0 - s));
Vec p = K(1.0-r) * (nodes[3] * s + nodes[0] * K(1.0 - s)) +
r * (nodes[2] * s + nodes[1] * K(1.0 - s));
auto [idx,key] = calcQuadKey(m,n,orders);
points[idx] = LP{p, key};
......@@ -262,7 +267,11 @@ template <class K>
template <class Points>
void LagrangePointSetBuilder<K,3>::operator() (GeometryType gt, unsigned int order, Points& points) const
{
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
std::size_t nPoints = numLagrangePoints(gt.id(), dim, order);
#else
std::size_t nPoints = numLagrangePoints(gt, order);
#endif
if (gt.isTetrahedron())
buildTetra(nPoints, order, points);
......@@ -410,8 +419,8 @@ void LagrangePointSetBuilder<K,3>::buildHex (std::size_t nPoints, int order, Poi
const K r = K(m) / orders[0];
const K s = K(n) / orders[1];
const K t = K(p) / orders[2];
Vec point = (1.0-r) * ((nodes[3] * (1.0-t) + nodes[7] * t) * s + (nodes[0] * (1.0-t) + nodes[4] * t) * (1.0-s)) +
r * ((nodes[2] * (1.0-t) + nodes[6] * t) * s + (nodes[1] * (1.0-t) + nodes[5] * t) * (1.0-s));
Vec point = K(1.0-r) * ((nodes[3] * K(1.0-t) + nodes[7] * t) * s + (nodes[0] * K(1.0-t) + nodes[4] * t) * K(1.0-s)) +
r * ((nodes[2] * K(1.0-t) + nodes[6] * t) * s + (nodes[1] * K(1.0-t) + nodes[5] * t) * K(1.0-s));
auto [idx,key] = calcHexKey(m,n,p,orders);
points[idx] = LP{point, key};
......
......@@ -2,6 +2,7 @@
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <dune/common/version.hh>
#include <dune/localfunctions/utility/field.hh>
#include <dune/localfunctions/utility/basisprint.hh>
......@@ -72,24 +73,43 @@ bool test(const Basis &basis, const Points &points, bool verbose)
return ret;
}
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
template <class Topology>
#else
template <Dune::GeometryType::Id geometryId>
#endif
bool test(unsigned int order, bool verbose = false)
{
typedef Dune::LagrangeBasisFactory<Dune::Vtk::LagrangePointSet,Topology::dimension,StorageField,ComputeField> BasisFactory;
typedef Dune::LagrangeCoefficientsFactory< Dune::Vtk::LagrangePointSet, Topology::dimension,double > LagrangeCoefficientsFactory;
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
typedef Dune::LagrangeBasisFactory<Dune::Vtk::LagrangePointSet, Topology::dimension, StorageField, ComputeField> BasisFactory;
typedef Dune::LagrangeCoefficientsFactory<Dune::Vtk::LagrangePointSet, Topology::dimension, double > LagrangeCoefficientsFactory;
#else
constexpr Dune::GeometryType geometry = geometryId;
typedef Dune::LagrangeBasisFactory<Dune::Vtk::LagrangePointSet, geometry.dim(), StorageField, ComputeField> BasisFactory;
typedef Dune::LagrangeCoefficientsFactory<Dune::Vtk::LagrangePointSet, geometry.dim(), double> LagrangeCoefficientsFactory;
#endif
bool ret = true;
for (unsigned int o = 0; o <= order; ++o)
{
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< Topology >( o );
#else
const typename LagrangeCoefficientsFactory::Object *pointsPtr = LagrangeCoefficientsFactory::template create< geometry >( o );
#endif
if ( pointsPtr == 0)
continue;
std::cout << "# Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl;
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
std::cout << "Testing " << Topology::name() << " in dimension " << Topology::dimension << " with order " << o << std::endl;
typename BasisFactory::Object &basis = *BasisFactory::template create<Topology>(o);
#else
std::cout << "Testing " << geometry << " with order " << o << std::endl;
typename BasisFactory::Object &basis = *BasisFactory::template create<geometry>(o);
#endif
ret |= test(basis,*pointsPtr,verbose);
......@@ -97,10 +117,17 @@ bool test(unsigned int order, bool verbose = false)
// derivatives in a human readabible form (aka LaTeX source)
#ifdef TEST_OUTPUT_FUNCTIONS
std::stringstream name;
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
name << "lagrange_" << Topology::name() << "_p" << o << ".basis";
std::ofstream out(name.str().c_str());
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField>(out,basis);
#else
name << "lagrange_" << geometry << "_p" << o << ".basis";
std::ofstream out(name.str().c_str());
Dune::basisPrint<0,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
Dune::basisPrint<1,BasisFactory,typename BasisFactory::StorageField,geometry>(out,basis);
#endif
#endif // TEST_OUTPUT_FUNCTIONS
LagrangeCoefficientsFactory::release( pointsPtr );
......@@ -149,19 +176,36 @@ int main ( int argc, char **argv )
bool tests = true;
#ifdef CHECKDIM1
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
tests &= test<Prism<Point> > (order); // line
tests &= test<Pyramid<Point> > (order); // line
#else
tests &= test<GeometryTypes::cube(1)>(order); // line
tests &= test<GeometryTypes::simplex(1)>(order); // line
#endif
#endif
#ifdef CHECKDIM2
//tests &= test<Prism<Prism<Point> > > (order); // quad
//tests &= test<Pyramid<Pyramid<Point> > >(order); // triangle
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
tests &= test<Prism<Prism<Point> > > (order); // quad
tests &= test<Pyramid<Pyramid<Point> > >(order); // triangle
#else
tests &= test<GeometryTypes::cube(2)>(order); // quad
tests &= test<GeometryTypes::simplex(2)>(order); // triangle
#endif
#endif
#ifdef CHECKDIM3
// tests &= test<Prism<Prism<Prism<Point> > > >(order); // hexahedron
#if DUNE_VERSION_LT(DUNE_LOCALFUNCTIONS,2,8)
tests &= test<Prism<Prism<Prism<Point> > > >(order); // hexahedron
// tests &= test<Prism<Pyramid<Pyramid<Point> > > >(order);
// tests &= test<Pyramid<Pyramid<Pyramid<Point> > > >(order); // tetrahedron
tests &= test<Pyramid<Pyramid<Pyramid<Point> > > >(order); // tetrahedron
#else
tests &= test<GeometryTypes::cube(3)>(order); // hexahedron
// tests &= test<GeometryTypes::prism>(order);
// tests &= test<GeometryTypes::pyramid>(order);
tests &= test<GeometryTypes::simplex(3)>(order);
#endif
#endif
return (tests ? 0 : 1);
......
Markdown is supported
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