src/HEXABLOCK/HexMatrix.hxx

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