src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdgeLin.cxx

   1 // Copyright (C) 2007-2020  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   const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
  31 }
  32
  33 SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):
  34         SameTypeEdgeIntersector(e1,e2)
  35 {
  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];
  40
  41   _determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
  42
  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];
  45
  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]))
  48     _ind=0;
  49   else
  50     _ind=1;
  51 }
  52
  53 /*!
  54  * Must be called when 'this' and 'other' have been detected to be at least colinear. Typically they are overlapped.
  55  */
  56 bool SegSegIntersector::haveTheySameDirection() const
  57 {
  58   return (_matrix[0]*_matrix[2]+_matrix[1]*_matrix[3])>0.;
  59 }
  60
  61 /*!
  62  * Precondition start and end must be so that there predecessor was in the same direction than 'e1'
  63  */
  64 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  65 {
  66   getCurveAbscisse(start,whereStart,commonNode);
  67   getCurveAbscisse(end,whereEnd,commonNode);
  68 }
  69
  70 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
  71 {
  72   bool obvious;
  73   obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
  74   if(obvious)
  75     return ;
  76   double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
  77   if(ret>0. && ret <1.)
  78     where=INSIDE;
  79   else if(ret<0.)
  80     where=OUT_BEFORE;
  81   else
  82     where=OUT_AFTER;
  83 }
  84
  85 /*!
  86  * areColinears method should be called before with a returned colinearity equal to false to avoid bad news.
  87  */
  88 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
  89 {
  90   std::list< IntersectElement > ret;
  91   if (_earlyInter)
  92     {
  93       // Intersection was already found: it is a common node shared by _e1 and _e2 - see areOverlappedOrOnlyColinears()
  94       ret.push_back(*_earlyInter);
  95       return ret;
  96     }
  97
  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);
 102   node->declareOn();
 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()));
 110   return ret;
 111 }
 112
 113 /*!
 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.
 118  */
 119 bool SegSegIntersector::areColinears() const
 120 {
 121   Bounds b1, b2;
 122   b1.prepareForAggregation();
 123   b2.prepareForAggregation();
 124   b1.aggregate(_e1.getBounds());
 125   b2.aggregate(_e2.getBounds());
 126   double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim());
 127
 128   // same criteria as in areOverlappedOrOnlyColinears, see comment below
 129   return fabs(_determinant)<dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision();
 130 }
 131
 132 /*!
 133  * Should be called \b once ! non const method.
 134  * \param obviousNoIntersection set to true if it is obvious that there is no intersection
 135  * \param areOverlapped if the two segs are colinears, this parameter looks if e1 and e2 are overlapped, i.e. is they lie on the same line (= this is different from
 136  * a true intersection, two segments can be in "overlap" mode, without intersecting)
 137  */
 138 void SegSegIntersector::areOverlappedOrOnlyColinears(bool& obviousNoIntersection, bool& areOverlapped)
 139 {
 140   Bounds b1, b2;
 141   b1.prepareForAggregation();
 142   b2.prepareForAggregation();
 143   b1.aggregate(_e1.getBounds());
 144   b2.aggregate(_e2.getBounds());
 145   double dimCharE1(b1.getCaracteristicDim()) ,dimCharE2(b2.getCaracteristicDim());
 146
 147   // Same criteria as in areColinears(), see doc.
 148   if(fabs(_determinant)>dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision())  // Non colinear vectors
 149     {
 150       areOverlapped=false;
 151       obviousNoIntersection=false;
 152
 153       // If they share one extremity, we can optimize since we already know where is the intersection:
 154       bool a,b,c,d;
 155       identifyEarlyIntersection(a,b,c,d);
 156     }
 157   else  // Colinear vectors
 158     {
 159       // Compute vectors joining tips of e1 and e2
 160       double xS=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
 161       double yS=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];
 162       double xE=(*(_e1.getEndNode()))[0]-(*(_e2.getEndNode()))[0];
 163       double yE=(*(_e1.getEndNode()))[1]-(*(_e2.getEndNode()))[1];
 164       double maxDimS(std::max(fabs(xS),fabs(yS))), maxDimE(std::max(fabs(xE), fabs(yE)));
 165       bool isS = (maxDimS > maxDimE), isE1 = (dimCharE1 >= dimCharE2);
 166       double x = isS ? xS : xE;
 167       double y = isS ? yS : yE;
 168       unsigned shift = isE1 ? 0 : 2;
 169       // test colinearity of the greatest tip-joining vector and greatest vector among {e1, e2}
 170       areOverlapped = fabs(x*_matrix[1+shift]-y*_matrix[0+shift]) < dimCharE1*dimCharE2*QuadraticPlanarPrecision::getPrecision();
 171       // explanation: if areOverlapped is true, we don't know yet if there will be an intersection (see meaning of areOverlapped in method doxy above)
 172       // if areOverlapped is false, we have two colinear vectors, not lying on the same line, so we're sure there is no intersec
 173       obviousNoIntersection = !areOverlapped;
 174     }
 175 }
 176
 177 EdgeLin::EdgeLin(std::istream& lineInXfig)
 178 {
 179   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
 180   lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
 181   _start=new Node(lineInXfig);
 182   _end=new Node(lineInXfig);
 183   updateBounds();
 184 }
 185
 186 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
 187 {
 188   updateBounds();
 189 }
 190
 191 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
 192 {
 193   updateBounds();
 194 }
 195
 196 EdgeLin::~EdgeLin()
 197 {
 198 }
 199
 200 /*!
 201  * Characteristic for edges is relative position btw 0.;1.
 202  */
 203 bool EdgeLin::isIn(double characterVal) const
 204 {
 205   return characterVal>0. && characterVal<1.;
 206 }
 207
 208 Node *EdgeLin::buildRepresentantOfMySelf() const
 209 {
 210   return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
 211 }
 212
 213 double EdgeLin::getCharactValue(const Node& node) const
 214 {
 215   return getCharactValueEng(node);
 216 }
 217
 218 double EdgeLin::getCharactValueBtw0And1(const Node& node) const
 219 {
 220   return getCharactValueEng(node);
 221 }
 222
 223 double EdgeLin::getDistanceToPoint(const double *pt) const
 224 {
 225   double loc=getCharactValueEng(pt);
 226   if(loc>0. && loc<1.)
 227     {
 228       double tmp[2];
 229       tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
 230       tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
 231       return Node::distanceBtw2Pt(pt,tmp);
 232     }
 233   else
 234     {
 235       double dist1=Node::distanceBtw2Pt(*_start,pt);
 236       double dist2=Node::distanceBtw2Pt(*_end,pt);
 237       return std::min(dist1,dist2);
 238     }
 239 }
 240
 241 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
 242 {
 243   double dBase=sqrt(_start->distanceWithSq(*_end));
 244   double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
 245   d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
 246   return Node::areDoubleEquals(dBase,d1);
 247 }
 248
 249 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 250 {
 251   stream << "2 1 0 1 ";
 252   fillXfigStreamForLoc(stream);
 253   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;
 254   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 255   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 256   stream << std::endl;
 257 }
 258
 259 void EdgeLin::update(Node *m)
 260 {
 261   updateBounds();
 262 }
 263
 264 double EdgeLin::getNormSq() const
 265 {
 266   return _start->distanceWithSq(*_end);
 267 }
 268
 269 /*!
 270  * This methods computes :
 271  * \f[
 272  * \int_{Current Edge} -ydx
 273  * \f]
 274  */
 275 double EdgeLin::getAreaOfZone() const
 276 {
 277   return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 278 }
 279
 280 void EdgeLin::getBarycenter(double *bary) const
 281 {
 282   bary[0]=((*_start)[0]+(*_end)[0])/2.;
 283   bary[1]=((*_start)[1]+(*_end)[1])/2.;
 284 }
 285
 286 /*!
 287  * \f[
 288  * bary[0]=\int_{Current Edge} -yxdx
 289  * \f]
 290  * \f[
 291  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 292  * \f]
 293  * To compute these 2 expressions in this class we have :
 294  * \f[
 295  * y=y_{1}+\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1})
 296  * \f]
 297  */
 298 void EdgeLin::getBarycenterOfZone(double *bary) const
 299 {
 300   double x1=(*_start)[0];
 301   double y1=(*_start)[1];
 302   double x2=(*_end)[0];
 303   double y2=(*_end)[1];
 304   bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
 305   //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
 306   //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
 307   bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
 308 }
 309
 310 /*!
 311  * Here \a this is not used (contrary to EdgeArcCircle class).
 312  */
 313 void EdgeLin::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
 314 {
 315   mid[0]=(p1[0]+p2[0])/2.;
 316   mid[1]=(p1[1]+p2[1])/2.;
 317 }
 318
 319 double EdgeLin::getCurveLength() const
 320 {
 321   double x=(*_start)[0]-(*_end)[0];
 322   double y=(*_start)[1]-(*_end)[1];
 323   return sqrt(x*x+y*y);
 324 }
 325
 326 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 327 {
 328   return new EdgeLin(start,end,direction);
 329 }
 330
 331 /*!
 332  * No precision should be introduced here. Just think as if precision was perfect.
 333  */
 334 void EdgeLin::updateBounds()
 335 {
 336   _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]));
 337 }
 338
 339 double EdgeLin::getCharactValueEng(const double *node) const
 340 {
 341   double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
 342   double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
 343   return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);
 344 }