Back to index

salome-med  6.5.0
InterpKernelGeo2DEdgeLin.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "InterpKernelGeo2DEdgeLin.hxx"
00021 #include "InterpKernelGeo2DNode.hxx"
00022 #include "InterpKernelException.hxx"
00023 #include "NormalizedUnstructuredMesh.hxx"
00024 
00025 using namespace INTERP_KERNEL;
00026 
00027 namespace INTERP_KERNEL
00028 {
00029   extern const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
00030 }
00031 
00032 SegSegIntersector::SegSegIntersector(const EdgeLin& e1, const EdgeLin& e2):SameTypeEdgeIntersector(e1,e2)
00033 {
00034   _matrix[0]=(*(e2.getStartNode()))[0]-(*(e2.getEndNode()))[0];
00035   _matrix[1]=(*(e1.getEndNode()))[0]-(*(e1.getStartNode()))[0];
00036   _matrix[2]=(*(e2.getStartNode()))[1]-(*(e2.getEndNode()))[1];
00037   _matrix[3]=(*(e1.getEndNode()))[1]-(*(e1.getStartNode()))[1];
00038   _col[0]=_matrix[3]*(*(e1.getStartNode()))[0]-_matrix[1]*(*(e1.getStartNode()))[1];
00039   _col[1]=-_matrix[2]*(*(e2.getStartNode()))[0]+_matrix[0]*(*(e2.getStartNode()))[1];
00040   //Little trick to avoid problems if 'e1' and 'e2' are colinears and along Ox or Oy axes.
00041   if(fabs(_matrix[3])>fabs(_matrix[1]))
00042     _ind=0;
00043   else
00044     _ind=1;
00045 }
00046 
00051 bool SegSegIntersector::haveTheySameDirection() const
00052 {
00053   return (_matrix[3]*_matrix[1]+_matrix[2]*_matrix[0])>0.;
00054   //return (_matrix[_ind?1:0]>0. && _matrix[_ind?3:2]>0.) || (_matrix[_ind?1:0]<0. && _matrix[_ind?3:2]<0.);
00055 }
00056 
00060 void SegSegIntersector::getPlacements(Node *start, Node *end, TypeOfLocInEdge& whereStart, TypeOfLocInEdge& whereEnd, MergePoints& commonNode) const
00061 {
00062   getCurveAbscisse(start,whereStart,commonNode);
00063   getCurveAbscisse(end,whereEnd,commonNode);
00064 }
00065 
00066 void SegSegIntersector::getCurveAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode) const
00067 {
00068   bool obvious;
00069   obviousCaseForCurvAbscisse(node,where,commonNode,obvious);
00070   if(obvious)
00071     return ;
00072   double ret=((*node)[!_ind]-(*_e1.getStartNode())[!_ind])/((*_e1.getEndNode())[!_ind]-(*_e1.getStartNode())[!_ind]);
00073   if(ret>0. && ret <1.)
00074     where=INSIDE;
00075   else if(ret<0.)
00076     where=OUT_BEFORE;
00077   else
00078     where=OUT_AFTER;
00079 }
00080 
00084 std::list< IntersectElement > SegSegIntersector::getIntersectionsCharacteristicVal() const
00085 {
00086   std::list< IntersectElement > ret;
00087   double x=_matrix[0]*_col[0]+_matrix[1]*_col[1];
00088   double y=_matrix[2]*_col[0]+_matrix[3]*_col[1];
00089   //Only one intersect point possible
00090   Node *node=new Node(x,y);
00091   node->declareOn();
00092   bool i_1S=_e1.getStartNode()->isEqual(*node);
00093   bool i_1E=_e1.getEndNode()->isEqual(*node);
00094   bool i_2S=_e2.getStartNode()->isEqual(*node);
00095   bool i_2E=_e2.getEndNode()->isEqual(*node);
00096   ret.push_back(IntersectElement(_e1.getCharactValue(*node),
00097                                  _e2.getCharactValue(*node),
00098                                  i_1S,i_1E,i_2S,i_2E,node,_e1,_e2,keepOrder()));
00099   return ret;
00100 }
00101 
00107 bool SegSegIntersector::areColinears() const
00108 {
00109   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
00110   return fabs(determinant)<QUADRATIC_PLANAR::_arc_detection_precision;
00111 }
00112 
00120 void SegSegIntersector::areOverlappedOrOnlyColinears(const Bounds *whereToFind, bool& colinearity, bool& areOverlapped)
00121 {
00122   double determinant=_matrix[0]*_matrix[3]-_matrix[1]*_matrix[2];
00123   if(fabs(determinant)>2.*QUADRATIC_PLANAR::_precision)//2*_precision due to max of offset on _start and _end
00124     {
00125       colinearity=false; areOverlapped=false;
00126       _matrix[0]/=determinant; _matrix[1]/=determinant; _matrix[2]/=determinant; _matrix[3]/=determinant;
00127     }
00128   else
00129     {
00130       colinearity=true;
00131       //retrieving initial matrix
00132       double tmp=_matrix[0]; _matrix[0]=_matrix[3]; _matrix[3]=tmp;
00133       _matrix[1]=-_matrix[1]; _matrix[2]=-_matrix[2];
00134       //
00135       double deno=sqrt(_matrix[0]*_matrix[0]+_matrix[1]*_matrix[1]);
00136       double x=(*(_e1.getStartNode()))[0]-(*(_e2.getStartNode()))[0];
00137       double y=(*(_e1.getStartNode()))[1]-(*(_e2.getStartNode()))[1];
00138       areOverlapped=fabs((_matrix[1]*y+_matrix[0]*x)/deno)<QUADRATIC_PLANAR::_precision;
00139     }
00140 }
00141 
00142 EdgeLin::EdgeLin(std::istream& lineInXfig)
00143 {
00144   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
00145   lineInXfig.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
00146   _start=new Node(lineInXfig);
00147   _end=new Node(lineInXfig);
00148   updateBounds();
00149 }
00150 
00151 EdgeLin::EdgeLin(Node *start, Node *end, bool direction):Edge(start,end,direction)
00152 {
00153   updateBounds();
00154 }
00155 
00156 EdgeLin::EdgeLin(double sX, double sY, double eX, double eY):Edge(sX,sY,eX,eY)
00157 {
00158   updateBounds();
00159 }
00160 
00161 EdgeLin::~EdgeLin()
00162 {
00163 }
00164 
00168 bool EdgeLin::isIn(double characterVal) const
00169 {
00170   return characterVal>0. && characterVal<1.;
00171 }
00172 
00173 Node *EdgeLin::buildRepresentantOfMySelf() const
00174 {
00175   return new Node(((*(_start))[0]+(*(_end))[0])/2.,((*(_start))[1]+(*(_end))[1])/2.);
00176 }
00177 
00178 double EdgeLin::getCharactValue(const Node& node) const
00179 {
00180   return getCharactValueEng(node);
00181 }
00182 
00183 double EdgeLin::getCharactValueBtw0And1(const Node& node) const
00184 {
00185   return getCharactValueEng(node);
00186 }
00187 
00188 double EdgeLin::getDistanceToPoint(const double *pt) const
00189 {
00190   double loc=getCharactValueEng(pt);
00191   if(loc>0. && loc<1.)
00192     {
00193       double tmp[2];
00194       tmp[0]=(*_start)[0]*(1-loc)+loc*(*_end)[0];
00195       tmp[1]=(*_start)[1]*(1-loc)+loc*(*_end)[1];
00196       return Node::distanceBtw2Pt(pt,tmp);
00197     }
00198   else
00199     {
00200       double dist1=Node::distanceBtw2Pt(*_start,pt);
00201       double dist2=Node::distanceBtw2Pt(*_end,pt);
00202       return std::min(dist1,dist2);
00203     }
00204 }
00205 
00206 bool EdgeLin::isNodeLyingOn(const double *coordOfNode) const
00207 {
00208   double dBase=sqrt(_start->distanceWithSq(*_end));
00209   double d1=Node::distanceBtw2Pt(*_start,coordOfNode);
00210   d1+=Node::distanceBtw2Pt(*_end,coordOfNode);
00211   return Node::areDoubleEquals(dBase,d1);
00212 }
00213 
00214 void EdgeLin::dumpInXfigFile(std::ostream& stream, bool direction, int resolution, const Bounds& box) const
00215 {
00216   stream << "2 1 0 1 ";
00217   fillXfigStreamForLoc(stream);
00218   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;
00219   direction?_start->dumpInXfigFile(stream,resolution,box):_end->dumpInXfigFile(stream,resolution,box);
00220   direction?_end->dumpInXfigFile(stream,resolution,box):_start->dumpInXfigFile(stream,resolution,box);
00221   stream << std::endl;
00222 }
00223 
00224 void EdgeLin::update(Node *m)
00225 {
00226   updateBounds();
00227 }
00228 
00229 double EdgeLin::getNormSq() const
00230 {
00231   return _start->distanceWithSq(*_end);
00232 }
00233 
00240 double EdgeLin::getAreaOfZone() const
00241 {
00242   return ((*_start)[0]-(*_end)[0])*((*_start)[1]+(*_end)[1])/2.;
00243 }
00244 
00245 void EdgeLin::getBarycenter(double *bary) const
00246 {
00247   bary[0]=((*_start)[0]+(*_end)[0])/2.;
00248   bary[1]=((*_start)[1]+(*_end)[1])/2.;
00249 }
00250 
00263 void EdgeLin::getBarycenterOfZone(double *bary) const
00264 {
00265   double x1=(*_start)[0];
00266   double y1=(*_start)[1];
00267   double x2=(*_end)[0];
00268   double y2=(*_end)[1];
00269   bary[0]=(x1-x2)*(y1*(2.*x1+x2)+y2*(2.*x2+x1))/6.;
00270   //bary[0]+=(y1-y2)*(x2*x2/3.-(x1*x2+x1*x1)/6.)+y1*(x1*x1-x2*x2)/2.;
00271   //bary[0]+=(y1-y2)*((x2*x2+x1*x2+x1*x1)/3.-(x2+x1)*x1/2.)+y1*(x1*x1-x2*x2)/2.;
00272   bary[1]=(x1-x2)*(y1*(y1+y2)+y2*y2)/6.;
00273 }
00274 
00275 double EdgeLin::getCurveLength() const
00276 {
00277   double x=(*_start)[0]-(*_end)[0];
00278   double y=(*_start)[1]-(*_end)[1];
00279   return sqrt(x*x+y*y);
00280 }
00281 
00282 Edge *EdgeLin::buildEdgeLyingOnMe(Node *start, Node *end, bool direction) const
00283 {
00284   return new EdgeLin(start,end,direction);
00285 }
00286 
00290 void EdgeLin::updateBounds()
00291 {
00292   _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]));
00293 }
00294 
00295 double EdgeLin::getCharactValueEng(const double *node) const
00296 {
00297   double car1_1x=node[0]-(*(_start))[0]; double car1_2x=(*(_end))[0]-(*(_start))[0];
00298   double car1_1y=node[1]-(*(_start))[1]; double car1_2y=(*(_end))[1]-(*(_start))[1];
00299   return (car1_1x*car1_2x+car1_1y*car1_2y)/(car1_2x*car1_2x+car1_2y*car1_2y);
00300 }
00301 
00302 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,
00303                                 std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
00304 {
00305   int tmp[2];
00306   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
00307   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
00308   if(direction)
00309     {
00310       edgesThis.push_back(tmp[0]);
00311       edgesThis.push_back(tmp[1]);
00312     }
00313   else
00314     {
00315       edgesThis.push_back(tmp[1]);
00316       edgesThis.push_back(tmp[0]);
00317     }
00318 }
00319 
00320 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,
00321                                  std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
00322 {
00323   _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
00324   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
00325 }