diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h
index f367de4e8bf2bdb5e1df85dedb828c1cf8b0abfa..8bc1bf7536ddef46bbddf03ff0715cabca2cabcb 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronGeneric.h
@@ -21,7 +21,7 @@ class GeoPolyhedronGeneric : public GeoShape {
GeoPolyhedronGeneric(){};
// Constructor
- GeoPolyhedronGeneric(const GeoShape*);
+ GeoPolyhedronGeneric(const GeoShape*, const GeoTrf::Transform3D& trf = GeoTrf::Transform3D::Identity(), int polyhedrize = 0);
// Copy constructor
GeoPolyhedronGeneric(const GeoPolyhedronGeneric& gpg);
@@ -84,6 +84,9 @@ class GeoPolyhedronGeneric : public GeoShape {
// check of the hermeticity of the polyhedron (all edge need to be mirrored), plus at least four faces required, false otherwise
bool isValid ( bool debugMode = false) const;
+ void setStatus( int status) { m_status = status;}
+ int getStatus() const { return m_status; }
+
// Returns coordinates of the specified vertex
GeoTrf::Vector3D getVertex(unsigned int i) const {
return m_Vertices.at(i);
@@ -114,7 +117,6 @@ class GeoPolyhedronGeneric : public GeoShape {
virtual bool isPolyhedron () const {
return true;
}
-
virtual ~GeoPolyhedronGeneric() = default;
protected:
@@ -122,6 +124,7 @@ protected:
private:
static const std::string s_classType;
static const ShapeType s_classTypeID;
+ GeoTrf::Transform3D m_transform{GeoTrf::Transform3D::Identity()};
std::vector<GeoTrf::Vector3D> m_Vertices{};
std::vector< std::vector<size_t> > m_Faces{};
@@ -140,7 +143,10 @@ protected:
double getRefSpan(int face = -1) const;
// for debugging purposes
- int m_origin{0};
+ int m_origin{-1};
+
+ // status : -1 unchecked, 0 valid (may be empty), 1 non-polyhedrizable component, 2 failure to construct
+ int m_status{-1};
};
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h
index c0b1fe90fcf5c84fd07f86e5f1e1361421fa3332..f0252d310d1896cf1f6a62805f7dd7dfa6ef6b55 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedronManipulator.h
@@ -12,8 +12,8 @@
*
*/
-#include "GeoModelKernel/GeoShape.h"
#include <GeoModelKernel/GeoDefinitions.h>
+#include "GeoModelKernel/GeoShape.h"
#include "GeoModelKernel/GeoShapeShift.h"
#include "GeoModelKernel/GeoPolyhedronGeneric.h"
#include <string>
@@ -21,13 +21,13 @@
class GeoShapeIntersection;
class GeoShapeUnion;
class GeoShapeSubtraction;
-class GeoPolyhedronGeneric;
struct ShiftedGeoShape {
const GeoShape* shape{nullptr};
const GeoPolyhedronGeneric* phed{nullptr};
GeoTrf::Transform3D transform;
double volume{-1}; // -1 not calculated/not analytical
+ int polyhedrize{16}; // >0 polyhedrize
// constructor
ShiftedGeoShape(const GeoShape* inputShape) {
shape = inputShape;
@@ -38,14 +38,15 @@ struct ShiftedGeoShape {
transform = shift->getX()*transform;
shape = shift->getOp();
}
- if (!(shape->type()=="Union" || shape->type()=="Subtraction" || shape->type()=="Intersection") ) {
+ if ( !(shape->type()=="Union" || shape->type()=="Subtraction" || shape->type()=="Intersection") ) {
volume = shape->volume();
- phed = new GeoPolyhedronGeneric(shape);
+ phed = new GeoPolyhedronGeneric(shape,GeoTrf::Transform3D::Identity(),polyhedrize);
if (!phed->isValid()) phed=nullptr; // temporary / to deal with non-polyhedrized shapes
}
}
};
+
struct NPL_vertex{
int index; bool shared{false}; int overlapping_face{-1}; bool used{false}; int onEdge{-1}; bool onEdgeDummy{false};
int shared_adjacent{-1}; bool oriented{false}; // int alternative_face;
@@ -72,7 +73,7 @@ struct NPL_edge{
struct NPL_ordered{
std::vector<size_t> vtx;
std::vector<int> edges;
- int top_overlap{-1}; bool shared{false}; bool oriented{false};
+ int top_overlap{-1}; bool shared{false}; bool oriented{false}; bool guess_orientation{false};
void reorder(size_t a) {
while ( vtx.front()!= a) {
vtx.push_back(vtx.front());
@@ -81,6 +82,17 @@ struct NPL_ordered{
edges.erase(edges.begin());
}
}
+ void reverse(std::vector<NPL_edge> plane_edges) {
+ std::vector<size_t> vtx_tmp;
+ std::vector<int> edges_tmp;
+ for (auto iv : vtx) vtx_tmp.insert(vtx_tmp.begin(),iv);
+ for (auto ie : edges) {
+ edges_tmp.insert(edges_tmp.begin(),plane_edges[ie].symmetrized);
+ }
+ edges_tmp.push_back(edges_tmp.front()); edges_tmp.erase(edges_tmp.begin()); // reorder to match vtx
+ vtx = vtx_tmp;
+ edges = edges_tmp;
+ }
};
struct NPH_plane{
@@ -94,6 +106,7 @@ struct NPH_plane{
//std::vector<NPL_ordered> shared_subface;
std::vector<NPL_ordered> subfaces;
std::vector<NPL_edge> plane_edges;
+ std::vector<NPL_edge> edge_loop_origin;
GeoTrf::Vector3D normal;
std::vector<int> overlapping;
NPH_plane(int index) { index_ph = index;}
@@ -103,7 +116,7 @@ struct NPH_vertex{
int index{-1};
int index_ph{-1};
bool shared{false};
- int onDoublEdge{-1}; int hitDoublEdge{-1};
+ //int onDoublEdge{-1}; //int hitDoublEdge{-1};
std::pair<int,int> hit_edge_edge;
std::vector<std::pair<int,int> > edges;
NPH_vertex() {hit_edge_edge = std::pair<int,int>(-1,-1);}
@@ -146,10 +159,18 @@ class GeoPolyhedronManipulator
{
public:
+ // steering routine for volume calculation with pre-defined precision
+ double volumeBoolean( const GeoShape* shape, double tolerance, int mode) const;
+
+ // resolve into generic polyhedron shape if possible, or deploy polyhedrization
+ GeoPolyhedronGeneric resolve( const GeoShape*, GeoTrf::Transform3D& trf, int polyhedrize = 0, int level = 0, bool debugMode = false );
+
GeoPolyhedronGeneric* intersection( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
GeoPolyhedronGeneric* subtraction( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
GeoPolyhedronGeneric* uni( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode);
- GeoPolyhedronGeneric* resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode);
+ GeoPolyhedronGeneric resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode);
+
+ void collectFaceNumber(const GeoShape*, int& nFaces);
private:
@@ -170,7 +191,7 @@ class GeoPolyhedronManipulator
void collectLines( std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D>& vtx, std::vector<NPH_vertex>& vlog,
std::vector<NPH_edge>& edges, bool debugMode ) const;
- std::vector<NPL_edge> evaluateVertexLogic(int plane_index, std::vector<NPH_vertex>& vlog, std::vector<NPH_edge> edges,
+ std::vector<NPL_edge> evaluateVertexLogic(int plane_index, std::vector<NPH_vertex> vlog, std::vector<NPH_edge> edges,
std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D> vtx) const;
bool insideFace(unsigned int plane_index, GeoTrf::Vector3D point, double precision = 1.e-5, bool debugMode = false) const {
@@ -206,6 +227,11 @@ class GeoPolyhedronManipulator
GeoTrf::Vector3D areaVec( std::vector<size_t> loop, std::vector< GeoTrf::Vector3D > vtx) const;
+ int nextOri(int input, std::vector<NPL_edge> plane_edges, std::vector<NPL_edge> edge_loop_origin, std::vector<int> availableEdges, bool previous = false) const;
+
+ // for debugging purpose only
+ std::vector<std::pair<int,int> > missing_edges( const GeoPolyhedronGeneric gpg) const;
+
};
#endif
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx
index 28e9b1825ca4e8bf3b5ed0dd510b8c1d460d0141..fb03ed4c47a3ea05dc2879ada3d9c97353e3e51d 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronGeneric.cxx
@@ -7,7 +7,10 @@
#include "GeoModelKernel/GeoBox.h"
#include "GeoModelKernel/GeoTrd.h"
#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoTubs.h"
+#include "GeoModelKernel/GeoPcon.h"
#include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include "GeoModelKernel/GeoShapeShift.h"
#include <cmath>
#include <stdexcept>
#include <iostream>
@@ -17,9 +20,10 @@
const std::string GeoPolyhedronGeneric::s_classType = "PolyhedronGeneric";
const ShapeType GeoPolyhedronGeneric::s_classTypeID = 0x25;
-GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoShape* shape) {
+GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoShape* shape, const GeoTrf::Transform3D& transform, int polyhedrize) {
std::vector<GeoTrf::Vector3D> vtx;
+ m_transform = transform;
if (shape->type() == "Box") {
const GeoBox* box = dynamic_cast<const GeoBox*> (shape);
@@ -78,9 +82,11 @@ GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoShape* shape) {
for (unsigned int i=0; i<getNFaces(); i++)
m_Normals.push_back(
((m_Vertices[m_Faces[i][0]]-m_Vertices[m_Faces[i][1]]).cross(m_Vertices[m_Faces[i][2]]-m_Vertices[m_Faces[i][1]])).normalized() );
- }
- if (shape->type() == "SimplePolygonBrep") {
+ m_status = isValid() ? 0 : 2;
+ //if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from Box/Trd not valid" << std::endl;
+
+ } else if (shape->type() == "SimplePolygonBrep") {
const GeoSimplePolygonBrep* spb = dynamic_cast<const GeoSimplePolygonBrep*> (shape);
double hz = spb->getDZ();
@@ -133,81 +139,184 @@ GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoShape* shape) {
}
}
}
-
- if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from SimplePolygonBrep not valid" << std::endl;
+ m_status = isValid() ? 0 : 2;
+ //if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from SimplePolygonBrep not valid" << std::endl;
+
m_origin = 3;
}
+ else if ( (shape->type() == "Tube" || shape->type() == "Tubs") ) {
- if (shape->type() == "Tube") { // polyhedrize : TODO : make the precision tunable
-
- const GeoTube* tube = dynamic_cast<const GeoTube*> (shape);
- double hz = tube->getZHalfLength();
- double rmin = tube->getRMin();
- double rmax = tube->getRMax();
+ m_origin = 4;
- int nsplit = 16;
+ if (polyhedrize == 0) {
+ m_status = 1; return;
+ }
+ const GeoTube* tube = (shape->type() == "Tube") ? dynamic_cast<const GeoTube*> (shape) : nullptr;
+ const GeoTubs* tubs = (shape->type() == "Tubs") ? dynamic_cast<const GeoTubs*> (shape) : nullptr;
+
+ double hz = tube ? tube->getZHalfLength() : tubs->getZHalfLength();
+ double rmin = tube ? tube->getRMin() : tubs->getRMin();
+ double rmax = tube ? tube->getRMax() : tubs->getRMax();
+ double sPhi = tubs ? tubs->getSPhi() : 0.;
+ double dPhi = tubs ? tubs->getDPhi() : 2*acos(-1.);
+ bool sector = dPhi < 2*acos(-1.);
+
+ int nsplit = dPhi>acos(-1.) ? polyhedrize : int(0.5*polyhedrize);
+ if ( nsplit < 2 ) nsplit = 2;
+
// build outer shell
- double ff = acos(-1.)/nsplit;
+ double ff = 0.5*dPhi/nsplit;
double rx = rmax * sqrt( ff / sin(ff) / cos(ff) );
- for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(-i*2*ff),rx*sin(-i*2*ff),hz));
- for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(-i*2*ff),rx*sin(-i*2*ff),-hz));
-
- std::vector<size_t> f0; std::vector<size_t> f1;
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(sPhi+dPhi-i*2*ff),rx*sin(sPhi+dPhi-i*2*ff),hz));
+ std::vector<size_t> f0;
for (unsigned int i=0; i < nsplit; i++) f0.push_back(i);
- for (unsigned int i=2*nsplit-1; i >= nsplit; i--) f1.push_back(i);
-
- if (rmin>0) {
-
- f0.push_back(0); f1.push_back(2*nsplit-1); // needed for combined loop
-
+
+ // sector needs the closing vertex,too
+ if ( sector ) {
+ f0.push_back(vtx.size()); vtx.push_back(GeoTrf::Vector3D(rx*cos(sPhi),rx*sin(sPhi),hz));
+ if (!(rmin>0) ) { // need the central point
+ f0.push_back(vtx.size()); vtx.push_back(GeoTrf::Vector3D(0.,0.,hz));
+ }
+ } else if (rmin>0) f0.push_back(0); // needed for combined loop
+
+ int nv = vtx.size(); // outer ring vertex count
+ if (rmin>0) { // inner ring mirrored
double rn = rmin * sqrt( ff / sin(ff) / cos(ff) );
- for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rn*cos(-i*2*ff),rn*sin(-i*2*ff),-hz));
- for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rn*cos(-i*2*ff),rn*sin(-i*2*ff),hz));
-
- for (unsigned int i=2*nsplit; i < 3*nsplit; i++) f1.push_back(i);
- for (unsigned int i=4*nsplit-1; i >= 3*nsplit; i--) f0.push_back(i);
-
- f0.push_back(4*nsplit-1); f1.push_back(2*nsplit); // needed for combined loop
+ for (int i = nv-1; i >= 0; i--) {
+ vtx.push_back(GeoTrf::Vector3D(rn/rx*vtx[i].x(),rn/rx*vtx[i].y(),hz));
+ }
+ f0.push_back(2*nv-1);
+ for (unsigned int i = nv; i< 2*nv; i++) f0.push_back(i);
}
-
addFace(f0); addNormal(GeoTrf::Vector3D(0.,0.,1.));
+
+ int mv = vtx.size();
+
+ // mirror the face
+ int ntop = f0.size();
+ for (unsigned int i=0; i<mv; i++) {
+ vtx.push_back(GeoTrf::Vector3D(vtx[i].x(),vtx[i].y(),-hz));
+ }
+ std::vector<size_t> f1;
+ for (int i=ntop-1; i>=0; i--) f1.push_back(f0[i]+mv);
addFace(f1); addNormal(GeoTrf::Vector3D(0.,0.,-1.));
addVertices(vtx);
- for (unsigned int i=0; i<nsplit; i++) {
+ for (unsigned int i=0; i<ntop; i++) {
+ if (rmin>0 && !sector && (i==nsplit || i==ntop-1) ) continue; // dummy lateral faces
std::vector<size_t> f; f.clear();
- f.push_back( i == nsplit-1 ? 0 : i+1); f.push_back(i); f.push_back(nsplit+i); f.push_back( i == nsplit-1 ? nsplit : nsplit+i+1);
- addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ f.push_back( i == ntop-1 ? f0[0] : f0[i+1]); f.push_back(f0[i]); f.push_back(f0[i]+mv); f.push_back( i == ntop-1 ? f0[0]+mv : f0[i+1]+mv);
+ addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
}
- if (rmin>0) {
- std::vector<size_t> f; f.clear();
- f.push_back(3*nsplit); f.push_back(0); f.push_back(nsplit); f.push_back(2*nsplit);
- addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ fillSpans();
- for (unsigned int i=0; i<nsplit; i++) {
+ m_status = isValid() ? 0 : 2;
+ //if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from Tube/Tubs not valid" << std::endl;
+
+ }
+ else if ( shape->type() == "Pcon") {
+
+ m_origin = 5;
+
+ if ( polyhedrize == 0 ) {
+ m_status = 1; return;
+ }
+
+ const GeoPcon* pcon = dynamic_cast<const GeoPcon*> (shape);
+
+ double sPhi = pcon->getSPhi();
+ double dPhi = pcon->getDPhi();
+ bool sector = dPhi < 2*acos(-1.);
+ unsigned int nplanes = pcon->getNPlanes();
+ double pord = ( pcon->getZPlane(0) < pcon->getZPlane(nplanes-1) ) ? 1. : -1.;
+
+ int nsplit = dPhi>acos(-1.) ? polyhedrize : int(0.5*polyhedrize);
+ if ( nsplit < 2 ) nsplit = 2;
+
+ double ff = 0.5*dPhi/nsplit;
+
+ std::vector<size_t> f; int nv = 0; int ntop = 0; int mv = 0;
+
+ for ( unsigned int iz = 0; iz < nplanes; iz++) {
+
+ double z = pcon->getZPlane(iz);
+ double rmin = pcon->getRMinPlane(iz);
+ double rmax = pcon->getRMaxPlane(iz);
+
+ // build outer shell
+ double rx = rmax * sqrt( ff / sin(ff) / cos(ff) );
+ for (unsigned int i=0; i < nsplit; i++) vtx.push_back(GeoTrf::Vector3D(rx*cos(sPhi+dPhi-i*2*ff),rx*sin(sPhi+dPhi-i*2*ff),z));
+
+ if (sector) {
+ vtx.push_back(GeoTrf::Vector3D(rx*cos(sPhi),rx*sin(sPhi),z)); // sector needs the closing vertex,too
+ if (!(rmin>0) ) vtx.push_back(GeoTrf::Vector3D(0.,0.,z)); // the central point needed
+ }
+
+ if (iz==0 || iz == nplanes-1 ) {
f.clear();
- f.push_back( i == 3*nsplit-1 ? 2*nsplit : 2*nsplit+i+1); f.push_back(2*nsplit+i); f.push_back(3*nsplit+i); f.push_back( i == 3*nsplit-1 ? 3*nsplit : 3*nsplit+i+1);
- addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ for (unsigned int i=0; i < nsplit; i++) f.push_back(vtx.size()-nsplit+i);
+ if ( sector ) {
+ f.push_back(vtx.size());
+ if (!(rmin>0) ) f.push_back(vtx.size());
+ } else if (rmin>0) f.push_back(vtx.size()-nsplit); // needed for combined loop
+ if (iz==0) nv = vtx.size(); // outer ring vertex count
+ }
+
+ // inner ring mirrored
+ if (rmin>0) {
+ double rn = rmin * sqrt( ff / sin(ff) / cos(ff) );
+ for (int i = vtx.size()-1; i >= vtx.size()-nv; i--) {
+ vtx.push_back(GeoTrf::Vector3D(rn/rx*vtx[i].x(),rn/rx*vtx[i].y(),z));
+ }
+
+ if (iz==0 || iz == nplanes-1) {
+ f.push_back(vtx.size()-1);
+ for (unsigned int i = vtx.size()-nv; i< vtx.size(); i++) f.push_back(i);
+ }
+ }
+
+ if (iz==0) {
+ ntop = f.size();
+ mv = vtx.size();
}
- f.clear();
- f.push_back(2*nsplit); f.push_back(nsplit); f.push_back(0); f.push_back(3*nsplit);
- addFace(f); addNormal( (m_Vertices[f[0]]-m_Vertices[f[1]]).cross(m_Vertices[f[2]]-m_Vertices[f[1]]).normalized() );
+ if (iz==0 || iz == nplanes-1) {
+ if ( (pord>0 && iz==nplanes-1) || (pord<0 && iz==0) ) {
+ addFace(f); addNormal(GeoTrf::Vector3D(0.,0.,+1.));
+ } else {
+ std::vector<size_t> fi;
+ for (int i=ntop-1; i>=0; i--) fi.push_back(f[i]);
+ addFace(fi); addNormal(GeoTrf::Vector3D(0.,0.,-1.));
+ }
+ }
+
+ if ( iz>0 ) {
+
+ int nvx = vtx.size();
+ for (unsigned int i=0; i<ntop; i++) {
+ if (rmin>0 && !sector && (i==nsplit || i==ntop-1) ) continue; // dummy lateral faces
+ f.clear();
+ f.push_back( i == ntop-1 ? nvx-ntop : nvx-ntop+i+1); f.push_back(nvx-ntop+i);
+ f.push_back(nvx-2*ntop+i); f.push_back( i == ntop-1 ? nvx-2*ntop : nvx-2*ntop+i+1);
+ addFace(f); addNormal( (vtx[f[0]]-vtx[f[1]]).cross(vtx[f[2]]-vtx[f[1]]).normalized() );
+ }
+ }
}
+
+ addVertices(vtx);
fillSpans();
- std::cout << "tube->polyhedron:rmin,rmax,hz:"<<rmin <<","<<rmax<<","<<hz<<":valid:" << this->isValid() <<":volume:" << this->volume() << std::endl;
-
- m_origin = 4;
-
- }
+ m_status = isValid() ? 0 : 2;
+ //if (!isValid() ) std::cout << " GeoPolyhedronGeneric constructor from Pcon not valid" << std::endl;
+ } else m_status = 1;
+
if (m_face_span.size()<m_Faces.size()) fillSpans();
}
@@ -218,6 +327,7 @@ GeoPolyhedronGeneric::GeoPolyhedronGeneric(const GeoPolyhedronGeneric& gpg) {
m_Faces = gpg.m_Faces;
m_Normals = gpg.m_Normals;
m_origin = gpg.m_origin; // only used for debugging of initial conversion from GeoShape
+ m_status = gpg.m_status;
if (isValid()) fillSpans();
@@ -345,6 +455,15 @@ bool GeoPolyhedronGeneric::contains ( GeoTrf::Vector3D point, double precision,
void GeoPolyhedronGeneric::addVertices(std::vector<GeoTrf::Vector3D> vertices)
{
m_Vertices = vertices;
+
+ bool isIdentity = true;
+ if (m_transform.translation().norm() > std::numeric_limits<float>::epsilon() ) isIdentity = false;
+ if (isIdentity) {
+ const GeoTrf::RotationMatrix3D& rot = m_transform.linear();
+ for (int i=0; i<3; i++) if ( std::abs(GeoTrf::Vector3D::Unit(i).dot(rot*GeoTrf::Vector3D::Unit(i)) - 1 ) > std::numeric_limits<float>::epsilon() ) isIdentity = false;
+ }
+
+ if (!isIdentity) for ( int i = 0; i < m_Vertices.size(); i++) m_Vertices[i] = m_transform * m_Vertices[i];
}
void GeoPolyhedronGeneric::addFace(std::vector<size_t> faceVtx)
@@ -530,7 +649,8 @@ std::vector< std::pair< double, std::pair<size_t,size_t> > > GeoPolyhedronGeneri
bool GeoPolyhedronGeneric::isValid( bool debugMode ) const
{
- if (getNFaces()<4) return false;
+ // if (getNFaces()<4) return false;
+ if ( m_status > 0 ) return false;
// loop over edges, check that each appears with its mirrored version ( attached faces )
for ( unsigned int ic=0; ic < getNFaces(); ic++) {
@@ -643,6 +763,8 @@ void GeoPolyhedronGeneric::cleanAfterBooleanOperation() {
bool GeoPolyhedronGeneric::construct(bool debugMode) {
fillSpans();
+ m_origin = 0;
+ m_status = isValid(debugMode) ? 0 : 2;
return isValid(debugMode);
}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx
index b1632ff9d48df1bb948d1a2579db2b30203a5792..3f0eb2b423af090bf505f297e60cb8b11fdbfa1e 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedronManipulator.cxx
@@ -4,25 +4,29 @@
#include "GeoModelKernel/GeoPolyhedronManipulator.h"
#include "GeoModelKernel/GeoShape.h"
-#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoShapeShift.h"
#include "GeoModelKernel/GeoShapeUnion.h"
#include "GeoModelKernel/GeoShapeSubtraction.h"
-#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoTubs.h"
+#include "GeoModelKernel/GeoPcon.h"
#include <cmath>
#include <cstdint>
#include <iostream>
#include <utility>
GeoPolyhedronGeneric* GeoPolyhedronManipulator::intersection( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
-{ return resolveBoolean( phA, phB, 0, debugMode); }
+{ return new GeoPolyhedronGeneric(resolveBoolean( phA, phB, 0, debugMode)); }
GeoPolyhedronGeneric* GeoPolyhedronManipulator::subtraction( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
-{ return resolveBoolean( phA, phB, 1, debugMode); }
+{ return new GeoPolyhedronGeneric(resolveBoolean( phA, phB, 1, debugMode)); }
GeoPolyhedronGeneric* GeoPolyhedronManipulator::uni( const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, bool debugMode)
-{ return resolveBoolean( phA, phB, 2, debugMode); }
+{ return new GeoPolyhedronGeneric(resolveBoolean( phA, phB, 2, debugMode)); }
-GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode) {
+GeoPolyhedronGeneric GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhedronGeneric* phA, const GeoPolyhedronGeneric* phB, int modeBoolean, bool debugMode) {
// modeBoolean : 0 intersection (A and B): shared(A,B) non-overlapping + shared A overlapping
// 1 subtraction (A and not B): non-shared(A) + inverted shared (B)
// 2 union (A or B) : non-shared(A,B) + shared A overlapping
@@ -39,12 +43,16 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
for (unsigned int iv=0; iv < nVA; iv++) {
vtx.push_back( phA->getVertex(iv) );
NPH_vertex phv; phv.index=vlog.size(); phv.index_ph=0; phv.shared = phB->contains( vtx.back(), 1.e-3 ); vlog.push_back( phv );
+
}
for (unsigned int ip=0; ip < phA->getNFaces(); ip++) {
planes.push_back(NPH_plane(0));
planes.back().vertices = phA->getFace(ip);
planes.back().normal = phA->faceNormal(ip);
- if (debugMode) std::cout <<"phA:face:" << ip <<":vtx:"; for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ if (debugMode) {
+ std::cout <<"phA:face:" << ip <<":vtx:";
+ for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ }
}
// phB + overlapping vertices + overlapping planes
std::vector<std::pair<int,int>> overlapping_vertices;
@@ -63,7 +71,10 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
planes.back().vertices.push_back( overlap>=0 ? overlap : vh+nVA);
}
planes.back().normal = phB->faceNormal(ip);
- if (debugMode) std::cout <<"phB:face:" << ip <<":vtx:"; for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ if (debugMode) {
+ std::cout <<"phB:face:" << ip <<":vtx:";
+ for ( auto jv : planes.back().vertices ) std::cout << jv<<","; std::cout <<"" << std::endl;
+ }
// overlapping planes
for (unsigned int ipA=0; ipA < nPA; ipA++) { NPH_plane plA = planes[ipA];
double nn = plA.normal.dot(phB->faceNormal(ip));
@@ -104,6 +115,8 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
if (debugMode) {
std::cout << "processing plane:" << ip <<":origin:" << plane.index_ph <<":overlapping planes:";
for (auto op : plane.overlapping) std::cout << op <<","; std::cout <<" " << std::endl;
+ std::cout << "original vertices:";
+ for (auto iv : plane.vertices) std::cout << iv<<","; std::cout << ""<< std::endl;
}
for ( auto ov : plane.overlapping ) { // add shared vertices TO DO : this loop seem obsolete / superflu
for ( auto iv : planes[ov%1000].vertices ) { if (!vlog[iv].shared) continue;
@@ -183,6 +196,7 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
plane_edges.push_back(ple);
}
}
+
// loop over edges to identify vertices with protected line connection
for (unsigned int ie = 0; ie< plane_edges.size()-1; ie++) {
if ( plane_edges[ie].shared != plane_edges[ie+1].shared ) planes[ip].hit_vertices[plane_edges[ie].endpoint_hitvx.second].onEdgeProtected = true;
@@ -196,13 +210,14 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
}
+
//
// ======================== prepare active lines ( mixed edges ) ====================================
std::vector<NPL_edge> resolved_lines = evaluateVertexLogic(ip, vlog, edges, planes, vtx);
if (debugMode) {
- std::cout << "resolved lines:" << resolved_lines.size() <<":"; for (auto il : resolved_lines) std::cout << il.adjacent_face <<",";
- std::cout <<""<<std::endl;
+ std::cout << "resolved lines:" << resolved_lines.size() <<":";
+ for (auto il : resolved_lines) std::cout << "adj:"<<il.adjacent_face <<":" << il.endpoints.first << "," <<il.endpoints.second<<std::endl;
}
plane = planes[ip];
@@ -247,8 +262,6 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
edgeInPlane = false; for ( auto ep : edge.inPlane ) if (ap == ep) edgeInPlane = true;
if (!edgeInPlane && inPlane(vlog[ple_next.endpoints.second],edges,ap) ) edgeInPlane = true;
if (!edgeInPlane) pruned_lines.push_back(ap);
- if (edgeInPlane && ( plane.hit_vertices[ple.endpoint_hitvx.first].onEdgeProtected || plane.hit_vertices[ple.endpoint_hitvx.second].onEdgeProtected) )
- std::cout << "pruning removed protected line:" << ap <<":endpoints:" << ple.endpoints.first <<"," << ple.endpoints.second << std::endl;
}
int nlines = pruned_lines.size();
if ( nlines == 1 ) {
@@ -293,12 +306,16 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
plane.hit_vertices[ple.endpoint_hitvx.second].faces.first = ple.adjacent_face ;
plane.hit_vertices[ple.endpoint_hitvx.second].faces.second = plane_edges[ie].endpoint_faces.second;
- plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.second = ple_next.adjacent_face ;
- plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.first = plane_edges[ie].endpoint_faces.second;
+ if ( ple_next.endpoint_hitvx.second >=0 ) {
+ plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.second = ple_next.adjacent_face ;
+ plane.hit_vertices[ple_next.endpoint_hitvx.second].faces.first = plane_edges[ie].endpoint_faces.second;
+ }
}
} // end non/shared vertices
} // end second loop over edges
+ planes[ip].edge_loop_origin = plane_edges;
+
if (debugMode) { // print overview of shared edges
for (unsigned int ie = 0; ie<plane_edges.size(); ie++) {
if (!plane_edges[ie].shared ) continue;
@@ -314,6 +331,7 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
<< plane.hit_vertices[ple.endpoint_hitvx.second].faces_shared_loop.second << std::endl;
}
}
+
// ============== merge edges and lines ================================================
// at this point, lines are not oriented and have not been symmetrized yet
// the ordering ensures oriented edges are recovered first
@@ -442,9 +460,24 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
// =================== shared loops : base ===========================================
// select shared non-overlapping first, then add all unused shared
+
+
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ if ( plane_edges[ie].generating_edge>=0 && plane_edges[ie].dummy_overlap>=0 ) {
+ for (unsigned int je=ie+1; je<plane_edges.size(); je++) {
+ if (plane_edges[je].endpoints.first == plane_edges[ie].endpoints.first &&
+ plane_edges[je].endpoints.second == plane_edges[ie].endpoints.second && plane_edges[je].generating_edge<0 )
+ plane_edges[je].dummy_overlap = plane_edges[ie].dummy_overlap;
+ if (plane_edges[je].endpoints.first == plane_edges[ie].endpoints.second &&
+ plane_edges[je].endpoints.second == plane_edges[ie].endpoints.first && plane_edges[je].generating_edge<0 )
+ plane_edges[je].dummy_overlap = plane_edges[ie].dummy_overlap;
+ }
+ }
+ }
+
int nedge = 0; std::vector<int> edgeSelect;
- for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
+ for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
if ( !plane_edges[ie].used && plane_edges[ie].shared && plane_edges[ie].overlap<0 ) {
edgeSelect.push_back(ie); nedge++;
}
@@ -463,54 +496,120 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
edgeSelect.push_back(ie); nedge++;
}
}
-
- if (edgeSelect.size()>1) {
- std::vector<int>::iterator it = edgeSelect.begin();
- NPL_ordered loop; loop.edges.push_back(*it); plane_edges[*it].used = true;
- loop.vtx.push_back(plane_edges[*it].endpoints.first);
- loop.vtx.push_back(plane_edges[*it].endpoints.second);
- loop.shared=true; loop.oriented = plane_edges[*it].generating_edge>=0 ? true : false; // after including check on dummy shared edges, oriented edges may appear in the middle of the loop
- it = edgeSelect.begin();
- while ( loop.vtx.front() != loop.vtx.back() && it != edgeSelect.end()) {
- if ( loop.edges.size() && plane_edges[loop.edges.back()].symmetrized == (*it) ) {
- it++;
- } else if (!plane_edges[*it].used && plane_edges[*it].endpoints.first == loop.vtx.back() ) {
- //
- loop.edges.push_back(*it); loop.vtx.push_back(plane_edges[*it].endpoints.second);
- plane_edges[*it].used = true; it = edgeSelect.begin(); // std::cout << ":adding:" << (*it) << std::endl;
-
- } else if (!plane_edges[*it].used && plane_edges[*it].generating_edge<0 && plane_edges[*it].symmetrized<0
- && plane_edges[*it].endpoints.second == loop.vtx.back()) { // need to symmetrize the line now
- NPL_edge symLine = symmetrize(plane_edges[*it]); plane_edges[*it].symmetrized=plane_edges.size(); symLine.used=true; symLine.symmetrized = (*it);
- (*it) = plane_edges.size(); plane_edges.push_back(symLine); // std::cout << ":adding symmetrized:" << (*it) << std::endl;
- loop.edges.push_back(*it); loop.vtx.push_back(symLine.endpoints.second); it = edgeSelect.begin();
-
- } else it++;
+ }
+
+ /*
+ std::cout << "shared edges:"<< std::endl;
+ for (auto ie : edgeSelect) std::cout << plane_edges[ie].endpoints.first <<":" << plane_edges[ie].endpoints.second
+ <<":adj:" << plane_edges[ie].adjacent_face<<":gen.edge:" << plane_edges[ie].generating_edge <<":sym:" << plane_edges[ie].symmetrized
+ <<":shared:" << plane_edges[ie].shared <<":hitvx:" << plane_edges[ie].endpoint_hitvx.first <<","
+ << plane_edges[ie].endpoint_hitvx.second <<":dummy overlap:" << plane_edges[ie].dummy_overlap << std::endl;
+ */
+ while ( nedge>2 ) {
+ if (shared_loops.size() && shared_loops.back().top_overlap<0) std::cout << "running for multiple shared loops" << std::endl;
+ int io = -1; int ien = -1;
+ //if (edgeSelect.size()>1) {
+ std::vector<int>::iterator it = edgeSelect.begin(); while (plane_edges[*it].used || plane_edges[*it].dummy_overlap>=1000) it++; if ( it == edgeSelect.end() ) break;
+ //if ( vlog[plane_edges[*it].endpoints.first].edges.size()>3)
+ // std::cout << "WARNING: ambiguous choice of shared loop start:" << plane_edges[*it].endpoints.first <<":" << plane_edges[*it].endpoint_hitvx.first << std::endl;
+ NPL_ordered loop; loop.edges.push_back(*it); plane_edges[*it].used = true;
+ loop.vtx.push_back(plane_edges[*it].endpoints.first);
+ loop.vtx.push_back(plane_edges[*it].endpoints.second);
+ loop.shared=true; loop.oriented = plane_edges[*it].generating_edge>=0 ? true : false; // after including check on dummy shared edges, oriented edges may appear in the middle of the loop
+ it = edgeSelect.begin();
+ if (plane_edges[*it].endpoint_hitvx.first>=0 && !planes[ip].hit_vertices[plane_edges[*it].endpoint_hitvx.first].active_lines.size() )
+ ien = nextOri(loop.edges.back(),plane_edges,planes[ip].edge_loop_origin,edgeSelect,true);
+ int inext=(*it);
+ while ( loop.vtx.front() != loop.vtx.back() && inext>=0 ) {
+ inext = -1;
+ // active loop end part of original loop -> fast selection track
+ if (plane_edges[loop.edges.back()].endpoint_hitvx.second>=0 &&
+ !planes[ip].hit_vertices[plane_edges[loop.edges.back()].endpoint_hitvx.second].active_lines.size() ) {
+ inext = nextOri(loop.edges.back(),plane_edges,planes[ip].edge_loop_origin,edgeSelect);
}
-
- if (loop.vtx.front() != loop.vtx.back() || loop.edges.size()<2 ) {
- for ( auto ie : loop.edges ) plane_edges[ie].used = false; // undo loop if not closed
- } else
- shared_loops.push_back(loop);
-
- int nseed = 0; // keep trace of non-symmetrized unused shared edges
- for (auto ie : edgeSelect ) if (!plane_edges[ie].used) {
- if (plane_edges[ie].symmetrized<0 ) nseed++;
- else if (!plane_edges[plane_edges[ie].symmetrized].used) nseed++;
+ if ( inext>=0 && !plane_edges[inext].used ) {
+ loop.edges.push_back(inext); loop.vtx.push_back(plane_edges[inext].endpoints.second);
+ plane_edges[inext].used = true;
+ } else { // deal with possiblity of multiple connections
+ std::vector<int>::iterator itn = edgeSelect.begin(); std::vector<int> edgeUpdate;
+ while (itn != edgeSelect.end() ) {
+ if ( loop.edges.size() && plane_edges[loop.edges.back()].symmetrized == (*itn) ) {
+ } else if (!plane_edges[*itn].used && plane_edges[*itn].endpoints.first == loop.vtx.back() ) {
+ if ( plane_edges[*itn].generating_edge<0 && plane_edges[*itn].symmetrized<0 ) {
+ NPL_edge symLine = symmetrize(plane_edges[*itn]); plane_edges[*itn].symmetrized=plane_edges.size();
+ symLine.symmetrized = (*itn); plane_edges.push_back(symLine); edgeUpdate.push_back(plane_edges.size()-1);
+ }
+ if (inext<0) inext = *itn;
+ else {
+ if (plane_edges[inext].generating_edge >= 0 && plane_edges[*itn].generating_edge < 0
+ && !plane_edges[plane_edges[*itn].symmetrized].used ) inext = *itn;
+ //std::cout <<"multiple choice:" << plane_edges[inext].generating_edge <<":"<<plane_edges[*itn].generating_edge << std::endl;
+ }
+ } else if (!plane_edges[*itn].used && plane_edges[*itn].generating_edge<0 && plane_edges[*itn].symmetrized<0
+ && plane_edges[*itn].endpoints.second == loop.vtx.back()) { // need to symmetrize the line now
+ NPL_edge symLine = symmetrize(plane_edges[*itn]); plane_edges[*itn].symmetrized=plane_edges.size();
+ symLine.symmetrized = (*itn); plane_edges.push_back(symLine); edgeUpdate.push_back(plane_edges.size()-1);
+ if (inext<0) inext = plane_edges.size()-1;
+ else {
+ if (plane_edges[inext].generating_edge >= 0 && !plane_edges[*itn].used) inext = plane_edges.size()-1;
+ //std::cout <<"multiple choice:" << plane_edges[inext].generating_edge <<":"<<plane_edges.back().generating_edge << std::endl;
+ }
+ }
+ itn++;
}
-
- if (debugMode) {
- std::cout << "base shared loop:selected:" << edgeSelect.size() <<":loop size:vtx:"<< loop.vtx.size()
- <<":loop size:edge:" << loop.edges.size() << std::endl;
- if (nseed>0) {
- std::cout <<"WARNING:remaining seeds:" << nseed << std::endl;
- for (auto iv : planes[ip].vertices) std::cout << "vtx ori:" <<iv<<":" << vtx[iv].x() <<"," << vtx[iv].y() <<"," << vtx[iv].z() << std::endl;
- for (auto iv : planes[ip].hit_vertices) std::cout << "vtx hit:" <<iv.index<<":" << vtx[iv.index].x() <<"," << vtx[iv.index].y() <<"," << vtx[iv.index].z() << std::endl;
+ if (inext>=0) {
+ //std::cout << ":next shared edge:" << inext <<":" << plane_edges[inext].endpoints.second << std::endl;
+ loop.edges.push_back(inext); loop.vtx.push_back(plane_edges[inext].endpoints.second);
+ plane_edges[inext].used = true;
+ }
+ if (edgeUpdate.size()) { edgeSelect.insert(edgeSelect.end(),edgeUpdate.begin(),edgeUpdate.end()); edgeUpdate.clear();}
+ }
+ // double-check closing conditions
+ if (ien>=0 && loop.vtx.front() == loop.vtx.back() && loop.edges.back()!=ien) {
+ int jnext = -1;
+ if (plane_edges[loop.edges.back()].endpoint_hitvx.second>=0 &&
+ !planes[ip].hit_vertices[plane_edges[loop.edges.back()].endpoint_hitvx.second].active_lines.size() ) {
+ jnext = nextOri(loop.edges.back(),plane_edges,planes[ip].edge_loop_origin,edgeSelect);
+ }
+ if ( jnext>=0 && !plane_edges[jnext].used ) {
+ loop.edges.push_back(jnext); loop.vtx.push_back(plane_edges[jnext].endpoints.second);
+ plane_edges[jnext].used = true; it = edgeSelect.begin();
+ //std::cout << "fast search applied in closing sequence" << std::endl;
}
}
}
- }
+ if (loop.vtx.front() != loop.vtx.back() || loop.edges.size()<2 ) {
+ for ( auto ie : loop.edges ) plane_edges[ie].used = false; // undo loop if not closed
+ if (debugMode) std::cout << "undo shared loop:"; for ( auto iv : loop.vtx ) std::cout <<iv<<","; std::cout <<"" << std::endl;
+ break;
+ } else {
+ shared_loops.push_back(loop);
+ if (debugMode) {
+ if (ien >= 0) {
+ for (auto iv : loop.vtx) std::cout << iv <<","; std::cout <<""<<std::endl;
+ for (auto ie : loop.edges) std::cout << plane_edges[ie].adjacent_face <<","; std::cout <<""<<std::endl;
+ std::cout << "closing edge:" << plane_edges[loop.edges.back()].endpoints.first << ":adj:" << plane_edges[loop.edges.back()].adjacent_face << std::endl;
+ if (plane_edges[loop.edges.back()].adjacent_face != planes[ip].edge_loop_origin[ien].adjacent_face) std::cout <<"CLOSING MISMATCH" << std::endl;
+ }
+ }
+ }
+
+ nedge = 0; // keep trace of non-symmetrized unused shared edges
+ for (auto ie : edgeSelect ) if (!plane_edges[ie].used) {
+ if (plane_edges[ie].symmetrized<0 ) nedge++;
+ else if (!plane_edges[plane_edges[ie].symmetrized].used) nedge++;
+ }
+
+ if (debugMode) {
+ std::cout << "base shared loop:selected:" << edgeSelect.size() <<":loop size:vtx:"<< loop.vtx.size()
+ <<":loop size:edge:" << loop.edges.size() << std::endl;
+ if (nedge>0) {
+ std::cout <<"WARNING:remaining seeds:" << nedge << std::endl;
+ }
+ }
+ }
+
//================ symmetrized remaining active lines =====================================
for (unsigned int ie=0; ie<plane_edges.size(); ie++ ) {
@@ -543,11 +642,11 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
loopi.shared = il.shared;
// resolve orientation
if ( areaVec(il.vtx,vtx).dot(plane.normal)>=0 ) {
- il.oriented = true; loopi.oriented = true;
+ il.oriented = true; loopi.oriented = true; il.guess_orientation = true;
planes[ip].subfaces.push_back(il);
shared_inverted_subface.push_back(loopi);
} else {
- il.oriented = true; loopi.oriented = true;
+ il.oriented = true; loopi.oriented = true; loopi.guess_orientation = true;
planes[ip].subfaces.push_back(loopi);
shared_inverted_subface.push_back(il);
}
@@ -679,8 +778,10 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
if ( non_shared_loops.size() ) {
// merge inverted shared with non-shared loops
- for (unsigned int il=0; il<shared_inverted_subface.size(); il++ ) {
- non_shared_loops[0] = mergeLoops(non_shared_loops[0], shared_inverted_subface[il], vtx);
+ for (unsigned int il=0; il<shared_inverted_subface.size(); il++ ) {
+ non_shared_loops[0] = mergeLoops(non_shared_loops[0], shared_inverted_subface[il], vtx);
+ //for (auto iv : non_shared_loops[0].vtx ) std::cout << iv <<","; std::cout <<""<< std::endl;
+ //for (auto ie : non_shared_loops[0].edges ) std::cout << ie <<","; std::cout <<""<< std::endl;
}
}
@@ -706,8 +807,75 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
*/
} // end loop over planes
+
+ // make final check of face orientation ( inner loops ... )
+ for ( unsigned int ip=0; ip<planes.size(); ip++ ) { NPH_plane plane = planes[ip];
+ for (unsigned int il=0; il<plane.subfaces.size(); il++) {
+ if (plane.subfaces[il].guess_orientation) {
+ bool swap = false; bool noswap = false;
+ for ( auto ie : plane.subfaces[il].edges ) {
+ if (ie<0) continue;
+ for (unsigned int jp=0; jp<planes.size(); jp++) {
+ bool overlap = false; for ( auto ov : planes[ip].overlapping ) if (jp == ov%1000) overlap = true;
+ if (overlap) continue;
+ for (unsigned int jf=0; jf<planes[jp].subfaces.size(); jf++) {
+ if (!planes[jp].subfaces[jf].shared) continue;
+ if (planes[jp].subfaces[jf].guess_orientation) continue;
+ for (auto je : planes[jp].subfaces[jf].edges) {
+ if (je<0) continue;
+ if (jp==ip && jf==il && je==ie) continue;
+ if (planes[jp].plane_edges[je].endpoints.first == plane.plane_edges[ie].endpoints.first &&
+ planes[jp].plane_edges[je].endpoints.second == plane.plane_edges[ie].endpoints.second) { swap = true; break; }
+ if (planes[jp].plane_edges[je].endpoints.first == plane.plane_edges[ie].endpoints.second &&
+ planes[jp].plane_edges[je].endpoints.second == plane.plane_edges[ie].endpoints.first) { noswap = true; break; }
+ }
+ if (swap || noswap) break;
+ }
+ if (swap || noswap) break;
+ }
+ if (swap || noswap) break;
+ }
+ if (noswap) { planes[ip].subfaces[il].guess_orientation = false; }
+ if (swap) { planes[ip].subfaces[il].guess_orientation = false;
+ NPL_ordered loop = planes[ip].subfaces[il];
+ planes[ip].subfaces[il].reverse(planes[ip].plane_edges);
+ /*
+ for (auto iv : planes[ip].subfaces[il].vtx ) std::cout << iv <<","; std::cout <<""<< std::endl;
+ for (auto ie : planes[ip].subfaces[il].edges ) {
+ if (ie<0) std::cout << "dummy edge link:" << std::endl;
+ else std::cout << planes[ip].plane_edges[ie].endpoints.first <<"->" << planes[ip].plane_edges[ie].endpoints.second << std::endl;
+ }
+ */
+ for ( unsigned int jl = 0; jl < plane.subfaces.size(); jl++) {
+ if ( plane.subfaces[jl].shared ) continue;
+ // by construction, the inverted loop sits at the beginning of non-shared loop
+ if (planes[ip].subfaces[jl].vtx.front() == planes[ip].subfaces[jl].vtx[loop.vtx.size()] &&
+ vlog[planes[ip].subfaces[jl].vtx.front()].shared) {
+ if ( planes[ip].subfaces[jl].vtx.front() != loop.vtx.front() ) loop.reorder(planes[ip].subfaces[jl].vtx.front());
+ //std::cout << "before replacement:" << std::endl;
+ //for ( auto iv : planes[ip].subfaces[jl].vtx ) std::cout << iv<<"("<<vlog[iv].shared<<"),"; std::cout <<"" <<std::endl;
+ //for ( auto ie : planes[ip].subfaces[jl].edges ) {
+ // if (ie<0) std::cout << "dummy edge link:" << std::endl;
+ // else std::cout << planes[ip].plane_edges[ie].endpoints.first <<"->" << planes[ip].plane_edges[ie].endpoints.second << std::endl;
+ //}
+ for (unsigned int i=0; i<loop.vtx.size(); i++) {
+ planes[ip].subfaces[jl].vtx[i] = loop.vtx[i];
+ planes[ip].subfaces[jl].edges[i] = loop.edges[i];
+ }
+ }
+ //
+ //for ( auto iv : planes[ip].subfaces[jl].vtx ) std::cout << iv<<"("<<vlog[iv].shared<<"),"; std::cout <<"" <<std::endl;
+ //for ( auto ie : planes[ip].subfaces[jl].edges ) {
+ // if (ie<0) std::cout << "dummy edge link:" << std::endl;
+ // else std::cout << planes[ip].plane_edges[ie].endpoints.first <<"->" << planes[ip].plane_edges[ie].endpoints.second << std::endl;
+ //}
+ }
+ }
+ }
+ }
+ }
- GeoPolyhedronGeneric* gpg = nullptr;
+ //GeoPolyhedronGeneric gpg;
// construct intersection
GeoPolyhedronGeneric intersection;
@@ -722,20 +890,43 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
}
- if (debugMode && intersection.getNFaces()>0 && !intersection.construct(true)) {
+ if (!intersection.construct(debugMode) && debugMode) {
+ //if (debugMode && intersection.getNFaces()>0 && !intersection.construct(true)) {
std::cout <<"VERIFY intersection:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin() << std::endl;
+ std::vector<std::pair<int,int> > missing = missing_edges(intersection);
for (unsigned int ii=0; ii<intersection.getNFaces(); ii++) {
std::vector<size_t> vtx = intersection.getFace(ii);
std::cout <<"face:"<<ii<<":"; for (auto vx : vtx) std::cout <<vx <<","; std::cout <<"" << std::endl;
}
+ std::vector<int> planes_to_check;
+ for (auto ii : missing) { int np = 0;
+ for (auto ip : planes) {
+ for (auto ie : ip.plane_edges) {
+ if (ie.endpoints.first == ii.first && ie.endpoints.second == ii.second ) {
+ std::cout << "missing edge:"<< ii.first <<":"<<ii.second <<":plane:"<< np << ":adj:" << ie.adjacent_face<<":gen.edge:"
+ << ie.generating_edge <<":protected:" <<ie.protectedLine << std::endl;
+ planes_to_check.push_back(np); planes_to_check.push_back(ie.adjacent_face);
+ }
+ }
+ np++;
+ }
+ }
for (unsigned int ip=0; ip<planes.size(); ip++) { NPH_plane plane = planes[ip];
+ bool failed = false; for ( auto jp : planes_to_check) if (ip==jp) failed = true;
+ if (!failed) continue;
+ std::cout << "plane check:" <<ip <<":original edge loop:";
+ for (auto ie : plane.edge_loop_origin) std::cout << ie.endpoints.first << "("<< vlog[ie.endpoints.first].shared<<")"; std::cout <<""<<std::endl;
+ for (auto ie : plane.edge_loop_origin) std::cout << ie.adjacent_face <<","; std::cout <<""<<std::endl;
+ for (auto il : plane.subfaces) { std::cout <<"subface(" << il.shared <<"):overlap:"<<il.top_overlap<<":";
+ for ( auto iv : il.vtx ) std::cout << iv <<","; std::cout <<"" << std::endl;
+ }
+ /*
int nShared = 0; for ( auto il : plane.subfaces ) if (il.shared) nShared++;
if ( plane.hit_vertices.size()>2 && nShared<1 ) std::cout <<"PROBLEM in plane:" << ip << std::endl;
+ */
}
}
- if ( modeBoolean == 0 && intersection.getNFaces() > 3 ) gpg = new GeoPolyhedronGeneric(intersection);
-
// construct subtraction
GeoPolyhedronGeneric subtraction;
subtraction.addVertices(vtx);
@@ -763,7 +954,8 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
}
- if (debugMode && subtraction.getNFaces()>0 && !subtraction.construct(true)) {
+ if (!subtraction.construct(debugMode) && debugMode) {
+ //if (debugMode && subtraction.getNFaces()>0 && !subtraction.construct(true)) {
std::cout <<"VERIFY subtraction:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin() << std::endl;
for (unsigned int ii=0; ii<subtraction.getNFaces(); ii++) {
std::vector<size_t> vtx = subtraction.getFace(ii);
@@ -771,8 +963,6 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
- if ( modeBoolean == 1 && subtraction.getNFaces()>3 ) gpg = new GeoPolyhedronGeneric(subtraction);
-
// construct union : non-shared + half of overlaps
GeoPolyhedronGeneric uni;
uni.addVertices(vtx);
@@ -793,7 +983,8 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
}
- if (debugMode && !uni.construct(true)) {
+ if (!uni.construct(debugMode) && debugMode) {
+ //if (debugMode && !uni.construct(true)) {
std::cout <<"VERIFY union:GPG origin:" << phA->getOrigin() <<":"<<phB->getOrigin()<<":valid? " << uni.isValid() << std::endl;
std::cout << "review input polyhedrons:A:vtx:faces:" << phA->getNVertices() <<":" << phA->getNFaces() << std::endl;
for (unsigned int iv=0 ; iv<phA->getNVertices(); iv++) { GeoTrf::Vector3D xx = phA->getVertex(iv);
@@ -822,7 +1013,6 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
}
}
- if ( modeBoolean == 2 ) gpg = new GeoPolyhedronGeneric(uni);
if (debugMode) {
if (uni.isValid() && intersection.isValid()) std::cout <<"VOLUME check uni:" << phA->volume()+phB->volume()-intersection.volume()
@@ -831,20 +1021,34 @@ GeoPolyhedronGeneric* GeoPolyhedronManipulator::resolveBoolean(const GeoPolyhed
<<"=?" << subtraction.volume() << std::endl;
}
- double v0 = (gpg && gpg->isValid()) ? gpg->volume() : 0.;
+ if ( modeBoolean == 0 ) {
+ intersection.cleanAfterBooleanOperation();
+ return intersection;
+ }
+ if ( modeBoolean == 1 ) {
+ subtraction.cleanAfterBooleanOperation();
+ return subtraction;
+ }
+ if ( modeBoolean == 2 ) {
+ uni.cleanAfterBooleanOperation();
+ return uni;
+ }
+
+ /*
+ double v0 = gpg.isValid() ? gpg.volume() : 0.;
- if (debugMode && gpg) std::cout << "test gpg before cleanup:" << gpg->isValid() <<":volume:" << v0 << std::endl;
+ if (debugMode) std::cout << "test gpg before cleanup:" << gpg.isValid() <<":volume:" << v0 << std::endl;
- if (gpg && gpg->isValid()) gpg->cleanAfterBooleanOperation();
+ if (gpg.isValid()) gpg.cleanAfterBooleanOperation();
if (debugMode) {
- if (gpg) { std::cout << "gpg conversion test:valid:" << gpg->isValid();
- if (gpg->isValid()) std::cout <<":volume:" << gpg->volume() <<":v0:"<<v0;
- std::cout<<""<< std::endl;
- }
+ std::cout << "gpg conversion test:valid:" << gpg.isValid();
+ if (gpg.isValid()) std::cout <<":volume:" << gpg.volume() <<":v0:"<<v0;
+ std::cout<<""<< std::endl;
}
-
- return gpg;
+ */
+
+ return GeoPolyhedronGeneric();
}
std::vector<NPH_edge> GeoPolyhedronManipulator::collectEdges(std::vector<NPH_plane>& planes, std::vector<NPH_vertex>& vlog, bool debugMode ) const {
@@ -903,6 +1107,18 @@ std::vector<NPH_edge> GeoPolyhedronManipulator::collectEdges(std::vector<NPH_pla
planes[iplane].vlog_ori.push_back(plv);
} // end loop vertices
} // end loop planes
+
+ // fourth loop to check presence of adjacent faces
+ for (unsigned int ie=0; ie<edges.size(); ie++) { if (edges[ie].faces.second<0) {
+ std::cout << "edge without adjacent face:" << edges[ie].endpoints.first <<"," << edges[ie].endpoints.second <<":" <<edges[ie].faces.first <<":" << edges[ie].faces.second << std::endl;
+ for (unsigned int je=0; je<edges.size(); je++) { if (je==ie) continue;
+ if (edges[je].endpoints.first== edges[ie].endpoints.first && edges[je].endpoints.second== edges[ie].endpoints.second ) {
+ std::cout << "duplicated edges" << edges[je].endpoints.first <<"," << edges[je].endpoints.second <<":" <<edges[je].faces.first <<":" << edges[je].faces.second << std::endl;
+ edges[ie].faces.second = edges[je].faces.second;
+ }
+ }
+ }
+ }
} // end loop ph
@@ -992,6 +1208,7 @@ void GeoPolyhedronManipulator::collectLines(std::vector<NPH_plane>& planes, std:
}
//
vlog[index].edges.push_back(std::make_pair<int>(ie,ip));
+ //std::cout <<"vlog.edges:" << index <<":"<<ie<<":" << ip <<":" << vlog[index].edges.size() << ":edges:" << edges.size() <<std::endl;
} // valid edge intersection with face
} else {
// special configuration: edge laying on face plane, need to calculate edge intersections
@@ -1088,7 +1305,7 @@ void GeoPolyhedronManipulator::intersectFace(unsigned int in_edge, std::vector<N
else {
GeoTrf::Vector3D hit = eInit + t*edge;
for (unsigned int jv = 0; jv<vtx.size(); jv++) if ( (hit - vtx[jv]).norm() < m_tolerance ) { index = jv ; break ;}
- if (index < 0) { vtx.push_back(hit); index = vtx.size()-1; NPH_vertex pv; pv.index = index; pv.shared = true; vlog.push_back(pv); }
+ if (index < 0) { vtx.push_back(hit); index = vtx.size()-1; NPH_vertex pv; pv.index = index; pv.shared = true; vlog.push_back(pv); }
}
// save on probe edge
NPH_vtx nvtx(index); nvtx.edge=in_edge ; nvtx.hitPlane = edges[in_edge].faces.second ; nvtx.distance = dist; nvtx.faces = std::pair<int,int>(edges[ie].faces);
@@ -1221,7 +1438,7 @@ double GeoPolyhedronManipulator::edgeFaceIntersection( GeoTrf::Vector3D eInit, G
return pen / en ;
}
-std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_index, std::vector<NPH_vertex>& vlog, std::vector<NPH_edge> edges, std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D> vtx) const {
+std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_index, std::vector<NPH_vertex> vlog, std::vector<NPH_edge> edges, std::vector<NPH_plane>& planes, std::vector<GeoTrf::Vector3D> vtx) const {
// loop over vertices to determine the level of ambiguity : results save as NPL_vertex properties
// default mode processes hit vertices ( intersections )
@@ -1230,7 +1447,8 @@ std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_in
std::vector<NPL_edge> lines; std::vector< std::pair< int, std::vector< int > > > active_intersected;
- for (unsigned int ih = 0; ih< plane.hit_vertices.size(); ih++) { NPL_vertex iplv = plane.hit_vertices[ih]; NPH_vertex iv = vlog[iplv.index];
+ for (unsigned int ih = 0; ih< plane.hit_vertices.size(); ih++) {
+ NPL_vertex iplv = plane.hit_vertices[ih]; NPH_vertex iv = vlog[iplv.index];
// overview of adjacent planes
std::vector<int> np; std::vector<int> native_planes; std::vector<int> intersected_planes; int pl; bool found;
for (auto ie : iv.edges) {
@@ -1286,10 +1504,6 @@ std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_in
if ( np.size() == 5) {
std::pair<int,int> doublEdg = hitDoubleEdge(iplv.index, plane_index, plane.index_ph, vlog, edges); // indicates adjacent face
- if (doublEdg.first>=0 && iv.index_ph==plane.index_ph) {
- vlog[iplv.index].hitDoublEdge = doublEdg.first;
- vlog[iplv.index].onDoublEdge = (planes[doublEdg.first].index_ph == plane.index_ph) ? 0 : 1;
- }
if (doublEdg.first>=0 && iv.index_ph!=plane.index_ph) {
planes[plane_index].hit_vertices[ih].faces.first = doublEdg.second;
// check if valid intersecting plane
@@ -1299,8 +1513,6 @@ std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_in
}
planes[plane_index].hit_vertices[ih].faces.second = validPlane ? doublEdg.first : doublEdg.second;
if (!validPlane) iplv.onEdgeDummy = true;
- vlog[iplv.index].hitDoublEdge = validPlane ? doublEdg.first : doublEdg.second ;
- vlog[iplv.index].onDoublEdge = validPlane ? 1 : 0;
}
}
}
@@ -1353,6 +1565,14 @@ std::vector<NPL_edge> GeoPolyhedronManipulator::evaluateVertexLogic(int plane_in
// create line edge
for (unsigned int jl =0; jl<vxord.size()-1; jl++) {
+ GeoTrf::Vector3D vlmid = 0.5 * (vtx[vxord[jl].first] + vtx[vxord[jl+1].first]);
+ //std::cout << "mid-line check:" << vxord[jl].first <<","<<vxord[jl+1].first <<":inside face:" << insideFace(plane_index,vlmid)
+ // << ":phA contains:"<< m_phA->contains(vlmid) <<":phB contains:" << m_phB->contains(vlmid) << std::endl;
+ if (!insideFace(plane_index,vlmid) || !m_phA->contains(vlmid) || !m_phB->contains(vlmid) ) continue;
+ bool foundLine = false;
+ for ( auto il : lines ) if ( vxord[jl].first == il.endpoints.first && vxord[jl+1].first == il.endpoints.second) foundLine = true;
+ for ( auto il : lines ) if ( vxord[jl].first == il.endpoints.second && vxord[jl+1].first == il.endpoints.first) foundLine = true;
+ if (foundLine) continue; // avoid duplicates
lines.push_back(createLine(vxord[jl].first,vxord[jl+1].first,ii.first, plane) );
if (!splitLine) lines.back().protectedLine = true;
}
@@ -1480,8 +1700,6 @@ void GeoPolyhedronManipulator::saveHitOnEdge(NPH_vtx nvtx, std::vector<NPH_edge
else break;
}
edges[nvtx.edge].evtx.insert(vIt,nvtx);
- //std::cout <<"inserting evtx(method):" << nvtx.index <<":at:" << nvtx.distance <<":endpoints:" <<
- // edges[nvtx.edge].endpoints.first << "," << edges[nvtx.edge].endpoints.second<<":hit plane:" << nvtx.hitPlane<< std::endl;
}
NPL_edge GeoPolyhedronManipulator::createLine(int iv1, int iv2, int adj, NPH_plane plane) const {
@@ -1545,8 +1763,12 @@ NPL_ordered GeoPolyhedronManipulator::mergeLoops( NPL_ordered loop1, NPL_ordered
merged.vtx.push_back(merged.vtx.front());
loop2.vtx.push_back(loop2.vtx.front());
+ merged.edges.push_back(-1);
+ loop2.edges.push_back(-1);
merged.vtx.insert(merged.vtx.begin(),loop2.vtx.begin(),loop2.vtx.end());
+ merged.edges.insert(merged.edges.begin(),loop2.edges.begin(),loop2.edges.end());
+
}
return merged;
@@ -1562,3 +1784,189 @@ GeoTrf::Vector3D GeoPolyhedronManipulator::areaVec( std::vector<size_t> loop, st
return area;
}
+
+int GeoPolyhedronManipulator::nextOri(int input, std::vector<NPL_edge> plane_edges,
+ std::vector<NPL_edge> edge_loop_origin, std::vector<int> availableEdges, bool previous) const {
+
+ int next = -1;
+
+ // find input edge in the original loop
+ int ori = -1;
+ for ( unsigned int ieo = 0; ieo < edge_loop_origin.size(); ieo++) {
+ if (plane_edges[input].endpoints.first == edge_loop_origin[ieo].endpoints.first &&
+ plane_edges[input].endpoints.second == edge_loop_origin[ieo].endpoints.second &&
+ plane_edges[input].adjacent_face == edge_loop_origin[ieo].adjacent_face ) { ori = ieo; break; }
+ }
+ // next edge in original set
+ int ori_next = ori<edge_loop_origin.size()-1 ? ori+1 : 0;
+ if (previous) ori_next = ori>0 ? ori-1 : edge_loop_origin.size()-1;
+
+ // find index in full set
+ for ( unsigned int ie = 0; ie < plane_edges.size(); ie++) {
+ if (plane_edges[ie].endpoints.first == edge_loop_origin[ori_next].endpoints.first &&
+ plane_edges[ie].endpoints.second == edge_loop_origin[ori_next].endpoints.second &&
+ plane_edges[ie].adjacent_face == edge_loop_origin[ori_next].adjacent_face ) { next = ori_next; break; }
+ }
+ // check is available
+ bool available = false;
+ for ( auto ie : availableEdges) if (next == ie) available = true;
+
+ if (available) return next;
+ else return -1;
+
+}
+
+std::vector<std::pair<int,int> > GeoPolyhedronManipulator::missing_edges( const GeoPolyhedronGeneric gpg) const {
+
+ std::vector<std::pair<int,int> > missing;
+
+ // loop over edges, check that each appears with its mirrored version ( attached faces )
+ for ( unsigned int ic=0; ic < gpg.getNFaces(); ic++) {
+ std::vector<size_t> vtx = gpg.getFace(ic);
+ if (vtx.size()<2) continue; // should not happen though
+ for (unsigned int iv=0; iv<vtx.size(); iv++) {
+ std::pair<size_t,size_t> edge(vtx[iv], iv<vtx.size()-1 ? vtx[iv+1] : vtx[0]); // probe
+ bool mirrored = false; // search mirrored edge
+ for ( unsigned int icm=0; icm < gpg.getNFaces(); icm++) {
+ std::vector<size_t> vtm = gpg.getFace(icm);
+ if (vtm.size()<2) continue; // should not happen though
+ for (unsigned int ivm=0; ivm<vtm.size(); ivm++) {
+ std::pair<size_t,size_t> edgm(vtm[ivm], ivm<vtm.size()-1 ? vtm[ivm+1] : vtm[0]); // test
+ if ( edgm.first == edge.second && edgm.second == edge.first) mirrored = true;
+ if (mirrored) break;
+ } // end loop over other edges
+ if (mirrored) break;
+ } // end loop over other faces
+ if (!mirrored) {
+ missing.push_back(std::pair<int,int>(edge.second,edge.first));
+ }
+ } // end loop over edges
+ } // end loop over faces
+
+ return missing;
+}
+
+GeoPolyhedronGeneric GeoPolyhedronManipulator::resolve(const GeoShape* shape, GeoTrf::Transform3D& trf, int polyhedrize, int level, bool debugMode) {
+
+ if (debugMode) std::cout << "recursively resolve:" << shape->type() << std::endl;
+
+ if (shape->type() == "Shift") {
+ const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*> (shape);
+ if (shift) {
+ GeoTrf::Transform3D trshift = trf*shift->getX();
+ GeoPolyhedronGeneric gp = resolve(shift->getOp(), trshift, polyhedrize, level+1 );
+ if (debugMode) {
+ std::cout << "resolved?:"<< gp.getStatus() <<":level:"<<level<<":" << shape->type() <<":number of faces:" << gp.getNFaces();
+ if (gp.getStatus()==0) std::cout<<":vol:"<< gp.volume()<< std::endl;
+ else std::cout << ""<< std::endl;
+ }
+ return gp;
+ }
+ }
+ else if (shape->type() == "Union") {
+ const GeoShapeUnion* uni = dynamic_cast<const GeoShapeUnion*> (shape);
+ GeoPolyhedronGeneric gpA = resolve(uni->getOpA(), trf, polyhedrize, level+1 );
+ GeoPolyhedronGeneric gpB = resolve(uni->getOpB(), trf, polyhedrize, level+1 );
+ if (gpA.getStatus() > 0 || gpB.getStatus() > 0 ) return ( gpA.getStatus() > gpB.getStatus() ? gpA : gpB );
+ if (debugMode) std::cout << "entering resolveBoolean for union:faces:" << gpA.getNFaces() <<"," << gpB.getNFaces() <<":status:" << gpA.getStatus()<<"," << gpB.getStatus()<< std::endl;
+ GeoPolyhedronGeneric gpp = (gpA.getNFaces()==0 || gpB.getNFaces()==0) ? (gpA.getNFaces()==0 ? gpB : gpA) : resolveBoolean(&gpA,&gpB,2,false);
+ if (debugMode) {
+ std::cout << "resolved ?:"<< gpp.getStatus()<<":level:" <<level<<":" << shape->type() <<":number of faces:" << gpp.getNFaces();
+ if (gpp.getStatus()==0) std::cout<<":vol:"<< gpp.volume()<< std::endl;
+ else std::cout << ""<< std::endl;
+ }
+ return gpp;
+ }
+ else if (shape->type() == "Subtraction") {
+ const GeoShapeSubtraction* sub = dynamic_cast<const GeoShapeSubtraction*> (shape);
+ GeoPolyhedronGeneric gpA = resolve(sub->getOpA(), trf, polyhedrize, level+1 );
+ GeoPolyhedronGeneric gpB = resolve(sub->getOpB(), trf, polyhedrize, level+1 );
+ if (gpA.getStatus() > 0 || gpB.getStatus() > 0 ) return ( gpA.getStatus() > gpB.getStatus() ? gpA : gpB );
+ if (debugMode) std::cout << "entering resolveBoolean for subtraction:faces:" << gpA.getNFaces() <<"," << gpB.getNFaces() <<":status:" << gpA.getStatus()<<"," << gpB.getStatus()<< std::endl;
+ GeoPolyhedronGeneric gpp = (gpA.getNFaces()==0 || gpB.getNFaces()==0) ? gpA : resolveBoolean(&gpA,&gpB,1,false);
+ if (debugMode) {
+ std::cout << "resolved?:"<< gpp.getStatus()<<":level:" <<level<<":" << shape->type() <<":number of faces:" << gpp.getNFaces();
+ if (gpp.getStatus()==0) std::cout<<":vol:"<< gpp.volume()<< std::endl;
+ else std::cout << ""<< std::endl;
+ }
+ return gpp;
+ }
+ else if (shape->type() == "Intersection") {
+ const GeoShapeIntersection* ix = dynamic_cast<const GeoShapeIntersection*> (shape);
+ GeoPolyhedronGeneric gpA = resolve(ix->getOpA(), trf, polyhedrize, level+1 );
+ GeoPolyhedronGeneric gpB = resolve(ix->getOpB(), trf, polyhedrize, level+1 );
+ if (gpA.getStatus() > 0 || gpB.getStatus() > 0 ) return ( gpA.getStatus() > gpB.getStatus() ? gpA : gpB );
+ if (debugMode) std::cout << "entering resolveBoolean for intersection:faces:" << gpA.getNFaces() <<"," << gpB.getNFaces() <<":status:" << gpA.getStatus()<<"," << gpB.getStatus()<< std::endl;
+ GeoPolyhedronGeneric gpp = (gpA.getNFaces()==0 || gpB.getNFaces()==0) ? (gpA.getNFaces()==0 ? gpA : gpB) : resolveBoolean(&gpA,&gpB,0,false);
+ if (debugMode) {
+ std::cout << "resolved?:"<< gpp.getStatus()<<":level:" <<level<<":" << shape->type() <<":number of faces:" << gpp.getNFaces();
+ if (gpp.getStatus()==0) std::cout<<":vol:"<< gpp.volume()<< std::endl;
+ else std::cout << ""<< std::endl;
+ }
+ return gpp;
+ } else {
+ GeoPolyhedronGeneric gp(shape,trf,polyhedrize);
+ if (debugMode) {
+ if (shape->type() == "PolyhedronGeneric") std::cout << "gpg on input:" << gp.isValid() << std::endl;
+ if (gp.isValid()) std::cout << "resolved:"<<level<<":" << shape->type() <<":number of faces:" << gp.getNFaces()<<":vol:"<<gp.volume() << std::endl;
+ }
+ if (gp.isValid()) return gp;
+ else if (debugMode) std::cout << "shape conversion ->gpg fails:" << shape->type() << std::endl;
+ }
+
+ GeoPolyhedronGeneric gpu; gpu.setStatus(1);
+ if (debugMode) std::cout << "unresolved:" <<level<<":"<<shape->type()<< std::endl;
+ return gpu;
+
+}
+
+void GeoPolyhedronManipulator::collectFaceNumber(const GeoShape* shape, int& nFaces ) {
+
+ if ( shape->type() == "Box" || shape->type() == "Trd") nFaces+=6;
+ else if ( shape->type() == "SimplePolygonBrep") {
+ const GeoSimplePolygonBrep* spb = dynamic_cast<const GeoSimplePolygonBrep*> (shape);
+ nFaces += 2+spb->getNVertices();
+ }
+ else if (shape->type() == "Shift") {
+ const GeoShapeShift* ix = dynamic_cast<const GeoShapeShift*> (shape);
+ collectFaceNumber( ix->getOp(), nFaces );
+ }
+ else if (shape->type() == "Intersection") {
+ const GeoShapeIntersection* ix = dynamic_cast<const GeoShapeIntersection*> (shape);
+ collectFaceNumber( ix->getOpA(), nFaces );
+ collectFaceNumber( ix->getOpB(), nFaces );
+ }
+ else if (shape->type() == "Subtraction") {
+ const GeoShapeSubtraction* ix = dynamic_cast<const GeoShapeSubtraction*> (shape);
+ collectFaceNumber( ix->getOpA(), nFaces );
+ collectFaceNumber( ix->getOpB(), nFaces );
+ }
+ else if (shape->type() == "Union") {
+ const GeoShapeUnion* ix = dynamic_cast<const GeoShapeUnion*> (shape);
+ collectFaceNumber( ix->getOpA(), nFaces );
+ collectFaceNumber( ix->getOpB(), nFaces );
+ }
+ else if (shape->type() == "Tube") {
+ const GeoTube* tub = dynamic_cast<const GeoTube*> (shape);
+ nFaces +=3;
+ if (tub->getRMin() > 0 ) nFaces +=1;
+ //if (tub->getDPhi() < acos(-1.) ) nFaces +=2;
+ }
+ else if (shape->type() == "Tubs") {
+ const GeoTubs* tubs = dynamic_cast<const GeoTubs*> (shape);
+ nFaces +=3;
+ if (tubs->getRMin() > 0 ) nFaces +=1;
+ if (tubs->getDPhi() < acos(-1.) ) nFaces +=2;
+ }
+ else if (shape->type() == "Pcon") {
+ const GeoPcon* pcon = dynamic_cast<const GeoPcon*> (shape);
+ int np = pcon->getNPlanes()-1;
+ nFaces +=2 + np;
+ if (pcon->getRMinPlane(0) > 0 ) nFaces += np;
+ if (pcon->getDPhi() < acos(-1.) ) nFaces += 2*np ;
+ }
+ else std::cout << "faces not collected:" << shape->type() << std::endl;
+
+
+
+}