diff --git a/applications/solvers/multiphase/interFoam/Allwclean b/applications/solvers/multiphase/interFoam/Allwclean index 350d4b268b070f686f65a6e959f2f5abdfde608d..2c54bfa3218478eef96fce9e9bbc8f23612add83 100755 --- a/applications/solvers/multiphase/interFoam/Allwclean +++ b/applications/solvers/multiphase/interFoam/Allwclean @@ -6,5 +6,6 @@ wclean wclean interDyMFoam wclean MRFInterFoam wclean porousInterFoam +wclean LTSInterFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interFoam/Allwmake b/applications/solvers/multiphase/interFoam/Allwmake index b25e4440e60bc084516b2829565268642a7ea2cc..8044426582eacdd44b4277674eb8d715c0365ed0 100755 --- a/applications/solvers/multiphase/interFoam/Allwmake +++ b/applications/solvers/multiphase/interFoam/Allwmake @@ -6,5 +6,6 @@ wmake wmake interDyMFoam wmake MRFInterFoam wmake porousInterFoam +wmake LTSInterFoam # ----------------------------------------------------------------- end-of-file diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C new file mode 100644 index 0000000000000000000000000000000000000000..9fcb19a7cc242f590087ea7b62ad268606516e43 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/LTSInterFoam.C @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +Application + interFoam + +Description + Solver for 2 incompressible, isothermal immiscible fluids using a VOF + (volume of fluid) phase-fraction based interface capturing approach. + + The momentum and other fluid properties are of the "mixture" and a single + momentum equation is solved. + + Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected. + + For a two-fluid approach see twoPhaseEulerFoam. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "MULES.H" +#include "subCycle.H" +#include "interfaceProperties.H" +#include "twoPhaseMixture.H" +#include "turbulenceModel.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readPISOControls.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "correctPhi.H" + #include "CourantNo.H" + #include "setInitialrDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readPISOControls.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + #include "setrDeltaT.H" + + twoPhaseProperties.correct(); + + #include "alphaEqnSubCycle.H" + turbulence->correct(); + #include "UEqn.H" + + // --- PISO loop + for (int corr=0; corr<nCorr; corr++) + { + #include "pEqn.H" + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C new file mode 100644 index 0000000000000000000000000000000000000000..e281ae575cbfd41e0000cb1f2c1d178417047c63 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.C @@ -0,0 +1,83 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "MULES.H" +#include "upwind.H" +#include "uncorrectedSnGrad.H" +#include "gaussConvectionScheme.H" +#include "gaussLaplacianScheme.H" +#include "uncorrectedSnGrad.H" +#include "surfaceInterpolate.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "syncTools.H" + +#include "fvCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void Foam::MULES::explicitLTSSolve +( + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +) +{ + explicitLTSSolve + ( + geometricOneField(), + psi, + phi, + phiPsi, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +void Foam::MULES::implicitSolve +( + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +) +{ + implicitSolve + ( + geometricOneField(), + psi, + phi, + phiPsi, + zeroField(), zeroField(), + psiMax, psiMin + ); +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H new file mode 100644 index 0000000000000000000000000000000000000000..729b041aa62d0b7870d28a613ab3f1a332d96bae --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULES.H @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +Global + MULES + +Description + MULES: Multidimensional universal limiter with explicit solution. + + Solve a convective-only transport equation using an explicit universal + multi-dimensional limiter. + + Parameters are the variable to solve, the normal convective flux and the + actual explicit flux of the variable which is also used to return limited + flux used in the bounded-solution. + +SourceFiles + MULES.C + +\*---------------------------------------------------------------------------*/ + +#ifndef MULES_H +#define MULES_H + +#include "volFields.H" +#include "surfaceFieldsFwd.H" +#include "primitiveFieldsFwd.H" +#include "zeroField.H" +#include "geometricOneField.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace MULES +{ + +template<class RhoType, class SpType, class SuType> +void explicitLTSSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void explicitLTSSolve +( + volScalarField& psi, + const surfaceScalarField& phiBD, + surfaceScalarField& phiPsi, + const scalar psiMax, + const scalar psiMin +); + +template<class RhoType, class SpType, class SuType> +void implicitSolve +( + const RhoType& rho, + volScalarField& gamma, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +); + +void implicitSolve +( + volScalarField& gamma, + const surfaceScalarField& phi, + surfaceScalarField& phiCorr, + const scalar psiMax, + const scalar psiMin +); + +template<class RhoType, class SpType, class SuType> +void limiter +( + scalarField& allLambda, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiBD, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +); + +} // End namespace MULES + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "MULESTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..a6d3a3f6b246e2463979604d84ad33dd108c5c01 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/MULESTemplates.C @@ -0,0 +1,602 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "MULES.H" +#include "upwind.H" +#include "uncorrectedSnGrad.H" +#include "gaussConvectionScheme.H" +#include "gaussLaplacianScheme.H" +#include "uncorrectedSnGrad.H" +#include "surfaceInterpolate.H" +#include "fvcSurfaceIntegrate.H" +#include "slicedSurfaceFields.H" +#include "syncTools.H" + +#include "fvCFD.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::explicitLTSSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + Info<< "MULES: Solving for " << psi.name() << endl; + + const fvMesh& mesh = psi.mesh(); + psi.correctBoundaryConditions(); + + surfaceScalarField phiBD = upwind<scalar>(psi.mesh(), phi).flux(psi); + + surfaceScalarField& phiCorr = phiPsi; + phiCorr -= phiBD; + + scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + limiter + ( + allLambda, + rho, + psi, + phiBD, + phiCorr, + Sp.field(), + Su.field(), + psiMax, + psiMin, + 3 + ); + + phiPsi = phiBD + lambda*phiCorr; + + scalarField& psiIf = psi; + const scalarField& psi0 = psi.oldTime(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT"); + + psiIf = 0.0; + fvc::surfaceIntegrate(psiIf, phiPsi); + + if (mesh.moving()) + { + psiIf = + ( + mesh.Vsc0()*rho.oldTime()*psi0*rDeltaT/mesh.Vsc() + + Su.field() + - psiIf + )/(rho*rDeltaT - Sp.field()); + } + else + { + psiIf = + ( + rho.oldTime()*psi0*rDeltaT + + Su.field() + - psiIf + )/(rho*rDeltaT - Sp.field()); + } + + psi.correctBoundaryConditions(); +} + + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::implicitSolve +( + const RhoType& rho, + volScalarField& psi, + const surfaceScalarField& phi, + surfaceScalarField& phiPsi, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin +) +{ + const fvMesh& mesh = psi.mesh(); + + const dictionary& MULEScontrols = mesh.solverDict(psi.name()); + + label maxIter + ( + readLabel(MULEScontrols.lookup("maxIter")) + ); + + label nLimiterIter + ( + readLabel(MULEScontrols.lookup("nLimiterIter")) + ); + + scalar maxUnboundedness + ( + readScalar(MULEScontrols.lookup("maxUnboundedness")) + ); + + scalar CoCoeff + ( + readScalar(MULEScontrols.lookup("CoCoeff")) + ); + + scalarField allCoLambda(mesh.nFaces()); + + { + surfaceScalarField Cof = + mesh.time().deltaT()*mesh.surfaceInterpolation::deltaCoeffs() + *mag(phi)/mesh.magSf(); + + slicedSurfaceScalarField CoLambda + ( + IOobject + ( + "CoLambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allCoLambda, + false // Use slices for the couples + ); + + CoLambda == 1.0/max(CoCoeff*Cof, scalar(1)); + } + + scalarField allLambda(allCoLambda); + //scalarField allLambda(mesh.nFaces(), 1.0); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + linear<scalar> CDs(mesh); + upwind<scalar> UDs(mesh, phi); + //fv::uncorrectedSnGrad<scalar> snGrads(mesh); + + fvScalarMatrix psiConvectionDiffusion + ( + fvm::ddt(rho, psi) + + fv::gaussConvectionScheme<scalar>(mesh, phi, UDs).fvmDiv(phi, psi) + //- fv::gaussLaplacianScheme<scalar, scalar>(mesh, CDs, snGrads) + //.fvmLaplacian(Dpsif, psi) + - fvm::Sp(Sp, psi) + - Su + ); + + surfaceScalarField phiBD = psiConvectionDiffusion.flux(); + + surfaceScalarField& phiCorr = phiPsi; + phiCorr -= phiBD; + + for (label i=0; i<maxIter; i++) + { + if (i != 0 && i < 4) + { + allLambda = allCoLambda; + } + + limiter + ( + allLambda, + rho, + psi, + phiBD, + phiCorr, + Sp.field(), + Su.field(), + psiMax, + psiMin, + nLimiterIter + ); + + solve + ( + psiConvectionDiffusion + fvc::div(lambda*phiCorr), + MULEScontrols + ); + + scalar maxPsiM1 = gMax(psi.internalField()) - 1.0; + scalar minPsi = gMin(psi.internalField()); + + scalar unboundedness = max(max(maxPsiM1, 0.0), -min(minPsi, 0.0)); + + if (unboundedness < maxUnboundedness) + { + break; + } + else + { + Info<< "MULES: max(" << psi.name() << " - 1) = " << maxPsiM1 + << " min(" << psi.name() << ") = " << minPsi << endl; + + phiBD = psiConvectionDiffusion.flux(); + + /* + word gammaScheme("div(phi,gamma)"); + word gammarScheme("div(phirb,gamma)"); + + const surfaceScalarField& phir = + mesh.lookupObject<surfaceScalarField>("phir"); + + phiCorr = + fvc::flux + ( + phi, + psi, + gammaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - psi, gammarScheme), + psi, + gammarScheme + ) + - phiBD; + */ + } + } + + phiPsi = psiConvectionDiffusion.flux() + lambda*phiCorr; +} + + +template<class RhoType, class SpType, class SuType> +void Foam::MULES::limiter +( + scalarField& allLambda, + const RhoType& rho, + const volScalarField& psi, + const surfaceScalarField& phiBD, + const surfaceScalarField& phiCorr, + const SpType& Sp, + const SuType& Su, + const scalar psiMax, + const scalar psiMin, + const label nLimiterIter +) +{ + const scalarField& psiIf = psi; + const volScalarField::GeometricBoundaryField& psiBf = psi.boundaryField(); + + const scalarField& psi0 = psi.oldTime(); + + const fvMesh& mesh = psi.mesh(); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighb = mesh.neighbour(); + tmp<volScalarField::DimensionedInternalField> tVsc = mesh.Vsc(); + const scalarField& V = tVsc(); + + const volScalarField& rDeltaT = + mesh.objectRegistry::lookupObject<volScalarField>("rSubDeltaT"); + + const scalarField& phiBDIf = phiBD; + const surfaceScalarField::GeometricBoundaryField& phiBDBf = + phiBD.boundaryField(); + + const scalarField& phiCorrIf = phiCorr; + const surfaceScalarField::GeometricBoundaryField& phiCorrBf = + phiCorr.boundaryField(); + + slicedSurfaceScalarField lambda + ( + IOobject + ( + "lambda", + mesh.time().timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimless, + allLambda, + false // Use slices for the couples + ); + + scalarField& lambdaIf = lambda; + surfaceScalarField::GeometricBoundaryField& lambdaBf = + lambda.boundaryField(); + + scalarField psiMaxn(psiIf.size(), psiMin); + scalarField psiMinn(psiIf.size(), psiMax); + + scalarField sumPhiBD(psiIf.size(), 0.0); + + scalarField sumPhip(psiIf.size(), VSMALL); + scalarField mSumPhim(psiIf.size(), VSMALL); + + forAll(phiCorrIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + psiMaxn[own] = max(psiMaxn[own], psiIf[nei]); + psiMinn[own] = min(psiMinn[own], psiIf[nei]); + + psiMaxn[nei] = max(psiMaxn[nei], psiIf[own]); + psiMinn[nei] = min(psiMinn[nei], psiIf[own]); + + sumPhiBD[own] += phiBDIf[facei]; + sumPhiBD[nei] -= phiBDIf[facei]; + + scalar phiCorrf = phiCorrIf[facei]; + + if (phiCorrf > 0.0) + { + sumPhip[own] += phiCorrf; + mSumPhim[nei] += phiCorrf; + } + else + { + mSumPhim[own] -= phiCorrf; + sumPhip[nei] -= phiCorrf; + } + } + + forAll(phiCorrBf, patchi) + { + const fvPatchScalarField& psiPf = psiBf[patchi]; + const scalarField& phiBDPf = phiBDBf[patchi]; + const scalarField& phiCorrPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + if (psiPf.coupled()) + { + scalarField psiPNf = psiPf.patchNeighbourField(); + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPNf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPNf[pFacei]); + } + } + else + { + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + psiMaxn[pfCelli] = max(psiMaxn[pfCelli], psiPf[pFacei]); + psiMinn[pfCelli] = min(psiMinn[pfCelli], psiPf[pFacei]); + } + } + + forAll(phiCorrPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + sumPhiBD[pfCelli] += phiBDPf[pFacei]; + + scalar phiCorrf = phiCorrPf[pFacei]; + + if (phiCorrf > 0.0) + { + sumPhip[pfCelli] += phiCorrf; + } + else + { + mSumPhim[pfCelli] -= phiCorrf; + } + } + } + + psiMaxn = min(psiMaxn, psiMax); + psiMinn = max(psiMinn, psiMin); + + //scalar smooth = 0.5; + //psiMaxn = min((1.0 - smooth)*psiIf + smooth*psiMaxn, psiMax); + //psiMinn = max((1.0 - smooth)*psiIf + smooth*psiMinn, psiMin); + + if (mesh.moving()) + { + tmp<volScalarField::DimensionedInternalField> V0 = mesh.Vsc0(); + + psiMaxn = + V*((rho*rDeltaT - Sp)*psiMaxn - Su) + - (V0()*rDeltaT)*rho.oldTime()*psi0 + + sumPhiBD; + + psiMinn = + V*(Su - (rho*rDeltaT - Sp)*psiMinn) + + (V0*rDeltaT)*rho.oldTime()*psi0 + - sumPhiBD; + } + else + { + psiMaxn = + V*((rho*rDeltaT - Sp)*psiMaxn - (rho.oldTime()*rDeltaT)*psi0 - Su) + + sumPhiBD; + + psiMinn = + V*((rho*rDeltaT)*psi0 - (rho.oldTime()*rDeltaT - Sp)*psiMinn + Su) + - sumPhiBD; + } + + scalarField sumlPhip(psiIf.size()); + scalarField mSumlPhim(psiIf.size()); + + for(int j=0; j<nLimiterIter; j++) + { + sumlPhip = 0.0; + mSumlPhim = 0.0; + + forAll(lambdaIf, facei) + { + label own = owner[facei]; + label nei = neighb[facei]; + + scalar lambdaPhiCorrf = lambdaIf[facei]*phiCorrIf[facei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[own] += lambdaPhiCorrf; + mSumlPhim[nei] += lambdaPhiCorrf; + } + else + { + mSumlPhim[own] -= lambdaPhiCorrf; + sumlPhip[nei] -= lambdaPhiCorrf; + } + } + + forAll(lambdaBf, patchi) + { + scalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + scalar lambdaPhiCorrf = lambdaPf[pFacei]*phiCorrfPf[pFacei]; + + if (lambdaPhiCorrf > 0.0) + { + sumlPhip[pfCelli] += lambdaPhiCorrf; + } + else + { + mSumlPhim[pfCelli] -= lambdaPhiCorrf; + } + } + } + + forAll (sumlPhip, celli) + { + sumlPhip[celli] = + max(min + ( + (sumlPhip[celli] + psiMaxn[celli])/mSumPhim[celli], + 1.0), 0.0 + ); + + mSumlPhim[celli] = + max(min + ( + (mSumlPhim[celli] + psiMinn[celli])/sumPhip[celli], + 1.0), 0.0 + ); + } + + const scalarField& lambdam = sumlPhip; + const scalarField& lambdap = mSumlPhim; + + forAll(lambdaIf, facei) + { + if (phiCorrIf[facei] > 0.0) + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdap[owner[facei]], lambdam[neighb[facei]]) + ); + } + else + { + lambdaIf[facei] = min + ( + lambdaIf[facei], + min(lambdam[owner[facei]], lambdap[neighb[facei]]) + ); + } + } + + + forAll(lambdaBf, patchi) + { + fvsPatchScalarField& lambdaPf = lambdaBf[patchi]; + const scalarField& phiCorrfPf = phiCorrBf[patchi]; + + const labelList& pFaceCells = mesh.boundary()[patchi].faceCells(); + + forAll(lambdaPf, pFacei) + { + label pfCelli = pFaceCells[pFacei]; + + if (phiCorrfPf[pFacei] > 0.0) + { + lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdap[pfCelli]); + } + else + { + lambdaPf[pFacei] = min(lambdaPf[pFacei], lambdam[pfCelli]); + } + } + } + + syncTools::syncFaceList(mesh, allLambda, minEqOp<scalar>()); + } +} + + +// ************************************************************************* // diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files new file mode 100644 index 0000000000000000000000000000000000000000..205812e978dd9dfe7ace29f84753167d1186ab39 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/files @@ -0,0 +1,4 @@ +LTSInterFoam.C +MULES.C + +EXE = $(FOAM_APPBIN)/LTSInterFoam diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options new file mode 100644 index 0000000000000000000000000000000000000000..24349f694e016e85732d7896b7b98c57083eaa66 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/Make/options @@ -0,0 +1,15 @@ +EXE_INC = \ + -I.. \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ + -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -ltwoPhaseInterfaceProperties \ + -lincompressibleTransportModels \ + -lincompressibleTurbulenceModel \ + -lincompressibleRASModels \ + -lincompressibleLESModels \ + -lfiniteVolume diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H new file mode 100644 index 0000000000000000000000000000000000000000..a4cb6e7688c168c92d3ac1aebd5cfb1cddf93382 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqn.H @@ -0,0 +1,36 @@ +{ + word alphaScheme("div(phi,alpha)"); + word alpharScheme("div(phirb,alpha)"); + + surfaceScalarField phic = mag(phi/mesh.magSf()); + phic = min(interface.cAlpha()*phic, max(phic)); + surfaceScalarField phir = phic*interface.nHatf(); + + for (int aCorr=0; aCorr<nAlphaCorr; aCorr++) + { + surfaceScalarField phiAlpha = + fvc::flux + ( + phi, + alpha1, + alphaScheme + ) + + fvc::flux + ( + -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme), + alpha1, + alpharScheme + ); + + MULES::explicitLTSSolve(alpha1, phi, phiAlpha, 1, 0); + //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0); + + rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2; + } + + Info<< "Liquid phase volume fraction = " + << alpha1.weightedAverage(mesh.V()).value() + << " Min(alpha1) = " << min(alpha1).value() + << " Max(alpha1) = " << max(alpha1).value() + << endl; +} diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H new file mode 100644 index 0000000000000000000000000000000000000000..9aae37a8bef50329d7e2b4cac649219d48dcc707 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H @@ -0,0 +1,35 @@ +label nAlphaCorr +( + readLabel(piso.lookup("nAlphaCorr")) +); + +label nAlphaSubCycles +( + readLabel(piso.lookup("nAlphaSubCycles")) +); + +if (nAlphaSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + surfaceScalarField rhoPhiSum = 0.0*rhoPhi; + + for + ( + subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles); + !(++alphaSubCycle).end(); + ) + { +# include "alphaEqn.H" + rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; + } + + rhoPhi = rhoPhiSum; +} +else +{ +# include "alphaEqn.H" +} + +interface.correct(); + +rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2; diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H new file mode 100644 index 0000000000000000000000000000000000000000..1f2f082fb36b8dc2f4538b2714990c1912482107 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setInitialrDeltaT.H @@ -0,0 +1,30 @@ +scalar maxDeltaT +( + piso.lookupOrDefault<scalar>("maxDeltaT", GREAT) +); + +volScalarField rDeltaT +( + IOobject + ( + "rDeltaT", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT) +); + +volScalarField rSubDeltaT +( + IOobject + ( + "rSubDeltaT", + runTime.timeName(), + mesh + ), + mesh, + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT) +); diff --git a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H new file mode 100644 index 0000000000000000000000000000000000000000..f5b8450ddd9d1032784bf7f96eed44932fc54e44 --- /dev/null +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -0,0 +1,109 @@ +{ + scalar maxCo + ( + piso.lookupOrDefault<scalar>("maxCo", 0.9) + ); + + scalar maxAlphaCo + ( + piso.lookupOrDefault<scalar>("maxAlphaCo", 0.2) + ); + + scalar rDeltaTSmoothingCoeff + ( + piso.lookupOrDefault<scalar>("rDeltaTSmoothingCoeff", 0.1) + ); + + label nAlphaSpreadIter + ( + piso.lookupOrDefault<label>("nAlphaSpreadIter", 1) + ); + + label nAlphaSweepIter + ( + piso.lookupOrDefault<label>("nAlphaSweepIter", 5) + ); + + scalar rDeltaTDampingCoeff + ( + piso.lookupOrDefault<scalar>("rDeltaTDampingCoeff", 1.0) + ); + + scalar maxDeltaT + ( + piso.lookupOrDefault<scalar>("maxDeltaT", GREAT) + ); + + volScalarField rDeltaT0 = rDeltaT; + + // Set the reciprocal time-step using an effective maximum Courant number + rDeltaT = max + ( + 1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT), + fvc::surfaceSum + ( + mag(rhoPhi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf()) + )/rho + ); + + // Limit the time-step further in the region of the interface + { + surfaceScalarField alphaf = fvc::interpolate(alpha1); + + surfaceScalarField SfUfbyDelta = + pos(alphaf - 0.01)*pos(0.99 - alphaf) + *mesh.surfaceInterpolation::deltaCoeffs()*mag(phi); + + rDeltaT = max + ( + rDeltaT, + fvc::surfaceSum(mag(SfUfbyDelta/(maxAlphaCo*mesh.magSf()))) + ); + } + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + if (rDeltaTSmoothingCoeff < 1.0) + { + fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); + } + + if (nAlphaSpreadIter > 0) + { + fvc::spread(rDeltaT, alpha1, nAlphaSpreadIter); + } + + if (nAlphaSweepIter > 0) + { + fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter); + } + + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + // Limit rate of change of time scale + // - reduce as much as required + // - only increase at a fraction of old time scale + if + ( + rDeltaTDampingCoeff < 1.0 + && runTime.timeIndex() > runTime.startTimeIndex() + 1 + ) + { + Info<< "Damping rDeltaT" << endl; + rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff); + } + + Info<< "Flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + + label nAlphaSubCycles + ( + readLabel(piso.lookup("nAlphaSubCycles")) + ); + + rSubDeltaT = rDeltaT*nAlphaSubCycles; +} diff --git a/applications/utilities/mesh/advanced/selectCells/selectCells.C b/applications/utilities/mesh/advanced/selectCells/selectCells.C index a298c4207df3792c4b09e34279d8de68a495f7c2..827449c02c91e24e5146c6d0e51d0e0b906dc406 100644 --- a/applications/utilities/mesh/advanced/selectCells/selectCells.C +++ b/applications/utilities/mesh/advanced/selectCells/selectCells.C @@ -55,6 +55,7 @@ Description #include "edgeStats.H" #include "treeDataTriSurface.H" #include "indexedOctree.H" +#include "globalMeshData.H" using namespace Foam; @@ -298,7 +299,7 @@ label selectOutsideCells mesh, outsideFaces.shrink(), outsideFacesInfo.shrink(), - mesh.nCells() // max iterations + mesh.globalData().nTotalCells()+1 // max iterations ); // Now regionCalc should hold info on cells that are reachable from diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L index 5fd94ef38dde348c9de1d82a08e8ccbb02e13303..6f0f7c976dff5c8215346adf0c726761612fabc1 100644 --- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L +++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L @@ -43,7 +43,7 @@ Description #include "polyMeshZipUpCells.H" #include "wallPolyPatch.H" #include "symmetryPolyPatch.H" -#include "cyclicPolyPatch.H" +#include "oldCyclicPolyPatch.H" #include "Swap.H" #include "IFstream.H" #include "readHexLabel.H" @@ -900,7 +900,7 @@ int main(int argc, char *argv[]) fluentToFoamType.insert("interface", polyPatch::typeName); fluentToFoamType.insert("internal", polyPatch::typeName); fluentToFoamType.insert("solid", polyPatch::typeName); - fluentToFoamType.insert("fan", cyclicPolyPatch::typeName); + fluentToFoamType.insert("fan", oldCyclicPolyPatch::typeName); fluentToFoamType.insert("radiator", polyPatch::typeName); fluentToFoamType.insert("porous-jump", polyPatch::typeName); diff --git a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C index ab764547a3582a9900246754f22ebf967da1a3a5..066fab386f3757c17abc4716ce4a8b0405a9225a 100644 --- a/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C +++ b/applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C @@ -41,7 +41,7 @@ Description #include "wallPolyPatch.H" #include "symmetryPolyPatch.H" #include "wedgePolyPatch.H" -#include "cyclicPolyPatch.H" +#include "oldCyclicPolyPatch.H" #include "unitConversion.H" using namespace Foam; diff --git a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H index 04b36020ea210bf0e3a5d88e8ea8e8297d2d36e2..1a559559e92509c4273e482d5ce937f9cfd4c083 100644 --- a/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H +++ b/applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H @@ -196,7 +196,7 @@ const word* kivaPatchTypes[nBCs] = &polyPatch::typeName, &polyPatch::typeName, &symmetryPolyPatch::typeName, - &cyclicPolyPatch::typeName + &oldCyclicPolyPatch::typeName }; enum patchTypeNames diff --git a/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C b/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C index 53dffd39c6bca64173bc5b5d6f7881e02d9e72ea..d01a114111ab7707bfe4c2569db8bc32fcffe910 100644 --- a/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C +++ b/applications/utilities/mesh/conversion/sammToFoam/readBoundary.C @@ -29,7 +29,7 @@ Description #include "sammMesh.H" #include "Time.H" #include "wallPolyPatch.H" -#include "cyclicPolyPatch.H" +#include "oldCyclicPolyPatch.H" #include "symmetryPolyPatch.H" #include "preservePatchTypes.H" #include "IFstream.H" @@ -208,7 +208,7 @@ void sammMesh::readBoundary() { // incorrect. should be cyclicPatch but this // requires info on connected faces. - patchTypes_[patchLabel] = cyclicPolyPatch::typeName; + patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName; } else { diff --git a/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C b/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C index a96ee6e96bb71a7df1ae627caf201cce00339f51..75a18600062522c62dd70f4bc143cea90f9a52c9 100644 --- a/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C +++ b/applications/utilities/mesh/conversion/star3ToFoam/readBoundary.C @@ -29,7 +29,7 @@ Description #include "starMesh.H" #include "Time.H" #include "wallPolyPatch.H" -#include "cyclicPolyPatch.H" +#include "oldCyclicPolyPatch.H" #include "symmetryPolyPatch.H" #include "preservePatchTypes.H" #include "IFstream.H" @@ -206,7 +206,7 @@ void starMesh::readBoundary() { // incorrect. should be cyclicPatch but this // requires info on connected faces. - patchTypes_[patchLabel] = cyclicPolyPatch::typeName; + patchTypes_[patchLabel] = oldCyclicPolyPatch::typeName; } else { diff --git a/applications/utilities/mesh/manipulation/subsetMesh/Make/options b/applications/utilities/mesh/manipulation/subsetMesh/Make/options index d27c95d033dd5d7b1995c8ff8dc406e35ca1f586..969020c4afaf5d784299462b9e1af282040ba6b4 100644 --- a/applications/utilities/mesh/manipulation/subsetMesh/Make/options +++ b/applications/utilities/mesh/manipulation/subsetMesh/Make/options @@ -4,4 +4,5 @@ EXE_INC = \ EXE_LIBS = \ -lfiniteVolume \ - -lmeshTools + -lmeshTools \ + -lgenericPatchFields diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C index 7ffa7a4ca56eadb0014219828e4b46ed74427453..54342097110b9fa455d5f449bf77556b27d7b6ba 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C +++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightField.C @@ -575,99 +575,99 @@ void ensightPointField } - label ensightPatchI = eMesh.patchPartOffset(); - - forAll(allPatchNames, patchi) - { - const word& patchName = allPatchNames[patchi]; - - eMesh.barrier(); - - if (patchNames.empty() || patchNames.found(patchName)) - { - const fvPatch& p = mesh.boundary()[patchi]; - if - ( - returnReduce(p.size(), sumOp<label>()) - > 0 - ) - { - // Renumber the patch points/faces into unique points - labelList pointToGlobal; - labelList uniqueMeshPointLabels; - autoPtr<globalIndex> globalPointsPtr = - mesh.globalData().mergePoints - ( - p.patch().meshPoints(), - p.patch().meshPointMap(), - pointToGlobal, - uniqueMeshPointLabels - ); - - if (Pstream::master()) - { - ensightFile.writePartHeader(ensightPatchI); - } - - writeField - ( - "coordinates", - Field<Type>(pf.internalField(), uniqueMeshPointLabels), - ensightFile - ); - - ensightPatchI++; - } - } - } - - // write faceZones, if requested - if (faceZoneNames.size()) - { - forAllConstIter(wordHashSet, faceZoneNames, iter) - { - const word& faceZoneName = iter.key(); + label ensightPatchI = eMesh.patchPartOffset(); - eMesh.barrier(); + forAll(allPatchNames, patchi) + { + const word& patchName = allPatchNames[patchi]; - label zoneID = mesh.faceZones().findZoneID(faceZoneName); + eMesh.barrier(); - const faceZone& fz = mesh.faceZones()[zoneID]; + if (patchNames.empty() || patchNames.found(patchName)) + { + const fvPatch& p = mesh.boundary()[patchi]; + if + ( + returnReduce(p.size(), sumOp<label>()) + > 0 + ) + { + // Renumber the patch points/faces into unique points + labelList pointToGlobal; + labelList uniqueMeshPointLabels; + autoPtr<globalIndex> globalPointsPtr = + mesh.globalData().mergePoints + ( + p.patch().meshPoints(), + p.patch().meshPointMap(), + pointToGlobal, + uniqueMeshPointLabels + ); - if (returnReduce(fz().nPoints(), sumOp<label>()) > 0) - { - // Renumber the faceZone points/faces into unique points - labelList pointToGlobal; - labelList uniqueMeshPointLabels; - autoPtr<globalIndex> globalPointsPtr = - mesh.globalData().mergePoints - ( - fz().meshPoints(), - fz().meshPointMap(), - pointToGlobal, - uniqueMeshPointLabels - ); + if (Pstream::master()) + { + ensightFile.writePartHeader(ensightPatchI); + } - if (Pstream::master()) - { - ensightFile.writePartHeader(ensightPatchI); - } + writeField + ( + "coordinates", + Field<Type>(pf.internalField(), uniqueMeshPointLabels), + ensightFile + ); - writeField - ( - "coordinates", - Field<Type> - ( - pf.internalField(), - uniqueMeshPointLabels - ), - ensightFile - ); + ensightPatchI++; + } + } + } - ensightPatchI++; - } - } - } + // write faceZones, if requested + if (faceZoneNames.size()) + { + forAllConstIter(wordHashSet, faceZoneNames, iter) + { + const word& faceZoneName = iter.key(); + + eMesh.barrier(); + + label zoneID = mesh.faceZones().findZoneID(faceZoneName); + + const faceZone& fz = mesh.faceZones()[zoneID]; + + if (returnReduce(fz().nPoints(), sumOp<label>()) > 0) + { + // Renumber the faceZone points/faces into unique points + labelList pointToGlobal; + labelList uniqueMeshPointLabels; + autoPtr<globalIndex> globalPointsPtr = + mesh.globalData().mergePoints + ( + fz().meshPoints(), + fz().meshPointMap(), + pointToGlobal, + uniqueMeshPointLabels + ); + + if (Pstream::master()) + { + ensightFile.writePartHeader(ensightPatchI); + } + + writeField + ( + "coordinates", + Field<Type> + ( + pf.internalField(), + uniqueMeshPointLabels + ), + ensightFile + ); + + ensightPatchI++; + } + } + } if (Pstream::master()) { diff --git a/bin/tools/foamConfigurePaths b/bin/tools/foamConfigurePaths old mode 100755 new mode 100644 index 823c816d51c4bbc1eb427112239e56383ada87a0..eb85ed9d972d0e439fd062684e843651bf377cdc --- a/bin/tools/foamConfigurePaths +++ b/bin/tools/foamConfigurePaths @@ -34,10 +34,10 @@ usage() { cat<<USAGE usage: ${0##*/} - --foamInstall dir specify installation directory (e.g. /opt) - --projectName name specify project name (e.g. openfoam170) - --archOption arch specify architecture option (only 32 or 64 applicable) - --paraviewInstall dir specify ParaView_DIR (e.g. /opt/paraviewopenfoam380) + --foamInstall dir specify installation directory (e.g. /opt) + --projectName name specify project name (e.g. openfoam170) + --archOption arch specify architecture option (only 32 or 64 applicable) + --paraviewInstall dir specify ParaView_DIR (e.g. /opt/paraviewopenfoam380) * hardcode paths to installation @@ -46,6 +46,17 @@ USAGE } +# Function to do replacement on file. Checks if any replacement has been done. +# inlineSed <file> <sedCommand> <description> +_inlineSed() +{ + backup=`tempfile` + cp $1 $backup + sed -i -e "$2" $1 + cmp $1 $backup || usage "Failed : $3" +} + + [ -f etc/bashrc -a -f etc/settings.sh ] || usage "Please run from top-level directory of installation" unset foamInstall projectName archOption paraviewInstall @@ -59,30 +70,40 @@ do ;; --foamInstall) [ "$#" -ge 2 ] || usage "'$1' option requires an argument" - foamInstall="$2" - echo "Replacing foamInstall setting by $foamInstall" - sed -i -e '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" etc/bashrc + foamInstall="$2" + echo "** foamInstall:$foamInstall" + + _inlineSed \ + etc/bashrc \ + '/^[^#]/s@foamInstall=.*@foamInstall='"$foamInstall@" \ + "Replacing foamInstall setting by $foamInstall" shift 2 ;; --projectName) [ "$#" -ge 2 ] || usage "'$1' option requires an argument" - projectName="$2" - echo "Replacing WM_PROJECT_DIR setting by $projectName" - sed -i -e '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" etc/bashrc + projectName="$2" + _inlineSed \ + etc/bashrc \ + '/^[^#]/s@WM_PROJECT_DIR=.*@WM_PROJECT_DIR=$WM_PROJECT_INST_DIR/'"$projectName@" \ + "Replacing WM_PROJECT_DIR setting by $projectName" shift 2 ;; --archOption) [ "$#" -ge 2 ] || usage "'$1' option requires an argument" - archOption="$2" - echo "Replacing WM_ARCH_OPTION setting by $archOption" - sed -i -e '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" etc/bashrc + archOption="$2" + _inlineSed \ + etc/bashrc \ + '/^[^#]/s@: ${WM_ARCH_OPTION:=64}@WM_ARCH_OPTION='"$archOption@" \ + "Replacing WM_ARCH_OPTION setting by $archOption" shift 2 ;; --paraviewInstall) [ "$#" -ge 2 ] || usage "'$1' option requires an argument" - paraviewInstall="$2" - echo "Replacing ParaView_DIR setting by $paraviewInstall" - sed -i -e '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" etc/apps/paraview3/bashrc + paraviewInstall="$2" + _inlineSed \ + etc/apps/paraview3/bashrc \ + '/^[^#]/s@ParaView_DIR=.*@ParaView_DIR='"$paraviewInstall@" \ + "Replacing ParaView_DIR setting by $paraviewInstall" shift 2 ;; *) @@ -97,11 +118,14 @@ done #sed -i -e 's@WM_PROJECT=.*@WM_PROJECT='"$projectName@" etc/bashrc # Replace the WM_MPLIB always -echo "Replacing WM_MPLIB setting by SYSTEMOPENMPI" -sed -i -e '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' etc/bashrc - +_inlineSed \ + etc/bashrc \ + '/^[^#]/s@: ${WM_MPLIB:=.*}@WM_MPLIB=SYSTEMOPENMPI@' \ + "Replacing WM_MPLIB setting by SYSTEMOPENMPI" # Replace the compilerInstall always -echo "Replacing compilerInstall setting by system" -sed -i -e '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' etc/settings.sh +_inlineSed \ + etc/settings.sh \ + '/^[^#]/s@: ${compilerInstall:=.*}@compilerInstall=system@' \ + "Replacing compilerInstall setting by system" #------------------------------------------------------------------------------ diff --git a/src/OSspecific/POSIX/fileMonitor.C b/src/OSspecific/POSIX/fileMonitor.C index 982d5e28dc0321496473f4d4cc1de8d59669781e..47791427785357d4940dc3f8da05f79ae874b96c 100644 --- a/src/OSspecific/POSIX/fileMonitor.C +++ b/src/OSspecific/POSIX/fileMonitor.C @@ -39,7 +39,7 @@ Class #else # include <sys/inotify.h> # include <sys/ioctl.h> - +# include <errno.h> # define EVENT_SIZE ( sizeof (struct inotify_event) ) # define EVENT_LEN (EVENT_SIZE + 16) # define EVENT_BUF_LEN ( 1024 * EVENT_LEN ) @@ -144,7 +144,7 @@ namespace Foam #else //- File descriptor for the inotify instance - int fd; + int inotifyFd_; //- Current watchIDs and corresponding directory id DynamicList<label> dirWatches_; @@ -153,23 +153,45 @@ namespace Foam //- initialise inotify inline fileMonitorWatcher(const label sz = 20) : - fd(inotify_init()), + inotifyFd_(inotify_init()), dirWatches_(sz), dirFiles_(sz) - {} + { + if (inotifyFd_ < 0) + { + static bool hasWarned = false; + if (!hasWarned) + { + hasWarned = true; + WarningIn("fileMonitorWatcher(const label)") + << "Failed allocating an inotify descriptor : " + << string(strerror(errno)) << endl + << " Please increase the number of allowable " + << "inotify instances" << endl + << " (/proc/sys/fs/inotify/max_user_instances" + << " on Linux)" << endl + << " or switch off runTimeModifiable." << endl + << " Continuing without additional file monitoring." + << endl; + } + } + } //- remove all watches inline ~fileMonitorWatcher() { - forAll(dirWatches_, i) + if (inotifyFd_ >= 0) { - if (dirWatches_[i] >= 0) + forAll(dirWatches_, i) { - if (inotify_rm_watch(fd, int(dirWatches_[i]))) + if (dirWatches_[i] >= 0) { - WarningIn("fileMonitor::~fileMonitor()") - << "Failed deleting directory watch " - << dirWatches_[i] << endl; + if (inotify_rm_watch(inotifyFd_, int(dirWatches_[i]))) + { + WarningIn("fileMonitor::~fileMonitor()") + << "Failed deleting directory watch " + << dirWatches_[i] << endl; + } } } } @@ -177,10 +199,15 @@ namespace Foam inline bool addWatch(const label watchFd, const fileName& fName) { + if (inotifyFd_ < 0) + { + return false; + } + // Add/retrieve watch on directory containing file label dirWatchID = inotify_add_watch ( - fd, + inotifyFd_, fName.path().c_str(), IN_CLOSE_WRITE ); @@ -189,7 +216,8 @@ namespace Foam { FatalErrorIn("addWatch(const label, const fileName&)") << "Failed adding watch " << watchFd - << " to directory " << fName + << " to directory " << fName << " due to " + << string(strerror(errno)) << exit(FatalError); } @@ -209,6 +237,11 @@ namespace Foam inline bool removeWatch(const label watchFd) { + if (inotifyFd_ < 0) + { + return false; + } + dirWatches_[watchFd] = -1; return true; } @@ -263,11 +296,11 @@ void Foam::fileMonitor::checkFiles() const fd_set fdSet; // Add notify descriptor to select fd_set FD_ZERO(&fdSet); - FD_SET(watcher_->fd, &fdSet); + FD_SET(watcher_->inotifyFd_, &fdSet); int ready = select ( - watcher_->fd+1, // num filedescriptors in fdSet + watcher_->inotifyFd_+1, // num filedescriptors in fdSet &fdSet, // fdSet with only inotifyFd NULL, // No writefds NULL, // No errorfds @@ -280,15 +313,15 @@ void Foam::fileMonitor::checkFiles() const << "Problem in issuing select." << abort(FatalError); } - else if (FD_ISSET(watcher_->fd, &fdSet)) + else if (FD_ISSET(watcher_->inotifyFd_, &fdSet)) { // Read events - ssize_t nBytes = read(watcher_->fd, buffer, EVENT_BUF_LEN); + ssize_t nBytes = read(watcher_->inotifyFd_, buffer, EVENT_BUF_LEN); if (nBytes < 0) { FatalErrorIn("fileMonitor::updateStates(const fileName&)") - << "read of " << watcher_->fd + << "read of " << watcher_->inotifyFd_ << " failed with " << label(nBytes) << abort(FatalError); } @@ -364,6 +397,7 @@ Foam::label Foam::fileMonitor::addWatch(const fileName& fName) label watchFd; label sz = freeWatchFds_.size(); + if (sz) { watchFd = freeWatchFds_[sz-1]; diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files index 97513b339d4a3a3ecf513e47f95b365a81a98b06..09fa5faa9468bfbfae097d1f3412e990b8c58ff7 100644 --- a/src/OpenFOAM/Make/files +++ b/src/OpenFOAM/Make/files @@ -346,6 +346,7 @@ $(basicPolyPatches)/generic/genericPolyPatch.C constraintPolyPatches = $(polyPatches)/constraint $(constraintPolyPatches)/cyclic/cyclicPolyPatch.C $(constraintPolyPatches)/cyclicSlip/cyclicSlipPolyPatch.C +$(constraintPolyPatches)/oldCyclic/oldCyclicPolyPatch.C $(constraintPolyPatches)/empty/emptyPolyPatch.C $(constraintPolyPatches)/nonuniformTransformCyclic/nonuniformTransformCyclicPolyPatch.C $(constraintPolyPatches)/processorCyclic/processorCyclicPolyPatch.C diff --git a/src/OpenFOAM/db/Time/Time.C b/src/OpenFOAM/db/Time/Time.C index 487576586bf025995778a59e2289efaca093f33f..5a5fbb8fa7f03d935dc16249fa2de156b09192ac 100644 --- a/src/OpenFOAM/db/Time/Time.C +++ b/src/OpenFOAM/db/Time/Time.C @@ -245,13 +245,14 @@ Foam::Time::Time readLibs_(controlDict_, "libs"), functionObjects_(*this) { + setControls(); + // Time objects not registered so do like objectRegistry::checkIn ourselves. if (runTimeModifiable_) { + monitorPtr_.reset(new fileMonitor()); controlDict_.watchIndex() = addWatch(controlDict_.filePath()); } - - setControls(); } @@ -307,14 +308,20 @@ Foam::Time::Time readLibs_(controlDict_, "libs"), functionObjects_(*this) { + setControls(); // Time objects not registered so do like objectRegistry::checkIn ourselves. if (runTimeModifiable_) { - controlDict_.watchIndex() = addWatch(controlDict_.filePath()); - } + monitorPtr_.reset(new fileMonitor()); - setControls(); + // File might not exist yet. + fileName f(controlDict_.filePath()); + if (f != fileName::null) + { + controlDict_.watchIndex() = addWatch(f); + } + } } @@ -388,18 +395,18 @@ Foam::Time::~Time() Foam::label Foam::Time::addWatch(const fileName& fName) const { - return monitor_.addWatch(fName); + return monitorPtr_().addWatch(fName); } bool Foam::Time::removeWatch(const label watchIndex) const { - return monitor_.removeWatch(watchIndex); + return monitorPtr_().removeWatch(watchIndex); } const Foam::fileName& Foam::Time::getFile(const label watchIndex) const { - return monitor_.getFile(watchIndex); + return monitorPtr_().getFile(watchIndex); } @@ -408,13 +415,13 @@ Foam::fileMonitor::fileState Foam::Time::getState const label watchFd ) const { - return monitor_.getState(watchFd); + return monitorPtr_().getState(watchFd); } void Foam::Time::setUnmodified(const label watchFd) const { - monitor_.setUnmodified(watchFd); + monitorPtr_().setUnmodified(watchFd); } diff --git a/src/OpenFOAM/db/Time/Time.H b/src/OpenFOAM/db/Time/Time.H index c644c757be1a971ec5f0f7f79ddb24834c66f8f1..f9e38b58eed5a14756f868be8c5e5a62954dd111 100644 --- a/src/OpenFOAM/db/Time/Time.H +++ b/src/OpenFOAM/db/Time/Time.H @@ -73,7 +73,7 @@ class Time // Private data //- file-change monitor for all registered files - mutable fileMonitor monitor_; + mutable autoPtr<fileMonitor> monitorPtr_; //- The controlDict IOdictionary controlDict_; diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C index e7c38a6aa6d3f7652a3ba99af3b1a9d65b3da61e..86db2e1c2f0ae40407fd95b9fe68d49b8405514c 100644 --- a/src/OpenFOAM/db/Time/TimeIO.C +++ b/src/OpenFOAM/db/Time/TimeIO.C @@ -211,10 +211,7 @@ void Foam::Time::readModifiedObjects() // valid filePath). // Note: requires same ordering in objectRegistries on different // processors! - monitor_.updateStates(Pstream::parRun()); - -//Pout<< "Time : runTimeModifiable_ and watchIndex:" -// << controlDict_.watchIndex() << endl; + monitorPtr_().updateStates(Pstream::parRun()); // Time handling is special since controlDict_ is the one dictionary // that is not registered to any database. diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C new file mode 100644 index 0000000000000000000000000000000000000000..218debfa5b0b0df7e7571e368309a5e2c209f00b --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C @@ -0,0 +1,1293 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "oldCyclicPolyPatch.H" +#include "addToRunTimeSelectionTable.H" +#include "polyBoundaryMesh.H" +#include "polyMesh.H" +#include "demandDrivenData.H" +#include "OFstream.H" +#include "patchZones.H" +#include "matchPoints.H" +#include "Time.H" +#include "transformList.H" +#include "cyclicPolyPatch.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(oldCyclicPolyPatch, 0); + + addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, word); + addToRunTimeSelectionTable(polyPatch, oldCyclicPolyPatch, dictionary); + + +template<> +const char* NamedEnum<oldCyclicPolyPatch::transformType, 3>::names[] = +{ + "unknown", + "rotational", + "translational" +}; + +const NamedEnum<oldCyclicPolyPatch::transformType, 3> + oldCyclicPolyPatch::transformTypeNames; +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::pointField Foam::oldCyclicPolyPatch::calcFaceCentres +( + const UList<face>& faces, + const pointField& points +) +{ + pointField ctrs(faces.size()); + + forAll(faces, faceI) + { + ctrs[faceI] = faces[faceI].centre(points); + } + + return ctrs; +} + + +Foam::pointField Foam::oldCyclicPolyPatch::getAnchorPoints +( + const UList<face>& faces, + const pointField& points +) +{ + pointField anchors(faces.size()); + + forAll(faces, faceI) + { + anchors[faceI] = points[faces[faceI][0]]; + } + + return anchors; +} + + +Foam::label Foam::oldCyclicPolyPatch::findMaxArea +( + const pointField& points, + const faceList& faces +) +{ + label maxI = -1; + scalar maxAreaSqr = -GREAT; + + forAll(faces, faceI) + { + scalar areaSqr = magSqr(faces[faceI].normal(points)); + + if (areaSqr > maxAreaSqr) + { + maxAreaSqr = areaSqr; + maxI = faceI; + } + } + return maxI; +} + + +// Get geometric zones of patch by looking at normals. +// Method 1: any edge with sharpish angle is edge between two halves. +// (this will handle e.g. wedge geometries). +// Also two fully disconnected regions will be handled this way. +// Method 2: sort faces into two halves based on face normal. +bool Foam::oldCyclicPolyPatch::getGeometricHalves +( + const primitivePatch& pp, + labelList& half0ToPatch, + labelList& half1ToPatch +) const +{ + // Calculate normals + const vectorField& faceNormals = pp.faceNormals(); + + // Find edges with sharp angles. + boolList regionEdge(pp.nEdges(), false); + + const labelListList& edgeFaces = pp.edgeFaces(); + + label nRegionEdges = 0; + + forAll(edgeFaces, edgeI) + { + const labelList& eFaces = edgeFaces[edgeI]; + + // Check manifold edges for sharp angle. + // (Non-manifold already handled by patchZones) + if (eFaces.size() == 2) + { + if ((faceNormals[eFaces[0]] & faceNormals[eFaces[1]])< featureCos_) + { + regionEdge[edgeI] = true; + + nRegionEdges++; + } + } + } + + + // For every face determine zone it is connected to (without crossing + // any regionEdge) + patchZones ppZones(pp, regionEdge); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::getGeometricHalves : " + << "Found " << nRegionEdges << " edges on patch " << name() + << " where the cos of the angle between two connected faces" + << " was less than " << featureCos_ << nl + << "Patch divided by these and by single sides edges into " + << ppZones.nZones() << " parts." << endl; + + + // Dumping zones to obj files. + + labelList nZoneFaces(ppZones.nZones()); + + for (label zoneI = 0; zoneI < ppZones.nZones(); zoneI++) + { + OFstream stream + ( + boundaryMesh().mesh().time().path() + /name()+"_zone_"+Foam::name(zoneI)+".obj" + ); + Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing zone " + << zoneI << " face centres to OBJ file " << stream.name() + << endl; + + labelList zoneFaces(findIndices(ppZones, zoneI)); + + forAll(zoneFaces, i) + { + writeOBJ(stream, pp[zoneFaces[i]].centre(pp.points())); + } + + nZoneFaces[zoneI] = zoneFaces.size(); + } + } + + + if (ppZones.nZones() == 2) + { + half0ToPatch = findIndices(ppZones, 0); + half1ToPatch = findIndices(ppZones, 1); + } + else + { + if (debug) + { + Pout<< "oldCyclicPolyPatch::getGeometricHalves :" + << " falling back to face-normal comparison" << endl; + } + label n0Faces = 0; + half0ToPatch.setSize(pp.size()); + + label n1Faces = 0; + half1ToPatch.setSize(pp.size()); + + // Compare to face 0 normal. + forAll(faceNormals, faceI) + { + if ((faceNormals[faceI] & faceNormals[0]) > 0) + { + half0ToPatch[n0Faces++] = faceI; + } + else + { + half1ToPatch[n1Faces++] = faceI; + } + } + half0ToPatch.setSize(n0Faces); + half1ToPatch.setSize(n1Faces); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::getGeometricHalves :" + << " Number of faces per zone:(" + << n0Faces << ' ' << n1Faces << ')' << endl; + } + } + + if (half0ToPatch.size() != half1ToPatch.size()) + { + fileName casePath(boundaryMesh().mesh().time().path()); + + // Dump halves + { + fileName nm0(casePath/name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, UIndirectList<face>(pp, half0ToPatch)(), pp.points()); + + fileName nm1(casePath/name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, UIndirectList<face>(pp, half1ToPatch)(), pp.points()); + } + + // Dump face centres + { + OFstream str0(casePath/name()+"_half0.obj"); + Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half0" + << " face centres to OBJ file " << str0.name() << endl; + + forAll(half0ToPatch, i) + { + writeOBJ(str0, pp[half0ToPatch[i]].centre(pp.points())); + } + + OFstream str1(casePath/name()+"_half1.obj"); + Pout<< "oldCyclicPolyPatch::getGeometricHalves : Writing half1" + << " face centres to OBJ file " << str1.name() << endl; + forAll(half1ToPatch, i) + { + writeOBJ(str1, pp[half1ToPatch[i]].centre(pp.points())); + } + } + + SeriousErrorIn + ( + "oldCyclicPolyPatch::getGeometricHalves" + "(const primitivePatch&, labelList&, labelList&) const" + ) << "Patch " << name() << " gets decomposed in two zones of" + << "inequal size: " << half0ToPatch.size() + << " and " << half1ToPatch.size() << endl + << "This means that the patch is either not two separate regions" + << " or one region where the angle between the different regions" + << " is not sufficiently sharp." << endl + << "Please adapt the featureCos setting." << endl + << "Continuing with incorrect face ordering from now on!" << endl; + + return false; + } + else + { + return true; + } +} + + +// Given a split of faces into left and right half calculate the centres +// and anchor points. Transform the left points so they align with the +// right ones. +void Foam::oldCyclicPolyPatch::getCentresAndAnchors +( + const primitivePatch& pp, + const faceList& half0Faces, + const faceList& half1Faces, + + pointField& ppPoints, + pointField& half0Ctrs, + pointField& half1Ctrs, + pointField& anchors0, + scalarField& tols +) const +{ + // Get geometric data on both halves. + half0Ctrs = calcFaceCentres(half0Faces, pp.points()); + anchors0 = getAnchorPoints(half0Faces, pp.points()); + half1Ctrs = calcFaceCentres(half1Faces, pp.points()); + + switch (transform_) + { + case ROTATIONAL: + { + label face0 = getConsistentRotationFace(half0Ctrs); + label face1 = getConsistentRotationFace(half1Ctrs); + + vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_); + vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_); + n0 /= mag(n0) + VSMALL; + n1 /= mag(n1) + VSMALL; + + if (debug) + { + Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :" + << " Specified rotation :" + << " n0:" << n0 << " n1:" << n1 << endl; + } + + // Rotation (around origin) + const tensor reverseT(rotationTensor(n0, -n1)); + + // Rotation + forAll(half0Ctrs, faceI) + { + half0Ctrs[faceI] = Foam::transform(reverseT, half0Ctrs[faceI]); + anchors0[faceI] = Foam::transform(reverseT, anchors0[faceI]); + } + + ppPoints = Foam::transform(reverseT, pp.points()); + + break; + } + //- Problem: usually specified translation is not accurate enough + //- to get proper match so keep automatic determination over here. + //case TRANSLATIONAL: + //{ + // // Transform 0 points. + // + // if (debug) + // { + // Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :" + // << "Specified translation : " << separationVector_ + // << endl; + // } + // + // half0Ctrs += separationVector_; + // anchors0 += separationVector_; + // break; + //} + default: + { + // Assumes that cyclic is planar. This is also the initial + // condition for patches without faces. + + // Determine the face with max area on both halves. These + // two faces are used to determine the transformation tensors + label max0I = findMaxArea(pp.points(), half0Faces); + vector n0 = half0Faces[max0I].normal(pp.points()); + n0 /= mag(n0) + VSMALL; + + label max1I = findMaxArea(pp.points(), half1Faces); + vector n1 = half1Faces[max1I].normal(pp.points()); + n1 /= mag(n1) + VSMALL; + + if (mag(n0 & n1) < 1-coupledPolyPatch::matchTol) + { + if (debug) + { + Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :" + << " Detected rotation :" + << " n0:" << n0 << " n1:" << n1 << endl; + } + + // Rotation (around origin) + const tensor reverseT(rotationTensor(n0, -n1)); + + // Rotation + forAll(half0Ctrs, faceI) + { + half0Ctrs[faceI] = Foam::transform + ( + reverseT, + half0Ctrs[faceI] + ); + anchors0[faceI] = Foam::transform + ( + reverseT, + anchors0[faceI] + ); + } + ppPoints = Foam::transform(reverseT, pp.points()); + } + else + { + // Parallel translation. Get average of all used points. + + primitiveFacePatch half0(half0Faces, pp.points()); + const pointField& half0Pts = half0.localPoints(); + const point ctr0(sum(half0Pts)/half0Pts.size()); + + primitiveFacePatch half1(half1Faces, pp.points()); + const pointField& half1Pts = half1.localPoints(); + const point ctr1(sum(half1Pts)/half1Pts.size()); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::getCentresAndAnchors :" + << " Detected translation :" + << " n0:" << n0 << " n1:" << n1 + << " ctr0:" << ctr0 << " ctr1:" << ctr1 << endl; + } + + half0Ctrs += ctr1 - ctr0; + anchors0 += ctr1 - ctr0; + ppPoints = pp.points() + ctr1 - ctr0; + } + break; + } + } + + + // Calculate typical distance per face + tols = calcFaceTol(half1Faces, pp.points(), half1Ctrs); +} + + +// Calculates faceMap and rotation. Returns true if all ok. +bool Foam::oldCyclicPolyPatch::matchAnchors +( + const bool report, + const primitivePatch& pp, + const labelList& half0ToPatch, + const pointField& anchors0, + + const labelList& half1ToPatch, + const faceList& half1Faces, + const labelList& from1To0, + + const scalarField& tols, + + labelList& faceMap, + labelList& rotation +) const +{ + // Set faceMap such that half0 faces get first and corresponding half1 + // faces last. + + forAll(half0ToPatch, half0FaceI) + { + // Label in original patch + label patchFaceI = half0ToPatch[half0FaceI]; + + faceMap[patchFaceI] = half0FaceI; + + // No rotation + rotation[patchFaceI] = 0; + } + + bool fullMatch = true; + + forAll(from1To0, half1FaceI) + { + label patchFaceI = half1ToPatch[half1FaceI]; + + // This face has to match the corresponding one on half0. + label half0FaceI = from1To0[half1FaceI]; + + label newFaceI = half0FaceI + pp.size()/2; + + faceMap[patchFaceI] = newFaceI; + + // Rotate patchFaceI such that its f[0] aligns with that of + // the corresponding face + // (which after shuffling will be at position half0FaceI) + + const point& wantedAnchor = anchors0[half0FaceI]; + + rotation[newFaceI] = getRotation + ( + pp.points(), + half1Faces[half1FaceI], + wantedAnchor, + tols[half1FaceI] + ); + + if (rotation[newFaceI] == -1) + { + fullMatch = false; + + if (report) + { + const face& f = half1Faces[half1FaceI]; + SeriousErrorIn + ( + "oldCyclicPolyPatch::matchAnchors(..)" + ) << "Patch:" << name() << " : " + << "Cannot find point on face " << f + << " with vertices:" + << UIndirectList<point>(pp.points(), f)() + << " that matches point " << wantedAnchor + << " when matching the halves of cyclic patch " << name() + << endl + << "Continuing with incorrect face ordering from now on!" + << endl; + } + } + } + return fullMatch; +} + + +Foam::label Foam::oldCyclicPolyPatch::getConsistentRotationFace +( + const pointField& faceCentres +) const +{ + const scalarField magRadSqr = + magSqr((faceCentres - rotationCentre_) ^ rotationAxis_); + scalarField axisLen = (faceCentres - rotationCentre_) & rotationAxis_; + axisLen = axisLen - min(axisLen); + const scalarField magLenSqr = magRadSqr + axisLen*axisLen; + + label rotFace = -1; + scalar maxMagLenSqr = -GREAT; + scalar maxMagRadSqr = -GREAT; + forAll(faceCentres, i) + { + if (magLenSqr[i] >= maxMagLenSqr) + { + if (magRadSqr[i] > maxMagRadSqr) + { + rotFace = i; + maxMagLenSqr = magLenSqr[i]; + maxMagRadSqr = magRadSqr[i]; + } + } + } + + if (debug) + { + Info<< "getConsistentRotationFace(const pointField&)" << nl + << " rotFace = " << rotFace << nl + << " point = " << faceCentres[rotFace] << endl; + } + + return rotFace; +} + + +// * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * * * // + +Foam::oldCyclicPolyPatch::oldCyclicPolyPatch +( + const word& name, + const label size, + const label start, + const label index, + const polyBoundaryMesh& bm +) +: + coupledPolyPatch(name, size, start, index, bm), + featureCos_(0.9), + transform_(UNKNOWN), + rotationAxis_(vector::zero), + rotationCentre_(point::zero), + separationVector_(vector::zero) +{} + + +Foam::oldCyclicPolyPatch::oldCyclicPolyPatch +( + const word& name, + const dictionary& dict, + const label index, + const polyBoundaryMesh& bm +) +: + coupledPolyPatch(name, dict, index, bm), + featureCos_(0.9), + transform_(UNKNOWN), + rotationAxis_(vector::zero), + rotationCentre_(point::zero), + separationVector_(vector::zero) +{ + if (dict.found("neighbourPatch")) + { + FatalIOErrorIn + ( + "oldCyclicPolyPatch::oldCyclicPolyPatch\n" + "(\n" + " const word& name,\n" + " const dictionary& dict,\n" + " const label index,\n" + " const polyBoundaryMesh& bm\n" + ")", + dict + ) << "Found \"neighbourPatch\" entry when reading cyclic patch " + << name << endl + << "Is this mesh already with split cyclics?" << endl + << "If so run a newer version that supports it" + << ", if not comment out the \"neighbourPatch\" entry and re-run" + << exit(FatalIOError); + } + + dict.readIfPresent("featureCos", featureCos_); + + if (dict.found("transform")) + { + transform_ = transformTypeNames.read(dict.lookup("transform")); + switch (transform_) + { + case ROTATIONAL: + { + dict.lookup("rotationAxis") >> rotationAxis_; + dict.lookup("rotationCentre") >> rotationCentre_; + break; + } + case TRANSLATIONAL: + { + dict.lookup("separationVector") >> separationVector_; + break; + } + default: + { + // no additional info required + } + } + } +} + + +Foam::oldCyclicPolyPatch::oldCyclicPolyPatch +( + const oldCyclicPolyPatch& pp, + const polyBoundaryMesh& bm +) +: + coupledPolyPatch(pp, bm), + featureCos_(pp.featureCos_), + transform_(pp.transform_), + rotationAxis_(pp.rotationAxis_), + rotationCentre_(pp.rotationCentre_), + separationVector_(pp.separationVector_) +{} + + +Foam::oldCyclicPolyPatch::oldCyclicPolyPatch +( + const oldCyclicPolyPatch& pp, + const polyBoundaryMesh& bm, + const label index, + const label newSize, + const label newStart +) +: + coupledPolyPatch(pp, bm, index, newSize, newStart), + featureCos_(pp.featureCos_), + transform_(pp.transform_), + rotationAxis_(pp.rotationAxis_), + rotationCentre_(pp.rotationCentre_), + separationVector_(pp.separationVector_) +{} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::oldCyclicPolyPatch::~oldCyclicPolyPatch() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::oldCyclicPolyPatch::initGeometry(PstreamBuffers& pBufs) +{ + polyPatch::initGeometry(pBufs); +} + + +void Foam::oldCyclicPolyPatch::calcGeometry +( + const primitivePatch& referPatch, + const UList<point>& thisCtrs, + const UList<point>& thisAreas, + const UList<point>& thisCc, + const UList<point>& nbrCtrs, + const UList<point>& nbrAreas, + const UList<point>& nbrCc +) +{} + + +void Foam::oldCyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs) +{} + + +void Foam::oldCyclicPolyPatch::initMovePoints +( + PstreamBuffers& pBufs, + const pointField& p +) +{ + polyPatch::initMovePoints(pBufs, p); +} + + +void Foam::oldCyclicPolyPatch::movePoints +( + PstreamBuffers& pBufs, + const pointField& p +) +{ + polyPatch::movePoints(pBufs, p); +} + + +void Foam::oldCyclicPolyPatch::initUpdateMesh(PstreamBuffers& pBufs) +{ + polyPatch::initUpdateMesh(pBufs); +} + + +void Foam::oldCyclicPolyPatch::updateMesh(PstreamBuffers& pBufs) +{ + polyPatch::updateMesh(pBufs); +} + + +void Foam::oldCyclicPolyPatch::initOrder +( + PstreamBuffers&, + const primitivePatch& pp +) const +{} + + +// Return new ordering. Ordering is -faceMap: for every face index +// the new face -rotation:for every new face the clockwise shift +// of the original face. Return false if nothing changes (faceMap +// is identity, rotation is 0) +bool Foam::oldCyclicPolyPatch::order +( + PstreamBuffers&, + const primitivePatch& pp, + labelList& faceMap, + labelList& rotation +) const +{ + faceMap.setSize(pp.size()); + faceMap = -1; + + rotation.setSize(pp.size()); + rotation = 0; + + if (pp.empty()) + { + // No faces, nothing to change. + return false; + } + + if (pp.size()&1) + { + FatalErrorIn("oldCyclicPolyPatch::order(..)") + << "Size of cyclic " << name() << " should be a multiple of 2" + << ". It is " << pp.size() << abort(FatalError); + } + + label halfSize = pp.size()/2; + + // Supplied primitivePatch already with new points. + // Cyclics are limited to one transformation tensor + // currently anyway (i.e. straight plane) so should not be too big a + // problem. + + + // Indices of faces on half0 + labelList half0ToPatch; + // Indices of faces on half1 + labelList half1ToPatch; + + + // 1. Test if already correctly oriented by starting from trivial ordering. + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + half0ToPatch = identity(halfSize); + half1ToPatch = half0ToPatch + halfSize; + + // Get faces + faceList half0Faces(UIndirectList<face>(pp, half0ToPatch)); + faceList half1Faces(UIndirectList<face>(pp, half1ToPatch)); + + // Get geometric quantities + pointField half0Ctrs, half1Ctrs, anchors0, ppPoints; + scalarField tols; + getCentresAndAnchors + ( + pp, + half0Faces, + half1Faces, + + ppPoints, + half0Ctrs, + half1Ctrs, + anchors0, + tols + ); + + // Geometric match of face centre vectors + labelList from1To0; + bool matchedAll = matchPoints + ( + half1Ctrs, + half0Ctrs, + tols, + false, + from1To0 + ); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::order : test if already ordered:" + << matchedAll << endl; + + // Dump halves + fileName nm0("match1_"+name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0Faces, ppPoints); + + fileName nm1("match1_"+name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1Faces, pp.points()); + + OFstream ccStr + ( + boundaryMesh().mesh().time().path() + /"match1_"+ name() + "_faceCentres.obj" + ); + Pout<< "oldCyclicPolyPatch::order : " + << "Dumping currently found cyclic match as lines between" + << " corresponding face centres to file " << ccStr.name() + << endl; + + // Recalculate untransformed face centres + //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points()); + label vertI = 0; + + forAll(half1Ctrs, i) + { + //if (from1To0[i] != -1) + { + // Write edge between c1 and c0 + //const point& c0 = rawHalf0Ctrs[from1To0[i]]; + //const point& c0 = half0Ctrs[from1To0[i]]; + const point& c0 = half0Ctrs[i]; + const point& c1 = half1Ctrs[i]; + writeOBJ(ccStr, c0, c1, vertI); + } + } + } + + + // 2. Ordered in pairs (so 0,1 coupled and 2,3 etc.) + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (!matchedAll) + { + label faceI = 0; + for (label i = 0; i < halfSize; i++) + { + half0ToPatch[i] = faceI++; + half1ToPatch[i] = faceI++; + } + + // And redo all matching + half0Faces = UIndirectList<face>(pp, half0ToPatch); + half1Faces = UIndirectList<face>(pp, half1ToPatch); + + getCentresAndAnchors + ( + pp, + half0Faces, + half1Faces, + + ppPoints, + half0Ctrs, + half1Ctrs, + anchors0, + tols + ); + + // Geometric match of face centre vectors + matchedAll = matchPoints + ( + half1Ctrs, + half0Ctrs, + tols, + false, + from1To0 + ); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::order : test if pairwise ordered:" + << matchedAll << endl; + // Dump halves + fileName nm0("match2_"+name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0Faces, ppPoints); + + fileName nm1("match2_"+name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1Faces, pp.points()); + + OFstream ccStr + ( + boundaryMesh().mesh().time().path() + /"match2_"+name()+"_faceCentres.obj" + ); + Pout<< "oldCyclicPolyPatch::order : " + << "Dumping currently found cyclic match as lines between" + << " corresponding face centres to file " << ccStr.name() + << endl; + + // Recalculate untransformed face centres + label vertI = 0; + + forAll(half1Ctrs, i) + { + if (from1To0[i] != -1) + { + // Write edge between c1 and c0 + const point& c0 = half0Ctrs[from1To0[i]]; + const point& c1 = half1Ctrs[i]; + writeOBJ(ccStr, c0, c1, vertI); + } + } + } + } + + + // 3. Baffles(coincident faces) converted into cyclics (e.g. jump) + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (!matchedAll) + { + label baffleI = 0; + + forAll(pp, faceI) + { + const face& f = pp.localFaces()[faceI]; + const labelList& pFaces = pp.pointFaces()[f[0]]; + + label matchedFaceI = -1; + + forAll(pFaces, i) + { + label otherFaceI = pFaces[i]; + + if (otherFaceI > faceI) + { + const face& otherF = pp.localFaces()[otherFaceI]; + + // Note: might pick up two similar oriented faces + // (but that is illegal anyway) + if (f == otherF) + { + matchedFaceI = otherFaceI; + break; + } + } + } + + if (matchedFaceI != -1) + { + half0ToPatch[baffleI] = faceI; + half1ToPatch[baffleI] = matchedFaceI; + baffleI++; + } + } + + if (baffleI == halfSize) + { + // And redo all matching + half0Faces = UIndirectList<face>(pp, half0ToPatch); + half1Faces = UIndirectList<face>(pp, half1ToPatch); + + getCentresAndAnchors + ( + pp, + half0Faces, + half1Faces, + + ppPoints, + half0Ctrs, + half1Ctrs, + anchors0, + tols + ); + + // Geometric match of face centre vectors + matchedAll = matchPoints + ( + half1Ctrs, + half0Ctrs, + tols, + false, + from1To0 + ); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::order : test if baffles:" + << matchedAll << endl; + // Dump halves + fileName nm0("match3_"+name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0Faces, ppPoints); + + fileName nm1("match3_"+name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1Faces, pp.points()); + + OFstream ccStr + ( + boundaryMesh().mesh().time().path() + /"match3_"+ name() + "_faceCentres.obj" + ); + Pout<< "oldCyclicPolyPatch::order : " + << "Dumping currently found cyclic match as lines between" + << " corresponding face centres to file " << ccStr.name() + << endl; + + // Recalculate untransformed face centres + label vertI = 0; + + forAll(half1Ctrs, i) + { + if (from1To0[i] != -1) + { + // Write edge between c1 and c0 + const point& c0 = half0Ctrs[from1To0[i]]; + const point& c1 = half1Ctrs[i]; + writeOBJ(ccStr, c0, c1, vertI); + } + } + } + } + } + + + // 4. Automatic geometric ordering + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + if (!matchedAll) + { + // Split faces according to feature angle or topology + bool okSplit = getGeometricHalves(pp, half0ToPatch, half1ToPatch); + + if (!okSplit) + { + // Did not split into two equal parts. + return false; + } + + // And redo all matching + half0Faces = UIndirectList<face>(pp, half0ToPatch); + half1Faces = UIndirectList<face>(pp, half1ToPatch); + + getCentresAndAnchors + ( + pp, + half0Faces, + half1Faces, + + ppPoints, + half0Ctrs, + half1Ctrs, + anchors0, + tols + ); + + // Geometric match of face centre vectors + matchedAll = matchPoints + ( + half1Ctrs, + half0Ctrs, + tols, + false, + from1To0 + ); + + if (debug) + { + Pout<< "oldCyclicPolyPatch::order : automatic ordering result:" + << matchedAll << endl; + // Dump halves + fileName nm0("match4_"+name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0Faces, ppPoints); + + fileName nm1("match4_"+name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1Faces, pp.points()); + + OFstream ccStr + ( + boundaryMesh().mesh().time().path() + /"match4_"+ name() + "_faceCentres.obj" + ); + Pout<< "oldCyclicPolyPatch::order : " + << "Dumping currently found cyclic match as lines between" + << " corresponding face centres to file " << ccStr.name() + << endl; + + // Recalculate untransformed face centres + label vertI = 0; + + forAll(half1Ctrs, i) + { + if (from1To0[i] != -1) + { + // Write edge between c1 and c0 + const point& c0 = half0Ctrs[from1To0[i]]; + const point& c1 = half1Ctrs[i]; + writeOBJ(ccStr, c0, c1, vertI); + } + } + } + } + + + if (!matchedAll || debug) + { + // Dump halves + fileName nm0(name()+"_half0_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half0" + << " faces to OBJ file " << nm0 << endl; + writeOBJ(nm0, half0Faces, pp.points()); + + fileName nm1(name()+"_half1_faces.obj"); + Pout<< "oldCyclicPolyPatch::order : Writing half1" + << " faces to OBJ file " << nm1 << endl; + writeOBJ(nm1, half1Faces, pp.points()); + + OFstream ccStr + ( + boundaryMesh().mesh().time().path() + /name() + "_faceCentres.obj" + ); + Pout<< "oldCyclicPolyPatch::order : " + << "Dumping currently found cyclic match as lines between" + << " corresponding face centres to file " << ccStr.name() + << endl; + + // Recalculate untransformed face centres + //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points()); + label vertI = 0; + + forAll(half1Ctrs, i) + { + if (from1To0[i] != -1) + { + // Write edge between c1 and c0 + //const point& c0 = rawHalf0Ctrs[from1To0[i]]; + const point& c0 = half0Ctrs[from1To0[i]]; + const point& c1 = half1Ctrs[i]; + writeOBJ(ccStr, c0, c1, vertI); + } + } + } + + + if (!matchedAll) + { + SeriousErrorIn + ( + "oldCyclicPolyPatch::order" + "(const primitivePatch&, labelList&, labelList&) const" + ) << "Patch:" << name() << " : " + << "Cannot match vectors to faces on both sides of patch" << endl + << " Perhaps your faces do not match?" + << " The obj files written contain the current match." << endl + << " Continuing with incorrect face ordering from now on!" + << endl; + + return false; + } + + + // Set faceMap such that half0 faces get first and corresponding half1 + // faces last. + matchAnchors + ( + true, // report if anchor matching error + pp, + half0ToPatch, + anchors0, + half1ToPatch, + half1Faces, + from1To0, + tols, + faceMap, + rotation + ); + + // Return false if no change neccesary, true otherwise. + + forAll(faceMap, faceI) + { + if (faceMap[faceI] != faceI || rotation[faceI] != 0) + { + return true; + } + } + + return false; +} + + +void Foam::oldCyclicPolyPatch::write(Ostream& os) const +{ + // Replacement of polyPatch::write to write 'cyclic' instead of type(): + os.writeKeyword("type") << cyclicPolyPatch::typeName + << token::END_STATEMENT << nl; + patchIdentifier::write(os); + os.writeKeyword("nFaces") << size() << token::END_STATEMENT << nl; + os.writeKeyword("startFace") << start() << token::END_STATEMENT << nl; + + + os.writeKeyword("featureCos") << featureCos_ << token::END_STATEMENT << nl; + switch (transform_) + { + case ROTATIONAL: + { + os.writeKeyword("transform") << transformTypeNames[transform_] + << token::END_STATEMENT << nl; + os.writeKeyword("rotationAxis") << rotationAxis_ + << token::END_STATEMENT << nl; + os.writeKeyword("rotationCentre") << rotationCentre_ + << token::END_STATEMENT << nl; + break; + } + case TRANSLATIONAL: + { + os.writeKeyword("transform") << transformTypeNames[transform_] + << token::END_STATEMENT << nl; + os.writeKeyword("separationVector") << separationVector_ + << token::END_STATEMENT << nl; + break; + } + default: + { + // no additional info to write + } + } + + WarningIn("oldCyclicPolyPatch::write(Ostream& os) const") + << "Please run foamUpgradeCyclics to convert these old-style" + << " cyclics into two separate cyclics patches." + << endl; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H new file mode 100644 index 0000000000000000000000000000000000000000..f8cfb169c0eaf609da95868c462f9df7b006eda1 --- /dev/null +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.H @@ -0,0 +1,318 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2010 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 <http://www.gnu.org/licenses/>. + +Class + Foam::oldCyclicPolyPatch + +Description + 'old' style cyclic polyPatch with all faces in single patch. Does ordering + but cannot be used to run. Writes 'type cyclic' so foamUpgradeCyclics + can be run afterwards. + Used to get cyclics from mesh converters that assume cyclics in single + patch (e.g. fluent3DMeshToFoam) + +SourceFiles + oldCyclicPolyPatch.C + +\*---------------------------------------------------------------------------*/ + +#ifndef oldCyclicPolyPatch_H +#define oldCyclicPolyPatch_H + +#include "coupledPolyPatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class oldCyclicPolyPatch Declaration +\*---------------------------------------------------------------------------*/ + +class oldCyclicPolyPatch +: + public coupledPolyPatch +{ +public: + + enum transformType + { + UNKNOWN, + ROTATIONAL, + TRANSLATIONAL + }; + static const NamedEnum<transformType, 3> transformTypeNames; + + +private: + + // Private data + + //- Morph:angle between normals of neighbouring faces. + // Used to split cyclic into halves. + scalar featureCos_; + + //- Type of transformation - rotational or translational + transformType transform_; + + // For rotation + + //- Axis of rotation for rotational cyclics + vector rotationAxis_; + + //- point on axis of rotation for rotational cyclics + point rotationCentre_; + + // For translation + + //- Translation vector + vector separationVector_; + + + // Private member functions + + //- Find amongst selected faces the one with the largest area + static label findMaxArea(const pointField&, const faceList&); + + void calcTransforms(); + + //- Calculate face centres + static pointField calcFaceCentres + ( + const UList<face>&, + const pointField& + ); + + //- Get f[0] for all faces + static pointField getAnchorPoints + ( + const UList<face>&, + const pointField& + ); + + // Face ordering + + //- Find the two parts of the faces of pp using feature edges. + // Returns true if successfull. + bool getGeometricHalves + ( + const primitivePatch&, + labelList&, + labelList& + ) const; + + //- Calculate geometric factors of the two halves. + void getCentresAndAnchors + ( + const primitivePatch&, + const faceList& half0Faces, + const faceList& half1Faces, + + pointField& ppPoints, + pointField& half0Ctrs, + pointField& half1Ctrs, + pointField& anchors0, + scalarField& tols + ) const; + + //- Given matched faces matches the anchor point. Sets faceMap, + // rotation. Returns true if all matched. + bool matchAnchors + ( + const bool report, + const primitivePatch&, + const labelList&, + const pointField&, + const labelList&, + const faceList&, + const labelList&, + const scalarField&, + + labelList& faceMap, + labelList& rotation + ) const; + + //- For rotational cases, try to find a unique face on each side + // of the cyclic. + label getConsistentRotationFace + ( + const pointField& faceCentres + ) const; + + +protected: + + // Protected Member functions + + //- Initialise the calculation of the patch geometry + virtual void initGeometry(PstreamBuffers&); + + //- Calculate the patch geometry + virtual void calcGeometry(PstreamBuffers&); + + //- Initialise the patches for moving points + virtual void initMovePoints(PstreamBuffers&, const pointField&); + + //- Correct patches after moving points + virtual void movePoints(PstreamBuffers&, const pointField&); + + //- Initialise the update of the patch topology + virtual void initUpdateMesh(PstreamBuffers&); + + //- Update of the patch topology + virtual void updateMesh(PstreamBuffers&); + +public: + + //- Runtime type information + TypeName("oldCyclic"); + + + // Constructors + + //- Construct from components + oldCyclicPolyPatch + ( + const word& name, + const label size, + const label start, + const label index, + const polyBoundaryMesh& bm + ); + + //- Construct from dictionary + oldCyclicPolyPatch + ( + const word& name, + const dictionary& dict, + const label index, + const polyBoundaryMesh& bm + ); + + //- Construct as copy, resetting the boundary mesh + oldCyclicPolyPatch(const oldCyclicPolyPatch&, const polyBoundaryMesh&); + + //- Construct given the original patch and resetting the + // face list and boundary mesh information + oldCyclicPolyPatch + ( + const oldCyclicPolyPatch& pp, + const polyBoundaryMesh& bm, + const label index, + const label newSize, + const label newStart + ); + + //- Construct and return a clone, resetting the boundary mesh + virtual autoPtr<polyPatch> clone(const polyBoundaryMesh& bm) const + { + return autoPtr<polyPatch>(new oldCyclicPolyPatch(*this, bm)); + } + + //- Construct and return a clone, resetting the face list + // and boundary mesh + virtual autoPtr<polyPatch> clone + ( + const polyBoundaryMesh& bm, + const label index, + const label newSize, + const label newStart + ) const + { + return autoPtr<polyPatch> + ( + new oldCyclicPolyPatch(*this, bm, index, newSize, newStart) + ); + } + + + // Destructor + + virtual ~oldCyclicPolyPatch(); + + + // Member Functions + + // Access + + //- Does this side own the patch ? + virtual bool owner() const + { + notImplemented("oldCyclicPolyPatch::owner()"); + return true; + } + + //- Transform a patch-based position from other side to this side + virtual void transformPosition(pointField& l) const + { + notImplemented("transformPosition(pointField&)"); + } + + //- Calculate the patch geometry + virtual void calcGeometry + ( + const primitivePatch& referPatch, + const UList<point>& thisCtrs, + const UList<point>& thisAreas, + const UList<point>& thisCc, + const UList<point>& nbrCtrs, + const UList<point>& nbrAreas, + const UList<point>& nbrCc + ); + + //- Initialize ordering for primitivePatch. Does not + // refer to *this (except for name() and type() etc.) + virtual void initOrder + ( + PstreamBuffers&, + const primitivePatch& + ) const; + + //- Return new ordering for primitivePatch. + // Ordering is -faceMap: for every face + // index of the new face -rotation:for every new face the clockwise + // shift of the original face. Return false if nothing changes + // (faceMap is identity, rotation is 0), true otherwise. + virtual bool order + ( + PstreamBuffers&, + const primitivePatch&, + labelList& faceMap, + labelList& rotation + ) const; + + //- Write the polyPatch data as a dictionary + virtual void write(Ostream&) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C index 5f5fb220afe231762290b679da2cc0808a9c8cc5..7024c2ffbc59bd8364e487462daed5c3ac892138 100644 --- a/src/dynamicMesh/meshCut/directions/directions.C +++ b/src/dynamicMesh/meshCut/directions/directions.C @@ -32,6 +32,7 @@ License #include "meshTools.H" #include "hexMatcher.H" #include "Switch.H" +#include "globalMeshData.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -206,7 +207,7 @@ Foam::vectorField Foam::directions::propagateDirection mesh, changedFaces, changedFacesInfo, - mesh.nCells() + mesh.globalData().nTotalCells()+1 ); const List<directionInfo>& cellInfo = directionCalc.allCellInfo(); diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C index 9482861bded80ba1e13807b32a0696549cb2b30e..c9919351fc169acf46485d1f87e54354c17ae5aa 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C @@ -2332,7 +2332,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement seedFacesInfo.clear(); // Iterate until no change. Now 2:1 face difference should be satisfied - levelCalc.iterate(mesh_.globalData().nTotalFaces()); + levelCalc.iterate(mesh_.globalData().nTotalFaces()+1); // Now check point-connected cells (face-connected cells already ok): @@ -2836,7 +2836,7 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement2 seedFacesInfo, allFaceInfo, allCellInfo, - mesh_.globalData().nTotalCells() + mesh_.globalData().nTotalCells()+1 ); diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files index 21b90fa409a32cf2acb74d8e419c8c6cc4133d60..dc325c1ddf89e2bab343bdf7dda360096099a462 100644 --- a/src/finiteVolume/Make/files +++ b/src/finiteVolume/Make/files @@ -337,6 +337,7 @@ $(laplacianSchemes)/laplacianScheme/laplacianSchemes.C $(laplacianSchemes)/gaussLaplacianScheme/gaussLaplacianSchemes.C finiteVolume/fvc/fvcMeshPhi.C +finiteVolume/fvc/fvcSmooth/fvcSmooth.C general = cfdTools/general $(general)/findRefCell/findRefCell.C diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C index 0010ab4a0a0b40b4f2611c509954d9a9ed32f610..c3949a6ac0715bf305a1256a8a98ce1c617f843a 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.C @@ -171,6 +171,23 @@ tmp<Field<Type> > cyclicFvPatchField<Type>::patchNeighbourField() const } +template<class Type> +const cyclicFvPatchField<Type>& cyclicFvPatchField<Type>::neighbourPatchField() +const +{ + const GeometricField<Type, fvPatchField, volMesh>& fld = + static_cast<const GeometricField<Type, fvPatchField, volMesh>&> + ( + this->internalField() + ); + + return refCast<const cyclicFvPatchField<Type> > + ( + fld.boundaryField()[this->cyclicPatch().neighbPatchID()] + ); +} + + template<class Type> void cyclicFvPatchField<Type>::updateInterfaceMatrix ( diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H index 654bd0831a560a95b6b1e86d0d5f9c8ad53e99e6..198dbab74138aaaf5fc3897316c49cc154b8a970 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H @@ -150,9 +150,12 @@ public: // Evaluation functions - //- Return neighbour coupled given internal cell data + //- Return neighbour coupled internal cell data tmp<Field<Type> > patchNeighbourField() const; + //- Return reference to neighbour patchField + const cyclicFvPatchField<Type>& neighbourPatchField() const; + //- Update result field based on interface functionality virtual void updateInterfaceMatrix ( diff --git a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H index 73014e048a2a17b6336760eab447e39e0e7a80da..31c51c999d4e9139cec0f556a519dec309a83d79 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/derived/fan/fanFvPatchField.H @@ -140,7 +140,17 @@ public: //- Return the "jump" across the patch. virtual tmp<Field<Type> > jump() const { - return jump_; + if (this->cyclicPatch().owner()) + { + return jump_; + } + else + { + return refCast<const fanFvPatchField<Type> > + ( + this->neighbourPatchField() + ).jump(); + } } diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C new file mode 100644 index 0000000000000000000000000000000000000000..c04d16dcc9eb4233cd75790b4e1665db071f76d5 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.C @@ -0,0 +1,318 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "fvcSmooth.H" +#include "volFields.H" +#include "FaceCellWave.H" +#include "smoothData.H" +#include "sweepData.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +Foam::scalar Foam::smoothData::maxRatio = 1.0; + +void Foam::fvc::smooth +( + volScalarField& field, + const scalar coeff +) +{ + const fvMesh& mesh = field.mesh(); + smoothData::maxRatio = 1 + coeff; + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<smoothData> changedFacesInfo(changedFaces.size()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + // Check if owner value much larger than neighbour value or vice versa + if (field[own] > smoothData::maxRatio*field[nbr]) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + else if (field[nbr] > smoothData::maxRatio*field[own]) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[nbr])); + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Set initial field on cells + List<smoothData> cellData(mesh.nCells()); + + forAll(field, celli) + { + cellData[celli] = field[celli]; + } + + // Set initial field on faces + List<smoothData> faceData(mesh.nFaces()); + + + // Propagate information over whole domain + FaceCellWave<smoothData > smoothData + ( + mesh, + changedFaces, + changedFacesInfo, + faceData, + cellData, + mesh.globalData().nTotalCells() // max iterations + ); + + forAll(field, celli) + { + field[celli] = cellData[celli].value(); + } + + field.correctBoundaryConditions(); +} + + +void Foam::fvc::spread +( + volScalarField& field, + const volScalarField& alpha, + const label nLayers +) +{ + const fvMesh& mesh = field.mesh(); + smoothData::maxRatio = 1; + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<smoothData> changedFacesInfo(changedFaces.size()); + + // Set initial field on cells + List<smoothData> cellData(mesh.nCells()); + + forAll(field, celli) + { + cellData[celli] = field[celli]; + } + + // Set initial field on faces + List<smoothData> faceData(mesh.nFaces()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + if + ( + (alpha[own] > 0.01 && alpha[own] < 0.99) + || (alpha[nbr] > 0.01 && alpha[nbr] < 0.99) + ) + { + if (mag(alpha[own] - alpha[nbr]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + smoothData(max(field[own], field[nbr])) + ); + } + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + scalarField alphapn = + alpha.boundaryField()[patchi].patchNeighbourField(); + + if + ( + (alpha[own] > 0.01 && alpha[own] < 0.99) + || (alphapn[patchFacei] > 0.01 && alphapn[patchFacei] < 0.99) + ) + { + if (mag(alpha[own] - alphapn[patchFacei]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append(smoothData(field[own])); + } + } + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Propagate information over whole domain + FaceCellWave<smoothData> smoothData + ( + mesh, + faceData, + cellData + ); + + smoothData.setFaceInfo(changedFaces, changedFacesInfo); + + smoothData.iterate(nLayers); + + forAll(field, celli) + { + field[celli] = cellData[celli].value(); + } + + field.correctBoundaryConditions(); +} + + +void Foam::fvc::sweep +( + volScalarField& field, + const volScalarField& alpha, + const label nLayers +) +{ + const fvMesh& mesh = field.mesh(); + + DynamicList<label> changedFaces(mesh.nFaces()/100 + 100); + DynamicList<sweepData> changedFacesInfo(changedFaces.size()); + + // Set initial field on cells + List<sweepData> cellData(mesh.nCells()); + + // Set initial field on faces + List<sweepData> faceData(mesh.nFaces()); + + const unallocLabelList& owner = mesh.owner(); + const unallocLabelList& neighbour = mesh.neighbour(); + const vectorField& Cf = mesh.faceCentres(); + + forAll(owner, facei) + { + const label own = owner[facei]; + const label nbr = neighbour[facei]; + + if (mag(alpha[own] - alpha[nbr]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + sweepData(max(field[own], field[nbr]), + Cf[facei]) + ); + } + } + + // Insert all faces of coupled patches - FaceCellWave will correct them + forAll(mesh.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh.boundaryMesh()[patchi]; + + if (patch.coupled()) + { + forAll(patch, patchFacei) + { + label facei = patch.start() + patchFacei; + label own = mesh.faceOwner()[facei]; + + scalarField alphapn = + alpha.boundaryField()[patchi].patchNeighbourField(); + + if (mag(alpha[own] - alphapn[patchFacei]) > 0.2) + { + changedFaces.append(facei); + changedFacesInfo.append + ( + sweepData(field[own], Cf[facei]) + ); + } + } + } + } + + changedFaces.shrink(); + changedFacesInfo.shrink(); + + // Propagate information over whole domain + FaceCellWave<sweepData> sweepData + ( + mesh, + faceData, + cellData + ); + + sweepData.setFaceInfo(changedFaces, changedFacesInfo); + + sweepData.iterate(nLayers); + + forAll(field, celli) + { + if (cellData[celli].valid()) + { + field[celli] = max(field[celli], cellData[celli].value()); + } + } + + field.correctBoundaryConditions(); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H new file mode 100644 index 0000000000000000000000000000000000000000..ced147e3dc5c2b5703e84e4f48bafabb3536cb20 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/fvcSmooth.H @@ -0,0 +1,73 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +function + Foam::fvc::smooth + +Description + Function that uses smoothData to apply spatial smoothing of a + volume field using the FaceCellWave algorithm + +SourceFiles + fvcSmooth.C + +\*---------------------------------------------------------------------------*/ + +#ifndef fvcSmooth_H +#define fvcSmooth_H + +#include "volFieldsFwd.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace fvc +{ + void smooth + ( + volScalarField& field, + const scalar coeff + ); + + void spread + ( + volScalarField& field, + const volScalarField& alpha, + const label nLayers + ); + + void sweep + ( + volScalarField& field, + const volScalarField& alpha, + const label nLayers + ); +} +} + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H new file mode 100644 index 0000000000000000000000000000000000000000..b55c93a52a48528afded661ffb6753cda9f211e2 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothData.H @@ -0,0 +1,204 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +Class + Foam::smoothData + +Description + Helper class used by the smoothVolField class + +Files + smoothData.H + smoothDataI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef smoothData_H +#define smoothData_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class smoothData Declaration +\*---------------------------------------------------------------------------*/ + +class smoothData +{ + scalar value_; + + // Private Member Functions + + //- Update - gets information from neighbouring face/cell and + // uses this to update itself (if necessary) and return true + inline bool update + ( + const smoothData& svf, + const scalar scale, + const scalar tol + ); + + +public: + + // Static data members + + //- Field fraction + static scalar maxRatio; + + + // Constructors + + //- Construct null + inline smoothData(); + + //- Construct from inverse field value + inline smoothData(const scalar value); + + + // Member Functions + + // Access + + //- Return value + scalar value() const + { + return value_; + } + + + // Needed by FaceCellWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value + inline bool valid() const; + + //- Check for identical geometrical data + // Used for cyclics checking + inline bool sameGeometry + ( + const polyMesh&, + const smoothData&, + const scalar + ) const; + + //- Convert any absolute coordinates into relative to + // (patch)face centre + inline void leaveDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Reverse of leaveDomain + inline void enterDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Apply rotation matrix to any coordinates + inline void transform + ( + const polyMesh&, + const tensor& + ); + + //- Influence of neighbouring face + inline bool updateCell + ( + const polyMesh&, + const label thisCellI, + const label neighbourFaceI, + const smoothData& svf, + const scalar tol + ); + + //- Influence of neighbouring cell + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const label neighbourCellI, + const smoothData& svf, + const scalar tol + ); + + //- Influence of different value on same face + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const smoothData& svf, + const scalar tol + ); + + + // Member Operators + + inline void operator=(const scalar value); + + // Needed for List IO + inline bool operator==(const smoothData&) const; + + inline bool operator!=(const smoothData&) const; + + + // IOstream Operators + + friend Ostream& operator<< + ( + Ostream& os, + const smoothData& svf + ) + { + return os << svf.value_; + } + + friend Istream& operator>>(Istream& is, smoothData& svf) + { + return is >> svf.value_; + } +}; + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "smoothDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H new file mode 100644 index 0000000000000000000000000000000000000000..3faaa694b617e5fbea9a744a2b67afde6a9992af --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/smoothDataI.H @@ -0,0 +1,192 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline bool Foam::smoothData::update +( + const smoothData::smoothData& svf, + const scalar scale, + const scalar tol +) +{ + if (!valid() || (value_ < VSMALL)) + { + // My value not set - take over neighbour + value_ = svf.value()/scale; + + // Something changed - let caller know + return true; + } + else if (svf.value() > (1 + tol)*scale*value_) + { + // Neighbour is too big for me - Up my value + value_ = svf.value()/scale; + + // Something changed - let caller know + return true; + } + else + { + // Neighbour is not too big for me or change is too small + // Nothing changed + return false; + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::smoothData::smoothData() +: + value_(-GREAT) +{} + + +inline Foam::smoothData::smoothData(const scalar value) +: + value_(value) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline bool Foam::smoothData::valid() const +{ + return value_ > -SMALL; +} + + +inline bool Foam::smoothData::sameGeometry +( + const polyMesh&, + const smoothData&, + const scalar +) const +{ + return true; +} + + +inline void Foam::smoothData::leaveDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& +) +{} + + +inline void Foam::smoothData::transform +( + const polyMesh&, + const tensor& +) +{} + + +inline void Foam::smoothData::enterDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& +) +{} + + +inline bool Foam::smoothData::updateCell +( + const polyMesh&, + const label, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over info from face if more than deltaRatio larger + return update(svf, maxRatio, tol); +} + + +inline bool Foam::smoothData::updateFace +( + const polyMesh&, + const label, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over information from cell without any scaling (scale = 1.0) + return update(svf, 1.0, tol); +} + + +// Update this (face) with coupled face information. +inline bool Foam::smoothData::updateFace +( + const polyMesh&, + const label, + const smoothData& svf, + const scalar tol +) +{ + // Take over information from coupled face without any scaling (scale = 1.0) + return update(svf, 1.0, tol); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline void Foam::smoothData::operator= +( + const scalar value +) +{ + value_ = value; +} + + +inline bool Foam::smoothData::operator== +( + const smoothData& rhs +) const +{ + return value_ == rhs.value(); +} + + +inline bool Foam::smoothData::operator!= +( + const smoothData& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H new file mode 100644 index 0000000000000000000000000000000000000000..ebee81c0f751db3cdfa1be56155a3cec2efea198 --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepData.H @@ -0,0 +1,205 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +Class + Foam::sweepData + +Description + Helper class used by fvcSmooth + +Files + sweepData.H + sweepDataI.H + +\*---------------------------------------------------------------------------*/ + +#ifndef sweepData_H +#define sweepData_H + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class sweepData Declaration +\*---------------------------------------------------------------------------*/ + +class sweepData +{ + scalar value_; + point origin_; + + // Private Member Functions + + //- Update - gets information from neighbouring face/cell and + // uses this to update itself (if necessary) and return true + inline bool update + ( + const sweepData& svf, + const point& position, + const scalar tol + ); + + +public: + + // Constructors + + //- Construct null + inline sweepData(); + + //- Construct from component + inline sweepData(const scalar value, const point& origin); + + + // Member Functions + + // Access + + //- Return value + scalar value() const + { + return value_; + } + + //- Return origin + const point& origin() const + { + return origin_; + } + + + // Needed by FaceCellWave + + //- Check whether origin has been changed at all or + // still contains original (invalid) value + inline bool valid() const; + + //- Check for identical geometrical data + // Used for cyclics checking + inline bool sameGeometry + ( + const polyMesh&, + const sweepData&, + const scalar + ) const; + + //- Convert any absolute coordinates into relative to + // (patch)face centre + inline void leaveDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Reverse of leaveDomain + inline void enterDomain + ( + const polyMesh&, + const polyPatch&, + const label patchFaceI, + const point& faceCentre + ); + + //- Apply rotation matrix to any coordinates + inline void transform + ( + const polyMesh&, + const tensor& + ); + + //- Influence of neighbouring face + inline bool updateCell + ( + const polyMesh&, + const label thisCellI, + const label neighbourFaceI, + const sweepData& svf, + const scalar tol + ); + + //- Influence of neighbouring cell + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const label neighbourCellI, + const sweepData& svf, + const scalar tol + ); + + //- Influence of different value on same face + inline bool updateFace + ( + const polyMesh&, + const label thisFaceI, + const sweepData& svf, + const scalar tol + ); + + + // Member Operators + + inline void operator=(const scalar value); + + // Needed for List IO + inline bool operator==(const sweepData&) const; + + inline bool operator!=(const sweepData&) const; + + + // IOstream Operators + + friend Ostream& operator<< + ( + Ostream& os, + const sweepData& svf + ) + { + return os << svf.value_ << svf.origin_; + } + + friend Istream& operator>>(Istream& is, sweepData& svf) + { + return is >> svf.value_ >> svf.origin_; + } +}; + + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#include "sweepDataI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H new file mode 100644 index 0000000000000000000000000000000000000000..651c5e4e5260252251da7cd445c056d419d8965c --- /dev/null +++ b/src/finiteVolume/finiteVolume/fvc/fvcSmooth/sweepDataI.H @@ -0,0 +1,208 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-2010 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 <http://www.gnu.org/licenses/>. + +\*---------------------------------------------------------------------------*/ + +#include "transform.H" + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +inline bool Foam::sweepData::update +( + const sweepData::sweepData& svf, + const point& position, + const scalar tol +) +{ + if (!valid()) + { + operator=(svf); + return true; + } + + scalar myDist2 = magSqr(position - origin()); + + if (myDist2 < SMALL) + { + if (svf.value() > value()) + { + operator=(svf); + return true; + } + else + { + return false; + } + } + + scalar dist2 = magSqr(position - svf.origin()); + + if (dist2 < myDist2) + { + operator=(svf); + return true; + } + + return false; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::sweepData::sweepData() +: + value_(-GREAT), + origin_(vector::max) +{} + + +inline Foam::sweepData::sweepData(const scalar value, const point& origin) +: + value_(value), + origin_(origin) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline bool Foam::sweepData::valid() const +{ + return value_ > -SMALL; +} + + +inline bool Foam::sweepData::sameGeometry +( + const polyMesh&, + const sweepData&, + const scalar +) const +{ + return true; +} + + +inline void Foam::sweepData::leaveDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& faceCentre +) +{ + origin_ -= faceCentre; +} + + +inline void Foam::sweepData::transform +( + const polyMesh&, + const tensor& rotTensor +) +{ + origin_ = Foam::transform(rotTensor, origin_); +} + + +inline void Foam::sweepData::enterDomain +( + const polyMesh&, + const polyPatch&, + const label, + const point& faceCentre +) +{ + // back to absolute form + origin_ += faceCentre; +} + + +inline bool Foam::sweepData::updateCell +( + const polyMesh& mesh, + const label thisCellI, + const label, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.cellCentres()[thisCellI], tol); +} + + +inline bool Foam::sweepData::updateFace +( + const polyMesh& mesh, + const label thisFaceI, + const label, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.faceCentres()[thisFaceI], tol); +} + + +// Update this (face) with coupled face information. +inline bool Foam::sweepData::updateFace +( + const polyMesh& mesh, + const label thisFaceI, + const sweepData& svf, + const scalar tol +) +{ + return update(svf, mesh.faceCentres()[thisFaceI], tol); +} + + +// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // + +inline void Foam::sweepData::operator= +( + const scalar value +) +{ + value_ = value; +} + + +inline bool Foam::sweepData::operator== +( + const sweepData& rhs +) const +{ + return origin() == rhs.origin(); +} + + +inline bool Foam::sweepData::operator!= +( + const sweepData& rhs +) const +{ + return !(*this == rhs); +} + + +// ************************************************************************* // diff --git a/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C b/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C index 72531d5cfd12427b5d43ecfeb90db1009a85b011..02db9ad941ab5838d4e031bba1cbcca56ab45e7c 100644 --- a/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C +++ b/src/fvMotionSolver/motionDiffusivity/inverseFaceDistance/inverseFaceDistanceDiffusivity.C @@ -115,7 +115,7 @@ void Foam::inverseFaceDistanceDiffusivity::correct() mesh, changedFaces, faceDist, - mesh.globalData().nTotalCells() // max iterations + mesh.globalData().nTotalCells()+1 // max iterations ); const List<wallPoint>& faceInfo = waveInfo.allFaceInfo(); diff --git a/src/meshTools/cellClassification/cellClassification.C b/src/meshTools/cellClassification/cellClassification.C index a7236654777262149c8ddff3918ee314f188ddd2..daf1de6ac0c03b6b1d9f7d505bc7ca2f59ec758e 100644 --- a/src/meshTools/cellClassification/cellClassification.C +++ b/src/meshTools/cellClassification/cellClassification.C @@ -322,7 +322,7 @@ void Foam::cellClassification::markCells changedFaces, // Labels of changed faces changedFacesInfo, // Information on changed faces cellInfoList, // Information on all cells - mesh_.globalData().nTotalCells() // max iterations + mesh_.globalData().nTotalCells()+1 // max iterations ); // Get information out of cellInfoList diff --git a/src/meshTools/cellDist/patchWave/patchDataWave.C b/src/meshTools/cellDist/patchWave/patchDataWave.C index bb6da1750dd6c518d8c873d9c6f48a258cb1be8b..988b8a0be5fdef061e978f98ef4f79e2e6193f2d 100644 --- a/src/meshTools/cellDist/patchWave/patchDataWave.C +++ b/src/meshTools/cellDist/patchWave/patchDataWave.C @@ -225,7 +225,7 @@ void Foam::patchDataWave<TransferType>::correct() mesh(), changedFaces, faceDist, - mesh().globalData().nTotalCells() // max iterations + mesh().globalData().nTotalCells()+1 // max iterations ); diff --git a/src/meshTools/cellDist/patchWave/patchWave.C b/src/meshTools/cellDist/patchWave/patchWave.C index 85118da58dbccc891d42dcb39048c6206772f82c..826430b00324a9788da1f71a34ee1784d0d40320 100644 --- a/src/meshTools/cellDist/patchWave/patchWave.C +++ b/src/meshTools/cellDist/patchWave/patchWave.C @@ -186,7 +186,7 @@ void Foam::patchWave::correct() mesh(), changedFaces, faceDist, - mesh().globalData().nTotalCells() // max iterations + mesh().globalData().nTotalCells()+1 // max iterations ); diff --git a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C index 32f86c82c83f6c08afcf840abdee572d22c06be6..04338dc46b9af2a45633cc0682ddbb1940b22110 100644 --- a/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C +++ b/src/parallel/decompose/decompositionMethods/structuredDecomp/structuredDecomp.C @@ -148,7 +148,7 @@ Foam::labelList Foam::structuredDecomp::decompose patchData, faceData, cellData, - mesh.globalData().nTotalCells() + mesh.globalData().nTotalCells()+1 ); // And extract diff --git a/src/parallel/decompose/ptscotchDecomp/Make/options b/src/parallel/decompose/ptscotchDecomp/Make/options index 488f5ef3bb46e3f4d6442264ce571cca8f4bc0f6..f3fa4b92b5b0a79978557f537a9519b87172d5ce 100644 --- a/src/parallel/decompose/ptscotchDecomp/Make/options +++ b/src/parallel/decompose/ptscotchDecomp/Make/options @@ -3,7 +3,7 @@ sinclude $(RULES)/mplib$(WM_MPLIB) EXE_INC = \ $(PFLAGS) $(PINC) \ - -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \ + -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \ -I/usr/include/scotch \ -I../decompositionMethods/lnInclude diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C index 98f231cf4858ce5b40e2c3a443e23ca848641783..ff262c5645d28388a71529f20f6c0a2dfaceed91 100644 --- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C @@ -180,17 +180,31 @@ Foam::label Foam::ptscotchDecomp::decompose // Info<< "Dumping Scotch graph file to " << str.name() << endl // << "Use this in combination with gpart." << endl; // -// label version = 0; +// // Distributed graph file (.grf) +// label version = 2; // str << version << nl; -// // Numer of vertices -// str << xadj.size()-1 << ' ' << adjncy.size() << nl; +// // Number of files + +// // Number of files (procglbnbr) +// str << Pstream::nProcs(); +// // My file number (procloc) +// str << ' ' << Pstream::myProcNo() << nl; +// +// // Total number of vertices (vertglbnbr) +// str << returnReduce(mesh.nCells(), sumOp<label>()); +// // Total number of connections (edgeglbnbr) +// str << ' ' << returnReduce(xadj[mesh.nCells()], sumOp<label>()) +// << nl; +// // Local number of vertices (vertlocnbr) +// str << mesh.nCells(); +// // Local number of connections (edgelocnbr) +// str << ' ' << xadj[mesh.nCells()] << nl; // // Numbering starts from 0 // label baseval = 0; -// // Has weights? -// label hasEdgeWeights = 0; -// label hasVertexWeights = 0; -// label numericflag = 10*hasEdgeWeights+hasVertexWeights; -// str << baseval << ' ' << numericflag << nl; +// // Start of my global vertices (procdsptab[proclocnum]) +// str << ' ' << globalCells.toGlobal(0); +// 100*hasVertlabels+10*hasEdgeWeights+1*hasVertWeighs +// str << ' ' << "0" << nl; // for (label cellI = 0; cellI < xadj.size()-1; cellI++) // { // label start = xadj[cellI]; diff --git a/src/parallel/decompose/scotchDecomp/Make/options b/src/parallel/decompose/scotchDecomp/Make/options index 2f01188ac6f6b08d684f4a3f9d4c953857d04c5f..10374b63c320be020896856f23a8cf21a2746e3a 100644 --- a/src/parallel/decompose/scotchDecomp/Make/options +++ b/src/parallel/decompose/scotchDecomp/Make/options @@ -1,5 +1,5 @@ EXE_INC = \ - -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \ + -I$(WM_THIRD_PARTY_DIR)/scotch_5.1.10a/include \ -I/usr/include/scotch \ -I../decompositionMethods/lnInclude diff --git a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C index af250a25f7b2048096672914ba92c18947379643..d97395d864ceaefb043f2c0447d7038af17c7f5f 100644 --- a/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C +++ b/src/turbulenceModels/LES/LESdeltas/smoothDelta/smoothDelta.C @@ -137,7 +137,7 @@ void smoothDelta::calcDelta() changedFacesInfo, faceDeltaData, cellDeltaData, - mesh_.globalData().nTotalCells() // max iterations + mesh_.globalData().nTotalCells()+1 // max iterations ); forAll(delta_, cellI) diff --git a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C index 486eccaa9813ab568ffe96912e42e52269507e37..4a76228d23f2c368a1b392510825f9fe69d18950 100644 --- a/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C +++ b/src/turbulenceModels/incompressible/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C @@ -172,7 +172,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe ( IOobject ( - "epsilon", + "nut", runTime_.timeName(), mesh_, IOobject::NO_READ, diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U index 51ef892fb8a5263548b1471887b3702fc8e41b88..f1006ceb5ed5db11c66679186c6c98a424d8190e 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/U @@ -15,47 +15,47 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 1 -1 0 0 0 0 ]; +dimensions [0 1 -1 0 0 0 0]; -internalField uniform ( 0 0 0 ); +internalField uniform (0 0 0); boundaryField { inlet { type pressureInletOutletVelocity; - value uniform ( 0 0 0 ); + value uniform (0 0 0); } outlet1 { type inletOutlet; - inletValue uniform ( 0 0 0 ); - value uniform ( 0 0 0 ); + inletValue uniform (0 0 0); + value uniform (0 0 0); } outlet2 { type inletOutlet; - inletValue uniform ( 0 0 0 ); - value uniform ( 0 0 0 ); + inletValue uniform (0 0 0); + value uniform (0 0 0); } baffles { type fixedValue; - value uniform ( 0 0 0 ); + value uniform (0 0 0); } fan_half0 { type cyclic; } - defaultFaces - { - type fixedValue; - value uniform ( 0 0 0 ); - } fan_half1 { type cyclic; } + defaultFaces + { + type fixedValue; + value uniform (0 0 0); + } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon index dc1eb9da3ec01e6e20b9eb1bfee53d5a810eeb34..5c6c53664db8773322bca5f71be4360714228dda 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/epsilon @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 2 -3 0 0 0 0 ]; +dimensions [0 2 -3 0 0 0 0]; internalField uniform 1; @@ -42,21 +42,27 @@ boundaryField baffles { type epsilonWallFunction; - value uniform 1; + Cmu 0.09; + kappa 0.41; + E 9.8; + value uniform 0; } fan_half0 { type cyclic; } + fan_half1 + { + type cyclic; + } defaultFaces { type epsilonWallFunction; + Cmu 0.09; + kappa 0.41; + E 9.8; value uniform 1; } - fan_half1 - { - type cyclic; - } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k index 06aeff0f1f5ab3f291870fd35e0fe6dea2557c93..18cb53d51920f57d6f114259326f42644171e03f 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/k @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 2 -2 0 0 0 0 ]; +dimensions [0 2 -2 0 0 0 0]; internalField uniform 1; @@ -42,21 +42,21 @@ boundaryField baffles { type kqRWallFunction; - value uniform 1; + value uniform 0; } fan_half0 { type cyclic; } + fan_half1 + { + type cyclic; + } defaultFaces { type kqRWallFunction; value uniform 1; } - fan_half1 - { - type cyclic; - } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda index 699a951041d71878bd5f2ac7a250c4eaf8ce4c36..ffc7ac906f4c311a7a4bf9ee078e66a7686c517d 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nuTilda @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 2 -1 0 0 0 0 ]; +dimensions [0 2 -1 0 0 0 0]; internalField uniform 0; @@ -41,14 +41,14 @@ boundaryField { type cyclic; } - defaultFaces - { - type zeroGradient; - } fan_half1 { type cyclic; } + defaultFaces + { + type zeroGradient; + } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut index 1a5afb9689fbb662c24a1cce95bdbfe3ba166031..6ed8f42dfc22716594c16f28ca52a8d0d38c29fb 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/nut @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 2 -1 0 0 0 0 ]; +dimensions [0 2 -1 0 0 0 0]; internalField uniform 0; @@ -39,21 +39,27 @@ boundaryField baffles { type nutkWallFunction; + Cmu 0.09; + kappa 0.41; + E 9.8; value uniform 0; } fan_half0 { type cyclic; } + fan_half1 + { + type cyclic; + } defaultFaces { type nutkWallFunction; + Cmu 0.09; + kappa 0.41; + E 9.8; value uniform 0; } - fan_half1 - { - type cyclic; - } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p index c77ad857259c08270564e1e72249631d59902d89..66f5c1d712ae70f900e35160f3dc300c688ebb67 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/0/p @@ -15,7 +15,7 @@ FoamFile } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -dimensions [ 0 2 -2 0 0 0 0 ]; +dimensions [0 2 -2 0 0 0 0]; internalField uniform 100000; @@ -24,14 +24,12 @@ boundaryField inlet { type timeVaryingTotalPressure; - p0 100040; - outOfBounds clamp; - fileName "$FOAM_CASE/constant/p0vsTime"; - U U; - phi phi; rho none; psi none; gamma 1; + p0 100040; + fileName "$FOAM_CASE/constant/p0vsTime"; + outOfBounds clamp; value uniform 100040; } outlet1 @@ -52,19 +50,19 @@ boundaryField { type fan; patchType cyclic; - f 2 ( 50 -0.1 ); - value $internalField; - } - defaultFaces - { - type zeroGradient; + f 2(100 -0.1); + value uniform 0; } fan_half1 { type fan; patchType cyclic; - f 2 ( 50 -0.1 ); - value $internalField; + f 2(100 -0.1); + value uniform 0; + } + defaultFaces + { + type zeroGradient; } } diff --git a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict index 9dc31d5561fb3c86c93e949d18a3dd5b39d41977..63b26690081d708fb88fc371777da73e5f6f5a7e 100644 --- a/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict +++ b/tutorials/incompressible/pimpleFoam/t-junction-with-fan/system/controlDict @@ -28,7 +28,6 @@ endTime 1; deltaT 0.001; writeControl adjustableRunTime; - writeInterval 0.1; purgeWrite 0; @@ -47,7 +46,7 @@ runTimeModifiable true; adjustTimeStep yes; -maxCo 5; +maxCo 3; libs (