diff --git a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H index 6c64056be2..4d7334c5f8 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/pEqn.H @@ -1,98 +1,102 @@ -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); - -volScalarField rAU(1.0/UEqn().A()); -volVectorField HbyA("HbyA", U); -HbyA = rAU*UEqn().H(); - -UEqn.clear(); - -bool closedVolume = false; - -if (simple.transonic()) -{ - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) - ); - - while (simple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvm::div(phid, p) - - fvm::laplacian(rho*rAU, p) - ); - - // Relax the pressure equation to ensure diagonal-dominance - pEqn.relax(); - - pEqn.setReference(pRefCell, pRefValue); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi == pEqn.flux(); - } - } -} -else -{ - surfaceScalarField phiHbyA - ( - "phiHbyA", - fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) - ); - - closedVolume = adjustPhi(phiHbyA, U, p); - - while (simple.correctNonOrthogonal()) - { - fvScalarMatrix pEqn - ( - fvc::div(phiHbyA) - - fvm::laplacian(rho*rAU, p) - ); - - pEqn.setReference(pRefCell, pRefValue); - - pEqn.solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi = phiHbyA + pEqn.flux(); - } - } -} - - -#include "incompressible/continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -p.relax(); - -U = HbyA - rAU*fvc::grad(p); -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); -} - -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); - -if (!simple.transonic()) { + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); rho.relax(); -} -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; + volScalarField rAU(1.0/UEqn().A()); + volVectorField HbyA("HbyA", U); + HbyA = rAU*UEqn().H(); + + UEqn.clear(); + + bool closedVolume = false; + + if (simple.transonic()) + { + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::div(phid, p) + - fvm::laplacian(rho*rAU, p) + ); + + // Relax the pressure equation to ensure diagonal-dominance + pEqn.relax(); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi == pEqn.flux(); + } + } + } + else + { + surfaceScalarField phiHbyA + ( + "phiHbyA", + fvc::interpolate(rho)*(fvc::interpolate(HbyA) & mesh.Sf()) + ); + + closedVolume = adjustPhi(phiHbyA, U, p); + + while (simple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) + ); + + pEqn.setReference(pRefCell, pRefValue); + + pEqn.solve(); + + if (simple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } + } + + + #include "incompressible/continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + U = HbyA - rAU*fvc::grad(p); + U.correctBoundaryConditions(); + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) + { + p += (initialMass - fvc::domainIntegrate(psi*p)) + /fvc::domainIntegrate(psi); + } + + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); + + if (!simple.transonic()) + { + rho.relax(); + } + + Info<< "rho max/min : " + << max(rho).value() << " " + << min(rho).value() << endl; +} diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H new file mode 100644 index 0000000000..96bdb2a122 --- /dev/null +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/createFields.H @@ -0,0 +1,61 @@ + Info<< "Reading thermophysical properties\n" << endl; + + autoPtr pThermo + ( + rhoThermo::New(mesh) + ); + rhoThermo& thermo = pThermo(); + + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + volScalarField& p = thermo.p(); + volScalarField& e = thermo.he(); + + Info<< "Reading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "compressibleCreatePhi.H" + + + label pRefCell = 0; + scalar pRefValue = 0.0; + setRefCell(p, simple.dict(), pRefCell, pRefValue); + + dimensionedScalar rhoMax(simple.dict().lookup("rhoMax")); + dimensionedScalar rhoMin(simple.dict().lookup("rhoMin")); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::RASModel::New + ( + rho, + U, + phi, + thermo + ) + ); + + dimensionedScalar initialMass = fvc::domainIntegrate(rho); diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H index 4791062d35..5d0f174623 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/eEqn.H @@ -1,4 +1,5 @@ { + // Kinetic + pressure energy volScalarField Ekp("Ekp", 0.5*magSqr(U) + p/rho); fvScalarMatrix eEqn @@ -10,7 +11,7 @@ fvc::div(phi)*Ekp - fvc::div(phi, Ekp) ); - //pZones.addEnergySource(thermo, rho, eEqn); + pZones.addEnergySource(thermo, rho, eEqn); eEqn.relax(); eEqn.solve(); diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H index c3e0a43a15..81b194823c 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/pEqn.H @@ -1,52 +1,24 @@ -volVectorField HbyA("HbyA", U); - -if (pressureImplicitPorosity) { - HbyA = trTU() & UEqn().H(); -} -else -{ - HbyA = trAU()*UEqn().H(); -} + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); + rho.relax(); -UEqn.clear(); + volVectorField HbyA("HbyA", U); -bool closedVolume = false; - -if (simple.transonic()) -{ - surfaceScalarField phid - ( - "phid", - fvc::interpolate(psi)*(fvc::interpolate(HbyA) & mesh.Sf()) - ); - mrfZones.relativeFlux(fvc::interpolate(psi), phid); - - while (simple.correctNonOrthogonal()) + if (pressureImplicitPorosity) { - tmp tpEqn; - - if (pressureImplicitPorosity) - { - tpEqn = (fvc::div(phid, p) - fvm::laplacian(rho*trTU(), p)); - } - else - { - tpEqn = (fvc::div(phid, p) - fvm::laplacian(rho*trAU(), p)); - } - - tpEqn().setReference(pRefCell, pRefValue); - - tpEqn().solve(); - - if (simple.finalNonOrthogonalIter()) - { - phi == tpEqn().flux(); - } + HbyA = trTU() & UEqn().H(); } -} -else -{ + else + { + HbyA = trAU()*UEqn().H(); + } + + UEqn.clear(); + + bool closedVolume = false; + surfaceScalarField phiHbyA ( "phiHbyA", @@ -79,34 +51,37 @@ else phi = phiHbyA - tpEqn().flux(); } } + + #include "incompressible/continuityErrs.H" + + // Explicitly relax pressure for momentum corrector + p.relax(); + + if (pressureImplicitPorosity) + { + U = HbyA - (trTU() & fvc::grad(p)); + } + else + { + U = HbyA - trAU()*fvc::grad(p); + } + + U.correctBoundaryConditions(); + + // For closed-volume cases adjust the pressure and density levels + // to obey overall mass continuity + if (closedVolume) + { + const volScalarField& psi = thermo.psi(); + p += (initialMass - fvc::domainIntegrate(psi*p)) + /fvc::domainIntegrate(psi); + } + + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); + rho.relax(); + Info<< "rho max/min : " + << max(rho).value() << " " + << min(rho).value() << endl; } - -#include "incompressible/continuityErrs.H" - -// Explicitly relax pressure for momentum corrector -p.relax(); - -if (pressureImplicitPorosity) -{ - U = HbyA - (trTU() & fvc::grad(p)); -} -else -{ - U = HbyA - trAU()*fvc::grad(p); -} - -U.correctBoundaryConditions(); - -// For closed-volume cases adjust the pressure and density levels -// to obey overall mass continuity -if (closedVolume) -{ - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); -} - -rho = thermo.rho(); -rho = max(rho, rhoMin); -rho = min(rho, rhoMax); -rho.relax(); -Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C index 9238106885..22a58da1cc 100644 --- a/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C +++ b/applications/solvers/compressible/rhoSimpleFoam/rhoPorousMRFSimpleFoam/rhoPorousMRFSimpleFoam.C @@ -32,7 +32,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "psiThermo.H" +#include "rhoThermo.H" #include "RASModel.H" #include "MRFZones.H" #include "thermalPorousZones.H" diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H index ca338975fb..9051805b87 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/porousFluid/pPorousFluidEqn.H @@ -16,7 +16,6 @@ porousPhi = for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix pEqn ( fvm::laplacian(porousRho*rAUPorous, porousP) == fvc::div(porousPhi) diff --git a/applications/test/BinSum/Make/files b/applications/test/BinSum/Make/files new file mode 100644 index 0000000000..34d265d30c --- /dev/null +++ b/applications/test/BinSum/Make/files @@ -0,0 +1,3 @@ +Test-BinSum.C + +EXE = $(FOAM_USER_APPBIN)/Test-BinSum diff --git a/applications/test/BinSum/Make/options b/applications/test/BinSum/Make/options new file mode 100644 index 0000000000..6a9e9810b3 --- /dev/null +++ b/applications/test/BinSum/Make/options @@ -0,0 +1,2 @@ +/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */ +/* EXE_LIBS = -lfiniteVolume */ diff --git a/applications/test/BinSum/Test-BinSum.C b/applications/test/BinSum/Test-BinSum.C new file mode 100644 index 0000000000..fb768f80c6 --- /dev/null +++ b/applications/test/BinSum/Test-BinSum.C @@ -0,0 +1,74 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 . + +Application + Test-BinSum + +Description + Test BinSum container + +\*---------------------------------------------------------------------------*/ + +#include "BinSum.H" +#include "IOstreams.H" +#include "Random.H" +#include "scalarField.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + Random rndGen(0); + + scalarField samples(10000000); + forAll(samples, i) + { + samples[i] = rndGen.scalar01(); + } + + const scalar min = 0; + const scalar max = 1; + const scalar delta = 0.1; + + BinSum count(min, max, delta); + BinSum sum(min, max, delta); + + forAll(samples, i) + { + count.add(samples[i], 1); + sum.add(samples[i], samples[i]); + } + + Info<< "sum : " << sum << endl; + Info<< "count : " << count << endl; + Info<< "average: " << sum/count << endl; + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H index cc22f92088..3be6dab81c 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H @@ -67,7 +67,6 @@ SourceFiles #include "transform.H" #include "volFields.H" #include "fvMesh.H" -#include "Histogram.H" #include "labelPair.H" #include "HashSet.H" #include "memInfo.H" diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 25405d2ff3..db67ba8a93 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -165,8 +165,12 @@ castellatedMeshControls } } + // Feature angle: + // - used if min and max refinement level of a surface differ + // - used if feature snapping (see snapControls below) is used resolveFeatureAngle 30; + // Region-wise refinement // ~~~~~~~~~~~~~~~~~~~~~~ diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSetDict b/applications/utilities/mesh/manipulation/topoSet/topoSetDict index 562c006ab1..be90a8fb80 100644 --- a/applications/utilities/mesh/manipulation/topoSet/topoSetDict +++ b/applications/utilities/mesh/manipulation/topoSet/topoSetDict @@ -171,6 +171,15 @@ FoamFile // insidePoint (1 2 3); // point inside region to select // } // +// // Cells underneath plane such that volume is reached. Can be used +// // in setFields. +// source targetVolumeToCell; +// sourceInfo +// { +// volume 2e-05; +// normal (1 1 1); +// } +// // // // faceSet diff --git a/src/OpenFOAM/containers/Lists/BinSum/BinSum.C b/src/OpenFOAM/containers/Lists/BinSum/BinSum.C new file mode 100644 index 0000000000..1e09961326 --- /dev/null +++ b/src/OpenFOAM/containers/Lists/BinSum/BinSum.C @@ -0,0 +1,113 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 "BinSum.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +Foam::BinSum::BinSum +( + const IndexType min, + const IndexType max, + const IndexType delta +) +: + List(ceil((max-min)/delta), pTraits::zero), + min_(min), + max_(max), + delta_(delta), + lowSum_(pTraits::zero), + highSum_(pTraits::zero) +{} + + +template +Foam::BinSum::BinSum +( + const IndexType min, + const IndexType max, + const IndexType delta, + const UList& indexVals, + const List& vals, + const CombineOp& cop +) +: + List(ceil((max-min)/delta), pTraits::zero), + min_(min), + max_(max), + delta_(delta), + lowSum_(pTraits::zero), + highSum_(pTraits::zero) +{ + forAll(indexVals, i) + { + add(indexVals[i], vals[i], cop); + } +} + + +// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * // + +template +void Foam::BinSum::add +( + const IndexType& indexVal, + const typename List::const_reference val, + const CombineOp& cop +) +{ + if (indexVal < min_) + { + cop(lowSum_, val); + } + else if (indexVal >= max_) + { + cop(highSum_, val); + } + else + { + label index = (indexVal-min_)/delta_; + cop(this->operator[](index), val); + } +} + + +template +void Foam::BinSum::add +( + const UList& indexVals, + const List& vals, + const CombineOp& cop +) +{ + forAll(indexVals, i) + { + add(indexVals[i], vals[i], cop); + } +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/containers/Lists/BinSum/BinSum.H b/src/OpenFOAM/containers/Lists/BinSum/BinSum.H new file mode 100644 index 0000000000..3cf09bc3a7 --- /dev/null +++ b/src/OpenFOAM/containers/Lists/BinSum/BinSum.H @@ -0,0 +1,149 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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::BinSum + +Description + Sums into bins + +SourceFiles + BinSum.C + +\*---------------------------------------------------------------------------*/ + +#ifndef BinSum_H +#define BinSum_H + +#include "ops.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + + +/*---------------------------------------------------------------------------*\ + Class BinSum Declaration +\*---------------------------------------------------------------------------*/ + +template +< + class IndexType, + class List, + class CombineOp = plusEqOp +> +class BinSum +: + public List +{ + // Private data + + const IndexType min_; + + const IndexType max_; + + const IndexType delta_; + + + //- Sum < lowest bin + typename List::value_type lowSum_; + + //- Sum of >= highest bin + typename List::value_type highSum_; + +public: + + // Constructors + + //- Construct given min, max, delta + BinSum + ( + const IndexType min, + const IndexType max, + const IndexType delta + ); + + //- Construct given min, max, delta and data + BinSum + ( + const IndexType min, + const IndexType max, + const IndexType delta, + const UList& indexVals, + const List& vals, + const CombineOp& cop = plusEqOp() + ); + + + // Access + + //- Return the delta + inline IndexType delta() const + { + return delta_; + } + + //- Return the sum of all added elements < min + inline const IndexType& lowSum() const + { + return lowSum_; + } + + //- Return the sum of all added elements >= max + inline const IndexType& highSum() const + { + return highSum_; + } + + void add + ( + const IndexType& indexVal, + const typename List::const_reference val, + const CombineOp& cop = plusEqOp() + ); + + void add + ( + const UList& indexVals, + const List& vals, + const CombineOp& cop = plusEqOp() + ); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "BinSum.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C index 35b18aaf69..6223661b81 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,33 +58,25 @@ void Foam::wedgePolyPatch::initTransforms() ); centreNormal_ /= mag(centreNormal_); - if - ( - mag(centreNormal_.x() + centreNormal_.y() + centreNormal_.z()) - < (1 - SMALL) - ) + const scalar cnCmptSum = + centreNormal_.x() + centreNormal_.y() + centreNormal_.z(); + + if (mag(cnCmptSum) < (1 - SMALL)) { - FatalErrorIn - ( - "wedgePolyPatch::wedgePolyPatch(const polyPatch&, " - "const fvBoundaryMesh&)" - ) << "wedge " << name() + FatalErrorIn("wedgePolyPatch::initTransforms()") + << "wedge " << name() << " centre plane does not align with a coordinate plane by " - << 1 - - mag(centreNormal_.x()+centreNormal_.y()+centreNormal_.z()) + << 1 - mag(cnCmptSum) << exit(FatalError); } axis_ = centreNormal_ ^ patchNormal_; scalar magAxis = mag(axis_); - axis_ /= magAxis; if (magAxis < SMALL) { - FatalErrorIn - ( - "wedgePolyPatch::initTransforms()" - ) << "wedge " << name() + FatalErrorIn("wedgePolyPatch::initTransforms()") + << "wedge " << name() << " plane aligns with a coordinate plane." << nl << " The wedge plane should make a small angle (~2.5deg)" " with the coordinate plane" << nl @@ -95,6 +87,8 @@ void Foam::wedgePolyPatch::initTransforms() << exit(FatalError); } + axis_ /= magAxis; + faceT_ = rotationTensor(centreNormal_, patchNormal_); cellT_ = faceT_ & faceT_; } diff --git a/src/OpenFOAM/primitives/quaternion/quaternion.C b/src/OpenFOAM/primitives/quaternion/quaternion.C index f7ff1904c9..cedbb56c32 100644 --- a/src/OpenFOAM/primitives/quaternion/quaternion.C +++ b/src/OpenFOAM/primitives/quaternion/quaternion.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -51,6 +51,31 @@ Foam::word Foam::name(const quaternion& q) } +Foam::quaternion Foam::slerp +( + const quaternion& qa, + const quaternion& qb, + const scalar t +) +{ + // Calculate angle between the quaternions + scalar cosHalfTheta = qa & qb; + + if (mag(cosHalfTheta) >= 1) + { + return qa; + } + + scalar halfTheta = acos(cosHalfTheta); + scalar sinHalfTheta = sqrt(1.0 - sqr(cosHalfTheta)); + + scalar wa = sin((1 - t)*halfTheta)/sinHalfTheta; + scalar wb = sin(t*halfTheta)/sinHalfTheta; + + return wa*qa + wb*qb; +} + + // * * * * * * * * * * * * * * * IOstream Operators * * * * * * * * * * * * // Foam::Istream& Foam::operator>>(Istream& is, quaternion& q) diff --git a/src/OpenFOAM/primitives/quaternion/quaternion.H b/src/OpenFOAM/primitives/quaternion/quaternion.H index 5a5e3fd697..dbe7ad15b5 100644 --- a/src/OpenFOAM/primitives/quaternion/quaternion.H +++ b/src/OpenFOAM/primitives/quaternion/quaternion.H @@ -106,6 +106,15 @@ public: // and angle theta inline quaternion(const vector& d, const scalar theta); + //- Construct a rotation quaternion given the direction d + // and cosine angle cosTheta and a if d is normalized + inline quaternion + ( + const vector& d, + const scalar cosTheta, + const bool normalized + ); + //- Construct given scalar part, the vector part = vector::zero inline explicit quaternion(const scalar w); @@ -140,6 +149,8 @@ public: //- The rotation tensor corresponding the quaternion inline tensor R() const; + inline quaternion normalized() const; + // Edit @@ -207,6 +218,14 @@ inline quaternion inv(const quaternion& q); //- Return a string representation of a quaternion word name(const quaternion&); +//- Spherical linear interpolation of quaternions +quaternion slerp +( + const quaternion& qa, + const quaternion& qb, + const scalar t +); + //- Data associated with quaternion type are contiguous template<> inline bool contiguous() {return true;} diff --git a/src/OpenFOAM/primitives/quaternion/quaternionI.H b/src/OpenFOAM/primitives/quaternion/quaternionI.H index 01c5b74bc2..1582f5c733 100644 --- a/src/OpenFOAM/primitives/quaternion/quaternionI.H +++ b/src/OpenFOAM/primitives/quaternion/quaternionI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -37,9 +37,27 @@ inline Foam::quaternion::quaternion(const scalar w, const vector& v) inline Foam::quaternion::quaternion(const vector& d, const scalar theta) : w_(cos(0.5*theta)), - v_((sin(0.5*theta)/magSqr(d))*d) + v_((sin(0.5*theta)/mag(d))*d) +{} + +inline Foam::quaternion::quaternion +( + const vector& d, + const scalar cosTheta, + const bool normalized +) { - normalize(); + scalar cosHalfTheta2 = 0.5*(cosTheta + 1); + w_ = sqrt(cosHalfTheta2); + + if (normalized) + { + v_ = sqrt(1 - cosHalfTheta2)*d; + } + else + { + v_ = (sqrt(1 - cosHalfTheta2)/mag(d))*d; + } } inline Foam::quaternion::quaternion(const scalar w) @@ -71,9 +89,10 @@ inline Foam::quaternion::quaternion const tensor& rotationTensor ) { - scalar trace = rotationTensor.xx() - + rotationTensor.yy() - + rotationTensor.zz(); + scalar trace = + rotationTensor.xx() + + rotationTensor.yy() + + rotationTensor.zz(); if (trace > 0) { @@ -168,6 +187,12 @@ inline Foam::vector& Foam::quaternion::v() } +inline Foam::quaternion Foam::quaternion::normalized() const +{ + return *this/mag(*this); +} + + inline void Foam::quaternion::normalize() { operator/=(mag(*this)); diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C index 348f61db2b..3a9773ce23 100644 --- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C +++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -1604,6 +1604,8 @@ Foam::autoPtr Foam::polyMeshAdder::add polyBoundaryMesh& allPatches = const_cast(mesh0.boundaryMesh()); allPatches.setSize(allPatchNames.size()); + labelList patchSizes(allPatches.size()); + labelList patchStarts(allPatches.size()); label startFaceI = nInternalFaces; @@ -1629,7 +1631,9 @@ Foam::autoPtr Foam::polyMeshAdder::add } else { - // Clone. + // Clone. Note dummy size and start. Gets overwritten later in + // resetPrimitives. This avoids getting temporarily illegal + // SubList construction in polyPatch. allPatches.set ( allPatchI, @@ -1637,10 +1641,12 @@ Foam::autoPtr Foam::polyMeshAdder::add ( allPatches, allPatchI, - nFaces[patch0], - startFaceI + 0, // dummy size + 0 // dummy start ) ); + patchSizes[allPatchI] = nFaces[patch0]; + patchStarts[allPatchI] = startFaceI; // Record new index in allPatches from0ToAllPatches[patch0] = allPatchI; @@ -1686,10 +1692,12 @@ Foam::autoPtr Foam::polyMeshAdder::add ( allPatches, allPatchI, - nFaces[uncompactAllPatchI], - startFaceI + 0, // dummy size + 0 // dummy start ) ); + patchSizes[allPatchI] = nFaces[uncompactAllPatchI]; + patchStarts[allPatchI] = startFaceI; // Record new index in allPatches from1ToAllPatches[patch1] = allPatchI; @@ -1702,6 +1710,8 @@ Foam::autoPtr Foam::polyMeshAdder::add } allPatches.setSize(allPatchI); + patchSizes.setSize(allPatchI); + patchStarts.setSize(allPatchI); // Construct map information before changing mesh0 primitives @@ -1740,9 +1750,6 @@ Foam::autoPtr Foam::polyMeshAdder::add // Now we have extracted all information from all meshes. // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - labelList patchSizes(getPatchSizes(allPatches)); - labelList patchStarts(getPatchStarts(allPatches)); - mesh0.resetMotion(); // delete any oldPoints. mesh0.resetPrimitives ( diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C index 475430bd83..97b84a6049 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporation/LiquidEvaporation.C @@ -216,7 +216,7 @@ Foam::scalar Foam::LiquidEvaporation::dh } case (parent::etEnthalpyDifference): { - scalar hc = this->owner().composition().carrier().Hs(idc, p, T); + scalar hc = this->owner().composition().carrier().Ha(idc, p, T); scalar hp = liquids_.properties()[idl].h(p, T); dh = hc - hp; diff --git a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C index 32821f8966..803b849fa4 100644 --- a/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C +++ b/src/lagrangian/intermediate/submodels/Reacting/PhaseChangeModel/LiquidEvaporationBoil/LiquidEvaporationBoil.C @@ -171,8 +171,8 @@ void Foam::LiquidEvaporationBoil::calculate forAll(this->owner().thermo().carrier().Y(), i) { scalar Yc = this->owner().thermo().carrier().Y()[i][cellI]; - Hc += Yc*this->owner().thermo().carrier().Hs(i, pc, Tc); - Hsc += Yc*this->owner().thermo().carrier().Hs(i, ps, Ts); + Hc += Yc*this->owner().thermo().carrier().Ha(i, pc, Tc); + Hsc += Yc*this->owner().thermo().carrier().Ha(i, ps, Ts); Cpc += Yc*this->owner().thermo().carrier().Cp(i, ps, Ts); kappac += Yc*this->owner().thermo().carrier().kappa(i, ps, Ts); } @@ -315,7 +315,7 @@ Foam::scalar Foam::LiquidEvaporationBoil::dh } case (parent::etEnthalpyDifference): { - scalar hc = this->owner().composition().carrier().Hs(idc, p, TDash); + scalar hc = this->owner().composition().carrier().Ha(idc, p, TDash); scalar hp = liquids_.properties()[idl].h(p, TDash); dh = hc - hp; diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 4abc57669e..2621adbc23 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -99,6 +99,7 @@ $(cellSources)/sphereToCell/sphereToCell.C $(cellSources)/cylinderToCell/cylinderToCell.C $(cellSources)/faceZoneToCell/faceZoneToCell.C $(cellSources)/cylinderAnnulusToCell/cylinderAnnulusToCell.C +$(cellSources)/targetVolumeToCell/targetVolumeToCell.C faceSources = sets/faceSources $(faceSources)/faceToFace/faceToFace.C diff --git a/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C new file mode 100644 index 0000000000..c64c29e1e6 --- /dev/null +++ b/src/meshTools/sets/cellSources/targetVolumeToCell/targetVolumeToCell.C @@ -0,0 +1,328 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation + \\/ 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 "targetVolumeToCell.H" +#include "polyMesh.H" +#include "globalMeshData.H" +#include "plane.H" + +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + +defineTypeNameAndDebug(targetVolumeToCell, 0); + +addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, word); + +addToRunTimeSelectionTable(topoSetSource, targetVolumeToCell, istream); + +} + + +Foam::topoSetSource::addToUsageTable Foam::targetVolumeToCell::usage_ +( + targetVolumeToCell::typeName, + "\n Usage: targetVolumeToCell (nx ny nz)\n\n" + " Adjust plane until obtained selected volume\n\n" +); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::scalar Foam::targetVolumeToCell::volumeOfSet +( + const PackedBoolList& selected +) const +{ + scalar sumVol = 0.0; + forAll(selected, cellI) + { + if (selected[cellI]) + { + sumVol += mesh_.cellVolumes()[cellI]; + } + } + return returnReduce(sumVol, sumOp()); +} + + +Foam::label Foam::targetVolumeToCell::selectCells +( + const scalar normalComp, + PackedBoolList& selected +) const +{ + selected.setSize(mesh_.nCells()); + selected = false; + + label nSelected = 0; + + forAll(mesh_.cellCentres(), cellI) + { + const point& cc = mesh_.cellCentres()[cellI]; + + if ((cc&n_) < normalComp) + { + selected[cellI] = true; + nSelected++; + } + } + return returnReduce(nSelected, sumOp