src/foam/primitives/VectorSpace/VectorSpaceI.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "error.H"
  27 #include "products.H"
  28 #include "VectorSpaceM.H"
  29 #include "ops.H"
  30 #include "pTraits.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 namespace Foam
  35 {
  36
  37 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  38
  39 template<class Form, class Cmpt, int nCmpt>
  40 inline VectorSpace<Form, Cmpt, nCmpt>::VectorSpace()
  41 {}
  42
  43
  44 template<class Form, class Cmpt, int nCmpt>
  45 inline VectorSpace<Form, Cmpt, nCmpt>::VectorSpace
  46 (
  47     const VectorSpace<Form, Cmpt, nCmpt>& vs
  48 )
  49 {
  50     VectorSpaceOps<nCmpt,0>::eqOp(*this, vs, eqOp<Cmpt>());
  51 }
  52
  53
  54 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  55
  56 template<class Form, class Cmpt, int nCmpt>
  57 inline label VectorSpace<Form, Cmpt, nCmpt>::size() const
  58 {
  59     return nCmpt;
  60 }
  61
  62
  63 template<class Form, class Cmpt, int nCmpt>
  64 inline const Cmpt& VectorSpace<Form, Cmpt, nCmpt>::component
  65 (
  66     const direction d
  67 ) const
  68 {
  69 #   ifdef FULLDEBUG
  70     if (d >= nCmpt)
  71     {
  72         FatalErrorIn
  73         (
  74             "VectorSpace<Form, Cmpt, nCmpt>::component(direction) const"
  75         )   << "index out of range"
  76             << abort(FatalError);
  77     }
  78 #   endif
  79
  80     return v_[d];
  81 }
  82
  83
  84 template<class Form, class Cmpt, int nCmpt>
  85 inline Cmpt& VectorSpace<Form, Cmpt, nCmpt>::component
  86 (
  87     const direction d
  88 )
  89 {
  90 #   ifdef FULLDEBUG
  91     if (d >= nCmpt)
  92     {
  93         FatalErrorIn("VectorSpace<Form, Cmpt, nCmpt>::component(direction)")
  94             << "index out of range"
  95             << abort(FatalError);
  96     }
  97 #   endif
  98
  99     return v_[d];
 100 }
 101
 102
 103 template<class Form, class Cmpt, int nCmpt>
 104 inline void VectorSpace<Form, Cmpt, nCmpt>::component
 105 (
 106     Cmpt& c,
 107     const direction d
 108 ) const
 109 {
 110 #   ifdef FULLDEBUG
 111     if (d >= nCmpt)
 112     {
 113         FatalErrorIn
 114         (
 115             "VectorSpace<Form, Cmpt, nCmpt>::component(Cmpt&, direction) const"
 116         )   << "index out of range"
 117             << abort(FatalError);
 118     }
 119 #   endif
 120
 121     c = v_[d];
 122 }
 123
 124
 125 template<class Form, class Cmpt, int nCmpt>
 126 inline void VectorSpace<Form, Cmpt, nCmpt>::replace
 127 (
 128     const direction d,
 129     const Cmpt& c
 130 )
 131 {
 132 #   ifdef FULLDEBUG
 133     if (d >= nCmpt)
 134     {
 135         FatalErrorIn
 136         (
 137             "VectorSpace<Form, Cmpt, nCmpt>::"
 138             "replace(direction, const Cmpt&) const"
 139         )   << "index out of range"
 140             << abort(FatalError);
 141     }
 142 #   endif
 143
 144     v_[d] = c;
 145 }
 146
 147
 148 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 149
 150 template<class Form, class Cmpt, int nCmpt>
 151 inline const Cmpt& VectorSpace<Form, Cmpt, nCmpt>::operator[]
 152 (
 153     const direction d
 154 ) const
 155 {
 156 #   ifdef FULLDEBUG
 157     if (d >= nCmpt)
 158     {
 159         FatalErrorIn
 160         (
 161             "VectorSpace<Form, Cmpt, nCmpt>::operator[](direction d) const"
 162         )   << "index out of range"
 163             << abort(FatalError);
 164     }
 165 #   endif
 166
 167     return v_[d];
 168 }
 169
 170
 171 template<class Form, class Cmpt, int nCmpt>
 172 inline Cmpt& VectorSpace<Form, Cmpt, nCmpt>::operator[]
 173 (
 174     const direction d
 175 )
 176 {
 177 #   ifdef FULLDEBUG
 178     if (d >= nCmpt)
 179     {
 180         FatalErrorIn("VectorSpace<Form, Cmpt, nCmpt>::operator[](direction d)")
 181             << "index out of range"
 182             << abort(FatalError);
 183     }
 184 #   endif
 185
 186     return v_[d];
 187 }
 188
 189
 190 template<class Form, class Cmpt, int nCmpt>
 191 inline void VectorSpace<Form, Cmpt, nCmpt>::operator=
 192 (
 193     const VectorSpace<Form, Cmpt, nCmpt>& vs
 194 )
 195 {
 196     VectorSpaceOps<nCmpt,0>::eqOp(*this, vs, eqOp<Cmpt>());
 197 }
 198
 199
 200 template<class Form, class Cmpt, int nCmpt>
 201 inline void VectorSpace<Form, Cmpt, nCmpt>::operator+=
 202 (
 203     const VectorSpace<Form, Cmpt, nCmpt>& vs
 204 )
 205 {
 206     VectorSpaceOps<nCmpt,0>::eqOp(*this, vs, plusEqOp<Cmpt>());
 207 }
 208
 209
 210 template<class Form, class Cmpt, int nCmpt>
 211 inline void VectorSpace<Form, Cmpt, nCmpt>::operator-=
 212 (
 213     const VectorSpace<Form, Cmpt, nCmpt>& vs
 214 )
 215 {
 216     VectorSpaceOps<nCmpt,0>::eqOp(*this, vs, minusEqOp<Cmpt>());
 217 }
 218
 219
 220 template<class Form, class Cmpt, int nCmpt>
 221 inline void VectorSpace<Form, Cmpt, nCmpt>::operator*=
 222 (
 223     const scalar s
 224 )
 225 {
 226     VectorSpaceOps<nCmpt,0>::eqOpS(*this, s, multiplyEqOp2<Cmpt, scalar>());
 227 }
 228
 229
 230 template<class Form, class Cmpt, int nCmpt>
 231 inline void VectorSpace<Form, Cmpt, nCmpt>::operator/=
 232 (
 233     const scalar s
 234 )
 235 {
 236     VectorSpaceOps<nCmpt,0>::eqOpS(*this, s, divideEqOp2<Cmpt, scalar>());
 237 }
 238
 239
 240 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 241
 242 template<class Form, class Cmpt, int nCmpt>
 243 inline Cmpt& setComponent
 244 (
 245     VectorSpace<Form, Cmpt, nCmpt>& vs,
 246     const direction d
 247 )
 248 {
 249     return vs.component(d);
 250 }
 251
 252
 253 template<class Form, class Cmpt, int nCmpt>
 254 inline const Cmpt& component
 255 (
 256     const VectorSpace<Form, Cmpt, nCmpt>& vs,
 257     const direction d
 258 )
 259 {
 260     return vs.component(d);
 261 }
 262
 263
 264 // Powers of a Form
 265 // Equivalent to outer-products between the Form and itself
 266 // Form^0 = 1.0
 267 template<class Form, class Cmpt, int nCmpt>
 268 inline typename powProduct<Form, 0>::type pow
 269 (
 270     const VectorSpace<Form, Cmpt, nCmpt>&,
 271     typename powProduct<Form, 0>::type
 272     = pTraits<typename powProduct<Form, 0>::type>::zero
 273 )
 274 {
 275     return 1.0;
 276 }
 277
 278 // Form^1 = Form
 279 template<class Form, class Cmpt, int nCmpt>
 280 inline typename powProduct<Form, 1>::type pow
 281 (
 282     const VectorSpace<Form, Cmpt, nCmpt>& v,
 283     typename powProduct<Form, 1>::type
 284   = pTraits<typename powProduct<Form, 1>::type>::zero
 285 )
 286 {
 287     return static_cast<const Form&>(v);
 288 }
 289
 290
 291 // Form^2 = sqr(Form)
 292 template<class Form, class Cmpt, int nCmpt>
 293 inline typename powProduct<Form, 2>::type pow
 294 (
 295     const VectorSpace<Form, Cmpt, nCmpt>& v,
 296     typename powProduct<Form, 2>::type
 297   = pTraits<typename powProduct<Form, 2>::type>::zero
 298 )
 299 {
 300     return sqr(static_cast<const Form&>(v));
 301 }
 302
 303
 304 template<class Form, class Cmpt, int nCmpt>
 305 inline scalar magSqr
 306 (
 307     const VectorSpace<Form, Cmpt, nCmpt>& vs
 308 )
 309 {
 310     scalar ms = magSqr(vs.v_[0]);
 311     VectorSpaceOps<nCmpt,1>::SeqOp(ms, vs, plusEqMagSqrOp2<scalar, Cmpt>());
 312     return ms;
 313 }
 314
 315
 316 template<class Form, class Cmpt, int nCmpt>
 317 inline scalar mag
 318 (
 319     const VectorSpace<Form, Cmpt, nCmpt>& vs
 320 )
 321 {
 322     return ::sqrt(magSqr(static_cast<const Form&>(vs)));
 323 }
 324
 325
 326 template<class Form, class Cmpt, int nCmpt>
 327 inline VectorSpace<Form, Cmpt, nCmpt> cmptMultiply
 328 (
 329     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 330     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 331 )
 332 {
 333     Form v;
 334     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, cmptMultiplyOp<Cmpt>());
 335     return v;
 336 }
 337
 338
 339 template<class Form, class Cmpt, int nCmpt>
 340 inline VectorSpace<Form, Cmpt, nCmpt> cmptDivide
 341 (
 342     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 343     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 344 )
 345 {
 346     Form v;
 347     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, cmptDivideOp<Cmpt>());
 348     return v;
 349 }
 350
 351
 352 template<class Form, class Cmpt, int nCmpt>
 353 inline Cmpt cmptSumMultiply
 354 (
 355     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 356     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 357 )
 358 {
 359     return cmptSum(cmptMultiply(vs1, vs2));
 360 }
 361
 362
 363 template<class Form, class Cmpt, int nCmpt>
 364 inline VectorSpace<Form, Cmpt, nCmpt> stabilise
 365 (
 366     const VectorSpace<Form, Cmpt, nCmpt>& vs,
 367     const Cmpt& small
 368 )
 369 {
 370     Form v;
 371     VectorSpaceOps<nCmpt,0>::opVS(v, vs, small, stabiliseOp<Cmpt>());
 372     return v;
 373 }
 374
 375
 376 template<class Form, class Cmpt, int nCmpt>
 377 inline Cmpt cmptMax
 378 (
 379     const VectorSpace<Form, Cmpt, nCmpt>& vs
 380 )
 381 {
 382     Cmpt cMax = vs.v_[0];
 383     VectorSpaceOps<nCmpt,1>::SeqOp(cMax, vs, maxEqOp<Cmpt>());
 384     return cMax;
 385 }
 386
 387
 388 template<class Form, class Cmpt, int nCmpt>
 389 inline Cmpt cmptMin
 390 (
 391     const VectorSpace<Form, Cmpt, nCmpt>& vs
 392 )
 393 {
 394     Cmpt cMin = vs.v_[0];
 395     VectorSpaceOps<nCmpt,1>::SeqOp(cMin, vs, minEqOp<Cmpt>());
 396     return cMin;
 397 }
 398
 399
 400 template<class Form, class Cmpt, int nCmpt>
 401 inline Cmpt cmptSum
 402 (
 403     const VectorSpace<Form, Cmpt, nCmpt>& vs
 404 )
 405 {
 406     Cmpt sum = vs.v_[0];
 407     VectorSpaceOps<nCmpt,1>::SeqOp(sum, vs, plusEqOp<Cmpt>());
 408     return sum;
 409 }
 410
 411
 412 template<class Form, class Cmpt, int nCmpt>
 413 inline Cmpt cmptAv
 414 (
 415     const VectorSpace<Form, Cmpt, nCmpt>& vs
 416 )
 417 {
 418     return cmptSum(vs)/nCmpt;
 419 }
 420
 421
 422 template<class Form, class Cmpt, int nCmpt>
 423 inline Form cmptMag
 424 (
 425     const VectorSpace<Form, Cmpt, nCmpt>& vs
 426 )
 427 {
 428     Form v;
 429     VectorSpaceOps<nCmpt,0>::eqOp(v, vs, eqMagOp<Cmpt>());
 430     return v;
 431 }
 432
 433
 434 template<class Form, class Cmpt, int nCmpt>
 435 inline Form cmptStabilise
 436 (
 437     const VectorSpace<Form, Cmpt, nCmpt>& vs,
 438     const scalar small,
 439     const scalar value
 440 )
 441 {
 442     Form v = vs;
 443     for (int i=0; i<nCmpt; i++)
 444     {
 445         if (mag(vs.v_[i]) < small)
 446         {
 447             v.v_[i] = sign(vs.v_[i])*value;
 448         }
 449     }
 450     return v;
 451 }
 452
 453
 454 template<class Form, class Cmpt, int nCmpt>
 455 inline Form max
 456 (
 457     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 458     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 459 )
 460 {
 461     Form v;
 462     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, maxOp<Cmpt>());
 463     return v;
 464 }
 465
 466
 467 template<class Form, class Cmpt, int nCmpt>
 468 inline Form min
 469 (
 470     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 471     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 472 )
 473 {
 474     Form v;
 475     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, minOp<Cmpt>());
 476     return v;
 477 }
 478
 479
 480 template<class Form, class Cmpt, int nCmpt>
 481 inline Form minMod
 482 (
 483     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 484     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 485 )
 486 {
 487     Form v;
 488     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, minModOp<Cmpt>());
 489     return v;
 490 }
 491
 492
 493 template<class Type>
 494 inline Type dot(const scalar s, const Type& t)
 495 {
 496     return s * t;
 497 }
 498
 499
 500 template<class Type>
 501 inline Type dot(const Type& t, const scalar s)
 502 {
 503     return t * s;
 504 }
 505
 506
 507 template
 508 <
 509     class Form1, class Cmpt1, int nCmpt1,
 510     class Form2, class Cmpt2, int nCmpt2
 511 >
 512 inline typename innerProduct<Form1, Form2>::type dot
 513 (
 514     const VectorSpace<Form1, Cmpt1, nCmpt1>& t1,
 515     const VectorSpace<Form2, Cmpt2, nCmpt2>& t2
 516 )
 517 {
 518     return static_cast<const Form1&>(t1) & static_cast<const Form2&>(t2);
 519 }
 520
 521
 522 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 523
 524 template<class Form, class Cmpt, int nCmpt>
 525 inline Form operator-
 526 (
 527     const VectorSpace<Form, Cmpt, nCmpt>& vs
 528 )
 529 {
 530     Form v;
 531     VectorSpaceOps<nCmpt,0>::eqOp(v, vs, eqMinusOp<Cmpt>());
 532     return v;
 533 }
 534
 535
 536 template<class Form, class Cmpt, int nCmpt>
 537 inline Form operator+
 538 (
 539     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 540     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 541 )
 542 {
 543     Form v;
 544     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, plusOp<Cmpt>());
 545     return v;
 546 }
 547
 548
 549 template<class Form, class Cmpt, int nCmpt>
 550 inline Form operator-
 551 (
 552     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 553     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 554 )
 555 {
 556     Form v;
 557     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, minusOp<Cmpt>());
 558     return v;
 559 }
 560
 561
 562 template<class Form, class Cmpt, int nCmpt>
 563 inline Form operator*
 564 (
 565     scalar s,
 566     const VectorSpace<Form, Cmpt, nCmpt>& vs
 567 )
 568 {
 569     Form v;
 570     VectorSpaceOps<nCmpt,0>::opSV(v, s, vs, multiplyOp3<Cmpt, scalar, Cmpt>());
 571     return v;
 572 }
 573
 574
 575 template<class Form, class Cmpt, int nCmpt>
 576 inline Form operator*
 577 (
 578     const VectorSpace<Form, Cmpt, nCmpt>& vs,
 579     scalar s
 580 )
 581 {
 582     Form v;
 583     VectorSpaceOps<nCmpt,0>::opVS(v, vs, s, multiplyOp3<Cmpt, Cmpt, scalar>());
 584     return v;
 585 }
 586
 587
 588 template<class Form, class Cmpt, int nCmpt>
 589 inline Form operator/
 590 (
 591     const VectorSpace<Form, Cmpt, nCmpt>& vs,
 592     scalar s
 593 )
 594 {
 595     Form v;
 596     VectorSpaceOps<nCmpt,0>::opVS(v, vs, s, divideOp3<Cmpt, Cmpt, scalar>());
 597     return v;
 598 }
 599
 600 /*
 601 template<class Form, class Cmpt, int nCmpt>
 602 inline Form operator/
 603 (
 604     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 605     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 606 )
 607 {
 608     Form v;
 609     VectorSpaceOps<nCmpt,0>::op(v, vs1, vs2, divideOp<Cmpt>());
 610     return v;
 611 }
 612
 613
 614 template<class Form, class Cmpt, int nCmpt>
 615 inline Form operator/
 616 (
 617     scalar s,
 618     const VectorSpace<Form, Cmpt, nCmpt>& vs
 619 )
 620 {
 621     Form v;
 622     VectorSpaceOps<nCmpt,0>::opSV(v, s, vs, divideOp2<scalar, Cmpt>());
 623     return v;
 624 }
 625 */
 626
 627
 628 template<class Form, class Cmpt, int nCmpt>
 629 inline Cmpt operator&&
 630 (
 631     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 632     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 633 )
 634 {
 635     Cmpt ddProd = vs1.v_[0]*vs2.v_[0];
 636     for (int i=1; i<nCmpt; ++i)
 637     {
 638         ddProd += vs1.v_[i]*vs2.v_[i];
 639     }
 640     return ddProd;
 641 }
 642
 643
 644 template<class Form, class Cmpt, int nCmpt>
 645 inline bool operator==
 646 (
 647     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 648     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 649 )
 650 {
 651     bool eq = true;
 652     for (int i=0; i<nCmpt; ++i)
 653     {
 654         if (!(eq &= (equal(vs1.v_[i], vs2.v_[i])))) break;
 655     }
 656     return eq;
 657 }
 658
 659
 660 template<class Form, class Cmpt, int nCmpt>
 661 inline bool operator!=
 662 (
 663     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 664     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 665 )
 666 {
 667     return !(vs1 == vs2);
 668 }
 669
 670
 671 template<class Form, class Cmpt, int nCmpt>
 672 inline bool operator>
 673 (
 674     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 675     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 676 )
 677 {
 678     bool gt = true;
 679     for (int i=0; i<nCmpt; ++i)
 680     {
 681         if (!(gt &= vs1.v_[i] > vs2.v_[i])) break;
 682     }
 683     return gt;
 684 }
 685
 686
 687 template<class Form, class Cmpt, int nCmpt>
 688 inline bool operator<
 689 (
 690     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 691     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 692 )
 693 {
 694     bool lt = true;
 695     for (int i=0; i<nCmpt; ++i)
 696     {
 697         if (!(lt &= vs1.v_[i] < vs2.v_[i])) break;
 698     }
 699     return lt;
 700 }
 701
 702
 703 template<class Form, class Cmpt, int nCmpt>
 704 inline bool operator>=
 705 (
 706     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 707     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 708 )
 709 {
 710     return !(vs1 < vs2);
 711 }
 712
 713
 714 template<class Form, class Cmpt, int nCmpt>
 715 inline bool operator<=
 716 (
 717     const VectorSpace<Form, Cmpt, nCmpt>& vs1,
 718     const VectorSpace<Form, Cmpt, nCmpt>& vs2
 719 )
 720 {
 721     return !(vs1 > vs2);
 722 }
 723
 724
 725 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 726
 727 } // End namespace Foam
 728
 729 // ************************************************************************* //