From ffe930843290f69f3d8227921db5e0c4862546f7 Mon Sep 17 00:00:00 2001 From: henry Date: Wed, 3 Sep 2008 08:34:04 +0100 Subject: [PATCH 01/39] Minor updates. --- .../multiphase/interDyMFoam/Make/options | 2 ++ .../solvers/multiphase/interDyMFoam/UEqn.H | 34 ------------------ .../multiphase/interDyMFoam/gammaEqn.H | 35 ------------------- .../interDyMFoam/gammaEqnSubCycle.H | 35 ------------------- .../multiphase/interDyMFoam/interDyMFoam.C | 1 - 5 files changed, 2 insertions(+), 105 deletions(-) delete mode 100644 applications/solvers/multiphase/interDyMFoam/UEqn.H delete mode 100644 applications/solvers/multiphase/interDyMFoam/gammaEqn.H delete mode 100644 applications/solvers/multiphase/interDyMFoam/gammaEqnSubCycle.H diff --git a/applications/solvers/multiphase/interDyMFoam/Make/options b/applications/solvers/multiphase/interDyMFoam/Make/options index 24d71bc7ee..6755d5bc59 100644 --- a/applications/solvers/multiphase/interDyMFoam/Make/options +++ b/applications/solvers/multiphase/interDyMFoam/Make/options @@ -1,4 +1,6 @@ EXE_INC = \ + -I../rasInterFoam \ + -I../interFoam \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/lnInclude \ -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \ diff --git a/applications/solvers/multiphase/interDyMFoam/UEqn.H b/applications/solvers/multiphase/interDyMFoam/UEqn.H deleted file mode 100644 index 4c14afe1a6..0000000000 --- a/applications/solvers/multiphase/interDyMFoam/UEqn.H +++ /dev/null @@ -1,34 +0,0 @@ -surfaceScalarField muEff -( - "muEff", - twoPhaseProperties.muf() - + fvc::interpolate(rho*turbulence->nut()) -); - -fvVectorMatrix UEqn -( - fvm::ddt(rho, U) - + fvm::div(rhoPhi, U) - - fvm::laplacian(muEff, U) - - (fvc::grad(U) & fvc::grad(muEff)) -//- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf())) -); - -UEqn.relax(); - -if (momentumPredictor) -{ - solve - ( - UEqn - == - fvc::reconstruct - ( - ( - fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma) - - ghf*fvc::snGrad(rho) - - fvc::snGrad(pd) - )*mesh.magSf() - ) - ); -} diff --git a/applications/solvers/multiphase/interDyMFoam/gammaEqn.H b/applications/solvers/multiphase/interDyMFoam/gammaEqn.H deleted file mode 100644 index 8978d1d293..0000000000 --- a/applications/solvers/multiphase/interDyMFoam/gammaEqn.H +++ /dev/null @@ -1,35 +0,0 @@ -{ - word gammaScheme("div(phi,gamma)"); - word gammarScheme("div(phirb,gamma)"); - - surfaceScalarField phic = mag(phi/mesh.magSf()); - phic = min(interface.cGamma()*phic, max(phic)); - surfaceScalarField phir = phic*interface.nHatf(); - - for (int gCorr=0; gCorr 1) -{ - dimensionedScalar totalDeltaT = runTime.deltaT(); - surfaceScalarField rhoPhiSum = 0.0*rhoPhi; - - for - ( - subCycle gammaSubCycle(gamma, nGammaSubCycles); - !(++gammaSubCycle).end(); - ) - { -# include "gammaEqn.H" - rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi; - } - - rhoPhi = rhoPhiSum; -} -else -{ -# include "gammaEqn.H" -} - -interface.correct(); - -rho == gamma*rho1 + (scalar(1) - gamma)*rho2; diff --git a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C index a126edaa92..915fee9a97 100644 --- a/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C +++ b/applications/solvers/multiphase/interDyMFoam/interDyMFoam.C @@ -41,7 +41,6 @@ Description #include "twoPhaseMixture.H" #include "incompressible/RASModel/RASModel.H" #include "probes.H" -#include "EulerDdtScheme.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // From 485a2549c821bd3b2261b8b2873a815fb2049f58 Mon Sep 17 00:00:00 2001 From: henry Date: Wed, 3 Sep 2008 08:35:44 +0100 Subject: [PATCH 02/39] Removed spaces. --- .../constant/polyMesh/blockMeshDict | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tutorials/boundaryFoam/boundaryWallFunctions/constant/polyMesh/blockMeshDict b/tutorials/boundaryFoam/boundaryWallFunctions/constant/polyMesh/blockMeshDict index 56b7dbe749..d5201ded84 100644 --- a/tutorials/boundaryFoam/boundaryWallFunctions/constant/polyMesh/blockMeshDict +++ b/tutorials/boundaryFoam/boundaryWallFunctions/constant/polyMesh/blockMeshDict @@ -16,7 +16,7 @@ FoamFile convertToMeters 0.05; -vertices +vertices ( (0 -1 0) (0 0 0) @@ -32,24 +32,24 @@ vertices (0.1 1 0.1) ); -blocks +blocks ( hex (0 3 4 1 6 9 10 7) (1 40 1) simpleGrading (1 1 1) hex (1 4 5 2 7 10 11 8) (1 40 1) simpleGrading (1 1 1) ); -edges +edges ( ); -patches +patches ( - wall lowerWall + wall lowerWall ( (0 3 9 6) ) - wall upperWall + wall upperWall ( (2 8 11 5) ) From 24fcef879db65db8302ecc801feb87a0691afb6e Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2008 10:55:04 +0100 Subject: [PATCH 03/39] empty surface fields --- .../postProcessing/dataConversion/foamToVTK/writeSurfFields.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C index 5e40ed5871..58107bc776 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C +++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/writeSurfFields.C @@ -28,7 +28,7 @@ License #include "OFstream.H" #include "floatScalar.H" #include "writeFuns.H" -#include "emptyFvPatchFields.H" +#include "emptyFvsPatchFields.H" #include "fvsPatchFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -100,7 +100,7 @@ void writeSurfFields const fvPatch& pp = mesh.boundary()[patchI]; - if (isA(pf)) + if (isA(pf)) { // Note: loop over polypatch size, not fvpatch size. forAll(pp.patch(), i) From ec6ba9c15bec00a508b6d8bb9116570f7f03ea0f Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2008 10:55:54 +0100 Subject: [PATCH 04/39] paraview .so library locations --- etc/apps/paraview3/bashrc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/etc/apps/paraview3/bashrc b/etc/apps/paraview3/bashrc index 56cdd8aacf..b8da8ba698 100644 --- a/etc/apps/paraview3/bashrc +++ b/etc/apps/paraview3/bashrc @@ -45,9 +45,9 @@ export ParaView_INST_DIR=$WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION export ParaView_DIR=$ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER if [ "$PYTHONPATH" ]; then - export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping + export PYTHONPATH=$PYTHONPATH:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3 else - export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping + export PYTHONPATH=$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3 fi From fe72c1c788bbfe2b0beb4edeb022a5631cea67ca Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2008 10:55:58 +0100 Subject: [PATCH 05/39] paraview .so library locations --- etc/apps/paraview3/cshrc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/etc/apps/paraview3/cshrc b/etc/apps/paraview3/cshrc index c9fd98f4a2..636b10eb75 100644 --- a/etc/apps/paraview3/cshrc +++ b/etc/apps/paraview3/cshrc @@ -45,9 +45,9 @@ setenv ParaView_INST_DIR $WM_THIRD_PARTY_DIR/ParaView$ParaView_VERSION setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER if ($?PYTHONPATH) then - setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping + setenv PYTHONPATH ${PYTHONPATH}:$ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3 else - setenv PYTHONPATH $ParaView_DIR/bin:$ParaView_DIR/Utilities/VTKPythonWrapping + setenv PYTHONPATH $ParaView_DIR/Utilities/VTKPythonWrapping:$ParaView_DIR/lib/paraview-3.3 endif if ( -r $ParaView_INST_DIR ) then From a8a3ea5475b3d8d6e141083b31df02bb1eaae10d Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2008 10:56:24 +0100 Subject: [PATCH 06/39] binary output --- src/OpenFOAM/meshes/boundBox/boundBox.C | 42 +++++++++++++++++++++++-- src/OpenFOAM/meshes/boundBox/boundBox.H | 23 ++++++++++++-- 2 files changed, 61 insertions(+), 4 deletions(-) diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C index c368e1a662..c87e552c85 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.C +++ b/src/OpenFOAM/meshes/boundBox/boundBox.C @@ -87,9 +87,47 @@ boundBox::boundBox(Istream& is) // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * // -Ostream& operator<<(Ostream& os, const boundBox& b) +Ostream& operator<<(Ostream& os, const boundBox& bb) { - return os << b.min() << token::SPACE << b.max(); + if (os.format() == IOstream::ASCII) + { + os << bb.min_ << token::SPACE << bb.max_; + } + else + { + os.write + ( + reinterpret_cast(&bb.min_), + sizeof(boundBox) + ); + } + + // Check state of Ostream + os.check("Ostream& operator<<(Ostream&, const boundBox&)"); + + return os; +} + + +Istream& operator>>(Istream& is, boundBox& bb) +{ + if (is.format() == IOstream::ASCII) + { + return is >> bb.min_ >> bb.max_; + } + else + { + is.read + ( + reinterpret_cast(&bb.min_), + sizeof(boundBox) + ); + } + + // Check state of Istream + is.check("Istream& operator>>(Istream&, boundBox&)"); + + return is; } diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H index 40287e53eb..1493c0fe2f 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.H +++ b/src/OpenFOAM/meshes/boundBox/boundBox.H @@ -153,12 +153,31 @@ public: } - // Ostream operator + // Friend Operators - friend Ostream& operator<<(Ostream& os, const boundBox& b); + friend bool operator==(const boundBox& a, const boundBox& b) + { + return (a.min_ == b.min_) && (a.max_ == b.max_); + } + + friend bool operator!=(const boundBox& a, const boundBox& b) + { + return !(a == b); + } + + + // IOstream operator + + friend Istream& operator>>(Istream& is, boundBox&); + friend Ostream& operator<<(Ostream& os, const boundBox&); }; +//- Specify data associated with boundBox type is contiguous +template<> +inline bool contiguous() {return contiguous();} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam From 5726847ddadab7bce3e0e3d0524649e0e9602528 Mon Sep 17 00:00:00 2001 From: mattijs Date: Wed, 3 Sep 2008 10:59:50 +0100 Subject: [PATCH 07/39] moved oriented surface out of refinementSurfaces --- .../autoHexMesh/shellSurfaces/shellSurfaces.C | 23 ++++++++++++++++++- 1 file changed, 22 insertions(+), 1 deletion(-) diff --git a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C index fe64ea81c3..1a577acecb 100644 --- a/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C +++ b/src/autoMesh/autoHexMesh/shellSurfaces/shellSurfaces.C @@ -30,6 +30,7 @@ License #include "triSurfaceMesh.H" #include "refinementSurfaces.H" #include "searchableSurfaces.H" +#include "orientedSurface.H" #include "pointIndexHit.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -211,7 +212,27 @@ void Foam::shellSurfaces::orient() refCast(s) ); - refinementSurfaces::orientSurface(outsidePt, shell); + // Flip surface so outsidePt is outside. + bool anyFlipped = orientedSurface::orient + ( + shell, + outsidePt, + true + ); + + if (anyFlipped) + { + // orientedSurface will have done a clearOut of the surface. + // we could do a clearout of the triSurfaceMeshes::trees() + // but these aren't affected by orientation + // (except for cached + // sideness which should not be set at this point. + // !!Should check!) + + Info<< "shellSurfaces : Flipped orientation of surface " + << s.name() + << " so point " << outsidePt << " is outside." << endl; + } } } } From 537b6074235af0a1299f8734344a601ee69e3da5 Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 4 Sep 2008 20:00:33 +0100 Subject: [PATCH 08/39] Fixed bug in ddtPhiCorr for compressible flow. --- .../CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C index dbbda95fcf..f591376bf5 100644 --- a/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C +++ b/src/finiteVolume/finiteVolume/ddtSchemes/CrankNicholsonDdtScheme/CrankNicholsonDdtScheme.C @@ -978,14 +978,14 @@ CrankNicholsonDdtScheme::fvcDdtPhiCorr ddt0_ > ( "ddt0(" + U.name() + ')', - rho.dimensions()*U.dimensions() + U.dimensions() ); DDt0Field& dphidt0 = ddt0_ ( "ddt0(" + phi.name() + ')', - phi.dimensions() + U.dimensions()*dimArea ); IOobject ddtIOobject From c2d3722bbded3d5afb2fd89421a3072cb07a1d57 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 4 Sep 2008 21:48:05 +0100 Subject: [PATCH 09/39] primitiveMesh instead of polyMesh --- src/meshTools/indexedOctree/treeDataCell.C | 6 +++--- src/meshTools/indexedOctree/treeDataCell.H | 10 +++++----- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/meshTools/indexedOctree/treeDataCell.C b/src/meshTools/indexedOctree/treeDataCell.C index cb41899306..131f0c69be 100644 --- a/src/meshTools/indexedOctree/treeDataCell.C +++ b/src/meshTools/indexedOctree/treeDataCell.C @@ -26,7 +26,7 @@ License #include "treeDataCell.H" #include "indexedOctree.H" -#include "polyMesh.H" +#include "primitiveMesh.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -71,7 +71,7 @@ Foam::treeBoundBox Foam::treeDataCell::calcCellBb(const label cellI) const Foam::treeDataCell::treeDataCell ( const bool cacheBb, - const polyMesh& mesh, + const primitiveMesh& mesh, const labelList& cellLabels ) : @@ -94,7 +94,7 @@ Foam::treeDataCell::treeDataCell Foam::treeDataCell::treeDataCell ( const bool cacheBb, - const polyMesh& mesh + const primitiveMesh& mesh ) : mesh_(mesh), diff --git a/src/meshTools/indexedOctree/treeDataCell.H b/src/meshTools/indexedOctree/treeDataCell.H index 1cb15db6fe..6fe5d86755 100644 --- a/src/meshTools/indexedOctree/treeDataCell.H +++ b/src/meshTools/indexedOctree/treeDataCell.H @@ -47,7 +47,7 @@ namespace Foam { // Forward declaration of classes -class polyMesh; +class primitiveMesh; template class indexedOctree; /*---------------------------------------------------------------------------*\ @@ -58,7 +58,7 @@ class treeDataCell { // Private data - const polyMesh& mesh_; + const primitiveMesh& mesh_; //- Subset of cells to work on const labelList cellLabels_; @@ -87,12 +87,12 @@ public: treeDataCell ( const bool cacheBb, - const polyMesh&, + const primitiveMesh&, const labelList& ); //- Construct from mesh. Uses all cells in mesh. - treeDataCell(const bool cacheBb, const polyMesh&); + treeDataCell(const bool cacheBb, const primitiveMesh&); // Member Functions @@ -104,7 +104,7 @@ public: return cellLabels_; } - const polyMesh& mesh() const + const primitiveMesh& mesh() const { return mesh_; } From fb99c7b63e6ecb3a090ecf9ef345560029225772 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 4 Sep 2008 21:48:34 +0100 Subject: [PATCH 10/39] indentation --- src/meshTools/Make/options | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options index 3e8506b54f..ef1033b0e8 100644 --- a/src/meshTools/Make/options +++ b/src/meshTools/Make/options @@ -3,5 +3,5 @@ EXE_INC = \ -I$(LIB_SRC)/lagrangian/basic/lnInclude LIB_LIBS = \ - -ltriSurface \ - -llagrangian + -ltriSurface \ + -llagrangian From 4773de4695a45cd81eeefaee5f45ff3b0718c4e2 Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 4 Sep 2008 21:49:11 +0100 Subject: [PATCH 11/39] \r handling --- .../mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L index f0c8e0a6c4..a8cb5dbb98 100644 --- a/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L +++ b/applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L @@ -697,7 +697,7 @@ endOfSection {space}")"{space} /* ------ Ignore remaining space and \n s. ------ */ -<*>{some_space}|\n { +<*>{some_space}|\n|\r { } From 987d716acf94539b2f1cc364acc0e07da34e9527 Mon Sep 17 00:00:00 2001 From: henry Date: Thu, 4 Sep 2008 22:39:00 +0100 Subject: [PATCH 12/39] Removed check for single component mass-fraction. --- .../mixtures/combustionMixture/combustionMixture.C | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/src/thermophysicalModels/combustion/mixtureThermos/mixtures/combustionMixture/combustionMixture.C b/src/thermophysicalModels/combustion/mixtureThermos/mixtures/combustionMixture/combustionMixture.C index 47a7cd3fbb..0b837ccd20 100644 --- a/src/thermophysicalModels/combustion/mixtureThermos/mixtures/combustionMixture/combustionMixture.C +++ b/src/thermophysicalModels/combustion/mixtureThermos/mixtures/combustionMixture/combustionMixture.C @@ -107,18 +107,6 @@ combustionMixture::combustionMixture ); } } - - // Error check for single component mixtures - should be uniform 1 - if (Y_.size() == 1) - { - if (average(Y_[0]).value() < 0.999) - { - FatalErrorIn("combustionMixture::combustionMixture") - << "Mass fraction for single component mixture must be unity" - << nl << "Please correct field: " << species_[0] - << " (Ydefault)" << nl << abort(FatalError); - } - } } From 1a9ec8954433585046b33c79b1e31b3007672a67 Mon Sep 17 00:00:00 2001 From: henry Date: Fri, 5 Sep 2008 23:28:27 +0100 Subject: [PATCH 13/39] Added standard High Re k-omega turbulence model. --- .../RAS/incompressible/Make/files | 1 + .../RAS/incompressible/kOmega/kOmega.C | 268 ++++++++++++++++++ .../RAS/incompressible/kOmega/kOmega.H | 205 ++++++++++++++ .../kOmega/kOmegaWallFunctionsI.H | 125 ++++++++ .../kOmega/kOmegaWallViscosityI.H | 70 +++++ .../RAS/incompressible/kOmega/wallOmegaI.H | 51 ++++ .../RAS/incompressible/kOmegaSST/kOmegaSST.H | 1 + 7 files changed, 721 insertions(+) create mode 100644 src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C create mode 100644 src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H create mode 100644 src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallFunctionsI.H create mode 100644 src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallViscosityI.H create mode 100644 src/turbulenceModels/RAS/incompressible/kOmega/wallOmegaI.H diff --git a/src/turbulenceModels/RAS/incompressible/Make/files b/src/turbulenceModels/RAS/incompressible/Make/files index c7746e7fd2..6584fc6113 100644 --- a/src/turbulenceModels/RAS/incompressible/Make/files +++ b/src/turbulenceModels/RAS/incompressible/Make/files @@ -5,6 +5,7 @@ laminar/laminar.C kEpsilon/kEpsilon.C RNGkEpsilon/RNGkEpsilon.C realizableKE/realizableKE.C +kOmega/kOmega.C kOmegaSST/kOmegaSST.C SpalartAllmaras/SpalartAllmaras.C LRR/LRR.C diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C new file mode 100644 index 0000000000..889f38961c --- /dev/null +++ b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.C @@ -0,0 +1,268 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +\*---------------------------------------------------------------------------*/ + +#include "kOmega.H" +#include "addToRunTimeSelectionTable.H" +#include "wallFvPatch.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +defineTypeNameAndDebug(kOmega, 0); +addToRunTimeSelectionTable(RASModel, kOmega, dictionary); + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from components +kOmega::kOmega +( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& lamTransportModel +) +: + RASModel(typeName, U, phi, lamTransportModel), + + Cmu_ + ( + dimensioned::lookupOrAddToDict + ( + "betaStar", + coeffDict_, + 0.09 + ) + ), + beta_ + ( + dimensioned::lookupOrAddToDict + ( + "beta", + coeffDict_, + 0.072 + ) + ), + alphaK_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaK", + coeffDict_, + 0.5 + ) + ), + alphaOmega_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaOmega", + coeffDict_, + 0.5 + ) + ), + + omega0_("omega0", dimless/dimTime, SMALL), + omegaSmall_("omegaSmall", dimless/dimTime, SMALL), + + k_ + ( + IOobject + ( + "k", + runTime_.timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + + omega_ + ( + IOobject + ( + "omega", + runTime_.timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh_ + ), + + nut_(k_/(omega_ + omegaSmall_)) +{ +# include "kOmegaWallViscosityI.H" + + printCoeffs(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +tmp kOmega::R() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "R", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + ((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)), + k_.boundaryField().types() + ) + ); +} + + +tmp kOmega::devReff() const +{ + return tmp + ( + new volSymmTensorField + ( + IOobject + ( + "devRhoReff", + runTime_.timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + -nuEff()*dev(twoSymm(fvc::grad(U_))) + ) + ); +} + + +tmp kOmega::divDevReff(volVectorField& U) const +{ + return + ( + - fvm::laplacian(nuEff(), U) + - fvc::div(nuEff()*dev(fvc::grad(U)().T())) + ); +} + + +bool kOmega::read() +{ + if (RASModel::read()) + { + Cmu_.readIfPresent(coeffDict_); + beta_.readIfPresent(coeffDict_); + alphaK_.readIfPresent(coeffDict_); + alphaOmega_.readIfPresent(coeffDict_); + + return true; + } + else + { + return false; + } +} + + +void kOmega::correct() +{ + transportModel_.correct(); + + if (!turbulence_) + { + return; + } + + RASModel::correct(); + + volScalarField G = nut_*2*magSqr(symm(fvc::grad(U_))); + +# include "kOmegaWallFunctionsI.H" + + // Turbulence specific dissipation rate equation + tmp omegaEqn + ( + fvm::ddt(omega_) + + fvm::div(phi_, omega_) + - fvm::Sp(fvc::div(phi_), omega_) + - fvm::laplacian(DomegaEff(), omega_) + == + G*omega_/k_ + - fvm::Sp(beta_*omega_, omega_) + ); + + omegaEqn().relax(); + +# include "wallOmegaI.H" + + solve(omegaEqn); + bound(omega_, omega0_); + + + // Turbulent kinetic energy equation + tmp kEqn + ( + fvm::ddt(k_) + + fvm::div(phi_, k_) + - fvm::Sp(fvc::div(phi_), k_) + - fvm::laplacian(DkEff(), k_) + == + G + - fvm::Sp(Cmu_*omega_, k_) + ); + + kEqn().relax(); + solve(kEqn); + bound(k_, k0_); + + + // Re-calculate viscosity + nut_ = k_/omega_; + +# include "kOmegaWallViscosityI.H" + +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace incompressible +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H new file mode 100644 index 0000000000..a524f64e40 --- /dev/null +++ b/src/turbulenceModels/RAS/incompressible/kOmega/kOmega.H @@ -0,0 +1,205 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Class + Foam::incompressible::RASModels::kOmega + +Description + Standard high Reynolds-number k-omega turbulence model for + incompressible flows. + + References: + @verbatim + "Turbulence Modeling for CFD" + D. C. Wilcox, + DCW Industries, Inc., La Canada, + California, 1998. + @endverbatim + + The default model coefficients correspond to the following: + @verbatim + kOmegaCoeffs + { + Cmu 0.09; // Equivalent to betaStar + beta 0.072; + alphak 0.5; + alphaOmega 0.5; + } + @endverbatim + +SourceFiles + kOmega.C + +\*---------------------------------------------------------------------------*/ + +#ifndef kOmega_H +#define kOmega_H + +#include "RASModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ +namespace incompressible +{ +namespace RASModels +{ + +/*---------------------------------------------------------------------------*\ + Class kOmega Declaration +\*---------------------------------------------------------------------------*/ + +class kOmega +: + public RASModel +{ + // Private data + + // Model coefficients + + dimensionedScalar Cmu_; + dimensionedScalar beta_; + dimensionedScalar alphaK_; + dimensionedScalar alphaOmega_; + + + dimensionedScalar omega0_; + dimensionedScalar omegaSmall_; + + + // Fields + + volScalarField k_; + volScalarField omega_; + volScalarField nut_; + + +public: + + //- Runtime type information + TypeName("kOmega"); + + // Constructors + + //- Construct from components + kOmega + ( + const volVectorField& U, + const surfaceScalarField& phi, + transportModel& transport + ); + + + // Destructor + + ~kOmega() + {} + + + // Member Functions + + //- Return the turbulence viscosity + tmp nut() const + { + return nut_; + } + + //- Return the effective diffusivity for k + tmp DkEff() const + { + return tmp + ( + new volScalarField("DkEff", alphaK_*nut_ + nu()) + ); + } + + //- Return the effective diffusivity for omega + tmp DomegaEff() const + { + return tmp + ( + new volScalarField("DomegaEff", alphaOmega_*nut_ + nu()) + ); + } + + //- Return the turbulence kinetic energy + tmp k() const + { + return k_; + } + + //- Return the turbulence specific dissipation rate + tmp omega() const + { + return omega_; + } + + //- Return the turbulence kinetic energy dissipation rate + tmp epsilon() const + { + return tmp + ( + new volScalarField + ( + IOobject + ( + "epsilon", + mesh_.time().timeName(), + mesh_ + ), + Cmu_*k_*omega_, + omega_.boundaryField().types() + ) + ); + } + + //- Return the Reynolds stress tensor + tmp R() const; + + //- Return the effective stress tensor including the laminar stress + tmp devReff() const; + + //- Return the source term for the momentum equation + tmp divDevReff(volVectorField& U) const; + + //- Solve the turbulence equations and correct the turbulence viscosity + void correct(); + + //- Read turbulenceProperties dictionary + bool read(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace RASModels +} // End namespace incompressible +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallFunctionsI.H b/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallFunctionsI.H new file mode 100644 index 0000000000..aca757443e --- /dev/null +++ b/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallFunctionsI.H @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Global + kOmegaWallFunctions + +Description + Calculate wall generation and frequency omega from wall-functions. + +\*---------------------------------------------------------------------------*/ + +{ + labelList cellBoundaryFaceCount(omega_.size(), 0); + + scalar Cmu25 = pow(Cmu_.value(), 0.25); + + const fvPatchList& patches = mesh_.boundary(); + + //- Initialise the near-wall omega and G fields to zero + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isType(curPatch)) + { + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + omega_[faceCelli] = 0.0; + G[faceCelli] = 0.0; + } + } + } + + //- Accumulate the wall face contributions to omega and G + // Increment cellBoundaryFaceCount for each face for averaging + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isType(curPatch)) + { +# include "checkkOmegaPatchFieldTypes.H" + + const scalarField& nuw = nu().boundaryField()[patchi]; + const scalarField& nutw = nut_.boundaryField()[patchi]; + + scalarField magFaceGradU = + mag(U_.boundaryField()[patchi].snGrad()); + + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + scalar yPlus = + Cmu25*y_[patchi][facei] + *sqrt(k_[faceCelli]) + /nuw[facei]; + + // For corner cells (with two boundary or more faces), + // omega and G in the near-wall cell are calculated + // as an average + + cellBoundaryFaceCount[faceCelli]++; + + omega_[faceCelli] += + sqrt(k_[faceCelli]) + /(Cmu25*kappa_.value()*y_[patchi][facei]); + + if (yPlus > yPlusLam_) + { + G[faceCelli] += + (nutw[facei] + nuw[facei]) + *magFaceGradU[facei] + *Cmu25*sqrt(k_[faceCelli]) + /(kappa_.value()*y_[patchi][facei]); + } + } + } + } + + + // Perform the averaging + + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isType(curPatch)) + { + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + omega_[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + G[faceCelli] /= cellBoundaryFaceCount[faceCelli]; + } + } + } +} + + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallViscosityI.H b/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallViscosityI.H new file mode 100644 index 0000000000..0f0571c93c --- /dev/null +++ b/src/turbulenceModels/RAS/incompressible/kOmega/kOmegaWallViscosityI.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Global + kOmegaWallViscosity + +Description + Calculate wall viscosity from wall-functions. + +\*---------------------------------------------------------------------------*/ + +{ + scalar Cmu25 = pow(Cmu_.value(), 0.25); + + const fvPatchList& patches = mesh_.boundary(); + + forAll(patches, patchi) + { + const fvPatch& curPatch = patches[patchi]; + + if (isType(curPatch)) + { + const scalarField& nuw = nu().boundaryField()[patchi]; + scalarField& nutw = nut_.boundaryField()[patchi]; + + forAll(curPatch, facei) + { + label faceCelli = curPatch.faceCells()[facei]; + + scalar yPlus = + Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei]; + + if (yPlus > yPlusLam_) + { + nutw[facei] = + nuw[facei] + *(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1); + } + else + { + nutw[facei] = 0.0; + } + } + } + } +} + + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/incompressible/kOmega/wallOmegaI.H b/src/turbulenceModels/RAS/incompressible/kOmega/wallOmegaI.H new file mode 100644 index 0000000000..eb1b3c5919 --- /dev/null +++ b/src/turbulenceModels/RAS/incompressible/kOmega/wallOmegaI.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 1991-2008 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 2 of the License, or (at your + option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM; if not, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Global + wallOmega + +Description + Set wall dissipation in the omega matrix + +\*---------------------------------------------------------------------------*/ + +{ + const fvPatchList& patches = mesh_.boundary(); + + forAll(patches, patchi) + { + const fvPatch& p = patches[patchi]; + + if (isType(p)) + { + omegaEqn().setValues + ( + p.faceCells(), + omega_.boundaryField()[patchi].patchInternalField() + ); + } + } +} + +// ************************************************************************* // diff --git a/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.H b/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.H index ef44b6a9ac..5564d84555 100644 --- a/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.H +++ b/src/turbulenceModels/RAS/incompressible/kOmegaSST/kOmegaSST.H @@ -238,6 +238,7 @@ public: return k_; } + //- Return the turbulence specific dissipation rate tmp omega() const { return omega_; From 44a19bc9034588e39ac5b585f318710bb2617bd8 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:34:11 +0100 Subject: [PATCH 14/39] built into user area --- applications/solvers/Lagrangian/kinematicParcelFoam/Make/files | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files b/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files index f98171fc93..a63fc64fa7 100644 --- a/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files +++ b/applications/solvers/Lagrangian/kinematicParcelFoam/Make/files @@ -1,3 +1,3 @@ kinematicParcelFoam.C -EXE = $(FOAM_USER_APPBIN)/kinematicParcelFoam +EXE = $(FOAM_APPBIN)/kinematicParcelFoam From 5b80d5d23e5b01b81405c57e79e17b445e69b492 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:34:24 +0100 Subject: [PATCH 15/39] added comment --- .../Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C b/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C index 290480583c..57b3dc5d13 100644 --- a/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C +++ b/applications/solvers/Lagrangian/kinematicParcelFoam/kinematicParcelFoam.C @@ -26,7 +26,8 @@ Application kinematicParcelFoam Description - Transient solver a single kinematicCloud. + Transient solver for a single kinematicCloud. Uses precalculated velocity + field to evolve a cloud. \*---------------------------------------------------------------------------*/ From 23d99c934dfb7218b036a93d58d3b99f4e201348 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:34:49 +0100 Subject: [PATCH 16/39] allow purgeWrite with variable time --- src/OpenFOAM/db/Time/TimeIO.C | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/OpenFOAM/db/Time/TimeIO.C b/src/OpenFOAM/db/Time/TimeIO.C index c74b674389..4299d05b5e 100644 --- a/src/OpenFOAM/db/Time/TimeIO.C +++ b/src/OpenFOAM/db/Time/TimeIO.C @@ -86,13 +86,6 @@ void Foam::Time::readDict() purgeWrite_ = 0; } - - if (writeControl_ != wcTimeStep && purgeWrite_ > 0) - { - FatalIOErrorIn("Time::readDict()", controlDict_) - << "writeControl must be set to timeStep for purgeWrite " - << exit(FatalIOError); - } } if (controlDict_.found("timeFormat")) From 00378a320280a587938a2b68205a3897ccd48edb Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:35:56 +0100 Subject: [PATCH 17/39] coupled edges did not get combined upon removal of faces on other procs --- .../polyTopoChange/removeFaces.C | 20 ++++++++----------- 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C index 28fe2a21ee..fab2b74609 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeFaces.C @@ -807,7 +807,7 @@ void Foam::removeFaces::setRefinement // Edges to remove labelHashSet edgesToRemove(faceLabels.size()); - // Per face the region it is. -1 for removed faces, -2 for regions + // Per face the region it is in. -1 for removed faces, -2 for regions // consisting of single face only. labelList faceRegion(mesh_.nFaces(), -1); @@ -1258,10 +1258,15 @@ void Foam::removeFaces::setRefinement // are only used by 2 unremoved edges. { // Usage of points by non-removed edges. - labelList nEdgesPerPoint(mesh_.nPoints(), labelMax); + labelList nEdgesPerPoint(mesh_.nPoints()); const labelListList& pointEdges = mesh_.pointEdges(); + forAll(pointEdges, pointI) + { + nEdgesPerPoint[pointI] = pointEdges[pointI].size(); + } + forAllConstIter(labelHashSet, edgesToRemove, iter) { // Edge will get removed. @@ -1269,16 +1274,7 @@ void Foam::removeFaces::setRefinement forAll(e, i) { - label pointI = e[i]; - - if (nEdgesPerPoint[pointI] == labelMax) - { - nEdgesPerPoint[pointI] = pointEdges[pointI].size()-1; - } - else - { - nEdgesPerPoint[pointI]--; - } + nEdgesPerPoint[e[i]]--; } } From 6093824d6ceeec1c0392eb9d9eeee14a46564be6 Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:36:17 +0100 Subject: [PATCH 18/39] better error message --- src/lagrangian/basic/Cloud/CloudIO.C | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lagrangian/basic/Cloud/CloudIO.C b/src/lagrangian/basic/Cloud/CloudIO.C index ea37ea27a9..0401955801 100644 --- a/src/lagrangian/basic/Cloud/CloudIO.C +++ b/src/lagrangian/basic/Cloud/CloudIO.C @@ -118,7 +118,8 @@ void Foam::Cloud::checkFieldIOobject "void Cloud::checkFieldIOobject" "(Cloud, IOField)" ) << "Size of " << data.name() - << " field does not match the number of particles" + << " field " << data.size() + << " does not match the number of particles " << c.size() << abort(FatalError); } } From ec1b7f7022c89e97a2d53a306950e2e89562f9cb Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:36:47 +0100 Subject: [PATCH 19/39] switch to indexedOctree --- src/meshTools/meshSearch/meshSearch.C | 391 +++++++++++++++++++------- src/meshTools/meshSearch/meshSearch.H | 120 +++++--- 2 files changed, 376 insertions(+), 135 deletions(-) diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C index 58c1d55e40..d71178eb86 100644 --- a/src/meshTools/meshSearch/meshSearch.C +++ b/src/meshTools/meshSearch/meshSearch.C @@ -22,14 +22,11 @@ License along with OpenFOAM; if not, write to the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA -Description - \*---------------------------------------------------------------------------*/ -#include "polyMesh.H" #include "meshSearch.H" -#include "triPointRef.H" -#include "octree.H" +#include "polyMesh.H" +#include "indexedOctree.H" #include "pointIndexHit.H" #include "DynamicList.H" #include "demandDrivenData.H" @@ -43,14 +40,75 @@ Foam::scalar Foam::meshSearch::tol_ = 1E-3; // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +bool Foam::meshSearch::findNearer +( + const point& sample, + const pointField& points, + label& nearestI, + scalar& nearestDistSqr +) +{ + bool nearer = false; + + forAll(points, pointI) + { + scalar distSqr = magSqr(points[pointI] - sample); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + nearestI = pointI; + nearer = true; + } + } + + return nearer; +} + + +bool Foam::meshSearch::findNearer +( + const point& sample, + const pointField& points, + const labelList& indices, + label& nearestI, + scalar& nearestDistSqr +) +{ + bool nearer = false; + + forAll(indices, i) + { + label pointI = indices[i]; + + scalar distSqr = magSqr(points[pointI] - sample); + + if (distSqr < nearestDistSqr) + { + nearestDistSqr = distSqr; + nearestI = i; + nearer = true; + } + } + + return nearer; +} + + // tree based searching Foam::label Foam::meshSearch::findNearestCellTree(const point& location) const { - treeBoundBox tightest(treeBoundBox::greatBox); + const indexedOctree& tree = cellCentreTree(); - scalar tightestDist = GREAT; + scalar span = mag(tree.bb().max() - tree.bb().min()); - return cellCentreTree().findNearest(location, tightest, tightestDist); + pointIndexHit info = tree.findNearest(location, Foam::sqr(span)); + + if (!info.hit()) + { + info = tree.findNearest(location, Foam::sqr(GREAT)); + } + return info.index(); } @@ -60,18 +118,15 @@ Foam::label Foam::meshSearch::findNearestCellLinear(const point& location) const const vectorField& centres = mesh_.cellCentres(); label nearestCelli = 0; - scalar minProximity = magSqr(centres[0] - location); + scalar minProximity = magSqr(centres[nearestCelli] - location); - forAll(centres, celli) - { - scalar proximity = magSqr(centres[celli] - location); - - if (proximity < minProximity) - { - nearestCelli = celli; - minProximity = proximity; - } - } + findNearer + ( + location, + centres, + nearestCelli, + minProximity + ); return nearestCelli; } @@ -105,30 +160,137 @@ Foam::label Foam::meshSearch::findNearestCellWalk do { - closer = false; - - // set the current list of neighbouring cells - const labelList& neighbours = cc[curCell]; - - forAll (neighbours, nI) - { - scalar curDistSqr = magSqr(centres[neighbours[nI]] - location); - - // search through all the neighbours. - // If the cell is closer, reset current cell and distance - if (curDistSqr < distanceSqr) - { - distanceSqr = curDistSqr; - curCell = neighbours[nI]; - closer = true; // a closer neighbour has been found - } - } + closer = findNearer + ( + location, + centres, + cc[curCell], + curCell, + distanceSqr + ); } while (closer); return curCell; } +// tree based searching +Foam::label Foam::meshSearch::findNearestFaceTree(const point& location) const +{ + // Search nearest cell centre. + const indexedOctree& tree = cellCentreTree(); + + scalar span = mag(tree.bb().max() - tree.bb().min()); + + // Search with decent span + pointIndexHit info = tree.findNearest(location, Foam::sqr(span)); + + if (!info.hit()) + { + // Search with desparate span + info = tree.findNearest(location, Foam::sqr(GREAT)); + } + + + // Now check any of the faces of the nearest cell + const vectorField& centres = mesh_.faceCentres(); + const cell& ownFaces = mesh_.cells()[info.index()]; + label nearestFaceI = ownFaces[0]; + scalar minProximity = magSqr(centres[nearestFaceI] - location); + + findNearer + ( + location, + centres, + ownFaces, + nearestFaceI, + minProximity + ); + + return nearestFaceI; +} + + +// linear searching +Foam::label Foam::meshSearch::findNearestFaceLinear(const point& location) const +{ + const vectorField& centres = mesh_.faceCentres(); + + label nearestFaceI = 0; + scalar minProximity = magSqr(centres[nearestFaceI] - location); + + findNearer + ( + location, + centres, + nearestFaceI, + minProximity + ); + + return nearestFaceI; +} + + +// walking from seed +Foam::label Foam::meshSearch::findNearestFaceWalk +( + const point& location, + const label seedFaceI +) const +{ + if (seedFaceI < 0) + { + FatalErrorIn + ( + "meshSearch::findNearestFaceWalk(const point&, const label)" + ) << "illegal seedFace:" << seedFaceI << exit(FatalError); + } + + const vectorField& centres = mesh_.faceCentres(); + + + // Walk in direction of face that decreases distance + + label curFace = seedFaceI; + scalar distanceSqr = magSqr(centres[curFace] - location); + + while (true) + { + label betterFace = curFace; + + findNearer + ( + location, + centres, + mesh_.cells()[mesh_.faceOwner()[curFace]], + betterFace, + distanceSqr + ); + + if (mesh_.isInternalFace(curFace)) + { + findNearer + ( + location, + centres, + mesh_.cells()[mesh_.faceNeighbour()[curFace]], + betterFace, + distanceSqr + ); + } + + if (betterFace == curFace) + { + break; + } + + curFace = betterFace; + } + + return curFace; +} + + Foam::label Foam::meshSearch::findCellLinear(const point& location) const { bool cellFound = false; @@ -180,12 +342,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk const face& f = mesh_.faces()[curFaceI]; - scalar minDist = - f.nearestPoint - ( - location, - mesh_.points() - ).distance(); + scalar minDist = f.nearestPoint + ( + location, + mesh_.points() + ).distance(); bool closer; @@ -202,8 +363,7 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk forAll (myEdges, myEdgeI) { - const labelList& neighbours = - mesh_.edgeFaces()[myEdges[myEdgeI]]; + const labelList& neighbours = mesh_.edgeFaces()[myEdges[myEdgeI]]; // Check any face which uses edge, is boundary face and // is not curFaceI itself. @@ -220,12 +380,11 @@ Foam::label Foam::meshSearch::findNearestBoundaryFaceWalk { const face& f = mesh_.faces()[faceI]; - pointHit curHit = - f.nearestPoint - ( - location, - mesh_.points() - ); + pointHit curHit = f.nearestPoint + ( + location, + mesh_.points() + ); // If the face is closer, reset current face and distance if (curHit.distance() < minDist) @@ -283,9 +442,11 @@ Foam::meshSearch::~meshSearch() clearOut(); } + // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::octree& Foam::meshSearch::boundaryTree() const +const Foam::indexedOctree& Foam::meshSearch::boundaryTree() + const { if (!boundaryTreePtr_) { @@ -294,17 +455,26 @@ const Foam::octree& Foam::meshSearch::boundaryTree() const // // all boundary faces (not just walls) - octreeDataFace shapes(mesh_); + labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces()); + forAll(bndFaces, i) + { + bndFaces[i] = mesh_.nInternalFaces() + i; + } treeBoundBox overallBb(mesh_.points()); - boundaryTreePtr_ = new octree + boundaryTreePtr_ = new indexedOctree ( - overallBb, // overall search domain - shapes, // all information needed to do checks on cells - 1, // min levels - 20.0, // maximum ratio of cubes v.s. cells - 10.0 + treeDataFace // all information needed to search faces + ( + false, // do not cache bb + mesh_, + bndFaces // boundary faces only + ), + overallBb, // overall search domain + 8, // maxLevel + 10, // leafsize + 3.0 // duplicity ); } @@ -312,7 +482,8 @@ const Foam::octree& Foam::meshSearch::boundaryTree() const } -const Foam::octree& Foam::meshSearch::cellTree() const +const Foam::indexedOctree& Foam::meshSearch::cellTree() + const { if (!cellTreePtr_) { @@ -320,17 +491,19 @@ const Foam::octree& Foam::meshSearch::cellTree() const // Construct tree // - octreeDataCell shapes(mesh_); - treeBoundBox overallBb(mesh_.points()); - cellTreePtr_ = new octree + cellTreePtr_ = new indexedOctree ( + treeDataCell + ( + false, // not cache bb + mesh_ + ), overallBb, // overall search domain - shapes, // all information needed to do checks on cells - 1, // min levels - 20.0, // maximum ratio of cubes v.s. cells - 2.0 + 8, // maxLevel + 10, // leafsize + 3.0 // duplicity ); } @@ -339,8 +512,8 @@ const Foam::octree& Foam::meshSearch::cellTree() const } -const Foam::octree& Foam::meshSearch::cellCentreTree() - const +const Foam::indexedOctree& + Foam::meshSearch::cellCentreTree() const { if (!cellCentreTreePtr_) { @@ -348,17 +521,14 @@ const Foam::octree& Foam::meshSearch::cellCentreTree() // Construct tree // - // shapes holds reference to cellCentres! - octreeDataPoint shapes(mesh_.cellCentres()); - treeBoundBox overallBb(mesh_.cellCentres()); - cellCentreTreePtr_ = new octree + cellCentreTreePtr_ = new indexedOctree ( + treeDataPoint(mesh_.cellCentres()), overallBb, // overall search domain - shapes, // all information needed to do checks on cells - 1, // min levels - 20.0, // maximum ratio of cubes v.s. cells + 10, // max levels + 10.0, // maximum ratio of cubes v.s. cells 100.0 // max. duplicity; n/a since no bounding boxes. ); } @@ -396,15 +566,14 @@ bool Foam::meshSearch::pointInCell(const point& p, label cellI) const forAll(f, fp) { - pointHit inter = - f.ray - ( - ctr, - dir, - mesh_.points(), - intersection::HALF_RAY, - intersection::VECTOR - ); + pointHit inter = f.ray + ( + ctr, + dir, + mesh_.points(), + intersection::HALF_RAY, + intersection::VECTOR + ); if (inter.hit()) { @@ -484,6 +653,31 @@ Foam::label Foam::meshSearch::findNearestCell } +Foam::label Foam::meshSearch::findNearestFace +( + const point& location, + const label seedFaceI, + const bool useTreeSearch +) const +{ + if (seedFaceI == -1) + { + if (useTreeSearch) + { + return findNearestFaceTree(location); + } + else + { + return findNearestFaceLinear(location); + } + } + else + { + return findNearestFaceWalk(location, seedFaceI); + } +} + + Foam::label Foam::meshSearch::findCell ( const point& location, @@ -607,19 +801,26 @@ Foam::label Foam::meshSearch::findNearestBoundaryFace { if (useTreeSearch) { - treeBoundBox tightest(treeBoundBox::greatBox); + const indexedOctree& tree = boundaryTree(); - scalar tightestDist = GREAT; + scalar span = mag(tree.bb().max() - tree.bb().min()); - label index = - boundaryTree().findNearest + pointIndexHit info = boundaryTree().findNearest + ( + location, + Foam::sqr(span) + ); + + if (!info.hit()) + { + info = boundaryTree().findNearest ( location, - tightest, - tightestDist + Foam::sqr(GREAT) ); + } - return boundaryTree().shapes().meshFaces()[index]; + return tree.shapes().faceLabels()[info.index()]; } else { @@ -670,7 +871,7 @@ Foam::pointIndexHit Foam::meshSearch::intersection if (curHit.hit()) { // Change index into octreeData into face label - curHit.setIndex(boundaryTree().shapes().meshFaces()[curHit.index()]); + curHit.setIndex(boundaryTree().shapes().faceLabels()[curHit.index()]); } return curHit; } @@ -733,7 +934,9 @@ Foam::List Foam::meshSearch::intersections bool Foam::meshSearch::isInside(const point& p) const { - return boundaryTree().getSampleType(p) == octree::INSIDE; + return + boundaryTree().getVolumeType(p) + == indexedOctree::INSIDE; } diff --git a/src/meshTools/meshSearch/meshSearch.H b/src/meshTools/meshSearch/meshSearch.H index c50c0b8bb0..38ca70fcd9 100644 --- a/src/meshTools/meshSearch/meshSearch.H +++ b/src/meshTools/meshSearch/meshSearch.H @@ -26,7 +26,8 @@ Class Foam::meshSearch Description - Various searches on polyMesh; uses (demand driven) octree to search. + Various (local, not parallel) searches on polyMesh; + uses (demand driven) octree to search. SourceFiles meshSearch.C @@ -36,11 +37,10 @@ SourceFiles #ifndef meshSearch_H #define meshSearch_H -#include "octreeDataCell.H" -#include "octreeDataFace.H" -#include "octreeDataPoint.H" +#include "treeDataCell.H" +#include "treeDataFace.H" +#include "treeDataPoint.H" #include "pointIndexHit.H" -#include "className.H" #include "Cloud.H" #include "passiveParticle.H" @@ -71,46 +71,77 @@ class meshSearch //- demand driven octrees - mutable octree* boundaryTreePtr_; - mutable octree* cellTreePtr_; - mutable octree* cellCentreTreePtr_; + mutable indexedOctree* boundaryTreePtr_; + mutable indexedOctree* cellTreePtr_; + mutable indexedOctree* cellCentreTreePtr_; // Private Member Functions - //- nearest cell centre using octree - label findNearestCellTree(const point& location) const; - - //- nearest cell centre going through all cells - label findNearestCellLinear(const point& location) const; - - //- walk from seed. Does not 'go around' boundary, just returns - // last cell before boundary. - label findNearestCellWalk + //- Updates nearestI, nearestDistSqr from any closer ones. + static bool findNearer ( - const point& location, - const label seedCellI - ) const; + const point& sample, + const pointField& points, + label& nearestI, + scalar& nearestDistSqr + ); - //- cell containing location. Linear search. - label findCellLinear(const point& location) const; - - //- walk from seed to find nearest boundary face. Gets stuck in - // local minimum. - label findNearestBoundaryFaceWalk + //- Updates nearestI, nearestDistSqr from any selected closer ones. + static bool findNearer ( - const point& location, - const label seedFaceI - ) const; + const point& sample, + const pointField& points, + const labelList& indices, + label& nearestI, + scalar& nearestDistSqr + ); - //- Calculate offset vector in direction dir with as length a fraction - // of the cell size (of the cell straddling boundary face) - vector offset - ( - const point& bPoint, - const label bFaceI, - const vector& dir - ) const; + + // Cells + + //- nearest cell centre using octree + label findNearestCellTree(const point&) const; + + //- nearest cell centre going through all cells + label findNearestCellLinear(const point&) const; + + //- walk from seed. Does not 'go around' boundary, just returns + // last cell before boundary. + label findNearestCellWalk(const point&, const label) const; + + //- cell containing location. Linear search. + label findCellLinear(const point&) const; + + + // Cells + + label findNearestFaceTree(const point&) const; + + label findNearestFaceLinear(const point&) const; + + label findNearestFaceWalk(const point&, const label) const; + + + + // Boundary faces + + //- walk from seed to find nearest boundary face. Gets stuck in + // local minimum. + label findNearestBoundaryFaceWalk + ( + const point& location, + const label seedFaceI + ) const; + + //- Calculate offset vector in direction dir with as length a + // fraction of the cell size (of the cell straddling boundary face) + vector offset + ( + const point& bPoint, + const label bFaceI, + const vector& dir + ) const; //- Disallow default bitwise copy construct @@ -154,13 +185,13 @@ public: //- Get (demand driven) reference to octree holding all // boundary faces - const octree& boundaryTree() const; + const indexedOctree& boundaryTree() const; //- Get (demand driven) reference to octree holding all cells - const octree& cellTree() const; + const indexedOctree& cellTree() const; //- Get (demand driven) reference to octree holding all cell centres - const octree& cellCentreTree() const; + const indexedOctree& cellCentreTree() const; // Queries @@ -181,6 +212,13 @@ public: const bool useTreeSearch = true ) const; + label findNearestFace + ( + const point& location, + const label seedFaceI, + const bool useTreeSearch + ) const; + //- Find cell containing (using pointInCell) location. // If seed provided walks and falls back to linear/tree search. // (so handles holes correctly)s @@ -206,7 +244,7 @@ public: //- Find first intersection of boundary in segment [pStart, pEnd] // (so inclusive of endpoints). Always octree for now pointIndexHit intersection(const point& pStart, const point& pEnd) - const; + const; //- Find all intersections of boundary within segment pStart .. pEnd // Always octree for now From a7fb0a36b3d319b088053a827f0c12f0d736642f Mon Sep 17 00:00:00 2001 From: mattijs Date: Tue, 9 Sep 2008 12:47:31 +0100 Subject: [PATCH 20/39] checking for hanging points; better debugging --- .../polyTopoChange/polyTopoChange/hexRef8.C | 354 ++++++++++++++---- 1 file changed, 284 insertions(+), 70 deletions(-) diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C index 1b3bf6e19b..0c567ee312 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8.C @@ -62,6 +62,7 @@ namespace Foam }; } + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // void Foam::hexRef8::reorder @@ -80,7 +81,7 @@ void Foam::hexRef8::reorder if (newI >= len) { - FatalErrorIn("hexRef8::reorder") << abort(FatalError); + FatalErrorIn("hexRef8::reorder(..)") << abort(FatalError); } if (newI >= 0) @@ -557,22 +558,11 @@ Foam::label Foam::hexRef8::getAnchorCell } // Problem. - - // Pick up points of the cell - const labelList cPoints(cellPoints(cellI)); - - Perr<< "cell:" << cellI << " points:" << endl; - forAll(cPoints, i) - { - label pointI = cPoints[i]; - - Perr<< " " << pointI << " coord:" << mesh_.points()[pointI] - << nl; - } + dumpCell(cellI); Perr<< "cell:" << cellI << " anchorPoints:" << cellAnchorPoints[cellI] << endl; - FatalErrorIn("hexRef8::getAnchorCell") + FatalErrorIn("hexRef8::getAnchorCell(..)") << "Could not find point " << pointI << " in the anchorPoints for cell " << cellI << endl << "Does your original mesh obey the 2:1 constraint and" @@ -690,9 +680,50 @@ Foam::label Foam::hexRef8::countAnchors } +void Foam::hexRef8::dumpCell(const label cellI) const +{ + OFstream str(mesh_.time().path()/"cell_" + Foam::name(cellI) + ".obj"); + Pout<< "hexRef8 : Dumping cell as obj to " << str.name() << endl; + + const cell& cFaces = mesh_.cells()[cellI]; + + Map