Commit 875a1bf2 authored by Dennis Gläser's avatar Dennis Gläser
Browse files

Merge branch 'feature/quad-face-intersections' into 'master'

Feature/quad face intersections

See merge request DennisGlaeser/frackit!34
parents 16f0ab8f aec074f0
......@@ -5,6 +5,8 @@ algo_cylsurface_quadrilateral.hh
algo_disk_disk.hh
algo_disk_line.hh
algo_face_disk.hh
algo_face_planargeom.hh
algo_face_quadrilateral.hh
algo_find_touching_points.hh
algo_planargeom_line.hh
algo_planargeom_planargeom.hh
......
......@@ -24,20 +24,8 @@
#ifndef FRACKIT_FACE_DISK_INTERSECTION_HH
#define FRACKIT_FACE_DISK_INTERSECTION_HH
#include <algorithm>
#include <stdexcept>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <frackit/geometry/point.hh>
#include <frackit/geometry/disk.hh>
#include <frackit/occ/breputilities.hh>
#include "algo_find_touching_points.hh"
#include "algo_face_planargeom.hh"
#include "intersectiontraits.hh"
#include "emptyintersection.hh"
namespace Frackit {
namespace IntersectionAlgorithms {
......@@ -54,102 +42,7 @@ Intersection< TopoDS_Face, Disk<ctype> >
intersect_face_disk(const TopoDS_Face& face,
const Disk<ctype>& disk,
ctype eps)
{
using ResultType = Intersection< TopoDS_Face, Disk<ctype> >;
std::vector<Point<ctype, 3>> pointIntersections;
std::vector<TopoDS_Edge> edgeIntersections;
std::vector<TopoDS_Face> faceIntersections;
// find possible face intersections
const auto diskShape = OCCUtilities::getShape(disk);
const auto diskIS = OCCUtilities::intersect(diskShape, face, eps);
const auto diskISFaces = OCCUtilities::getFaces(diskIS);
for (auto& f : diskISFaces)
faceIntersections.emplace_back(std::move(f));
// find possible edge intersections
const auto diskCut = OCCUtilities::cut(diskShape, face, eps);
const auto diskCutFaces = OCCUtilities::getFaces(diskCut);
for (unsigned int i = 0; i < diskCutFaces.size(); ++i)
{
// edges of the fragments coinciding with the face
const auto cutFaceWires = OCCUtilities::getWires(diskCutFaces[i]);
for (const auto& wire : cutFaceWires)
{
const auto wireIS = OCCUtilities::intersect(wire, face, eps);
const auto edgesOnFace = OCCUtilities::getEdges(wireIS);
for (auto& edge : edgesOnFace)
{
// only use this edge if it is not on a previously found face
bool isContained = false;
for (const auto& faceIS : faceIntersections)
{
const auto edgeFaceIS = OCCUtilities::intersect(edge, faceIS, eps);
const auto edgeFaceISEdges = OCCUtilities::getEdges(edgeFaceIS);
if (!edgeFaceISEdges.empty()) { isContained = true; break; }
}
if (isContained)
continue;
if (std::none_of(edgeIntersections.begin(),
edgeIntersections.end(),
[&] (const auto& e) { return edge.IsSame(e); }))
edgeIntersections.emplace_back(edge);
}
}
}
// Find possible touching points. We might find points
// again that resulted from intersections in faces or edges.
// Thus, we have to check for duplicates here and avoid all
// points that have been detected so far. The newly detected
// ones are true touching points. Note that this would not be
// necessary if we restricted the face to be planar. But, on
// nonplanar and periodic faces you can have all kinds of
// intersections. To this end, collect all points found so far.
std::vector<Point<ctype, 3>> curIsPoints;
for (const auto& faceIS : faceIntersections)
for (const auto& v : OCCUtilities::getVertices(faceIS))
{
const auto p = OCCUtilities::point(v);
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
curIsPoints.emplace_back(p);
}
for (const auto& edgeIS : edgeIntersections)
for (const auto& v : OCCUtilities::getVertices(edgeIS))
{
const auto p = OCCUtilities::point(v);
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
curIsPoints.emplace_back(p);
}
const auto touchPoints = find_touching_points(diskShape, face, eps);
for (auto&& p : touchPoints)
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); })
&& std::none_of(pointIntersections.begin(),
pointIntersections.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
pointIntersections.emplace_back(std::move(p));
ResultType result;
for (auto& p : pointIntersections) result.emplace_back(std::move(p));
for (auto& e : edgeIntersections) result.emplace_back(std::move(e));
for (auto& f : faceIntersections) result.emplace_back(std::move(f));
if (result.empty())
return ResultType({EmptyIntersection<3>()});
return result;
}
{ return intersect_face_planarGeometry(face, disk, eps); }
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between a face shape and a planar
* two-dimensional geometry in 3 space.
*/
#ifndef FRACKIT_FACE_PLANAR_GEOM_INTERSECTION_HH
#define FRACKIT_FACE_PLANAR_GEOM_INTERSECTION_HH
#include <algorithm>
#include <stdexcept>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <frackit/geometry/point.hh>
#include <frackit/occ/breputilities.hh>
#include <frackit/common/extractctype.hh>
#include "algo_find_touching_points.hh"
#include "intersectiontraits.hh"
#include "emptyintersection.hh"
namespace Frackit {
namespace IntersectionAlgorithms {
/*!
* \brief Intersect a face shape and a planar 2d geometry in 3d space.
* The result can be:
* - points
* - edge shapes
* - face shapes
* Multiple of the above are possible since the face shape might be curved.
* \param face The (possible curved) face shape
* \param faceGeom The planar geometry (planar bounded 2d surface in 3d space)
* \param eps Tolerance to be used for boolean operations
*/
template<class PlanarGeometry>
Intersection< TopoDS_Face, PlanarGeometry >
intersect_face_planarGeometry(const TopoDS_Face& face,
const PlanarGeometry& faceGeom,
typename CoordinateTypeTraits<PlanarGeometry>::type eps)
{
using ResultType = Intersection< TopoDS_Face, PlanarGeometry >;
using ctype = typename CoordinateTypeTraits<PlanarGeometry>::type;
std::vector<Point<ctype, 3>> pointIntersections;
std::vector<TopoDS_Edge> edgeIntersections;
std::vector<TopoDS_Face> faceIntersections;
// find possible face intersections
const auto faceGeomShape = OCCUtilities::getShape(faceGeom);
const auto faceGeomIS = OCCUtilities::intersect(faceGeomShape, face, eps);
const auto faceGeomISFaces = OCCUtilities::getFaces(faceGeomIS);
for (auto& f : faceGeomISFaces)
faceIntersections.emplace_back(std::move(f));
// find possible edge intersections
const auto faceGeomCut = OCCUtilities::cut(faceGeomShape, face, eps);
const auto faceGeomCutFaces = OCCUtilities::getFaces(faceGeomCut);
for (unsigned int i = 0; i < faceGeomCutFaces.size(); ++i)
{
// edges of the fragments coinciding with the face
const auto cutFaceWires = OCCUtilities::getWires(faceGeomCutFaces[i]);
for (const auto& wire : cutFaceWires)
{
const auto wireIS = OCCUtilities::intersect(wire, face, eps);
const auto edgesOnFace = OCCUtilities::getEdges(wireIS);
for (auto& edge : edgesOnFace)
{
// only use this edge if it is not on a previously found face
bool isContained = false;
for (const auto& faceIS : faceIntersections)
{
const auto edgeFaceIS = OCCUtilities::intersect(edge, faceIS, eps);
const auto edgeFaceISEdges = OCCUtilities::getEdges(edgeFaceIS);
if (!edgeFaceISEdges.empty()) { isContained = true; break; }
}
if (isContained)
continue;
if (std::none_of(edgeIntersections.begin(),
edgeIntersections.end(),
[&] (const auto& e) { return edge.IsSame(e); }))
edgeIntersections.emplace_back(edge);
}
}
}
// Find possible touching points. We might find points
// again that resulted from intersections in faces or edges.
// Thus, we have to check for duplicates here and avoid all
// points that have been detected so far. The newly detected
// ones are true touching points. Note that this would not be
// necessary if we restricted the face to be planar. But, on
// nonplanar and periodic faces you can have all kinds of
// intersections. To this end, collect all points found so far.
std::vector<Point<ctype, 3>> curIsPoints;
for (const auto& faceIS : faceIntersections)
for (const auto& v : OCCUtilities::getVertices(faceIS))
{
const auto p = OCCUtilities::point(v);
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
curIsPoints.emplace_back(p);
}
for (const auto& edgeIS : edgeIntersections)
for (const auto& v : OCCUtilities::getVertices(edgeIS))
{
const auto p = OCCUtilities::point(v);
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
curIsPoints.emplace_back(p);
}
const auto touchPoints = find_touching_points(faceGeomShape, face, eps);
for (auto&& p : touchPoints)
if (std::none_of(curIsPoints.begin(),
curIsPoints.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); })
&& std::none_of(pointIntersections.begin(),
pointIntersections.end(),
[&] (const auto& isP) { return isP.isEqual(p, eps); }))
pointIntersections.emplace_back(std::move(p));
ResultType result;
for (auto& p : pointIntersections) result.emplace_back(std::move(p));
for (auto& e : edgeIntersections) result.emplace_back(std::move(e));
for (auto& f : faceIntersections) result.emplace_back(std::move(f));
if (result.empty())
return ResultType({EmptyIntersection<3>()});
return result;
}
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_FACE_PLANAR_GEOM_INTERSECTION_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between a face shape and a quadrilateral in 3d space.
*/
#ifndef FRACKIT_FACE_QUADRILATERAL_INTERSECTION_HH
#define FRACKIT_FACE_QUADRILATERAL_INTERSECTION_HH
#include "algo_face_planargeom.hh"
#include "intersectiontraits.hh"
namespace Frackit {
namespace IntersectionAlgorithms {
//! Intersect a face shape and a quadrilateral in 3d space.
//! The result can be composed of:
//! - points
//! - edges
//! - faces
//! Multiple of the above are possible
//! since the face shape might be curved.
template<class ctype>
Intersection< TopoDS_Face, Quadrilateral<ctype, 3> >
intersect_face_quadrilateral(const TopoDS_Face& face,
const Quadrilateral<ctype, 3>& quad,
ctype eps)
{ return intersect_face_planarGeometry(face, quad, eps); }
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_FACE_QUADRILATERAL_INTERSECTION_HH
......@@ -50,6 +50,7 @@
#include "algo_cylsurface_quadrilateral.hh"
#include "algo_shell_disk.hh"
#include "algo_face_disk.hh"
#include "algo_face_quadrilateral.hh"
namespace Frackit {
......@@ -308,6 +309,28 @@ Intersection< TopoDS_Face, Disk<ctype> >
intersect(const TopoDS_Face& face, const Disk<ctype>& disk, ctype eps)
{ return intersect(disk, face, eps); }
/*!
* \brief Intersect a quadrilateral and a face shape.
* \param quad The quadrilateral
* \param face The face shape
* \param eps Tolerance to be used for floating point comparisons
*/
template<class ctype>
Intersection< Quadrilateral<ctype, 3>, TopoDS_Face >
intersect(const Quadrilateral<ctype, 3>& quad, const TopoDS_Face& face, ctype eps)
{ return IntersectionAlgorithms::intersect_face_quadrilateral(face, quad, eps); }
/*!
* \brief Intersect a face shape and a quadrilateral.
* \param face The face shape
* \param quad The quadrilateral
* \param eps Tolerance to be used for floating point comparisons
*/
template<class ctype>
Intersection< TopoDS_Face, Quadrilateral<ctype, 3> >
intersect(const TopoDS_Face& face, const Quadrilateral<ctype, 3>& quad, ctype eps)
{ return intersect(quad, face, eps); }
} // end namespace Frackit
#endif // FRACKIT_INTERSECT_HH
......@@ -227,16 +227,18 @@ struct IntersectionTraits< Disk<ctype>, TopoDS_Shell >
: public IntersectionTraits< TopoDS_Shell, Disk<ctype> >
{};
namespace IntersectionDetail {
template<class ctype>
using FaceIntersectionWith2d = std::vector< std::variant< Point<ctype, 3>,
TopoDS_Edge,
TopoDS_Face,
EmptyIntersection<3> > >;
} // end namespace IntersectionDetail
//! Result type of the intersection of a face and a disk
template<class ctype>
struct IntersectionTraits< TopoDS_Face, Disk<ctype> >
{
using BaseType = std::variant< Point<ctype, 3>,
TopoDS_Edge,
TopoDS_Face,
EmptyIntersection<3> >;
using type = std::vector<BaseType>;
};
{ using type = IntersectionDetail::FaceIntersectionWith2d<ctype>; };
//! Result type of the intersection of a disk and a face
template<class ctype>
......@@ -244,6 +246,17 @@ struct IntersectionTraits< Disk<ctype>, TopoDS_Face >
: public IntersectionTraits< TopoDS_Face, Disk<ctype> >
{};
//! Result type of the intersection of a quadrilateral and a face
template<class ctype>
struct IntersectionTraits< Quadrilateral<ctype, 3>, TopoDS_Face >
{ using type = IntersectionDetail::FaceIntersectionWith2d<ctype>; };
//! Result type of the intersection of a face and a quadrilateral
template<class ctype>
struct IntersectionTraits< TopoDS_Face, Quadrilateral<ctype, 3> >
: public IntersectionTraits< Quadrilateral<ctype, 3>, TopoDS_Face >
{};
} // end namespace Frackit
#endif // FRACKIT_INTERSECTION_TRAITS_HH
......@@ -52,6 +52,11 @@ frackit_add_test(NAME test_disk_face_intersection
SOURCES test_disk_face_intersection.cc
LABELS intersection)
# face - quadrilateral
frackit_add_test(NAME test_quad_face_intersection
SOURCES test_quad_face_intersection.cc
LABELS intersection)
# shell - disk
frackit_add_test(NAME test_disk_shell_intersection
SOURCES test_disk_shell_intersection.cc
......
#include <TopExp.hxx>
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/quadrilateral.hh>
#include <frackit/magnitude/length.hh>
#include <frackit/magnitude/area.hh>
#include <frackit/precision/precision.hh>
#include <frackit/occ/breputilities.hh>
#include <frackit/intersection/intersect.hh>
enum IntersectionType { point, edge, face, empty };
template<class G>
void checkResultGeometry(const G& geometry, IntersectionType expected)
{
throw std::runtime_error("Unexpected intersection geometry");
}
template<int wd>
void checkResultGeometry(const Frackit::EmptyIntersection<wd>& empty, IntersectionType expected)
{
std::cout << "Found empty intersection" << std::endl;
if (expected != IntersectionType::empty)
throw std::runtime_error("Got an unexpected empty intersection");
}
template<class CT, int wd>
void checkResultGeometry(const Frackit::Point<CT, wd>& p, IntersectionType expected)
{
std::cout << "Found intersection point at " << p << std::endl;
if (expected != IntersectionType::point)
throw std::runtime_error("Got an unexpected point intersection");
}
void checkResultGeometry(const TopoDS_Edge& edge, IntersectionType expected)
{
std::cout << "Found intersection edge with corners "
<< Frackit::OCCUtilities::point(TopExp::FirstVertex(edge)) << " - "
<< Frackit::OCCUtilities::point(TopExp::LastVertex(edge))
<< " and length " << Frackit::computeLength(edge) << std::endl;
if (expected != IntersectionType::edge)
throw std::runtime_error("Got an unexpected segment intersection");
}
void checkResultGeometry(const TopoDS_Face& face, IntersectionType expected)
{
std::cout << "Found intersection face with area "
<< Frackit::computeArea(face) << std::endl;
if (expected != IntersectionType::face)
throw std::runtime_error("Got an unexpected face intersection");
}
//! test disk-face intersections
int main()
{
using ctype = double;
using Quad = Frackit::Quadrilateral<ctype, 3>;
using Point = typename Quad::Point;
std::vector<ctype> scales({1.0e-3, 1.0, 1.0e3});
for (auto f : scales)
{
std::cout << "Checking scale factor " << f << std::endl;
const auto eps = Frackit::Precision<ctype>::confusion()*f;
// Make quad
Quad quadGeom(Point(0.0, 0.0, 0.0), Point(f, 0.0, 0.0),
Point(0.0, f, 0.0), Point(f, f, 0.0));
const auto quadFace = Frackit::OCCUtilities::getShape(quadGeom);
// create a quad that intersects in an edge of length f
Quad quad1(Point(-f, 0.0, -f), Point(f, 0.0, -f),
Point(-f, 0.0, f), Point(f, 0.0, f));
auto result = intersect(quad1, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::edge); }, result[0]);
using std::abs;
if ( abs(Frackit::computeLength(std::get<TopoDS_Edge>(result[0])) - f) > eps )
throw std::runtime_error("Unexpected intersection edge length");
std::cout << "Test 1 passed" << std::endl;
// create a quad that intersects in an edge of length 0.5*f
Quad quad2(Point(-0.5*f, 0.0, -f), Point(0.5*f, 0.0, -f),
Point(-0.5*f, 0.0, f), Point(0.5*f, 0.0, f));
result = intersect(quad2, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
for (const auto& isection : result)
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::edge); }, isection);
using std::abs;
if ( abs(Frackit::computeLength(std::get<TopoDS_Edge>(result[0])) - 0.5*f) > eps )
throw std::runtime_error("Unexpected intersection edge 1 length");
std::cout << "Test 2 passed" << std::endl;
// 1 touching point
Quad quad3(Point(-0.5*f, 0.0, 0.5*f), Point(0.0, 0.0, 0.0),
Point( 0.0, 0.0, f), Point(0.5*f, 0.0, 0.5*f));
result = intersect(quad3, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
for (const auto& isection : result)
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::point); }, isection);
std::cout << "Test 3 passed" << std::endl;
// 1 touching point (in-plane)
Quad quad4(Point(-0.5*f, 0.0, 0.0), Point(0.0, 0.5*f, 0.0),
Point(-0.5*f, 0.5*f, 0.0), Point(-0.25*f, 1.0*f, 0.0));
result = intersect(quad4, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
for (const auto& isection : result)
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::point); }, isection);
std::cout << "Test 4 passed" << std::endl;
// no touching point
Quad quad5(Point(-0.5*f, 0.0, 0.5*f), Point(0.0, 0.0, 1e-3*f),
Point( 0.0, 0.0, f), Point(0.5*f, 0.0, 0.5*f));
result = intersect(quad5, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::empty); }, result[0]);
std::cout << "Test 5 passed" << std::endl;
// intersection face that has an area of f^2
Quad quad6(Point(-2.0*f, -2.0*f, 0.0), Point(2.0*f, -2.0*f, 0.0),
Point(-2.0*f, 2.0*f, 0.0), Point(2.0*f, 2.0*f, 0.0));
result = intersect(quad6, quadFace);
if (result.size() != 1)
throw std::runtime_error("Unexpected intersection size");
std::visit([&] (auto&& is) { checkResultGeometry(is, IntersectionType::face); }, result[0]);
if ( abs(Frackit::computeArea(std::get<TopoDS_Face>(result[0])) - f*f) > eps*eps )
throw std::runtime_error("Unexpected intersection face area");
std::cout << "Test 5 passed" << std::endl;
std::cout << "All tests passed" << std::endl;
}
return 0;
}