FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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 <OpenFOAM/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 unallocLabelList& l = lduAddr().lowerAddr();
40  const unallocLabelList& u = lduAddr().upperAddr();
41 
42  for (register 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 unallocLabelList& l = lduAddr().lowerAddr();
57  const unallocLabelList& u = lduAddr().upperAddr();
58 
59  for (register 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 unallocLabelList& l = lduAddr().lowerAddr();
76  const unallocLabelList& u = lduAddr().upperAddr();
77 
78  for (register 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  FatalErrorIn("lduMatrix::operator+=(const lduMatrix& A)")
206  << "Unknown matrix type combination"
207  << abort(FatalError);
208  }
209 }
210 
211 
213 {
214  if (A.diagPtr_)
215  {
216  diag() -= A.diag();
217  }
218 
219  if (symmetric() && A.symmetric())
220  {
221  upper() -= A.upper();
222  }
223  else if (symmetric() && A.asymmetric())
224  {
225  if (upperPtr_)
226  {
227  lower();
228  }
229  else
230  {
231  upper();
232  }
233 
234  upper() -= A.upper();
235  lower() -= A.lower();
236  }
237  else if (asymmetric() && A.symmetric())
238  {
239  if (A.upperPtr_)
240  {
241  lower() -= A.upper();
242  upper() -= A.upper();
243  }
244  else
245  {
246  lower() -= A.lower();
247  upper() -= A.lower();
248  }
249 
250  }
251  else if (asymmetric() && A.asymmetric())
252  {
253  lower() -= A.lower();
254  upper() -= A.upper();
255  }
256  else if (diagonal())
257  {
258  if (A.upperPtr_)
259  {
260  upper() = -A.upper();
261  }
262 
263  if (A.lowerPtr_)
264  {
265  lower() = -A.lower();
266  }
267  }
268  else if (A.diagonal())
269  {
270  }
271  else
272  {
273  FatalErrorIn("lduMatrix::operator-=(const lduMatrix& A)")
274  << "Unknown matrix type combination"
275  << abort(FatalError);
276  }
277 }
278 
279 
281 {
282  if (diagPtr_)
283  {
284  *diagPtr_ *= sf;
285  }
286 
287  if (upperPtr_)
288  {
289  scalarField& upper = *upperPtr_;
290 
291  const unallocLabelList& l = lduAddr().lowerAddr();
292 
293  for (register label face=0; face<upper.size(); face++)
294  {
295  upper[face] *= sf[l[face]];
296  }
297  }
298 
299  if (lowerPtr_)
300  {
301  scalarField& lower = *lowerPtr_;
302 
303  const unallocLabelList& u = lduAddr().upperAddr();
304 
305  for (register label face=0; face<lower.size(); face++)
306  {
307  lower[face] *= sf[u[face]];
308  }
309  }
310 }
311 
312 
314 {
315  if (diagPtr_)
316  {
317  *diagPtr_ *= s;
318  }
319 
320  if (upperPtr_)
321  {
322  *upperPtr_ *= s;
323  }
324 
325  if (lowerPtr_)
326  {
327  *lowerPtr_ *= s;
328  }
329 }
330 
331 
333 {
335  (
336  new scalarField(lduAddr().size(), 0.0)
337  );
338 
339  if (lowerPtr_ || upperPtr_)
340  {
341  scalarField& H1_ = tH1();
342 
343  scalar* __restrict__ H1Ptr = H1_.begin();
344 
345  const label* __restrict__ uPtr = lduAddr().upperAddr().begin();
346  const label* __restrict__ lPtr = lduAddr().lowerAddr().begin();
347 
348  const scalar* __restrict__ lowerPtr = lower().begin();
349  const scalar* __restrict__ upperPtr = upper().begin();
350 
351  register const label nFaces = upper().size();
352 
353  for (register label face=0; face<nFaces; face++)
354  {
355  H1Ptr[uPtr[face]] -= lowerPtr[face];
356  H1Ptr[lPtr[face]] -= upperPtr[face];
357  }
358  }
359 
360  return tH1;
361 }
362 
363 
364 // ************************ vim: set sw=4 sts=4 et: ************************ //