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   _col[0]=_matrix[3]*(*(e1.getStartNode()))[0]-_matrix[1]*(*(e1.getStartNode()))[1];
  40   _col[1]=-_matrix[2]*(*(e2.getStartNode()))[0]+_matrix[0]*(*(e2.getStartNode()))[1];
  41   //Little trick to avoid problems if 'e1' and 'e2' are colinears and along Ox or Oy axes.
  42   if(fabs(_matrix[3])>fabs(_matrix[1]))
  43     _ind=0;
  44   else
  45     _ind=1;
  46 }
  47
  48 /*!
  49  * Must be called when 'this' and 'other' have been detected to be at least colinear. Typically they are overlapped.
  50  * Must be called after call of areOverlappedOrOnlyColinears.
  51  */
  52 bool SegSegIntersector::haveTheySameDirection() const
  53 {
  54   return (_matrix[3]*_matrix[1]+_matrix[2]*_matrix[0])>0.;
  55   //return (_matrix[_ind?1:0]>0. && _matrix[_ind?3:2]>0.) || (_matrix[_ind?1:0]<0. && _matrix[_ind?3:2]<0.);
  56 }
  57
  58 /*!
  59  * Precondition start and end must be so that there predecessor was in the same direction than 'e1'
  60  */
  61 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
  62 {
  63   getCurveAbscisse(start,whereStart,commonNode);
  64   getCurveAbscisse(end,whereEnd,commonNode);
  65 }
  66
  67 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
  68 {
  69   bool obvious;
  70   obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
  71   if(obvious)
  72     return ;
  73   double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
  74   if(ret>0. && ret <1.)
  75     where=INSIDE;
  76   else if(ret<0.)
  77     where=OUT_BEFORE;
  78   else
  79     where=OUT_AFTER;
  80 }
  81
  82 /*!
  83  * areColinears method should be called before with a returned colinearity equal to false to avoid bad news.
  84  */
  85 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
  86 {
  87   std::list< IntersectElement > ret;
  88   double x=_matrix[0]*_col[0]+_matrix[1]*_col[1];
  89   double y=_matrix[2]*_col[0]+_matrix[3]*_col[1];
  90   //Only one intersect point possible
  91   Node *node=new Node(x,y);
  92   node->declareOn();
  93   bool i_1S=_e1.getStartNode()->isEqual(*node);
  94   bool i_1E=_e1.getEndNode()->isEqual(*node);
  95   bool i_2S=_e2.getStartNode()->isEqual(*node);
  96   bool i_2E=_e2.getEndNode()->isEqual(*node);
  97   ret.push_back(IntersectElement(_e1.getCharactValue(*node),
  98       _e2.getCharactValue(*node),
  99       i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder()));
 100   return ret;
 101 }
 102
 103 /*!
 104  * retrieves if segs are colinears.
 105  * WARNING !!! Contrary to areOverlappedOrOnlyColinears method, this method use an
 106  * another precision to detect colinearity !
 107  */
 108 bool SegSegIntersector::areColinears() const
 109 {
 110   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 111   return fabs(determinant)<QUADRATIC_PLANAR::_arc_detection_precision;
 112 }
 113
 114 /*!
 115  * Should be called \b once ! non const method.
 116  * \param whereToFind specifies the box where final seek should be done. Essentially it is used for caracteristic reason.
 117  * \param colinearity returns if regarding QUADRATIC_PLANAR::_precision ; e1 and e2 are colinears
 118  *                    If true 'this' is modified ! So this method be called once above all if true is returned for this parameter.
 119  * \param areOverlapped if colinearity if true, this parameter looks if e1 and e2 are overlapped.
 120  */
 121 void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& colinearity, bool& areOverlapped)
 122 {
 123   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
 124   if(fabs(determinant)>2.*QUADRATIC_PLANAR::_precision)//2*_precision due to max of offset on _start and _end
 125     {
 126       colinearity=false; areOverlapped=false;
 127       _matrix[0]/=determinant; _matrix[1]/=determinant; _matrix[2]/=determinant; _matrix[3]/=determinant;
 128     }
 129   else
 130     {
 131       colinearity=true;
 132       //retrieving initial matrix
 133       double tmp=_matrix[0]; _matrix[0]=_matrix[3]; _matrix[3]=tmp;
 134       _matrix[1]=-_matrix[1]; _matrix[2]=-_matrix[2];
 135       //
 136       double deno=sqrt(_matrix[0]*_matrix[0]+_matrix[1]*_matrix[1]);
 137       double x=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
 138       double y=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];
 139       areOverlapped=fabs((_matrix[1]*y+_matrix[0]*x)/deno)<QUADRATIC_PLANAR::_precision;
 140     }
 141 }
 142
 143 EdgeLin::EdgeLin(std::istream& lineInXfig)
 144 {
 145   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
 146   lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
 147   _start=new Node(lineInXfig);
 148   _end=new Node(lineInXfig);
 149   updateBounds();
 150 }
 151
 152 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
 153 {
 154   updateBounds();
 155 }
 156
 157 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
 158 {
 159   updateBounds();
 160 }
 161
 162 EdgeLin::~EdgeLin()
 163 {
 164 }
 165
 166 /*!
 167  * Characteristic for edges is relative position btw 0.;1.
 168  */
 169 bool EdgeLin::isIn(double characterVal) const
 170 {
 171   return characterVal>0. && characterVal<1.;
 172 }
 173
 174 Node *EdgeLin::buildRepresentantOfMySelf() const
 175 {
 176   return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
 177 }
 178
 179 double EdgeLin::getCharactValue(const Node& node) const
 180 {
 181   return getCharactValueEng(node);
 182 }
 183
 184 double EdgeLin::getCharactValueBtw0And1(const Node& node) const
 185 {
 186   return getCharactValueEng(node);
 187 }
 188
 189 double EdgeLin::getDistanceToPoint(const double *pt) const
 190 {
 191   double loc=getCharactValueEng(pt);
 192   if(loc>0. && loc<1.)
 193     {
 194       double tmp[2];
 195       tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
 196       tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
 197       return Node::distanceBtw2Pt(pt,tmp);
 198     }
 199   else
 200     {
 201       double dist1=Node::distanceBtw2Pt(*_start,pt);
 202       double dist2=Node::distanceBtw2Pt(*_end,pt);
 203       return std::min(dist1,dist2);
 204     }
 205 }
 206
 207 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
 208 {
 209   double dBase=sqrt(_start->distanceWithSq(*_end));
 210   double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
 211   d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
 212   return Node::areDoubleEquals(dBase,d1);
 213 }
 214
 215 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
 216 {
 217   stream << "2 1 0 1 ";
 218   fillXfigStreamForLoc(stream);
 219   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;
 220   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
 221   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
 222   stream << std::endl;
 223 }
 224
 225 void EdgeLin::update(Node *m)
 226 {
 227   updateBounds();
 228 }
 229
 230 double EdgeLin::getNormSq() const
 231 {
 232   return _start->distanceWithSq(*_end);
 233 }
 234
 235 /*!
 236  * This methods computes :
 237  * \f[
 238  * \int_{Current Edge} -ydx
 239  * \f]
 240  */
 241 double EdgeLin::getAreaOfZone() const
 242 {
 243   return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
 244 }
 245
 246 void EdgeLin::getBarycenter(double *bary) const
 247 {
 248   bary[0]=((*_start)[0]+(*_end)[0])/2.;
 249   bary[1]=((*_start)[1]+(*_end)[1])/2.;
 250 }
 251
 252 /*!
 253  * \f[
 254  * bary[0]=\int_{Current Edge} -yxdx
 255  * \f]
 256  * \f[
 257  * bary[1]=\int_{Current Edge} -\frac{y^{2}}{2}dx
 258  * \f]
 259  * To compute these 2 expressions in this class we have :
 260  * \f[
 261  * y=y_{1}+\frac{y_{2}-y_{1}}{x_{2}-x_{1}}(x-x_{1})
 262  * \f]
 263  */
 264 void EdgeLin::getBarycenterOfZone(double *bary) const
 265 {
 266   double x1=(*_start)[0];
 267   double y1=(*_start)[1];
 268   double x2=(*_end)[0];
 269   double y2=(*_end)[1];
 270   bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
 271   //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
 272   //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
 273   bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
 274 }
 275
 276 /*!
 277  * Here \a this is not used (contrary to EdgeArcCircle class).
 278  */
 279 void EdgeLin::getMiddleOfPoints(const double *p1, const double *p2, double *mid) const
 280 {
 281   mid[0]=(p1[0]+p2[0])/2.;
 282   mid[1]=(p1[1]+p2[1])/2.;
 283 }
 284
 285 double EdgeLin::getCurveLength() const
 286 {
 287   double x=(*_start)[0]-(*_end)[0];
 288   double y=(*_start)[1]-(*_end)[1];
 289   return sqrt(x*x+y*y);
 290 }
 291
 292 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
 293 {
 294   return new EdgeLin(start,end,direction);
 295 }
 296
 297 /*!
 298  * No precision should be introduced here. Just think as if precision was perfect.
 299  */
 300 void EdgeLin::updateBounds()
 301 {
 302   _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]));
 303 }
 304
 305 double EdgeLin::getCharactValueEng(const double *node) const
 306 {
 307   double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
 308   double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
 309   return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);
 310 }
 311
 312 void EdgeLin::fillGlobalInfoAbs(bool direction, const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
 313                                 std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
 314 {
 315   int tmp[2];
 316   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
 317   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
 318   if(direction)
 319     {
 320       edgesThis.push_back(tmp[0]);
 321       edgesThis.push_back(tmp[1]);
 322     }
 323   else
 324     {
 325       edgesThis.push_back(tmp[1]);
 326       edgesThis.push_back(tmp[0]);
 327     }
 328 }
 329
 330 void EdgeLin::fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
 331                                  std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
 332 {
 333   _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
 334   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
 335 }