src/foam/primitives/transform/transform.H

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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
 150 (
 151     const tensor& L,
 152     const SymmTensor4thOrder<Cmpt>& C
 153 )
 154 {
 155     //- represent fourth order tensors in 6x6 notation.  Rotation is given by
 156     //- C_rotated_af = A_ba * C_cd * A_ef
 157     //- where A is a function of L
 158
 159     const scalar s = ::sqrt(2);
 160     const scalar A[6][6] =
 161     {
 162         {
 163             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() },
 164         { 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() },
 165         { 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() },
 166         { 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()) },
 167         { 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()) },
 168         { 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()) }
 169     };
 170
 171     return symmTensor4thOrder
 172     (
 173         // xxxx
 174         A[0][0] * C.xxxx() * A[0][0] +
 175         A[1][0] * C.xxyy() * A[0][0] +
 176         A[2][0] * C.xxzz() * A[0][0] +
 177         A[0][0] * C.xxyy() * A[1][0] +
 178         A[1][0] * C.yyyy() * A[1][0] +
 179         A[2][0] * C.yyzz() * A[1][0] +
 180         A[0][0] * C.xxzz() * A[2][0] +
 181         A[1][0] * C.yyzz() * A[2][0] +
 182         A[2][0] * C.zzzz() * A[2][0] +
 183         A[3][0] * C.xyxy() * A[3][0] +
 184         A[4][0] * C.yzyz() * A[4][0] +
 185         A[5][0] * C.zxzx() * A[5][0],
 186
 187         // xxyy
 188         A[0][0] * C.xxxx() * A[0][1] +
 189         A[1][0] * C.xxyy() * A[0][1] +
 190         A[2][0] * C.xxzz() * A[0][1] +
 191         A[0][0] * C.xxyy() * A[1][1] +
 192         A[1][0] * C.yyyy() * A[1][1] +
 193         A[2][0] * C.yyzz() * A[1][1] +
 194         A[0][0] * C.xxzz() * A[2][1] +
 195         A[1][0] * C.yyzz() * A[2][1] +
 196         A[2][0] * C.zzzz() * A[2][1] +
 197         A[3][0] * C.xyxy() * A[3][1] +
 198         A[4][0] * C.yzyz() * A[4][1] +
 199         A[5][0] * C.zxzx() * A[5][1],
 200
 201         // xzz
 202         A[0][0] * C.xxxx() * A[0][2] +
 203         A[1][0] * C.xxyy() * A[0][2] +
 204         A[2][0] * C.xxzz() * A[0][2] +
 205         A[0][0] * C.xxyy() * A[1][2] +
 206         A[1][0] * C.yyyy() * A[1][2] +
 207         A[2][0] * C.yyzz() * A[1][2] +
 208         A[0][0] * C.xxzz() * A[2][2] +
 209         A[1][0] * C.yyzz() * A[2][2] +
 210         A[2][0] * C.zzzz() * A[2][2] +
 211         A[3][0] * C.xyxy() * A[3][2] +
 212         A[4][0] * C.yzyz() * A[4][2] +
 213         A[5][0] * C.zxzx() * A[5][2],
 214
 215         // yyyy
 216         A[0][1] * C.xxxx() * A[0][1] +
 217         A[1][1] * C.xxyy() * A[0][1] +
 218         A[2][1] * C.xxzz() * A[0][1] +
 219         A[0][1] * C.xxyy() * A[1][1] +
 220         A[1][1] * C.yyyy() * A[1][1] +
 221         A[2][1] * C.yyzz() * A[1][1] +
 222         A[0][1] * C.xxzz() * A[2][1] +
 223         A[1][1] * C.yyzz() * A[2][1] +
 224         A[2][1] * C.zzzz() * A[2][1] +
 225         A[3][1] * C.xyxy() * A[3][1] +
 226         A[4][1] * C.yzyz() * A[4][1] +
 227         A[5][1] * C.zxzx() * A[5][1],
 228
 229         // yyzz
 230         A[0][1] * C.xxxx() * A[0][2] +
 231         A[1][1] * C.xxyy() * A[0][2] +
 232         A[2][1] * C.xxzz() * A[0][2] +
 233         A[0][1] * C.xxyy() * A[1][2] +
 234         A[1][1] * C.yyyy() * A[1][2] +
 235         A[2][1] * C.yyzz() * A[1][2] +
 236         A[0][1] * C.xxzz() * A[2][2] +
 237         A[1][1] * C.yyzz() * A[2][2] +
 238         A[2][1] * C.zzzz() * A[2][2] +
 239         A[3][1] * C.xyxy() * A[3][2] +
 240         A[4][1] * C.yzyz() * A[4][2] +
 241         A[5][1] * C.zxzx() * A[5][2],
 242
 243         // zzzz
 244         A[0][2] * C.xxxx() * A[0][2] +
 245         A[1][2] * C.xxyy() * A[0][2] +
 246         A[2][2] * C.xxzz() * A[0][2] +
 247         A[0][2] * C.xxyy() * A[1][2] +
 248         A[1][2] * C.yyyy() * A[1][2] +
 249         A[2][2] * C.yyzz() * A[1][2] +
 250         A[0][2] * C.xxzz() * A[2][2] +
 251         A[1][2] * C.yyzz() * A[2][2] +
 252         A[2][2] * C.zzzz() * A[2][2] +
 253         A[3][2] * C.xyxy() * A[3][2] +
 254         A[4][2] * C.yzyz() * A[4][2] +
 255         A[5][2] * C.zxzx() * A[5][2],
 256
 257         // xyxy
 258         A[0][3] * C.xxxx() * A[0][3] +
 259         A[1][3] * C.xxyy() * A[0][3] +
 260         A[2][3] * C.xxzz() * A[0][3] +
 261         A[0][3] * C.xxyy() * A[1][3] +
 262         A[1][3] * C.yyyy() * A[1][3] +
 263         A[2][3] * C.yyzz() * A[1][3] +
 264         A[0][3] * C.xxzz() * A[2][3] +
 265         A[1][3] * C.yyzz() * A[2][3] +
 266         A[2][3] * C.zzzz() * A[2][3] +
 267         A[3][3] * C.xyxy() * A[3][3] +
 268         A[4][3] * C.yzyz() * A[4][3] +
 269         A[5][3] * C.zxzx() * A[5][3],
 270
 271         // yzyz
 272         A[0][4] * C.xxxx() * A[0][4] +
 273         A[1][4] * C.xxyy() * A[0][4] +
 274         A[2][4] * C.xxzz() * A[0][4] +
 275         A[0][4] * C.xxyy() * A[1][4] +
 276         A[1][4] * C.yyyy() * A[1][4] +
 277         A[2][4] * C.yyzz() * A[1][4] +
 278         A[0][4] * C.xxzz() * A[2][4] +
 279         A[1][4] * C.yyzz() * A[2][4] +
 280         A[2][4] * C.zzzz() * A[2][4] +
 281         A[3][4] * C.xyxy() * A[3][4] +
 282         A[4][4] * C.yzyz() * A[4][4] +
 283         A[5][4] * C.zxzx() * A[5][4],
 284
 285         // zxzx
 286         A[0][5] * C.xxxx() * A[0][5] +
 287         A[1][5] * C.xxyy() * A[0][5] +
 288         A[2][5] * C.xxzz() * A[0][5] +
 289         A[0][5] * C.xxyy() * A[1][5] +
 290         A[1][5] * C.yyyy() * A[1][5] +
 291         A[2][5] * C.yyzz() * A[1][5] +
 292         A[0][5] * C.xxzz() * A[2][5] +
 293         A[1][5] * C.yyzz() * A[2][5] +
 294         A[2][5] * C.zzzz() * A[2][5] +
 295         A[3][5] * C.xyxy() * A[3][5] +
 296         A[4][5] * C.yzyz() * A[4][5] +
 297         A[5][5] * C.zxzx() * A[5][5]
 298     );
 299 }
 300
 301
 302 template<class Cmpt>
 303 inline DiagTensor<Cmpt> transform
 304 (
 305     const tensor& tt,
 306     const DiagTensor<Cmpt>& st
 307 )
 308 {
 309     notImplemented
 310     (
 311         "transform.H\n"
 312         "template<>\n"
 313         "inline DiagTensor<Cmpt> transform\n"
 314         "(\n"
 315         "const tensor& tt,\n"
 316         "const DiagTensor<Cmpt>& st\n"
 317         ")\n"
 318         "not implemented"
 319     );
 320
 321     return st;
 322 }
 323
 324
 325 template<class Cmpt>
 326 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
 327 {
 328     return SymmTensor<Cmpt>
 329     (
 330         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.xx()
 331       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.xy()
 332       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.xz(),
 333
 334         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.yx()
 335       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.yy()
 336       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.yz(),
 337
 338         (tt.xx()*st.xx() + tt.xy()*st.xy() + tt.xz()*st.xz())*tt.zx()
 339       + (tt.xx()*st.xy() + tt.xy()*st.yy() + tt.xz()*st.yz())*tt.zy()
 340       + (tt.xx()*st.xz() + tt.xy()*st.yz() + tt.xz()*st.zz())*tt.zz(),
 341
 342         (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.yx()
 343       + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.yy()
 344       + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.yz(),
 345
 346         (tt.yx()*st.xx() + tt.yy()*st.xy() + tt.yz()*st.xz())*tt.zx()
 347       + (tt.yx()*st.xy() + tt.yy()*st.yy() + tt.yz()*st.yz())*tt.zy()
 348       + (tt.yx()*st.xz() + tt.yy()*st.yz() + tt.yz()*st.zz())*tt.zz(),
 349
 350         (tt.zx()*st.xx() + tt.zy()*st.xy() + tt.zz()*st.xz())*tt.zx()
 351       + (tt.zx()*st.xy() + tt.zy()*st.yy() + tt.zz()*st.yz())*tt.zy()
 352       + (tt.zx()*st.xz() + tt.zy()*st.yz() + tt.zz()*st.zz())*tt.zz()
 353     );
 354 }
 355
 356
 357 template<class Type1, class Type2>
 358 inline Type1 transformMask(const Type2& t)
 359 {
 360     return t;
 361 }
 362
 363
 364 template<>
 365 inline sphericalTensor transformMask<sphericalTensor>(const tensor& t)
 366 {
 367     return sph(t);
 368 }
 369
 370
 371 template<>
 372 inline symmTensor transformMask<symmTensor>(const tensor& t)
 373 {
 374     return symm(t);
 375 }
 376
 377
 378 template<>
 379 inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const tensor& t)
 380 {
 381   notImplemented("transform.H\n"
 382          "template<>\n"
 383          "inline symmTensor4thOrder transformMask<symmTensor4thOrder>(const tensor& t)\n"
 384          "not implemented");
 385
 386   return symmTensor4thOrder::zero;
 387 }
 388
 389
 390 template<>
 391 inline diagTensor transformMask<diagTensor>(const tensor& t)
 392 {
 393   return diagTensor(t.xx(), t.yy(), t.zz());
 394 }
 395
 396
 397 //- Estimate angle of vec in coordinate system (e0, e1, e0^e1).
 398 //  Is guaranteed to return increasing number but is not correct
 399 //  angle. Used for sorting angles.  All input vectors need to be normalized.
 400 //
 401 // Calculates scalar which increases with angle going from e0 to vec in
 402 // the coordinate system e0, e1, e0^e1
 403 //
 404 // Jumps from 2PI -> 0 at -SMALL so parallel vectors with small rounding errors
 405 // should hopefully still get the same quadrant.
 406 //
 407 inline scalar pseudoAngle
 408 (
 409     const vector& e0,
 410     const vector& e1,
 411     const vector& vec
 412 )
 413 {
 414     scalar cos = vec & e0;
 415     scalar sin = vec & e1;
 416
 417     if (sin < -SMALL)
 418     {
 419         return (3.0 + cos)*mathematicalConstant::piByTwo;
 420     }
 421     else
 422     {
 423         return (1.0 - cos)*mathematicalConstant::piByTwo;
 424     }
 425 }
 426
 427
 428 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 429
 430 } // End namespace Foam
 431
 432 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 433
 434 #endif
 435
 436 // ************************************************************************* //