View source for OpenFOAM guide/Matrices in OpenFOAM

Linear equations resulting from the [[discretization]] method are solved in [[OpenFOAM]] using matrix methods.  [[OpenFOAM]] has a specialized strategy for handling these matrices.

== Overview ==

The [[discretization]] process results in the equation:

<math>\mathbf{A} \boldsymbol \Psi = \mathbf B</math>

Where:
* '''A''' is the discretization matrix, which holds the [[matrix coefficients]];
* <math>\boldsymbol \Psi</math> is the vector array of variables being solved for; and
* '''B''' is the discretization source term.

The '''A''' component is an ''N'' x ''N'' square matrix, where ''N'' is the number of equations produced by the [[discretization]] process.  [[OpenFOAM]] has developed a specialized strategy to handle this matrix.  Important components of this strategy include:

* lduMatrix (lower - diagonal - upper - matrix storage class); and
* lduAddressing (an indexing class);
* fvMatrix (contains components specific to the [[finite volume method]]).

== lduMatrix ==
The lduMatrix class stores the [[matrix coefficients]] in three arrays:

 scalarField *lowerPtr_;
 scalarField *diagPtr_;
 scalarField *upperPtr_;

These arrays are Fields that contain all the non-zero entries in the [[matrix coefficients|coefficient matrix]].

The access of the coefficients goes over the functions:

 scalarField& upper()
 scalarField& lower()
 scalarField& diag()

The size of the arrays lowerPtr_ and upperPtr_ are equal to the number of internal faces. The size of the array diagPtr_ is equal to the number of cells. 

=== Diagonals ===
Every cell has one and only one diagonal coefficient, none of which are zero.  Therefore the diagonal list contains ''N'' values, where ''N'' is the number of cells in the mesh.  These are denoted:

* As [[discretization coefficients]], <math>A_p \neq 0\,\!</math>; or
* As [[matrix coefficients]], <math>A_{ii} \neq 0\,\!</math>.

=== Upper and lower triangles ===

The upper and lower triangles have (''N ''- 1)''N'' possible values, however most are zero.  The non-zero values correspond to cell-pairs that influence one another.  [[OpenFOAM]] uses a small computational molecule, therefore only adjacent cells influence one another.  In other words, there are two non-zero coefficients for every ''internal face'' in the mesh: one that appears in the lower triangle; and one that appears in the upper triangle.

For instance, the next figure shows a simple 3 x 3 mesh: 

[[File:FigureGridLduMatrix.png]]

It consists in 9 cell centers, 12 internal faces and 12 external (boundary) faces.  The numbers of the cell centers are shown in black and the number of the faces are shown in blue.

The non-zero coefficients for [[OpenFOAM]] would be:
{| class="wikitable" border="1"
|+ Non-zero coefficients
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| X || N || || || || N || || ||
|- align="center"
! -2-
| O || X || N || || N || || || ||
|- align="center"
! -3-
| || O || X || N || || || || ||
|- align="center"
! -4-
| || || O || X || N || || || || N
|- align="center"
! -5-
| || O || || O || X || N || || N ||
|- align="center"
! -6-
| O || || || || O || X || N || ||
|- align="center"
! -7-
| || || || || || O || X || N ||
|- align="center"
! -8-
| || || || || O || || O || X || N
|- align="center"
! -9-
| || || || O || || || || O || X
|}

