src/HEXABLOCK/HexMatrix.hxx

   1
   2 // class : Les matrices
   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
  23 #ifndef __MATRIX_H
  24 #define __MATRIX_H
  25
  26 #include "Hex_defines.hxx"
  27 #include "HexVertex.hxx"
  28 #include "HexVector.hxx"
  29
  30 #include <cmath>
  31
  32 BEGIN_NAMESPACE_HEXA
  33
  34 class HEXABLOCKENGINE_EXPORT Matrix
  35 {
  36 public:
  37     Matrix ();
  38     int defTranslation   (Vector* depl);
  39     int defTranslation   (double* depl);
  40     int defScale         (Vertex* center, double scale);
  41     int defScale         (double* center, double scale);
  42     int defRotation      (Vertex* center, Vector* axe, double degres);
  43     int defRotation      (Vertex* center, double* axe, double degres);
  44     int defSymmetryPoint (Vertex* center);
  45     int defSymmetryLine  (Vertex* center, Vector* dir);
  46     int defSymmetryPlane (Vertex* center, Vector* normale);
  47
  48     void define (double* orig, double* iprim, double* jprim, double* kprim);
  49
  50     void perform (Vertex* noeud);
  51     void perform (double* point);
  52     void perform (double* point, double* result);
  53
  54     void getCoeff (double& a11, double& a12, double& a13, double& a14,
  55                    double& a21, double& a22, double& a23, double& a24,
  56                    double& a31, double& a32, double& a33, double& a34);
  57 private:
  58     void erase();
  59 private:
  60     double mat11, mat12, mat13, mat14;
  61     double mat21, mat22, mat23, mat24;
  62     double mat31, mat32, mat33, mat34;
  63 };
  64 // ------------------------------------------- Inlining
  65 // ========================================================= Constructeur
  66 inline Matrix::Matrix ()
  67 {
  68    erase ();
  69 }
  70 // ========================================================= erase
  71 inline void Matrix::erase ()
  72 {
  73    mat12 = mat13 = mat14 = 0.0;
  74    mat21 = mat23 = mat24 = 0.0;
  75    mat31 = mat32 = mat34 = 0.0;
  76    mat11 = mat22 = mat33 = 1.0;
  77 }
  78 // ========================================================= perform
  79 inline void Matrix::perform (Vertex* noeud)
  80 {
  81    double px, py, pz;
  82    px = mat11*noeud->getX()+mat12*noeud->getY()+mat13*noeud->getZ()+mat14;
  83    py = mat21*noeud->getX()+mat22*noeud->getY()+mat23*noeud->getZ()+mat24;
  84    pz = mat31*noeud->getX()+mat32*noeud->getY()+mat33*noeud->getZ()+mat34;
  85
  86    noeud->setCoord (px, py, pz);
  87 }
  88 // ========================================================= defTranslation
  89 inline int Matrix::defTranslation (Vector* boulevard)
  90 {
  91    erase();
  92    mat11 = mat22 = mat33 = 1.0;
  93
  94    mat14 = boulevard->getDx ();
  95    mat24 = boulevard->getDy ();
  96    mat34 = boulevard->getDz ();
  97
  98    return HOK;
  99 }
 100 // ========================================================= defTranslation (2)
 101 inline int Matrix::defTranslation (double* decal)
 102 {
 103    erase();
 104    mat11 = mat22 = mat33 = 1.0;
 105
 106    mat14 = decal [dir_x];
 107    mat24 = decal [dir_y];
 108    mat34 = decal [dir_z];
 109    return HOK;
 110 }
 111 // ========================================================= defScale
 112 inline int Matrix::defScale (double* center, double scale)
 113 {
 114    erase();
 115    mat11 = mat22 = mat33 = scale;
 116
 117    mat14 = (1-scale) * center[dir_x];
 118    mat24 = (1-scale) * center[dir_y];
 119    mat34 = (1-scale) * center[dir_z];
 120
 121    return HOK;
 122 }
 123 // ========================================================= defScale
 124 inline int Matrix::defScale (Vertex* center, double scale)
 125 {
 126    if (center==NULL)
 127       {
 128       erase();
 129       return HERR;
 130       }
 131
 132    Real3 coord;
 133    int ier = defScale (center->getPoint (coord ), scale);
 134    return ier;
 135 }
 136 // ========================================================= defRotation
 137 inline int Matrix::defRotation (Vertex* center, double* axe, double degres)
 138 {
 139    erase();
 140    if (BadElement (center))
 141       return HERR;
 142
 143    double norme = calc_norme (axe);
 144    if (norme< 1e-30)
 145       return HERR;
 146
 147    double ux = axe [dir_x] / norme;
 148    double uy = axe [dir_y] / norme;
 149    double uz = axe [dir_z] / norme;
 150
 151    double cx = center->getX ();
 152    double cy = center->getY ();
 153    double cz = center->getZ ();
 154
 155    double cost = cos (degres*M_PI/180);
 156    double sint = sin (degres*M_PI/180);
 157
 158    mat11 = ux*ux * (1-cost) + cost;
 159    mat12 = ux*uy * (1-cost) - uz*sint;
 160    mat13 = ux*uz * (1-cost) + uy*sint;
 161
 162    mat21 = ux*uy * (1-cost) + uz*sint;
 163    mat22 = uy*uy * (1-cost) + cost;
 164    mat23 = uy*uz * (1-cost) - ux*sint;
 165
 166    mat31 = ux*uz * (1-cost) - uy*sint;
 167    mat32 = uy*uz * (1-cost) + ux*sint;
 168    mat33 = uz*uz * (1-cost) + cost;
 169
 170    mat14 = cx - mat11*cx - mat12*cy - mat13*cz;
 171    mat24 = cy - mat21*cx - mat22*cy - mat23*cz;
 172    mat34 = cz - mat31*cx - mat32*cy - mat33*cz;
 173
 174    return HOK;
 175 }
 176 // ========================================================= defRotation
 177 inline int Matrix::defRotation (Vertex* center, Vector* dir, double degres)
 178 {
 179    if (BadElement (dir))
 180       return HERR;
 181
 182    Real3 axe;
 183    dir->getCoord (axe);
 184    int ier = defRotation (center, axe, degres);
 185    return ier;
 186 }
 187 // ========================================================= defSymmetryPoint
 188 inline int Matrix::defSymmetryPoint (Vertex* center)
 189 {
 190    erase();
 191
 192    mat11 = mat22 = mat33 = -1;
 193
 194    mat14 = 2 * center->getX();
 195    mat24 = 2 * center->getY();
 196    mat34 = 2 * center->getZ();
 197
 198    return HOK;
 199 }
 200 // ========================================================= defSymmetryLine
 201 //     MH.d = 0        (1)
 202 //     CH  = lambda*d  (2)
 203 //     MM' = 2MH       (3)
 204 //
 205 // (1) et (2) => lambda = ((x-xc)*xd + (y-yc)*yd + (z-zc)*zd) / norme(d)
 206 //
 207 //     MM' = 2MH (3)
 208 // <=> MO + OM' =  2 (MO + OC + CH)
 209 // <=>      OM' =  MO + 2.OC + 2.CH
 210 // <=>      OM' = -OM + 2.OC + 2.lambda.d   (2) et (3)
 211 //
 212 //           x' = -x  + 2*xc + 2*xd*((x-xc)*xd + (y-yc)*yd + (z-zc)*zd)
 213 //           y' = -y  + 2*yc + 2*yd*((x-xc)*xd + (y-yc)*yd + (z-zc)*zd)
 214 //           z' = -z  + 2*zc + 2*zd*((x-xc)*xd + (y-yc)*yd + (z-zc)*zd)
 215 //
 216 inline int Matrix::defSymmetryLine (Vertex* center, Vector* dir)
 217 {
 218    erase ();
 219
 220    double normed = dir->getNorm ();
 221    if (normed< 1e-30)
 222       return HERR;
 223
 224    double xc =  center->getX();
 225    double yc =  center->getY();
 226    double zc =  center->getZ();
 227
 228    double xd = dir->getDx() / normed;
 229    double yd = dir->getDy() / normed;
 230    double zd = dir->getDz() / normed;
 231
 232    mat11 = 2*xd*xd -1;
 233    mat12 = 2*xd*yd;
 234    mat13 = 2*xd*zd;
 235    mat14 = 2*(xc - xd*(xc*xd + yc*yd + zc*zd));
 236
 237    mat21 = 2*yd*xd;
 238    mat22 = 2*yd*yd - 1;
 239    mat23 = 2*yd*zd;
 240    mat24 = 2*(yc - yd*(xc*xd + yc*yd + zc*zd));
 241
 242    mat31 = 2*zd*xd;
 243    mat32 = 2*zd*yd;
 244    mat33 = 2*zd*zd - 1;
 245    mat34 = 2*(zc - zd*(xc*xd + yc*yd + zc*zd));
 246
 247    return HOK;
 248 }
 249 // ========================================================= defSymmetryPlane
 250 //     CH.n = 0         (1)
 251 //     MH   = lambda*n  (2)
 252 //     MM'  = 2MH       (3)
 253 //
 254 // (1) et (2) => lambda = ((x-xc)*xn + (y-yc)*yn + (z-zc)*zn) / norme(n)
 255 //
 256 //     MM' = 2MH (3)
 257 // <=> MO + OM' =  2.lambda.n
 258 // <=>      OM' =  OM + 2.lambda.n
 259 //
 260 //           x' = x + 2*lambda*xn
 261 //           y' = y + 2*lambda*yn
 262 //           z' = z + 2*lambda*zn
 263 //
 264 //           x' = x + 2*xn*((x-xc)*xn + (y-yc)*yn + (z-zc)*zn)
 265 //           y' = y + 2*yn*((x-xc)*xn + (y-yc)*yn + (z-zc)*zn)
 266 //           z' = z + 2*zn*((x-xc)*xn + (y-yc)*yn + (z-zc)*zn)
 267 //
 268 inline int Matrix::defSymmetryPlane (Vertex* center, Vector* normale)
 269 {
 270    erase ();
 271
 272    double normed = normale->getNorm ();
 273    if (normed< 1e-30)
 274       return HERR;
 275
 276    double xc =  center->getX();
 277    double yc =  center->getY();
 278    double zc =  center->getZ();
 279
 280    double xn = normale->getDx() / normed;
 281    double yn = normale->getDy() / normed;
 282    double zn = normale->getDz() / normed;
 283
 284    mat11 = -2*xn*xn + 1;
 285    mat12 = -2*xn*yn;
 286    mat13 = -2*xn*zn;
 287    mat14 =  2*xn*(xc*xn + yc*yn + zc*zn);
 288
 289    mat21 = -2*yn*xn;
 290    mat22 = -2*yn*yn + 1;
 291    mat23 = -2*yn*zn;
 292    mat24 =  2*yn*(xc*xn + yc*yn + zc*zn);
 293
 294    mat31 = -2*zn*xn;
 295    mat32 = -2*zn*yn;
 296    mat33 = -2*zn*zn + 1;
 297    mat34 =  2*zn*(xc*xn + yc*yn + zc*zn);
 298
 299    return HOK;
 300 }
 301 // ========================================================= define
 302 inline void Matrix::define (double* omega, double* iprim, double* jprim,
 303                                                           double* kprim)
 304 {
 305    mat11 = iprim [dir_x];
 306    mat21 = iprim [dir_y];
 307    mat31 = iprim [dir_z];
 308
 309    mat12 = jprim [dir_x];
 310    mat22 = jprim [dir_y];
 311    mat32 = jprim [dir_z];
 312
 313    mat13 = kprim [dir_x];
 314    mat23 = kprim [dir_y];
 315    mat33 = kprim [dir_z];
 316
 317    mat14 = omega [dir_x];
 318    mat24 = omega [dir_y];
 319    mat34 = omega [dir_z];
 320 }
 321 // ========================================================= getCoeff
 322 inline void Matrix::getCoeff(double& a11, double& a12, double& a13, double& a14,
 323                              double& a21, double& a22, double& a23, double& a24,
 324                              double& a31, double& a32, double& a33, double& a34)
 325 {
 326    a11 = mat11;
 327    a12 = mat12;
 328    a13 = mat13;
 329    a14 = mat14;
 330
 331    a21 = mat21;
 332    a22 = mat22;
 333    a23 = mat23;
 334    a24 = mat24;
 335
 336    a31 = mat31;
 337    a32 = mat32;
 338    a33 = mat33;
 339    a34 = mat34;
 340 }
 341 // ========================================================= perform
 342 inline void Matrix::perform (double* point)
 343 {
 344    double px = point [dir_x];
 345    double py = point [dir_y];
 346    double pz = point [dir_z];
 347
 348    point [dir_x] = mat11*px + mat12*py + mat13*pz + mat14;
 349    point [dir_y] = mat21*px + mat22*py + mat23*pz + mat24;
 350    point [dir_z] = mat31*px + mat32*py + mat33*pz + mat34;
 351 }
 352 // ========================================================= perform
 353 inline void Matrix::perform (double* point, double* result)
 354 {
 355    double px = point [dir_x];
 356    double py = point [dir_y];
 357    double pz = point [dir_z];
 358
 359    result [dir_x] = mat11*px + mat12*py + mat13*pz + mat14;
 360    result [dir_y] = mat21*px + mat22*py + mat23*pz + mat24;
 361    result [dir_z] = mat31*px + mat32*py + mat33*pz + mat34;
 362 }
 363 END_NAMESPACE_HEXA
 364 #endif