src/HEXABLOCK/HexCylinder.cxx

   1
   2 // C++ : Gestion des cylindres
   3
   4 // Copyright (C) 2009-2013  CEA/DEN, EDF R&D
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #include "HexCylinder.hxx"
  23 #include "HexVertex.hxx"
  24 #include "HexVector.hxx"
  25 #include "HexXmlWriter.hxx"
  26 #include "HexDocument.hxx"
  27
  28 #include <cmath>
  29
  30 BEGIN_NAMESPACE_HEXA
  31
  32 bool is_out (double val, double orig, double lg);
  33
  34 // ======================================================== Constructeur
  35 Cylinder::Cylinder(Vertex* b, Vector* v, double r, double h)
  36         : EltBase (b->dad(), EL_CYLINDER)
  37 {
  38     c_base   = b;
  39     c_dir    = v;
  40     c_radius = r;
  41     c_height = h;
  42
  43     if (el_root->debug ())
  44        {
  45        Echo ("---------------------------------- AddCylinder");
  46        HexDisplay (c_radius);
  47        HexDisplay (c_height);
  48        HexDump (c_base);
  49        HexDump (c_dir);
  50        cout << endl;
  51        }
  52
  53     if (c_base==NULL || c_dir==NULL || el_root==NULL)
  54        setError ();
  55     else
  56        {
  57        double tol   = el_root->getTolerance ();
  58        double norme = c_dir  ->getNorm ();
  59        if (c_radius <= tol || c_height <= tol || norme <= tol)
  60            setError ();
  61        }
  62 }
  63 // ======================================================== rdiffers
  64 bool rdiffers (double v1, double v2)
  65 {
  66    double dd = v1-v2;
  67    return dd*dd > 1e-12;
  68 }
  69 // ======================================================== norme
  70 double norme (double px, double py, double pz)
  71 {
  72    double res = sqrt (px*px + py*py + pz*pz);
  73    return res;
  74 }
  75 // ======================================================== interCylinder
  76 ///  Intersection de 2 droites (c1,v1) et (c2,v2)
  77 ///  Le point M(x,yz) est solutions de 2 equations :
  78 ///            (c1,m) = k * v1
  79 ///            (c2,m) vectoriel v2 = 0
  80 ///
  81 ///  x = xc1 + k * xv1                  (1)
  82 ///  y = yc1 + k * yv1                  (2)
  83 ///  z = zc1 + k * zv1                  (3)
  84 ///  (x-xc2) * yv2 - (y-yc2) * xv2 = 0  (4)
  85 ///  (z-zc2) * yv2 - (y-yc2) * zv2 = 0  (5)
  86 ///  (x-xc2) * zv2 - (z-zc2) * xv2 = 0  (6)
  87 ///
  88 ///  En substituant les (x,y,z)
  89 ///
  90 ///  (k*xv1+xc1-xc2) * yv2 - (k*yv1+yc1-yc2) * xv2 = 0  (4)
  91 ///  (k*zv1+zc1-zc2) * yv2 - (k*yv1+yc1-yc2) * zv2 = 0  (5)
  92 ///  (k*xv1+xc1-xc2) * zv2 - (k*zv1+zc1-zc2) * xv2 = 0  (6)
  93 ///
  94 ///  k * (xv1*yv2 - yv1*xv2) + (xc1-xc2) * yv2 - (yc1-yc2) * xv2 = 0 (4)
  95 ///  k * dxy - nxy = 0      (4)
  96 ///  k * dzy - nzy = 0      (5)
  97 ///  k * dxz - nxz = 0      (6)
  98 ///
  99 ///  nxy = (xc2-xc1) * yv2 - (yc2-yc1) * xv2
 100 ///  dxy =  xv1*yv2 - yv1*xv2
 101 ///
 102 // ======================================================== interCylinder
 103 Vertex* Cylinder::interCylinder (Cylinder* small, bool& left, bool& right)
 104 {
 105    Real3 orig;
 106    int ier = interCylinder (small, left, right, orig);
 107    if (ier!=HOK)
 108        return NULL;
 109
 110     Vertex* sol = new Vertex (c_base->dad(), orig[0],  orig[1],  orig[2]);
 111     return sol;
 112 }
 113 ///
 114 // ======================================================== interCylinder
 115 int Cylinder::interCylinder (Cylinder* small, bool& left, bool& right,
 116                              double* orig)
 117 {
 118     left = right = true;
 119     if (el_root->debug ())
 120        {
 121        Echo ("---------------------------------- interCylinders");
 122        PutData (this->c_radius);
 123        PutData (small->c_radius);
 124        PutData (this->c_height);
 125        PutData (small->c_height);
 126        HexDump (this->c_base);
 127        HexDump (small->c_base);
 128        HexDump (this->c_dir);
 129        HexDump (small->c_dir);
 130        }
 131
 132     double norm1 = c_dir->getNorm();
 133     double norm2 = small->c_dir->getNorm();
 134
 135     double xc1 = c_base->getX();
 136     double yc1 = c_base->getY();
 137     double zc1 = c_base->getZ();
 138
 139     double xc2 = small->c_base->getX();
 140     double yc2 = small->c_base->getY();
 141     double zc2 = small->c_base->getZ();
 142
 143     double xv1 = c_dir->getDx()/norm1;
 144     double yv1 = c_dir->getDy()/norm1;
 145     double zv1 = c_dir->getDz()/norm1;
 146
 147     double xv2 = small->c_dir->getDx()/norm2;
 148     double yv2 = small->c_dir->getDy()/norm2;
 149     double zv2 = small->c_dir->getDz()/norm2;
 150
 151     double nxyz [DIM3] = { (xc2-xc1) * yv2 - (yc2-yc1) * xv2,
 152                            (zc2-zc1) * yv2 - (yc2-yc1) * zv2,
 153                            (xc2-xc1) * zv2 - (zc2-zc1) * xv2 };
 154
 155     double dxyz [DIM3] = { xv1*yv2 - yv1*xv2,
 156                            zv1*yv2 - yv1*zv2,
 157                            xv1*zv2 - zv1*xv2 };
 158     double lambda = 0;
 159     bool   prems  = true;
 160     for (int dd=0 ; dd<DIM3 ; dd++)
 161         {
 162                                       // dxyz!=0 : calcul de lambda
 163         if (rdiffers (dxyz [dd], ZEROR))
 164            {
 165            double kk = nxyz [dd] / dxyz[dd];
 166                                       // 1er lambda trouve
 167            if (prems)
 168               {
 169               prems  = false;
 170               lambda = kk;
 171               }
 172                                       // Solutions incompatibles
 173            else if (rdiffers (kk, lambda))
 174               {
 175               cout << "*** InterCylinders : Solutions incompatibles " << endl;
 176               return HERR;
 177               }
 178                                       // equation O*lamda = 0 : ignoree
 179         // else ....
 180            }
 181                                       // dxyz=0 et nxyz!=0 : Pas de solution
 182         else if (rdiffers (nxyz [dd], ZEROR))
 183            {
 184            cout << "*** InterCylinders : intersection vide " << endl;
 185            return HERR;
 186            }
 187         }
 188
 189                                       // Aucune solution trouvee
 190     if (prems)
 191        {
 192        HexDisplay (prems);
 193        cout << "*** InterCylinders : intersection non trouvee " << endl;
 194        return HERR;
 195        }
 196                                     // Intersection des droites
 197     orig [dir_x] = xc1 + lambda*xv1;
 198     orig [dir_y] = yc1 + lambda*yv1;
 199     orig [dir_z] = zc1 + lambda*zv1;
 200
 201     Real3 base1, dir1, extr1;
 202     Real3 base2, dir2, extr2;
 203
 204     this ->c_base->getPoint (base1);
 205     small->c_base->getPoint (base2);
 206
 207     this ->c_dir->getCoord (dir1);
 208     small->c_dir->getCoord (dir2);
 209
 210     for (int dd=dir_x ; dd<=dir_z ; dd++)
 211         {
 212         extr1 [dd] = base1[dd] + dir1[dd]*this ->c_height/norm1;
 213         extr2 [dd] = base2[dd] + dir2[dd]*small->c_height/norm2;
 214         }
 215
 216     double dbase1 = calc_distance (orig, base1);
 217     double dbase2 = calc_distance (orig, base2);
 218     double dextr1 = calc_distance (orig, extr1);
 219     double dextr2 = calc_distance (orig, extr2);
 220
 221     double dmin2  = std::min (dbase2, dextr2);
 222     double dmax2  = std::max (dbase2, dextr2);
 223     double dmax1  = std::max (dbase1, dextr1);
 224     double coeff  = 1.1;
 225                                     // Appartenance axe petit cylindre
 226     if (dmax2 > small->c_height + c_radius)
 227        {
 228        cout << "*** InterCylinders : Petit cylindre trop court" << endl;
 229        cout << "*** Intersection = (" << orig[0] << ", "  << orig[1]
 230                              << ", "  << orig[2] << " )" << endl;
 231        cout << "*** distance = "  << dmin2
 232             << " > gros rayon + longueur "
 233             << small->c_height + c_radius << endl;
 234        cout << endl;
 235        PutCoord (base2);
 236        PutCoord (extr2);
 237        PutData(dbase2);
 238        PutData(dextr2);
 239        PutData(dextr2);
 240        return HERR;
 241        }
 242
 243     left  = dbase2 > c_radius*coeff;
 244     right = dextr2 > c_radius*coeff;
 245                                     // Le gros cylindre doit depasser le petit
 246     if (dmax1 > c_height-coeff*small->c_radius)
 247        {
 248        cout << "*** InterCylinders : Gros cylindre trop court" << endl;
 249        cout << "*** Intersection = (" << orig[0] << ", "  << orig[1]
 250                              << ", "  << orig[2] << " )" << endl;
 251        cout << "*** distance maximale = " << dmax1 << endl;
 252        cout << endl;
 253        return HERR;
 254        }
 255
 256     if (el_root->debug ())
 257        {
 258        PutCoord   (orig);
 259        HexDisplay (left);
 260        HexDisplay (right);
 261        }
 262     return  HOK;
 263 }
 264 // ======================================================== is_out
 265 bool is_out (double val, double v1, double v2)
 266 {
 267    if (v1 < v2)
 268       return (val < v1-Epsil || val > v2+Epsil) ;
 269
 270    else if (v1 > v2)
 271       return (val < v2-Epsil || val > v1+Epsil) ;
 272
 273    else
 274       return fabs (val - v1)>Epsil;
 275 }
 276 // ========================================================= saveXml
 277 void Cylinder::saveXml  (XmlWriter* xml)
 278 {
 279    char buffer[12];
 280
 281    xml->openMark     ("Cylinder");
 282    xml->addAttribute ("c_base", c_base->getName (buffer));
 283    xml->addAttribute ("c_dir",  c_dir->getName  (buffer));
 284    xml->addAttribute ("c_radius", c_radius);
 285    xml->addAttribute ("c_height", c_height);
 286    xml->closeMark ();
 287 }
 288 END_NAMESPACE_HEXA
 289