{| class="wikitable" border="1"
|+ Non-zero coefficients
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| <math>A_{1,1}</math> || <math>A_{1,2}</math> || || || || <math>A_{1,6}</math> || || ||
|- align="center"
! -2-
| <math>A_{2,1}</math> || <math>A_{2,2}</math> || <math>A_{2,3}</math> || || <math>A_{2,5}</math> || || || ||
|- align="center"
! -3-
| || <math>A_{3,2}</math> || <math>A_{3,3}</math> || <math>A_{3,4}</math> || || || || ||
|- align="center"
! -4-
| || || <math>A_{4,3}</math> || <math>A_{4,4}</math> || <math>A_{4,5}</math> || || || || <math>A_{4,9}</math>
|- align="center"
! -5-
| || <math>A_{5,2}</math> || || <math>A_{5,4}</math> || <math>A_{5,5}</math>|| <math>A_{5,6}</math> || || <math>A_{5,8}</math> ||
|- align="center"
! -6-
| <math>A_{6,1}</math> || || || || <math>A_{6,5}</math> || <math>A_{6,6}</math> || <math>A_{6,7}</math> || ||
|- align="center"
! -7-
| || || || || || <math>A_{7,6}</math> || <math>A_{7,7}</math>|| <math>A_{7,8}</math> ||
|- align="center"
! -8-
| || || || || <math>A_{8,5}</math> || || <math>A_{8,7}</math> || <math>A_{8,8}</math> || <math>A_{8,9}</math>
|- align="center"
! -9-
| || || || <math>A_{9,4}</math> || || || || <math>A_{9,8}</math> || <math>A_{9,9}</math>
|}
Some mathematical operations require distinguishing between the two coefficients as they are sometimes opposite in sign (for instance when fluxes are involved).  [[OpenFOAM]]'s mesh design anticipated this requirement, which is why every face has an ''owner'' cell and a ''neighbour'' cell.  ''Owner'', ''neighbour'' and ''diagonal'' coefficients are defined as:

<math>
A_{ij} =
\begin{cases}
\text{owner}, & i > j \\
\text{diagonal}, & i = j \\
\text{neighbour}, & i < j
\end{cases}
</math>

Where:
* ''i'' is the matrix row number (corresponds to cell index); and
* ''j'' is the matrix column number (corresponds to cell index).

With this definition, the ''lower'' triangle has the ''owner'' coefficients; and the ''upper'' triangle has the ''neighbour'' coefficients.  They are denoted:

* As [[discretization coefficients]], <math>A_{r,i}\,\!</math>, where ''r'' indicates "related cells"<ref>What I'm calling ''related cells'' is conventionally called ''neighbours''.  But [[OpenFOAM]] has a different meaning for ''neighbours'', so the term ''related cells'' is used.</ref>;
* As [[matrix coefficients]]:
** Owner coefficients are: <math>\left \lbrack A_{ij}\right \rbrack_{i>j} = A_o\,\!</math>; and
** Neighbour coefficients are: <math>\left \lbrack A_{ij}\right \rbrack_{i<j} = A_n\,\!</math>.

With the above example the entries of the lowerPtr_, upperPtr_ and diagPtr_ array are:

<math>  lower() =  A_{2,1}, A_{3,2}, A_{6,1}, A_{5,2}, A_{4,3}, A_{6,5}, A_{5,4}, A_{7,6}, A_{8,5}, A_{9,4}, A_{8,7}, A_{9,8} </math>

<math>  upper() =  A_{1,2}, A_{2,3}, A_{1,6}, A_{2,5}, A_{3,4}, A_{5,6}, A_{4,5}, A_{6,7}, A_{5,8}, A_{4,9}, A_{7,8}, A_{8,9} </math>

<math>  diag() =  A_{1,1}, A_{2,2}, A_{3,3}, A_{4,4}, A_{5,5}, A_{6,6}, A_{7,7}, A_{8,8}, A_{9,9} </math>

Note that the ordering of the entries of the arrays lowerPtr_ and upperPtr_ are stored according to the face number and of the array diagPtr_ according to the cell number.  

With this consistently applied across the mesh, a direction for positive surface-normal vectors is defined: a positive surface-normal vector (such as fluxes) points away from the ''owner'' cell, and towards the ''neighbour'' cell.



=== lduAddressing ===
The lower, diagonal and upper coefficients are all scalarFields (basically a list), but they are indexed differently.  The diagonal coefficients are indexed by ''cell index''; whereas the upper and lower triangles are indexed by ''face index''.  In order to get the indices to match, an addressing array is provided for the lower and upper triangles to translate their ''face index'' into a corresponding ''cell index''.  This addressing array is called an ''lduAddressing'' array.

