#ifndef NETGEN_CSG2D_HPP_INCLUDED #define NETGEN_CSG2D_HPP_INCLUDED #include #include "geometry2d.hpp" namespace netgen { using namespace ngcore; using netgen::Point; using netgen::Vec; using Spline = SplineSeg3<2>; using netgen::Box; 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]); } // compute weight of spline such that p lies on it void ComputeWeight( Spline & s, Point<2> p ); 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 }; inline constexpr const double MAXH_DEFAULT{1e99}; inline const string BC_DEFAULT{""}; inline const string MAT_DEFAULT{""}; struct EdgeInfo { optional> control_point = nullopt; // for spline segments double maxh = MAXH_DEFAULT; string bc = BC_DEFAULT; EdgeInfo() = default; EdgeInfo(Point<2> p) : control_point(p) {} EdgeInfo(double h) : maxh(h) {} EdgeInfo(string s) : bc(s) {} EdgeInfo(optional> p, double h, string s) : control_point(p), maxh(h), bc(s) {} void Assign( EdgeInfo other ) { if(other.control_point != nullopt) control_point = other.control_point; if(other.bc != BC_DEFAULT) bc = other.bc; if(other.maxh != MAXH_DEFAULT) maxh = min(maxh, other.maxh); } }; struct Vertex : Point<2> { Vertex (Point<2> p) : Point<2>(p) {} Vertex (const Vertex & v) : Point<2>(v) { spline = v.spline; info = v.info; is_source = true; } Vertex * prev = nullptr; Vertex * next = nullptr; unique_ptr 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" // In case the edge this - next is curved, store the spline information here optional spline = nullopt; EdgeInfo info; DLL_HEADER 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 Loop { unique_ptr first = nullptr; unique_ptr> bbox = nullptr; Loop() = default; Loop(const Loop & p) : first(nullptr) { for(auto v : p.Vertices(ALL)) AppendVertex(*v); } Loop(Loop && p) = default; Loop & operator=(Loop && p) = default; Loop & operator=(const Loop & p) { // static Timer t("Loop::operator="); RegionTimer rt(t); first = nullptr; if(p.first) { size_t n = p.Size(); Array> new_verts(n); { size_t i = 0; for(const auto v : p.Vertices(ALL)) new_verts[i++] = make_unique(*v); } for(auto i : IntRange(n-1)) { Vertex * v = new_verts[i].get(); Vertex * vn = new_verts[i+1].get(); v->next = vn; vn->prev = v; } Vertex * vfirst = new_verts[0].get(); Vertex * vlast = new_verts[n-1].get(); vfirst->prev = vlast; vlast->next = vfirst; for(auto i : IntRange(1,n)) new_verts[n-1-i]->pnext = std::move(new_verts[n-i]); first = std::move(new_verts[0]); } bbox = nullptr; return *this; } void Clear() { first = nullptr; } Vertex & AppendVertex(const Vertex & v) { auto & vnew = Append( static_cast>(v), true ); vnew.info = v.info; if(v.spline) vnew.spline = *v.spline; if(bbox) bbox->Add(v); return vnew; } Vertex & Append(Point<2> p, bool source = false) { Vertex * vnew; if(first==nullptr) { first = make_unique(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; if(bbox) bbox->Add(p); 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); bbox.reset(); } bool IsInside( Point<2> r ) const; bool IsLeftInside( const Vertex & p0 ); bool IsRightInside( const Vertex & p0 ); 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->info.bc = bc; } size_t Size() const { if(first==nullptr) return 0; size_t cnt = 0; for(auto v : Vertices(ALL)) cnt++; return cnt; } const Box<2> & GetBoundingBox() { if(bbox==nullptr) { static Timer tall("Loop::GetBoundingBox"); RegionTimer rt(tall); bbox = make_unique>(Box<2>::EMPTY_BOX); for(auto v : Vertices(ALL)) bbox->Add(*v); } return *bbox; } }; struct Solid2d { Array polys; string name = MAT_DEFAULT; Solid2d() = default; Solid2d(string name_) : name(name_) {} DLL_HEADER Solid2d(const Array, EdgeInfo>> & points, string name_=MAT_DEFAULT, string bc_=BC_DEFAULT); Solid2d(Solid2d && other) = default; Solid2d(const Solid2d & other) = default; DLL_HEADER Solid2d operator+(const Solid2d & other) const; DLL_HEADER Solid2d operator*(const Solid2d & other) const; DLL_HEADER Solid2d operator-(const Solid2d & other) const; Solid2d& operator=(Solid2d && other) = default; Solid2d& operator=(const Solid2d & other) = default; DLL_HEADER Solid2d& operator+=(const Solid2d & other); DLL_HEADER Solid2d& operator*=(const Solid2d & other); DLL_HEADER Solid2d& operator-=(const Solid2d & other); void Append( const Loop & poly ) { polys.Append(poly); } bool IsInside( Point<2> r ) const; bool IsLeftInside( const Vertex & p0 ); bool IsRightInside( const Vertex & p0 ); template Solid2d & Transform( const TFunc & func ) { for(auto & poly : polys) for(auto v : poly.Vertices(ALL)) { auto p = func(*v); (*v)[0] = p[0]; (*v)[1] = p[1]; if(v->spline) { auto &s = *v->spline; auto pmid = func(s.GetPoint(0.5)); s = Spline(func(s.StartPI()), func(s.TangentPoint()), func(s.EndPI())); ComputeWeight(s, pmid); } } return *this; } Solid2d & Move( Vec<2> v ); Solid2d & Scale( double sx, double sy=0.0 ); Solid2d & RotateRad( double ang, Point<2> center = {0,0} ); Solid2d & RotateDeg( double ang, Point<2> center = {0,0} ) { return RotateRad( ang/180.*M_PI ); } Solid2d & BC(string bc) { for(auto & p : polys) for(auto v : p.Vertices(ALL)) v->info.bc = bc; return *this; } Solid2d & Maxh(double maxh) { for(auto & p : polys) for(auto v : p.Vertices(ALL)) v->info.maxh = maxh; return *this; } Solid2d & Mat(string mat) { name = mat; return *this; } Box<2> GetBoundingBox() const; }; class CSG2d { public: Array solids; void Add ( Solid2d s ) { solids.Append(s); } DLL_HEADER shared_ptr GenerateSplineGeometry(); DLL_HEADER shared_ptr GenerateMesh(MeshingParameters & mp); }; DLL_HEADER Solid2d Circle( Point<2> center, double r, string name="", string bc=""); DLL_HEADER Solid2d Rectangle( Point<2> p0, Point<2> p1, string mat=MAT_DEFAULT, string bc=BC_DEFAULT ); DLL_HEADER void AddIntersectionPoints ( Solid2d & s1, Solid2d & s2 ); DLL_HEADER Solid2d ClipSolids ( const Solid2d & s1, const Solid2d & s2, char op); DLL_HEADER Solid2d ClipSolids ( const Solid2d & s1, Solid2d && s2, char op); DLL_HEADER Solid2d ClipSolids ( Solid2d && s1, const Solid2d & s2, char op); DLL_HEADER Solid2d ClipSolids ( Solid2d && s1, Solid2d && s2, char op); DLL_HEADER IntersectionType intersect(const Point<2> P1, const Point<2> P2, const Point<2> Q1, const Point<2> Q2, double& alpha, double& beta); } #endif // NETGEN_CSG2D_HPP_INCLUDED