LduMatrixOperations.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-2016 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  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
34  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
35  Field<DType>& Diag = diag();
36 
37  const unallocLabelList& l = lduAddr().lowerAddr();
38  const unallocLabelList& u = lduAddr().upperAddr();
39 
40  for (label face=0; face<l.size(); face++)
41  {
42  Diag[l[face]] += Lower[face];
43  Diag[u[face]] += Upper[face];
44  }
45 }
46 
47 
48 template<class Type, class DType, class LUType>
50 {
51  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
52  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
53  Field<DType>& Diag = diag();
54 
55  const unallocLabelList& l = lduAddr().lowerAddr();
56  const unallocLabelList& u = lduAddr().upperAddr();
57 
58  for (label face=0; face<l.size(); face++)
59  {
60  Diag[l[face]] -= Lower[face];
61  Diag[u[face]] -= Upper[face];
62  }
63 }
64 
65 
66 template<class Type, class DType, class LUType>
68 (
69  Field<LUType>& sumOff
70 ) const
71 {
72  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
73  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
74 
75  const unallocLabelList& l = lduAddr().lowerAddr();
76  const unallocLabelList& u = lduAddr().upperAddr();
77 
78  for (label face = 0; face < l.size(); face++)
79  {
80  sumOff[u[face]] += cmptMag(Lower[face]);
81  sumOff[l[face]] += cmptMag(Upper[face]);
82  }
83 }
84 
85 
86 template<class Type, class DType, class LUType>
89 {
90  tmp<Field<Type>> tHpsi
91  (
92  new Field<Type>(lduAddr().size(), Zero)
93  );
94 
95  if (lowerPtr_ || upperPtr_)
96  {
97  Field<Type> & Hpsi = tHpsi();
98 
99  Type* __restrict__ HpsiPtr = Hpsi.begin();
100 
101  const Type* __restrict__ psiPtr = psi.begin();
102 
103  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
104  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
105 
106  const LUType* __restrict__ lowerPtr = lower().begin();
107  const LUType* __restrict__ upperPtr = upper().begin();
108 
109  const label nFaces = upper().size();
110 
111  for (label face=0; face<nFaces; face++)
112  {
113  HpsiPtr[uPtr[face]] -= lowerPtr[face]*psiPtr[lPtr[face]];
114  HpsiPtr[lPtr[face]] -= upperPtr[face]*psiPtr[uPtr[face]];
115  }
116  }
117 
118  return tHpsi;
119 }
120 
121 template<class Type, class DType, class LUType>
124 {
125  tmp<Field<Type>> tHpsi(H(tpsi()));
126  tpsi.clear();
127  return tHpsi;
128 }
129 
130 
131 template<class Type, class DType, class LUType>
134 {
135  const Field<LUType>& Lower = const_cast<const LduMatrix&>(*this).lower();
136  const Field<LUType>& Upper = const_cast<const LduMatrix&>(*this).upper();
137 
138  // Take refereces to addressing
139  const unallocLabelList& l = lduAddr().lowerAddr();
140  const unallocLabelList& u = lduAddr().upperAddr();
141 
142  tmp<Field<Type>> tfaceHpsi(new Field<Type> (Lower.size()));
143  Field<Type> & faceHpsi = tfaceHpsi();
144 
145  for (label face=0; face<l.size(); face++)
146  {
147  faceHpsi[face] = Upper[face]*psi[u[face]] - Lower[face]*psi[l[face]];
148  }
149 
150  return tfaceHpsi;
151 }
152 
153 
154 template<class Type, class DType, class LUType>
157 {
158  tmp<Field<Type>> tfaceHpsi(faceH(tpsi()));
159  tpsi.clear();
160  return tfaceHpsi;
161 }
162 
163 
164 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
165 
166 template<class Type, class DType, class LUType>
168 {
169  if (this == &A)
170  {
172  << "attempted assignment to self"
173  << abort(FatalError);
174  }
175 
176  if (A.diagPtr_)
177  {
178  diag() = A.diag();
179  }
180 
181  if (A.upperPtr_)
182  {
183  upper() = A.upper();
184  }
185  else if (upperPtr_)
186  {
187  delete upperPtr_;
188  upperPtr_ = nullptr;
189  }
190 
191  if (A.lowerPtr_)
192  {
193  lower() = A.lower();
194  }
195  else if (lowerPtr_)
196  {
197  delete lowerPtr_;
198  lowerPtr_ = nullptr;
199  }
200 
201  if (A.sourcePtr_)
202  {
203  source() = A.source();
204  }
205 
206  interfacesUpper_ = A.interfacesUpper_;
207  interfacesLower_ = A.interfacesLower_;
208 }
209 
210 
211 template<class Type, class DType, class LUType>
213 {
214  if (diagPtr_)
215  {
216  diagPtr_->negate();
217  }
218 
219  if (upperPtr_)
220  {
221  upperPtr_->negate();
222  }
223 
224  if (lowerPtr_)
225  {
226  lowerPtr_->negate();
227  }
228 
229  if (sourcePtr_)
230  {
231  sourcePtr_->negate();
232  }
233 
234  negate(interfacesUpper_);
235  negate(interfacesLower_);
236 }
237 
238 
239 template<class Type, class DType, class LUType>
241 {
242  if (A.diagPtr_)
243  {
244  diag() += A.diag();
245  }
246 
247  if (A.sourcePtr_)
248  {
249  source() += A.source();
250  }
251 
252  if (symmetric() && A.symmetric())
253  {
254  upper() += A.upper();
255  }
256  else if (symmetric() && A.asymmetric())
257  {
258  if (upperPtr_)
259  {
260  lower();
261  }
262  else
263  {
264  upper();
265  }
266 
267  upper() += A.upper();
268  lower() += A.lower();
269  }
270  else if (asymmetric() && A.symmetric())
271  {
272  if (A.upperPtr_)
273  {
274  lower() += A.upper();
275  upper() += A.upper();
276  }
277  else
278  {
279  lower() += A.lower();
280  upper() += A.lower();
281  }
282 
283  }
284  else if (asymmetric() && A.asymmetric())
285  {
286  lower() += A.lower();
287  upper() += A.upper();
288  }
289  else if (diagonal())
290  {
291  if (A.upperPtr_)
292  {
293  upper() = A.upper();
294  }
295 
296  if (A.lowerPtr_)
297  {
298  lower() = A.lower();
299  }
300  }
301  else if (A.diagonal())
302  {
303  }
304  else
305  {
307  << "Unknown matrix type combination"
308  << abort(FatalError);
309  }
310 
311  interfacesUpper_ += A.interfacesUpper_;
312  interfacesLower_ += A.interfacesLower_;
313 }
314 
315 
316 template<class Type, class DType, class LUType>
318 {
319  if (A.diagPtr_)
320  {
321  diag() -= A.diag();
322  }
323 
324  if (A.sourcePtr_)
325  {
326  source() -= A.source();
327  }
328 
329  if (symmetric() && A.symmetric())
330  {
331  upper() -= A.upper();
332  }
333  else if (symmetric() && A.asymmetric())
334  {
335  if (upperPtr_)
336  {
337  lower();
338  }
339  else
340  {
341  upper();
342  }
343 
344  upper() -= A.upper();
345  lower() -= A.lower();
346  }
347  else if (asymmetric() && A.symmetric())
348  {
349  if (A.upperPtr_)
350  {
351  lower() -= A.upper();
352  upper() -= A.upper();
353  }
354  else
355  {
356  lower() -= A.lower();
357  upper() -= A.lower();
358  }
359 
360  }
361  else if (asymmetric() && A.asymmetric())
362  {
363  lower() -= A.lower();
364  upper() -= A.upper();
365  }
366  else if (diagonal())
367  {
368  if (A.upperPtr_)
369  {
370  upper() = -A.upper();
371  }
372 
373  if (A.lowerPtr_)
374  {
375  lower() = -A.lower();
376  }
377  }
378  else if (A.diagonal())
379  {
380  }
381  else
382  {
384  << "Unknown matrix type combination"
385  << abort(FatalError);
386  }
387 
388  interfacesUpper_ -= A.interfacesUpper_;
389  interfacesLower_ -= A.interfacesLower_;
390 }
391 
392 
393 template<class Type, class DType, class LUType>
394 void Foam::LduMatrix<Type, DType, LUType>::operator*=
395 (
396  const scalarField& sf
397 )
398 {
399  if (diagPtr_)
400  {
401  *diagPtr_ *= sf;
402  }
403 
404  if (sourcePtr_)
405  {
406  *sourcePtr_ *= sf;
407  }
408 
409  // Non-uniform scaling causes a symmetric matrix
410  // to become asymmetric
411  if (symmetric() || asymmetric())
412  {
413  Field<LUType>& upper = this->upper();
414  Field<LUType>& lower = this->lower();
415 
416  const unallocLabelList& l = lduAddr().lowerAddr();
417  const unallocLabelList& u = lduAddr().upperAddr();
418 
419  for (label face=0; face<upper.size(); face++)
420  {
421  upper[face] *= sf[l[face]];
422  }
423 
424  for (label face=0; face<lower.size(); face++)
425  {
426  lower[face] *= sf[u[face]];
427  }
428  }
429 
431  << "Scaling a matrix by scalarField is not currently supported\n"
432  "because scaling interfacesUpper_ and interfacesLower_ "
433  "require special transfers"
434  << abort(FatalError);
435 
436  //interfacesUpper_ *= ;
437  //interfacesLower_ *= sf;
438 }
439 
440 
441 template<class Type, class DType, class LUType>
443 {
444  if (diagPtr_)
445  {
446  *diagPtr_ *= s;
447  }
448 
449  if (sourcePtr_)
450  {
451  *sourcePtr_ *= s;
452  }
453 
454  if (upperPtr_)
455  {
456  *upperPtr_ *= s;
457  }
458 
459  if (lowerPtr_)
460  {
461  *lowerPtr_ *= s;
462  }
463 
464  interfacesUpper_ *= s;
465  interfacesLower_ *= s;
466 }
467 
468 
469 // ************************************************************************* //
bool symmetric() const
Definition: LduMatrix.H:581
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
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:319
Field< Type > & source()
Definition: LduMatrix.C:248
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
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
void operator*=(const scalarField &)
Field< LUType > & upper()
Definition: LduMatrix.C:202
static const zero Zero
Definition: zero.H:91
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:61
void operator-=(const LduMatrix< Type, DType, LUType > &)
void cmptMag(FieldField< Field, Type > &cf, const FieldField< Field, Type > &f)
volScalarField sf(fieldObject, mesh)
bool asymmetric() const
Definition: LduMatrix.H:586
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
const volScalarField & psi
A class for managing temporary objects.
Definition: PtrList.H:53
label size() const
Return the number of elements in the UList.
Definition: UListI.H:299
void operator=(const LduMatrix< Type, DType, LUType > &)