1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License, or (at your option) any later version.
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
11 // Lesser General Public License for more details.
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #include "InterpKernelGeo2DEdgeLin.hxx"
22 #include "InterpKernelGeo2DNode.hxx"
23 #include "InterpKernelException.hxx"
24 #include "NormalizedUnstructuredMesh.hxx"
26 using namespace INTERP_KERNEL;
28 namespace INTERP_KERNEL
30 extern const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
33 SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):
34 SameTypeEdgeIntersector(e1,e2)
36 _matrix[0]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0];
37 _matrix[1]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1];
38 _matrix[2]=(*(e2.getEndNode()))[0]-(*(e2.getStartNode()))[0];
39 _matrix[3]=(*(e2.getEndNode()))[1]-(*(e2.getStartNode()))[1];
41 _determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
43 _col[0]=_matrix[1]*(*(e1.getStartNode()))[0]-_matrix[0]*(*(e1.getStartNode()))[1];
44 _col[1]=_matrix[3]*(*(e2.getStartNode()))[0]-_matrix[2]*(*(e2.getStartNode()))[1];
46 //Little trick to avoid problems if 'e1' and 'e2' are colinears and along Ox or Oy axes.
47 if(fabs(_matrix[1])>fabs(_matrix[0]))
54 * Must be called when 'this' and 'other' have been detected to be at least colinear. Typically they are overlapped.
56 bool SegSegIntersector::haveTheySameDirection() const
58 return (_matrix[0]*_matrix[2]+_matrix[1]*_matrix[3])>0.;
62 * Precondition start and end must be so that there predecessor was in the same direction than 'e1'
64 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
66 getCurveAbscisse(start,whereStart,commonNode);
67 getCurveAbscisse(end,whereEnd,commonNode);
70 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
73 obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
76 double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
86 * areColinears method should be called before with a returned colinearity equal to false to avoid bad news.
88 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
90 std::list< IntersectElement > ret;
93 // Intersection was already found: it is a common node shared by _e1 and _e2 - see areOverlappedOrOnlyColinears()
94 ret.push_back(*_earlyInter);
98 double x= (-_matrix[2]*_col[0]+_matrix[0]*_col[1]) / _determinant;
99 double y= (-_matrix[3]*_col[0]+_matrix[1]*_col[1]) / _determinant;
100 //Only one intersect point possible
101 Node *node=new Node(x,y);
103 bool i_1S=_e1.getStartNode()->isEqual(*node);
104 bool i_1E=_e1.getEndNode()->isEqual(*node);
105 bool i_2S=_e2.getStartNode()->isEqual(*node);
106 bool i_2E=_e2.getEndNode()->isEqual(*node);
107 ret.push_back(IntersectElement(_e1.getCharactValue(*node),
108 _e2.getCharactValue(*node),
109 i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder()));
114 * Retrieves if segs are colinears.
115 * Same philosophy as in other intersectors: we use epsilon as an absolute distance.
116 * If one puts the two vectors starting at the origin, determinant/dimChar is a close representative of the absolute distance between the tip of one vector
117 * to the other vector.
119 bool SegSegIntersector::areColinears() const
122 b1.prepareForAggregation();
123 b2.prepareForAggregation();
124 b1.aggregate(_e1.getBounds());
125 b2.aggregate(_e2.getBounds());
126 double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim());
128 // same criteria as in areOverlappedOrOnlyColinears, see comment below
129 return fabs(_determinant)<dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision();
133 * Should be called \b once ! non const method.
134 * \param whereToFind specifies the box where final seek should be done. Essentially it is used for caracteristic reason.
135 * \param colinearity returns if regarding QuadraticPlanarPrecision::getPrecision() ; e1 and e2 are colinears
136 * If true 'this' is modified ! So this method be called once above all if true is returned for this parameter.
137 * \param areOverlapped if colinearity if true, this parameter looks if e1 and e2 are overlapped, i.e. is they lie on the same line (= this is different from
138 * a true intersection, two segments can be in "overlap" mode, without intersecting)
140 void SegSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
143 b1.prepareForAggregation();
144 b2.prepareForAggregation();
145 b1.aggregate(_e1.getBounds());
146 b2.aggregate(_e2.getBounds());
147 double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim());
149 // Same criteria as in areColinears(), see doc.
150 if(fabs(_determinant)>dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision()) // Non colinear vectors
153 obviousNoIntersection=false;
155 // If they share one extremity, we can optimize since we already know where is the intersection:
157 identifyEarlyIntersection(a,b,c,d);
159 else // Colinear vectors
161 // Compute vectors joining tips of e1 and e2
162 double xS=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
163 double yS=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];
164 double xE=(*(_e1.getEndNode()))[0]-(*(_e2.getEndNode()))[0];
165 double yE=(*(_e1.getEndNode()))[1]-(*(_e2.getEndNode()))[1];
166 double maxDimS(std::max(fabs(xS),fabs(yS))), maxDimE(std::max(fabs(xE), fabs(yE)));
167 bool isS = (maxDimS > maxDimE), isE1 = (dimCharE1 >= dimCharE2);
168 double x = isS ? xS : xE;
169 double y = isS ? yS : yE;
170 unsigned shift = isE1 ? 0 : 2;
171 // test colinearity of the greatest tip-joining vector and greatest vector among {e1, e2}
172 areOverlapped = fabs(x*_matrix[1+shift]-y*_matrix[0+shift]) < dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision();
173 // explanation: if areOverlapped is true, we don't know yet if there will be an intersection (see meaning of areOverlapped in method doxy above)
174 // if areOverlapped is false, we have two colinear vectors, not lying on the same line, so we're sure there is no intersec
175 obviousNoIntersection = !areOverlapped;
179 EdgeLin::EdgeLin(std::istream& lineInXfig)
181 char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
182 lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
183 _start=new Node(lineInXfig);
184 _end=new Node(lineInXfig);
188 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
193 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
203 * Characteristic for edges is relative position btw 0.;1.
205 bool EdgeLin::isIn(double characterVal) const
207 return characterVal>0. && characterVal<1.;
210 Node *EdgeLin::buildRepresentantOfMySelf() const
212 return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
215 double EdgeLin::getCharactValue(const Node& node) const
217 return getCharactValueEng(node);
220 double EdgeLin::getCharactValueBtw0And1(const Node& node) const
222 return getCharactValueEng(node);
225 double EdgeLin::getDistanceToPoint(const double *pt) const
227 double loc=getCharactValueEng(pt);
231 tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
232 tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
233 return Node::distanceBtw2Pt(pt,tmp);
237 double dist1=Node::distanceBtw2Pt(*_start,pt);
238 double dist2=Node::distanceBtw2Pt(*_end,pt);
239 return std::min(dist1,dist2);
243 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
245 double dBase=sqrt(_start->distanceWithSq(*_end));
246 double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
247 d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
248 return Node::areDoubleEquals(dBase,d1);
251 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
253 stream << "2 1 0 1 ";
254 fillXfigStreamForLoc(stream);
255 stream << " 7 50 -1 -1 0.000 0 0 -1 1 0 2" << std::endl << "1 1 1.00 60.00 120.00" << std::endl;
256 direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
257 direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
261 void EdgeLin::update(Node *m)
266 double EdgeLin::getNormSq() const
268 return _start->distanceWithSq(*_end);
272 * This methods computes :
274 * \int_{Current Edge} -ydx
277 double EdgeLin::getAreaOfZone() const
279 return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
282 void EdgeLin::getBarycenter(double *bary) const
284 bary[0]=((*_start)[0]+(*_end)[0])/2.;
285 bary[1]=((*_start)[1]+(*_end)[1])/2.;
290 * bary[0]=\int_{Current Edge} -yxdx
293 * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
295 * To compute these 2 expressions in this class we have :
297 * y=y_{1}+\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1})
300 void EdgeLin::getBarycenterOfZone(double *bary) const
302 double x1=(*_start)[0];
303 double y1=(*_start)[1];
304 double x2=(*_end)[0];
305 double y2=(*_end)[1];
306 bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
307 //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
308 //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
309 bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
313 * Here \a this is not used (contrary to EdgeArcCircle class).
315 void EdgeLin::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
317 mid[0]=(p1[0]+p2[0])/2.;
318 mid[1]=(p1[1]+p2[1])/2.;
321 double EdgeLin::getCurveLength() const
323 double x=(*_start)[0]-(*_end)[0];
324 double y=(*_start)[1]-(*_end)[1];
325 return sqrt(x*x+y*y);
328 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
330 return new EdgeLin(start,end,direction);
334 * No precision should be introduced here. Just think as if precision was perfect.
336 void EdgeLin::updateBounds()
338 _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]));
341 double EdgeLin::getCharactValueEng(const double *node) const
343 double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
344 double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
345 return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);