The lduAdressing class (see https://www.openfoam.com/documentation/guides/latest/api/classFoam_1_1lduAddressing.html) contains two function which return the label list containing the owner (lowerAddr()) and neighbor  (upperAddr())  of each face. For our example above the two arrays look like:

 lowerAddr() = 1, 2, 1, 2, 3, 5, 4, 6, 5, 4, 7, 8
 upperAddr() = 2, 3, 6, 5, 4, 6, 5, 7, 8, 9, 8, 9

The above arrays are ordered according to the face number and contain the owner (lowerAddr()) and neighbor (upperAddr()) of the face. 

==== When do you need it? ====
''lduAddressing'' is needed anytime you need to perform an operation on all rows of an <tt>lduMatrix</tt> that involves both ''diagonal'' and ''off-diagonal'' coefficients.

=== How do you use it? ===
To use ''lduAddressing'':
* separate operations involving the lower triangle from operations involving the upper triangle;
* loop over all the faces in the mesh;
* upper and lower triangle coefficients use the loop index - e.g. <tt>upper[facei] = ...</tt>
* diagonal coefficients use the lduAddressing coefficient of the face index:
** <tt>diag[l[facei]] -= lower[facei]</tt>, where <tt>l</tt> is the ''lduAddressing'' for the lower triangle; and
** <tt>diag[u[facei]] -= upper[facei]</tt>, where <tt>u</tt> is the ''lduAddressing'' for the upper triangle.

==== Example ====
For example, the <tt>[[lduMatrix;;negSumDiag|negSumDiag]]</tt> operation is frequently required.  This is where the diagonal is the negative sum of its neighbours:

<math>A_p = -\sum\limits_r A_r</math>

The guts of its code is:

 for (register label face=0; face<l.size(); face++)
 {
     Diag[l[face]] -= Lower[face];
     Diag[u[face]] -= Upper[face];
 }
Where:
* <tt>Diag</tt> is the diagonal coefficient;
* <tt>l</tt> is the ''lduAddressing'' for the lower triangle;
* <tt>u</tt> is the ''lduAddressing'' for the upper triangle;
* <tt>Lower</tt> is the lower triangle; and
* <tt>Upper</tt> is the upper triangle.


=== Diagonal, symmetric, and asymmetric matrices ===
It is mandatory to define the diagonal; the others are not required.  If only one of the lower or upper triangles are defined, the matrix is assumed to be symmetric.  If both are defined, it will assume the matrix is asymmetric, even if they are equal.

{| class="wikitable" border="1"
|-
!  lowerPtr !!  diagPtr !!  upperPtr !!  Result
|-align=center
|  Yes/No || No || Yes/No || Invalid matrix - error thrown
|-align=center
|  No || Yes || No || Diagonal matrix
|-align=center
|  Yes || Yes || No || Symmetric matrix
|-align=center
|  No || Yes || Yes || Symmetric matrix
|-align=center
|  Yes || Yes || Yes || Asymmetric matrix
|}
:'''Caution:''' Accessing the opposite pointer without a <tt>const</tt> modifier will convert the matrix to an asymmetric matrix.  To modify the off-diagonal of a symmetric matrix, first test which pointer is active using <tt>hasUpper()</tt> and <tt>hasLower()</tt>.

=== OpenFOAM code example ===
For example, in calculating the diffusion of a quantity <math>\phi\,\!</math>, the governing equation is:

<math>\boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf U = 0</math>

This is the Laplacian, which is calculated in the [[finite volume method]] using fvm::laplacian.  The [[discretization coefficients]] are:

<math>A_r = - \frac{\Gamma S_f}{\delta}</math>

<math>A_p = \sum\limits_r A_r</math>

To implement this in [[OpenFOAM]]:

 GeometricField<scalar, fvsPatchField, surfaceMesh> gammaMagSf ( gamma*mesh.magSf() );
 fvMatrix<Type>& fvm;
 fvm.upper() = [[surfaceInterpolation;;deltaCoeffs|deltaCoeffs]].internalField() * gamma*mesh.magSf().internalField();
 fvm.[[lduMatrix;;negSumDiag|negSumDiag]]();
:'''Note:'''
:*<tt>fvm.lower</tt> is never mentioned - this is because it is equal to <tt>fvm.upper</tt> and a symmetric matrix is being created.
:*<tt>fvm.upper</tt> is opposite in sign to <math>A_r</math> above - this is because <tt>fvm.lower</tt> are [[matrix coefficients]], whereas <math>A_r</math> are [[discretization coefficients]].

=== Examples how the matrix is build ===

==== Laplacian ====
In the following example we show how the coefficients for the matrix resulting from the discretization of the Laplacian equation are derived. Similar information can be found in the very good book of <ref> Moukalled, F., L. Mangani, and M. Darwish. "The finite volume method in computational fluid dynamics." An Advanced Introduction with OpenFOAM and Matlab (2016): </ref> which describes the details of the discretization in OpenFOAM. 

The Laplacian equation reads:   

<math>\boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf \phi = 0</math>

Note the for the explanation here we assume a constant <math>\Gamma</math>. 

Being <math> \phi</math> the variable solved for and <math> \Gamma</math> the diffusivity.

The above equation is integrated over the control volume and together with the Gauss theorem we get

<math>\int \boldsymbol \nabla \bullet \Gamma \boldsymbol \nabla \mathbf \phi dV = \int \Gamma \boldsymbol \nabla \phi \mathbf \bullet \boldsymbol n dS = 0</math>

Substituting the integral over the surface of the volume with a summation over the cell faces we get:

<math> \sum_f  (\Gamma \boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol  dS = 0</math>

The the scalar product of the gradient of <math> \phi </math> with the face Area vector <math> \boldsymbol  dS = 0</math> is discretized as follows:

<math>  (\Gamma \boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol  dS =  \Gamma \frac{\phi_N - \phi_O}{\delta}S_f</math>

The face normal vector in OpenFoam is pointing from the face owner (the cell with the lower index) to the face neighbor (the cell with the higher index). This means that if the value <math> \phi_N </math> is higher than  <math> \phi_O </math> the scalar product <math>  ( \Gamma\boldsymbol \nabla \phi )_f\mathbf \bullet \boldsymbol  dS =  \Gamma \frac{\phi_N - \phi_O}{\delta}S_f</math> is positive and vice versa. 

Finally we get: 

<math> \sum_f \Gamma \frac{\phi_{N,f }- \phi_O}{\delta_f}S_f  = 0</math>


Were <math> \phi_{N,f } , \phi_O, \delta_f, S_f</math> are the variable stored at the face neighbor cell center, the variable stored at the owner cell center, the distance from the owner cell center to the neighbor cell center and the face are, respectively.  Finally we get for the matrix coefficients for the case where the value of <math>\phi</math> at the boundary is specified:


{| class="wikitable" border="1"
|+ Non-zero coefficients
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| <math>- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_3}{\delta_3}+ \frac{S_{14}}{\delta_{14}}   + \frac{S_{25}}{\delta_{25}}  )</math> || <math> \Gamma\frac{S_1}{\delta_1}</math> || || || || <math> \Gamma\frac{S_3}{\delta_3}</math> || || ||
|- align="center"
! -2-
| <math> \Gamma\frac{S_1}{\delta_1}</math> || <math>- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_2}{\delta_2} + \frac{S_4}{\delta_4}+ \frac{S_{15}}{\delta_{15}}    )</math> || <math> \Gamma\frac{S_2}{\delta_2}</math> || || <math> \Gamma\frac{S_4}{\delta_4}</math> || || || ||
|- align="center"
! -3-
| || <math> \Gamma\frac{S_2}{\delta_2}</math> || <math>- \Gamma(\frac{S_2}{\delta_2}+ \frac{S_5}{\delta_5} +\frac{S_{16}}{\delta_{16}}  +\frac{S_{17}}{\delta_{17}}    )</math> || <math> \Gamma\frac{S_5}{\delta_5}</math> || || || || ||
|- align="center"
! -4-
| || || <math> \Gamma\frac{S_5}{\delta_5}</math> || <math>- \Gamma(\frac{S_5}{\delta_5}+ \frac{S_7}{\delta_7} + \frac{S_{10}}{\delta_{10}} + \frac{S_{18}}{\delta_{18}}    )</math> || <math> \Gamma\frac{S_7}{\delta_7}</math> || || || || <math> \Gamma\frac{S_{10}}{\delta_{10}}</math>
|- align="center"
! -5-
| || <math> \Gamma\frac{S_4}{\delta_4}</math> || || <math> \Gamma\frac{S_7}{\delta_7}</math> || <math>- \Gamma(\frac{S_4}{\delta_4}+ \frac{S_6}{\delta_6} + \frac{S_7}{\delta_7} + \frac{S_{9}}{\delta_{9}}  )</math>|| <math> \Gamma\frac{S_6}{\delta_6}</math> || || <math> \Gamma\frac{S_9}{\delta_9}</math> ||
|- align="center"
! -6-
| <math> \Gamma\frac{S_3}{\delta_3}</math> || || || || <math> \Gamma\frac{S_6}{\delta_6}</math> || <math>- \Gamma(\frac{S_3}{\delta_3}+ \frac{S_6}{\delta_6} + \frac{S_{8}}{\delta_{8}} +\frac{S_{24}}{\delta_{24}}   )</math> || <math> \Gamma\frac{S_8}{\delta_8}</math> || ||
|- align="center"
! -7-
| || || || || || <math> \Gamma\frac{S_8}{\delta_8}</math> || <math>- \Gamma( \frac{S_8}{\delta_8} + \frac{S_{11}}{\delta_{11}} +\frac{S_{22}}{\delta_{22}}  +\frac{S_{23}}{\delta_{23}}   )</math>|| <math> \Gamma\frac{S_{11}}{\delta_{11}}</math> ||
|- align="center"
! -8-
| || || || || <math> \Gamma\frac{S_9}{\delta_9}</math> || || <math> \Gamma\frac{S_{11}}{\delta_{11}}</math> || <math>- \Gamma(\frac{S_9}{\delta_9}+ \frac{S_{11}}{\delta_{11}} + \frac{S_{12}}{\delta_{12}} + \frac{S_{21}}{\delta_{21}}   )</math> || <math> \Gamma\frac{S_{12}}{\delta_{12}}</math>
|- align="center"
! -9-
| || || ||<math> \Gamma\frac{S_{10}}{\delta_{10}}</math>|| || || || <math> \Gamma\frac{S_{12}}{\delta_{12}}</math> || <math>- \Gamma(\frac{S_{10}}{\delta_{10}} + \frac{S_{12}}{\delta_{12}}  + \frac{S_{19}}{\delta_{19}}  + \frac{S_{20}}{\delta_{20}}  )</math>
|}

