src/tetFiniteElement/tetFemMatrix/tetFemMatrix.C

   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 Description
  25     Tetrahedral Finite Element matrix member functions and operators
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "PstreamReduceOps.H"
  30
  31 #include "tetFemMatrix.H"
  32 #include "tetPointFields.H"
  33
  34 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  35
  36 namespace Foam
  37 {
  38
  39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  40
  41 template<class Type>
  42 const label tetFemMatrix<Type>::fixFillIn = 4;
  43
  44
  45 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  46
  47 template<class Type>
  48 tetFemMatrix<Type>::tetFemMatrix
  49 (
  50     GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi,
  51     const dimensionSet& ds
  52 )
  53 :
  54     lduMatrix(psi.mesh()),
  55     psi_(psi),
  56     dimensions_(ds),
  57     source_(psi.size(), pTraits<Type>::zero),
  58     boundaryConditionsSet_(false),
  59     fixedEqns_(psi.mesh().lduAddr().size()/fixFillIn),
  60     solvingComponent(0)
  61 {
  62     if (debug)
  63     {
  64         Info<< "tetFemMatrix<Type>(GeometricField<Type, tetPolyPatchField, "
  65             << "tetPointMesh>&, const dimensionSet&) : "
  66             << "constructing tetFemMatrix<Type> for field " << psi_.name()
  67             << endl;
  68     }
  69 }
  70
  71
  72 template<class Type>
  73 tetFemMatrix<Type>::tetFemMatrix(const tetFemMatrix<Type>& tetFem)
  74 :
  75     refCount(),
  76     lduMatrix(tetFem),
  77     psi_(tetFem.psi_),
  78     dimensions_(tetFem.dimensions_),
  79     source_(tetFem.source_),
  80     boundaryConditionsSet_(false),
  81     fixedEqns_(psi_.mesh().lduAddr().size()/fixFillIn),
  82     solvingComponent(0)
  83 {
  84     if (debug)
  85     {
  86         Info<< "tetFemMatrix<Type>::tetFemMatrix(const tetFemMatrix<Type>&) : "
  87             << "copying tetFemMatrix<Type> for field " << psi_.name()
  88             << endl;
  89     }
  90 }
  91
  92
  93 template<class Type>
  94 tetFemMatrix<Type>::tetFemMatrix
  95 (
  96     GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi,
  97     Istream& is
  98 )
  99 :
 100     lduMatrix(psi.mesh()),
 101     psi_(psi),
 102     dimensions_(is),
 103     source_(is),
 104     boundaryConditionsSet_(false),
 105     fixedEqns_(psi.mesh().lduAddr().size()/fixFillIn),
 106     solvingComponent(0)
 107 {
 108     if (debug)
 109     {
 110         Info<< "tetFemMatrix<Type>(GeometricField<Type, tetPolyPatchField, "
 111             << "tetPointMesh>&, Istream&) : "
 112             << "constructing tetFemMatrix<Type> for field " << psi_.name()
 113             << endl;
 114     }
 115 }
 116
 117
 118 template<class Type>
 119 tetFemMatrix<Type>::~tetFemMatrix()
 120 {
 121     if (debug)
 122     {
 123         Info<< "tetFemMatrix<Type>::~tetFemMatrix<Type>() : "
 124             << "destroying tetFemMatrix<Type> for field " << psi_.name()
 125             << endl;
 126     }
 127 }
 128
 129
 130 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 131
 132 // Does the matrix need a reference level for solution
 133 // template<class Type>
 134 // bool tetFemMatrix<Type>::needReference()
 135 // {
 136 //     // Search all boundary conditions, if any are
 137 //     // fixed-value or mixed (Robin) do not set reference level for solution.
 138
 139 //     static bool searched = false;
 140 //     static bool needRef = true;
 141
 142 //     if (!searched)
 143 //     {
 144 //         const BoundaryField<Type>& patchFields = psi_.boundaryField();
 145
 146 //         forAll(patchFields, patchi)
 147 //         {
 148 //             if (patchFields[patchi].fixesValue())
 149 //             {
 150 //                 needRef = false;
 151 //             }
 152 //         }
 153
 154 //         reduce(needRef, andOp<bool>());
 155
 156 //         searched = true;
 157 //     }
 158
 159 //     return needRef;
 160 // }
 161
 162
 163 // Set constrained equation
 164 template<class Type>
 165 void tetFemMatrix<Type>::addConstraint
 166 (
 167     const label vertex,
 168     const Type& value
 169 )
 170 {
 171     constraint<Type> cp(vertex, value);
 172
 173     if (!fixedEqns_.found(vertex))
 174     {
 175         fixedEqns_.insert(vertex, cp);
 176     }
 177     else
 178     {
 179         fixedEqns_[vertex].combine(cp);
 180     }
 181 }
 182
 183
 184 template<class Type>
 185 void tetFemMatrix<Type>::relax(const scalar alpha)
 186 {
 187     if (alpha <= 0)
 188     {
 189         return;
 190     }
 191
 192     Field<Type>& S = source();
 193     scalarField& D = diag();
 194
 195     // Store the current unrelaxed diagonal for use in updating the source
 196     scalarField D0(D);
 197
 198     // Calculate the sum-mag off-diagonal from the interior faces
 199     scalarField sumOff(D.size(), 0.0);
 200     sumMagOffDiag(sumOff);
 201
 202     // Some treatment of coupled boundaries may be added here
 203     // HJ, 3/Sep/2007
 204
 205     // Ensure the matrix is diagonally dominant...
 206     max(D, D, sumOff);
 207
 208     // ... then relax
 209     D /= alpha;
 210
 211     S += (D - D0)*psi_.internalField();
 212 }
 213
 214
 215 template<class Type>
 216 void tetFemMatrix<Type>::relax()
 217 {
 218     scalar alpha = 0;
 219
 220     if (psi_.mesh().solutionDict().relax(psi_.name()))
 221     {
 222         alpha = psi_.mesh().solutionDict().relaxationFactor(psi_.name());
 223     }
 224
 225     if (alpha > 0)
 226     {
 227         relax(alpha);
 228     }
 229 }
 230
 231
 232 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 233
 234 template<class Type>
 235 void tetFemMatrix<Type>::operator=(const tetFemMatrix<Type>& tetFem)
 236 {
 237     if (this == &tetFem)
 238     {
 239         FatalErrorIn
 240         (
 241             "tetFemMatrix<Type>::operator=(const tetFemMatrix<Type>&)"
 242         )   << "attempted assignment to self"
 243             << abort(FatalError);
 244     }
 245
 246     if (&psi_ != &(tetFem.psi_))
 247     {
 248         FatalErrorIn
 249         (
 250             "tetFemMatrix<Type>::operator=(const tetFemMatrix<Type>&)"
 251         )   << "different fields"
 252             << abort(FatalError);
 253     }
 254
 255     lduMatrix::operator=(tetFem);
 256     source_ = tetFem.source_;
 257     boundaryConditionsSet_ = false;
 258     fixedEqns_.clear();
 259 }
 260
 261
 262 template<class Type>
 263 void tetFemMatrix<Type>::operator=(const tmp<tetFemMatrix<Type> >& ttetFem)
 264 {
 265     operator=(ttetFem());
 266     ttetFem.clear();
 267 }
 268
 269
 270 template<class Type>
 271 void tetFemMatrix<Type>::negate()
 272 {
 273     lduMatrix::negate();
 274     source_.negate();
 275 }
 276
 277
 278 template<class Type>
 279 void tetFemMatrix<Type>::operator+=(const tetFemMatrix<Type>& tetFem)
 280 {
 281     checkMethod(*this, tetFem, "+=");
 282
 283     dimensions_ += tetFem.dimensions_;
 284     lduMatrix::operator+=(tetFem);
 285     source_ += tetFem.source_;
 286 }
 287
 288
 289 template<class Type>
 290 void tetFemMatrix<Type>::operator+=(const tmp<tetFemMatrix<Type> >& ttetFem)
 291 {
 292     operator+=(ttetFem());
 293     ttetFem.clear();
 294 }
 295
 296
 297 template<class Type>
 298 void tetFemMatrix<Type>::operator+=
 299 (
 300     const GeometricField<Type, elementPatchField, elementMesh>& su
 301 )
 302 {
 303     checkMethod(*this, su, "+=");
 304     source() -= distributeSource(su.internalField());
 305 }
 306
 307
 308 template<class Type>
 309 void tetFemMatrix<Type>::operator+=
 310 (
 311     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 312 )
 313 {
 314     operator+=(tsu());
 315     tsu.clear();
 316 }
 317
 318
 319 template<class Type>
 320 void tetFemMatrix<Type>::operator-=(const tetFemMatrix<Type>& tetFem)
 321 {
 322     checkMethod(*this, tetFem, "+=");
 323
 324     dimensions_ -= tetFem.dimensions_;
 325     lduMatrix::operator-=(tetFem);
 326     source_ -= tetFem.source_;
 327 }
 328
 329
 330 template<class Type>
 331 void tetFemMatrix<Type>::operator-=(const tmp<tetFemMatrix<Type> >& ttetFem)
 332 {
 333     operator-=(ttetFem());
 334     ttetFem.clear();
 335 }
 336
 337
 338 template<class Type>
 339 void tetFemMatrix<Type>::operator-=
 340 (
 341     const GeometricField<Type, elementPatchField, elementMesh>& su
 342 )
 343 {
 344     checkMethod(*this, su, "-=");
 345     source() += distributeSource(su.internalField());
 346 }
 347
 348
 349 template<class Type>
 350 void tetFemMatrix<Type>::operator-=
 351 (
 352     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 353 )
 354 {
 355     operator-=(tsu());
 356     tsu.clear();
 357 }
 358
 359
 360 template<class Type>
 361 void tetFemMatrix<Type>::operator+=
 362 (
 363     const dimensioned<Type>& su
 364 )
 365 {
 366     checkMethod(*this, su, "+=");
 367     source() -= distributeSource(Field<Type>(psi_.mesh().nCells(), su.value()));
 368 }
 369
 370
 371 template<class Type>
 372 void tetFemMatrix<Type>::operator-=
 373 (
 374     const dimensioned<Type>& su
 375 )
 376 {
 377     checkMethod(*this, su, "-=");
 378     source() += distributeSource(Field<Type>(psi_.mesh().nCells(), su.value()));
 379 }
 380
 381
 382 template<class Type>
 383 void tetFemMatrix<Type>::operator*=
 384 (
 385     const dimensioned<scalar>& ds
 386 )
 387 {
 388     dimensions_ *= ds.dimensions();
 389     lduMatrix::operator*=(ds.value());
 390     source_ *= ds.value();
 391 }
 392
 393
 394 // * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
 395
 396 template<class Type>
 397 void checkMethod
 398 (
 399     const tetFemMatrix<Type>& tetFem1,
 400     const tetFemMatrix<Type>& tetFem2,
 401     const char* op
 402 )
 403 {
 404     if (&tetFem1.psi() != &tetFem2.psi())
 405     {
 406         FatalErrorIn
 407         (
 408             "checkMethod(const tetFemMatrix<Type>&, "
 409             "const tetFemMatrix<Type>&) : "
 410         )   << "incompatible fields for operation "
 411             << endl << "    "
 412             << "[" << tetFem1.psi().name() << "] "
 413             << op
 414             << " [" << tetFem1.psi().name() << "]"
 415             << abort(FatalError);
 416     }
 417
 418     if (dimensionSet::debug && tetFem1.dimensions() != tetFem2.dimensions())
 419     {
 420         FatalErrorIn
 421         (
 422             "checkMethod(const tetFemMatrix<Type>&, "
 423             "const tetFemMatrix<Type>&) : "
 424         )   << "incompatible dimensions for operation "
 425             << endl << "    "
 426             << "[" << tetFem1.psi().name() << tetFem1.dimensions()/dimVolume
 427             << " ] "
 428             << op
 429             << " [" << tetFem1.psi().name() << tetFem2.dimensions()/dimVolume
 430             << " ]"
 431             << abort(FatalError);
 432     }
 433 }
 434
 435
 436 template<class Type>
 437 void checkMethod
 438 (
 439     const tetFemMatrix<Type>& tetFem,
 440     const GeometricField<Type, elementPatchField, elementMesh>& vf,
 441     const char* op
 442 )
 443 {
 444     if
 445     (
 446         dimensionSet::debug
 447      && tetFem.dimensions() != vf.dimensions()
 448     )
 449     {
 450         FatalErrorIn
 451         (
 452             "checkMethod(const tetFemMatrix<Type>&, "
 453             "const GeometricField<Type, elementPatchField, "
 454             "elementMesh>&) : "
 455         )   << "incompatible dimensions for operation "
 456             << endl << "    "
 457             << "[" << tetFem.psi().name() << tetFem.dimensions()/dimVolume
 458             << " ] "
 459             << op
 460             << " [" << tetFem.psi().name() << vf.dimensions() << " ]"
 461             << abort(FatalError);
 462     }
 463 }
 464
 465
 466 template<class Type>
 467 void checkMethod
 468 (
 469     const tetFemMatrix<Type>& tetFem,
 470     const dimensioned<Type>& dt,
 471     const char* op
 472 )
 473 {
 474     if
 475     (
 476         dimensionSet::debug
 477      && tetFem.dimensions() != dt.dimensions()
 478     )
 479     {
 480         FatalErrorIn
 481         (
 482             "checkMethod(const tetFemMatrix<Type>&, "
 483             "const dimensioned<Type>&) : "
 484         )   << "incompatible dimensions for operation "
 485             << endl << "    "
 486             << "[" << tetFem.psi().name() << tetFem.dimensions()/dimVolume
 487             << " ] "
 488             << op
 489             << " [" << dt.name() << dt.dimensions() << " ]"
 490             << abort(FatalError);
 491     }
 492 }
 493
 494
 495 template<class Type>
 496 lduSolverPerformance solve
 497 (
 498     tetFemMatrix<Type>& tetFem,
 499     Istream& solverControls
 500 )
 501 {
 502     return tetFem.solve(solverControls);
 503 }
 504
 505 template<class Type>
 506 lduSolverPerformance solve
 507 (
 508     const tmp<tetFemMatrix<Type> >& ttetFem,
 509     Istream& solverControls
 510 )
 511 {
 512     lduSolverPerformance solverPerf =
 513         const_cast<tetFemMatrix<Type>&>(ttetFem()).solve(solverControls);
 514
 515     ttetFem.clear();
 516
 517     return solverPerf;
 518 }
 519
 520
 521 template<class Type>
 522 lduSolverPerformance solve(tetFemMatrix<Type>& tetFem)
 523 {
 524     return tetFem.solve();
 525 }
 526
 527
 528 template<class Type>
 529 lduSolverPerformance solve(const tmp<tetFemMatrix<Type> >& ttetFem)
 530 {
 531     lduSolverPerformance solverPerf =
 532         const_cast<tetFemMatrix<Type>&>(ttetFem()).solve();
 533
 534     ttetFem.clear();
 535     return solverPerf;
 536 }
 537
 538
 539 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 540
 541 template<class Type>
 542 tmp<tetFemMatrix<Type> > operator+
 543 (
 544     const tetFemMatrix<Type>& A,
 545     const tetFemMatrix<Type>& B
 546 )
 547 {
 548     checkMethod(A, B, "+");
 549     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 550     tC() += B;
 551     return tC;
 552 }
 553
 554
 555 template<class Type>
 556 tmp<tetFemMatrix<Type> > operator+
 557 (
 558     const tmp<tetFemMatrix<Type> >& tA,
 559     const tetFemMatrix<Type>& B
 560 )
 561 {
 562     checkMethod(tA(), B, "+");
 563     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 564     tC() += B;
 565     return tC;
 566 }
 567
 568
 569 template<class Type>
 570 tmp<tetFemMatrix<Type> > operator+
 571 (
 572     const tetFemMatrix<Type>& A,
 573     const tmp<tetFemMatrix<Type> >& tB
 574 )
 575 {
 576     checkMethod(A, tB(), "+");
 577     tmp<tetFemMatrix<Type> > tC(tB.ptr());
 578     tC() += A;
 579     return tC;
 580 }
 581
 582
 583 template<class Type>
 584 tmp<tetFemMatrix<Type> > operator+
 585 (
 586     const tmp<tetFemMatrix<Type> >& tA,
 587     const tmp<tetFemMatrix<Type> >& tB
 588 )
 589 {
 590     checkMethod(tA(), tB(), "+");
 591     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 592     tC() += tB();
 593     tB.clear();
 594     return tC;
 595 }
 596
 597
 598 template<class Type>
 599 tmp<tetFemMatrix<Type> > operator-
 600 (
 601     const tetFemMatrix<Type>& A
 602 )
 603 {
 604     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 605     tC().negate();
 606     return tC;
 607 }
 608
 609
 610 template<class Type>
 611 tmp<tetFemMatrix<Type> > operator-
 612 (
 613     const tmp<tetFemMatrix<Type> >& tA
 614 )
 615 {
 616     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 617     tC().negate();
 618     return tC;
 619 }
 620
 621
 622 template<class Type>
 623 tmp<tetFemMatrix<Type> > operator-
 624 (
 625     const tetFemMatrix<Type>& A,
 626     const tetFemMatrix<Type>& B
 627 )
 628 {
 629     checkMethod(A, B, "-");
 630     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 631     tC() -= B;
 632     return tC;
 633 }
 634
 635
 636 template<class Type>
 637 tmp<tetFemMatrix<Type> > operator-
 638 (
 639     const tmp<tetFemMatrix<Type> >& tA,
 640     const tetFemMatrix<Type>& B
 641 )
 642 {
 643     checkMethod(tA(), B, "-");
 644     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 645     tC() -= B;
 646     return tC;
 647 }
 648
 649
 650 template<class Type>
 651 tmp<tetFemMatrix<Type> > operator-
 652 (
 653     const tetFemMatrix<Type>& A,
 654     const tmp<tetFemMatrix<Type> >& tB
 655 )
 656 {
 657     checkMethod(A, tB(), "-");
 658     tmp<tetFemMatrix<Type> > tC(tB.ptr());
 659     tC() -= A;
 660     tC().negate();
 661     return tC;
 662 }
 663
 664
 665 template<class Type>
 666 tmp<tetFemMatrix<Type> > operator-
 667 (
 668     const tmp<tetFemMatrix<Type> >& tA,
 669     const tmp<tetFemMatrix<Type> >& tB
 670 )
 671 {
 672     checkMethod(tA(), tB(), "-");
 673     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 674     tC() -= tB();
 675     tB.clear();
 676     return tC;
 677 }
 678
 679
 680 template<class Type>
 681 tmp<tetFemMatrix<Type> > operator==
 682 (
 683     const tetFemMatrix<Type>& A,
 684     const tetFemMatrix<Type>& B
 685 )
 686 {
 687     checkMethod(A, B, "==");
 688     return (A - B);
 689 }
 690
 691
 692 template<class Type>
 693 tmp<tetFemMatrix<Type> > operator==
 694 (
 695     const tmp<tetFemMatrix<Type> >& tA,
 696     const tetFemMatrix<Type>& B
 697 )
 698 {
 699     checkMethod(tA(), B, "==");
 700     return (tA - B);
 701 }
 702
 703
 704 template<class Type>
 705 tmp<tetFemMatrix<Type> > operator==
 706 (
 707     const tetFemMatrix<Type>& A,
 708     const tmp<tetFemMatrix<Type> >& tB
 709 )
 710 {
 711     checkMethod(A, tB(), "==");
 712     return (A - tB);
 713 }
 714
 715
 716 template<class Type>
 717 tmp<tetFemMatrix<Type> > operator==
 718 (
 719     const tmp<tetFemMatrix<Type> >& tA,
 720     const tmp<tetFemMatrix<Type> >& tB
 721 )
 722 {
 723     checkMethod(tA(), tB(), "==");
 724     return (tA - tB);
 725 }
 726
 727
 728 template<class Type>
 729 tmp<tetFemMatrix<Type> > operator+
 730 (
 731     const tetFemMatrix<Type>& A,
 732     const GeometricField<Type, elementPatchField, elementMesh>& su
 733 )
 734 {
 735     checkMethod(A, su, "+");
 736     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 737     tC() -= su;
 738     return tC;
 739 }
 740
 741 template<class Type>
 742 tmp<tetFemMatrix<Type> > operator+
 743 (
 744     const tmp<tetFemMatrix<Type> >& tA,
 745     const GeometricField<Type, elementPatchField, elementMesh>& su
 746 )
 747 {
 748     checkMethod(tA(), su, "+");
 749     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 750     tC() -= su;
 751     return tC;
 752 }
 753
 754 template<class Type>
 755 tmp<tetFemMatrix<Type> > operator+
 756 (
 757     const tetFemMatrix<Type>& A,
 758     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 759 )
 760 {
 761     checkMethod(A, tsu(), "+");
 762     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 763     tC() -= tsu();
 764     tsu.clear();
 765     return tC;
 766 }
 767
 768
 769 template<class Type>
 770 tmp<tetFemMatrix<Type> > operator+
 771 (
 772     const tmp<tetFemMatrix<Type> >& tA,
 773     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 774 )
 775 {
 776     checkMethod(tA(), tsu(), "+");
 777     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 778     tC() -= tsu();
 779     tsu.clear();
 780     return tC;
 781 }
 782
 783 template<class Type>
 784 tmp<tetFemMatrix<Type> > operator+
 785 (
 786     const GeometricField<Type, elementPatchField, elementMesh>& su,
 787     const tetFemMatrix<Type>& A
 788 )
 789 {
 790     checkMethod(A, su, "+");
 791     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 792     tC() -= su;
 793     return tC;
 794 }
 795
 796 template<class Type>
 797 tmp<tetFemMatrix<Type> > operator+
 798 (
 799     const GeometricField<Type, elementPatchField, elementMesh>& su,
 800     const tmp<tetFemMatrix<Type> >& tA
 801 )
 802 {
 803     checkMethod(tA(), su, "+");
 804     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 805     tC() -= su;
 806     return tC;
 807 }
 808
 809 template<class Type>
 810 tmp<tetFemMatrix<Type> > operator+
 811 (
 812     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu,
 813     const tetFemMatrix<Type>& A
 814 )
 815 {
 816     checkMethod(A, tsu(), "+");
 817     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 818     tC() -= tsu();
 819     tsu.clear();
 820     return tC;
 821 }
 822
 823 template<class Type>
 824 tmp<tetFemMatrix<Type> > operator+
 825 (
 826     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu,
 827     const tmp<tetFemMatrix<Type> >& tA
 828 )
 829 {
 830     checkMethod(tA(), tsu(), "+");
 831     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 832     tC() -= tsu();
 833     tsu.clear();
 834     return tC;
 835 }
 836
 837
 838 template<class Type>
 839 tmp<tetFemMatrix<Type> > operator-
 840 (
 841     const tetFemMatrix<Type>& A,
 842     const GeometricField<Type, elementPatchField, elementMesh>& su
 843 )
 844 {
 845     checkMethod(A, su, "-");
 846     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 847     tC() += su;
 848     return tC;
 849 }
 850
 851 template<class Type>
 852 tmp<tetFemMatrix<Type> > operator-
 853 (
 854     const tmp<tetFemMatrix<Type> >& tA,
 855     const GeometricField<Type, elementPatchField, elementMesh>& su
 856 )
 857 {
 858     checkMethod(tA(), su, "-");
 859     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 860     tC() += su;
 861     return tC;
 862 }
 863
 864 template<class Type>
 865 tmp<tetFemMatrix<Type> > operator-
 866 (
 867     const tetFemMatrix<Type>& A,
 868     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 869 )
 870 {
 871     checkMethod(A, tsu(), "-");
 872     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 873     tC() += tsu();
 874     tsu.clear();
 875     return tC;
 876 }
 877
 878 template<class Type>
 879 tmp<tetFemMatrix<Type> > operator-
 880 (
 881     const tmp<tetFemMatrix<Type> >& tA,
 882     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
 883 )
 884 {
 885     checkMethod(tA(), tsu(), "-");
 886     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 887     tC() += tsu();
 888     tsu.clear();
 889     return tC;
 890 }
 891
 892
 893 template<class Type>
 894 tmp<tetFemMatrix<Type> > operator-
 895 (
 896     const GeometricField<Type, elementPatchField, elementMesh>& su,
 897     const tetFemMatrix<Type>& A
 898 )
 899 {
 900     checkMethod(A, su, "-");
 901     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 902     tC().negate();
 903     tC() -= su;
 904     return tC;
 905 }
 906
 907
 908 template<class Type>
 909 tmp<tetFemMatrix<Type> > operator-
 910 (
 911     const GeometricField<Type, elementPatchField, elementMesh>& su,
 912     const tmp<tetFemMatrix<Type> >& tA
 913 )
 914 {
 915     checkMethod(tA(), su, "-");
 916     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 917     tC().negate();
 918     tC() -= su;
 919     return tC;
 920 }
 921
 922 template<class Type>
 923 tmp<tetFemMatrix<Type> > operator-
 924 (
 925     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu,
 926     const tetFemMatrix<Type>& A
 927 )
 928 {
 929     checkMethod(A, tsu(), "-");
 930     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 931     tC().negate();
 932     tC() -= tsu();
 933     tsu.clear();
 934     return tC;
 935 }
 936
 937
 938 template<class Type>
 939 tmp<tetFemMatrix<Type> > operator-
 940 (
 941     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu,
 942     const tmp<tetFemMatrix<Type> >& tA
 943 )
 944 {
 945     checkMethod(tA(), tsu(), "-");
 946     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 947     tC().negate();
 948     tC() -= tsu();
 949     tsu.clear();
 950     return tC;
 951 }
 952
 953
 954 template<class Type>
 955 tmp<tetFemMatrix<Type> > operator+
 956 (
 957     const tetFemMatrix<Type>& A,
 958     const dimensioned<Type>& su
 959 )
 960 {
 961     checkMethod(A, su, "+");
 962     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 963     tC() -= su;
 964     return tC;
 965 }
 966
 967
 968 template<class Type>
 969 tmp<tetFemMatrix<Type> > operator+
 970 (
 971     const tmp<tetFemMatrix<Type> >& tA,
 972     const dimensioned<Type>& su
 973 )
 974 {
 975     checkMethod(tA(), su, "+");
 976     tmp<tetFemMatrix<Type> > tC(tA.ptr());
 977     tC() -= su;
 978     return tC;
 979 }
 980
 981
 982 template<class Type>
 983 tmp<tetFemMatrix<Type> > operator+
 984 (
 985     const dimensioned<Type>& su,
 986     const tetFemMatrix<Type>& A
 987 )
 988 {
 989     checkMethod(A, su, "+");
 990     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
 991     tC() -= su;
 992     return tC;
 993 }
 994
 995
 996 template<class Type>
 997 tmp<tetFemMatrix<Type> > operator+
 998 (
 999     const dimensioned<Type>& su,
1000     const tmp<tetFemMatrix<Type> >& tA
1001 )
1002 {
1003     checkMethod(tA(), su, "+");
1004     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1005     tC() -= su;
1006     return tC;
1007 }
1008
1009
1010 template<class Type>
1011 tmp<tetFemMatrix<Type> > operator-
1012 (
1013     const tetFemMatrix<Type>& A,
1014     const dimensioned<Type>& su
1015 )
1016 {
1017     checkMethod(A, su, "-");
1018     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1019     tC() += su;
1020     return tC;
1021 }
1022
1023
1024 template<class Type>
1025 tmp<tetFemMatrix<Type> > operator-
1026 (
1027     const tmp<tetFemMatrix<Type> >& tA,
1028     const dimensioned<Type>& su
1029 )
1030 {
1031     checkMethod(tA(), su, "-");
1032     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1033     tC() += su;
1034     return tC;
1035 }
1036
1037
1038 template<class Type>
1039 tmp<tetFemMatrix<Type> > operator-
1040 (
1041     const dimensioned<Type>& su,
1042     const tetFemMatrix<Type>& A
1043 )
1044 {
1045     checkMethod(A, su, "-");
1046     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1047     tC().negate();
1048     tC() -= su;
1049     return tC;
1050 }
1051
1052
1053 template<class Type>
1054 tmp<tetFemMatrix<Type> > operator-
1055 (
1056     const dimensioned<Type>& su,
1057     const tmp<tetFemMatrix<Type> >& tA
1058 )
1059 {
1060     checkMethod(tA(), su, "-");
1061     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1062     tC().negate();
1063     tC() -= su;
1064     return tC;
1065 }
1066
1067
1068 template<class Type>
1069 tmp<tetFemMatrix<Type> > operator==
1070 (
1071     const tetFemMatrix<Type>& A,
1072     const GeometricField<Type, elementPatchField, elementMesh>& su
1073 )
1074 {
1075     checkMethod(A, su, "==");
1076     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1077     tC() += su;
1078     return tC;
1079 }
1080
1081 template<class Type>
1082 tmp<tetFemMatrix<Type> > operator==
1083 (
1084     const tmp<tetFemMatrix<Type> >& tA,
1085     const GeometricField<Type, elementPatchField, elementMesh>& su
1086 )
1087 {
1088     checkMethod(tA(), su, "==");
1089     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1090     tC() += su;
1091     return tC;
1092 }
1093
1094 template<class Type>
1095 tmp<tetFemMatrix<Type> > operator==
1096 (
1097     const tetFemMatrix<Type>& A,
1098     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
1099 )
1100 {
1101     checkMethod(A, tsu(), "==");
1102     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1103     tC() += tsu();
1104     tsu.clear();
1105     return tC;
1106 }
1107
1108 template<class Type>
1109 tmp<tetFemMatrix<Type> > operator==
1110 (
1111     const tmp<tetFemMatrix<Type> >& tA,
1112     const tmp<GeometricField<Type, elementPatchField, elementMesh> >& tsu
1113 )
1114 {
1115     checkMethod(tA(), tsu(), "==");
1116     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1117     tC() += tsu();
1118     tsu.clear();
1119     return tC;
1120 }
1121
1122
1123 template<class Type>
1124 tmp<tetFemMatrix<Type> > operator==
1125 (
1126     const tetFemMatrix<Type>& A,
1127     const dimensioned<Type>& su
1128 )
1129 {
1130     checkMethod(A, su, "==");
1131     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1132     tC() += su;
1133     return tC;
1134 }
1135
1136
1137 template<class Type>
1138 tmp<tetFemMatrix<Type> > operator==
1139 (
1140     const tmp<tetFemMatrix<Type> >& tA,
1141     const dimensioned<Type>& su
1142 )
1143 {
1144     checkMethod(tA(), su, "==");
1145     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1146     tC() += su.value();
1147     return tC;
1148 }
1149
1150
1151 template<class Type>
1152 tmp<tetFemMatrix<Type> > operator*
1153 (
1154     const dimensioned<scalar>& ds,
1155     const tetFemMatrix<Type>& A
1156 )
1157 {
1158     tmp<tetFemMatrix<Type> > tC(new tetFemMatrix<Type>(A));
1159     tC() *= ds;
1160     return tC;
1161 }
1162
1163
1164 template<class Type>
1165 tmp<tetFemMatrix<Type> > operator*
1166 (
1167     const dimensioned<scalar>& ds,
1168     const tmp<tetFemMatrix<Type> >& tA
1169 )
1170 {
1171     tmp<tetFemMatrix<Type> > tC(tA.ptr());
1172     tC() *= ds;
1173     return tC;
1174 }
1175
1176
1177 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
1178
1179 template<class Type>
1180 Ostream& operator<<(Ostream& os, const tetFemMatrix<Type>& tetFem)
1181 {
1182     os  << static_cast<const lduMatrix&>(tetFem) << nl
1183         << tetFem.dimensions_ << nl
1184         << tetFem.source_ << endl;
1185
1186     os.check("Ostream& operator<<(Ostream&, tetFemMatrix<Type>&");
1187
1188     return os;
1189 }
1190
1191
1192 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1193
1194 } // End namespace Foam
1195
1196 // ************************************************************************* //