Commit 08c8e3c6 authored by graeser's avatar graeser

Merge branch 'feature/update-refinedsimplexgeometry' into 'master'

Modernize makeRefinedSimplexGeometry implementation and test

See merge request !68
parents 7a3c7508 408d93d3
Pipeline #28151 canceled with stage
in 5 minutes
......@@ -8,6 +8,10 @@
#include <dune/geometry/type.hh>
#include <dune/geometry/affinegeometry.hh>
// Borrow internal helper from dune-localfunctions
// for transition of function interface
#include <dune/localfunctions/common/localinterpolation.hh>
template <class ct, int localdim, int worlddim>
class RefinedSimplexGeometry {
public:
......@@ -258,12 +262,22 @@ class RefinedSimplexGeometry<ct,2,3> {
}
};
template <class Element, class Parameterization>
RefinedSimplexGeometry<typename Element::Geometry::ctype,2,3> makeRefinedSimplexGeometry(const Element& element, const Parameterization& parameterization)
/**
* \brief Create RefinedSimplexGeometry with transformed geometry
*
* \param element The element whose geometry should be transformed
* \param transform The transformation mapping global element coordinated to transformed global coordinates
*
* This will only transform edge mid points but leave the vertices in place.
*/
template <class Element, class Transformation>
RefinedSimplexGeometry<typename Element::Geometry::ctype,2,3> makeRefinedSimplexGeometry(const Element& element, const Transformation& transformation)
{
if(not element.type().isTriangle())
DUNE_THROW(Dune::NotImplemented, "Free method makeRefinedSimplexGeometry(.,.) is only implemented for 2D simplex elements");
using GlobalCoordinate = typename Element::Geometry::GlobalCoordinate;
auto&& transformationF = Dune::Impl::makeFunctionWithCallOperator<GlobalCoordinate>(transformation);
std::array<typename Element::Geometry::GlobalCoordinate,6> movedVirtualVertices;
typename Element::Geometry eltGeometry = element.geometry();
......@@ -273,8 +287,7 @@ RefinedSimplexGeometry<typename Element::Geometry::ctype,2,3> makeRefinedSimplex
for (std::size_t i=0; i<(std::size_t)element.subEntities(Element::dimension - 1); ++i)
{
typename Element::Geometry::GlobalCoordinate y;
parameterization.evaluate(element.template subEntity<Element::dimension -1>(i).geometry().center(), y);
auto y = transformationF(element.template subEntity<Element::dimension -1>(i).geometry().center());
movedVirtualVertices[eltGeometry.corners() + i] = y;
}
......
......@@ -6,6 +6,7 @@
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/test/testsuite.hh>
#if HAVE_DUNE_ALUGRID
#include <dune/grid/common/gridfactory.hh>
......@@ -25,21 +26,12 @@ class MakeRefinedSimplexGeometryTestSuite
{}
bool check()
Dune::TestSuite test()
{
return check_makeRefinedSimplexGeometry();
}
private:
bool check_makeRefinedSimplexGeometry()
{
bool passed = true;
Dune::TestSuite test;
#if HAVE_DUNE_ALUGRID
typedef Dune::ALUGrid<2,3,Dune::simplex, Dune::conforming> GridType;
GridType* grid;
Dune::GridFactory<GridType> factory;
std::vector<WorldCoords> vertices(3);
......@@ -61,66 +53,55 @@ class MakeRefinedSimplexGeometryTestSuite
factory.insertElement(Dune::GeometryType(Dune::GeometryType::simplex,2),face);
}
grid = factory.createGrid();
auto grid = factory.createGrid();
auto eltIt = grid->leafGridView().begin<0>();
MapToSphere<WorldCoords,WorldCoords> parameterization;
auto transformation = [] (const WorldCoords& x) {
WorldCoords y = x;
y /= x.two_norm();
return y;
};
RefinedSimplexGeometry<typename GridType::ctype, 2,3> refinedGeometry = makeRefinedSimplexGeometry(*eltIt,parameterization);
RefinedSimplexGeometry<typename GridType::ctype, 2,3> refinedGeometry = makeRefinedSimplexGeometry(*eltIt,transformation);
{
WorldCoords testPoint{-sqrt_two/4, 0 , sqrt_two/2}, sollResult(testPoint);
sollResult /= testPoint.two_norm();
passed = passed and
isCloseDune(sollResult,
refinedGeometry.global(LocalCoords{0, 0.5}));
auto xLocal = LocalCoords{0.0, 0.5};
auto xGlobal = eltIt->geometry().global(xLocal);
auto xMappedGlobal = transformation(xGlobal);
test.check(isCloseDune(xMappedGlobal, refinedGeometry.global(xLocal)))
<< "Edge center (" << xLocal << ") has wrong global coordinate. Should be (" << xMappedGlobal << ") but is (" << refinedGeometry.global(xLocal) << ").";
}
{
auto xLocal = LocalCoords{0.5, 0.5};
auto xGlobal = eltIt->geometry().global(xLocal);
auto xMappedGlobal = transformation(xGlobal);
test.check(isCloseDune(xMappedGlobal, refinedGeometry.global(xLocal)))
<< "Edge center (" << xLocal << ") has wrong global coordinate. Should be (" << xMappedGlobal << ") but is (" << refinedGeometry.global(xLocal) << ").";
}
{
WorldCoords testPoint{sqrt_two/8, sqrt_six/8 , sqrt_two/2}, sollResult(testPoint);
sollResult /= testPoint.two_norm();
passed = passed and
isCloseDune(sollResult,
refinedGeometry.global(LocalCoords{0.5, 0.5}));
auto xLocal = LocalCoords{0.5, 0.0};
auto xGlobal = eltIt->geometry().global(xLocal);
auto xMappedGlobal = transformation(xGlobal);
test.check(isCloseDune(xMappedGlobal, refinedGeometry.global(xLocal)))
<< "Edge center (" << xLocal << ") has wrong global coordinate. Should be (" << xMappedGlobal << ") but is (" << refinedGeometry.global(xLocal) << ").";
}
{
WorldCoords testPoint{sqrt_two/8, -sqrt_six/8 , sqrt_two/2}, sollResult(testPoint);
sollResult /= testPoint.two_norm();
passed = passed and
isCloseDune(sollResult,
refinedGeometry.global(LocalCoords{0.5, 0}));
WorldCoords testPoint2{0, 0 , sqrt_two/2}, sollResult2(sollResult);
sollResult2[0] = sollResult2[1] = 0;
passed = passed and isCloseDune(sollResult2,
refinedGeometry.global(LocalCoords{
1.0 / 3.0, 1.0 / 3.0}));
auto xLocal = LocalCoords{1.0/3.0, 1.0/3.0};
// Center of element should be the center of refined sub-element
// spanned by the transfomed edge midpoints.
auto xMappedGlobal = WorldCoords{0.0,0.0,0.0};
xMappedGlobal += transformation(eltIt->geometry().global(LocalCoords{0.0, 0.5}));
xMappedGlobal += transformation(eltIt->geometry().global(LocalCoords{0.5, 0.5}));
xMappedGlobal += transformation(eltIt->geometry().global(LocalCoords{0.5, 0.0}));
xMappedGlobal /= 3.0;
test.check(isCloseDune(xMappedGlobal, refinedGeometry.global(xLocal)))
<< "Element center (" << xLocal << ") has wrong global coordinate. Should be (" << xMappedGlobal << ") but is (" << refinedGeometry.global(xLocal) << ").";
}
#endif
return passed;
return test;
}
template <class DomainVector, class RangeVector>
class MapToSphere
{
private:
const double radius_;
const DomainVector center_;
public:
MapToSphere(const double radius=1.0, const DomainVector& center=DomainVector{0,0,0}):
radius_(radius),
center_(center)
{}
void evaluate( const DomainVector &x, RangeVector &y ) const
{
y = x;
y -= center_;
y /= x.two_norm()/radius_;
y += center_;
}
};
};
int main (int argc, char *argv[])
......@@ -128,7 +109,7 @@ int main (int argc, char *argv[])
Dune::MPIHelper::instance(argc, argv);
MakeRefinedSimplexGeometryTestSuite testSuite;
#if HAVE_DUNE_ALUGRID
return testSuite.check() ? 0 : 1;
return testSuite.test().exit();
#else
return 77;
#endif
......
......@@ -5,6 +5,7 @@
#include <limits>
#include <dune/common/exceptions.hh>
#include <dune/common/test/testsuite.hh>
#include <dune/fufem/test/common.hh>
#include <dune/fufem/geometry/refinedsimplexgeometry.hh>
......@@ -40,24 +41,37 @@ class RefinedSimplexGeometryTestSuite
}
bool check()
Dune::TestSuite test()
{
return true and
check_affine() and
check_type() and
check_corners() and
check_corner() and
check_center() and
check_global() and
check_local() and
check_integrationElement() and
check_volume() and
check_jacobianTransposed() and
check_jacobianInverseTransposed();
Dune::TestSuite test;
test.check(check_affine())
<< "Geometry is not affine";
test.check(check_type())
<< "type() is not implemented";
test.check(check_corners())
<< "Incorrect number of corners";
test.check(check_corner())
<< "corner() is inconsistent";
test.check(check_center())
<< "center() is not implemented";
test.check(check_global())
<< "global() is inconsistent";
test.check(check_local())
<< "local() is inconsistent";
test.check(check_integrationElement())
<< "integrationElement() is inconsistent";
test.check(check_volume())
<< "volume() is not implemented";
test.check(check_jacobianTransposed())
<< "jacobianTransposed() is inconsistent";
test.check(check_jacobianInverseTransposed())
<< "jacobianInverseTransposed() is inconsistent";
return test;
}
private:
bool check_affine()
{
return not testGeometry.affine();
......@@ -198,5 +212,5 @@ class RefinedSimplexGeometryTestSuite
int main (int argc, char *argv[])
{
RefinedSimplexGeometryTestSuite testSuite;
return not testSuite.check();
return testSuite.test().exit();
}
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