CSG for 2D

This commit is contained in:
Matthias Hochsteger 2020-08-19 16:46:32 +02:00
parent 9e105c48ea
commit 12b2e073ac
9 changed files with 2172 additions and 272 deletions

View File

@ -1,5 +1,5 @@
add_definitions(-DNGLIB_EXPORTS)
add_library(geom2d ${NG_LIB_TYPE} genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
add_library(geom2d ${NG_LIB_TYPE} csg2d.cpp genmesh2d.cpp geometry2d.cpp python_geom2d.cpp )
if(APPLE)
set_target_properties( geom2d PROPERTIES SUFFIX ".so")
endif(APPLE)
@ -20,6 +20,6 @@ endif(USE_GUI)
install(FILES
geometry2d.hpp spline2d.hpp
vsgeom2d.hpp
vsgeom2d.hpp csg2d.hpp
DESTINATION ${NG_INSTALL_DIR_INCLUDE}/geom2d COMPONENT netgen_devel
)

1549
libsrc/geom2d/csg2d.cpp Normal file

File diff suppressed because it is too large Load Diff

583
libsrc/geom2d/csg2d.hpp Normal file
View File

@ -0,0 +1,583 @@
#ifndef NETGEN_CSG2D_HPP_INCLUDED
#define NETGEN_CSG2D_HPP_INCLUDED
#include "geometry2d.hpp"
namespace netgen
{
using namespace ngcore;
using netgen::Point;
using netgen::Vec;
inline double Area(const Point<2>& P, const Point<2>& Q, const Point<2>& R)
{
return (Q[0]-P[0]) * (R[1]-P[1]) - (Q[1]-P[1]) * (R[0]-P[0]);
}
enum IntersectionType
{ // types of intersection (detected in the first phase)
NO_INTERSECTION = 0,
X_INTERSECTION,
T_INTERSECTION_Q,
T_INTERSECTION_P,
V_INTERSECTION,
X_OVERLAP,
T_OVERLAP_Q,
T_OVERLAP_P,
V_OVERLAP
};
enum IntersectionLabel
{ // for the classification of intersection vertices in the second phase
NONE,
CROSSING,
BOUNCING,
LEFT_ON,
RIGHT_ON,
ON_ON,
ON_LEFT,
ON_RIGHT,
DELAYED_CROSSING,
DELAYED_BOUNCING
};
enum EntryExitLabel
{ // for marking intersection vertices as "entry" or "exit"
EXIT,
ENTRY,
NEITHER
};
enum IteratorType
{
SOURCE,
INTERSECTION,
CROSSING_INTERSECTION,
ALL
};
using Spline = SplineSeg3<2>;
struct Vertex : Point<2>
{
Vertex (Point<2> p) : Point<2>(p) {}
Vertex * prev = nullptr;
Vertex * next = nullptr;
unique_ptr<Vertex> pnext = nullptr;
Vertex * neighbour = nullptr; // same vertex in other polygon (at intersections)
double lam = -1.0;
bool is_intersection = false;
bool is_source = false;
IntersectionLabel label = NONE; // type of intersection vertex
EntryExitLabel enex = NEITHER; // entry/exit "flag"
string bc = "";
// In case the edge this - next is curved, store the spline information here
optional<Spline> spline = nullopt;
Vertex * Insert(Point<2> p, double lam = -1.0);
void Link( Vertex * v )
{
neighbour = v;
v->neighbour = this;
is_intersection = true;
v->is_intersection = true;
}
};
struct VertexIterator
{
struct iterator
{
iterator(Vertex* root, IteratorType IterType) :
root(root), V(NULL), iterType(IterType)
{
if (root == NULL)
return;
if (nextVertex() == NULL) // no (source/intersection) vertex found
root = V = NULL; // -> mark iterator as "end"
}
const iterator& operator++()
{
nextVertex();
return *this;
}
Vertex* operator*()
{
return V;
}
bool operator!=(const iterator& other) const
{
return (root != other.root) || (V != other.V);
}
private:
Vertex* root;
Vertex* V;
IteratorType iterType;
//
// find the next vertex
// if iterType is ALL, then it is just the next vertex
// if iterType is SOURCE, then it is the next source vertex
// if iterType is INTERSECTION, then it is the next intersection vertex
// if iterType is CROSSING_INTERSECTION, then it is the next intersection vertex with CROSSING label
//
Vertex* nextVertex()
{
bool nextFound = false;
if (V == NULL)
{ // find first (source/intersection) vertex
V = root;
switch(iterType)
{
case ALL:
nextFound = true;
break;
case SOURCE:
if (V->is_source)
nextFound = true;
break;
case INTERSECTION:
if (V->is_intersection)
nextFound = true;
break;
case CROSSING_INTERSECTION:
if (V->is_intersection && (V->label == CROSSING))
nextFound = true;
break;
}
}
while (!nextFound)
{ // find next (source/intersection) vertex
switch(iterType)
{
case ALL:
V = V->next;
break;
case SOURCE:
do {
V = V->next;
} while (!V->is_source && V != root);
break;
case INTERSECTION:
do {
V = V->next;
} while (!V->is_intersection && V != root);
break;
case CROSSING_INTERSECTION:
do {
V = V->next;
} while ( ( !V->is_intersection || (V->label != CROSSING) ) && V != root);
break;
}
if (V == root)
{ // back at the root vertex?
root = V = NULL; // -> mark iterator as "end"
return(V);
}
switch(iterType)
{
case ALL:
nextFound = true;
break;
case SOURCE:
if (V->is_source)
nextFound = true;
break;
case INTERSECTION:
if (V->is_intersection)
nextFound = true;
break;
case CROSSING_INTERSECTION:
if (V->is_intersection && (V->label == CROSSING))
nextFound = true;
break;
}
}
return(V);
}
};
public:
VertexIterator() : root(NULL) {};
iterator begin() { return iterator(root, iterType); }
iterator end() { return iterator(NULL, iterType); }
Vertex* root;
IteratorType iterType;
};
struct Edge
{
Vertex * v0 = nullptr;
Vertex * v1 = nullptr;
Edge (Vertex* v, Vertex* w) : v0(v), v1(w) { };
};
struct EdgeIterator
{
struct iterator
{
iterator(Vertex* root, IteratorType IterType) :
root(root), one(NULL), two(NULL), iterType(IterType)
{
if (root == NULL)
return;
if (nextEdge() == NULL) // no source edge found
root = one = two = NULL; // -> mark iterator as "end"
}
const iterator& operator++() { nextEdge(); return *this; }
Edge operator*()
{
return Edge(one,two);
}
bool operator!=(const iterator& other) const
{
return (root != other.root) || (one != other.one) || (two != other.two);
}
private:
Vertex* root;
Vertex* one;
Vertex* two;
IteratorType iterType;
//
// find the next vertex, starting at curr
// if iterType is ALL, then it is just the next vertex
// if iterType is SOURCE, then it is the next source vertex
//
Vertex* nextVertex(Vertex* curr)
{
if (curr == NULL)
return(NULL);
switch(iterType)
{
case ALL:
curr = curr->next;
break;
case SOURCE:
do {
curr = curr->next;
} while (!curr->is_source);
break;
default:
;
}
return(curr);
}
//
// find the next edge
//
Vertex* nextEdge()
{
if (root == NULL) // empty polygon?
return (NULL);
if (one == NULL)
{ // find one (source) vertex
one = root; // note: root is always a (source) vertex
two = nextVertex(one);
if (two == one) // just one (source) vertex
return(NULL); // -> no (source) edges
return(one);
}
if (two == root)
{ // back at the root vertex?
root = one = two = NULL; // -> mark iterator as "end"
return(NULL);
}
one = two;
two = nextVertex(one);
return (one);
}
};
public:
EdgeIterator() : root(NULL) {};
iterator begin() { return iterator(root, iterType); }
iterator end() { return iterator(NULL, iterType); }
Vertex* root;
IteratorType iterType;
};
inline int CalcSide( const Point<2> & p0, const Point<2> & p1, const Point<2> & r )
{
if ( (p0[1] < r[1]) != (p1[1] < r[1]) )
{
if (p0[0] >= r[0])
{
if (p1[0] > r[0])
return 2 * (p1[1] > p0[1]) - 1;
else
if ( (Area(p0,p1,r) > 0) == (p1[1] > p0[1]) )
return 2 * (p1[1] > p0[1]) - 1;
}
else
{
if (p1[0] > r[0])
if ( (Area(p0,p1,r) > 0) == (p1[1] > p0[1]) )
return 2 * (p1[1] > p0[1]) - 1;
}
}
return 0;
}
struct Polygon2d
{
unique_ptr<Vertex> first = nullptr;
Polygon2d() = default;
Polygon2d(const Polygon2d & p)
: first(nullptr)
{
for(auto v : p.Vertices(ALL))
AppendVertex(*v);
}
Polygon2d & operator=(const Polygon2d & p)
{
first = nullptr;
if(p.first)
for(const auto v : p.Vertices(ALL))
AppendVertex(*v);
return *this;
}
Vertex & AppendVertex(const Vertex & v)
{
auto & vnew = Append( static_cast<Point<2>>(v), true );
vnew.bc = v.bc;
if(v.spline)
vnew.spline = *v.spline;
return vnew;
}
Vertex & Append(Point<2> p, bool source = false)
{
Vertex * vnew;
if(first==nullptr)
{
first = make_unique<Vertex>(p);
first->next = first.get();
first->prev = first.get();
vnew = first.get();
}
else
{
vnew = first->prev->Insert(p);
}
vnew->is_source = source;
// cout << "size after " << Size() << endl;
return *vnew;
}
void Remove (Vertex* v)
{
v->prev->next = v->next;
v->next->prev = v->prev;
if(first.get() == v)
first = std::move(v->pnext);
else
v->prev->pnext = std::move(v->pnext);
}
bool IsInside( Point<2> r ) const
{
int w = 0;
for(auto e : Edges(ALL))
w += CalcSide(*e.v0, *e.v1, r);
return ( (w % 2) != 0 );
}
EdgeIterator Edges(IteratorType iterType) const
{
EdgeIterator it;
it.iterType = iterType;
it.root = first.get();
return it;
}
VertexIterator Vertices(IteratorType iterType, Vertex* first_ = nullptr) const
{
VertexIterator it;
it.iterType = iterType;
it.root = (first_ == nullptr) ? first.get() : first_;
return it;
}
//
// check, if all vertices have the ON_ON label
//
bool allOnOn()
{
for (Vertex* v : Vertices(ALL))
if (v->label != ON_ON)
return(false);
return(true);
}
//
// check, if the polygon does not contain any crossing intersection vertex
// or crossing intersection chain or (if we want to compute the union instead
// of the intersection) a bouncing vertex or a bouncing intersection chain
//
bool noCrossingVertex(bool union_case = false)
{
for (Vertex* v : Vertices(ALL))
if (v->is_intersection)
{
if ( (v->label == CROSSING) || (v->label == DELAYED_CROSSING) )
return(false);
if (union_case && ( (v->label == BOUNCING) || (v->label == DELAYED_BOUNCING) ) )
return(false);
}
return(true);
}
//
// return a non-intersection point
//
Point<2> getNonIntersectionPoint()
{
for (Vertex* v : Vertices(ALL))
if (!v->is_intersection)
return *v;
// no non-intersection vertex found -> find suitable edge midpoint
for (Vertex* v : Vertices(ALL))
// make sure that edge from V to V->next is not collinear with other polygon
if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
// return edge midpoint
return Center(*v, *v->next);
throw Exception("no point found");
}
//
// return and insert a non-intersection vertex
//
Vertex* getNonIntersectionVertex()
{
for (Vertex* v : Vertices(ALL))
if (!v->is_intersection)
return(v);
// no non-intersection vertex found -> generate and return temporary vertex
for (Vertex* v : Vertices(ALL))
// make sure that edge from V to V->next is not collinear with other polygon
if ( (v->next->neighbour != v->neighbour->prev) && (v->next->neighbour != v->neighbour->next) )
{
// add edge midpoint as temporary vertex
auto p = Center(*v, *v->next);
return v->Insert(p);
}
return(NULL);
}
void SetBC(string bc)
{
for(auto v : Vertices(ALL))
v->bc = bc;
}
size_t Size() const
{
if(first==nullptr) return 0;
size_t cnt = 0;
for(auto v : Vertices(ALL))
cnt++;
return cnt;
}
};
struct Solid2d
{
Array<Polygon2d> polys;
string name = "";
Solid2d() = default;
Solid2d(string name_) : name(name_) {}
Solid2d operator+(Solid2d & other);
Solid2d operator*(Solid2d & other);
Solid2d operator-(Solid2d other);
void Append( const Polygon2d & poly )
{
polys.Append(poly);
}
bool IsInside( Point<2> r ) const;
bool IsLeftInside( const Vertex & p0 );
bool IsRightInside( const Vertex & p0 );
void SetBC(string bc)
{
for(auto & p : polys)
for(auto v : p.Vertices(ALL))
v->bc = bc;
}
};
class CSG2d
{
public:
Array<Solid2d> solids;
void Add ( Solid2d s )
{
solids.Append(s);
}
shared_ptr<netgen::SplineGeometry2d> GenerateSplineGeometry();
};
Solid2d Circle(double x, double y, double r, string name="", string bc="");
Solid2d Rectangle(double x0, double x1, double y0, double y1, string name="", string bc="");
Solid2d AddIntersectionPoints ( Solid2d s1, Solid2d s2 );
Solid2d ClipSolids ( Solid2d s1, Solid2d s2, bool intersect=true );
}
#endif // NETGEN_CSG2D_HPP_INCLUDED

View File

@ -9,6 +9,7 @@
#include <myadt.hpp>
#include <gprim.hpp>
#include <meshing.hpp>
// #include "../gprim/spline.hpp"

View File

@ -6,6 +6,7 @@
#include <meshing.hpp>
#include <geometry2d.hpp>
#include <csg2d.hpp>
using namespace netgen;
@ -265,6 +266,7 @@ DLL_HEADER void ExportGeom2d(py::module &m)
{
double len = self.splines[i]->Length();
int n = floor(len/(0.05*min(xdist,ydist)));
n = max(3, n);
lst.push_back(self.splines[i]->StartPI());
for (int j = 1; j < n; j++){
lst.push_back(self.splines[i]->GetPoint(j*1./n));
@ -395,6 +397,36 @@ DLL_HEADER void ExportGeom2d(py::module &m)
.def("_SetDomainTensorMeshing", &SplineGeometry2d::SetDomainTensorMeshing)
;
py::class_<Solid2d>(m, "Solid2d")
.def(py::init<>())
.def(py::init<std::string>())
.def_readwrite("name", &Solid2d::name)
.def("__mul__", [](Solid2d & self, Solid2d & other) { return self*other; })
.def("__add__", [](Solid2d & self, Solid2d & other) { return self+other; })
.def("__sub__", [](Solid2d & self, Solid2d & other) { return self-other; })
.def("Append", &Solid2d::Append)
.def("SetBC", &Solid2d::SetBC)
;
py::class_<Polygon2d>(m, "Polygon2d")
.def(py::init<>())
.def("SetBC", &Polygon2d::SetBC)
.def("Append", [](Polygon2d & self, double x, double y)
{
self.Append({x,y});
})
;
m.def("Rectangle", &Rectangle);
m.def("Circle", &Circle);
py::class_<CSG2d>(m, "CSG2d")
.def(py::init<>())
.def("GenerateSplineGeometry", &CSG2d::GenerateSplineGeometry)
.def("Add", &CSG2d::Add)
;
}
PYBIND11_MODULE(libgeom2d, m) {

View File

@ -251,239 +251,4 @@ int PTRIANGLE2D :: IsIn (const Point2d & p) const
}
#endif
Polygon2d :: Polygon2d ()
{
;
}
Polygon2d :: ~Polygon2d ()
{
;
}
void Polygon2d :: AddPoint (const Point2d & p)
{
points.Append(p);
}
double Polygon2d :: HArea () const
{
int i;
double ar = 0;
for (i = 1; i <= points.Size(); i++)
{
const Point2d & p1 = points.Get(i);
const Point2d & p2 = points.Get(i%points.Size()+1);
ar +=
(p2.X()-p1.X()) * p1.Y() -
(p2.Y()-p1.Y()) * p1.X();
}
return ar/2;
/*
CURSOR c;
double ar = 0;
Point2d * p1, * p2, p0 = Point2d(0, 0);
Vec2d v1, v2 = Vec2d(1, 0);
p2 = points[points.Last()];
for (c = points.First(); c != points.Head(); c++)
{
p1 = p2;
p2 = points[c];
ar += Cross ( (*p2-*p1), (*p1 - p0));
}
return ar / 2;
*/
}
int Polygon2d :: IsOn (const Point2d & p) const
{
int i;
for (i = 1; i <= points.Size(); i++)
{
const Point2d & p1 = points.Get(i);
const Point2d & p2 = points.Get(i%points.Size()+1);
if (IsOnLine (Line2d(p1, p2), p)) return 1;
}
return 0;
/*
CURSOR c;
Point2d * p1, * p2;
p2 = points[points.Last()];
for (c = points.First(); c != points.Head(); c++)
{
p1 = p2;
p2 = points[c];
if (IsOnLine (Line2d(*p1, *p2), p)) return 1;
}
return 0;
*/
}
int Polygon2d :: IsIn (const Point2d & p) const
{
int i;
double sum = 0, ang;
for (i = 1; i <= points.Size(); i++)
{
const Point2d & p1 = points.Get(i);
const Point2d & p2 = points.Get(i%points.Size()+1);
ang = Angle ( (p1 - p), (p2 - p) );
if (ang > M_PI) ang -= 2 * M_PI;
sum += ang;
}
return fabs(sum) > M_PI;
/*
CURSOR c;
Point2d * p1, * p2;
double sum = 0, ang;
p2 = points[points.Last()];
for (c = points.First(); c != points.Head(); c++)
{
p1 = p2;
p2 = points[c];
ang = Angle ( (*p1 - p), (*p2 - p) );
if (ang > M_PI) ang -= 2 * M_PI;
sum += ang;
}
return fabs(sum) > M_PI;
*/
}
int Polygon2d :: IsConvex () const
{
/*
Point2d *p, *pold, *pnew;
char cw;
CURSOR c;
if (points.Length() < 3) return 0;
c = points.Last();
p = points[c];
c--;
pold = points[c];
pnew = points[points.First()];
cw = ::CW (*pold, *p, *pnew);
for (c = points.First(); c != points.Head(); c++)
{
pnew = points[c];
if (cw != ::CW (*pold, *p, *pnew))
return 0;
pold = p;
p = pnew;
}
*/
return 0;
}
int Polygon2d :: IsStarPoint (const Point2d & p) const
{
/*
Point2d *pnew, *pold;
char cw;
CURSOR c;
if (points.Length() < 3) return 0;
pold = points[points.Last()];
pnew = points[points.First()];
cw = ::CW (p, *pold, *pnew);
for (c = points.First(); c != points.Head(); c++)
{
pnew = points[c];
if (cw != ::CW (p, *pold, *pnew))
return 0;
pold = pnew;
}
return 1;
*/
return 0;
}
Point2d Polygon2d :: Center () const
{
/*
double ai, a = 0, x = 0, y = 0;
Point2d * p, *p2;
Point2d p0 = Point2d(0, 0);
CURSOR c;
p2 = points[points.Last()];
for (c = points.First(); c != points.Head(); c++)
{
p = points[c];
ai = Cross (*p2 - p0, *p - p0);
x += ai / 3 * (p2->X() + p->X());
y += ai / 3 * (p2->Y() + p->Y());
a+= ai;
p2 = p;
}
if (a != 0)
return Point2d (x / a, y / a);
else
return Point2d (0, 0);
*/
return Point2d (0, 0);
}
Point2d Polygon2d :: EqualAreaPoint () const
{
/*
double a11 = 0, a12 = 0, a21= 0, a22 = 0;
double b1 = 0, b2 = 0, dx, dy;
double det;
Point2d * p, *p2;
CURSOR c;
p = points[points.Last()];
for (c = points.First(); c != points.Head(); c++)
{
p2 = p;
p = points[c];
dx = p->X() - p2->X();
dy = p->Y() - p2->Y();
a11 += sqr (dy);
a12 -= dx * dy;
a21 -= dx * dy;
a22 += sqr (dx);
b1 -= dy * (p->X() * p2->Y() - p2->X() * p->Y());
b2 -= dx * (p->Y() * p2->X() - p2->Y() * p->X());
}
det = a11 * a22 - a21 * a12;
if (det != 0)
return Point2d ( (b1 * a22 - b2 * a12) / det,
(a11 * b2 - a21 * b1) / det);
else
return Point2d (0, 0);
*/
return Point2d (0, 0);
}
}

View File

@ -609,40 +609,6 @@ namespace netgen
#endif
class Polygon2d
{
protected:
NgArray<Point2d> points;
public:
Polygon2d ();
~Polygon2d ();
void AddPoint (const Point2d & p);
int GetNP() const { return points.Size(); }
void GetPoint (int i, Point2d & p) const
{ p = points.Get(i); }
void GetLine (int i, Point2d & p1, Point2d & p2) const
{ p1 = points.Get(i); p2 = points.Get(i%points.Size()+1); }
double Area () const { return fabs (HArea()); }
int CW () const { return HArea() > 0; }
int CCW () const { return HArea() < 0; }
int IsOn (const Point2d & p) const;
int IsIn (const Point2d & p) const;
int IsConvex () const;
int IsStarPoint (const Point2d & p) const;
Point2d Center() const;
Point2d EqualAreaPoint () const;
private:
double HArea () const;
};
/** Cheap approximation to atan2.
A monotone function of atan2(x,y) is computed.
*/

View File

@ -196,6 +196,10 @@ namespace netgen
{
ar & p1 & p2 & p3 & weight & proj_latest_t;
}
///
double GetWeight () const { return weight; }
void SetWeight (double w) { weight = w; }
///
virtual Point<D> GetPoint (double t) const;
///
virtual Vec<D> GetTangent (const double t) const;

View File

@ -1,4 +1,4 @@
from .libngpy._geom2d import SplineGeometry
from .libngpy._geom2d import SplineGeometry, Solid2d, Polygon2d, CSG2d, Rectangle, Circle
from .meshing import meshsize
unit_square = SplineGeometry()