Skip to content

GitLab

Commit 16f0ab8f authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/isection-angle-quads' into 'master' 

Feature/isection angle quads

See merge request DennisGlaeser/frackit!33
parents e044e12c 1a02817f
Hide whitespace changes
Inline Side-by-side
frackit/distance/distancetoboundary.hh
View file @ 16f0ab8f
......@@ -24,6 +24,8 @@
#ifndef FRACKIT_DISTANCE_TO_BOUNDARY_HH
#define FRACKIT_DISTANCE_TO_BOUNDARY_HH
#include <cmath>
#include <limits>
#include <algorithm>
#include <Extrema_ExtAlgo.hxx>
......@@ -32,6 +34,7 @@
#include <frackit/precision/precision.hh>
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/cylindersurface.hh>
#include <frackit/geometry/name.hh>
......@@ -59,7 +62,8 @@ computeDistanceToBoundary(const Geo1& geo1,
std::string msg = "Distance to boundary computation not implemented for ";
msg += "\"" + geometryName(geo1) + "\"";
msg += " and ";
msg += "\"" + geometryName(geo1) + "\"";
msg += "\"" + geometryName(geo2) + "\"";
throw std::runtime_error(msg);
}
/*!
......@@ -87,6 +91,26 @@ computeDistanceToBoundary(const Geo& geo,
extAlgo);
}
/*!
* \brief Compute the distance of a point
* to the bounding of a quadrilateral.
* \param p The point
* \param quad The quadrilateral
*/
template<class ctype1, class ctype2>
PromotedType<ctype1, ctype2>
computeDistanceToBoundary(const Point<ctype1, 3>& p,
const Quadrilateral<ctype2, 3>& quad)
{
using std::min;
using ctype = PromotedType<ctype1, ctype2>;
ctype distance = std::numeric_limits<ctype>::max();
for (unsigned int i = 0; i < quad.numEdges(); ++i)
distance = min(distance, computeDistance(quad.edge(i), p));
return distance;
}
/*!
* \brief Compute the distance of a geometry
* to the bounding circles of a cylinder surface.
......
frackit/distance/pointonboundary.hh
View file @ 16f0ab8f
......@@ -50,7 +50,7 @@ namespace Frackit {
template<class ctype1, int wd, class Geo, class ctype2>
bool pointOnGeometryBoundary(const Point<ctype1, wd>& p, const Geo& geo, ctype2 eps)
{
static_assert(wd == DimensionalityTraits<Geo>::geomDim, "World dimension mismatch");
static_assert(wd == DimensionalityTraits<Geo>::worldDimension(), "World dimension mismatch");
std::string msg = "Point on boundary not implemented for ";
msg += "\"" + geometryName(geo) + "\" \n";
throw std::runtime_error(msg);
......@@ -68,6 +68,24 @@ bool pointOnGeometryBoundary(const Point<ctype1, wd>& p,
ctype2 eps)
{ return pointOnGeometry(p, disk.boundingEllipse(), eps); }
/*!
* \brief Evaluate if a point lies on the boundary of
* a quadrilateral in 3d space.
* \param p The point
* \param quad The quadrilateral
* \param eps Epsilon value to be used for the check
*/
template<class ctype1, class ctype2, class ctype3>
bool pointOnGeometryBoundary(const Point<ctype1, 3>& p,
const Quadrilateral<ctype2, 3>& quad,
ctype3 eps)
{
for (unsigned int i = 0; i < quad.numEdges(); ++i)
if (pointOnGeometry(p, quad.edge(i), eps))
return true;
return false;
}
/*!
* \brief Evaluate if a point lies on the boundary of a cylinder surface.
* \param p The point
......
frackit/entitynetwork/constraints.hh
View file @ 16f0ab8f
......@@ -94,6 +94,7 @@ public:
* \note This is only available if the engines are default constructible
*/
EntityNetworkConstraints()
: EntityNetworkConstraints(AngleComputationEngine())
{
static_assert(std::is_default_constructible<AE>::value,
"Angle computation engine not default constructible. "
......
frackit/intersection/intersectionangle.hh
View file @ 16f0ab8f
......@@ -49,7 +49,9 @@
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/line.hh>
#include <frackit/geometry/cylindersurface.hh>
#include <frackit/geometry/name.hh>
#include <frackit/common/extractdimension.hh>
#include <frackit/occ/gputilities.hh>
#include <frackit/occ/geomutilities.hh>
......@@ -68,6 +70,13 @@ namespace Frackit {
template<class ctype = double>
class IntersectionAngle
{
template<class G>
struct IsPlanarGeometry
{
static constexpr bool value = DimensionalityTraits<G>::geometryDimension() == 2
&& DimensionalityTraits<G>::worldDimension() == 3;
};
public:
/*!
......@@ -127,55 +136,67 @@ public:
}
/*!
* \brief Returns the angle in which two disks touch in a point
* \param disk1 The first plane
* \param disk2 The second plane
* \brief Returns the angle in which two planar 2-dimensional
* geometries embedded in 3d space touch in a point
* \param geo1 The first planar geometry
* \param geo2 The first planar geometry
* \param isPoint The touching point
*/
ctype operator() (const Disk<ctype>& disk1,
const Disk<ctype>& disk2,
template<class Geo1, class Geo2,
std::enable_if_t<IsPlanarGeometry<Geo1>::value
&& IsPlanarGeometry<Geo2>::value, int> = 0>
ctype operator() (const Geo1& geo1,
const Geo2& geo2,
const Point<ctype, 3>& isPoint)
{ return (*this)(disk1.supportingPlane(), disk2.supportingPlane()); }
{ return (*this)(geo1.supportingPlane(), geo2.supportingPlane()); }
/*!
* \brief Returns the angle in which two disks intersect in a segment
* \param disk1 The first plane
* \param disk2 The second plane
* \brief Returns the angle in which two planar 2-dimensional
* geometries embedded in 3d space touch intersect in a segment
* \param geo1 The first planar geometry
* \param geo2 The first planar geometry
* \param isSeg The intersection segment
*/
ctype operator() (const Disk<ctype>& disk1,
const Disk<ctype>& disk2,
template<class Geo1, class Geo2,
std::enable_if_t<IsPlanarGeometry<Geo1>::value
&& IsPlanarGeometry<Geo2>::value, int> = 0>
ctype operator() (const Geo1& geo1,
const Geo2& geo2,
const Segment<ctype, 3>& isSeg)
{ return (*this)(disk1.supportingPlane(), disk2.supportingPlane()); }
{ return (*this)(geo1.supportingPlane(), geo2.supportingPlane()); }
/*!
* \brief Returns the angle in which a disk and a
* cylinder surface touch in a point.
* \param disk The disk
* \brief Returns the angle in which a planar
* 2-dimensional geometry in 2d space
* and a cylinder surface touch in a point.
* \param geo The planar geometry
* \param cylSurface The cylinder surface
* \param isPoint The touching point
*/
ctype operator() (const Disk<ctype>& disk,
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const Geo& geo,
const CylinderSurface<ctype>& cylSurface,
const Point<ctype, 3>& isPoint)
{ return (*this)(disk.supportingPlane(), cylSurface.getTangentPlane(isPoint)); }
{ return (*this)(geo.supportingPlane(), cylSurface.getTangentPlane(isPoint)); }
/*!
* \brief Returns the angle in which a cylinder
* surface and a disk touch in a point.
* surface and planar 2-dimensional geometry
* in 2d space touch in a point.
* \param cylSurface The cylinder surface
* \param disk The disk
* \param geo The planar geometry
* \param isPoint The touching point
*/
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const CylinderSurface<ctype>& cylSurface,
const Disk<ctype>& disk,
const Geo& geo,
const Point<ctype, 3>& isPoint)
{ return(*this)(disk, cylSurface, isPoint); }
{ return(*this)(geo, cylSurface, isPoint); }
/*!
* \brief Returns the angle in which a disk
* \brief Returns the angle in which a planar 2-dimensional geometry
* and a cylinder surface intersect in a segment.
* \param disk The disk
* \param geo The planar geometry
* \param cylSurface The cylinder surface
* \param isSeg The intersection segment
* \note An intersection segment on the cylinder surface can only
......@@ -184,16 +205,17 @@ public:
* the intersection segment, and we compute the angle for an
* arbitrary point on the segment; here: the first corner.
*/
ctype operator() (const Disk<ctype>& disk,
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const Geo& geo,
const CylinderSurface<ctype>& cylSurface,
const Segment<ctype, 3>& isSeg)
{ return (*this)(disk.supportingPlane(), cylSurface.getTangentPlane(isSeg.source())); }
{ return (*this)(geo.supportingPlane(), cylSurface.getTangentPlane(isSeg.source())); }
/*!
* \brief Returns the angle in which a cylinder surface
* and a disk intersect in a segment.
* \param disk The disk
* \brief Returns the angle in which a cylinder surface and
* a planar 2-dimensional geometry intersect in a segment.
* \param cylSurface The cylinder surface
* \param geo The planar geometry
* \param isSeg The intersection segment
* \note An intersection segment on the cylinder surface can only
* occur if the intersection is parallel to the cylinder axis.
......@@ -201,65 +223,69 @@ public:
* the intersection segment, and we compute the angle for an
* arbitrary point on the segment; here: the first corner.
*/
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const CylinderSurface<ctype>& cylSurface,
const Disk<ctype>& disk,
const Geo& geo,
const Segment<ctype, 3>& isSeg)
{ return (*this)(disk, cylSurface, isSeg); }
{ return (*this)(geo, cylSurface, isSeg); }
/*!
* \brief Returns the angle in which a cylinder surface
* and a disk intersect in an ellipse arc.
* \param disk The disk
* \brief Returns the angle in which a planar 2-dimensional
* geometry and a cylinder surface intersect in an ellipse arc.
* \param geo The planar geometry
* \param cylSurface The cylinder surface
* \param isArc The intersection arc
*/
ctype operator() (const Disk<ctype>& disk,
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const Geo& geo,
const CylinderSurface<ctype>& cylSurface,
const EllipseArc<ctype, 3>& isArc)
{
// use the minimum angle between the disk plane and the
// use the minimum angle between the geometry plane and the
// tangent plane on the surface at four sample points
std::array<ctype, 5> params({0.0, 0.25, 0.5, 0.75, 1.0});
ctype resultAngle = std::numeric_limits<ctype>::max();
const auto& diskPlane = disk.supportingPlane();
const auto& geoPlane = geo.supportingPlane();
using std::min;
for (auto param : params)
{
const auto p = isArc.getPoint(param);
const auto tangentPlane = cylSurface.getTangentPlane(p);
resultAngle = min(resultAngle, (*this)(diskPlane, tangentPlane));
resultAngle = min(resultAngle, (*this)(geoPlane, tangentPlane));
}
return resultAngle;
}
/*!
* \brief Returns the angle in which a disk
* and a cylinder surface intersect in an ellipse arc.
* \brief Returns the angle in which a cyliner surface and
* a planar 2-dimensional geometry intersect in an ellipse arc.
* \param cylSurface The cylinder surface
* \param disk The disk
* \param geo The planar geometry
* \param isArc The intersection arc
*/
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const CylinderSurface<ctype>& cylSurface,
const Disk<ctype>& disk,
const Geo& geo,
const EllipseArc<ctype, 3>& isArc)
{ return (*this)(disk, cylSurface, isArc); }
{ return (*this)(geo, cylSurface, isArc); }
/*!
* \brief Returns the angle in which a disk
* \brief Returns the angle in which a planar 2-dimensional geometry
* and a cylinder surface intersect in an ellipse.
* \param cylSurface The cylinder surface
* \param disk The disk
* \param geo The planar geometry
* \param isEllipse The intersection ellipse
*/
ctype operator() (const Disk<ctype>& disk,
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const Geo& geo,
const CylinderSurface<ctype>& cylSurface,
const Ellipse<ctype, 3>& isEllipse)
{
// use the minimum angle between the disk plane and the
// use the minimum angle between the geometry plane and the
// tangent plane on the surface at eight sample points
const auto& diskPlane = disk.supportingPlane();
const auto& geoPlane = geo.supportingPlane();
ctype resultAngle = std::numeric_limits<ctype>::max();
std::array<ctype, 9> params({0.0, 0.125, 0.25, 0.375,
0.5, 0.625, 0.75, 0.875,
......@@ -270,7 +296,7 @@ public:
{
const auto p = isEllipse.getPoint(param);
const auto tangentPlane = cylSurface.getTangentPlane(p);
resultAngle = min(resultAngle, (*this)(diskPlane, tangentPlane));
resultAngle = min(resultAngle, (*this)(geoPlane, tangentPlane));
}
return resultAngle;
......@@ -278,15 +304,16 @@ public:
/*!
* \brief Returns the angle in which a cylinder surface
* and a disk intersect in an ellipse.
* \param disk The disk
* and a planar 2-dimensional geometry intersect in an ellipse.
* \param geo The planar geometry
* \param cylSurface The cylinder surface
* \param isEllipse The intersection ellipse
*/
template<class Geo, std::enable_if_t<IsPlanarGeometry<Geo>::value, int> = 0>
ctype operator() (const CylinderSurface<ctype>& cylSurface,
const Disk<ctype>& disk,
const Geo& geo,
const Ellipse<ctype, 3>& isEllipse)
{ return (*this)(disk, cylSurface, isEllipse); }
{ return (*this)(geo, cylSurface, isEllipse); }
/*!
* \brief Returns the angle in which a disk
......
test/entitynetwork/CMakeLists.txt
View file @ 16f0ab8f
......@@ -7,6 +7,11 @@ frackit_add_test(NAME test_constraints_disk
SOURCES test_constraints_disk.cc
LABELS entitynetwork constraints)
# test network constaints (quads)
frackit_add_test(NAME test_constraints_quads
SOURCES test_constraints_quads.cc
LABELS entitynetwork constraints)
# test network constaints (cylinder surface)
frackit_add_test(NAME test_constraints_cylindersurface
SOURCES test_constraints_cylindersurface.cc
......
test/entitynetwork/test_constraints_disk.cc
View file @ 16f0ab8f
......@@ -3,6 +3,7 @@
#include <stdexcept>
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/cylinder.hh>
#include <frackit/entitynetwork/constraints.hh>
......@@ -12,6 +13,7 @@ int main()
using ctype = double;
using Disk = Frackit::Disk<ctype>;
using Quad = Frackit::Quadrilateral<ctype, 3>;
using Point = typename Disk::Point;
using Direction = typename Disk::Direction;
using Vector = typename Direction::Vector;
......@@ -87,6 +89,27 @@ int main()
throw std::runtime_error("Did not detect intersection angle violation");
std::cout << "Test 9 passed" << std::endl;
// test both intersection angles also with quadrilaterals
{
// intersection angle just ok
Quad quad1(Point(-1.0, -1.0, -1.0 - 0.5e-3),
Point( 1.0, -1.0, -1.0 - 0.5e-3),
Point(-1.0, 1.0, 1.0 + 0.5e-3),
Point( 1.0, 1.0, 1.0 + 0.5e-3));
if (!constraints.evaluate(mainDisk, quad1))
throw std::runtime_error("False positive intersection angle violation");
std::cout << "Test 8 with quad passed" << std::endl;
// intersection angle too small
Quad quad2(Point(-1.0, -1.0, -1.0 + 0.5e-3),
Point( 1.0, -1.0, -1.0 + 0.5e-3),
Point(-1.0, 1.0, 1.0 - 0.5e-3),
Point( 1.0, 1.0, 1.0 - 0.5e-3));
if (constraints.evaluate(mainDisk, quad2))
throw std::runtime_error("Did not detect intersection angle violation");
std::cout << "Test 9 with quad passed" << std::endl;
}
// Test constraints w.r.t. cylinder
Frackit::Cylinder<ctype> cylinder(0.5, 1.0);
Disk disk10(Point(0.0, 0.0, 0.5), e1, e2, 2.0, 2.0);
......
test/entitynetwork/test_constraints_quads.cc 0 → 100644
View file @ 16f0ab8f
#include <cmath>
#include <string>
#include <stdexcept>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/geometry/cylinder.hh>
#include <frackit/entitynetwork/constraints.hh>
//! test the constraints for entity networks of quadrilaterals
int main()
{
using ctype = double;
using Quad = Frackit::Quadrilateral<ctype, 3>;
using Point = typename Quad::Point;
using Vector = typename Frackit::Vector<ctype, 3>;
// Basis Vectors
const Vector e1(1.0, 0.0, 0.0);
const Vector e2(0.0, 1.0, 0.0);
const Vector e3(0.0, 0.0, 1.0);
// Define constraints
Frackit::EntityNetworkConstraints<ctype> constraints;
constraints.setMinDistance(0.1);
constraints.setMinIntersectingAngle(M_PI/4.0);
constraints.setMinIntersectionMagnitude(0.05);
constraints.setMinIntersectionDistance(0.05);
// disk to test the others against
Quad mainQuad(Point(-1.0, -1.0, 0.0),
Point( 1.0, -1.0, 0.0),
Point(-1.0, 1.0, 0.0),
Point( 1.0, 1.0, 0.0));
// violates distance constraint
Quad quad1(Point(-1.0, -1.0, 0.1 - 1.0e-3),
Point( 1.0, -1.0, 0.1 - 1.0e-3),
Point(-1.0, 1.0, 0.1 - 1.0e-3),
Point( 1.0, 1.0, 0.1 - 1.0e-3));
if (constraints.evaluate(mainQuad, quad1))
throw std::runtime_error("Did not detect distance violation");
std::cout << "Test 1 passed" << std::endl;
// just doesn't violates distance constraint
Quad quad2(Point(-1.0, -1.0, 0.1 + 1.0e-3),
Point( 1.0, -1.0, 0.1 + 1.0e-3),
Point(-1.0, 1.0, 0.1 + 1.0e-3),
Point( 1.0, 1.0, 0.1 + 1.0e-3));
if (!constraints.evaluate(mainQuad, quad2))
throw std::runtime_error("Detected false positive distance violation");
std::cout << "Test 2 passed" << std::endl;
// violates distance constraint (orthogonal)
Quad quad3(Point(-1.0, 0.0, 0.1 - 1.0e-3),
Point( 1.0, 0.0, 0.1 - 1.0e-3),
Point(-1.0, 0.0, 1.1 - 1.0e-3),
Point( 1.0, 0.0, 1.1 - 1.0e-3));
if (constraints.evaluate(mainQuad, quad3))
throw std::runtime_error("Did not detect distance violation");
std::cout << "Test 3 passed" << std::endl;
// just doesn't violates distance constraint (orthogonal)
Quad quad4(Point(-1.0, 0.0, 0.1 + 1.0e-3),
Point( 1.0, 0.0, 0.1 + 1.0e-3),
Point(-1.0, 0.0, 1.1 + 1.0e-3),
Point( 1.0, 0.0, 1.1 + 1.0e-3));
if (!constraints.evaluate(mainQuad, quad4))
throw std::runtime_error("Detected false positive distance violation");
std::cout << "Test 4 passed" << std::endl;
// too small intersection
Quad quad5(Point(-0.05+1e-3, 0.0, 1.0),
Point( 0.0 , 0.0, -1.0),
Point( 0.0 , 0.0, 1.5),
Point( 0.05-1e-3, 0.0, 1.0));
if (constraints.evaluate(mainQuad, quad5))
throw std::runtime_error("Did not detect intersection magnitude violation");
std::cout << "Test 5 passed" << std::endl;
// intersection magnitude ok
Quad quad6(Point(-0.05-1e-3, 0.0, 1.0),
Point( 0.0 , 0.0, -1.0),
Point( 0.0 , 0.0, 1.5),
Point( 0.05+1e-3, 0.0, 1.0));
if (!constraints.evaluate(mainQuad, quad6))
throw std::runtime_error("False positive intersection magnitude violation");
std::cout << "Test 6 passed" << std::endl;
// intersection magnitude ok, but distance to boundary is too small
Quad quad7(Point(-2.0 + 0.0999, 0.0, 1.0),
Point( 0.0, 0.0, -1.0),
Point( 0.0, 0.0, 1.5),
Point( 2.0 - 0.0999, 0.0, 1.0));
if (constraints.evaluate(mainQuad, quad7))
throw std::runtime_error("Did not detect intersection distance violation");
std::cout << "Test 7 passed" << std::endl;
// intersection magnitude ok, distance to boundary is just ok
Quad quad8(Point(-2.0 + 0.1001, 0.0, 1.0),
Point( 0.0, 0.0, -1.0),
Point( 0.0, 0.0, 1.5),
Point( 2.0 - 0.1001, 0.0, 1.0));
if (!constraints.evaluate(mainQuad, quad8))
throw std::runtime_error("False positive intersection distance violation");
std::cout << "Test 8 passed" << std::endl;
// intersection angle just ok
Quad quad9(Point(-0.5, -0.5, -0.5 - 1e-3),
Point( 0.5, -0.5, -0.5 - 1e-3),
Point(-0.5, 0.5, 0.5 + 1e-3),
Point( 0.5, 0.5, 0.5 + 1e-3));
if (!constraints.evaluate(mainQuad, quad9))
throw std::runtime_error("False positive intersection angle violation");
std::cout << "Test 9 passed" << std::endl;
// intersection angle too small
Quad quad10(Point(-0.5, -0.5, -0.5 + 1e-3),
Point( 0.5, -0.5, -0.5 + 1e-3),
Point(-0.5, 0.5, 0.5 - 1e-3),
Point( 0.5, 0.5, 0.5 - 1e-3));
if (constraints.evaluate(mainQuad, quad10))
throw std::runtime_error("Did not detect intersection angle violation");
std::cout << "Test 10 passed" << std::endl;
// Test constraints w.r.t. cylinder
Frackit::Cylinder<ctype> cylinder(0.5, 1.0);
Quad quad11(Point(-0.6, -0.0, 1.0 - 0.05111),
Point( 0.0, -0.6, 1.0 - 0.05111),
Point( 0.6, 0.0, 1.0 - 0.05111),
Point( 0.0, 0.6, 1.0 - 0.05111));
if (!constraints.evaluate(cylinder.lateralFace(), quad11))
throw std::runtime_error("False positive intersection distance violation");
std::cout << "Test 11 passed" << std::endl;
// violates intersection distance constraint
Quad quad12(Point(-0.6, -0.0, 1.0 - 0.04999),
Point( 0.0, -0.6, 1.0 - 0.04999),
Point( 0.6, 0.0, 1.0 - 0.04999),
Point( 0.0, 0.6, 1.0 - 0.04999));
if (constraints.evaluate(cylinder.lateralFace(), quad12))
throw std::runtime_error("Did not detect intersection distance violation");
std::cout << "Test 12 passed" << std::endl;
// just doesn't violate distance constraint
Quad quad13(Point(-0.40001, -0.0, 01.5),
Point( 0.0, -0.40001, 01.5),
Point( 0.40001, 0.0, 01.5),
Point( 0.0, 0.40001, 01.5));
if (!constraints.evaluate(cylinder.lateralFace(), quad13))
throw std::runtime_error("False positive distance violation");
std::cout << "Test 13 passed" << std::endl;
// just violates distance constraint
Quad quad14(Point(-0.39999, -0.0, 0.5),
Point( 0.0, -0.39999, 0.5),
Point( 0.39999, 0.0, 0.5),
Point( 0.0, 0.39999, 0.5));
if (constraints.evaluate(cylinder.lateralFace(), quad14))
throw std::runtime_error("Did not detect distance violation");
std::cout << "Test 14 passed" << std::endl;
std::cout << "All tests passed" << std::endl;
return 0;
}
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