src/tetDecompositionFiniteElement/tetFemMatrix/tetFemMatrix.C

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