src/tetDecompositionFiniteElement/tetFemMatrix/tetFemMatrix.H

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