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-2020 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 Description
25  lduMatrix member operations.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "lduMatrix.H"
30 
31 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
32 
34 {
35  if (!lowerPtr_ && !upperPtr_)
36  {
37  return;
38  }
39 
40  const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
41  const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
42  scalarField& Diag = diag();
43 
44  const labelUList& l = lduAddr().lowerAddr();
45  const labelUList& u = lduAddr().upperAddr();
46 
47  for (label face=0; face<l.size(); face++)
48  {
49  Diag[l[face]] += Lower[face];
50  Diag[u[face]] += Upper[face];
51  }
52 }
53 
54 
56 {
57  if (!lowerPtr_ && !upperPtr_)
58  {
59  return;
60  }
61 
62  const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
63  const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
64  scalarField& Diag = diag();
65 
66  const labelUList& l = lduAddr().lowerAddr();
67  const labelUList& u = lduAddr().upperAddr();
68 
69  for (label face=0; face<l.size(); face++)
70  {
71  Diag[l[face]] -= Lower[face];
72  Diag[u[face]] -= Upper[face];
73  }
74 }
75 
76 
78 (
79  scalarField& sumOff
80 ) const
81 {
82  if (!lowerPtr_ && !upperPtr_)
83  {
84  return;
85  }
86 
87  const scalarField& Lower = const_cast<const lduMatrix&>(*this).lower();
88  const scalarField& Upper = const_cast<const lduMatrix&>(*this).upper();
89 
90  const labelUList& l = lduAddr().lowerAddr();
91  const labelUList& u = lduAddr().upperAddr();
92 
93  for (label face = 0; face < l.size(); face++)
94  {
95  sumOff[u[face]] += mag(Lower[face]);
96  sumOff[l[face]] += mag(Upper[face]);
97  }
98 }
99 
100 
101 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
102 
104 {
105  if (this == &A)
106  {
107  FatalError
108  << "lduMatrix::operator=(const lduMatrix&) : "
109  << "attempted assignment to self"
110  << abort(FatalError);
111  }
112 
113  if (A.lowerPtr_)
114  {
115  lower() = A.lower();
116  }
117  else if (lowerPtr_)
118  {
119  delete lowerPtr_;
120  lowerPtr_ = nullptr;
121  }
122 
123  if (A.upperPtr_)
124  {
125  upper() = A.upper();
126  }
127  else if (upperPtr_)
128  {
129  delete upperPtr_;
130  upperPtr_ = nullptr;
131  }
132 
133  if (A.diagPtr_)
134  {
135  diag() = A.diag();
136  }
137 }
138 
139 
141 {
142  if (lowerPtr_)
143  {
144  lowerPtr_->negate();
145  }
146 
147  if (upperPtr_)
148  {
149  upperPtr_->negate();
150  }
151 
152  if (diagPtr_)
153  {
154  diagPtr_->negate();
155  }
156 }
157 
158 
160 {
161  if (A.diagPtr_)
162  {
163  diag() += A.diag();
164  }
165 
166  if (symmetric() && A.symmetric())
167  {
168  upper() += A.upper();
169  }
170  else if (symmetric() && A.asymmetric())
171  {
172  if (upperPtr_)
173  {
174  lower();
175  }
176  else
177  {
178  upper();
179  }
180 
181  upper() += A.upper();
182  lower() += A.lower();
183  }
184  else if (asymmetric() && A.symmetric())
185  {
186  if (A.upperPtr_)
187  {
188  lower() += A.upper();
189  upper() += A.upper();
190  }
191  else
192  {
193  lower() += A.lower();
194  upper() += A.lower();
195  }
196 
197  }
198  else if (asymmetric() && A.asymmetric())
199  {
200  lower() += A.lower();
201  upper() += A.upper();
202  }
203  else if (diagonal())
204  {
205  if (A.upperPtr_)
206  {
207  upper() = A.upper();
208  }
209 
210  if (A.lowerPtr_)
211  {
212  lower() = A.lower();
213  }
214  }
215  else if (A.diagonal())
216  {
217  }
218  else
219  {
220  if (debug > 1)
221  {
223  << "Unknown matrix type combination" << nl
224  << " this :"
225  << " diagonal:" << diagonal()
226  << " symmetric:" << symmetric()
227  << " asymmetric:" << asymmetric() << nl
228  << " A :"
229  << " diagonal:" << A.diagonal()
230  << " symmetric:" << A.symmetric()
231  << " asymmetric:" << A.asymmetric()
232  << endl;
233  }
234  }
235 }
236 
237 
239 {
240  if (A.diagPtr_)
241  {
242  diag() -= A.diag();
243  }
244 
245  if (symmetric() && A.symmetric())
246  {
247  upper() -= A.upper();
248  }
249  else if (symmetric() && A.asymmetric())
250  {
251  if (upperPtr_)
252  {
253  lower();
254  }
255  else
256  {
257  upper();
258  }
259 
260  upper() -= A.upper();
261  lower() -= A.lower();
262  }
263  else if (asymmetric() && A.symmetric())
264  {
265  if (A.upperPtr_)
266  {
267  lower() -= A.upper();
268  upper() -= A.upper();
269  }
270  else
271  {
272  lower() -= A.lower();
273  upper() -= A.lower();
274  }
275 
276  }
277  else if (asymmetric() && A.asymmetric())
278  {
279  lower() -= A.lower();
280  upper() -= A.upper();
281  }
282  else if (diagonal())
283  {
284  if (A.upperPtr_)
285  {
286  upper() = -A.upper();
287  }
288 
289  if (A.lowerPtr_)
290  {
291  lower() = -A.lower();
292  }
293  }
294  else if (A.diagonal())
295  {
296  }
297  else
298  {
299  if (debug > 1)
300  {
302  << "Unknown matrix type combination" << nl
303  << " this :"
304  << " diagonal:" << diagonal()
305  << " symmetric:" << symmetric()
306  << " asymmetric:" << asymmetric() << nl
307  << " A :"
308  << " diagonal:" << A.diagonal()
309  << " symmetric:" << A.symmetric()
310  << " asymmetric:" << A.asymmetric()
311  << endl;
312  }
313  }
314 }
315 
316 
318 {
319  if (diagPtr_)
320  {
321  *diagPtr_ *= sf;
322  }
323 
324  // Non-uniform scaling causes a symmetric matrix
325  // to become asymmetric
326  if (symmetric() || asymmetric())
327  {
328  scalarField& upper = this->upper();
329  scalarField& lower = this->lower();
330 
331  const labelUList& l = lduAddr().lowerAddr();
332  const labelUList& u = lduAddr().upperAddr();
333 
334  for (label face=0; face<upper.size(); face++)
335  {
336  upper[face] *= sf[l[face]];
337  }
338 
339  for (label face=0; face<lower.size(); face++)
340  {
341  lower[face] *= sf[u[face]];
342  }
343  }
344 }
345 
346 
348 {
349  if (diagPtr_)
350  {
351  *diagPtr_ *= s;
352  }
353 
354  if (upperPtr_)
355  {
356  *upperPtr_ *= s;
357  }
358 
359  if (lowerPtr_)
360  {
361  *lowerPtr_ *= s;
362  }
363 }
364 
365 
367 {
368  if (diagPtr_)
369  {
370  *diagPtr_ /= sf;
371  }
372 
373  // Non-uniform scaling causes a symmetric matrix
374  // to become asymmetric
375  if (symmetric() || asymmetric())
376  {
377  scalarField& upper = this->upper();
378  scalarField& lower = this->lower();
379 
380  const labelUList& l = lduAddr().lowerAddr();
381  const labelUList& u = lduAddr().upperAddr();
382 
383  for (label face=0; face<upper.size(); face++)
384  {
385  upper[face] /= sf[l[face]];
386  }
387 
388  for (label face=0; face<lower.size(); face++)
389  {
390  lower[face] /= sf[u[face]];
391  }
392  }
393 }
394 
395 
397 {
398  if (diagPtr_)
399  {
400  *diagPtr_ /= s;
401  }
402 
403  if (upperPtr_)
404  {
405  *upperPtr_ /= s;
406  }
407 
408  if (lowerPtr_)
409  {
410  *lowerPtr_ /= s;
411  }
412 }
413 
414 
415 // ************************************************************************* //
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
void negate()
Negate this field.
Definition: Field.C:446
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
scalarField & upper()
Definition: lduMatrix.C:197
virtual const labelUList & lowerAddr() const =0
Return lower addressing.
void sumMagOffDiag(scalarField &sumOff) const
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))
virtual const labelUList & upperAddr() const =0
Return upper addressing.
errorManip< error > abort(error &err)
Definition: errorManip.H:131
void operator=(const lduMatrix &)
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
bool asymmetric() const
Definition: lduMatrix.H:623
static const char nl
Definition: Ostream.H:260
void operator*=(const scalarField &)
void operator/=(const scalarField &)
volScalarField sf(fieldObject, mesh)
lduMatrix is a general matrix class in which the coefficients are stored as three arrays...
Definition: lduMatrix.H:79
#define WarningInFunction
Report a warning using Foam::Warning.
const lduAddressing & lduAddr() const
Return the LDU addressing.
Definition: lduMatrix.H:550
scalarField & lower()
Definition: lduMatrix.C:168
bool diagonal() const
Definition: lduMatrix.H:595
dimensioned< scalar > mag(const dimensioned< Type > &)
void operator-=(const lduMatrix &)
void operator+=(const lduMatrix &)
scalarField & diag()
Definition: lduMatrix.C:186
label size() const
Return the number of elements in the UList.
Definition: UListI.H:311
bool symmetric() const
Definition: lduMatrix.H:609