The source vector resulting from the fixed value boundary condition will look like:

{| class="wikitable" border="1"
|+ source vector for the fixed value boundary conditons
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| <math>\Gamma\frac{S_{14}}{\delta_{14}}\phi_{14}   + \Gamma\frac{S_{25}}{\delta_{25}}\phi_{25}</math> ||<math> \Gamma\frac{S_{15}}{\delta_{15}}\phi_{15}</math>  ||<math> \Gamma\frac{S_{16}}{\delta_{16}}\phi_{16}   + \Gamma\frac{S_{17}}{\delta_{17}}\phi_{17}</math> ||<math> \Gamma\frac{S_{18}}{\delta_{18}}\phi_{18}  </math> ||<math> 0</math>|| <math> \Gamma\frac{S_{24}}{\delta_{24}}\phi_{24}</math> ||<math>\Gamma\frac{S_{22}}{\delta_{22}}\phi_{22}   + \Gamma\frac{S_{23}}{\delta_{23}}\phi_{23}</math> ||<math> \Gamma\frac{S_{21}}{\delta_{21}}\phi_{21}</math> || <math>\Gamma\frac{S_{19}}{\delta_{19}}\phi_{19}   + \Gamma\frac{S_{20}}{\delta_{20}}\phi_{20}</math>

|}

