src/foam/primitives/transform/transform.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 InNamespace
  25     Foam
  26
  27 Description
  28     3D tensor transformation operations.
  29
  30 \*---------------------------------------------------------------------------*/
  31
  32 #ifndef transform_H
  33 #define transform_H
  34
  35 #include "tensorField.H"
  36 #include "symmTensor4thOrder.H"
  37 #include "diagTensor.H"
  38 #include "mathematicalConstants.H"
  39
  40 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  41
  42 namespace Foam
  43 {
  44
  45 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  46
  47 inline tensor rotationTensor
  48 (
  49     const vector& n1,
  50     const vector& n2
  51 )
  52 {
  53     return
  54         (n1 & n2)*I
  55       + (1 - (n1 & n2))*sqr(n1 ^ n2)/(magSqr(n1 ^ n2) + VSMALL)
  56       + (n2*n1 - n1*n2);
  57 }
  58
  59
  60 inline tmp<tensorField> rotationTensor
  61 (
  62     const vectorField& n1,
  63     const vectorField& n2
  64 )
  65 {
  66     return
  67         (n1 & n2)*I
  68       + (1 - (n1 & n2))*sqr(n1 ^ n2)/(magSqr(n1 ^ n2) + VSMALL)
  69       + (n2*n1 - n1*n2);
  70 }
  71
  72
  73 inline label transform(const tensor&, const label i)
  74 {
  75     return i;
  76 }
  77
  78
  79 inline scalar transform(const tensor&, const scalar s)
  80 {
  81     return s;
  82 }
  83
  84
  85 template<class Cmpt>
  86 inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
  87 {
  88     return tt & v;
  89 }
  90
  91
  92 template<class Cmpt>
  93 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
  94 {
  95     //return tt & t & tt.T();
  96     return Tensor<Cmpt>
  97     (
  98         (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.xx()
  99       + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.xy()
 100       + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.xz(),
 101
 102         (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.yx()
 103       + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.yy()
 104       + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.yz(),
 105
 106         (tt.xx()*t.xx() + tt.xy()*t.yx() + tt.xz()*t.zx())*tt.zx()
 107       + (tt.xx()*t.xy() + tt.xy()*t.yy() + tt.xz()*t.zy())*tt.zy()
 108       + (tt.xx()*t.xz() + tt.xy()*t.yz() + tt.xz()*t.zz())*tt.zz(),
 109
 110         (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.xx()
 111       + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.xy()
 112       + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.xz(),
 113
 114         (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.yx()
 115       + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.yy()
 116       + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.yz(),
 117
 118         (tt.yx()*t.xx() + tt.yy()*t.yx() + tt.yz()*t.zx())*tt.zx()
 119       + (tt.yx()*t.xy() + tt.yy()*t.yy() + tt.yz()*t.zy())*tt.zy()
 120       + (tt.yx()*t.xz() + tt.yy()*t.yz() + tt.yz()*t.zz())*tt.zz(),
 121
 122         (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.xx()
 123       + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.xy()
 124       + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.xz(),
 125
 126         (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.yx()
 127       + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.yy()
 128       + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.yz(),
 129
 130         (tt.zx()*t.xx() + tt.zy()*t.yx() + tt.zz()*t.zx())*tt.zx()
 131       + (tt.zx()*t.xy() + tt.zy()*t.yy() + tt.zz()*t.zy())*tt.zy()
 132       + (tt.zx()*t.xz() + tt.zy()*t.yz() + tt.zz()*t.zz())*tt.zz()
 133     );
 134 }
 135
 136
 137 template<class Cmpt>
 138 inline SphericalTensor<Cmpt> transform
 139 (
 140     const tensor& tt,
 141     const SphericalTensor<Cmpt>& st
 142 )
 143 {
 144     return st;
 145 }
 146
 147
 148 template<class Cmpt>
 149 inline SymmTensor4thOrder<Cmpt> transform(const tensor& L, const SymmTensor4thOrder<Cmpt>& C)
 150 {
 151   //- represent fourth order tensors in 6x6 notation then rotation is given by
 152   //- C_rotated_af = A_ba * C_cd * A_ef
 153   //- where A is a function of L
 154
 155   const double s = ::sqrt(2);
 156   const double A[6][6] =
 157     {
 158       { L.xx()*L.xx(), L.xy()*L.xy(), L.xz()*L.xz(), s*L.xx()*L.xy(), s*L.xy()*L.xz(), s*L.xz()*L.xx() },
 159       { L.yx()*L.yx(), L.yy()*L.yy(), L.yz()*L.yz(), s*L.yx()*L.yy(), s*L.yy()*L.yz(), s*L.yz()*L.yx() },
 160       { L.zx()*L.zx(), L.zy()*L.zy(), L.zz()*L.zz(), s*L.zx()*L.zy(), s*L.zy()*L.zz(), s*L.zz()*L.zx() },
 161       { s*L.xx()*L.yx(), s*L.xy()*L.yy(), s*L.xz()*L.yz(), (L.xy()*L.yx()+L.xx()*L.yy()), (L.xz()*L.yy()+L.xy()*L.yz()), (L.xx()*L.yz()+L.xz()*L.yx()) },
 162       { s*L.yx()*L.zx(), s*L.yy()*L.zy(), s*L.yz()*L.zz(), (L.yy()*L.zx()+L.yx()*L.zy()), (L.yz()*L.zy()+L.yy()*L.zz()), (L.yx()*L.zz()+L.yz()*L.zx()) },
 163       { s*L.zx()*L.xx(), s*L.zy()*L.xy(), s*L.zz()*L.xz(), (L.zy()*L.xx()+L.zx()*L.xy()), (L.zz()*L.xy()+L.zy()*L.xz()), (L.zx()*L.xz()+L.zz()*L.xx()) }
 164     };
 165
 166   return symmTensor4thOrder
 167     (
 168      // xxxx
 169      A[0][0] * C.xxxx() * A[0][0] +
 170      A[1][0] * C.xxyy() * A[0][0] +
 171      A[2][0] * C.xxzz() * A[0][0] +
 172      A[0][0] * C.xxyy() * A[1][0] +
 173      A[1][0] * C.yyyy() * A[1][0] +
 174      A[2][0] * C.yyzz() * A[1][0] +
 175      A[0][0] * C.xxzz() * A[2][0] +
 176      A[1][0] * C.yyzz() * A[2][0] +
 177      A[2][0] * C.zzzz() * A[2][0] +
 178      A[3][0] * C.xyxy() * A[3][0] +
 179      A[4][0] * C.yzyz() * A[4][0] +
 180      A[5][0] * C.zxzx() * A[5][0],
 181
 182      // xxyy
 183      A[0][0] * C.xxxx() * A[0][1] +
 184      A[1][0] * C.xxyy() * A[0][1] +
 185      A[2][0] * C.xxzz() * A[0][1] +
 186      A[0][0] * C.xxyy() * A[1][1] +
 187      A[1][0] * C.yyyy() * A[1][1] +
 188      A[2][0] * C.yyzz() * A[1][1] +
 189      A[0][0] * C.xxzz() * A[2][1] +
 190      A[1][0] * C.yyzz() * A[2][1] +
 191      A[2][0] * C.zzzz() * A[2][1] +
 192      A[3][0] * C.xyxy() * A[3][1] +
 193      A[4][0] * C.yzyz() * A[4][1] +
 194      A[5][0] * C.zxzx() * A[5][1],
 195
 196      // xzz
 197      A[0][0] * C.xxxx() * A[0][2] +
 198      A[1][0] * C.xxyy() * A[0][2] +
 199      A[2][0] * C.xxzz() * A[0][2] +
 200      A[0][0] * C.xxyy() * A[1][2] +
 201      A[1][0] * C.yyyy() * A[1][2] +
 202      A[2][0] * C.yyzz() * A[1][2] +
 203      A[0][0] * C.xxzz() * A[2][2] +
 204      A[1][0] * C.yyzz() * A[2][2] +
 205      A[2][0] * C.zzzz() * A[2][2] +
 206      A[3][0] * C.xyxy() * A[3][2] +
 207      A[4][0] * C.yzyz() * A[4][2] +
 208      A[5][0] * C.zxzx() * A[5][2],
 209
 210      // yyyy
 211      A[0][1] * C.xxxx() * A[0][1] +
 212      A[1][1] * C.xxyy() * A[0][1] +
 213      A[2][1] * C.xxzz() * A[0][1] +
 214      A[0][1] * C.xxyy() * A[1][1] +
 215      A[1][1] * C.yyyy() * A[1][1] +
 216      A[2][1] * C.yyzz() * A[1][1] +
 217      A[0][1] * C.xxzz() * A[2][1] +
 218      A[1][1] * C.yyzz() * A[2][1] +
 219      A[2][1] * C.zzzz() * A[2][1] +
 220      A[3][1] * C.xyxy() * A[3][1] +
 221      A[4][1] * C.yzyz() * A[4][1] +
 222      A[5][1] * C.zxzx() * A[5][1],
 223
 224      // yyzz
 225      A[0][1] * C.xxxx() * A[0][2] +
 226      A[1][1] * C.xxyy() * A[0][2] +
 227      A[2][1] * C.xxzz() * A[0][2] +
 228      A[0][1] * C.xxyy() * A[1][2] +
 229      A[1][1] * C.yyyy() * A[1][2] +
 230      A[2][1] * C.yyzz() * A[1][2] +
 231      A[0][1] * C.xxzz() * A[2][2] +
 232      A[1][1] * C.yyzz() * A[2][2] +
 233      A[2][1] * C.zzzz() * A[2][2] +
 234      A[3][1] * C.xyxy() * A[3][2] +
 235      A[4][1] * C.yzyz() * A[4][2] +
 236      A[5][1] * C.zxzx() * A[5][2],
 237
 238      // zzzz
 239      A[0][2] * C.xxxx() * A[0][2] +
 240      A[1][2] * C.xxyy() * A[0][2] +
 241      A[2][2] * C.xxzz() * A[0][2] +
 242      A[0][2] * C.xxyy() * A[1][2] +
 243      A[1][2] * C.yyyy() * A[1][2] +
 244      A[2][2] * C.yyzz() * A[1][2] +
 245      A[0][2] * C.xxzz() * A[2][2] +
 246      A[1][2] * C.yyzz() * A[2][2] +
 247      A[2][2] * C.zzzz() * A[2][2] +
 248      A[3][2] * C.xyxy() * A[3][2] +
 249      A[4][2] * C.yzyz() * A[4][2] +
 250      A[5][2] * C.zxzx() * A[5][2],
 251
 252      // xyxy
 253      A[0][3] * C.xxxx() * A[0][3] +
 254      A[1][3] * C.xxyy() * A[0][3] +
 255      A[2][3] * C.xxzz() * A[0][3] +
 256      A[0][3] * C.xxyy() * A[1][3] +
 257      A[1][3] * C.yyyy() * A[1][3] +
 258      A[2][3] * C.yyzz() * A[1][3] +
 259      A[0][3] * C.xxzz() * A[2][3] +
 260      A[1][3] * C.yyzz() * A[2][3] +
 261      A[2][3] * C.zzzz() * A[2][3] +
 262      A[3][3] * C.xyxy() * A[3][3] +
 263      A[4][3] * C.yzyz() * A[4][3] +
 264      A[5][3] * C.zxzx() * A[5][3],
 265
 266      // yzyz
 267      A[0][4] * C.xxxx() * A[0][4] +
 268      A[1][4] * C.xxyy() * A[0][4] +
 269      A[2][4] * C.xxzz() * A[0][4] +
 270      A[0][4] * C.xxyy() * A[1][4] +
 271      A[1][4] * C.yyyy() * A[1][4] +
 272      A[2][4] * C.yyzz() * A[1][4] +
 273      A[0][4] * C.xxzz() * A[2][4] +
 274      A[1][4] * C.yyzz() * A[2][4] +
 275      A[2][4] * C.zzzz() * A[2][4] +
 276      A[3][4] * C.xyxy() * A[3][4] +
 277      A[4][4] * C.yzyz() * A[4][4] +
 278      A[5][4] * C.zxzx() * A[5][4],
 279
 280      // zxzx
 281      A[0][5] * C.xxxx() * A[0][5] +
 282      A[1][5] * C.xxyy() * A[0][5] +
 283      A[2][5] * C.xxzz() * A[0][5] +
 284      A[0][5] * C.xxyy() * A[1][5] +
 285      A[1][5] * C.yyyy() * A[1][5] +
 286      A[2][5] * C.yyzz() * A[1][5] +
 287      A[0][5] * C.xxzz() * A[2][5] +
 288      A[1][5] * C.yyzz() * A[2][5] +
 289      A[2][5] * C.zzzz() * A[2][5] +
 290      A[3][5] * C.xyxy() * A[3][5] +
 291      A[4][5] * C.yzyz() * A[4][5] +
 292      A[5][5] * C.zxzx() * A[5][5]
 293      );
 294 }
 295
 296
 297 template<class Cmpt>
 298 inline DiagTensor<Cmpt> transform
 299 (
 300  const tensor& tt,
 301  const DiagTensor<Cmpt>& st
 302  )
 303 {
 304   notImplemented("transform.H\n"
 305                  "template<>\n"
 306                  "inline DiagTensor<Cmpt> transform\n"
 307                  "(\n"
 308                  "const tensor& tt,\n"
 309                  "const DiagTensor<Cmpt>& st\n"
 310                  ")\n"
 311                  "not implemented");
 312
 313   return st;
 314 }
 315
 316
 317 template<class Cmpt>
 318 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
 319 {
 320     return SymmTensor<Cmpt>
 321     (
 322         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.xx()
 323       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.xy()
 324       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.xz(),
 325
 326         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.yx()
 327       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.yy()
 328       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.yz(),
 329
 330         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.zx()
 331       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.zy()
 332       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.zz(),
 333
 334         (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.yx()
 335       + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.yy()
 336       + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.yz(),
 337
 338         (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
 339       + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
 340       + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
 341
 342         (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
 343       + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
 344       + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
 345     );
 346 }
 347
 348
 349 template<class Type1, class Type2>
 350 inline Type1 transformMask(const Type2& t)
 351 {
 352     return t;
 353 }
 354
 355
 356 template<>
 357 inline sphericalTensor transformMask<sphericalTensor>(const tensor& t)
 358 {
 359     return sph(t);
 360 }
 361
 362
 363 template<>
 364 inline symmTensor transformMask<symmTensor>(const tensor& t)
 365 {
 366     return symm(t);
 367 }
 368
 369
 370 template<>
 371 inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const tensor& t)
 372 {
 373   notImplemented("transform.H\n"
 374                  "template<>\n"
 375                  "inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const tensor& t)\n"
 376                  "not implemented");
 377
 378   return symmTensor4thOrder::zero;
 379 }
 380
 381
 382 template<>
 383 inline diagTensor transformMask<diagTensor>(const tensor& t)
 384 {
 385   return diagTensor(t.xx(), t.yy(), t.zz());
 386 }
 387
 388
 389 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
 390 //  Is guaranteed to return increasing number but is not correct
 391 //  angle. Used for sorting angles.  All input vectors need to be normalized.
 392 //
 393 // Calculates scalar which increases with angle going from e0 to vec in
 394 // the coordinate system e0, e1, e0^e1
 395 //
 396 // Jumps from 2PI -> 0 at -SMALL so parallel vectors with small rounding errors
 397 // should hopefully still get the same quadrant.
 398 //
 399 inline scalar pseudoAngle
 400 (
 401     const vector& e0,
 402     const vector& e1,
 403     const vector& vec
 404 )
 405 {
 406     scalar cos = vec & e0;
 407     scalar sin = vec & e1;
 408
 409     if (sin < -SMALL)
 410     {
 411         return (3.0 + cos)*mathematicalConstant::piByTwo;
 412     }
 413     else
 414     {
 415         return (1.0 - cos)*mathematicalConstant::piByTwo;
 416     }
 417 }
 418
 419
 420 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 421
 422 } // End namespace Foam
 423
 424 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 425
 426 #endif
 427
 428 // ************************************************************************* //