src/tetFiniteElement/tetFemMatrix/tetFemMatrix.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 Class
  25     tetFemMatrix<Type>
  26
  27 Description
  28     Tetrahedral Finite Element matrix.
  29
  30 SourceFiles
  31     tetFemMatrix.C
  32     tetFemMatrixSolve.C
  33
  34 \*---------------------------------------------------------------------------*/
  35
  36 #ifndef tetFemMatrix_H
  37 #define tetFemMatrix_H
  38
  39 #include "tetPointFields.H"
  40 #include "elementFields.H"
  41 #include "lduMatrix.H"
  42 #include "tmp.H"
  43 #include "autoPtr.H"
  44 #include "dimensionedTypes.H"
  45 #include "Map.H"
  46 #include "constraints.H"
  47
  48 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  49
  50 namespace Foam
  51 {
  52
  53 // * * * * * * Forward declaration of template friend fuctions * * * * * * * //
  54
  55 template<class Type>
  56 class tetFemMatrix;
  57
  58 template<class Type>
  59 Ostream& operator<<(Ostream&, const tetFemMatrix<Type>&);
  60
  61
  62 /*---------------------------------------------------------------------------*\
  63                            Class tetFemMatrix Declaration
  64 \*---------------------------------------------------------------------------*/
  65
  66 template<class Type>
  67 class tetFemMatrix
  68 :
  69     public refCount,
  70     public lduMatrix
  71 {
  72     // Private data
  73
  74         //- Reference to GeometricField
  75         GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi_;
  76
  77         //- Dimension set
  78         dimensionSet dimensions_;
  79
  80         //- Source term
  81         Field<Type> source_;
  82
  83         //- Are boundary conditions set?
  84         mutable bool boundaryConditionsSet_;
  85
  86         //- Equation triangle map
  87         mutable Map<constraint<Type> > fixedEqns_;
  88
  89         direction solvingComponent;
  90
  91
  92     // Private static data
  93
  94         //- Matrix constraint fill-in
  95         //  Equals to the estimated fraction of fixed nodes in the matrix
  96         static const label fixFillIn;
  97
  98     // Private Member Functions
  99
 100         //- Add boundary source and diagonal for gradient-type conditions
 101         void addBoundarySourceDiag();
 102
 103         //- Store boundary coefficients
 104         void storeBoundaryCoeffs() const;
 105
 106         //- Add coupling coefficients
 107         //  In parallel computations, the consequence of cell-by-cell
 108         //  matrix assembly is the fact that the "other part" of all the
 109         //  edge coefficients on the coupled edge and the diagonal
 110         //  coefficient is assembled on the neighbouring processor.  The
 111         //  matrix assembly is therefore completed only once the matrix is
 112         //  completed, by adding the contribution from the neighbouring
 113         //  domain.
 114         //  Note: There will be an additional component of
 115         //  vector-matrix multiply, for the edges originating on the
 116         //  parallel patch and protruding into the neighbouring domain.
 117         //  This is treated by doing a "local" multiplication on the
 118         //  neighbouring side and then adding the result.
 119         //  HJ, 13/Nov/2001
 120         void addCouplingCoeffs();
 121
 122         void eliminateCouplingCoeffs();
 123
 124         void addCouplingSource(scalarField&) const;
 125
 126         //- Set component boundary conditions
 127         void setComponentBoundaryConditions
 128         (
 129             const direction,
 130             scalarField& psiCmpt,
 131             scalarField& sourceCmpt
 132         );
 133
 134         //- Reconstruct matrix
 135         void reconstructMatrix();
 136
 137         //- Distribute source.  Given an element field construct a source
 138         //  by distributing the volume integral over the element points
 139         tmp<Field<Type> > distributeSource(const Field<Type>&) const;
 140
 141
 142 public:
 143
 144     // No tetFemSolver: boundary coefficients are not cached
 145     // HJ, 6/Jul/2007
 146
 147     ClassName("tetFemMatrix");
 148
 149     typedef constraint<Type> ConstraintType;
 150
 151
 152     // Constructors
 153
 154         //- Construct given a field to solve for
 155         tetFemMatrix
 156         (
 157             GeometricField<Type, tetPolyPatchField, tetPointMesh>&,
 158             const dimensionSet&
 159         );
 160
 161         //- Construct as copy
 162         tetFemMatrix(const tetFemMatrix<Type>&);
 163
 164         //- Construct from Istream given field to solve for
 165         tetFemMatrix
 166         (
 167             GeometricField<Type, tetPolyPatchField, tetPointMesh>&,
 168             Istream&
 169         );
 170
 171
 172     // Destructor
 173
 174         virtual ~tetFemMatrix();
 175
 176
 177     // Member functions
 178
 179         // Access
 180
 181             const GeometricField<Type, tetPolyPatchField, tetPointMesh>&
 182             psi() const
 183             {
 184                 return psi_;
 185             }
 186
 187             GeometricField<Type, tetPolyPatchField, tetPointMesh>& psi()
 188             {
 189                 return psi_;
 190             }
 191
 192             const dimensionSet& dimensions() const
 193             {
 194                 return dimensions_;
 195             }
 196
 197             Field<Type>& source()
 198             {
 199                 return source_;
 200             }
 201
 202             const Field<Type>& source() const
 203             {
 204                 return source_;
 205             }
 206
 207             //- Does the matrix need a reference level for solution
 208             //bool needReference();
 209
 210
 211         // Operations
 212
 213             //- Set constrained value in a prescribed point
 214             void addConstraint
 215             (
 216                 const label vertex,
 217                 const Type& value
 218             );
 219
 220             //- Relax matrix (for steady-state solution).
 221             //  alpha = 1 : diagonally equal
 222             //  alpha < 1 :    ,,      dominant
 223             //  alpha = 0 : do nothing
 224             void relax(const scalar alpha);
 225
 226             //- Relax matrix (for steady-state solution).
 227             //  alpha is read from controlDict
 228             void relax();
 229
 230             //- Check matrix for diagonal dominance
 231             void check();
 232
 233             //- Solve returning the solution statistics.
 234             //  Convergence tolerance read from dictionary
 235             lduSolverPerformance solve(const dictionary&);
 236
 237             //- Solve returning the solution statistics.
 238             //  Solver controls read from tetFemSolution
 239             lduSolverPerformance solve();
 240
 241             //- Return the matrix residual
 242             tmp<Field<Type> > residual();
 243
 244
 245     // Member operators
 246
 247         void operator=(const tetFemMatrix<Type>&);
 248         void operator=(const tmp<tetFemMatrix<Type> >&);
 249
 250         void negate();
 251
 252         void operator+=(const tetFemMatrix<Type>&);
 253         void operator+=(const tmp<tetFemMatrix<Type> >&);
 254
 255         void operator-=(const tetFemMatrix<Type>&);
 256         void operator-=(const tmp<tetFemMatrix<Type> >&);
 257
 258         void operator+=
 259         (
 260             const GeometricField<Type, elementPatchField, elementMesh>&
 261         );
 262         void operator+=
 263         (
 264             const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 265         );
 266
 267         void operator-=
 268         (
 269             const GeometricField<Type, elementPatchField, elementMesh>&
 270         );
 271
 272         void operator-=
 273         (
 274             const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 275         );
 276
 277         void operator+=(const dimensioned<Type>&);
 278         void operator-=(const dimensioned<Type>&);
 279
 280         void operator*=(const dimensioned<scalar>&);
 281
 282     // Ostream operator
 283
 284         friend Ostream& operator<< <Type>
 285         (
 286             Ostream&,
 287             const tetFemMatrix<Type>&
 288         );
 289 };
 290
 291
 292 // * * * * * * * * * * * * * * * Global functions  * * * * * * * * * * * * * //
 293
 294 template<class Type>
 295 void checkMethod
 296 (
 297     const tetFemMatrix<Type>&,
 298     const tetFemMatrix<Type>&,
 299     const char*
 300 );
 301
 302
 303 template<class Type>
 304 void checkMethod
 305 (
 306     const tetFemMatrix<Type>&,
 307     const GeometricField<Type, elementPatchField, elementMesh>&,
 308     const char*
 309 );
 310
 311
 312 template<class Type>
 313 void checkMethod
 314 (
 315     const tetFemMatrix<Type>&,
 316     const dimensioned<Type>&,
 317     const char*
 318 );
 319
 320
 321 //- Solve returning the solution statistics given convergence tolerance
 322 //  Solver controls read from Istream
 323 template<class Type>
 324 lduSolverPerformance solve(tetFemMatrix<Type>&, Istream&);
 325
 326
 327 //- Solve returning the solution statistics given convergence tolerance,
 328 //  deleting temporary matrix after solution.
 329 //  Solver controls read from Istream
 330 template<class Type>
 331 lduSolverPerformance solve(const tmp<tetFemMatrix<Type> >&, Istream&);
 332
 333
 334 //- Solve returning the solution statistics given convergence tolerance
 335 //  Convergence tolerance read from controlDict
 336 template<class Type>
 337 lduSolverPerformance solve(tetFemMatrix<Type>&);
 338
 339 //- Solve returning the solution statistics given convergence tolerance,
 340 //  deleting temporary matrix after solution.
 341 //  Convergence tolerance read from controlDict
 342 template<class Type>
 343 lduSolverPerformance solve(const tmp<tetFemMatrix<Type> >&);
 344
 345
 346 // * * * * * * * * * * * * * * * Global operators  * * * * * * * * * * * * * //
 347
 348 template<class Type>
 349 tmp<tetFemMatrix<Type> > operator-
 350 (
 351     const tetFemMatrix<Type>&
 352 );
 353
 354 template<class Type>
 355 tmp<tetFemMatrix<Type> > operator-
 356 (
 357     const tmp<tetFemMatrix<Type> >&
 358 );
 359
 360 template<class Type>
 361 tmp<tetFemMatrix<Type> > operator+
 362 (
 363     const tetFemMatrix<Type>&,
 364     const tetFemMatrix<Type>&
 365 );
 366
 367 template<class Type>
 368 tmp<tetFemMatrix<Type> > operator+
 369 (
 370     const tmp<tetFemMatrix<Type> >&,
 371     const tetFemMatrix<Type>&
 372 );
 373
 374 template<class Type>
 375 tmp<tetFemMatrix<Type> > operator+
 376 (
 377     const tetFemMatrix<Type>&,
 378     const tmp<tetFemMatrix<Type> >&
 379 );
 380
 381 template<class Type>
 382 tmp<tetFemMatrix<Type> > operator+
 383 (
 384     const tmp<tetFemMatrix<Type> >&,
 385     const tmp<tetFemMatrix<Type> >&
 386 );
 387
 388 template<class Type>
 389 tmp<tetFemMatrix<Type> > operator-
 390 (
 391     const tetFemMatrix<Type>&,
 392     const tetFemMatrix<Type>&
 393 );
 394
 395 template<class Type>
 396 tmp<tetFemMatrix<Type> > operator-
 397 (
 398     const tmp<tetFemMatrix<Type> >&,
 399     const tetFemMatrix<Type>&
 400 );
 401
 402 template<class Type>
 403 tmp<tetFemMatrix<Type> > operator-
 404 (
 405     const tetFemMatrix<Type>&,
 406     const tmp<tetFemMatrix<Type> >&
 407 );
 408
 409 template<class Type>
 410 tmp<tetFemMatrix<Type> > operator-
 411 (
 412     const tmp<tetFemMatrix<Type> >&,
 413     const tmp<tetFemMatrix<Type> >&
 414 );
 415
 416 template<class Type>
 417 tmp<tetFemMatrix<Type> > operator==
 418 (
 419     const tetFemMatrix<Type>&,
 420     const tetFemMatrix<Type>&
 421 );
 422
 423 template<class Type>
 424 tmp<tetFemMatrix<Type> > operator==
 425 (
 426     const tmp<tetFemMatrix<Type> >&,
 427     const tetFemMatrix<Type>&
 428 );
 429
 430 template<class Type>
 431 tmp<tetFemMatrix<Type> > operator==
 432 (
 433     const tetFemMatrix<Type>&,
 434     const tmp<tetFemMatrix<Type> >&
 435 );
 436
 437 template<class Type>
 438 tmp<tetFemMatrix<Type> > operator==
 439 (
 440     const tmp<tetFemMatrix<Type> >&,
 441     const tmp<tetFemMatrix<Type> >&
 442 );
 443
 444 template<class Type>
 445 tmp<tetFemMatrix<Type> > operator+
 446 (
 447     const tetFemMatrix<Type>&,
 448     const GeometricField<Type, elementPatchField, elementMesh>&
 449 );
 450
 451 template<class Type>
 452 tmp<tetFemMatrix<Type> > operator+
 453 (
 454     const tmp<tetFemMatrix<Type> >&,
 455     const GeometricField<Type, elementPatchField, elementMesh>&
 456 );
 457
 458 template<class Type>
 459 tmp<tetFemMatrix<Type> > operator+
 460 (
 461     const tetFemMatrix<Type>&,
 462     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 463 );
 464
 465 template<class Type>
 466 tmp<tetFemMatrix<Type> > operator+
 467 (
 468     const tmp<tetFemMatrix<Type> >&,
 469     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 470 );
 471
 472 template<class Type>
 473 tmp<tetFemMatrix<Type> > operator+
 474 (
 475     const GeometricField<Type, elementPatchField, elementMesh>&,
 476     const tetFemMatrix<Type>&
 477 );
 478
 479 template<class Type>
 480 tmp<tetFemMatrix<Type> > operator+
 481 (
 482     const GeometricField<Type, elementPatchField, elementMesh>&,
 483     const tmp<tetFemMatrix<Type> >&
 484 );
 485
 486 template<class Type>
 487 tmp<tetFemMatrix<Type> > operator+
 488 (
 489     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&,
 490     const tetFemMatrix<Type>&
 491 );
 492
 493 template<class Type>
 494 tmp<tetFemMatrix<Type> > operator+
 495 (
 496     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&,
 497     const tmp<tetFemMatrix<Type> >&
 498 );
 499
 500 template<class Type>
 501 tmp<tetFemMatrix<Type> > operator-
 502 (
 503     const tetFemMatrix<Type>&,
 504     const GeometricField<Type, elementPatchField, elementMesh>&
 505 );
 506
 507 template<class Type>
 508 tmp<tetFemMatrix<Type> > operator-
 509 (
 510     const tmp<tetFemMatrix<Type> >&,
 511     const GeometricField<Type, elementPatchField, elementMesh>&
 512 );
 513
 514 template<class Type>
 515 tmp<tetFemMatrix<Type> > operator-
 516 (
 517     const tetFemMatrix<Type>&,
 518     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 519 );
 520
 521 template<class Type>
 522 tmp<tetFemMatrix<Type> > operator-
 523 (
 524     const tmp<tetFemMatrix<Type> >&,
 525     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 526 );
 527
 528 template<class Type>
 529 tmp<tetFemMatrix<Type> > operator-
 530 (
 531     const GeometricField<Type, elementPatchField, elementMesh>&,
 532     const tetFemMatrix<Type>&
 533 );
 534
 535 template<class Type>
 536 tmp<tetFemMatrix<Type> > operator-
 537 (
 538     const GeometricField<Type, elementPatchField, elementMesh>&,
 539     const tmp<tetFemMatrix<Type> >&
 540 );
 541
 542 template<class Type>
 543 tmp<tetFemMatrix<Type> > operator-
 544 (
 545     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&,
 546     const tetFemMatrix<Type>&
 547 );
 548
 549 template<class Type>
 550 tmp<tetFemMatrix<Type> > operator-
 551 (
 552     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&,
 553     const tmp<tetFemMatrix<Type> >&
 554 );
 555
 556 template<class Type>
 557 tmp<tetFemMatrix<Type> > operator+
 558 (
 559     const tmp<tetFemMatrix<Type> >&,
 560     const dimensioned<Type>&
 561 );
 562
 563 template<class Type>
 564 tmp<tetFemMatrix<Type> > operator+
 565 (
 566     const dimensioned<Type>&,
 567     const tmp<tetFemMatrix<Type> >&
 568 );
 569
 570 template<class Type>
 571 tmp<tetFemMatrix<Type> > operator-
 572 (
 573     const tmp<tetFemMatrix<Type> >&,
 574     const dimensioned<Type>&
 575 );
 576
 577 template<class Type>
 578 tmp<tetFemMatrix<Type> > operator-
 579 (
 580     const dimensioned<Type>&,
 581     const tmp<tetFemMatrix<Type> >&
 582 );
 583
 584 template<class Type>
 585 tmp<tetFemMatrix<Type> > operator==
 586 (
 587     const tetFemMatrix<Type>&,
 588     const GeometricField<Type, elementPatchField, elementMesh>&
 589 );
 590
 591 template<class Type>
 592 tmp<tetFemMatrix<Type> > operator==
 593 (
 594     const tmp<tetFemMatrix<Type> >&,
 595     const GeometricField<Type, elementPatchField, elementMesh>&
 596 );
 597
 598 template<class Type>
 599 tmp<tetFemMatrix<Type> > operator==
 600 (
 601     const tetFemMatrix<Type>&,
 602     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 603 );
 604
 605 template<class Type>
 606 tmp<tetFemMatrix<Type> > operator==
 607 (
 608     const tmp<tetFemMatrix<Type> >&,
 609     const tmp<GeometricField<Type, elementPatchField, elementMesh> >&
 610 );
 611
 612 template<class Type>
 613 tmp<tetFemMatrix<Type> > operator==
 614 (
 615     const tetFemMatrix<Type>&,
 616     const dimensioned<Type>&
 617 );
 618
 619 template<class Type>
 620 tmp<tetFemMatrix<Type> > operator==
 621 (
 622     const tmp<tetFemMatrix<Type> >&,
 623     const dimensioned<Type>&
 624 );
 625
 626
 627 template<class Type>
 628 tmp<tetFemMatrix<Type> > operator*
 629 (
 630     const dimensioned<scalar>&,
 631     const tetFemMatrix<Type>&
 632 );
 633
 634 template<class Type>
 635 tmp<tetFemMatrix<Type> > operator*
 636 (
 637     const dimensioned<scalar>&,
 638     const tmp<tetFemMatrix<Type> >&
 639 );
 640
 641
 642 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 643
 644 } // End namespace Foam
 645
 646 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 647
 648 #ifdef NoRepository
 649 #   include "tetFemMatrix.C"
 650 #   include "tetFemMatrixSolve.C"
 651 #   include "tetFemMatrixTools.C"
 652 #   include "tetFemMatrixCheck.C"
 653 #endif
 654
 655 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 656
 657 #endif
 658
 659 // ************************************************************************* //