From b7e8d9794a3319b8816ae34ceebdb9617cbd28eb Mon Sep 17 00:00:00 2001 From: vbd Date: Fri, 18 Apr 2008 08:24:32 +0000 Subject: [PATCH] adding Anthony modification for integration in V 5.0 --- .../Geometric2D/AbstractEdge.hxx | 3 +- src/INTERP_KERNEL/Geometric2D/Bounds.cxx | 16 ++++ src/INTERP_KERNEL/Geometric2D/Bounds.hxx | 4 + .../Geometric2D/ComposedEdge.cxx | 20 +++- .../Geometric2D/ComposedEdge.hxx | 3 +- src/INTERP_KERNEL/Geometric2D/Edge.hxx | 3 +- .../Geometric2D/EdgeArcCircle.cxx | 95 ++++++++++++++++--- .../Geometric2D/EdgeArcCircle.hxx | 4 +- src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx | 31 ++++-- src/INTERP_KERNEL/Geometric2D/EdgeLin.hxx | 4 +- .../Geometric2D/ElementaryEdge.cxx | 10 +- .../Geometric2D/ElementaryEdge.hxx | 3 +- src/INTERP_KERNEL/Geometric2D/Makefile.am | 1 + src/INTERP_KERNEL/Geometric2D/Node.cxx | 8 +- src/INTERP_KERNEL/Geometric2D/Node.hxx | 6 +- src/INTERP_KERNEL/Geometric2D/Precision.cxx | 15 +++ src/INTERP_KERNEL/Geometric2D/Precision.hxx | 10 +- .../Geometric2D/QuadraticPolygon.cxx | 69 ++++++++++---- .../Geometric2D/QuadraticPolygon.hxx | 5 +- 19 files changed, 243 insertions(+), 67 deletions(-) create mode 100644 src/INTERP_KERNEL/Geometric2D/Precision.cxx diff --git a/src/INTERP_KERNEL/Geometric2D/AbstractEdge.hxx b/src/INTERP_KERNEL/Geometric2D/AbstractEdge.hxx index a30eaf23b..ec2e21e5b 100644 --- a/src/INTERP_KERNEL/Geometric2D/AbstractEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/AbstractEdge.hxx @@ -57,6 +57,7 @@ namespace INTERP_KERNEL virtual double getAreaOfZone() const = 0; virtual void fillBounds(Bounds& output) const = 0; virtual void getAllNodes(std::set& output) const = 0; + virtual void getBarycenter(double *bary, double& weigh) const = 0; virtual ElementaryEdge* &getLastElementary(IteratorOnComposedEdge::ItOnFixdLev &delta) = 0; virtual ElementaryEdge* &getFirstElementary(IteratorOnComposedEdge::ItOnFixdLev &delta) = 0; virtual const AbstractEdge *&operator[](IteratorOnComposedEdge::ItOnFixdLev i) const = 0; @@ -65,9 +66,9 @@ namespace INTERP_KERNEL virtual Node *getStartNode() const = 0; virtual bool changeStartNodeWith(Node *node) const = 0; virtual bool changeEndNodeWith(Node *node) const = 0; - virtual void dumpInXfigFile(std::ostream& stream) const = 0; virtual bool intresicEqual(const AbstractEdge *other) const = 0; virtual bool intresicEqualDirSensitive(const AbstractEdge *other) const = 0; + virtual void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const = 0; virtual bool intresincEqCoarse(const Edge *other) const = 0; virtual bool getDirection() const = 0; }; diff --git a/src/INTERP_KERNEL/Geometric2D/Bounds.cxx b/src/INTERP_KERNEL/Geometric2D/Bounds.cxx index 74bfc7ab6..bd6297221 100644 --- a/src/INTERP_KERNEL/Geometric2D/Bounds.cxx +++ b/src/INTERP_KERNEL/Geometric2D/Bounds.cxx @@ -43,6 +43,22 @@ void Bounds::prepareForAggregation() _xMin=1e200; _xMax=-1e200; _yMin=1e200; _yMax=-1e200; } +double Bounds::fitXForXFigD(double val, int res) const +{ + double delta=fmax(_xMax-_xMin,_yMax-_yMin)/2.; + double ret=val-(_xMax+_xMin)/2.+delta; + delta=11.1375*res/(2.*delta); + return ret*delta; +} + +double Bounds::fitYForXFigD(double val, int res) const +{ + double delta=fmax(_xMax-_xMin,_yMax-_yMin)/2.; + double ret=val-(_yMax+_yMin)/2.+delta; + delta=11.1375*res/(2.*delta); + return ret*delta; +} + Bounds *Bounds::nearlyAmIIntersectingWith(const Bounds& other) const { if( (other._xMin > _xMax+QUADRATIC_PLANAR::_precision) || (other._xMax < _xMin-QUADRATIC_PLANAR::_precision) || (other._yMin > _yMax+QUADRATIC_PLANAR::_precision) diff --git a/src/INTERP_KERNEL/Geometric2D/Bounds.hxx b/src/INTERP_KERNEL/Geometric2D/Bounds.hxx index 89ccf8242..6b3553786 100644 --- a/src/INTERP_KERNEL/Geometric2D/Bounds.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Bounds.hxx @@ -24,6 +24,10 @@ namespace INTERP_KERNEL Bounds(double xMin, double xMax, double yMin, double yMax):_xMin(xMin),_xMax(xMax),_yMin(yMin),_yMax(yMax) { } void setValues(double xMin, double xMax, double yMin, double yMax) { _xMin=xMin; _xMax=xMax; _yMin=yMin; _yMax=yMax; } void prepareForAggregation(); + int fitXForXFig(double val, int res) const { return (int)fitXForXFigD(val,res); } + int fitYForXFig(double val, int res) const { return (int)fitYForXFigD(val,res); } + double fitXForXFigD(double val, int res) const; + double fitYForXFigD(double val, int res) const; Bounds *nearlyAmIIntersectingWith(const Bounds& other) const; Bounds *amIIntersectingWith(const Bounds& other) const; //! No approximations. diff --git a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx index ed0f94c3f..5d7f1eebd 100644 --- a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx @@ -103,10 +103,11 @@ double ComposedEdge::getAreaOfZone() const return ret; } -void ComposedEdge::dumpInXfigFile(std::ostream& stream) const +void ComposedEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const { + stream.precision(10); for(std::vector::const_iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++) - (*iter)->dumpInXfigFile(stream); + (*iter)->dumpInXfigFile(stream,resolution,box); } Node *ComposedEdge::getEndNode() const @@ -186,6 +187,21 @@ void ComposedEdge::getAllNodes(std::set& output) const (*iter)->getAllNodes(output); } +void ComposedEdge::getBarycenter(double *bary, double& weigh) const +{ + weigh=0.; bary[0]=0.; bary[1]=0.; + double tmp1,tmp2[2]; + for(vector::const_iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++) + { + (*iter)->getBarycenter(tmp2,tmp1); + weigh+=tmp1; + bary[0]+=tmp1*tmp2[0]; + bary[1]+=tmp1*tmp2[1]; + } + bary[0]/=weigh; + bary[1]/=weigh; +} + bool ComposedEdge::isInOrOut(Node *nodeToTest) const { Bounds b; b.prepareForAggregation(); diff --git a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx index 60c93e1e9..16a24adf1 100644 --- a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx @@ -21,6 +21,7 @@ namespace INTERP_KERNEL double getAreaOfZone() const; void fillBounds(Bounds& output) const; void getAllNodes(std::set& output) const; + void getBarycenter(double *bary, double& weigh) const; bool completed() const { return getEndNode()==getStartNode(); } ElementaryEdge * &getLastElementary(IteratorOnComposedEdge::ItOnFixdLev &delta); ElementaryEdge * &getFirstElementary(IteratorOnComposedEdge::ItOnFixdLev &delta); @@ -42,7 +43,7 @@ namespace INTERP_KERNEL bool addEdgeIfIn(ElementaryEdge *edge); bool changeEndNodeWith(Node *node) const; bool changeStartNodeWith(Node *node) const; - void dumpInXfigFile(std::ostream& stream) const; + void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; bool intresicEqual(const AbstractEdge *other) const; bool intresicEqualDirSensitive(const AbstractEdge *other) const; bool isInOrOut(Node *nodeToTest) const; diff --git a/src/INTERP_KERNEL/Geometric2D/Edge.hxx b/src/INTERP_KERNEL/Geometric2D/Edge.hxx index 608fed38a..2a72b6c39 100644 --- a/src/INTERP_KERNEL/Geometric2D/Edge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Edge.hxx @@ -189,6 +189,7 @@ namespace INTERP_KERNEL //! returns area between this and axe Ox delimited along Ox by _start and _end. virtual double getAreaOfZone() const = 0; virtual double getCurveLength() const = 0; + virtual void getBarycenter(double *bary) const = 0; //! Retrieves a point that is owning to this, well placed for IN/OUT detection of this. Typically midlle of this is returned. virtual Node *buildRepresentantOfMySelf() const = 0; //! Given a magnitude specified by sub-type returns if in or not. See getCharactValue method. @@ -203,7 +204,7 @@ namespace INTERP_KERNEL const EdgeArcCircle * &arcSeg) const = 0; bool intersectWith(const Edge *other, MergePoints& commonNode, ComposedEdge& outVal1, ComposedEdge& outVal2) const; - virtual void dumpInXfigFile(std::ostream& stream, bool direction) const = 0; + virtual void dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const = 0; protected: Edge():_cnt(1),_loc(FULL_UNKNOWN),_start(0),_end(0) { } virtual ~Edge(); diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx index bd8ad250e..9b91e860e 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx @@ -122,8 +122,71 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi double u[2];//u is normalized vector from center1 to center2. u[0]=(center2[0]-center1[0])/_dist; u[1]=(center2[1]-center1[1])/_dist; double angle0_1=EdgeArcCircle::safeAcos(u[0]); - double angle0_2; + if(u[1]<0.) + angle0_1=-angle0_1; + double d1_1y=EdgeArcCircle::safeSqrt(radius1*radius1-d1_1*d1_1); + double angleE1=normalizeAngle(getE1().getAngle0()+getE1().getAngle()); + double angleE2=normalizeAngle(getE2().getAngle0()+getE2().getAngle()); + if(!Node::areDoubleEquals(d1_1y,0)) + { + //2 intersections + double v1[2],v2[2]; + v1[0]=u[0]*d1_1-u[1]*d1_1y; v1[1]=u[1]*d1_1+u[0]*d1_1y; + v2[0]=u[0]*d1_1+u[1]*d1_1y; v2[1]=u[1]*d1_1-u[0]*d1_1y; + Node *node1=new Node(center1[0]+v1[0],center1[1]+v1[1]); node1->declareOn(); + Node *node2=new Node(center1[0]+v2[0],center1[1]+v2[1]); node2->declareOn(); + double angle1_1=EdgeArcCircle::safeAcos(v1[0]/radius1); + if(v1[1]<0.) + angle1_1=-angle1_1; + double angle2_1=EdgeArcCircle::safeAcos(v2[0]/radius1); + if(v2[1]<0.) + angle2_1=-angle2_1; + double v3[2],v4[2]; + v3[0]=center1[0]-center2[0]+v1[0]; v3[1]=center1[1]-center2[1]+v1[1]; + v4[0]=center1[0]-center2[0]+v2[0]; v4[1]=center1[1]-center2[1]+v2[1]; + double angle1_2=EdgeArcCircle::safeAcos(v3[0]/radius2); + if(v3[1]<0.) + angle1_2=-angle1_2; + double angle2_2=EdgeArcCircle::safeAcos(v4[0]/radius2); + if(v4[1]<0.) + angle2_2=-angle2_2; + // + bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1); + bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1); + bool e1_2S=Node::areDoubleEqualsWP(angle1_2,getE2().getAngle0(),radius1); + bool e1_2E=Node::areDoubleEqualsWP(angle1_2,angleE2,radius1); + // + bool e2_1S=Node::areDoubleEqualsWP(angle2_1,getE1().getAngle0(),radius2); + bool e2_1E=Node::areDoubleEqualsWP(angle2_1,angleE1,radius2); + bool e2_2S=Node::areDoubleEqualsWP(angle2_2,getE2().getAngle0(),radius2); + bool e2_2E=Node::areDoubleEqualsWP(angle2_2,angleE2,radius2); + ret.push_back(IntersectElement(angle1_1,angle1_2,e1_1S,e1_1E,e1_2S,e1_2E,node1,_e1,_e2,keepOrder())); + ret.push_back(IntersectElement(angle2_1,angle2_2,e2_1S,e2_1E,e2_2S,e2_2E,node2,_e1,_e2,keepOrder())); + } + else + { + //tangent intersection + double v1[2],v2[2]; + v1[0]=d1_1*u[0]; v1[1]=d1_1*u[1]; + v2[0]=center1[0]-center2[0]+v1[0]; v2[1]=center1[1]-center2[1]+v1[1]; + double angle0_1=EdgeArcCircle::safeAcos(v1[0]/radius1); + if(v1[1]<0.) + angle0_1=-angle0_1; + double angle0_2=EdgeArcCircle::safeAcos(v2[0]/radius2); + if(v2[1]<0.) + angle0_2=-angle0_2; + bool e0_1S=Node::areDoubleEqualsWP(angle0_1,getE1().getAngle0(),radius1); + bool e0_1E=Node::areDoubleEqualsWP(angle0_1,angleE1,radius1); + bool e0_2S=Node::areDoubleEqualsWP(angle0_2,getE2().getAngle0(),radius2); + bool e0_2E=Node::areDoubleEqualsWP(angle0_2,angleE2,radius2); + Node *node=new Node(center1[0]+d1_1*u[0],center1[1]+d1_1*u[1]); node->declareOnTangent(); + ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder())); + } + return ret; +} + /*double angle0_2; double signDeltaAngle2; + double d1_2; if(u[1]<0.) angle0_1=-angle0_1; if(d1_1>=0.) @@ -151,8 +214,7 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi double angleE2=normalizeAngle(getE2().getAngle0()+getE2().getAngle()); if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) ) { - //2 intersections - double d1_2=_dist-d1_1; + //2 intersections double deltaAngle1=EdgeArcCircle::safeAcos(fabs(d1_1)/radius1); //owns to 0;Pi/2 by construction double deltaAngle2=EdgeArcCircle::safeAcos(fabs(d1_2)/radius2); //owns to 0;Pi/2 by construction double angle1_1=normalizeAngle(angle0_1+deltaAngle1);// Intersection 1 seen for _e1 @@ -184,8 +246,7 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi Node *node=new Node(center1[0]+radius1*cos(angle0_1),center1[0]+radius1*sin(angle0_1)); node->declareOnTangent(); ret.push_back(IntersectElement(angle0_1,angle0_2,e0_1S,e0_1E,e0_2S,e0_2E,node,_e1,_e2,keepOrder())); } - return ret; -} + return ret;*/ /*! * Idem isAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[. @@ -257,8 +318,7 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic const double *center=getE1().getCenter(); if(!(fabs(_determinant)<(QUADRATIC_PLANAR::_precision*10.)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx)) { - double determinant=fmax(_determinant,0.); - determinant=sqrt(_determinant); + double determinant=EdgeArcCircle::safeSqrt(_determinant); double x1=(_cross*_dy+Node::sign(_dy)*_dx*determinant)/_drSq+center[0]; double y1=(-_cross*_dx+fabs(_dy)*determinant)/_drSq+center[1]; Node *intersect1=new Node(x1,y1); intersect1->declareOn(); @@ -327,6 +387,7 @@ EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, doubl { _center[0]=center[0]; _center[1]=center[1]; + updateBounds(); } void EdgeArcCircle::changeMiddle(Node *newMiddle) @@ -342,7 +403,7 @@ Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) double ex=((*end)[0]-_center[0])/_radius; double ey=((*end)[1]-_center[1])/_radius; double angle0=safeAcos(direction?sx:ex); - angle0=direction?sy:ey>0.?angle0:-angle0; + angle0=(direction?sy:ey)>0.?angle0:-angle0; double deltaAngle=safeAcos(sx*ex+sy*ey); deltaAngle=sx*ey-sy*ex>0.?deltaAngle:-deltaAngle; if(deltaAngle>0. && _angle<0.) @@ -361,7 +422,7 @@ void EdgeArcCircle::getArcOfCirclePassingThru(const double *start, const double double b2=(end[1]*end[1]+end[0]*end[0]-middle[0]*middle[0]-middle[1]*middle[1])/2; center[0]=((end[1]-middle[1])*b1+(start[1]-middle[1])*b2)/delta; center[1]=((middle[0]-end[0])*b1+(middle[0]-start[0])*b2)/delta; - radius=sqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1])); + radius=safeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1])); angleInRad0=safeAcos((start[0]-center[0])/radius); if((start[1]-center[1])<0.) angleInRad0=-angleInRad0; @@ -375,7 +436,7 @@ void EdgeArcCircle::getArcOfCirclePassingThru(const double *start, const double angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI; } -void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction) const +void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const { stream << "5 1 0 1 "; fillXfigStreamForLoc(stream); @@ -385,12 +446,12 @@ void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction) const else stream << '1';//'1' stream << " 0 0 "; - stream << (int)(_center[0]*Node::CST_FOR_XFIG) << " " << (int)(_center[1]*Node::CST_FOR_XFIG) << " "; - direction?_start->dumpInXfigFile(stream):_end->dumpInXfigFile(stream); + stream << box.fitXForXFigD(_center[0],resolution) << " " << box.fitYForXFigD(_center[1],resolution) << " "; + direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box); Node *middle=buildRepresentantOfMySelf(); - middle->dumpInXfigFile(stream); + middle->dumpInXfigFile(stream,resolution,box); middle->decrRef(); - direction?_end->dumpInXfigFile(stream):_start->dumpInXfigFile(stream); + direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box); stream << endl; } @@ -410,6 +471,12 @@ double EdgeArcCircle::getCurveLength() const return fabs(_angle*_radius); } +void EdgeArcCircle::getBarycenter(double *bary) const +{ + bary[0]=((*_start)[0]+(*_end)[0])/2.; + bary[1]=((*_start)[1]+(*_end)[1])/2.; +} + /*! * Characteristic value used is angle in ]_Pi;Pi[ from axe 0x. */ diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx index d787a6f7e..d8868d231 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx @@ -57,10 +57,11 @@ namespace INTERP_KERNEL EdgeArcCircle(Node *start, Node *end, const double *center, double radius, double angle0, double deltaAngle, bool direction=true); //! for tests void changeMiddle(Node *newMiddle); - void dumpInXfigFile(std::ostream& stream, bool direction) const; + void dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const; void update(Node *m); double getAreaOfZone() const; double getCurveLength() const; + void getBarycenter(double *bary) const; bool isIn(double characterVal) const; Node *buildRepresentantOfMySelf() const; bool isLower(double val1, double val2) const; @@ -78,6 +79,7 @@ namespace INTERP_KERNEL static void getArcOfCirclePassingThru(const double *start, const double *middle, const double *end, double *center, double& radius, double& angleInRad, double& angleInRad0); //! To avoid in aggressive optimizations nan. + static double safeSqrt(double val) { double ret=fmax(val,0.); return sqrt(ret); } static double safeAcos(double cosAngle) { double ret=fmin(cosAngle,1.); ret=fmax(ret,-1.); return acos(ret); } protected: void updateBounds(); diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx b/src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx index 749eb08e1..fd4b86178 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeLin.cxx @@ -1,5 +1,6 @@ #include "EdgeLin.hxx" #include "Node.hxx" +#include "InterpolationUtils.hxx" using namespace std; using namespace INTERP_KERNEL; @@ -78,6 +79,17 @@ std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicV return ret; } +/*! + * retrieves if segs are colinears. + * WARNING !!! Contrary to areOverlappedOrOnlyColinears method, this method use an + * another precision to detect colinearity ! + */ +bool SegSegIntersector::areColinears() const +{ + double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2]; + return fabs(determinant)dumpInXfigFile(stream):_end->dumpInXfigFile(stream); - direction?_end->dumpInXfigFile(stream):_start->dumpInXfigFile(stream); + direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box); + direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box); stream << endl; } @@ -182,6 +187,12 @@ double EdgeLin::getAreaOfZone() const return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.; } +void EdgeLin::getBarycenter(double *bary) const +{ + bary[0]=((*_start)[0]+(*_end)[0])/2.; + bary[1]=((*_start)[1]+(*_end)[1])/2.; +} + double EdgeLin::getCurveLength() const { double x=(*_start)[0]-(*_end)[0]; diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeLin.hxx b/src/INTERP_KERNEL/Geometric2D/EdgeLin.hxx index ec8d5f101..23cde4982 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeLin.hxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeLin.hxx @@ -10,6 +10,7 @@ namespace INTERP_KERNEL friend class Edge; public: SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2); + bool areColinears() const; bool haveTheySameDirection() const; void getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const; void areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped); @@ -32,11 +33,12 @@ namespace INTERP_KERNEL EdgeLin(double sX, double sY, double eX, double eY); ~EdgeLin(); TypeOfFunction getTypeOfFunc() const { return SEG; } - void dumpInXfigFile(std::ostream& stream, bool direction) const; + void dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const; void update(Node *m); double getNormSq() const; double getAreaOfZone() const; double getCurveLength() const; + void getBarycenter(double *bary) const; bool isIn(double characterVal) const; Node *buildRepresentantOfMySelf() const; double getCharactValue(const Node& node) const; diff --git a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.cxx b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.cxx index 797feee48..be84a880b 100644 --- a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.cxx @@ -30,6 +30,12 @@ void ElementaryEdge::getAllNodes(std::set& output) const output.insert(_ptr->getEndNode()); } +void ElementaryEdge::getBarycenter(double *bary, double& weigh) const +{ + _ptr->getBarycenter(bary); + weigh=_ptr->getCurveLength(); +} + AbstractEdge *ElementaryEdge::clone() const { return new ElementaryEdge(*this); @@ -145,9 +151,9 @@ bool ElementaryEdge::changeStartNodeWith(Node *node) const return _ptr->changeEndNodeWith(node); } -void ElementaryEdge::dumpInXfigFile(std::ostream& stream) const +void ElementaryEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const { - _ptr->dumpInXfigFile(stream,_direction); + _ptr->dumpInXfigFile(stream,_direction,resolution,box); } bool ElementaryEdge::intresicEqual(const AbstractEdge *other) const diff --git a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx index 364d22f16..ae2dab936 100644 --- a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx @@ -26,6 +26,7 @@ namespace INTERP_KERNEL double getAreaOfZone() const { return getAreaOfZoneFast(); } void fillBounds(Bounds& output) const; void getAllNodes(std::set& output) const; + void getBarycenter(double *bary, double& weigh) const; double getAreaOfZoneFast() const { double ret=_ptr->getAreaOfZone(); return _direction?ret:-ret; } AbstractEdge *clone() const; int recursiveSize() const; @@ -39,9 +40,9 @@ namespace INTERP_KERNEL double getCurveLength() const { return _ptr->getCurveLength(); } bool changeEndNodeWith(Node *node) const; bool changeStartNodeWith(Node *node) const; - void dumpInXfigFile(std::ostream& stream) const; bool intresicEqual(const AbstractEdge *other) const; bool intresicEqualDirSensitive(const AbstractEdge *other) const; + void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; bool getDirection() const { return _direction; } bool intresincEqCoarse(const Edge *other) const; private: diff --git a/src/INTERP_KERNEL/Geometric2D/Makefile.am b/src/INTERP_KERNEL/Geometric2D/Makefile.am index 75551d0cd..65a017b53 100644 --- a/src/INTERP_KERNEL/Geometric2D/Makefile.am +++ b/src/INTERP_KERNEL/Geometric2D/Makefile.am @@ -32,6 +32,7 @@ noinst_LTLIBRARIES = libInterpGeometric2DAlg.la dist_libInterpGeometric2DAlg_la_SOURCES = \ AbstractEdge.cxx \ Bounds.cxx \ +Precision.cxx \ ComposedEdge.cxx \ EdgeArcCircle.cxx \ Edge.cxx \ diff --git a/src/INTERP_KERNEL/Geometric2D/Node.cxx b/src/INTERP_KERNEL/Geometric2D/Node.cxx index ed60e38bf..1f88b6d94 100644 --- a/src/INTERP_KERNEL/Geometric2D/Node.cxx +++ b/src/INTERP_KERNEL/Geometric2D/Node.cxx @@ -4,8 +4,6 @@ using namespace std; using namespace INTERP_KERNEL; -const double Node::CST_FOR_XFIG=1e4; - Node::Node(double x, double y):_isToDel(true),_cnt(1),_loc(UNKNOWN) { const unsigned SPACEDIM=2; @@ -25,7 +23,7 @@ Node::Node(std::istream& stream):_isToDel(true),_cnt(1),_loc(UNKNOWN) { int tmp; stream >> tmp; - _coords[i]=((double) tmp)/CST_FOR_XFIG; + _coords[i]=((double) tmp)/1e4; } } @@ -76,9 +74,9 @@ bool Node::isEqualAndKeepTrack(const Node& other, std::vector& track) co return ret; } -void Node::dumpInXfigFile(std::ostream& stream) const +void Node::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const { - stream << (int)(_coords[0]*CST_FOR_XFIG) << " " << (int)(_coords[1]*CST_FOR_XFIG) << " "; + stream << box.fitXForXFig(_coords[0],resolution) << " " << box.fitYForXFig(_coords[1],resolution) << " "; } double Node::distanceWithSq(const Node& other) const diff --git a/src/INTERP_KERNEL/Geometric2D/Node.hxx b/src/INTERP_KERNEL/Geometric2D/Node.hxx index 45c4025c1..5c5c9fadf 100644 --- a/src/INTERP_KERNEL/Geometric2D/Node.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Node.hxx @@ -18,6 +18,8 @@ namespace INTERP_KERNEL OUT_1 = 10, UNKNOWN = 11 } TypeOfLocInPolygon; + + class Bounds; /*! * As nodes can be shared between edges it is dealed with ref counting. @@ -40,7 +42,7 @@ namespace INTERP_KERNEL bool isEqual(const Node& other) const; double getSlope(const Node& other) const; bool isEqualAndKeepTrack(const Node& other, std::vector& track) const; - void dumpInXfigFile(std::ostream& stream) const; + void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; double distanceWithSq(const Node& other) const; double operator[](int i) const { return _coords[i]; } //!for tests only ! @@ -57,8 +59,6 @@ namespace INTERP_KERNEL mutable unsigned char _cnt; mutable TypeOfLocInPolygon _loc; double *_coords; - public: - static const double CST_FOR_XFIG; }; } diff --git a/src/INTERP_KERNEL/Geometric2D/Precision.cxx b/src/INTERP_KERNEL/Geometric2D/Precision.cxx new file mode 100644 index 000000000..cb67b2016 --- /dev/null +++ b/src/INTERP_KERNEL/Geometric2D/Precision.cxx @@ -0,0 +1,15 @@ +#include "Precision.hxx" + +double INTERP_KERNEL::QUADRATIC_PLANAR::_precision=1e-14; + +double INTERP_KERNEL::QUADRATIC_PLANAR::_arcDetectionPrecision=1e-14; + +void INTERP_KERNEL::QUADRATIC_PLANAR::setPrecision(double precision) +{ + INTERP_KERNEL::QUADRATIC_PLANAR::_precision=precision; +} + +void INTERP_KERNEL::QUADRATIC_PLANAR::setArcDetectionPrecision(double precision) +{ + INTERP_KERNEL::QUADRATIC_PLANAR::_arcDetectionPrecision=precision; +} diff --git a/src/INTERP_KERNEL/Geometric2D/Precision.hxx b/src/INTERP_KERNEL/Geometric2D/Precision.hxx index 5eb791359..9eb3d275e 100644 --- a/src/INTERP_KERNEL/Geometric2D/Precision.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Precision.hxx @@ -1,7 +1,11 @@ namespace INTERP_KERNEL { - namespace QUADRATIC_PLANAR + class QUADRATIC_PLANAR { - static double _precision=1e-14; - } + public: + static double _precision; + static double _arcDetectionPrecision; + static void setPrecision(double precision); + static void setArcDetectionPrecision(double precision); + }; } diff --git a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx index 712eb8c0f..c12a58678 100644 --- a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx @@ -65,19 +65,25 @@ QuadraticPolygon *QuadraticPolygon::buildArcCirclePolygon(std::vector& n { EdgeLin *e1,*e2; e1=new EdgeLin(nodes[i],nodes[i+size/2]); - e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%size]); + e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]); SegSegIntersector inters(*e1,*e2); - bool colinearity,areOverlapped; - inters.areOverlappedOrOnlyColinears(0,colinearity,areOverlapped); + bool colinearity=inters.areColinears(); + delete e1; delete e2; if(colinearity) - ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size])); + ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)])); else - ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%size])); + ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)])); nodes[i]->decrRef(); nodes[i+size/2]->decrRef(); } return ret; } +void QuadraticPolygon::closeMe() const +{ + if(!front()->changeStartNodeWith(back()->getEndNode())) + throw(Exception("big error: not closed polygon...")); +} + void QuadraticPolygon::circularPermute() { vector::iterator iter1=_subEdges.begin(); @@ -112,24 +118,45 @@ double QuadraticPolygon::getPerimeterFast() const return ret; } -double QuadraticPolygon::getHydroulicDiameter() const +double QuadraticPolygon::getHydraulicDiameter() const { return 4*fabs(getAreaFast())/getPerimeterFast(); } +void QuadraticPolygon::dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const +{ + ofstream file(fileName); + const int resolution=1200; + Bounds box; + box.prepareForAggregation(); + fillBounds(box); + other.fillBounds(box); + dumpInXfigFile(file,resolution,box); + other.ComposedEdge::dumpInXfigFile(file,resolution,box); +} + void QuadraticPolygon::dumpInXfigFile(const char *fileName) const { ofstream file(fileName); - file << "#FIG 3.2 Produced by xfig version 3.2.5-alpha5" << endl; - file << "Landscape" << endl; - file << "Center" << endl; - file << "Inches" << endl; - file << "Letter" << endl; - file << "100.00" << endl; - file << "Single" << endl; - file << "-2" << endl; - file << "1200 2" << endl; - ComposedEdge::dumpInXfigFile(file); + const int resolution=1200; + Bounds box; + box.prepareForAggregation(); + fillBounds(box); + dumpInXfigFile(file,resolution,box); +} + +void QuadraticPolygon::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const +{ + stream << "#FIG 3.2 Produced by xfig version 3.2.5-alpha5" << endl; + stream << "Landscape" << endl; + stream << "Center" << endl; + stream << "Metric" << endl; + stream << "Letter" << endl; + stream << "100.00" << endl; + stream << "Single" << endl; + stream << "-2" << endl; + stream << resolution << " 2" << endl; + ComposedEdge::dumpInXfigFile(stream,resolution,box); } double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const @@ -138,7 +165,7 @@ double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const vector polygs=intersectMySelfWith(other); for(vector::iterator iter=polygs.begin();iter!=polygs.end();iter++) { - ret+=fabs((*iter)->getAreaFast()); + ret+=fabs((*iter)->getAreaOfZone()); delete *iter; } return ret; @@ -197,10 +224,10 @@ void QuadraticPolygon::splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticP c2=new ComposedEdgeWithIt; } else - { - updateNeighbours(merge,it1,it2,curE1,curE2); - it1.next(); - } + { + updateNeighbours(merge,it1,it2,curE1,curE2); + it1.next(); + } } } ComposedEdge::Delete(c1); delete c2; diff --git a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx index 3232b914f..9096db02d 100644 --- a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx @@ -19,12 +19,14 @@ namespace INTERP_KERNEL static QuadraticPolygon *buildLinearPolygon(std::vector& nodes); static QuadraticPolygon *buildArcCirclePolygon(std::vector& nodes); ~QuadraticPolygon(); + void closeMe() const; void circularPermute(); //! warning : use it if and only if this is composed of ElementaryEdges only : typical case. double getAreaFast() const; double getPerimeterFast() const; - double getHydroulicDiameter() const; + double getHydraulicDiameter() const; void dumpInXfigFile(const char *fileName) const; + void dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const; double intersectWith(const QuadraticPolygon& other) const; std::vector intersectMySelfWith(const QuadraticPolygon& other) const; public://Only public for tests reasons @@ -34,6 +36,7 @@ namespace INTERP_KERNEL bool amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction); protected: std::list zipConsecutiveInSegments() const; + void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; void closePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, std::vector& results) const; static void updateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, const AbstractEdge *e1, const AbstractEdge *e2); -- 2.39.2