src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx

   1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "InterpKernelGeo2DEdgeLin.hxx"
  22 #include "InterpKernelGeo2DNode.hxx"
  23 #include "InterpKernelException.hxx"
  24 #include "NormalizedUnstructuredMesh.hxx"
  25
  26 using namespace INTERP_KERNEL;
  27
  28 namespace INTERP_KERNEL
  29 {
  30   extern const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
  31 }
  32
  33 SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):SameTypeEdgeIntersector(e1,e2)
  34 {
  35 //  _matrix[0]=(*(e2.getStartNode()))[0]-(*(e2.getEndNode()))[0];
  36 //  _matrix[1]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0];
  37 //  _matrix[2]=(*(e2.getStartNode()))[1]-(*(e2.getEndNode()))[1];
  38 //  _matrix[3]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1];
  39
  40   _matrix[0]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0];
  41   _matrix[1]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1];
  42   _matrix[2]=(*(e2.getEndNode()))[0]-(*(e2.getStartNode()))[0];
  43   _matrix[3]=(*(e2.getEndNode()))[1]-(*(e2.getStartNode()))[1];
  44
  45
  46 //  _col[0]=_matrix[3]*(*(e1.getStartNode()))[0]-_matrix[1]*(*(e1.getStartNode()))[1];
  47 //  _col[1]=-_matrix[2]*(*(e2.getStartNode()))[0]+_matrix[0]*(*(e2.getStartNode()))[1];
  48   _col[0]=_matrix[1]*(*(e1.getStartNode()))[0]-_matrix[0]*(*(e1.getStartNode()))[1];
  49   _col[1]=_matrix[3]*(*(e2.getStartNode()))[0]-_matrix[2]*(*(e2.getStartNode()))[1];
  50
  51   //Little trick to avoid problems if 'e1' and 'e2' are colinears and along Ox or Oy axes.
  52   if(fabs(_matrix[1])>fabs(_matrix[0]))
  53     _ind=0;
  54   else
  55     _ind=1;
  56 }
  57
  58 /*!
  59  * Must be called when 'this' and 'other' have been detected to be at least colinear. Typically they are overlapped.
  60  */
  61 bool SegSegIntersector::haveTheySameDirection() const
  62 {
  63   return (_matrix[0]*_matrix[2]+_matrix[1]*_matrix[3])>0.;
  64 }
  65
  66 /*!
  67  * Precondition start and end must be so that there predecessor was in the same direction than 'e1'
  68  */
  69 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  70 {
  71   getCurveAbscisse(start,whereStart,commonNode);
  72   getCurveAbscisse(end,whereEnd,commonNode);
  73 }
  74
  75 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
  76 {
  77   bool obvious;
  78   obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
  79   if(obvious)
  80     return ;
  81   double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
  82   if(ret>0. && ret <1.)
  83     where=INSIDE;
  84   else if(ret<0.)
  85     where=OUT_BEFORE;
  86   else
  87     where=OUT_AFTER;
  88 }
  89
  90 /*!
  91  * areColinears method should be called before with a returned colinearity equal to false to avoid bad news.
  92  */
  93 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
  94 {
  95   std::list< IntersectElement > ret;
  96   double x=-_matrix[2]*_col[0]+_matrix[0]*_col[1];
  97   double y=-_matrix[3]*_col[0]+_matrix[1]*_col[1];
  98   //Only one intersect point possible
  99   Node *node=new Node(x,y);
 100   node->declareOn();
 101   bool i_1S=_e1.getStartNode()->isEqual(*node);
 102   bool i_1E=_e1.getEndNode()->isEqual(*node);
 103   bool i_2S=_e2.getStartNode()->isEqual(*node);
 104   bool i_2E=_e2.getEndNode()->isEqual(*node);
 105   ret.push_back(IntersectElement(_e1.getCharactValue(*node),
 106       _e2.getCharactValue(*node),
 107       i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder()));
 108   return ret;
 109 }
 110
 111 /*!
 112  * Retrieves if segs are colinears.
 113  * Same philosophy as in other intersectors: we use epsilon as an absolute distance.
 114  * 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
 115  * to the other vector.
 116  */
 117 bool SegSegIntersector::areColinears() const
 118 {
 119   Bounds b;
 120   b.prepareForAggregation();
 121   b.aggregate(_e1.getBounds());
 122   b.aggregate(_e2.getBounds());
 123   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 124   double dimChar=b.getCaracteristicDim();
 125
 126   return fabs(determinant)< 2.*dimChar*QuadraticPlanarPrecision::getPrecision(); // same criteria as in areOverlappedOrOnlyColinears, see comment below
 127 }
 128
 129 /*!
 130  * Should be called \b once ! non const method.
 131  * \param whereToFind specifies the box where final seek should be done. Essentially it is used for caracteristic reason.
 132  * \param colinearity returns if regarding QuadraticPlanarPrecision::getPrecision() ; e1 and e2 are colinears
 133  *                    If true 'this' is modified ! So this method be called once above all if true is returned for this parameter.
 134  * \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
 135  * a true intersection, two segments can be in "overlap" mode, without intersecting)
 136  */
 137 void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& obviousNoIntersection, bool& areOverlapped)
 138 {
 139   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 140   Bounds b;
 141   b.prepareForAggregation();
 142   b.aggregate(_e1.getBounds());
 143   b.aggregate(_e2.getBounds());
 144   double dimChar=b.getCaracteristicDim();
 145
 146   // Same criteria as in areColinears(), see doc.
 147   // [ABN] the 2 is not really justified, but the initial tests from Tony were written so closely to precision that I can't bother to change all of them ...
 148   if(fabs(determinant)>2.*dimChar*QuadraticPlanarPrecision::getPrecision())
 149     {
 150       obviousNoIntersection=false; areOverlapped=false;
 151       _matrix[0]/=determinant; _matrix[1]/=determinant; _matrix[2]/=determinant; _matrix[3]/=determinant;
 152     }
 153   else  // colinear vectors
 154     {
 155       double x=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
 156       double y=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];   // (x,y) is the vector between the two start points of e1 and e2
 157       areOverlapped = fabs(-_matrix[0]*y+_matrix[1]*x) < dimChar*QuadraticPlanarPrecision::getPrecision(); // test colinearity of (x,y) with e1
 158       // explanation: if areOverlapped is true, we don't know yet if there will be an intersection (see meaning of areOverlapped in method doxy above)
 159       // if areOverlapped is false, we have two colinear vectors, not lying on the same line, so we're sure there is no intersec
 160       obviousNoIntersection = !areOverlapped;
 161     }
 162 }
 163
 164 EdgeLin::EdgeLin(std::istream& lineInXfig)
 165 {
 166   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
 167   lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
 168   _start=new Node(lineInXfig);
 169   _end=new Node(lineInXfig);
 170   updateBounds();
 171 }
 172
 173 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
 174 {
 175   updateBounds();
 176 }
 177
 178 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
 179 {
 180   updateBounds();
 181 }
 182
 183 EdgeLin::~EdgeLin()
 184 {
 185 }
 186
 187 /*!
 188  * Characteristic for edges is relative position btw 0.;1.
 189  */
 190 bool EdgeLin::isIn(double characterVal) const
 191 {
 192   return characterVal>0. && characterVal<1.;
 193 }
 194
 195 Node *EdgeLin::buildRepresentantOfMySelf() const
 196 {
 197   return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
 198 }
 199
 200 double EdgeLin::getCharactValue(const Node& node) const
 201 {
 202   return getCharactValueEng(node);
 203 }
 204
 205 double EdgeLin::getCharactValueBtw0And1(const Node& node) const
 206 {
 207   return getCharactValueEng(node);
 208 }
 209
 210 double EdgeLin::getDistanceToPoint(const double *pt) const
 211 {
 212   double loc=getCharactValueEng(pt);
 213   if(loc>0. && loc<1.)
 214     {
 215       double tmp[2];
 216       tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
 217       tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
 218       return Node::distanceBtw2Pt(pt,tmp);
 219     }
 220   else
 221     {
 222       double dist1=Node::distanceBtw2Pt(*_start,pt);
 223       double dist2=Node::distanceBtw2Pt(*_end,pt);
 224       return std::min(dist1,dist2);
 225     }
 226 }
 227
 228 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
 229 {
 230   double dBase=sqrt(_start->distanceWithSq(*_end));
 231   double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
 232   d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
 233   return Node::areDoubleEquals(dBase,d1);
 234 }
 235
 236 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 237 {
 238   stream << "2 1 0 1 ";
 239   fillXfigStreamForLoc(stream);
 240   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;
 241   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 242   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 243   stream << std::endl;
 244 }
 245
 246 void EdgeLin::update(Node *m)
 247 {
 248   updateBounds();
 249 }
 250
 251 double EdgeLin::getNormSq() const
 252 {
 253   return _start->distanceWithSq(*_end);
 254 }
 255
 256 /*!
 257  * This methods computes :
 258  * \f[
 259  * \int_{Current Edge} -ydx
 260  * \f]
 261  */
 262 double EdgeLin::getAreaOfZone() const
 263 {
 264   return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 265 }
 266
 267 void EdgeLin::getBarycenter(double *bary) const
 268 {
 269   bary[0]=((*_start)[0]+(*_end)[0])/2.;
 270   bary[1]=((*_start)[1]+(*_end)[1])/2.;
 271 }
 272
 273 /*!
 274  * \f[
 275  * bary[0]=\int_{Current Edge} -yxdx
 276  * \f]
 277  * \f[
 278  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 279  * \f]
 280  * To compute these 2 expressions in this class we have :
 281  * \f[
 282  * y=y_{1}+\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1})
 283  * \f]
 284  */
 285 void EdgeLin::getBarycenterOfZone(double *bary) const
 286 {
 287   double x1=(*_start)[0];
 288   double y1=(*_start)[1];
 289   double x2=(*_end)[0];
 290   double y2=(*_end)[1];
 291   bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
 292   //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
 293   //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
 294   bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
 295 }
 296
 297 /*!
 298  * Here \a this is not used (contrary to EdgeArcCircle class).
 299  */
 300 void EdgeLin::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
 301 {
 302   mid[0]=(p1[0]+p2[0])/2.;
 303   mid[1]=(p1[1]+p2[1])/2.;
 304 }
 305
 306 double EdgeLin::getCurveLength() const
 307 {
 308   double x=(*_start)[0]-(*_end)[0];
 309   double y=(*_start)[1]-(*_end)[1];
 310   return sqrt(x*x+y*y);
 311 }
 312
 313 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 314 {
 315   return new EdgeLin(start,end,direction);
 316 }
 317
 318 /*!
 319  * No precision should be introduced here. Just think as if precision was perfect.
 320  */
 321 void EdgeLin::updateBounds()
 322 {
 323   _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]));
 324 }
 325
 326 double EdgeLin::getCharactValueEng(const double *node) const
 327 {
 328   double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
 329   double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
 330   return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);
 331 }