FreeFOAM The Cross-Platform CFD Toolkit
MRFZone.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 \*---------------------------------------------------------------------------*/
25 
26 #include "MRFZone.H"
27 #include <finiteVolume/fvMesh.H>
28 #include <finiteVolume/volFields.H>
31 #include <OpenFOAM/syncTools.H>
32 #include <meshTools/faceSet.H>
34 
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
36 
38 
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::MRFZone::setMRFFaces()
43 {
44  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
45 
46  // Type per face:
47  // 0:not in zone
48  // 1:moving with frame
49  // 2:other
50  labelList faceType(mesh_.nFaces(), 0);
51 
52  // Determine faces in cell zone
53  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54  // (without constructing cells)
55 
56  const labelList& own = mesh_.faceOwner();
57  const labelList& nei = mesh_.faceNeighbour();
58 
59  // Cells in zone
60  boolList zoneCell(mesh_.nCells(), false);
61 
62  if (cellZoneID_ != -1)
63  {
64  const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
65  forAll(cellLabels, i)
66  {
67  zoneCell[cellLabels[i]] = true;
68  }
69  }
70 
71 
72  label nZoneFaces = 0;
73 
74  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
75  {
76  if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
77  {
78  faceType[faceI] = 1;
79  nZoneFaces++;
80  }
81  }
82 
83 
84  labelHashSet excludedPatches(excludedPatchLabels_);
85 
86  forAll(patches, patchI)
87  {
88  const polyPatch& pp = patches[patchI];
89 
90  if (pp.coupled() || excludedPatches.found(patchI))
91  {
92  forAll(pp, i)
93  {
94  label faceI = pp.start()+i;
95 
96  if (zoneCell[own[faceI]])
97  {
98  faceType[faceI] = 2;
99  nZoneFaces++;
100  }
101  }
102  }
103  else if (!isA<emptyPolyPatch>(pp))
104  {
105  forAll(pp, i)
106  {
107  label faceI = pp.start()+i;
108 
109  if (zoneCell[own[faceI]])
110  {
111  faceType[faceI] = 1;
112  nZoneFaces++;
113  }
114  }
115  }
116  }
117 
118  // Now we have for faceType:
119  // 0 : face not in cellZone
120  // 1 : internal face or normal patch face
121  // 2 : coupled patch face or excluded patch face
122 
123  // Sort into lists per patch.
124 
125  internalFaces_.setSize(mesh_.nFaces());
126  label nInternal = 0;
127 
128  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
129  {
130  if (faceType[faceI] == 1)
131  {
132  internalFaces_[nInternal++] = faceI;
133  }
134  }
135  internalFaces_.setSize(nInternal);
136 
137  labelList nIncludedFaces(patches.size(), 0);
138  labelList nExcludedFaces(patches.size(), 0);
139 
140  forAll(patches, patchi)
141  {
142  const polyPatch& pp = patches[patchi];
143 
144  forAll(pp, patchFacei)
145  {
146  label faceI = pp.start()+patchFacei;
147 
148  if (faceType[faceI] == 1)
149  {
150  nIncludedFaces[patchi]++;
151  }
152  else if (faceType[faceI] == 2)
153  {
154  nExcludedFaces[patchi]++;
155  }
156  }
157  }
158 
159  includedFaces_.setSize(patches.size());
160  excludedFaces_.setSize(patches.size());
161  forAll(nIncludedFaces, patchi)
162  {
163  includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
164  excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
165  }
166  nIncludedFaces = 0;
167  nExcludedFaces = 0;
168 
169  forAll(patches, patchi)
170  {
171  const polyPatch& pp = patches[patchi];
172 
173  forAll(pp, patchFacei)
174  {
175  label faceI = pp.start()+patchFacei;
176 
177  if (faceType[faceI] == 1)
178  {
179  includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
180  }
181  else if (faceType[faceI] == 2)
182  {
183  excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
184  }
185  }
186  }
187 
188 
189  if (debug)
190  {
191  faceSet internalFaces(mesh_, "internalFaces", internalFaces_);
192  Pout<< "Writing " << internalFaces.size()
193  << " internal faces in MRF zone to faceSet "
194  << internalFaces.name() << endl;
195  internalFaces.write();
196 
197  faceSet MRFFaces(mesh_, "includedFaces", 100);
198  forAll(includedFaces_, patchi)
199  {
200  forAll(includedFaces_[patchi], i)
201  {
202  label patchFacei = includedFaces_[patchi][i];
203  MRFFaces.insert(patches[patchi].start()+patchFacei);
204  }
205  }
206  Pout<< "Writing " << MRFFaces.size()
207  << " patch faces in MRF zone to faceSet "
208  << MRFFaces.name() << endl;
209  MRFFaces.write();
210 
211  faceSet excludedFaces(mesh_, "excludedFaces", 100);
212  forAll(excludedFaces_, patchi)
213  {
214  forAll(excludedFaces_[patchi], i)
215  {
216  label patchFacei = excludedFaces_[patchi][i];
217  excludedFaces.insert(patches[patchi].start()+patchFacei);
218  }
219  }
220  Pout<< "Writing " << excludedFaces.size()
221  << " faces in MRF zone with special handling to faceSet "
222  << excludedFaces.name() << endl;
223  excludedFaces.write();
224  }
225 }
226 
227 
228 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
229 
230 Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
231 :
232  mesh_(mesh),
233  name_(is),
234  dict_(is),
235  cellZoneID_(mesh_.cellZones().findZoneID(name_)),
236  excludedPatchNames_
237  (
238  dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
239  ),
240  origin_(dict_.lookup("origin")),
241  axis_(dict_.lookup("axis")),
242  omega_(dict_.lookup("omega")),
243  Omega_("Omega", omega_*axis_)
244 {
245  if (dict_.found("patches"))
246  {
247  WarningIn("MRFZone(const fvMesh&, Istream&)")
248  << "Ignoring entry 'patches'\n"
249  << " By default all patches within the rotating region rotate.\n"
250  << " Optionally supply excluded patches using 'nonRotatingPatches'."
251  << endl;
252  }
253 
254  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
255 
256  axis_ = axis_/mag(axis_);
257  Omega_ = omega_*axis_;
258 
259  excludedPatchLabels_.setSize(excludedPatchNames_.size());
260 
261  forAll(excludedPatchNames_, i)
262  {
263  excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
264 
265  if (excludedPatchLabels_[i] == -1)
266  {
268  (
269  "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
270  ) << "cannot find MRF patch " << excludedPatchNames_[i]
271  << exit(FatalError);
272  }
273  }
274 
275 
276  bool cellZoneFound = (cellZoneID_ != -1);
277  reduce(cellZoneFound, orOp<bool>());
278 
279  if (!cellZoneFound)
280  {
282  (
283  "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
284  ) << "cannot find MRF cellZone " << name_
285  << exit(FatalError);
286  }
287 
288  setMRFFaces();
289 }
290 
291 
292 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
293 
295 {
296  if (cellZoneID_ == -1)
297  {
298  return;
299  }
300 
301  const labelList& cells = mesh_.cellZones()[cellZoneID_];
302  const scalarField& V = mesh_.V();
303  vectorField& Usource = UEqn.source();
304  const vectorField& U = UEqn.psi();
305  const vector& Omega = Omega_.value();
306 
307  forAll(cells, i)
308  {
309  label celli = cells[i];
310  Usource[celli] -= V[celli]*(Omega ^ U[celli]);
311  }
312 }
313 
314 
316 (
317  const volScalarField& rho,
319 ) const
320 {
321  if (cellZoneID_ == -1)
322  {
323  return;
324  }
325 
326  const labelList& cells = mesh_.cellZones()[cellZoneID_];
327  const scalarField& V = mesh_.V();
328  vectorField& Usource = UEqn.source();
329  const vectorField& U = UEqn.psi();
330  const vector& Omega = Omega_.value();
331 
332  forAll(cells, i)
333  {
334  label celli = cells[i];
335  Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
336  }
337 }
338 
339 
341 {
342  const volVectorField& C = mesh_.C();
343 
344  const vector& origin = origin_.value();
345  const vector& Omega = Omega_.value();
346 
347  const labelList& cells = mesh_.cellZones()[cellZoneID_];
348 
349  forAll(cells, i)
350  {
351  label celli = cells[i];
352  U[celli] -= (Omega ^ (C[celli] - origin));
353  }
354 
355  // Included patches
356  forAll(includedFaces_, patchi)
357  {
358  forAll(includedFaces_[patchi], i)
359  {
360  label patchFacei = includedFaces_[patchi][i];
361  U.boundaryField()[patchi][patchFacei] = vector::zero;
362  }
363  }
364 
365  // Excluded patches
366  forAll(excludedFaces_, patchi)
367  {
368  forAll(excludedFaces_[patchi], i)
369  {
370  label patchFacei = excludedFaces_[patchi][i];
371  U.boundaryField()[patchi][patchFacei] -=
372  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
373  }
374  }
375 }
376 
377 
379 {
380  const volVectorField& C = mesh_.C();
381 
382  const vector& origin = origin_.value();
383  const vector& Omega = Omega_.value();
384 
385  const labelList& cells = mesh_.cellZones()[cellZoneID_];
386 
387  forAll(cells, i)
388  {
389  label celli = cells[i];
390  U[celli] += (Omega ^ (C[celli] - origin));
391  }
392 
393  // Included patches
394  forAll(includedFaces_, patchi)
395  {
396  forAll(includedFaces_[patchi], i)
397  {
398  label patchFacei = includedFaces_[patchi][i];
399  U.boundaryField()[patchi][patchFacei] =
400  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
401  }
402  }
403 
404  // Excluded patches
405  forAll(excludedFaces_, patchi)
406  {
407  forAll(excludedFaces_[patchi], i)
408  {
409  label patchFacei = excludedFaces_[patchi][i];
410  U.boundaryField()[patchi][patchFacei] +=
411  (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
412  }
413  }
414 }
415 
416 
418 {
419  relativeRhoFlux(geometricOneField(), phi);
420 }
421 
422 
424 (
425  const surfaceScalarField& rho,
427 ) const
428 {
429  relativeRhoFlux(rho, phi);
430 }
431 
432 
434 {
435  absoluteRhoFlux(geometricOneField(), phi);
436 }
437 
438 
440 (
441  const surfaceScalarField& rho,
442  surfaceScalarField& phi
443 ) const
444 {
445  absoluteRhoFlux(rho, phi);
446 }
447 
448 
450 {
451  const vector& origin = origin_.value();
452  const vector& Omega = Omega_.value();
453 
454  // Included patches
455  forAll(includedFaces_, patchi)
456  {
457  const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
458 
459  vectorField pfld(U.boundaryField()[patchi]);
460 
461  forAll(includedFaces_[patchi], i)
462  {
463  label patchFacei = includedFaces_[patchi][i];
464 
465  pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
466  }
467 
468  U.boundaryField()[patchi] == pfld;
469  }
470 }
471 
472 
473 // ************************ vim: set sw=4 sts=4 et: ************************ //