For a constant gradient <math>q</math> boundary conditions the matrix looks like:


{| class="wikitable" border="1"
|+ Non-zero coefficients
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| <math>- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_3}{\delta_3}  )</math> || <math> \Gamma\frac{S_1}{\delta_1}</math> || || || || <math> \Gamma\frac{S_3}{\delta_3}</math> || || ||
|- align="center"
! -2-
| <math> \Gamma\frac{S_1}{\delta_1}</math> || <math>- \Gamma(\frac{S_1}{\delta_1}+ \frac{S_2}{\delta_2} + \frac{S_4}{\delta_4}   )</math> || <math> \Gamma\frac{S_2}{\delta_2}</math> || || <math> \Gamma\frac{S_4}{\delta_4}</math> || || || ||
|- align="center"
! -3-
| || <math> \Gamma\frac{S_2}{\delta_2}</math> || <math>- \Gamma(\frac{S_2}{\delta_2}+ \frac{S_5}{\delta_5}   )</math> || <math> \Gamma\frac{S_5}{\delta_5}</math> || || || || ||
|- align="center"
! -4-
| || || <math> \Gamma\frac{S_5}{\delta_5}</math> || <math>- \Gamma(\frac{S_5}{\delta_5}+ \frac{S_7}{\delta_7} + \frac{S_{10}}{\delta_{10}}    )</math> || <math> \Gamma\frac{S_7}{\delta_7}</math> || || || || <math> \Gamma\frac{S_{10}}{\delta_{10}}</math>
|- align="center"
! -5-
| || <math> \Gamma\frac{S_4}{\delta_4}</math> || || <math> \Gamma\frac{S_7}{\delta_7}</math> || <math>- \Gamma(\frac{S_4}{\delta_4}+ \frac{S_6}{\delta_6} + \frac{S_7}{\delta_7} + \frac{S_{9}}{\delta_{9}}  )</math>|| <math> \Gamma\frac{S_6}{\delta_6}</math> || || <math> \Gamma\frac{S_9}{\delta_9}</math> ||
|- align="center"
! -6-
| <math> \Gamma\frac{S_3}{\delta_3}</math> || || || || <math> \Gamma\frac{S_6}{\delta_6}</math> || <math>- \Gamma(\frac{S_3}{\delta_3}+ \frac{S_6}{\delta_6} + \frac{S_{8}}{\delta_{8}}   )</math> || <math> \Gamma\frac{S_8}{\delta_8}</math> || ||
|- align="center"
! -7-
| || || || || || <math> \Gamma\frac{S_8}{\delta_8}</math> || <math>- \Gamma( \frac{S_8}{\delta_8} + \frac{S_{11}}{\delta_{11}}   )</math>|| <math> \Gamma\frac{S_{11}}{\delta_{11}}</math> ||
|- align="center"
! -8-
| || || || || <math> \Gamma\frac{S_9}{\delta_9}</math> || || <math> \Gamma\frac{S_{11}}{\delta_{11}}</math> || <math>- \Gamma(\frac{S_9}{\delta_9}+ \frac{S_{11}}{\delta_{11}} + \frac{S_{12}}{\delta_{12}}   )</math> || <math> \Gamma\frac{S_{12}}{\delta_{12}}</math>
|- align="center"
! -9-
| || || ||<math> \Gamma\frac{S_{10}}{\delta_{10}}</math>|| || || || <math> \Gamma\frac{S_{12}}{\delta_{12}}</math> || <math>- \Gamma(\frac{S_{10}}{\delta_{10}} + \frac{S_{12}}{\delta_{12}}   )</math>
|}

