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);
}
}
}
{
if(obvious1)
{
- if(EdgeArcCircle::isIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
+ if(EdgeArcCircle::IsIn2Pi(getE1().getAngle0(),getE1().getAngle(),angleInRadEnd))
whereEnd=INSIDE;
else
whereEnd=OUT_AFTER;
}
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;
*/
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)
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.);
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
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);
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);
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);
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();
_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();
}
void EdgeArcCircle::changeMiddle(Node *newMiddle)
{
- getArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
+ GetArcOfCirclePassingThru(*_start,*newMiddle,*_end,_center,_radius,_angle,_angle0);
updateBounds();
}
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.)
_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<int>& newConn, std::vector<double>& 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;i<nbOfSubDiv;i++,offset2++)
+ {
+ double angle=_angle0+i*eps*signOfAngle;
+ newConn.push_back(offset2);
+ addCoo.push_back(_center[0]+_radius*cos(angle)); addCoo.push_back(_center[1]+_radius*sin(angle));
+ }
+ newConn.push_back(conn[1]);
+}
+
+EdgeArcCircle *EdgeArcCircle::BuildFromNodes(Node *start, Node *middle, Node *end)
+{
+ EdgeLin *e1,*e2;
+ e1=new EdgeLin(start,middle);
+ e2=new EdgeLin(middle,end);
+ SegSegIntersector inters(*e1,*e2);
+ bool colinearity=inters.areColinears();
+ delete e1; delete e2;
+ if(colinearity)
+ {
+ start->decrRef(); 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);
}
/*!
* 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;
}
else
{
- double ret=safeAsin(uy);
+ double ret=SafeAsin(uy);
if(ux>0.)
return ret;
if(ret>0.)
}
}
-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]);
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;
}
void EdgeArcCircle::update(Node *m)
{
- getArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
+ GetArcOfCirclePassingThru(*_start,*m,*_end,_center,_radius,_angle,_angle0);
updateBounds();
}
*/
bool EdgeArcCircle::isIn(double characterVal) const
{
- return isIn2Pi(_angle0,_angle,characterVal);
+ return IsIn2Pi(_angle0,_angle,characterVal);
}
Node *EdgeArcCircle::buildRepresentantOfMySelf() 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
{
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.)
/*!
* 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.)
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;
}
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<int>& newConn, std::vector<double>& 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;
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
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);
}
/*!
{
}
-QuadraticPolygon *QuadraticPolygon::buildLinearPolygon(std::vector<Node *>& nodes)
+QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
{
QuadraticPolygon *ret=new QuadraticPolygon;
std::size_t size=nodes.size();
return ret;
}
-QuadraticPolygon *QuadraticPolygon::buildArcCirclePolygon(std::vector<Node *>& nodes)
+QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
{
QuadraticPolygon *ret=new QuadraticPolygon;
std::size_t size=nodes.size();
return ret;
}
-void QuadraticPolygon::buildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
+void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
{
std::ofstream file(fileName);
file << std::setprecision(16);
{
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
}
else
{
- updateNeighbours(merge,it1,it2,curE1,curE2);
+ UpdateNeighbours(merge,it1,it2,curE1,curE2);
it1.next();
}
}
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);
isColinear = false;
perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
QuadraticPolygon cpyOfThis(*this);
QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
- splitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
+ SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
performLocatingOperation(cpyOfOther);
other.performLocatingOperation(cpyOfThis);
cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
QuadraticPolygon tmp;
tmp.pushBack(curE1->clone());
int tmp2;
- splitPolygonsEachOther(tmp,cpyOfOther,tmp2);
+ SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
other.performLocatingOperation(tmp);
tmp.dispatchPerimeter(polThis[edgeId]);
}
QuadraticPolygon tmp;
tmp.pushBack(curE2->clone());
int tmp2;
- splitPolygonsEachOther(tmp,cpyOfThis,tmp2);
+ SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
performLocatingOperation(tmp);
tmp.dispatchPerimeter(polOther[edgeId]);
}
QuadraticPolygon tmp;
tmp.pushBack(curE1->clone());
int tmp2;
- splitPolygonsEachOther(tmp,cpyOfOther,tmp2);
+ SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
}
}
{
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);
* 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;
{
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
}
else
{
- updateNeighbours(merge,it1,it2,curE1,curE2);
+ UpdateNeighbours(merge,it1,it2,curE1,curE2);
it1.next();
}
}
pushBack(tmp);
nodeToTest=tmp->getEndNode();
direction?it.nextLoop():it.previousLoop();
- ret=checkInList(nodeToTest,iStart,iEnd);
+ ret=CheckInList(nodeToTest,iStart,iEnd);
if(completed())
return iEnd;
}
return ret;
}
-std::list<QuadraticPolygon *>::iterator QuadraticPolygon::checkInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
+std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
std::list<QuadraticPolygon *>::iterator iEnd)
{
for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
QuadraticPolygon() { }
QuadraticPolygon(const QuadraticPolygon& other):ComposedEdge(other) { }
QuadraticPolygon(const char *fileName);
- static QuadraticPolygon *buildLinearPolygon(std::vector<Node *>& nodes);
- static QuadraticPolygon *buildArcCirclePolygon(std::vector<Node *>& nodes);
- static void buildDbgFile(const std::vector<Node *>& nodes, const char *fileName);
+ static QuadraticPolygon *BuildLinearPolygon(std::vector<Node *>& nodes);
+ static QuadraticPolygon *BuildArcCirclePolygon(std::vector<Node *>& nodes);
+ static void BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName);
~QuadraticPolygon();
void closeMe() const;
void circularPermute();
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<QuadraticPolygon *> buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const;
bool amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted, const QuadraticPolygon& pol2NotSplitted, bool& direction);
protected:
void dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const;
void closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1, std::vector<QuadraticPolygon *>& results) const;
template<class EDGES>
- 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<QuadraticPolygon *>::iterator fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
std::list<QuadraticPolygon *>::iterator iStart,
std::list<QuadraticPolygon *>::iterator iEnd,
bool direction);
- static std::list<QuadraticPolygon *>::iterator checkInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
+ static std::list<QuadraticPolygon *>::iterator CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
std::list<QuadraticPolygon *>::iterator iEnd);
};
}
namespace INTERP_KERNEL
{
template<class EDGES>
- 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();
std::vector<Node *> nodes2(nbOfSourceNodes);
for(int i=0;i<nbOfSourceNodes;i++)
nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
- QuadraticPolygon *p1=QuadraticPolygon::buildLinearPolygon(nodes);
+ QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
QuadraticPolygon *p2;
if(!isSourceQuad)
- p2=QuadraticPolygon::buildLinearPolygon(nodes2);
+ p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
else
- p2=QuadraticPolygon::buildArcCirclePolygon(nodes2);
+ p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
double ret=p1->intersectWithAbs(*p2);
delete p1; delete p2;
return ret;
std::vector<Node *> nodes2(nbOfSourceNodes);
for(int i=0;i<nbOfSourceNodes;i++)
nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
- QuadraticPolygon *p1=QuadraticPolygon::buildLinearPolygon(nodes);
- QuadraticPolygon *p2=QuadraticPolygon::buildLinearPolygon(nodes2);
+ QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
+ QuadraticPolygon *p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
double ret=p1->intersectWithAbs(*p2);
delete p1; delete p2;
return ret;
std::vector<Node *> nodes2(nbOfTargetNodes);
for(int i=0;i<nbOfTargetNodes;i++)
nodes2[i]=new Node(targetCell[i*SPACEDIM],targetCell[i*SPACEDIM+1]);
- QuadraticPolygon *p1=QuadraticPolygon::buildLinearPolygon(nodes);
+ QuadraticPolygon *p1=QuadraticPolygon::BuildLinearPolygon(nodes);
QuadraticPolygon *p2;
if(!targetCellQuadratic)
- p2=QuadraticPolygon::buildLinearPolygon(nodes2);
+ p2=QuadraticPolygon::BuildLinearPolygon(nodes2);
else
- p2=QuadraticPolygon::buildArcCirclePolygon(nodes2);
+ p2=QuadraticPolygon::BuildArcCirclePolygon(nodes2);
double barycenter[2];
double ret=p1->intersectWithAbs(*p2,barycenter);
delete p1; delete p2;
for(int i=0;i<nbNodes;i++)
nodes[i]=new Node(coords[i*SPACEDIM],coords[i*SPACEDIM+1]);
if(!CellModel::GetCellModel(type).isQuadratic())
- return QuadraticPolygon::buildLinearPolygon(nodes);
+ return QuadraticPolygon::BuildLinearPolygon(nodes);
else
- return QuadraticPolygon::buildArcCirclePolygon(nodes);
+ return QuadraticPolygon::BuildArcCirclePolygon(nodes);
}
INTERSECTOR_TEMPLATE
for(int i=0;i<nbOfPoints;i++)
nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsT+OTT<ConnType,numPol>::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
for(int i=0;i<nbOfPoints;i++)
nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsS+OTT<ConnType,numPol>::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);
}
}
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;
std::vector<Node *> nodes2(nbOfSourceNodes);
for(int i=0;i<nbOfSourceNodes;i++)
nodes2[i]=new Node(sourceCoords[i*SPACEDIM],sourceCoords[i*SPACEDIM+1]);
- QuadraticPolygon *p=QuadraticPolygon::buildLinearPolygon(nodes2);
+ QuadraticPolygon *p=QuadraticPolygon::BuildLinearPolygon(nodes2);
double bary[SPACEDIM];
p->getBarycenterGeneral(bary);
delete p;
for(int i=0;i<nbNodes;i++)
nodes[i]=new Node(coords[i*SPACEDIM],coords[i*SPACEDIM+1]);
if(!CellModel::GetCellModel(type).isQuadratic())
- return QuadraticPolygon::buildLinearPolygon(nodes);
+ return QuadraticPolygon::BuildLinearPolygon(nodes);
else
- return QuadraticPolygon::buildArcCirclePolygon(nodes);
+ return QuadraticPolygon::BuildArcCirclePolygon(nodes);
}
INTERSECTOR_TEMPLATE
for(int i=0;i<nbOfPoints;i++)
nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsT+OTT<ConnType,numPol>::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
for(int i=0;i<nbOfPoints;i++)
nodes[i]=new Node(PlanarIntersector<MyMeshType,MyMatrix>::_coordsS+OTT<ConnType,numPol>::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);
}
}
}
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;
}
}
}
-
- /*void MEDCouplingUMeshAssignOnLoc(const INTERP_KERNEL::QuadraticPolygon& pol1, INTERP_KERNEL::QuadraticPolygon& pol2,
- const int *desc1Bg, const int *desc1End,const std::vector<std::vector<int> >& intesctEdges1,
- const int *desc2Bg, const int *desc2End,const std::vector<std::vector<int> >& intesctEdges2,
- const std::vector<std::vector<int> >& colinear1)
- {
- for(const int *desc1=desc1Bg;desc1!=desc1End;desc1++)
- {
- int eltId1=abs(*desc1)-1;
- for(std::vector<int>::const_iterator it1=intesctEdges1[eltId1].begin();it1!=intesctEdges1[eltId1].end();it1++)
- {
- std::map<int,INTERP_KERNEL::Node *>::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
const double *p0=i+1<nbOfLevsInVec?begin:third;
const double *p1=i+1<nbOfLevsInVec?end:begin;
const double *p2=i+1<nbOfLevsInVec?third:end;
- INTERP_KERNEL::EdgeArcCircle::getArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
+ INTERP_KERNEL::EdgeArcCircle::GetArcOfCirclePassingThru(p0,p1,p2,tmp3,radius,alpha,alpha0);
double cosangle=i+1<nbOfLevsInVec?(p0[0]-tmp3[0])*(p1[0]-tmp3[0])+(p0[1]-tmp3[1])*(p1[1]-tmp3[1]):(p2[0]-tmp3[0])*(p1[0]-tmp3[0])+(p2[1]-tmp3[1])*(p1[1]-tmp3[1]);
double angle=acos(cosangle/(radius*radius));
tmp->rotate(end,0,angle);
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+1<nbOfLevsInVec?(p0r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p0r[1]-tmp3[1])*(p1r[1]-tmp3[1]):(p2r[0]-tmp3[0])*(p1r[0]-tmp3[0])+(p2r[1]-tmp3[1])*(p1r[1]-tmp3[1]);
double angle=acos(cosangle/(radius*radius));
tmp->rotate(end,vecPlane,angle);
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<std::numeric_limits<double>::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<DataArrayInt> desc1=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> descIndx1=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDesc1=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> revDescIndx1=DataArrayInt::New();
+ MEDCouplingAutoRefCountObjectPtr<MEDCouplingUMesh> 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<std::numeric_limits<double>::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<double> addCoo;
+ std::vector<int> newConn;
+ MEDCouplingAutoRefCountObjectPtr<DataArrayInt> newConnI=DataArrayInt::New();
+ newConnI->alloc(nbCells+1,1);
+ int *newConnIPtr=newConnI->getPointer();
+ *newConnIPtr=0;
+ int tmp1[3];
+ INTERP_KERNEL::Node *tmp2[3];
+ std::set<INTERP_KERNEL::NormalizedCellType> types;
+ for(int i=0;i<nbCells;i++,newConnIPtr++)
+ {
+ const INTERP_KERNEL::CellModel& cm=INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[connI[i]]);
+ if(cm.isQuadratic())
+ {//assert(connI[i+1]-connI[i]-1==3)
+ tmp1[0]=conn[connI[i]+1+0]; tmp1[1]=conn[connI[i]+1+1]; tmp1[2]=conn[connI[i]+1+2];
+ tmp2[0]=new INTERP_KERNEL::Node(coords[2*tmp1[0]],coords[2*tmp1[0]+1]);
+ tmp2[1]=new INTERP_KERNEL::Node(coords[2*tmp1[1]],coords[2*tmp1[1]+1]);
+ tmp2[2]=new INTERP_KERNEL::Node(coords[2*tmp1[2]],coords[2*tmp1[2]+1]);
+ INTERP_KERNEL::EdgeArcCircle *eac=INTERP_KERNEL::EdgeArcCircle::BuildFromNodes(tmp2[0],tmp2[2],tmp2[1]);
+ if(eac)
+ {
+ eac->tesselate(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<DataArrayInt> newConnArr=DataArrayInt::New();
+ newConnArr->alloc((int)newConn.size(),1);
+ std::copy(newConn.begin(),newConn.end(),newConnArr->getPointer());
+ DataArrayInt::SetArrayIn(newConnArr,_nodal_connec);
+ MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> 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.
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;i<nbOfCells;i++,connI++)
+ {
+ int offset=descIndex[i];
+ int nbOfEdges=descIndex[i+1]-offset;
+ //
+ bool ddirect=desc[offset+nbOfEdges-1]>0;
+ 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;j<nbOfEdges;j++)
+ {
+ bool direct=desc[offset+j]>0;
+ 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<DataArrayInt> newConn=DataArrayInt::New();
+ newConn->alloc(newConnLgth,1);
+ int *work=newConn->getPointer();
+ for(int i=0;i<nbOfCells;i++)
+ {
+ *work++=INTERP_KERNEL::NORM_POLYGON;
+ int offset=descIndex[i];
+ int nbOfEdges=descIndex[i+1]-offset;
+ for(int j=0;j<nbOfEdges;j++)
+ {
+ bool direct=desc[offset+j]>0;
+ 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<const int *> 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
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 << " <Piece NumberOfPoints=\"" << getNumberOfNodes() << "\" NumberOfCells=\"" << getNumberOfCells() << "\">\n";
ofs << " <PointData>\n" << pointData << std::endl;
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<int> > 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;
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);
//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<int>& cellIdsKept) const;
MEDCouplingUMesh *buildExtrudedMeshFromThisLowLev(int nbOfNodesOf1Lev, bool isQuad) const;
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);
