LduMatrixATmul.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) 2011-2015 OpenFOAM Foundation
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 "LduMatrix.H"
28 
29 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33 
34 template<class Type, class LUType>
36 :
38 {
39  const Field<LUType>& A_;
40 
41 public:
42 
44  :
45  A_(A)
46  {}
47 
48  virtual ~Amultiplier()
49  {}
50 
51  virtual void addAmul(Field<Type>& Apsi, const Field<Type>& psi) const
52  {
53  Apsi += A_*psi;
54  }
55 };
56 
57 }
58 
59 
60 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
61 
62 template<class Type, class DType, class LUType>
64 (
65  Field<Type>& Apsi,
66  const tmp<Field<Type> >& tpsi
67 ) const
68 {
69  Type* __restrict__ ApsiPtr = Apsi.begin();
70 
71  const Field<Type>& psi = tpsi();
72  const Type* const __restrict__ psiPtr = psi.begin();
73 
74  const DType* const __restrict__ diagPtr = diag().begin();
75 
76  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
77  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
78 
79  const LUType* const __restrict__ upperPtr = upper().begin();
80  const LUType* const __restrict__ lowerPtr = lower().begin();
81 
82  // Initialise the update of interfaced interfaces
83  initMatrixInterfaces
84  (
85  interfacesUpper_,
86  psi,
87  Apsi
88  );
89 
90  const label nCells = diag().size();
91  for (label cell=0; cell<nCells; cell++)
92  {
93  ApsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
94  }
95 
96 
97  const label nFaces = upper().size();
98  for (label face=0; face<nFaces; face++)
99  {
100  ApsiPtr[uPtr[face]] += dot(lowerPtr[face], psiPtr[lPtr[face]]);
101  ApsiPtr[lPtr[face]] += dot(upperPtr[face], psiPtr[uPtr[face]]);
102  }
103 
104  // Update interface interfaces
105  updateMatrixInterfaces
106  (
107  interfacesUpper_,
108  psi,
109  Apsi
110  );
111 
112  tpsi.clear();
113 }
114 
115 
116 template<class Type, class DType, class LUType>
118 (
119  Field<Type>& Tpsi,
120  const tmp<Field<Type> >& tpsi
121 ) const
122 {
123  Type* __restrict__ TpsiPtr = Tpsi.begin();
124 
125  const Field<Type>& psi = tpsi();
126  const Type* const __restrict__ psiPtr = psi.begin();
127 
128  const DType* const __restrict__ diagPtr = diag().begin();
129 
130  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
131  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
132 
133  const LUType* const __restrict__ lowerPtr = lower().begin();
134  const LUType* const __restrict__ upperPtr = upper().begin();
135 
136  // Initialise the update of interfaced interfaces
137  initMatrixInterfaces
138  (
139  interfacesLower_,
140  psi,
141  Tpsi
142  );
143 
144  const label nCells = diag().size();
145  for (label cell=0; cell<nCells; cell++)
146  {
147  TpsiPtr[cell] = dot(diagPtr[cell], psiPtr[cell]);
148  }
149 
150  const label nFaces = upper().size();
151  for (label face=0; face<nFaces; face++)
152  {
153  TpsiPtr[uPtr[face]] += dot(upperPtr[face], psiPtr[lPtr[face]]);
154  TpsiPtr[lPtr[face]] += dot(lowerPtr[face], psiPtr[uPtr[face]]);
155  }
156 
157  // Update interface interfaces
158  updateMatrixInterfaces
159  (
160  interfacesLower_,
161  psi,
162  Tpsi
163  );
164 
165  tpsi.clear();
166 }
167 
168 
169 template<class Type, class DType, class LUType>
171 (
172  Field<Type>& sumA
173 ) const
174 {
175  Type* __restrict__ sumAPtr = sumA.begin();
176 
177  const DType* __restrict__ diagPtr = diag().begin();
178 
179  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
180  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
181 
182  const LUType* __restrict__ lowerPtr = lower().begin();
183  const LUType* __restrict__ upperPtr = upper().begin();
184 
185  const label nCells = diag().size();
186  const label nFaces = upper().size();
187 
188  for (label cell=0; cell<nCells; cell++)
189  {
190  sumAPtr[cell] = dot(diagPtr[cell], pTraits<Type>::one);
191  }
192 
193  for (label face=0; face<nFaces; face++)
194  {
195  sumAPtr[uPtr[face]] += dot(lowerPtr[face], pTraits<Type>::one);
196  sumAPtr[lPtr[face]] += dot(upperPtr[face], pTraits<Type>::one);
197  }
198 
199  // Add the interface internal coefficients to diagonal
200  // and the interface boundary coefficients to the sum-off-diagonal
201  forAll(interfaces_, patchI)
202  {
203  if (interfaces_.set(patchI))
204  {
205  const unallocLabelList& pa = lduAddr().patchAddr(patchI);
206  const Field<LUType>& pCoeffs = interfacesUpper_[patchI];
207 
208  forAll(pa, face)
209  {
210  sumAPtr[pa[face]] -= dot(pCoeffs[face], pTraits<Type>::one);
211  }
212  }
213  }
214 }
215 
216 
217 template<class Type, class DType, class LUType>
219 (
220  Field<Type>& rA,
221  const Field<Type>& psi
222 ) const
223 {
224  Type* __restrict__ rAPtr = rA.begin();
225 
226  const Type* const __restrict__ psiPtr = psi.begin();
227  const DType* const __restrict__ diagPtr = diag().begin();
228  const Type* const __restrict__ sourcePtr = source().begin();
229 
230  const label* const __restrict__ uPtr = lduAddr().upperAddr().begin();
231  const label* const __restrict__ lPtr = lduAddr().lowerAddr().begin();
232 
233  const LUType* const __restrict__ upperPtr = upper().begin();
234  const LUType* const __restrict__ lowerPtr = lower().begin();
235 
236  // Parallel boundary initialisation.
237  // Note: there is a change of sign in the coupled
238  // interface update to add the contibution to the r.h.s.
239 
240  FieldField<Field, LUType> mBouCoeffs(interfacesUpper_.size());
241 
242  forAll(mBouCoeffs, patchi)
243  {
244  if (interfaces_.set(patchi))
245  {
246  mBouCoeffs.set(patchi, -interfacesUpper_[patchi]);
247  }
248  }
249 
250  // Initialise the update of interfaced interfaces
251  initMatrixInterfaces
252  (
253  mBouCoeffs,
254  psi,
255  rA
256  );
257 
258  const label nCells = diag().size();
259  for (label cell=0; cell<nCells; cell++)
260  {
261  rAPtr[cell] = sourcePtr[cell] - dot(diagPtr[cell], psiPtr[cell]);
262  }
263 
264 
265  const label nFaces = upper().size();
266  for (label face=0; face<nFaces; face++)
267  {
268  rAPtr[uPtr[face]] -= dot(lowerPtr[face], psiPtr[lPtr[face]]);
269  rAPtr[lPtr[face]] -= dot(upperPtr[face], psiPtr[uPtr[face]]);
270  }
271 
272  // Update interface interfaces
273  updateMatrixInterfaces
274  (
275  mBouCoeffs,
276  psi,
277  rA
278  );
279 }
280 
281 
282 template<class Type, class DType, class LUType>
284 (
285  const Field<Type>& psi
286 ) const
287 {
288  tmp<Field<Type> > trA(new Field<Type>(psi.size()));
289  residual(trA(), psi);
290  return trA;
291 }
292 
293 
294 // ************************************************************************* //
Generic field type.
Definition: FieldField.H:51
An abstract base class for implicitly-coupled interface fields e.g. processor and cyclic patch fields...
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
virtual ~Amultiplier()
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
void Tmul(Field< Type > &, const tmp< Field< Type > > &) const
Matrix transpose multiplication.
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
Namespace for OpenFOAM.
Amultiplier(const Field< LUType > &A)
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
iterator begin()
Return an iterator to begin traversing the UList.
void dot(FieldField< Field1, typename innerProduct< Type1, Type2 >::type > &f, const FieldField< Field1, Type1 > &f1, const FieldField< Field2, Type2 > &f2)
void sumA(Field< Type > &) const
Sum the coefficients on each row of the matrix.
#define forAll(list, i)
Definition: UList.H:421
label patchi
void Amul(Field< Type > &, const tmp< Field< Type > > &) const
Matrix multiplication.
const volScalarField & psi
Definition: createFields.H:24
Traits class for primitives.
Definition: pTraits.H:50
virtual void addAmul(Field< Type > &Apsi, const Field< Type > &psi) const
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:60
void residual(Field< Type > &rA, const Field< Type > &psi) const
A class for managing temporary objects.
Definition: PtrList.H:118
List of coupled interface fields to be used in coupling.