The source vector resulting from the fixed gradient boundary condition will look like:

{| class="wikitable" border="1"
|+ source vector for the fixed value boundary conditons
! i !! -1- !! -2- !! -3- !! -4- !! -5- !! -6- !! -7- !! -8- !! -9-
|- align="center"
! -1-
| <math>\Gamma S_{14}q   + \Gamma S_{25}q</math> ||<math> \Gamma S_{15}q</math>  ||<math> \Gamma S_{16}q   + \Gamma S_{17}q</math> ||<math> \Gamma S_{18}q  </math> ||<math> 0</math>|| <math> \Gamma S_{24} q </math> ||<math>\Gamma S_{22}q   + \Gamma S_{23}q</math> ||<math> \Gamma S_{21}q</math> || <math>\Gamma S_{19}q   + \Gamma S_{20}q</math>

|}

<br><cpp>
template<class Type, class GType>
tmp<fvMatrix<Type>>
gaussLaplacianScheme<Type, GType>::fvmLaplacian
(
    const GeometricField<GType, fvsPatchField, surfaceMesh>& gamma,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    const fvMesh& mesh = this->mesh();

    const surfaceVectorField Sn(mesh.Sf()/mesh.magSf());

    const surfaceVectorField SfGamma(mesh.Sf() & gamma);
    const GeometricField<scalar, fvsPatchField, surfaceMesh> SfGammaSn
    (
        SfGamma & Sn
    );
    const surfaceVectorField SfGammaCorr(SfGamma - SfGammaSn*Sn);

    tmp<fvMatrix<Type>> tfvm = fvmLaplacianUncorrected
    (
        SfGammaSn,
        this->tsnGradScheme_().deltaCoeffs(vf),
        vf
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tfaceFluxCorrection
        = gammaSnGradCorr(SfGammaCorr, vf);

    if (this->tsnGradScheme_().corrected())
    {
        tfaceFluxCorrection.ref() +=
            SfGammaSn*this->tsnGradScheme_().correction(vf);
    }

    fvm.source() -= mesh.V()*fvc::div(tfaceFluxCorrection())().primitiveField();

    if (mesh.fluxRequired(vf.name()))
    {
        fvm.faceFluxCorrectionPtr() = tfaceFluxCorrection.ptr();
    }

    return tfvm;
}



template<class Type, class GType>
tmp<fvMatrix<Type>>
gaussLaplacianScheme<Type, GType>::fvmLaplacianUncorrected
(
    const surfaceScalarField& gammaMagSf,
    const surfaceScalarField& deltaCoeffs,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type>> tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            deltaCoeffs.dimensions()*gammaMagSf.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm.ref();

    fvm.upper() = deltaCoeffs.primitiveField()*gammaMagSf.primitiveField();
    fvm.negSumDiag();

    forAll(vf.boundaryField(), patchi)
    {
        const fvPatchField<Type>& pvf = vf.boundaryField()[patchi];
        const fvsPatchScalarField& pGamma = gammaMagSf.boundaryField()[patchi];
        const fvsPatchScalarField& pDeltaCoeffs =
            deltaCoeffs.boundaryField()[patchi];

        if (pvf.coupled())
        {
            fvm.internalCoeffs()[patchi] =
                pGamma*pvf.gradientInternalCoeffs(pDeltaCoeffs);
            fvm.boundaryCoeffs()[patchi] =
               -pGamma*pvf.gradientBoundaryCoeffs(pDeltaCoeffs);
        }
        else
        {
            fvm.internalCoeffs()[patchi] = pGamma*pvf.gradientInternalCoeffs();
            fvm.boundaryCoeffs()[patchi] = -pGamma*pvf.gradientBoundaryCoeffs();
        }
    }

    return tfvm;
}

</cpp><br>

<br><cpp>

template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedValueFvPatchField<Type>::gradientInternalCoeffs() const
{
    return -pTraits<Type>::one*this->patch().deltaCoeffs();
}


template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedValueFvPatchField<Type>::gradientBoundaryCoeffs() const
{
    return this->patch().deltaCoeffs()*(*this);
}

</cpp><br>

<br><cpp>

template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedGradientFvPatchField<Type>::gradientInternalCoeffs() const
{
    return tmp<Field<Type>>::New(this->size(), Zero);
}


template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::fixedGradientFvPatchField<Type>::gradientBoundaryCoeffs() const
{
    return gradient();
}

</cpp><br>

==== Divergence ====

== fvMatrix ==
Since [[OpenFOAM]] implements multiple numerical strategies such as finite difference, finite area, finite element, and so on, everything that is specific to any one of these strategies is stripped out of the base class lduMatrix.  This facilitates code reuse.

fvMatrix is the finite volume extension of lduMatrix.  It includes the source, '''B''', as well as a reference to the field, <math>\boldsymbol \Psi</math> from <math>\mathbf A \boldsymbol \Psi = \mathbf B</math>.  It handles what happens on the boundaries, and also implements specific operations that occur frequently in finite volume algorithms, such as [[fvMatrix;;H]].

== Notes ==
<references/>
[[Category:OpenFOAM overview]]