LduMatrixOperations.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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"
27 
28 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
29 
30 template<class Type, class DType, class LUType>
32 {
33  if (!lowerPtr_ && !upperPtr_)
34  {
35  return;
36  }
37 
38  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
39  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
40  Field<DType>& Diag = diag();
41 
42  const unallocLabelList& l = lduAddr().lowerAddr();
43  const unallocLabelList& u = lduAddr().upperAddr();
44 
45  for (label face=0; face<l.size(); face++)
46  {
47  Diag[l[face]] += Lower[face];
48  Diag[u[face]] += Upper[face];
49  }
50 }
51 
52 
53 template<class Type, class DType, class LUType>
55 {
56  if (!lowerPtr_ && !upperPtr_)
57  {
58  return;
59  }
60 
61  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
62  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
63  Field<DType>& Diag = diag();
64 
65  const unallocLabelList& l = lduAddr().lowerAddr();
66  const unallocLabelList& u = lduAddr().upperAddr();
67 
68  for (label face=0; face<l.size(); face++)
69  {
70  Diag[l[face]] -= Lower[face];
71  Diag[u[face]] -= Upper[face];
72  }
73 }
74 
75 
76 template<class Type, class DType, class LUType>
78 (
79  Field<LUType>& sumOff
80 ) const
81 {
82  if (!lowerPtr_ && !upperPtr_)
83  {
84  return;
85  }
86 
87  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
88  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
89 
90  const unallocLabelList& l = lduAddr().lowerAddr();
91  const unallocLabelList& u = lduAddr().upperAddr();
92 
93  for (label face = 0; face < l.size(); face++)
94  {
95  sumOff[u[face]] += cmptMag(Lower[face]);
96  sumOff[l[face]] += cmptMag(Upper[face]);
97  }
98 }
99 
100 
101 template<class Type, class DType, class LUType>
104 {
105  tmp<Field<Type>> tHpsi
106  (
107  new Field<Type>(lduAddr().size(), Zero)
108  );
109 
110  if (lowerPtr_ || upperPtr_)
111  {
112  Field<Type> & Hpsi = tHpsi();
113 
114  Type* __restrict__ HpsiPtr = Hpsi.begin();
115 
116  const Type* __restrict__ psiPtr = psi.begin();
117 
118  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
119  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
120 
121  const LUType* __restrict__ lowerPtr = lower().begin();
122  const LUType* __restrict__ upperPtr = upper().begin();
123 
124  const label nFaces = upper().size();
125 
126  for (label face=0; face<nFaces; face++)
127  {
128  HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
129  HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
130  }
131  }
132 
133  return tHpsi;
134 }
135 
136 template<class Type, class DType, class LUType>
139 {
140  tmp<Field<Type>> tHpsi(H(tpsi()));
141  tpsi.clear();
142  return tHpsi;
143 }
144 
145 
146 template<class Type, class DType, class LUType>
149 {
150  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
151  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
152 
153  // Take references to addressing
154  const unallocLabelList& l = lduAddr().lowerAddr();
155  const unallocLabelList& u = lduAddr().upperAddr();
156 
157  tmp<Field<Type>> tfaceHpsi(new Field<Type> (Lower.size()));
158  Field<Type> & faceHpsi = tfaceHpsi();
159 
160  for (label face=0; face<l.size(); face++)
161  {
162  faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
163  }
164 
165  return tfaceHpsi;
166 }
167 
168 
169 template<class Type, class DType, class LUType>
172 {
173  tmp<Field<Type>> tfaceHpsi(faceH(tpsi()));
174  tpsi.clear();
175  return tfaceHpsi;
176 }
177 
178 
179 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
180 
181 template<class Type, class DType, class LUType>
183 {
184  if (this == &A)
185  {
187  << "attempted assignment to self"
188  << abort(FatalError);
189  }
190 
191  if (A.diagPtr_)
192  {
193  diag() = A.diag();
194  }
195 
196  if (A.upperPtr_)
197  {
198  upper() = A.upper();
199  }
200  else if (upperPtr_)
201  {
202  delete upperPtr_;
203  upperPtr_ = nullptr;
204  }
205 
206  if (A.lowerPtr_)
207  {
208  lower() = A.lower();
209  }
210  else if (lowerPtr_)
211  {
212  delete lowerPtr_;
213  lowerPtr_ = nullptr;
214  }
215 
216  if (A.sourcePtr_)
217  {
218  source() = A.source();
219  }
220 
221  interfacesUpper_ = A.interfacesUpper_;
222  interfacesLower_ = A.interfacesLower_;
223 }
224 
225 
226 template<class Type, class DType, class LUType>
228 {
229  if (diagPtr_)
230  {
231  diagPtr_->negate();
232  }
233 
234  if (upperPtr_)
235  {
236  upperPtr_->negate();
237  }
238 
239  if (lowerPtr_)
240  {
241  lowerPtr_->negate();
242  }
243 
244  if (sourcePtr_)
245  {
246  sourcePtr_->negate();
247  }
248 
249  negate(interfacesUpper_);
250  negate(interfacesLower_);
251 }
252 
253 
254 template<class Type, class DType, class LUType>
256 {
257  if (A.diagPtr_)
258  {
259  diag() += A.diag();
260  }
261 
262  if (A.sourcePtr_)
263  {
264  source() += A.source();
265  }
266 
267  if (symmetric() && A.symmetric())
268  {
269  upper() += A.upper();
270  }
271  else if (symmetric() && A.asymmetric())
272  {
273  if (upperPtr_)
274  {
275  lower();
276  }
277  else
278  {
279  upper();
280  }
281 
282  upper() += A.upper();
283  lower() += A.lower();
284  }
285  else if (asymmetric() && A.symmetric())
286  {
287  if (A.upperPtr_)
288  {
289  lower() += A.upper();
290  upper() += A.upper();
291  }
292  else
293  {
294  lower() += A.lower();
295  upper() += A.lower();
296  }
297 
298  }
299  else if (asymmetric() && A.asymmetric())
300  {
301  lower() += A.lower();
302  upper() += A.upper();
303  }
304  else if (diagonal())
305  {
306  if (A.upperPtr_)
307  {
308  upper() = A.upper();
309  }
310 
311  if (A.lowerPtr_)
312  {
313  lower() = A.lower();
314  }
315  }
316  else if (A.diagonal())
317  {
318  }
319  else
320  {
322  << "Unknown matrix type combination"
323  << abort(FatalError);
324  }
325 
326  interfacesUpper_ += A.interfacesUpper_;
327  interfacesLower_ += A.interfacesLower_;
328 }
329 
330 
331 template<class Type, class DType, class LUType>
333 {
334  if (A.diagPtr_)
335  {
336  diag() -= A.diag();
337  }
338 
339  if (A.sourcePtr_)
340  {
341  source() -= A.source();
342  }
343 
344  if (symmetric() && A.symmetric())
345  {
346  upper() -= A.upper();
347  }
348  else if (symmetric() && A.asymmetric())
349  {
350  if (upperPtr_)
351  {
352  lower();
353  }
354  else
355  {
356  upper();
357  }
358 
359  upper() -= A.upper();
360  lower() -= A.lower();
361  }
362  else if (asymmetric() && A.symmetric())
363  {
364  if (A.upperPtr_)
365  {
366  lower() -= A.upper();
367  upper() -= A.upper();
368  }
369  else
370  {
371  lower() -= A.lower();
372  upper() -= A.lower();
373  }
374 
375  }
376  else if (asymmetric() && A.asymmetric())
377  {
378  lower() -= A.lower();
379  upper() -= A.upper();
380  }
381  else if (diagonal())
382  {
383  if (A.upperPtr_)
384  {
385  upper() = -A.upper();
386  }
387 
388  if (A.lowerPtr_)
389  {
390  lower() = -A.lower();
391  }
392  }
393  else if (A.diagonal())
394  {
395  }
396  else
397  {
399  << "Unknown matrix type combination"
400  << abort(FatalError);
401  }
402 
403  interfacesUpper_ -= A.interfacesUpper_;
404  interfacesLower_ -= A.interfacesLower_;
405 }
406 
407 
408 template<class Type, class DType, class LUType>
409 void Foam::LduMatrix<Type, DType, LUType>::operator*=
410 (
411  const scalarField& sf
412 )
413 {
414  if (diagPtr_)
415  {
416  *diagPtr_ *= sf;
417  }
418 
419  if (sourcePtr_)
420  {
421  *sourcePtr_ *= sf;
422  }
423 
424  // Non-uniform scaling causes a symmetric matrix
425  // to become asymmetric
426  if (symmetric() || asymmetric())
427  {
428  Field<LUType>& upper = this->upper();
429  Field<LUType>& lower = this->lower();
430 
431  const unallocLabelList& l = lduAddr().lowerAddr();
432  const unallocLabelList& u = lduAddr().upperAddr();
433 
434  for (label face=0; face<upper.size(); face++)
435  {
436  upper[face] *= sf[l[face]];
437  }
438 
439  for (label face=0; face<lower.size(); face++)
440  {
441  lower[face] *= sf[u[face]];
442  }
443  }
444 
446  << "Scaling a matrix by scalarField is not currently supported\n"
447  "because scaling interfacesUpper_ and interfacesLower_ "
448  "require special transfers"
449  << abort(FatalError);
450 
451  // interfacesUpper_ *= ;
452  // interfacesLower_ *= sf;
453 }
454 
455 
456 template<class Type, class DType, class LUType>
458 {
459  if (diagPtr_)
460  {
461  *diagPtr_ *= s;
462  }
463 
464  if (sourcePtr_)
465  {
466  *sourcePtr_ *= s;
467  }
468 
469  if (upperPtr_)
470  {
471  *upperPtr_ *= s;
472  }
473 
474  if (lowerPtr_)
475  {
476  *lowerPtr_ *= s;
477  }
478 
479  interfacesUpper_ *= s;
480  interfacesLower_ *= s;
481 }
482 
483 
484 // ************************************************************************* //
bool symmetric() const
Definition: LduMatrix.H:590
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
Field< Type > & source()
Definition: LduMatrix.C:248
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
void sumMagOffDiag(Field< LUType > &sumOff) const
Field< LUType > & lower()
Definition: LduMatrix.C:225
void diag(pointPatchField< vector > &, const pointPatchField< tensor > &)
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
tmp< Field< Type > > H(const Field< Type > &) const
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
dimensionSet cmptMag(const dimensionSet &)
Definition: dimensionSet.C:275
const volScalarField & psi
void operator*=(const scalarField &)
Field< LUType > & upper()
Definition: LduMatrix.C:202
static const zero Zero
Definition: zero.H:97
errorManip< error > abort(error &err)
Definition: errorManip.H:131
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 operator-=(const LduMatrix< Type, DType, LUType > &)
volScalarField sf(fieldObject, mesh)
bool asymmetric() const
Definition: LduMatrix.H:604
bool diagonal() const
Definition: LduMatrix.H:576
Field< DType > & diag()
Definition: LduMatrix.C:190
void negate(FieldField< Field, Type > &res, const FieldField< Field, Type > &f)
LduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: LduMatrix.H:69
void operator+=(const LduMatrix< Type, DType, LUType > &)
tmp< Field< Type > > faceH(const Field< Type > &) const
A class for managing temporary objects.
Definition: PtrList.H:53
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
void operator=(const LduMatrix< Type, DType, LUType > &)