1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
28 #include <finiteVolume/fvMesh.H>
29 #include <finiteVolume/volFields.H>
30 #include <finiteVolume/surfaceFields.H>
31 #include <finiteVolume/fvMatrices.H>
32 #include <OpenFOAM/syncTools.H>
33 #include <meshTools/faceSet.H>
35 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 defineTypeNameAndDebug(Foam::MRFZone, 0);
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 void Foam::MRFZone::setMRFFaces()
44 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
48 // 1:moving with frame
50 labelList faceType(mesh_.nFaces(), 0);
52 // Determine faces in cell zone
53 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
54 // (without constructing cells)
56 const labelList& own = mesh_.faceOwner();
57 const labelList& nei = mesh_.faceNeighbour();
60 boolList zoneCell(mesh_.nCells(), false);
62 if (cellZoneID_ != -1)
64 const labelList& cellLabels = mesh_.cellZones()[cellZoneID_];
67 zoneCell[cellLabels[i]] = true;
74 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
76 if (zoneCell[own[faceI]] || zoneCell[nei[faceI]])
84 labelHashSet excludedPatches(excludedPatchLabels_);
86 forAll(patches, patchI)
88 const polyPatch& pp = patches[patchI];
90 if (pp.coupled() || excludedPatches.found(patchI))
94 label faceI = pp.start()+i;
96 if (zoneCell[own[faceI]])
103 else if (!isA<emptyPolyPatch>(pp))
107 label faceI = pp.start()+i;
109 if (zoneCell[own[faceI]])
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
123 // Sort into lists per patch.
125 internalFaces_.setSize(mesh_.nFaces());
128 for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
130 if (faceType[faceI] == 1)
132 internalFaces_[nInternal++] = faceI;
135 internalFaces_.setSize(nInternal);
137 labelList nIncludedFaces(patches.size(), 0);
138 labelList nExcludedFaces(patches.size(), 0);
140 forAll(patches, patchi)
142 const polyPatch& pp = patches[patchi];
144 forAll(pp, patchFacei)
146 label faceI = pp.start()+patchFacei;
148 if (faceType[faceI] == 1)
150 nIncludedFaces[patchi]++;
152 else if (faceType[faceI] == 2)
154 nExcludedFaces[patchi]++;
159 includedFaces_.setSize(patches.size());
160 excludedFaces_.setSize(patches.size());
161 forAll(nIncludedFaces, patchi)
163 includedFaces_[patchi].setSize(nIncludedFaces[patchi]);
164 excludedFaces_[patchi].setSize(nExcludedFaces[patchi]);
169 forAll(patches, patchi)
171 const polyPatch& pp = patches[patchi];
173 forAll(pp, patchFacei)
175 label faceI = pp.start()+patchFacei;
177 if (faceType[faceI] == 1)
179 includedFaces_[patchi][nIncludedFaces[patchi]++] = patchFacei;
181 else if (faceType[faceI] == 2)
183 excludedFaces_[patchi][nExcludedFaces[patchi]++] = patchFacei;
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();
197 faceSet MRFFaces(mesh_, "includedFaces", 100);
198 forAll(includedFaces_, patchi)
200 forAll(includedFaces_[patchi], i)
202 label patchFacei = includedFaces_[patchi][i];
203 MRFFaces.insert(patches[patchi].start()+patchFacei);
206 Pout<< "Writing " << MRFFaces.size()
207 << " patch faces in MRF zone to faceSet "
208 << MRFFaces.name() << endl;
211 faceSet excludedFaces(mesh_, "excludedFaces", 100);
212 forAll(excludedFaces_, patchi)
214 forAll(excludedFaces_[patchi], i)
216 label patchFacei = excludedFaces_[patchi][i];
217 excludedFaces.insert(patches[patchi].start()+patchFacei);
220 Pout<< "Writing " << excludedFaces.size()
221 << " faces in MRF zone with special handling to faceSet "
222 << excludedFaces.name() << endl;
223 excludedFaces.write();
228 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
230 Foam::MRFZone::MRFZone(const fvMesh& mesh, Istream& is)
235 cellZoneID_(mesh_.cellZones().findZoneID(name_)),
238 dict_.lookupOrDefault("nonRotatingPatches", wordList(0))
240 origin_(dict_.lookup("origin")),
241 axis_(dict_.lookup("axis")),
242 omega_(dict_.lookup("omega")),
243 Omega_("Omega", omega_*axis_)
245 if (dict_.found("patches"))
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'."
254 const polyBoundaryMesh& patches = mesh_.boundaryMesh();
256 axis_ = axis_/mag(axis_);
257 Omega_ = omega_*axis_;
259 excludedPatchLabels_.setSize(excludedPatchNames_.size());
261 forAll(excludedPatchNames_, i)
263 excludedPatchLabels_[i] = patches.findPatchID(excludedPatchNames_[i]);
265 if (excludedPatchLabels_[i] == -1)
269 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
270 ) << "cannot find MRF patch " << excludedPatchNames_[i]
276 bool cellZoneFound = (cellZoneID_ != -1);
277 reduce(cellZoneFound, orOp<bool>());
283 "Foam::MRFZone::MRFZone(const fvMesh&, Istream&)"
284 ) << "cannot find MRF cellZone " << name_
292 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
294 void Foam::MRFZone::addCoriolis(fvVectorMatrix& UEqn) const
296 if (cellZoneID_ == -1)
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();
309 label celli = cells[i];
310 Usource[celli] -= V[celli]*(Omega ^ U[celli]);
315 void Foam::MRFZone::addCoriolis
317 const volScalarField& rho,
321 if (cellZoneID_ == -1)
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();
334 label celli = cells[i];
335 Usource[celli] -= V[celli]*rho[celli]*(Omega ^ U[celli]);
340 void Foam::MRFZone::relativeVelocity(volVectorField& U) const
342 const volVectorField& C = mesh_.C();
344 const vector& origin = origin_.value();
345 const vector& Omega = Omega_.value();
347 const labelList& cells = mesh_.cellZones()[cellZoneID_];
351 label celli = cells[i];
352 U[celli] -= (Omega ^ (C[celli] - origin));
356 forAll(includedFaces_, patchi)
358 forAll(includedFaces_[patchi], i)
360 label patchFacei = includedFaces_[patchi][i];
361 U.boundaryField()[patchi][patchFacei] = vector::zero;
366 forAll(excludedFaces_, patchi)
368 forAll(excludedFaces_[patchi], i)
370 label patchFacei = excludedFaces_[patchi][i];
371 U.boundaryField()[patchi][patchFacei] -=
372 (Omega ^ (C.boundaryField()[patchi][patchFacei] - origin));
378 void Foam::MRFZone::relativeFlux(surfaceScalarField& phi) const
380 const surfaceVectorField& Cf = mesh_.Cf();
381 const surfaceVectorField& Sf = mesh_.Sf();
383 const vector& origin = origin_.value();
384 const vector& Omega = Omega_.value();
387 forAll(internalFaces_, i)
389 label facei = internalFaces_[i];
390 phi[facei] -= (Omega ^ (Cf[facei] - origin)) & Sf[facei];
394 forAll(includedFaces_, patchi)
396 forAll(includedFaces_[patchi], i)
398 label patchFacei = includedFaces_[patchi][i];
400 phi.boundaryField()[patchi][patchFacei] = 0.0;
405 forAll(excludedFaces_, patchi)
407 forAll(excludedFaces_[patchi], i)
409 label patchFacei = excludedFaces_[patchi][i];
411 phi.boundaryField()[patchi][patchFacei] -=
412 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
413 & Sf.boundaryField()[patchi][patchFacei];
419 void Foam::MRFZone::absoluteFlux(surfaceScalarField& phi) const
421 const surfaceVectorField& Cf = mesh_.Cf();
422 const surfaceVectorField& Sf = mesh_.Sf();
424 const vector& origin = origin_.value();
425 const vector& Omega = Omega_.value();
428 forAll(internalFaces_, i)
430 label facei = internalFaces_[i];
431 phi[facei] += (Omega ^ (Cf[facei] - origin)) & Sf[facei];
435 forAll(includedFaces_, patchi)
437 forAll(includedFaces_[patchi], i)
439 label patchFacei = includedFaces_[patchi][i];
441 phi.boundaryField()[patchi][patchFacei] +=
442 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
443 & Sf.boundaryField()[patchi][patchFacei];
448 forAll(excludedFaces_, patchi)
450 forAll(excludedFaces_[patchi], i)
452 label patchFacei = excludedFaces_[patchi][i];
454 phi.boundaryField()[patchi][patchFacei] +=
455 (Omega ^ (Cf.boundaryField()[patchi][patchFacei] - origin))
456 & Sf.boundaryField()[patchi][patchFacei];
462 void Foam::MRFZone::correctBoundaryVelocity(volVectorField& U) const
464 const vector& origin = origin_.value();
465 const vector& Omega = Omega_.value();
468 forAll(includedFaces_, patchi)
470 const vectorField& patchC = mesh_.Cf().boundaryField()[patchi];
472 vectorField pfld(U.boundaryField()[patchi]);
474 forAll(includedFaces_[patchi], i)
476 label patchFacei = includedFaces_[patchi][i];
478 pfld[patchFacei] = (Omega ^ (patchC[patchFacei] - origin));
481 U.boundaryField()[patchi] == pfld;
486 // ************************ vim: set sw=4 sts=4 et: ************************ //