From 3738dbc05d50d62eb003079c74af0fcab18274e9 Mon Sep 17 00:00:00 2001 From: ageay Date: Mon, 22 Sep 2008 14:16:32 +0000 Subject: [PATCH] *** empty log message *** --- src/INTERP_KERNEL/Geometric2D/Bounds.cxx | 6 ++ src/INTERP_KERNEL/Geometric2D/Bounds.hxx | 3 +- .../Geometric2D/ComposedEdge.cxx | 56 +++++++++++++++++++ .../Geometric2D/ComposedEdge.hxx | 5 ++ src/INTERP_KERNEL/Geometric2D/Edge.cxx | 6 ++ src/INTERP_KERNEL/Geometric2D/Edge.hxx | 2 + .../Geometric2D/EdgeArcCircle.cxx | 6 ++ .../Geometric2D/EdgeArcCircle.hxx | 1 + .../Geometric2D/ElementaryEdge.hxx | 1 + src/INTERP_KERNEL/Geometric2D/Node.cxx | 6 ++ src/INTERP_KERNEL/Geometric2D/Node.hxx | 1 + .../Geometric2D/QuadraticPolygon.cxx | 23 ++++++++ .../Geometric2D/QuadraticPolygon.hxx | 1 + .../Test/QuadraticPlanarInterpTest4.cxx | 23 +++++++- 14 files changed, 138 insertions(+), 2 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/Bounds.cxx b/src/INTERP_KERNEL/Geometric2D/Bounds.cxx index dbce54efc..3def5ad51 100644 --- a/src/INTERP_KERNEL/Geometric2D/Bounds.cxx +++ b/src/INTERP_KERNEL/Geometric2D/Bounds.cxx @@ -46,6 +46,12 @@ double Bounds::getDiagonal() const return sqrt(a*a+b*b); } +void Bounds::getBarycenter(double& xBary, double& yBary) const +{ + xBary=(_xMin+_xMax)/2.; + yBary=(_yMax+_yMin)/2.; +} + void Bounds::prepareForAggregation() { _xMin=1e200; _xMax=-1e200; _yMin=1e200; _yMax=-1e200; diff --git a/src/INTERP_KERNEL/Geometric2D/Bounds.hxx b/src/INTERP_KERNEL/Geometric2D/Bounds.hxx index 379f30b36..85d6378cd 100644 --- a/src/INTERP_KERNEL/Geometric2D/Bounds.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Bounds.hxx @@ -17,10 +17,11 @@ namespace INTERP_KERNEL class Bounds { public: - Bounds():_xMin(0),_xMax(0.),_yMin(0.),_yMax(0.) { } + Bounds():_xMin(0.),_xMax(0.),_yMin(0.),_yMax(0.) { } double &operator[](int i); const double& operator[](int i) const; double getDiagonal() const; + void getBarycenter(double& xBary, double& yBary) const; Bounds& operator=(const Bounds& other) { _xMin=other._xMin; _xMax=other._xMax; _yMin=other._yMin; _yMax=other._yMax; return *this; } 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; } diff --git a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx index b866ee643..6b69bad8c 100644 --- a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.cxx @@ -129,6 +129,20 @@ double ComposedEdge::getHydraulicDiameter() const return 4*fabs(getArea())/getPerimeter(); } +double ComposedEdge::normalize(ComposedEdge *other) +{ + Bounds b; + b.prepareForAggregation(); + fillBounds(b); + other->fillBounds(b); + double dimChar=b.getDiagonal(); + double xBary,yBary; + b.getBarycenter(xBary,yBary); + applySimilarity(xBary,yBary,dimChar); + other->applySimilarity(xBary,yBary,dimChar); + return dimChar; +} + void ComposedEdge::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const { stream.precision(10); @@ -162,6 +176,12 @@ void ComposedEdge::fillBounds(Bounds& output) const (*iter)->fillBounds(output); } +void ComposedEdge::applySimilarity(double xBary, double yBary, double dimChar) +{ + for(list::iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++) + (*iter)->applySimilarity(xBary,yBary,dimChar); +} + /*! * @param part1 INOUT parameter. If an edge in 'this' is found with an id in 'ids1', 'part1' is \b incremented. * @param part2 INOUT parameter. If an edge in 'this' is found with an id in 'ids2', 'part2' is \b incremented. @@ -188,6 +208,42 @@ void ComposedEdge::dispatchPerimeter(const std::set& ids1, const std::set& result) const +{ + double ret=0; + for(list::const_iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++) + ret+=father.checkFatherHood(*iter,result); + return ret; +} + +/*! + * Givent edge 'edge' is one of subedge is father of this edge AND NOT ON on the corresponding rank in 'this' the curvelength is added. + * If this sub edge is son and ON the corresponding length is returned by return param. + */ +double ComposedEdge::checkFatherHood(ElementaryEdge *edge, std::vector& result) const +{ + double ret=0.; + int rank=0; + int refId=edge->getPtr()->getId(); + for(list::const_iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++,rank++) + { + int curId=(*iter)->getPtr()->getId(); + if(curId==refId) + { + if((*iter)->getLoc()!=FULL_ON_1) + result[rank]+=edge->getCurveLength(); + else + ret+=edge->getCurveLength(); + break; + } + } + return ret; +} + void ComposedEdge::fillAllEdgeIds(std::set& ids) const { for(list::const_iterator iter=_subEdges.begin();iter!=_subEdges.end();iter++) diff --git a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx index fd3e9c8be..b02d8f6ff 100644 --- a/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/ComposedEdge.hxx @@ -3,6 +3,7 @@ #include #include +#include namespace INTERP_KERNEL { @@ -29,9 +30,13 @@ namespace INTERP_KERNEL double getArea() const; double getPerimeter() const; double getHydraulicDiameter() const; + double normalize(ComposedEdge *other); void fillBounds(Bounds& output) const; + void applySimilarity(double xBary, double yBary, double dimChar); void dispatchPerimeter(const std::set& ids1, const std::set& ids2, double& part1, double& part2, double& commonPart) const; + double dispatchPerimeterAdv(const ComposedEdge& father, std::vector& result) const; + double checkFatherHood(ElementaryEdge *edge, std::vector& result) const; void fillAllEdgeIds(std::set& ids) const; void getAllNodes(std::set& output) const; void getBarycenter(double *bary, double& weigh) const; diff --git a/src/INTERP_KERNEL/Geometric2D/Edge.cxx b/src/INTERP_KERNEL/Geometric2D/Edge.cxx index bdf640f5f..708873564 100644 --- a/src/INTERP_KERNEL/Geometric2D/Edge.cxx +++ b/src/INTERP_KERNEL/Geometric2D/Edge.cxx @@ -644,6 +644,12 @@ Intersector *Edge::buildIntersectorWith(const Edge *e1, const Edge *e2) return ret; } +void Edge::applySimilarity(double xBary, double yBary, double dimChar) +{ + _start->applySimilarity(xBary,yBary,dimChar); + _end->applySimilarity(xBary,yBary,dimChar); +} + bool Edge::intersect(const Edge *f1, const Edge *f2, Intersector *intersector, const Bounds *whereToFind, MergePoints& commonNode, ComposedEdge& outValForF1, ComposedEdge& outValForF2) { diff --git a/src/INTERP_KERNEL/Geometric2D/Edge.hxx b/src/INTERP_KERNEL/Geometric2D/Edge.hxx index 2a0b02a2b..0f95773bc 100644 --- a/src/INTERP_KERNEL/Geometric2D/Edge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Edge.hxx @@ -193,6 +193,8 @@ namespace INTERP_KERNEL virtual void update(Node *m) = 0; //! returns area between this and axe Ox delimited along Ox by _start and _end. virtual double getAreaOfZone() const = 0; + //! apply a similiraty transformation on 'this' + virtual void applySimilarity(double xBary, double yBary, double dimChar); //! return the length of arc. Value is always > 0. ! virtual double getCurveLength() const = 0; virtual void getBarycenter(double *bary) const = 0; diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx index 188aa66c4..3d94d2021 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.cxx @@ -460,6 +460,12 @@ Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) return new EdgeArcCircle(start,end,_center,_radius,angle0,deltaAngle,direction); } +void EdgeArcCircle::applySimilarity(double xBary, double yBary, double dimChar) +{ + Edge::applySimilarity(xBary,yBary,dimChar); + _radius/=dimChar; +} + double EdgeArcCircle::getAbsoluteAngle(const double *vect, double& normVect) { normVect=Node::norm(vect); diff --git a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx index 307d5d875..90bdbcf18 100644 --- a/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/EdgeArcCircle.hxx @@ -75,6 +75,7 @@ namespace INTERP_KERNEL const double *getCenter() const { return _center; } void getCenter(double *center) const { center[0]=_center[0]; center[1]=_center[1]; } bool doIHaveSameDirectionAs(const Edge& other) const { return false; } + void applySimilarity(double xBary, double yBary, double dimChar); double getAngle0() const { return _angle0; } double getRadius() const { return _radius; } double getAngle() const { return _angle; } diff --git a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx index 860c764f1..8657632de 100644 --- a/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx +++ b/src/INTERP_KERNEL/Geometric2D/ElementaryEdge.hxx @@ -25,6 +25,7 @@ namespace INTERP_KERNEL bool isNodeIn(Node *n) const; double getAreaOfZone() const { double ret=_ptr->getAreaOfZone(); return _direction?ret:-ret; } void fillBounds(Bounds& output) const; + void applySimilarity(double xBary, double yBary, double dimChar) { _ptr->applySimilarity(xBary,yBary,dimChar); } void getAllNodes(std::set& output) const; void getBarycenter(double *bary, double& weigh) const; ElementaryEdge *clone() const; diff --git a/src/INTERP_KERNEL/Geometric2D/Node.cxx b/src/INTERP_KERNEL/Geometric2D/Node.cxx index 4afeb667f..5a0aaf022 100644 --- a/src/INTERP_KERNEL/Geometric2D/Node.cxx +++ b/src/INTERP_KERNEL/Geometric2D/Node.cxx @@ -109,3 +109,9 @@ double Node::computeAngle(const double *pt1, const double *pt2) else return -ret; } + +void Node::applySimilarity(double xBary, double yBary, double dimChar) +{ + _coords[0]=(_coords[0]+xBary)/dimChar; + _coords[1]=(_coords[1]+yBary)/dimChar; +} diff --git a/src/INTERP_KERNEL/Geometric2D/Node.hxx b/src/INTERP_KERNEL/Geometric2D/Node.hxx index 01e7b839b..961090050 100644 --- a/src/INTERP_KERNEL/Geometric2D/Node.hxx +++ b/src/INTERP_KERNEL/Geometric2D/Node.hxx @@ -53,6 +53,7 @@ namespace INTERP_KERNEL static double computeSlope(const double *pt1, const double *pt2); //returns an angle in -Pi;Pi static double computeAngle(const double *pt1, const double *pt2); + void applySimilarity(double xBary, double yBary, double dimChar); static double dot(const double *vect1, const double *vect2) { return vect1[0]*vect2[0]+vect1[1]*vect2[1]; } static double sign(double val) { if(val>=0) return 1.; else return -1.; } static double norm(const double *vect) { return sqrt(vect[0]*vect[0]+vect[1]*vect[1]); } diff --git a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx index 0ab60b1ac..3f0c2e890 100644 --- a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.cxx @@ -173,6 +173,29 @@ void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, doub } } +/*! + * polThis.size()==this->size() and polOther.size()==other.size(). + * For each ElementaryEdge of 'this', the corresponding UNIQUE contribution in resulting polygon is in 'polThis'. + * For each ElementaryEdge of 'other', the corresponding UNIQUE contribution in resulting polygon is in 'polOther'. + * The returned value is the sum of common edge lengths. + */ +double QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther, double& area) const +{ + area=0.; + double ret=0.; + polThis.resize(size()); + polOther.resize(other.size()); + vector polygs=intersectMySelfWith(other); + for(vector::iterator iter=polygs.begin();iter!=polygs.end();iter++) + { + area+=fabs((*iter)->getArea()); + ret+=(*iter)->dispatchPerimeterAdv(*this,polThis); + ret+=(*iter)->dispatchPerimeterAdv(other,polOther); + delete *iter; + } + return ret; +} + std::vector QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const { QuadraticPolygon cpyOfThis(*this); diff --git a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx index 41c4556ea..301c991a8 100644 --- a/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/QuadraticPolygon.hxx @@ -30,6 +30,7 @@ namespace INTERP_KERNEL double intersectWith(const QuadraticPolygon& other) const; std::vector intersectMySelfWith(const QuadraticPolygon& other) const; void intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart, double& area) const; + double intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther, double& area) const; public://Only public for tests reasons void performLocatingOperation(QuadraticPolygon& pol2) const; void splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits) const; diff --git a/src/INTERP_KERNEL/Test/QuadraticPlanarInterpTest4.cxx b/src/INTERP_KERNEL/Test/QuadraticPlanarInterpTest4.cxx index 0adf79e1b..5dd6639f3 100644 --- a/src/INTERP_KERNEL/Test/QuadraticPlanarInterpTest4.cxx +++ b/src/INTERP_KERNEL/Test/QuadraticPlanarInterpTest4.cxx @@ -23,7 +23,7 @@ void QuadraticPlanarInterpTest::checkPolygonsIntersection1() // vector result; for(int k=0;k<2;k++) - for(int i=0;i<1;i++) + for(int i=0;i<3;i++) { for(int j=0;j<1;j++) { @@ -37,6 +37,21 @@ void QuadraticPlanarInterpTest::checkPolygonsIntersection1() CPPUNIT_ASSERT_EQUAL(3,result[0]->recursiveSize()); double tmp1=0.,tmp2=0.,tmp3=0.,tmp4=0.; pol1.intersectForPerimeter(pol2,tmp1,tmp2,tmp3,tmp4); + vector v1,v2; + double area2; + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,pol1.intersectForPerimeterAdvanced(pol2,v1,v2,area2),1.e-14);//no common edge + CPPUNIT_ASSERT_EQUAL(3,(int)v1.size()); + CPPUNIT_ASSERT_EQUAL(3,(int)v2.size()); + if(k==0) + { + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.7,v1[(3-i)%3],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,v1[(4-i)%3],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,v1[(5-i)%3],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,v2[0],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.78262379212492639,v2[1],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.78262379212492639,v2[2],1.e-14); + } + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.245,area2,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.7,tmp1,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(1.5652475842498528,tmp2,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,tmp3,1.e-14);//no common edge @@ -124,6 +139,12 @@ void QuadraticPlanarInterpTest::checkPolygonsIntersection1() delete result[0]; double tmp1=0.,tmp2=0.,tmp3=0.,tmp4=0.; pol7.intersectForPerimeter(pol8,tmp1,tmp2,tmp3,tmp4); + vector v1,v2; + double area2; + CPPUNIT_ASSERT_DOUBLES_EQUAL(3.2360679774997898,pol7.intersectForPerimeterAdvanced(pol8,v1,v2,area2),1.e-14);//only common edges. + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,v1[0]+v1[1]+v1[2],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,v2[0]+v2[1]+v2[2],1.e-14); + CPPUNIT_ASSERT_DOUBLES_EQUAL(0.5,area2,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,tmp1,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(0.,tmp2,1.e-14); CPPUNIT_ASSERT_DOUBLES_EQUAL(3.2360679774997898,tmp3,1.e-14); -- 2.39.2