From 5eef5b0d951c208c109afee14987d0de8ebee8df Mon Sep 17 00:00:00 2001 From: ageay Date: Thu, 16 Feb 2012 13:44:00 +0000 Subject: [PATCH] Tesselate. --- .../Geometric2D/InterpKernelGeo2DBounds.cxx | 8 +- .../InterpKernelGeo2DEdgeArcCircle.cxx | 169 +++++++++----- .../InterpKernelGeo2DEdgeArcCircle.hxx | 20 +- .../Geometric2D/InterpKernelGeo2DNode.cxx | 4 +- .../InterpKernelGeo2DQuadraticPolygon.cxx | 32 +-- .../InterpKernelGeo2DQuadraticPolygon.hxx | 16 +- src/INTERP_KERNEL/Geometric2DIntersector.txx | 28 +-- .../PointLocator2DIntersector.txx | 18 +- src/MEDCoupling/MEDCouplingPointSet.cxx | 4 +- src/MEDCoupling/MEDCouplingUMesh.cxx | 218 ++++++++++++++++-- src/MEDCoupling/MEDCouplingUMesh.hxx | 3 + src/MEDCoupling_Swig/MEDCoupling.i | 2 + 12 files changed, 373 insertions(+), 149 deletions(-) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx index 64dd24818..d450601da 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DBounds.cxx @@ -113,18 +113,18 @@ void Bounds::getInterceptedArc(const double *center, double radius, double& intr double v1[2],v2[2],w1[2],w2[2]; v1[0]=_x_min-center[0]; v1[1]=_y_max-center[1]; v2[0]=_x_max-center[0]; v2[1]=_y_min-center[1]; w1[0]=v1[0]; w1[1]=_y_min-center[1]; w2[0]=v2[0]; w2[1]=_y_max-center[1]; - double delta1=EdgeArcCircle::safeAsin(v1[0]*v2[1]-v1[1]*v2[0]); - double delta2=EdgeArcCircle::safeAsin(w1[0]*w2[1]-w1[1]*w2[0]); + double delta1=EdgeArcCircle::SafeAsin(v1[0]*v2[1]-v1[1]*v2[0]); + double delta2=EdgeArcCircle::SafeAsin(w1[0]*w2[1]-w1[1]*w2[0]); double tmp; if(fabs(delta1)>fabs(delta2)) { intrcptArcDelta=delta1; - intrcptArcAngle0=EdgeArcCircle::getAbsoluteAngle(v1,tmp); + intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngle(v1,tmp); } else { intrcptArcDelta=delta2; - intrcptArcAngle0=EdgeArcCircle::getAbsoluteAngle(w1,tmp); + intrcptArcAngle0=EdgeArcCircle::GetAbsoluteAngle(w1,tmp); } } } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx index e884131f2..9f5b846c7 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.cxx @@ -53,7 +53,7 @@ void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& { if(obvious1) { - if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) + if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) whereEnd=INSIDE; else whereEnd=OUT_AFTER; @@ -61,31 +61,31 @@ void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& } else { - if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart)) + if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart)) whereStart=INSIDE; else whereStart=OUT_BEFORE; return ; } } - if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart)) + if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadStart)) { whereStart=INSIDE; - if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) + if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) whereEnd=INSIDE; else whereEnd=OUT_AFTER; } else {//we are out in start. - if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) + if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd)) { whereStart=OUT_BEFORE; whereEnd=INSIDE; } else { - if(EdgeArcCircle::isIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0())) + if(EdgeArcCircle::IsIn2Pi(getE2().getAngle0(),getE2().getAngle(),getE1().getAngle0())) {//_e2 contains stictly _e1 whereStart=OUT_BEFORE; whereEnd=OUT_AFTER; @@ -104,7 +104,7 @@ void ArcCArcCIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& */ double ArcCArcCIntersector::getAngle(Node *node) const { - return EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius()); + return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(((*node)[0]-getE1().getCenter()[0])/getE1().getRadius(),((*node)[1]-getE1().getCenter()[1])/getE1().getRadius()); } bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeArcCircle& a2) @@ -142,17 +142,17 @@ bool ArcCArcCIntersector::areArcsOverlapped(const EdgeArcCircle& a1, const EdgeA else return false; } - double phi=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp); + double phi=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect((centerL[0]-centerB[0])/tmp,(centerL[1]-centerB[1])/tmp); double cst2=2*radiusL*tmp/(radiusB*radiusB); double cmpContainer[4]; int sizeOfCmpContainer=2; cmpContainer[0]=cst+cst2*cos(phi-angle0L); cmpContainer[1]=cst+cst2*cos(phi-angle0L+angleL); - double a=EdgeArcCircle::normalizeAngle(phi-angle0L); - if(EdgeArcCircle::isIn2Pi(angle0L,angleL,a)) + double a=EdgeArcCircle::NormalizeAngle(phi-angle0L); + if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a)) cmpContainer[sizeOfCmpContainer++]=cst+cst2; - a=EdgeArcCircle::normalizeAngle(phi-angle0L+M_PI); - if(EdgeArcCircle::isIn2Pi(angle0L,angleL,a)) + a=EdgeArcCircle::NormalizeAngle(phi-angle0L+M_PI); + if(EdgeArcCircle::IsIn2Pi(angle0L,angleL,a)) cmpContainer[sizeOfCmpContainer++]=cst-cst2; a=*std::max_element(cmpContainer,cmpContainer+sizeOfCmpContainer); return Node::areDoubleEqualsWP(a,1.,2.); @@ -189,9 +189,9 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi double d1_1=(_dist*_dist-radius2*radius2+radius1*radius1)/(2.*_dist); 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 d1_1y=EdgeArcCircle::safeSqrt(radius1*radius1-d1_1*d1_1); - double angleE1=EdgeArcCircle::normalizeAngle(getE1().getAngle0()+getE1().getAngle()); - double angleE2=EdgeArcCircle::normalizeAngle(getE2().getAngle0()+getE2().getAngle()); + double d1_1y=EdgeArcCircle::SafeSqrt(radius1*radius1-d1_1*d1_1); + double angleE1=EdgeArcCircle::NormalizeAngle(getE1().getAngle0()+getE1().getAngle()); + double angleE2=EdgeArcCircle::NormalizeAngle(getE2().getAngle0()+getE2().getAngle()); if(!Node::areDoubleEquals(d1_1y,0)) { //2 intersections @@ -200,13 +200,13 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi 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::getAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); - double angle2_1=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1); + double angle1_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); + double angle2_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius1,v2[1]/radius1); 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::getAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2); - double angle2_2=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2); + double angle1_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v3[0]/radius2,v3[1]/radius2); + double angle2_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v4[0]/radius2,v4[1]/radius2); // bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1); bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1); @@ -226,8 +226,8 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi 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::getAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); - double angle0_2=EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2); + double angle0_1=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v1[0]/radius1,v1[1]/radius1); + double angle0_2=EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(v2[0]/radius2,v2[1]/radius2); 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); @@ -261,19 +261,19 @@ std::list< IntersectElement > ArcCArcCIntersector::getIntersectionsCharacteristi angle0_2=angle0_1; signDeltaAngle2=1.; } - angle0_1=normalizeAngle(angle0_1); - angle0_2=normalizeAngle(angle0_2); - double angleE1=normalizeAngle(getE1().getAngle0()+getE1().getAngle()); - double angleE2=normalizeAngle(getE2().getAngle0()+getE2().getAngle()); + angle0_1=NormalizeAngle(angle0_1); + angle0_2=NormalizeAngle(angle0_2); + double angleE1=NormalizeAngle(getE1().getAngle0()+getE1().getAngle()); + double angleE2=NormalizeAngle(getE2().getAngle0()+getE2().getAngle()); if(!(Node::areDoubleEquals(d1_1,radius1) || Node::areDoubleEquals(d1_1,-radius1)) ) { //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 - double angle2_1=normalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1 - double angle1_2=normalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2 - double angle2_2=normalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2 + 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 + double angle2_1=NormalizeAngle(angle0_1-deltaAngle1);// Intersection 2 seen for _e1 + double angle1_2=NormalizeAngle(angle0_2+signDeltaAngle2*deltaAngle2);// Intersection 1 seen for _e2 + double angle2_2=NormalizeAngle(angle0_2-signDeltaAngle2*deltaAngle2);// Intersection 2 seen for _e2 // bool e1_1S=Node::areDoubleEqualsWP(angle1_1,getE1().getAngle0(),radius1); bool e1_1E=Node::areDoubleEqualsWP(angle1_1,angleE1,radius1); @@ -333,7 +333,7 @@ std::list< IntersectElement > ArcCSegIntersector::getIntersectionsCharacteristic const double *center=getE1().getCenter(); if(!(fabs(_determinant)<(2.*QUADRATIC_PLANAR::_precision)))//QUADRATIC_PLANAR::_precision*QUADRATIC_PLANAR::_precision*_drSq*_drSq/(2.*_dx*_dx)) { - double determinant=EdgeArcCircle::safeSqrt(_determinant); + double determinant=EdgeArcCircle::SafeSqrt(_determinant); double x1=(_cross*_dy/_drSq+Node::sign(_dy)*_dx*determinant)+center[0]; double y1=(-_cross*_dx/_drSq+fabs(_dy)*determinant)+center[1]; Node *intersect1=new Node(x1,y1); intersect1->declareOn(); @@ -375,21 +375,21 @@ EdgeArcCircle::EdgeArcCircle(std::istream& lineInXfig) _start=new Node(lineInXfig); Node *middle=new Node(lineInXfig); _end=new Node(lineInXfig); - getArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0); + GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0); middle->decrRef(); updateBounds(); } EdgeArcCircle::EdgeArcCircle(Node *start, Node *middle, Node *end, bool direction):Edge(start,end, direction) { - getArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0); + GetArcOfCirclePassingThru(*_start,*middle,*_end,_center,_radius,_angle,_angle0); updateBounds(); } EdgeArcCircle::EdgeArcCircle(double sX, double sY, double mX, double mY, double eX, double eY):Edge(sX,sY,eX,eY) { double middle[2]; middle[0]=mX; middle[1]=mY; - getArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0); + GetArcOfCirclePassingThru(*_start,middle,*_end,_center,_radius,_angle,_angle0); updateBounds(); } @@ -407,7 +407,7 @@ EdgeArcCircle::EdgeArcCircle(Node *start, Node *end, const double *center, doubl void EdgeArcCircle::changeMiddle(Node *newMiddle) { - getArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0); + GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0); updateBounds(); } @@ -417,8 +417,8 @@ Edge *EdgeArcCircle::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) double sy=((*start)[1]-_center[1])/_radius; double ex=((*end)[0]-_center[0])/_radius; double ey=((*end)[1]-_center[1])/_radius; - double angle0=getAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey); - double deltaAngle=getAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex); + double angle0=GetAbsoluteAngleOfNormalizedVect(direction?sx:ex,direction?sy:ey); + double deltaAngle=GetAbsoluteAngleOfNormalizedVect(sx*ex+sy*ey,sx*ey-sy*ex); if(deltaAngle>0. && _angle<0.) deltaAngle-=2.*M_PI; else if(deltaAngle<0. && _angle>0.) @@ -443,14 +443,61 @@ void EdgeArcCircle::unApplySimilarity(double xBary, double yBary, double dimChar _center[1]=_center[1]*dimChar+yBary; } +/*! + * 'eps' is expected to be > 0. + * 'conn' is of size 3. conn[0] is start id, conn[1] is end id and conn[2] is middle id. + * 'offset' is typically the number of nodes already existing in global 2D curve mesh. Additionnal coords 'addCoo' ids will be put after the already existing. + */ +void EdgeArcCircle::tesselate(const int *conn, int offset, double eps, std::vector& newConn, std::vector& addCoo) const +{ + newConn.push_back(INTERP_KERNEL::NORM_POLYL); + int nbOfSubDiv=fabs(_angle)/eps; + if(nbOfSubDiv<=2) + { + newConn.push_back(conn[0]); newConn.push_back(conn[2]); newConn.push_back(conn[1]); + return ; + } + double signOfAngle=_angle>0.?1.:-1.; + int offset2=offset+((int)addCoo.size())/2; + newConn.push_back(conn[0]); + for(int i=1;idecrRef(); middle->decrRef(); end->decrRef(); + return 0; + } + else + { + EdgeArcCircle *ret=new EdgeArcCircle(start,middle,end); + start->decrRef(); middle->decrRef(); end->decrRef(); + return ret; + } +} + /*! * Given an \b NON normalized vector 'vect', returns its norm 'normVect' and its * angle in ]-Pi,Pi] relative to Ox axe. */ -double EdgeArcCircle::getAbsoluteAngle(const double *vect, double& normVect) +double EdgeArcCircle::GetAbsoluteAngle(const double *vect, double& normVect) { normVect=Node::norm(vect); - return getAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect); + return GetAbsoluteAngleOfNormalizedVect(vect[0]/normVect,vect[1]/normVect); } /*! @@ -460,12 +507,12 @@ double EdgeArcCircle::getAbsoluteAngle(const double *vect, double& normVect) * It is NOT ALWAYS possible to do that only in one call of acos. Sometimes call to asin is necessary * due to imperfection of acos near 0. and Pi (cos x ~ 1-x*x/2.) */ -double EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(double ux, double uy) +double EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(double ux, double uy) { //When arc is lower than 0.707 Using Asin if(fabs(ux)<0.707) { - double ret=safeAcos(ux); + double ret=SafeAcos(ux); if(uy>0.) return ret; ret=-ret; @@ -473,7 +520,7 @@ double EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(double ux, double uy) } else { - double ret=safeAsin(uy); + double ret=SafeAsin(uy); if(ux>0.) return ret; if(ret>0.) @@ -483,7 +530,7 @@ double EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(double ux, double uy) } } -void EdgeArcCircle::getArcOfCirclePassingThru(const double *start, const double *middle, const double *end, +void EdgeArcCircle::GetArcOfCirclePassingThru(const double *start, const double *middle, const double *end, double *center, double& radius, double& angleInRad, double& angleInRad0) { double delta=(middle[0]-start[0])*(end[1]-middle[1])-(end[0]-middle[0])*(middle[1]-start[1]); @@ -491,12 +538,12 @@ 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=safeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1])); - angleInRad0=getAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius); - double angleInRadM=getAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius); - angleInRad=getAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius), + radius=SafeSqrt((start[0]-center[0])*(start[0]-center[0])+(start[1]-center[1])*(start[1]-center[1])); + angleInRad0=GetAbsoluteAngleOfNormalizedVect((start[0]-center[0])/radius,(start[1]-center[1])/radius); + double angleInRadM=GetAbsoluteAngleOfNormalizedVect((middle[0]-center[0])/radius,(middle[1]-center[1])/radius); + angleInRad=GetAbsoluteAngleOfNormalizedVect(((start[0]-center[0])*(end[0]-center[0])+(start[1]-center[1])*(end[1]-center[1]))/(radius*radius), ((start[0]-center[0])*(end[1]-center[1])-(start[1]-center[1])*(end[0]-center[0]))/(radius*radius)); - if(isAngleNotIn(angleInRad0,angleInRad,angleInRadM)) + if(IsAngleNotIn(angleInRad0,angleInRad,angleInRadM)) angleInRad=angleInRad<0?2*M_PI+angleInRad:angleInRad-2*M_PI; } @@ -521,7 +568,7 @@ void EdgeArcCircle::dumpInXfigFile(std::ostream& stream, bool direction, int res void EdgeArcCircle::update(Node *m) { - getArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0); + GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0); updateBounds(); } @@ -587,7 +634,7 @@ void EdgeArcCircle::getBarycenterOfZone(double *bary) const */ bool EdgeArcCircle::isIn(double characterVal) const { - return isIn2Pi(_angle0,_angle,characterVal); + return IsIn2Pi(_angle0,_angle,characterVal); } Node *EdgeArcCircle::buildRepresentantOfMySelf() const @@ -624,13 +671,13 @@ double EdgeArcCircle::getCharactValue(const Node& node) const { double dx=(node[0]-_center[0])/_radius; double dy=(node[1]-_center[1])/_radius; - return getAbsoluteAngleOfNormalizedVect(dx,dy); + return GetAbsoluteAngleOfNormalizedVect(dx,dy); } double EdgeArcCircle::getDistanceToPoint(const double *pt) const { double angle=Node::computeAngle(_center,pt); - if(isIn2Pi(_angle0,_angle,angle)) + if(IsIn2Pi(_angle0,_angle,angle)) return fabs(Node::distanceBtw2Pt(_center,pt)-_radius); else { @@ -646,17 +693,17 @@ bool EdgeArcCircle::isNodeLyingOn(const double *coordOfNode) const if(Node::areDoubleEquals(dist,_radius)) { double angle=Node::computeAngle(_center,coordOfNode); - return isIn2Pi(_angle0,_angle,angle); + return IsIn2Pi(_angle0,_angle,angle); } else return false; } /*! - * Idem isAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[. + * Idem IsAngleNotIn except that here 'start' in ]-Pi;Pi[ and delta in ]-2*Pi;2Pi[. * @param angleIn in ]-Pi;Pi[. */ -bool EdgeArcCircle::isIn2Pi(double start, double delta, double angleIn) +bool EdgeArcCircle::IsIn2Pi(double start, double delta, double angleIn) { double myDelta=angleIn-start; if(delta>0.) @@ -674,7 +721,7 @@ bool EdgeArcCircle::isIn2Pi(double start, double delta, double angleIn) /*! * Given the arc 'a' defined by 'start' angle and a 'delta' [-Pi;Pi] states for the angle 'angleIn' [-Pi;Pi] if it owns or not 'a'. */ -bool EdgeArcCircle::isAngleNotIn(double start, double delta, double angleIn) +bool EdgeArcCircle::IsAngleNotIn(double start, double delta, double angleIn) { double tmp=start; if(tmp<0.) @@ -693,13 +740,13 @@ bool EdgeArcCircle::isAngleNotIn(double start, double delta, double angleIn) void EdgeArcCircle::updateBounds() { _bounds.setValues(std::min((*_start)[0],(*_end)[0]),std::max((*_start)[0],(*_end)[0]),std::min((*_start)[1],(*_end)[1]),std::max((*_start)[1],(*_end)[1])); - if(isIn2Pi(_angle0,_angle,M_PI/2)) + if(IsIn2Pi(_angle0,_angle,M_PI/2)) _bounds[3]=_center[1]+_radius; - if(isIn2Pi(_angle0,_angle,-M_PI/2)) + if(IsIn2Pi(_angle0,_angle,-M_PI/2)) _bounds[2]=_center[1]-_radius; - if(isIn2Pi(_angle0,_angle,0.)) + if(IsIn2Pi(_angle0,_angle,0.)) _bounds[1]=_center[0]+_radius; - if(isIn2Pi(_angle0,_angle,M_PI)) + if(IsIn2Pi(_angle0,_angle,M_PI)) _bounds[0]=_center[0]-_radius; } diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx index fb2efdc41..9657a816a 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeArcCircle.hxx @@ -95,20 +95,22 @@ namespace INTERP_KERNEL double getAngle0() const { return _angle0; } double getRadius() const { return _radius; } double getAngle() const { return _angle; } - static double getAbsoluteAngle(const double *vect, double& normVect); - static double getAbsoluteAngleOfNormalizedVect(double ux, double uy); - static void getArcOfCirclePassingThru(const double *start, const double *middle, const double *end, + void tesselate(const int *conn, int offset, double eps, std::vector& newConn, std::vector& addCoo) const; + static EdgeArcCircle *BuildFromNodes(Node *start, Node *middle, Node *end); + static double GetAbsoluteAngle(const double *vect, double& normVect); + static double GetAbsoluteAngleOfNormalizedVect(double ux, double uy); + 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=std::max(val,0.); return sqrt(ret); } - static double safeAcos(double cosAngle) { double ret=std::min(cosAngle,1.); ret=std::max(ret,-1.); return acos(ret); } - static double safeAsin(double sinAngle) { double ret=std::min(sinAngle,1.); ret=std::max(ret,-1.); return asin(ret); } + static double SafeSqrt(double val) { double ret=std::max(val,0.); return sqrt(ret); } + static double SafeAcos(double cosAngle) { double ret=std::min(cosAngle,1.); ret=std::max(ret,-1.); return acos(ret); } + static double SafeAsin(double sinAngle) { double ret=std::min(sinAngle,1.); ret=std::max(ret,-1.); return asin(ret); } //! @param start and @param angleIn in ]-Pi;Pi] and @param delta in ]-2*Pi,2*Pi[ - static bool isIn2Pi(double start, double delta, double angleIn); + static bool IsIn2Pi(double start, double delta, double angleIn); //! 'delta' 'start' in ]-Pi;Pi[ - static bool isAngleNotIn(double start, double delta, double angleIn); + static bool IsAngleNotIn(double start, double delta, double angleIn); //! for an angle 'angle' in ]-3*Pi;3*Pi[ returns angle in ]-Pi;Pi[ - static double normalizeAngle(double angle) { if(angle>M_PI) return angle-2.*M_PI; if(angle<-M_PI) return angle+2.*M_PI; return angle; } + static double NormalizeAngle(double angle) { if(angle>M_PI) return angle-2.*M_PI; if(angle<-M_PI) return angle+2.*M_PI; return angle; } protected: void updateBounds(); Edge *buildEdgeLyingOnMe(Node *start, Node *end, bool direction=true) const; diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx index 6ebe9a53f..c45939aa6 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DNode.cxx @@ -100,7 +100,7 @@ double Node::computeSlope(const double *pt1, const double *pt2) double x=pt2[0]-pt1[0]; double y=pt2[1]-pt1[1]; double norm=sqrt(x*x+y*y); - double ret=EdgeArcCircle::safeAcos(fabs(x)/norm); + double ret=EdgeArcCircle::SafeAcos(fabs(x)/norm); if( (x>=0. && y>=0.) || (x<0. && y<0.) ) return ret; else @@ -116,7 +116,7 @@ double Node::computeAngle(const double *pt1, const double *pt2) double x=pt2[0]-pt1[0]; double y=pt2[1]-pt1[1]; double norm=sqrt(x*x+y*y); - return EdgeArcCircle::getAbsoluteAngleOfNormalizedVect(x/norm,y/norm); + return EdgeArcCircle::GetAbsoluteAngleOfNormalizedVect(x/norm,y/norm); } /*! diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx index 0e4fc6181..5618cc78a 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx @@ -68,7 +68,7 @@ QuadraticPolygon::~QuadraticPolygon() { } -QuadraticPolygon *QuadraticPolygon::buildLinearPolygon(std::vector& nodes) +QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector& nodes) { QuadraticPolygon *ret=new QuadraticPolygon; std::size_t size=nodes.size(); @@ -80,7 +80,7 @@ QuadraticPolygon *QuadraticPolygon::buildLinearPolygon(std::vector& node return ret; } -QuadraticPolygon *QuadraticPolygon::buildArcCirclePolygon(std::vector& nodes) +QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector& nodes) { QuadraticPolygon *ret=new QuadraticPolygon; std::size_t size=nodes.size(); @@ -101,7 +101,7 @@ QuadraticPolygon *QuadraticPolygon::buildArcCirclePolygon(std::vector& n return ret; } -void QuadraticPolygon::buildDbgFile(const std::vector& nodes, const char *fileName) +void QuadraticPolygon::BuildDbgFile(const std::vector& nodes, const char *fileName) { std::ofstream file(fileName); file << std::setprecision(16); @@ -253,7 +253,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, const std::mapgetDirection()) c1->reverse(); if(!curE2->getDirection()) c2->reverse(); - updateNeighbours(merge,it1,it2,c1,c2); + UpdateNeighbours(merge,it1,it2,c1,c2); //Substitution of simple edge by sub-edges. delete curE1; // <-- destroying simple edge coming from pol1 delete curE2; // <-- destroying simple edge coming from pol2 @@ -269,7 +269,7 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other, const std::mapclone()); int tmp2; - splitPolygonsEachOther(tmp,cpyOfOther,tmp2); + SplitPolygonsEachOther(tmp,cpyOfOther,tmp2); other.performLocatingOperation(tmp); tmp.dispatchPerimeter(polThis[edgeId]); } @@ -664,7 +664,7 @@ void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& oth QuadraticPolygon tmp; tmp.pushBack(curE2->clone()); int tmp2; - splitPolygonsEachOther(tmp,cpyOfThis,tmp2); + SplitPolygonsEachOther(tmp,cpyOfThis,tmp2); performLocatingOperation(tmp); tmp.dispatchPerimeter(polOther[edgeId]); } @@ -689,7 +689,7 @@ void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vec QuadraticPolygon tmp; tmp.pushBack(curE1->clone()); int tmp2; - splitPolygonsEachOther(tmp,cpyOfOther,tmp2); + SplitPolygonsEachOther(tmp,cpyOfOther,tmp2); numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1; } } @@ -703,7 +703,7 @@ std::vector QuadraticPolygon::intersectMySelfWith(const Quad { QuadraticPolygon cpyOfThis(*this); QuadraticPolygon cpyOfOther(other); int nbOfSplits=0; - splitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits); + SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits); //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done. performLocatingOperation(cpyOfOther); return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther); @@ -714,7 +714,7 @@ std::vector QuadraticPolygon::intersectMySelfWith(const Quad * This method perform the minimal splitting so that at the end each edges constituting pol1 are fully either IN or OUT or ON. * @param pol1 IN/OUT param that is equal to 'this' when called. */ -void QuadraticPolygon::splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits) +void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits) { IteratorOnComposedEdge it1(&pol1),it2(&pol2); MergePoints merge; @@ -736,7 +736,7 @@ void QuadraticPolygon::splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticP { if(!curE1->getDirection()) c1->reverse(); if(!curE2->getDirection()) c2->reverse(); - updateNeighbours(merge,it1,it2,c1,c2); + UpdateNeighbours(merge,it1,it2,c1,c2); //Substitution of simple edge by sub-edges. delete curE1; // <-- destroying simple edge coming from pol1 delete curE2; // <-- destroying simple edge coming from pol2 @@ -752,7 +752,7 @@ void QuadraticPolygon::splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticP } else { - updateNeighbours(merge,it1,it2,curE1,curE2); + UpdateNeighbours(merge,it1,it2,curE1,curE2); it1.next(); } } @@ -950,7 +950,7 @@ std::list::iterator QuadraticPolygon::fillAsMuchAsPossibleWi pushBack(tmp); nodeToTest=tmp->getEndNode(); direction?it.nextLoop():it.previousLoop(); - ret=checkInList(nodeToTest,iStart,iEnd); + ret=CheckInList(nodeToTest,iStart,iEnd); if(completed()) return iEnd; } @@ -958,7 +958,7 @@ std::list::iterator QuadraticPolygon::fillAsMuchAsPossibleWi return ret; } -std::list::iterator QuadraticPolygon::checkInList(Node *n, std::list::iterator iStart, +std::list::iterator QuadraticPolygon::CheckInList(Node *n, std::list::iterator iStart, std::list::iterator iEnd) { for(std::list::iterator iter=iStart;iter!=iEnd;iter++) diff --git a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx index 931ba1d9d..e24ee3b12 100644 --- a/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx +++ b/src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.hxx @@ -40,9 +40,9 @@ namespace INTERP_KERNEL QuadraticPolygon() { } QuadraticPolygon(const QuadraticPolygon& other):ComposedEdge(other) { } QuadraticPolygon(const char *fileName); - static QuadraticPolygon *buildLinearPolygon(std::vector& nodes); - static QuadraticPolygon *buildArcCirclePolygon(std::vector& nodes); - static void buildDbgFile(const std::vector& nodes, const char *fileName); + static QuadraticPolygon *BuildLinearPolygon(std::vector& nodes); + static QuadraticPolygon *BuildArcCirclePolygon(std::vector& nodes); + static void BuildDbgFile(const std::vector& nodes, const char *fileName); ~QuadraticPolygon(); void closeMe() const; void circularPermute(); @@ -76,7 +76,7 @@ namespace INTERP_KERNEL void intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const; public://Only public for tests reasons void performLocatingOperation(QuadraticPolygon& pol2) const; - static void splitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits); + static void SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits); std::vector buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const; bool amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction); protected: @@ -84,13 +84,13 @@ namespace INTERP_KERNEL void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const; void closePolygons(std::list& pol2Zip, const QuadraticPolygon& pol1, std::vector& results) const; template - static void updateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, - const EDGES *e1, const EDGES *e2); + static void UpdateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, + const EDGES *e1, const EDGES *e2); std::list::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted, std::list::iterator iStart, std::list::iterator iEnd, bool direction); - static std::list::iterator checkInList(Node *n, std::list::iterator iStart, + static std::list::iterator CheckInList(Node *n, std::list::iterator iStart, std::list::iterator iEnd); }; } @@ -98,7 +98,7 @@ namespace INTERP_KERNEL namespace INTERP_KERNEL { template - void QuadraticPolygon::updateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, + void QuadraticPolygon::UpdateNeighbours(const MergePoints& merger, IteratorOnComposedEdge it1, IteratorOnComposedEdge it2, const EDGES *e1, const EDGES *e2) { it1.previousLoop(); it2.previousLoop(); diff --git a/src/INTERP_KERNEL/Geometric2DIntersector.txx b/src/INTERP_KERNEL/Geometric2DIntersector.txx index d1fd6b09c..4c0e1768b 100644 --- a/src/INTERP_KERNEL/Geometric2DIntersector.txx +++ b/src/INTERP_KERNEL/Geometric2DIntersector.txx @@ -96,12 +96,12 @@ namespace INTERP_KERNEL std::vector nodes2(nbOfSourceNodes); for(int i=0;iintersectWithAbs(*p2); delete p1; delete p2; return ret; @@ -119,8 +119,8 @@ namespace INTERP_KERNEL std::vector nodes2(nbOfSourceNodes); for(int i=0;iintersectWithAbs(*p2); delete p1; delete p2; return ret; @@ -150,12 +150,12 @@ namespace INTERP_KERNEL std::vector nodes2(nbOfTargetNodes); for(int i=0;iintersectWithAbs(*p2,barycenter); delete p1; delete p2; @@ -186,9 +186,9 @@ namespace INTERP_KERNEL for(int i=0;i::_coordsT+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); if(CellModel::GetCellModel(type).isQuadratic()) - return QuadraticPolygon::buildLinearPolygon(nodes); + return QuadraticPolygon::BuildLinearPolygon(nodes); else - return QuadraticPolygon::buildArcCirclePolygon(nodes); + return QuadraticPolygon::BuildArcCirclePolygon(nodes); } INTERSECTOR_TEMPLATE @@ -238,9 +238,9 @@ namespace INTERP_KERNEL for(int i=0;i::_coordsS+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); if(type!=NORM_TRI6 && type!=NORM_QUAD8) - return QuadraticPolygon::buildLinearPolygon(nodes); + return QuadraticPolygon::BuildLinearPolygon(nodes); else - return QuadraticPolygon::buildArcCirclePolygon(nodes); + return QuadraticPolygon::BuildArcCirclePolygon(nodes); } } diff --git a/src/INTERP_KERNEL/PointLocator2DIntersector.txx b/src/INTERP_KERNEL/PointLocator2DIntersector.txx index 1af497b71..ab470829d 100644 --- a/src/INTERP_KERNEL/PointLocator2DIntersector.txx +++ b/src/INTERP_KERNEL/PointLocator2DIntersector.txx @@ -72,9 +72,9 @@ namespace INTERP_KERNEL nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]); QuadraticPolygon *p2; if(!isSourceQuad) - p2=QuadraticPolygon::buildLinearPolygon(nodes2); + p2=QuadraticPolygon::BuildLinearPolygon(nodes2); else - p2=QuadraticPolygon::buildArcCirclePolygon(nodes2); + p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2); double bary[SPACEDIM]; p2->getBarycenter(bary); delete p2; @@ -92,7 +92,7 @@ namespace INTERP_KERNEL std::vector nodes2(nbOfSourceNodes); for(int i=0;igetBarycenterGeneral(bary); delete p; @@ -129,9 +129,9 @@ namespace INTERP_KERNEL for(int i=0;i::_coordsT+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); if(CellModel::GetCellModel(type).isQuadratic()) - return QuadraticPolygon::buildLinearPolygon(nodes); + return QuadraticPolygon::BuildLinearPolygon(nodes); else - return QuadraticPolygon::buildArcCirclePolygon(nodes); + return QuadraticPolygon::BuildArcCirclePolygon(nodes); } INTERSECTOR_TEMPLATE @@ -155,9 +155,9 @@ namespace INTERP_KERNEL for(int i=0;i::_coordsS+OTT::coo2C(startOfCellNodeConn[i])*SPACEDIM); if(type!=NORM_TRI6 && type!=NORM_QUAD8) - return QuadraticPolygon::buildLinearPolygon(nodes); + return QuadraticPolygon::BuildLinearPolygon(nodes); else - return QuadraticPolygon::buildArcCirclePolygon(nodes); + return QuadraticPolygon::BuildArcCirclePolygon(nodes); } } diff --git a/src/MEDCoupling/MEDCouplingPointSet.cxx b/src/MEDCoupling/MEDCouplingPointSet.cxx index 3b7f0f55b..d76908c76 100644 --- a/src/MEDCoupling/MEDCouplingPointSet.cxx +++ b/src/MEDCoupling/MEDCouplingPointSet.cxx @@ -902,9 +902,9 @@ bool MEDCouplingPointSet::isButterfly2DCell(const std::vector& res, bool } INTERP_KERNEL::QuadraticPolygon *pol=0; if(isQuad) - pol=INTERP_KERNEL::QuadraticPolygon::buildArcCirclePolygon(nodes); + pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes); else - pol=INTERP_KERNEL::QuadraticPolygon::buildLinearPolygon(nodes); + pol=INTERP_KERNEL::QuadraticPolygon::BuildLinearPolygon(nodes); bool ret=pol->isButterfly(); delete pol; return ret; diff --git a/src/MEDCoupling/MEDCouplingUMesh.cxx b/src/MEDCoupling/MEDCouplingUMesh.cxx index 11f7bab04..cca5a8fe1 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.cxx +++ b/src/MEDCoupling/MEDCouplingUMesh.cxx @@ -2614,27 +2614,6 @@ namespace ParaMEDMEM } } } - - /*void MEDCouplingUMeshAssignOnLoc(const INTERP_KERNEL::QuadraticPolygon& pol1, INTERP_KERNEL::QuadraticPolygon& pol2, - const int *desc1Bg, const int *desc1End,const std::vector >& intesctEdges1, - const int *desc2Bg, const int *desc2End,const std::vector >& intesctEdges2, - const std::vector >& colinear1) - { - for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++) - { - int eltId1=abs(*desc1)-1; - for(std::vector::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++) - { - std::map::const_iterator it=mappRev.find(*it1); - if(it==mappRev.end()) - { - INTERP_KERNEL::Node *node=MEDCouplingUMeshBuildQPNode(*it1,coo1,offset1,coo2,offset2,addCoo); - mapp[node]=*it1; - mappRev[*it1]=node; - } - } - } - }*/ } /// @endcond @@ -2944,7 +2923,7 @@ DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation2D(con const double *p0=i+1rotate(end,0,angle); @@ -3011,7 +2990,7 @@ DataArrayDouble *MEDCouplingUMesh::fillExtCoordsUsingTranslAndAutoRotation3D(con double p0r[3]={m[0][0]*p0[0]+m[0][1]*p0[1]+m[0][2]*p0[2], m[1][0]*p0[0]+m[1][1]*p0[1]+m[1][2]*p0[2], m[2][0]*p0[0]+m[2][1]*p0[1]+m[2][2]*p0[2]}; double p1r[3]={m[0][0]*p1[0]+m[0][1]*p1[1]+m[0][2]*p1[2], m[1][0]*p1[0]+m[1][1]*p1[1]+m[1][2]*p1[2], m[2][0]*p1[0]+m[2][1]*p1[1]+m[2][2]*p1[2]}; double p2r[3]={m[0][0]*p2[0]+m[0][1]*p2[1]+m[0][2]*p2[2], m[1][0]*p2[0]+m[1][1]*p2[1]+m[1][2]*p2[2], m[2][0]*p2[0]+m[2][1]*p2[1]+m[2][2]*p2[2]}; - INTERP_KERNEL::EdgeArcCircle::getArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0); + INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0r,p1r,p2r,tmp3,radius,alpha,alpha0); double cosangle=i+1rotate(end,vecPlane,angle); @@ -3167,6 +3146,112 @@ void MEDCouplingUMesh::convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exce setConnectivity(newConn,newConnI,false); } +/*! + * This method tessallates 'this' so that the number of cells remains the same. + * This method works only for meshes with spaceDim equal to 2 and meshDim equal to 2. + * If no cells are quadratic in 'this' (INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QPOLYG ) this method will remain unchanged. + * + * \b WARNING this method can lead to a uge amount of nodes if eps is very low. + * @param eps specifies the maximal angle (in radian) between 2 subedges of polylinized edge constituting the input polygon. + */ +void MEDCouplingUMesh::tessellate2D(double eps) throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=2 || getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2D works on umeshes with meshdim equal to 2 and spaceDim equal to 2 too!"); + double epsa=fabs(eps); + if(epsa::min()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !"); + MEDCouplingAutoRefCountObjectPtr desc1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr descIndx1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr revDesc1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr revDescIndx1=DataArrayInt::New(); + MEDCouplingAutoRefCountObjectPtr mDesc=buildDescendingConnectivity2(desc1,descIndx1,revDesc1,revDescIndx1); + revDesc1=0; revDescIndx1=0; + mDesc->tessellate2DCurve(eps); + subDivide2DMesh(mDesc->_nodal_connec->getConstPointer(),mDesc->_nodal_connec_index->getConstPointer(),desc1->getConstPointer(),descIndx1->getConstPointer()); + setCoords(mDesc->getCoords()); +} + +/*! + * This method tessallates 'this' so that the number of cells remains the same. + * This method works only for meshes with spaceDim equal to 2 and meshDim equal to 1. + * If no cells are quadratic in 'this' (INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_QPOLYG ) this method will remain unchanged. + * + * \b WARNING this method can lead to a uge amount of nodes if eps is very low. + * @param eps specifies the maximal angle (in radian) between 2 subedges of polylinized edge constituting the input polygon. + */ +void MEDCouplingUMesh::tessellate2DCurve(double eps) throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=1 || getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve works on umeshes with meshdim equal to 1 and spaceDim equal to 2 too!"); + double epsa=fabs(eps); + if(epsa::min()) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2DCurve : epsilon is null ! Please specify a higher epsilon. If too tiny it can lead to a huge amount of nodes and memory !"); + INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=1.e-10; + int nbCells=getNumberOfCells(); + int nbNodes=getNumberOfNodes(); + const int *conn=_nodal_connec->getConstPointer(); + const int *connI=_nodal_connec_index->getConstPointer(); + const double *coords=_coords->getConstPointer(); + std::vector addCoo; + std::vector newConn; + MEDCouplingAutoRefCountObjectPtr newConnI=DataArrayInt::New(); + newConnI->alloc(nbCells+1,1); + int *newConnIPtr=newConnI->getPointer(); + *newConnIPtr=0; + int tmp1[3]; + INTERP_KERNEL::Node *tmp2[3]; + std::set types; + for(int i=0;itesselate(tmp1,nbNodes,epsa,newConn,addCoo); + types.insert((INTERP_KERNEL::NormalizedCellType)newConn[newConnIPtr[0]]); + delete eac; + newConnIPtr[1]=(int)newConn.size(); + } + else + { + types.insert(INTERP_KERNEL::NORM_SEG2); + newConn.push_back(INTERP_KERNEL::NORM_SEG2); + newConn.insert(newConn.end(),conn+connI[i]+1,conn+connI[i]+3); + newConnIPtr[1]=newConnIPtr[0]+3; + } + } + else + { + types.insert((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]); + newConn.insert(newConn.end(),conn+connI[i],conn+connI[i+1]); + newConnIPtr[1]=newConnIPtr[0]+3; + } + } + if(addCoo.empty() && newConn.size()==_nodal_connec->getNumberOfTuples())//nothing happens during tasselation : no update needed + return ; + _types=types; + DataArrayInt::SetArrayIn(newConnI,_nodal_connec_index); + MEDCouplingAutoRefCountObjectPtr newConnArr=DataArrayInt::New(); + newConnArr->alloc((int)newConn.size(),1); + std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer()); + DataArrayInt::SetArrayIn(newConnArr,_nodal_connec); + MEDCouplingAutoRefCountObjectPtr newCoords=DataArrayDouble::New(); + newCoords->alloc(nbNodes+((int)addCoo.size())/2,2); + double *work=std::copy(_coords->begin(),_coords->end(),newCoords->getPointer()); + std::copy(addCoo.begin(),addCoo.end(),work); + DataArrayDouble::SetArrayIn(newCoords,_coords); + updateTime(); +} + /*! * This methods modify this by converting each cells into simplex cell, that is too say triangle for meshdim==2 or tetra for meshdim==3. * This cut into simplex is performed following the parameter 'policy'. This method so typically increases the number of cells of this. @@ -3318,6 +3403,87 @@ DataArrayInt *MEDCouplingUMesh::simplexizePol1() throw(INTERP_KERNEL::Exception) return ret; } +/*! + * This private method is used to subdivide edges of a mesh with meshdim==2. If 'this' has no a meshdim equal to 2 an exception will be thrown. + * This method completly ignore coordinates. + * @param nodeSubdived is the nodal connectivity of subdivision of edges + * @param nodeIndxSubdived is the nodal connectivity index of subdivision of edges + * @param desc is descending connectivity in format specified in MEDCouplingUMesh::buildDescendingConnectivity2 + * @param descIndex is descending connectivity index in format specified in MEDCouplingUMesh::buildDescendingConnectivity2 + */ +void MEDCouplingUMesh::subDivide2DMesh(const int *nodeSubdived, const int *nodeIndxSubdived, const int *desc, const int *descIndex) throw(INTERP_KERNEL::Exception) +{ + checkFullyDefined(); + if(getMeshDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : works only on umesh with meshdim==2 !"); + int nbOfCells=getNumberOfCells(); + int *connI=_nodal_connec_index->getPointer(); + int newConnLgth=0; + for(int i=0;i0; + int eedgeId=std::abs(desc[offset+nbOfEdges-1])-1; + int ref=ddirect?nodeSubdived[nodeIndxSubdived[eedgeId+1]-1]:nodeSubdived[nodeIndxSubdived[eedgeId]+1]; + for(int j=0;j0; + int edgeId=std::abs(desc[offset+j])-1; + if(!INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)nodeSubdived[nodeIndxSubdived[edgeId]]).isQuadratic()) + { + int id1=nodeSubdived[nodeIndxSubdived[edgeId]+1]; + int id2=nodeSubdived[nodeIndxSubdived[edgeId+1]-1]; + int ref2=direct?id1:id2; + if(ref==ref2) + { + int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1; + newConnLgth+=nbOfSubNodes-1; + ref=direct?id2:id1; + } + else + { + std::ostringstream oss; oss << "MEDCouplingUMesh::subDivide2DMesh : On polygon #" << i << " edgeid #" << j << " subedges mismatch : end subedge k!=start subedge k+1 !"; + throw INTERP_KERNEL::Exception(oss.str().c_str()); + } + } + else + { + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::subDivide2DMesh : this method only subdivides into linear edges !"); + } + } + newConnLgth++;//+1 is for cell type + connI[1]=newConnLgth; + } + // + MEDCouplingAutoRefCountObjectPtr newConn=DataArrayInt::New(); + newConn->alloc(newConnLgth,1); + int *work=newConn->getPointer(); + for(int i=0;i0; + int edgeId=std::abs(desc[offset+j])-1; + if(direct) + work=std::copy(nodeSubdived+nodeIndxSubdived[edgeId]+1,nodeSubdived+nodeIndxSubdived[edgeId+1]-1,work); + else + { + int nbOfSubNodes=nodeIndxSubdived[edgeId+1]-nodeIndxSubdived[edgeId]-1; + std::reverse_iterator it(nodeSubdived+nodeIndxSubdived[edgeId+1]); + work=std::copy(it,it+nbOfSubNodes-1,work); + } + } + } + DataArrayInt::SetArrayIn(newConn,_nodal_connec); + _types.clear(); + if(nbOfCells>0) + _types.insert(INTERP_KERNEL::NORM_POLYGON); +} /*! * This method converts all degenerated cells to simpler cells. For example a NORM_QUAD4 cell consituted from 2 same node id in its @@ -4784,7 +4950,7 @@ void MEDCouplingUMesh::FillInCompact3DMode(int spaceDim, int nbOfNodesInCell, co void MEDCouplingUMesh::writeVTKLL(std::ostream& ofs, const std::string& cellData, const std::string& pointData) const throw(INTERP_KERNEL::Exception) { - static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,-1,23,-1,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,-1,-1,-1,25,42,-1}; + static const int PARAMEDMEM2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,-1,23,-1,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,-1,-1,-1,25,42,-1,4}; ofs << " <" << getVTKDataSetType() << ">\n"; ofs << " \n"; ofs << " \n" << pointData << std::endl; @@ -4826,6 +4992,10 @@ std::string MEDCouplingUMesh::getVTKDataSetType() const throw(INTERP_KERNEL::Exc MEDCouplingUMesh *MEDCouplingUMesh::Intersect2DMeshes(const MEDCouplingUMesh *m1, const MEDCouplingUMesh *m2, double eps, DataArrayInt *&cellNb1, DataArrayInt *&cellNb2) throw(INTERP_KERNEL::Exception) { + m1->checkFullyDefined(); + m2->checkFullyDefined(); + if(m1->getMeshDimension()!=2 || m1->getSpaceDimension()!=2 || m2->getMeshDimension()!=2 || m2->getSpaceDimension()!=2) + throw INTERP_KERNEL::Exception("MEDCouplingUMesh::Intersect2DMeshes works on umeshes m1 AND m2 with meshdim equal to 2 and spaceDim equal to 2 too!"); std::vector< std::vector > intersectEdge1, colinear2, subDiv2; MEDCouplingUMesh *m1Desc=0,*m2Desc=0; DataArrayInt *desc1=0,*descIndx1=0,*revDesc1=0,*revDescIndx1=0,*desc2=0,*descIndx2=0,*revDesc2=0,*revDescIndx2=0; diff --git a/src/MEDCoupling/MEDCouplingUMesh.hxx b/src/MEDCoupling/MEDCouplingUMesh.hxx index 611f7bb14..d61087018 100644 --- a/src/MEDCoupling/MEDCouplingUMesh.hxx +++ b/src/MEDCoupling/MEDCouplingUMesh.hxx @@ -139,6 +139,8 @@ namespace ParaMEDMEM MEDCOUPLING_EXPORT bool isFullyQuadratic() const; MEDCOUPLING_EXPORT bool isPresenceOfQuadratic() const; MEDCOUPLING_EXPORT void convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void tessellate2D(double eps) throw(INTERP_KERNEL::Exception); + MEDCOUPLING_EXPORT void tessellate2DCurve(double eps) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT DataArrayInt *simplexize(int policy) throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT bool areOnlySimplexCells() const throw(INTERP_KERNEL::Exception); MEDCOUPLING_EXPORT void convertDegeneratedCells() throw(INTERP_KERNEL::Exception); @@ -194,6 +196,7 @@ namespace ParaMEDMEM //tools DataArrayInt *simplexizePol0() throw(INTERP_KERNEL::Exception); DataArrayInt *simplexizePol1() throw(INTERP_KERNEL::Exception); + void subDivide2DMesh(const int *nodeSubdived, const int *nodeIndxSubdived, const int *desc, const int *descIndex) throw(INTERP_KERNEL::Exception); void renumberNodesInConn(const int *newNodeNumbers); void fillCellIdsToKeepFromNodeIds(const int *begin, const int *end, bool fullyIn, std::vector& cellIdsKept) const; MEDCouplingUMesh *buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const; diff --git a/src/MEDCoupling_Swig/MEDCoupling.i b/src/MEDCoupling_Swig/MEDCoupling.i index 2c0366dac..16238d612 100644 --- a/src/MEDCoupling_Swig/MEDCoupling.i +++ b/src/MEDCoupling_Swig/MEDCoupling.i @@ -1040,6 +1040,8 @@ namespace ParaMEDMEM bool isPresenceOfQuadratic() const throw(INTERP_KERNEL::Exception); MEDCouplingFieldDouble *buildDirectionVectorField() const throw(INTERP_KERNEL::Exception); bool isContiguous1D() const throw(INTERP_KERNEL::Exception); + void tessellate2D(double eps) throw(INTERP_KERNEL::Exception); + void tessellate2DCurve(double eps) throw(INTERP_KERNEL::Exception); void convertQuadraticCellsToLinear() throw(INTERP_KERNEL::Exception); void convertDegeneratedCells() throw(INTERP_KERNEL::Exception); bool areOnlySimplexCells() const throw(INTERP_KERNEL::Exception); -- 2.39.2