FreeFOAM The Cross-Platform CFD Toolkit
scalarMatricesTemplates.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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
13  the Free Software Foundation, either version 3 of the License, or
14  (at your 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, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "scalarMatrices.H"
27 #include <OpenFOAM/Swap.H>
28 
29 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
30 
31 template<class Type>
32 void Foam::solve
33 (
34  scalarSquareMatrix& tmpMatrix,
35  Field<Type>& sourceSol
36 )
37 {
38  label n = tmpMatrix.n();
39 
40  // Elimination
41  for (register label i=0; i<n; i++)
42  {
43  label iMax = i;
44  scalar largestCoeff = mag(tmpMatrix[iMax][i]);
45 
46  // Swap entries around to find a good pivot
47  for (register label j=i+1; j<n; j++)
48  {
49  if (mag(tmpMatrix[j][i]) > largestCoeff)
50  {
51  iMax = j;
52  largestCoeff = mag(tmpMatrix[iMax][i]);
53  }
54  }
55 
56  if (i != iMax)
57  {
58  //Info<< "Pivoted on " << i << " " << iMax << endl;
59 
60  for (register label k=i; k<n; k++)
61  {
62  Swap(tmpMatrix[i][k], tmpMatrix[iMax][k]);
63  }
64  Swap(sourceSol[i], sourceSol[iMax]);
65  }
66 
67  // Check that the system of equations isn't singular
68  if (mag(tmpMatrix[i][i]) < 1e-20)
69  {
70  FatalErrorIn("solve(scalarSquareMatrix&, Field<Type>& sourceSol)")
71  << "Singular Matrix"
72  << exit(FatalError);
73  }
74 
75  // Reduce to upper triangular form
76  for (register label j=i+1; j<n; j++)
77  {
78  sourceSol[j] -= sourceSol[i]*(tmpMatrix[j][i]/tmpMatrix[i][i]);
79 
80  for (register label k=n-1; k>=i; k--)
81  {
82  tmpMatrix[j][k] -=
83  tmpMatrix[i][k]*tmpMatrix[j][i]/tmpMatrix[i][i];
84  }
85  }
86  }
87 
88  // Back-substitution
89  for (register label j=n-1; j>=0; j--)
90  {
91  Type ntempvec = pTraits<Type>::zero;
92 
93  for (register label k=j+1; k<n; k++)
94  {
95  ntempvec += tmpMatrix[j][k]*sourceSol[k];
96  }
97 
98  sourceSol[j] = (sourceSol[j] - ntempvec)/tmpMatrix[j][j];
99  }
100 }
101 
102 
103 template<class Type>
104 void Foam::solve
105 (
106  Field<Type>& psi,
107  const scalarSquareMatrix& matrix,
108  const Field<Type>& source
109 )
110 {
111  scalarSquareMatrix tmpMatrix = matrix;
112  psi = source;
113  solve(tmpMatrix, psi);
114 }
115 
116 
117 template<class Type>
119 (
120  const scalarSquareMatrix& luMatrix,
121  const labelList& pivotIndices,
122  Field<Type>& sourceSol
123 )
124 {
125  label n = luMatrix.n();
126 
127  label ii = 0;
128 
129  for (register label i=0; i<n; i++)
130  {
131  label ip = pivotIndices[i];
132  Type sum = sourceSol[ip];
133  sourceSol[ip] = sourceSol[i];
134  const scalar* __restrict__ luMatrixi = luMatrix[i];
135 
136  if (ii != 0)
137  {
138  for (label j=ii-1; j<i; j++)
139  {
140  sum -= luMatrixi[j]*sourceSol[j];
141  }
142  }
143  else if (sum != pTraits<Type>::zero)
144  {
145  ii = i+1;
146  }
147 
148  sourceSol[i] = sum;
149  }
150 
151  for (register label i=n-1; i>=0; i--)
152  {
153  Type sum = sourceSol[i];
154  const scalar* __restrict__ luMatrixi = luMatrix[i];
155 
156  for (register label j=i+1; j<n; j++)
157  {
158  sum -= luMatrixi[j]*sourceSol[j];
159  }
160 
161  sourceSol[i] = sum/luMatrixi[i];
162  }
163 }
164 
165 
166 template<class Type>
167 void Foam::LUsolve
168 (
169  scalarSquareMatrix& matrix,
170  Field<Type>& sourceSol
171 )
172 {
173  labelList pivotIndices(matrix.n());
174  LUDecompose(matrix, pivotIndices);
175  LUBacksubstitute(matrix, pivotIndices, sourceSol);
176 }
177 
178 
179 // ************************ vim: set sw=4 sts=4 et: ************************ //