src/HEXABLOCK/HexCylinder.cxx

   1
   2 // C++ : Gestion des cylindres
   3
   4 // Copyright (C) 2009-2012  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
  27 #include <cmath>
  28
  29 BEGIN_NAMESPACE_HEXA
  30
  31 bool is_out (double val, double orig, double lg);
  32
  33 // ======================================================== Constructeur
  34 Cylinder::Cylinder(Vertex* b, Vector* v, double r, double h)
  35         : EltBase (b->dad())
  36 {
  37     c_base   = b;
  38     c_dir    = v;
  39     c_radius = r;
  40     c_height = h;
  41
  42     if (el_root->debug ())
  43        {
  44        Echo ("---------------------------------- AddCylinder");
  45        HexDisplay (c_radius);
  46        HexDisplay (c_height);
  47        HexDump (c_base);
  48        HexDump (c_dir);
  49        cout << endl;
  50        }
  51 }
  52 // ======================================================== rdiffers
  53 bool rdiffers (double v1, double v2)
  54 {
  55    double dd = v1-v2;
  56    return dd*dd > 1e-12;
  57 }
  58 // ======================================================== norme
  59 double norme (double px, double py, double pz)
  60 {
  61    double res = sqrt (px*px + py*py + pz*pz);
  62    return res;
  63 }
  64 // ======================================================== interCylinder
  65 ///  Intersection de 2 droites (c1,v1) et (c2,v2)
  66 ///  Le point M(x,yz) est solutions de 2 equations :
  67 ///            (c1,m) = k * v1
  68 ///            (c2,m) vectoriel v2 = 0
  69 ///
  70 ///  x = xc1 + k * xv1                  (1)
  71 ///  y = yc1 + k * yv1                  (2)
  72 ///  z = zc1 + k * zv1                  (3)
  73 ///  (x-xc2) * yv2 - (y-yc2) * xv2 = 0  (4)
  74 ///  (z-zc2) * yv2 - (y-yc2) * zv2 = 0  (5)
  75 ///  (x-xc2) * zv2 - (z-zc2) * xv2 = 0  (6)
  76 ///
  77 ///  En substituant les (x,y,z)
  78 ///
  79 ///  (k*xv1+xc1-xc2) * yv2 - (k*yv1+yc1-yc2) * xv2 = 0  (4)
  80 ///  (k*zv1+zc1-zc2) * yv2 - (k*yv1+yc1-yc2) * zv2 = 0  (5)
  81 ///  (k*xv1+xc1-xc2) * zv2 - (k*zv1+zc1-zc2) * xv2 = 0  (6)
  82 ///
  83 ///  k * (xv1*yv2 - yv1*xv2) + (xc1-xc2) * yv2 - (yc1-yc2) * xv2 = 0 (4)
  84 ///  k * dxy - nxy = 0      (4)
  85 ///  k * dzy - nzy = 0      (5)
  86 ///  k * dxz - nxz = 0      (6)
  87 ///
  88 ///  nxy = (xc2-xc1) * yv2 - (yc2-yc1) * xv2
  89 ///  dxy =  xv1*yv2 - yv1*xv2
  90 ///
  91 // ======================================================== interCylinder
  92 Vertex* Cylinder::interCylinder (Cylinder* small, bool& left, bool& right)
  93 {
  94    Real3 orig;
  95    int ier = interCylinder (small, left, right, orig);
  96    if (ier!=HOK)
  97        return NULL;
  98
  99     Vertex* sol = new Vertex (c_base->dad(), orig[0],  orig[1],  orig[2]);
 100     return sol;
 101 }
 102 ///
 103 // ======================================================== interCylinder
 104 int Cylinder::interCylinder (Cylinder* small, bool& left, bool& right,
 105                              double* orig)
 106 {
 107     left = right = true;
 108     if (el_root->debug ())
 109        {
 110        Echo ("---------------------------------- interCylinders");
 111        PutData (this->c_radius);
 112        PutData (small->c_radius);
 113        PutData (this->c_height);
 114        PutData (small->c_height);
 115        HexDump (this->c_base);
 116        HexDump (small->c_base);
 117        HexDump (this->c_dir);
 118        HexDump (small->c_dir);
 119        }
 120
 121     double norm1 = c_dir->getNorm();
 122     double norm2 = small->c_dir->getNorm();
 123
 124     double xc1 = c_base->getX();
 125     double yc1 = c_base->getY();
 126     double zc1 = c_base->getZ();
 127
 128     double xc2 = small->c_base->getX();
 129     double yc2 = small->c_base->getY();
 130     double zc2 = small->c_base->getZ();
 131
 132     double xv1 = c_dir->getDx()/norm1;
 133     double yv1 = c_dir->getDy()/norm1;
 134     double zv1 = c_dir->getDz()/norm1;
 135
 136     double xv2 = small->c_dir->getDx()/norm2;
 137     double yv2 = small->c_dir->getDy()/norm2;
 138     double zv2 = small->c_dir->getDz()/norm2;
 139
 140     double nxyz [DIM3] = { (xc2-xc1) * yv2 - (yc2-yc1) * xv2,
 141                            (zc2-zc1) * yv2 - (yc2-yc1) * zv2,
 142                            (xc2-xc1) * zv2 - (zc2-zc1) * xv2 };
 143
 144     double dxyz [DIM3] = { xv1*yv2 - yv1*xv2,
 145                            zv1*yv2 - yv1*zv2,
 146                            xv1*zv2 - zv1*xv2 };
 147     double lambda = 0;
 148     bool   prems  = true;
 149     for (int dd=0 ; dd<DIM3 ; dd++)
 150         {
 151                                       // dxyz!=0 : calcul de lambda
 152         if (rdiffers (dxyz [dd], ZEROR))
 153            {
 154            double kk = nxyz [dd] / dxyz[dd];
 155                                       // 1er lambda trouve
 156            if (prems)
 157               {
 158               prems  = false;
 159               lambda = kk;
 160               }
 161                                       // Solutions incompatibles
 162            else if (rdiffers (kk, lambda))
 163               {
 164               cout << "*** InterCylinders : Solutions incompatibles " << endl;
 165               return HERR;
 166               }
 167                                       // equation O*lamda = 0 : ignoree
 168         // else ....
 169            }
 170                                       // dxyz=0 et nxyz!=0 : Pas de solution
 171         else if (rdiffers (nxyz [dd], ZEROR))
 172            {
 173            cout << "*** InterCylinders : intersection vide " << endl;
 174            return HERR;
 175            }
 176         }
 177
 178                                       // Aucune solution trouvee
 179     if (prems)
 180        {
 181        HexDisplay (prems);
 182        cout << "*** InterCylinders : intersection non trouvee " << endl;
 183        return HERR;
 184        }
 185                                     // Intersection des droites
 186     orig [dir_x] = xc1 + lambda*xv1;
 187     orig [dir_y] = yc1 + lambda*yv1;
 188     orig [dir_z] = zc1 + lambda*zv1;
 189
 190     Real3 base1, dir1, extr1;
 191     Real3 base2, dir2, extr2;
 192
 193     this ->c_base->getPoint (base1);
 194     small->c_base->getPoint (base2);
 195
 196     this ->c_dir->getCoord (dir1);
 197     small->c_dir->getCoord (dir2);
 198
 199     for (int dd=dir_x ; dd<=dir_z ; dd++)
 200         {
 201         extr1 [dd] = base1[dd] + dir1[dd]*this ->c_height/norm1;
 202         extr2 [dd] = base2[dd] + dir2[dd]*small->c_height/norm2;
 203         }
 204
 205     double dbase1 = calc_distance (orig, base1);
 206     double dbase2 = calc_distance (orig, base2);
 207     double dextr1 = calc_distance (orig, extr1);
 208     double dextr2 = calc_distance (orig, extr2);
 209
 210     double dmin2  = std::min (dbase2, dextr2);
 211     double dmax2  = std::max (dbase2, dextr2);
 212     double dmax1  = std::max (dbase1, dextr1);
 213     double coeff  = 1.1;
 214                                     // Appartenance axe petit cylindre
 215     if (dmax2 > small->c_height + c_radius)
 216        {
 217        cout << "*** InterCylinders : Petit cylindre trop court" << endl;
 218        cout << "*** Intersection = (" << orig[0] << ", "  << orig[1]
 219                              << ", "  << orig[2] << " )" << endl;
 220        cout << "*** distance = "  << dmin2
 221             << " > gros rayon + longueur "
 222             << small->c_height + c_radius << endl;
 223        cout << endl;
 224        PutCoord (base2);
 225        PutCoord (extr2);
 226        PutData(dbase2);
 227        PutData(dextr2);
 228        PutData(dextr2);
 229        return HERR;
 230        }
 231
 232     left  = dbase2 > c_radius*coeff;
 233     right = dextr2 > c_radius*coeff;
 234                                     // Le gros cylindre doit depasser le petit
 235     if (dmax1 > c_height-coeff*small->c_radius)
 236        {
 237        cout << "*** InterCylinders : Gros cylindre trop court" << endl;
 238        cout << "*** Intersection = (" << orig[0] << ", "  << orig[1]
 239                              << ", "  << orig[2] << " )" << endl;
 240        cout << "*** distance maximale = " << dmax1 << endl;
 241        cout << endl;
 242        return HERR;
 243        }
 244
 245     if (el_root->debug ())
 246        {
 247        PutCoord   (orig);
 248        HexDisplay (left);
 249        HexDisplay (right);
 250        }
 251     return  HOK;
 252 }
 253 // ======================================================== is_out
 254 bool is_out (double val, double v1, double v2)
 255 {
 256    if (v1 < v2)
 257       return (val < v1-Epsil || val > v2+Epsil) ;
 258
 259    else if (v1 > v2)
 260       return (val < v2-Epsil || val > v1+Epsil) ;
 261
 262    else
 263       return fabs (val - v1)>Epsil;
 264 }
 265 // ========================================================= saveXml
 266 void Cylinder::saveXml  (XmlWriter* xml)
 267 {
 268    char buffer[12];
 269
 270    xml->openMark     ("Cylinder");
 271    xml->addAttribute ("c_base", c_base->getName (buffer));
 272    xml->addAttribute ("c_dir",  c_dir->getName  (buffer));
 273    xml->addAttribute ("c_radius", c_radius);
 274    xml->addAttribute ("c_height", c_height);
 275    xml->closeMark ();
 276 }
 277 END_NAMESPACE_HEXA
 278