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