From 7e3b60a9ea9f3b8679c0ef19da24b886fa8423e0 Mon Sep 17 00:00:00 2001 From: sergio Date: Mon, 7 Mar 2011 12:58:02 +0000 Subject: [PATCH 1/9] ENH: add of activePressureForceBaffleVelocity BC --- src/finiteVolume/Make/files | 1 + ...ureForceBaffleVelocityFvPatchVectorField.C | 385 ++++++++++++++++++ ...ureForceBaffleVelocityFvPatchVectorField.H | 207 ++++++++++ 3 files changed, 593 insertions(+) create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C create mode 100644 src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index e6f6c26148..9eaab89c40 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -112,6 +112,7 @@ $(constraintFvPatchFields)/wedge/wedgeFvPatchScalarField.C derivedFvPatchFields = $(fvPatchFields)/derived $(derivedFvPatchFields)/activeBaffleVelocity/activeBaffleVelocityFvPatchVectorField.C +$(derivedFvPatchFields)/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C $(derivedFvPatchFields)/advective/advectiveFvPatchFields.C $(derivedFvPatchFields)/codedFixedValue/codedFixedValueFvPatchScalarField.C diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C new file mode 100644 index 0000000000..245bf5662d --- /dev/null +++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.C @@ -0,0 +1,385 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "activePressureForceBaffleVelocityFvPatchVectorField.H" +#include "addToRunTimeSelectionTable.H" +#include "volFields.H" +#include "surfaceFields.H" +#include "cyclicFvPatch.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +activePressureForceBaffleVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF +) +: + fixedValueFvPatchVectorField(p, iF), + pName_("p"), + cyclicPatchName_(), + cyclicPatchLabel_(-1), + orientation_(1), + initWallSf_(0), + initCyclicSf_(0), + nbrCyclicSf_(0), + openFraction_(0), + openingTime_(0), + maxOpenFractionDelta_(0), + curTimeIndex_(-1), + minThresholdValue_(0), + fBased_(1), + baffleActivated_(0) +{} + + +Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +activePressureForceBaffleVelocityFvPatchVectorField +( + const activePressureForceBaffleVelocityFvPatchVectorField& ptf, + const fvPatch& p, + const DimensionedField& iF, + const fvPatchFieldMapper& mapper +) +: + fixedValueFvPatchVectorField(ptf, p, iF, mapper), + pName_(ptf.pName_), + cyclicPatchName_(ptf.cyclicPatchName_), + cyclicPatchLabel_(ptf.cyclicPatchLabel_), + orientation_(ptf.orientation_), + initWallSf_(ptf.initWallSf_), + initCyclicSf_(ptf.initCyclicSf_), + nbrCyclicSf_(ptf.nbrCyclicSf_), + openFraction_(ptf.openFraction_), + openingTime_(ptf.openingTime_), + maxOpenFractionDelta_(ptf.maxOpenFractionDelta_), + curTimeIndex_(-1), + minThresholdValue_(ptf.minThresholdValue_), + fBased_(ptf.fBased_), + baffleActivated_(ptf.baffleActivated_) +{} + + +Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +activePressureForceBaffleVelocityFvPatchVectorField +( + const fvPatch& p, + const DimensionedField& iF, + const dictionary& dict +) +: + fixedValueFvPatchVectorField(p, iF), + pName_("p"), + cyclicPatchName_(dict.lookup("cyclicPatch")), + cyclicPatchLabel_(p.patch().boundaryMesh().findPatchID(cyclicPatchName_)), + orientation_(readLabel(dict.lookup("orientation"))), + initWallSf_(0), + initCyclicSf_(0), + nbrCyclicSf_(0), + openFraction_(readScalar(dict.lookup("openFraction"))), + openingTime_(readScalar(dict.lookup("openingTime"))), + maxOpenFractionDelta_(readScalar(dict.lookup("maxOpenFractionDelta"))), + curTimeIndex_(-1), + minThresholdValue_(readScalar(dict.lookup("minThresholdValue"))), + fBased_(readBool(dict.lookup("forceBased"))), + baffleActivated_(0) +{ + fvPatchVectorField::operator=(vector::zero); + + if (p.size() > 0) + { + initWallSf_ = p.Sf(); + initCyclicSf_ = p.boundaryMesh()[cyclicPatchLabel_].Sf(); + nbrCyclicSf_ = refCast + ( + p.boundaryMesh()[cyclicPatchLabel_] + ).neighbFvPatch().Sf(); + } + + if (dict.found("p")) + { + dict.lookup("p") >> pName_; + } +} + + +Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +activePressureForceBaffleVelocityFvPatchVectorField +( + const activePressureForceBaffleVelocityFvPatchVectorField& ptf +) +: + fixedValueFvPatchVectorField(ptf), + pName_(ptf.pName_), + cyclicPatchName_(ptf.cyclicPatchName_), + cyclicPatchLabel_(ptf.cyclicPatchLabel_), + orientation_(ptf.orientation_), + initWallSf_(ptf.initWallSf_), + initCyclicSf_(ptf.initCyclicSf_), + nbrCyclicSf_(ptf.nbrCyclicSf_), + openFraction_(ptf.openFraction_), + openingTime_(ptf.openingTime_), + maxOpenFractionDelta_(ptf.maxOpenFractionDelta_), + curTimeIndex_(-1), + minThresholdValue_(ptf.minThresholdValue_), + fBased_(ptf.fBased_), + baffleActivated_(ptf.baffleActivated_) +{} + + +Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +activePressureForceBaffleVelocityFvPatchVectorField +( + const activePressureForceBaffleVelocityFvPatchVectorField& ptf, + const DimensionedField& iF +) +: + fixedValueFvPatchVectorField(ptf, iF), + pName_(ptf.pName_), + cyclicPatchName_(ptf.cyclicPatchName_), + cyclicPatchLabel_(ptf.cyclicPatchLabel_), + orientation_(ptf.orientation_), + initWallSf_(ptf.initWallSf_), + initCyclicSf_(ptf.initCyclicSf_), + nbrCyclicSf_(ptf.nbrCyclicSf_), + openFraction_(ptf.openFraction_), + openingTime_(ptf.openingTime_), + maxOpenFractionDelta_(ptf.maxOpenFractionDelta_), + curTimeIndex_(-1), + minThresholdValue_(ptf.minThresholdValue_), + fBased_(ptf.fBased_), + baffleActivated_(ptf.baffleActivated_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::activePressureForceBaffleVelocityFvPatchVectorField::autoMap +( + const fvPatchFieldMapper& m +) +{ + fixedValueFvPatchVectorField::autoMap(m); + + //- Note: cannot map field from cyclic patch anyway so just recalculate + // Areas should be consistent when doing autoMap except in case of + // topo changes. + //- Note: we don't want to use Sf here since triggers rebuilding of + // fvMesh::S() which will give problems when mapped (since already + // on new mesh) + forAll (patch().boundaryMesh().mesh().faceAreas(), i) + { + if (mag(patch().boundaryMesh().mesh().faceAreas()[i]) == 0) + { + Info << "faceArea[active] "<< i << endl; + } + } + if (patch().size() > 0) + { + const vectorField& areas = patch().boundaryMesh().mesh().faceAreas(); + initWallSf_ = patch().patchSlice(areas); + initCyclicSf_ = patch().boundaryMesh() + [ + cyclicPatchLabel_ + ].patchSlice(areas); + nbrCyclicSf_ = refCast + ( + patch().boundaryMesh() + [ + cyclicPatchLabel_ + ] + ).neighbFvPatch().patch().patchSlice(areas); + } +} + +void Foam::activePressureForceBaffleVelocityFvPatchVectorField::rmap +( + const fvPatchVectorField& ptf, + const labelList& addr +) +{ + fixedValueFvPatchVectorField::rmap(ptf, addr); + + // See autoMap. + const vectorField& areas = patch().boundaryMesh().mesh().faceAreas(); + initWallSf_ = patch().patchSlice(areas); + initCyclicSf_ = patch().boundaryMesh() + [ + cyclicPatchLabel_ + ].patchSlice(areas); + nbrCyclicSf_ = refCast + ( + patch().boundaryMesh() + [ + cyclicPatchLabel_ + ] + ).neighbFvPatch().patch().patchSlice(areas); +} + + +void Foam::activePressureForceBaffleVelocityFvPatchVectorField::updateCoeffs() +{ + if (updated()) + { + return; + } + // Execute the change to the openFraction only once per time-step + if (curTimeIndex_ != this->db().time().timeIndex()) + { + const volScalarField& p = db().lookupObject + ( + pName_ + ); + + const fvPatch& cyclicPatch = patch().boundaryMesh()[cyclicPatchLabel_]; + const labelList& cyclicFaceCells = cyclicPatch.patch().faceCells(); + const fvPatch& nbrPatch = refCast + ( + cyclicPatch + ).neighbFvPatch(); + + const labelList& nbrFaceCells = nbrPatch.patch().faceCells(); + + scalar valueDiff = 0; + + if (fBased_) + { + // Add this side + forAll(cyclicFaceCells, facei) + { + valueDiff +=p[cyclicFaceCells[facei]]*mag(initCyclicSf_[facei]); + } + + // Remove other side + forAll(nbrFaceCells, facei) + { + valueDiff -=p[nbrFaceCells[facei]]*mag(initCyclicSf_[facei]); + } + } + else //pressure based + { + forAll(cyclicFaceCells, facei) + { + valueDiff += p[cyclicFaceCells[facei]]; + } + + forAll(nbrFaceCells, facei) + { + valueDiff -= p[nbrFaceCells[facei]]; + } + } + + if ((mag(valueDiff) > mag(minThresholdValue_) || baffleActivated_)) + { + openFraction_ = + max( + min( + openFraction_ + + max + ( + this->db().time().deltaT().value()/openingTime_, + maxOpenFractionDelta_ + )*(orientation_), + 1 - 1e-6 + ), + 1e-6 + ); + + baffleActivated_ = true; + } + else + { + openFraction_ = max(min(1 - 1e-6, openFraction_), 1e-6); + } + + Info<< "Open fraction = " << openFraction_ << endl; + Info<< "Pressure difference = " << valueDiff << endl; + + vectorField::subField Sfw = patch().patch().faceAreas(); + vectorField newSfw = (1 - openFraction_)*initWallSf_; + forAll(Sfw, facei) + { + Sfw[facei] = newSfw[facei]; + } + const_cast(patch().magSf()) = mag(patch().Sf()); + + // Update owner side of cyclic + const_cast(cyclicPatch.Sf()) = + openFraction_*initCyclicSf_; + + const_cast(cyclicPatch.magSf()) = + mag(cyclicPatch.Sf()); + + // Update neighbour side of cyclic + const_cast(nbrPatch.Sf()) = + openFraction_*nbrCyclicSf_; + + const_cast(nbrPatch.magSf()) = + mag(nbrPatch.Sf()); + + curTimeIndex_ = this->db().time().timeIndex(); + } + + fixedValueFvPatchVectorField::updateCoeffs(); +} + + +void Foam::activePressureForceBaffleVelocityFvPatchVectorField:: +write(Ostream& os) const +{ + fvPatchVectorField::write(os); + os.writeKeyword("cyclicPatch") + << cyclicPatchName_ << token::END_STATEMENT << nl; + os.writeKeyword("orientation") + << orientation_ << token::END_STATEMENT << nl; + os.writeKeyword("openingTime") + << openingTime_ << token::END_STATEMENT << nl; + os.writeKeyword("maxOpenFractionDelta") + << maxOpenFractionDelta_ << token::END_STATEMENT << nl; + os.writeKeyword("openFraction") + << openFraction_ << token::END_STATEMENT << nl; + os.writeKeyword("p") + << pName_ << token::END_STATEMENT << nl; + os.writeKeyword("minThresholdValue") + << minThresholdValue_ << token::END_STATEMENT << nl; + os.writeKeyword("forceBased") + << fBased_ << token::END_STATEMENT << nl; + writeEntry("value", os); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + makePatchTypeField + ( + fvPatchVectorField, + activePressureForceBaffleVelocityFvPatchVectorField + ); +} + +// ************************************************************************* // diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H new file mode 100644 index 0000000000..aea26d662b --- /dev/null +++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H @@ -0,0 +1,207 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::activePressureForceBaffleVelocityFvPatchVectorField + +Description + Bounday which emulates the operation of a release pressure panel. + + The boundary condition modifies mesh areas based on difference + of pressure or force face beween both sides of the panel. Once opened the + panel continues to open at a fixed rate. + +SourceFiles + activePressureForceBaffleVelocityFvPatchVectorField.C + +\*---------------------------------------------------------------------------*/ + +#ifndef activePressureForceBaffleVelocityFvPatchVectorField_H +#define activePressureForceBaffleVelocityFvPatchVectorField_H + +#include "fvPatchFields.H" +#include "fixedValueFvPatchFields.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class activePressureForceBaffleVelocityFvPatch Declaration +\*---------------------------------------------------------------------------*/ + +class activePressureForceBaffleVelocityFvPatchVectorField +: + public fixedValueFvPatchVectorField +{ + // Private data + + //- Name of the pressure field used to calculate the force + // on the active baffle + word pName_; + + //- Name of the cyclic patch used when the active baffle is open + word cyclicPatchName_; + + //- Index of the cyclic patch used when the active baffle is open + label cyclicPatchLabel_; + + //- Orientation (1 or -1) of the active baffle patch. + // Used to change the direction of opening without the need for + // reordering the patch faces + label orientation_; + + //- Initial wall patch areas + vectorField initWallSf_; + + //- Initial cyclic patch areas + vectorField initCyclicSf_; + + //- Initial neighbour-side cyclic patch areas + vectorField nbrCyclicSf_; + + //- Current fraction of the active baffle which is open + scalar openFraction_; + + //- Time taken for the active baffle to open + scalar openingTime_; + + //- Maximum fractional change to the active baffle openness + // per time-step + scalar maxOpenFractionDelta_; + + label curTimeIndex_; + + //- Minimum value for the active baffle to start opening + scalar minThresholdValue_; + + //- Force based active baffle + bool fBased_; + + //- Baffle is activated + bool baffleActivated_; + +public: + + //- Runtime type information + TypeName("activePressureForceBaffleVelocity"); + + + // Constructors + + //- Construct from patch and internal field + activePressureForceBaffleVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField& + ); + + //- Construct from patch, internal field and dictionary + activePressureForceBaffleVelocityFvPatchVectorField + ( + const fvPatch&, + const DimensionedField&, + const dictionary& + ); + + //- Construct by mapping given activePressureForceBaffleVelocityFvPatchVectorField + // onto a new patch + activePressureForceBaffleVelocityFvPatchVectorField + ( + const activePressureForceBaffleVelocityFvPatchVectorField&, + const fvPatch&, + const DimensionedField&, + const fvPatchFieldMapper& + ); + + //- Construct as copy + activePressureForceBaffleVelocityFvPatchVectorField + ( + const activePressureForceBaffleVelocityFvPatchVectorField& + ); + + //- Construct and return a clone + virtual tmp clone() const + { + return tmp + ( + new activePressureForceBaffleVelocityFvPatchVectorField(*this) + ); + } + + //- Construct as copy setting internal field reference + activePressureForceBaffleVelocityFvPatchVectorField + ( + const activePressureForceBaffleVelocityFvPatchVectorField&, + const DimensionedField& + ); + + //- Construct and return a clone setting internal field reference + virtual tmp clone + ( + const DimensionedField& iF + ) const + { + return tmp + ( + new activePressureForceBaffleVelocityFvPatchVectorField(*this, iF) + ); + } + + + // Member functions + + // Mapping functions + + //- Map (and resize as needed) from self given a mapping object + virtual void autoMap + ( + const fvPatchFieldMapper& + ); + + //- Reverse map the given fvPatchField onto this fvPatchField + virtual void rmap + ( + const fvPatchVectorField&, + const labelList& + ); + + + //- Update the coefficients associated with the patch field + virtual void updateCoeffs(); + + //- Write + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // From b3adb680f841076181a4eb1d482829192794c462 Mon Sep 17 00:00:00 2001 From: sergio Date: Mon, 7 Mar 2011 13:00:31 +0000 Subject: [PATCH 2/9] STY: max lines checks --- ...activePressureForceBaffleVelocityFvPatchVectorField.H | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H index aea26d662b..538049b34f 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/activePressureForceBaffleVelocity/activePressureForceBaffleVelocityFvPatchVectorField.H @@ -125,8 +125,7 @@ public: const dictionary& ); - //- Construct by mapping given activePressureForceBaffleVelocityFvPatchVectorField - // onto a new patch + //- Construct by mapping activePressureForceBaffleVelocityFvPatchVectorField ( const activePressureForceBaffleVelocityFvPatchVectorField&, @@ -165,7 +164,11 @@ public: { return tmp ( - new activePressureForceBaffleVelocityFvPatchVectorField(*this, iF) + new activePressureForceBaffleVelocityFvPatchVectorField + ( + *this, + iF + ) ); } From e95597a2c559dc746100e940e6d653ea9d06032f Mon Sep 17 00:00:00 2001 From: sergio Date: Mon, 7 Mar 2011 13:19:37 +0000 Subject: [PATCH 3/9] ENH: PDRMesh --- .../mesh/advanced/PDRMesh/Make/files | 3 + .../mesh/advanced/PDRMesh/Make/options | 11 + .../utilities/mesh/advanced/PDRMesh/PDRMesh.C | 1181 +++++++++++++++++ .../mesh/advanced/PDRMesh/PDRMeshDict | 40 + 4 files changed, 1235 insertions(+) create mode 100644 applications/utilities/mesh/advanced/PDRMesh/Make/files create mode 100644 applications/utilities/mesh/advanced/PDRMesh/Make/options create mode 100644 applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C create mode 100644 applications/utilities/mesh/advanced/PDRMesh/PDRMeshDict diff --git a/applications/utilities/mesh/advanced/PDRMesh/Make/files b/applications/utilities/mesh/advanced/PDRMesh/Make/files new file mode 100644 index 0000000000..998d5028f2 --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/Make/files @@ -0,0 +1,3 @@ +PDRMesh.C + +EXE = $(FOAM_APPBIN)/PDRMesh diff --git a/applications/utilities/mesh/advanced/PDRMesh/Make/options b/applications/utilities/mesh/advanced/PDRMesh/Make/options new file mode 100644 index 0000000000..6ab51e876a --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/Make/options @@ -0,0 +1,11 @@ +EXE_INC = \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -lmeshTools \ + -ldynamicMesh \ + -lfiniteVolume \ + -lcompressibleRASModels + diff --git a/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C new file mode 100644 index 0000000000..c69e102bc2 --- /dev/null +++ b/applications/utilities/mesh/advanced/PDRMesh/PDRMesh.C @@ -0,0 +1,1181 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd. + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software; you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by the + Free Software Foundation; either version 3 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Description + Mesh and field preparation utility for PDR type simulations. + + Reads + - cellSet giving blockedCells + - faceSets giving blockedFaces and the patch they should go into + + NOTE: To avoid exposing wrong fields values faceSets should include + faces contained in the blockedCells cellset. + + - coupledFaces reads coupledFacesSet to introduces mixe-coupled baffles + + Subsets out the blocked cells and splits the blockedFaces and updates + fields. + + The face splitting is done by duplicating the faces. No duplication of + points for now (so checkMesh will show a lot of 'duplicate face' messages) + +\*---------------------------------------------------------------------------*/ + +#include "fvMeshSubset.H" +#include "argList.H" +#include "cellSet.H" +#include "IOobjectList.H" +#include "volFields.H" +#include "mapPolyMesh.H" +#include "faceSet.H" +#include "cellSet.H" +#include "syncTools.H" +#include "polyTopoChange.H" +#include "polyModifyFace.H" +#include "polyAddFace.H" +#include "regionSplit.H" +#include "Tuple2.H" +#include "cyclicFvPatch.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void modifyOrAddFace +( + polyTopoChange& meshMod, + const face& f, + const label faceI, + const label own, + const bool flipFaceFlux, + const label newPatchI, + const label zoneID, + const bool zoneFlip, + + PackedBoolList& modifiedFace +) +{ + if (!modifiedFace[faceI]) + { + // First usage of face. Modify. + meshMod.setAction + ( + polyModifyFace + ( + f, // modified face + faceI, // label of face + own, // owner + -1, // neighbour + flipFaceFlux, // face flip + newPatchI, // patch for face + false, // remove from zone + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + modifiedFace[faceI] = 1; + } + else + { + // Second or more usage of face. Add. + meshMod.setAction + ( + polyAddFace + ( + f, // modified face + own, // owner + -1, // neighbour + -1, // master point + -1, // master edge + faceI, // master face + flipFaceFlux, // face flip + newPatchI, // patch for face + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + } +} + + +template +void subsetVolFields +( + const fvMeshSubset& subsetter, + const IOobjectList& objectsList, + const label patchI, + const Type& exposedValue, + const word GeomVolType, + PtrList >& subFields +) +{ + const fvMesh& baseMesh = subsetter.baseMesh(); + + label i = 0; + + forAllConstIter(IOobjectList , objectsList, iter) + { + if (iter()->headerClassName() == GeomVolType) + { + const word fieldName = iter()->name(); + + Info<< "Subsetting field " << fieldName << endl; + + GeometricField volField + ( + *iter(), + baseMesh + ); + + subFields.set(i, subsetter.interpolate(volField)); + + // Explicitly set exposed faces (in patchI) to exposedValue. + if (patchI >= 0) + { + fvPatchField& fld = + subFields[i++].boundaryField()[patchI]; + + label newStart = fld.patch().patch().start(); + + label oldPatchI = subsetter.patchMap()[patchI]; + + if (oldPatchI == -1) + { + // New patch. Reset whole value. + fld = exposedValue; + } + else + { + // Reset those faces that originate from different patch + // or internal faces. + label oldSize = volField.boundaryField()[oldPatchI].size(); + label oldStart = volField.boundaryField() + [ + oldPatchI + ].patch().patch().start(); + + forAll(fld, j) + { + label oldFaceI = subsetter.faceMap()[newStart+j]; + + if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize) + { + fld[j] = exposedValue; + } + } + } + } + } + } +} + + +template +void subsetSurfaceFields +( + const fvMeshSubset& subsetter, + const IOobjectList& objectsList, + const label patchI, + const Type& exposedValue, + const word GeomSurfType, + PtrList >& subFields +) +{ + const fvMesh& baseMesh = subsetter.baseMesh(); + + label i(0); + + forAllConstIter(IOobjectList , objectsList, iter) + { + if (iter()->headerClassName() == GeomSurfType) + { + const word& fieldName = iter.key(); + + Info<< "Subsetting field " << fieldName << endl; + + GeometricField volField + ( + *iter(), + baseMesh + ); + + subFields.set(i, subsetter.interpolate(volField)); + + + // Explicitly set exposed faces (in patchI) to exposedValue. + if (patchI >= 0) + { + fvsPatchField& fld = + subFields[i++].boundaryField()[patchI]; + + label newStart = fld.patch().patch().start(); + + label oldPatchI = subsetter.patchMap()[patchI]; + + if (oldPatchI == -1) + { + // New patch. Reset whole value. + fld = exposedValue; + } + else + { + // Reset those faces that originate from different patch + // or internal faces. + label oldSize = volField.boundaryField()[oldPatchI].size(); + label oldStart = volField.boundaryField() + [ + oldPatchI + ].patch().patch().start(); + + forAll(fld, j) + { + label oldFaceI = subsetter.faceMap()[newStart+j]; + + if(oldFaceI < oldStart || oldFaceI >= oldStart+oldSize) + { + fld[j] = exposedValue; + } + } + } + } + } + } +} + + +// Faces introduced into zero-sized patches don't get a value at all. +// This is hack to set them to an initial value. +template +void initCreatedPatches +( + const fvMesh& mesh, + const mapPolyMesh& map, + const typename GeoField::value_type initValue +) +{ + HashTable fields + ( + mesh.objectRegistry::lookupClass() + ); + + for + ( + typename HashTable:: + iterator fieldIter = fields.begin(); + fieldIter != fields.end(); + ++fieldIter + ) + { + GeoField& field = const_cast(*fieldIter()); + + forAll(field.boundaryField(), patchi) + { + if (map.oldPatchSizes()[patchi] == 0) + { + // Not mapped. + field.boundaryField()[patchi] = initValue; + + if (field.boundaryField()[patchi].fixesValue()) + { + field.boundaryField()[patchi] == initValue; + } + } + } + } +} + + +void createCoupledBaffles +( + fvMesh& mesh, + const labelList& coupledWantedPatch, + polyTopoChange& meshMod, + PackedBoolList& modifiedFace +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + + forAll(coupledWantedPatch, faceI) + { + if (coupledWantedPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + // Use owner side of face + modifyOrAddFace + ( + meshMod, + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + false, // face flip + coupledWantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + // Use neighbour side of face + modifyOrAddFace + ( + meshMod, + f.reverseFace(), // modified face + faceI, // label of face + mesh.faceNeighbour()[faceI],// owner + false, // face flip + coupledWantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + } + } + } +} + + +void createCyclicCoupledBaffles +( + fvMesh& mesh, + const labelList& cyclicMasterPatch, + const labelList& cyclicSlavePatch, + polyTopoChange& meshMod, + PackedBoolList& modifiedFace +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + + forAll(cyclicMasterPatch, faceI) + { + if (cyclicMasterPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + modifyOrAddFace + ( + meshMod, + f.reverseFace(), // modified face + faceI, // label of face + mesh.faceNeighbour()[faceI], // owner + false, // face flip + cyclicMasterPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add + ); + } + } + + forAll(cyclicSlavePatch, faceI) + { + if (cyclicSlavePatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + // Use owner side of face + modifyOrAddFace + ( + meshMod, + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + false, // face flip + cyclicSlavePatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip, // face flip in zone + modifiedFace // modify or add status + ); + } + } + } +} + + +void createBaffles +( + fvMesh& mesh, + const labelList& wantedPatch, + polyTopoChange& meshMod +) +{ + const faceZoneMesh& faceZones = mesh.faceZones(); + Info << "faceZone:createBaffle " << faceZones << endl; + forAll(wantedPatch, faceI) + { + if (wantedPatch[faceI] != -1) + { + const face& f = mesh.faces()[faceI]; + + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + meshMod.setAction + ( + polyModifyFace + ( + f, // modified face + faceI, // label of face + mesh.faceOwner()[faceI], // owner + -1, // neighbour + false, // face flip + wantedPatch[faceI], // patch for face + false, // remove from zone + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + + if (mesh.isInternalFace(faceI)) + { + label zoneID = faceZones.whichZone(faceI); + bool zoneFlip = false; + + if (zoneID >= 0) + { + const faceZone& fZone = faceZones[zoneID]; + zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)]; + } + + meshMod.setAction + ( + polyAddFace + ( + f.reverseFace(), // modified face + mesh.faceNeighbour()[faceI],// owner + -1, // neighbour + -1, // masterPointID + -1, // masterEdgeID + faceI, // masterFaceID, + false, // face flip + wantedPatch[faceI], // patch for face + zoneID, // zone for face + zoneFlip // face flip in zone + ) + ); + } + } + } +} + + +// Wrapper around find patch. Also makes sure same patch in parallel. +label findPatch(const polyBoundaryMesh& patches, const word& patchName) +{ + label patchI = patches.findPatchID(patchName); + + if (patchI == -1) + { + FatalErrorIn("findPatch(const polyBoundaryMesh&, const word&)") + << "Illegal patch " << patchName + << nl << "Valid patches are " << patches.names() + << exit(FatalError); + } + + // Check same patch for all procs + { + label newPatch = patchI; + reduce(newPatch, minOp