diff --git a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H index 5e919efcc4..5cd98cd9d6 100644 --- a/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H +++ b/applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H @@ -124,11 +124,6 @@ public: // Member Functions - tmp mut() const - { - return mut_; - } - //- Return the effective diffusivity for k tmp DkEff() const { @@ -147,41 +142,44 @@ public: ); } - //- Return the effective turbulent thermal diffusivity - tmp alphaEff() const + //- Return the turbulence viscosity + virtual tmp mut() const { - return tmp - ( - new volScalarField("alphaEff", alphat_ + alpha()) - ); + return mut_; + } + + //- Return the turbulence thermal diffusivity + virtual tmp alphat() const + { + return alphat_; } //- Return the turbulence kinetic energy - tmp k() const + virtual tmp k() const { return k_; } //- Return the turbulence kinetic energy dissipation rate - tmp epsilon() const + virtual tmp epsilon() const { return epsilon_; } //- Return the Reynolds stress tensor - tmp R() const; + virtual tmp R() const; //- Return the effective stress tensor including the laminar stress - tmp devRhoReff() const; + virtual tmp devRhoReff() const; //- Return the source term for the momentum equation - tmp divDevRhoReff(volVectorField& U) const; + virtual tmp divDevRhoReff(volVectorField& U) const; //- Solve the turbulence equations and correct the turbulence viscosity - void correct(); + virtual void correct(); //- Read turbulenceProperties dictionary - bool read(); + virtual bool read(); }; diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index c7f2d5b6f3..783469ee7d 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -40,6 +40,10 @@ geometry { type triSurfaceMesh; + //tolerance 1E-6; // optional:non-default tolerance on intersections + //maxTreeDepth 10; // optional:depth of octree. Decrease only in case + // of memory limitations. + // Per region the patchname. If not provided will be _. regions { diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C index 7833b4ec24..cf24a1c447 100644 --- a/applications/utilities/mesh/manipulation/setSet/setSet.C +++ b/applications/utilities/mesh/manipulation/setSet/setSet.C @@ -94,7 +94,7 @@ void backup { if (fromSet.size()) { - Pout<< " Backing up " << fromName << " into " << toName << endl; + Info<< " Backing up " << fromName << " into " << toName << endl; topoSet::New(setType, mesh, toName, fromSet)().write(); } @@ -525,7 +525,7 @@ bool doCommand { topoSet& currentSet = currentSetPtr(); - Pout<< " Set:" << currentSet.name() + Info<< " Set:" << currentSet.name() << " Size:" << currentSet.size() << " Action:" << actionName << endl; @@ -622,7 +622,7 @@ bool doCommand + ".vtk" ); - Pout<< " Writing " << currentSet.name() + Info<< " Writing " << currentSet.name() << " (size " << currentSet.size() << ") to " << currentSet.instance()/currentSet.local() /currentSet.name() @@ -634,7 +634,7 @@ bool doCommand } else { - Pout<< " Writing " << currentSet.name() + Info<< " Writing " << currentSet.name() << " (size " << currentSet.size() << ") to " << currentSet.instance()/currentSet.local() /currentSet.name() << endl << endl; @@ -683,7 +683,7 @@ enum commandStatus void printMesh(const Time& runTime, const polyMesh& mesh) { - Pout<< "Time:" << runTime.timeName() + Info<< "Time:" << runTime.timeName() << " cells:" << mesh.nCells() << " faces:" << mesh.nFaces() << " points:" << mesh.nPoints() @@ -703,19 +703,19 @@ commandStatus parseType { if (setType.empty()) { - Pout<< "Type 'help' for usage information" << endl; + Info<< "Type 'help' for usage information" << endl; return INVALID; } else if (setType == "help") { - printHelp(Pout); + printHelp(Info); return INVALID; } else if (setType == "list") { - printAllSets(mesh, Pout); + printAllSets(mesh, Info); return INVALID; } @@ -726,7 +726,7 @@ commandStatus parseType label nearestIndex = Time::findClosestTimeIndex(Times, requestedTime); - Pout<< "Changing time from " << runTime.timeName() + Info<< "Changing time from " << runTime.timeName() << " to " << Times[nearestIndex].name() << endl; @@ -737,24 +737,24 @@ commandStatus parseType { case polyMesh::UNCHANGED: { - Pout<< " mesh not changed." << endl; + Info<< " mesh not changed." << endl; break; } case polyMesh::POINTS_MOVED: { - Pout<< " points moved; topology unchanged." << endl; + Info<< " points moved; topology unchanged." << endl; break; } case polyMesh::TOPO_CHANGE: { - Pout<< " topology changed; patches unchanged." << nl + Info<< " topology changed; patches unchanged." << nl << " "; printMesh(runTime, mesh); break; } case polyMesh::TOPO_PATCH_CHANGE: { - Pout<< " topology changed and patches changed." << nl + Info<< " topology changed and patches changed." << nl << " "; printMesh(runTime, mesh); @@ -773,7 +773,7 @@ commandStatus parseType } else if (setType == "quit") { - Pout<< "Quitting ..." << endl; + Info<< "Quitting ..." << endl; return QUIT; } @@ -864,7 +864,7 @@ int main(int argc, char *argv[]) printMesh(runTime, mesh); // Print current sets - printAllSets(mesh, Pout); + printAllSets(mesh, Info); @@ -874,7 +874,7 @@ int main(int argc, char *argv[]) { fileName batchFile(args.option("batch")); - Pout<< "Reading commands from file " << batchFile << endl; + Info<< "Reading commands from file " << batchFile << endl; // we cannot handle .gz files if (!isFile(batchFile, false)) @@ -888,11 +888,11 @@ int main(int argc, char *argv[]) #if READLINE != 0 else if (!read_history(historyFile)) { - Pout<< "Successfully read history from " << historyFile << endl; + Info<< "Successfully read history from " << historyFile << endl; } #endif - Pout<< "Please type 'help', 'quit' or a set command after prompt." << endl; + Info<< "Please type 'help', 'quit' or a set command after prompt." << endl; bool ok = true; @@ -916,7 +916,7 @@ int main(int argc, char *argv[]) { if (!fileStreamPtr->good()) { - Pout<< "End of batch file" << endl; + Info<< "End of batch file" << endl; break; } @@ -924,7 +924,7 @@ int main(int argc, char *argv[]) if (rawLine.size()) { - Pout<< "Doing:" << rawLine << endl; + Info<< "Doing:" << rawLine << endl; } } else @@ -945,7 +945,7 @@ int main(int argc, char *argv[]) } # else { - Pout<< "Command>" << flush; + Info<< "Command>" << flush; std::getline(std::cin, rawLine); } # endif @@ -992,7 +992,7 @@ int main(int argc, char *argv[]) delete fileStreamPtr; } - Pout<< "\nEnd\n" << endl; + Info<< "\nEnd\n" << endl; return 0; } diff --git a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C index 70927ca8f3..553ea07138 100644 --- a/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C +++ b/applications/utilities/mesh/manipulation/stitchMesh/stitchMesh.C @@ -96,6 +96,7 @@ void checkPatch(const polyBoundaryMesh& bMesh, const word& name) int main(int argc, char *argv[]) { argList::noParallel(); +# include "addRegionOption.H" argList::validArgs.append("masterPatch"); argList::validArgs.append("slavePatch"); @@ -109,7 +110,7 @@ int main(int argc, char *argv[]) # include "setRootCase.H" # include "createTime.H" runTime.functionObjects().off(); -# include "createMesh.H" +# include "createNamedMesh.H" const word oldInstance = mesh.pointsInstance(); diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C index 7e4363ffb0..cd51599ca1 100644 --- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C +++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C @@ -31,6 +31,14 @@ Description Mainly used to convert binary mesh/field files to ASCII. + Problem: any zero-size List written binary gets written as '0'. When + reading the file as a dictionary this is interpreted as a label. This + is (usually) not a problem when doing patch fields since these get the + 'uniform', 'nonuniform' prefix. However zone contents are labelLists + not labelFields and these go wrong. For now hacked a solution where + we detect the keywords in zones and redo the dictionary entries + to be labelLists. + \*---------------------------------------------------------------------------*/ #include "argList.H" @@ -56,6 +64,82 @@ namespace Foam } +// Hack to do zones which have Lists in them. See above. +bool writeZones(const word& name, Time& runTime) +{ + IOobject io + ( + name, + runTime.timeName(), + polyMesh::meshSubDir, + runTime, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ); + + bool writeOk = false; + + if (io.headerOk()) + { + Info<< " Reading " << io.headerClassName() + << " : " << name << endl; + + // Switch off type checking (for reading e.g. faceZones as + // generic list of dictionaries). + const word oldTypeName = IOPtrList::typeName; + const_cast(IOPtrList::typeName) = word::null; + + IOPtrList meshObject(io); + + forAll(meshObject, i) + { + if (meshObject[i].isDict()) + { + dictionary& d = meshObject[i].dict(); + + if (d.found("faceLabels")) + { + d.set("faceLabels", labelList(d.lookup("faceLabels"))); + } + + if (d.found("flipMap")) + { + d.set("flipMap", boolList(d.lookup("flipMap"))); + } + + if (d.found("cellLabels")) + { + d.set("cellLabels", labelList(d.lookup("cellLabels"))); + } + + if (d.found("pointLabels")) + { + d.set("pointLabels", labelList(d.lookup("pointLabels"))); + } + } + } + + const_cast(IOPtrList::typeName) = oldTypeName; + // Fake type back to what was in field + const_cast(meshObject.type()) = io.headerClassName(); + + Info<< " Writing " << name << endl; + + // Force writing as ascii + writeOk = meshObject.regIOobject::writeObject + ( + IOstream::ASCII, + IOstream::currentVersion, + runTime.writeCompression() + ); + } + + return writeOk; +} + + + // Main program: int main(int argc, char *argv[]) @@ -76,9 +160,19 @@ int main(int argc, char *argv[]) writeMeshObject("neighbour", runTime); writeMeshObject("faces", runTime); writeMeshObject("points", runTime); - writeMeshObject >("cellZones", runTime); - writeMeshObject >("faceZones", runTime); - writeMeshObject >("pointZones", runTime); + writeMeshObject("pointProcAddressing", runTime); + writeMeshObject("faceProcAddressing", runTime); + writeMeshObject("cellProcAddressing", runTime); + writeMeshObject("boundaryProcAddressing", runTime); + + if (runTime.writeFormat() == IOstream::ASCII) + { + // Only do zones when converting from binary to ascii + // The other way gives problems since working on dictionary level. + writeZones("cellZones", runTime); + writeZones("faceZones", runTime); + writeZones("pointZones", runTime); + } // Get list of objects from the database IOobjectList objects(runTime, runTime.timeName()); diff --git a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C index 6c1065a489..44384b7efe 100644 --- a/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C +++ b/applications/utilities/postProcessing/dataConversion/foamToTecplot360/foamToTecplot360.C @@ -1077,83 +1077,92 @@ int main(int argc, char *argv[]) { const faceZone& pp = zones[zoneI]; - const indirectPrimitivePatch ipp - ( - IndirectList(mesh.faces(), pp), - mesh.points() - ); + if (pp.size() > 0) + { + const indirectPrimitivePatch ipp + ( + IndirectList(mesh.faces(), pp), + mesh.points() + ); - writer.writePolygonalZone - ( - pp.name(), - strandID++, //1+patchIDs.size()+zoneI, //strandID, - ipp, - allVarLocation - ); + writer.writePolygonalZone + ( + pp.name(), + strandID++, //1+patchIDs.size()+zoneI, //strandID, + ipp, + allVarLocation + ); - // Write coordinates - writer.writeField(ipp.localPoints().component(0)()); - writer.writeField(ipp.localPoints().component(1)()); - writer.writeField(ipp.localPoints().component(2)()); + // Write coordinates + writer.writeField(ipp.localPoints().component(0)()); + writer.writeField(ipp.localPoints().component(1)()); + writer.writeField(ipp.localPoints().component(2)()); - // Write all volfields - forAll(vsf, i) - { - writer.writeField - ( - writer.getFaceField + // Write all volfields + forAll(vsf, i) + { + writer.writeField ( - linearInterpolate(vsf[i])(), - pp - )() - ); - } - forAll(vvf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vsf[i])(), + pp + )() + ); + } + forAll(vvf, i) + { + writer.writeField ( - linearInterpolate(vvf[i])(), - pp - )() - ); - } - forAll(vSpheretf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vvf[i])(), + pp + )() + ); + } + forAll(vSpheretf, i) + { + writer.writeField ( - linearInterpolate(vSpheretf[i])(), - pp - )() - ); - } - forAll(vSymmtf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vSpheretf[i])(), + pp + )() + ); + } + forAll(vSymmtf, i) + { + writer.writeField ( - linearInterpolate(vSymmtf[i])(), - pp - )() - ); - } - forAll(vtf, i) - { - writer.writeField - ( - writer.getFaceField + writer.getFaceField + ( + linearInterpolate(vSymmtf[i])(), + pp + )() + ); + } + forAll(vtf, i) + { + writer.writeField ( - linearInterpolate(vtf[i])(), - pp - )() - ); - } + writer.getFaceField + ( + linearInterpolate(vtf[i])(), + pp + )() + ); + } - writer.writeConnectivity(ipp); + writer.writeConnectivity(ipp); + } + else + { + Info<< " Skipping zero sized faceZone " << zoneI + << "\t" << pp.name() + << nl << endl; + } } } diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C index 12785fec9b..f3bcb4599b 100644 --- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C +++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C @@ -202,6 +202,7 @@ void massPropertiesSolid J *= density; } + void massPropertiesShell ( const pointField& pts, @@ -272,7 +273,7 @@ tensor applyParallelAxisTheorem vector d = (refPt - cM); - return (J + m*((d & d)*I - d*d)); + return J + m*((d & d)*I - d*d); } @@ -389,7 +390,7 @@ int main(int argc, char *argv[]) << eVec.x() << ' ' << eVec.y() << ' ' << eVec.z() << endl; - if(calcAroundRefPt) + if (calcAroundRefPt) { Info << "Inertia tensor relative to " << refPt << " = " << applyParallelAxisTheorem(m, cM, J, refPt) diff --git a/etc/apps/paraview3/bashrc b/etc/apps/paraview3/bashrc index 7209dc0da8..bbeb6c0ebd 100644 --- a/etc/apps/paraview3/bashrc +++ b/etc/apps/paraview3/bashrc @@ -68,7 +68,7 @@ fi if [ -r $ParaView_DIR ] then export PATH=$ParaView_DIR/bin:$PATH - export PV_PLUGIN_PATH=$FOAM_LIBBIN + export PV_PLUGIN_PATH=$FOAM_LIBBIN/paraview fi unset cmake ParaView_PYTHON_DIR diff --git a/etc/apps/paraview3/cshrc b/etc/apps/paraview3/cshrc index 9160edf0f8..33c863215e 100644 --- a/etc/apps/paraview3/cshrc +++ b/etc/apps/paraview3/cshrc @@ -62,7 +62,7 @@ endif if ( -r $ParaView_INST_DIR ) then set path=($ParaView_DIR/bin $path) - setenv PV_PLUGIN_PATH $FOAM_LIBBIN + setenv PV_PLUGIN_PATH $FOAM_LIBBIN/paraview endif unset cmake paraviewMajor paraviewPython diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H index 5e28ac1e14..10aa63f2c9 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraint.H @@ -77,6 +77,9 @@ public: //- Apply and accumulate the effect of the given constraint direction inline void applyConstraint(const vector& cd); + //- Combine constraints + inline void combine(const pointConstraint&); + //- Return the accumulated constraint transformation tensor inline tensor constraintTransformation() const; }; diff --git a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H index a0af05b5af..c488bf603e 100644 --- a/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H +++ b/src/OpenFOAM/fields/pointPatchFields/pointPatchField/pointConstraint/pointConstraintI.H @@ -69,6 +69,47 @@ void Foam::pointConstraint::applyConstraint(const vector& cd) } +void Foam::pointConstraint::combine(const pointConstraint& pc) +{ + if (first() == 0) + { + operator=(pc); + } + else if (first() == 1) + { + // Save single normal + vector n = second(); + // Apply to supplied point constaint + operator=(pc); + applyConstraint(n); + } + else if (first() == 2) + { + if (pc.first() == 0) + {} + else if (pc.first() == 1) + { + applyConstraint(pc.second()); + } + else if (pc.first() == 2) + { + // Both constrained to line. Same (+-)direction? + if (mag(second() & pc.second()) <= (1.0-1e-3)) + { + // Different directions + first() = 3; + second() = vector::zero; + } + } + else + { + first() = 3; + second() = vector::zero; + } + } +} + + Foam::tensor Foam::pointConstraint::constraintTransformation() const { if (first() == 0) diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H index a77af7af30..ec5ae840e4 100644 --- a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H +++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H @@ -132,7 +132,7 @@ public: // associated with any faces virtual const labelList& loneMeshPoints() const; - //- Return point unit normals. Not impelemented. + //- Return point unit normals. Not implemented. virtual const vectorField& pointNormals() const; }; diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C index ee2d966836..8e18f68c77 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.C @@ -199,13 +199,20 @@ Foam::scalarField Foam::coupledPolyPatch::calcFaceTol const face& f = faces[faceI]; + // 1. calculate a typical size of the face. Use maximum distance + // to face centre scalar maxLenSqr = -GREAT; + // 2. as measure of truncation error when comparing two coordinates + // use SMALL * maximum component + scalar maxCmpt = -GREAT; forAll(f, fp) { - maxLenSqr = max(maxLenSqr, magSqr(points[f[fp]] - cc)); + const point& pt = points[f[fp]]; + maxLenSqr = max(maxLenSqr, magSqr(pt - cc)); + maxCmpt = max(maxCmpt, cmptMax(cmptMag(pt))); } - tols[faceI] = matchTol * Foam::sqrt(maxLenSqr); + tols[faceI] = max(SMALL*maxCmpt, matchTol*Foam::sqrt(maxLenSqr)); } return tols; } diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C index b50069d606..b77d2cefa4 100644 --- a/src/dynamicMesh/meshCut/directions/directions.C +++ b/src/dynamicMesh/meshCut/directions/directions.C @@ -446,19 +446,4 @@ Foam::directions::directions } -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // - - -// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // - - // ************************************************************************* // diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index f42f2a4ef4..f777081d03 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -824,7 +824,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness // (i.e. treat as if mirror cell on other side) vector faceNormal = faceAreas[faceI]; - faceNormal /= mag(faceNormal) + VSMALL; + faceNormal /= mag(faceNormal) + ROOTVSMALL; vector dOwn = faceCentres[faceI] - cellCentres[own[faceI]]; @@ -834,7 +834,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness scalar skewness = mag(faceCentres[faceI] - faceIntersection) - /(2*mag(dWall) + VSMALL); + /(2*mag(dWall) + ROOTVSMALL); // Check if the skewness vector is greater than the PN vector. // This does not cause trouble but is a good indication of a poor diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H index 062ffa1d86..5e179da0f9 100644 --- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H +++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclic/cyclicFvPatchField.H @@ -142,7 +142,7 @@ public: // Access - //- Retirn local reference cast into the cyclic patch + //- Return local reference cast into the cyclic patch const cyclicFvPatch& cyclicPatch() const { return cyclicPatch_; diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H index 7fa48ccecb..f8b94c8d5f 100644 --- a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H +++ b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/fullStencils/cellToCellStencil.H @@ -95,7 +95,7 @@ protected: class unionEqOp { public: - void operator()( labelList& x, const labelList& y ) const; + void operator()(labelList& x, const labelList& y) const; }; //- Collect cell neighbours of faces in global numbering diff --git a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H index f02475a316..3a8d69c2ef 100644 --- a/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H +++ b/src/lagrangian/intermediate/clouds/derived/BasicReactingMultiphaseCloud/BasicReactingMultiphaseCloud.H @@ -42,13 +42,31 @@ Description namespace Foam { - typedef ReactingMultiphaseCloud > + typedef ReactingMultiphaseCloud + < + BasicReactingMultiphaseParcel + < + constGasThermoPhysics + > + > constThermoReactingMultiphaseCloud; - typedef ReactingMultiphaseCloud > + typedef ReactingMultiphaseCloud + < + BasicReactingMultiphaseParcel + < + gasThermoPhysics + > + > thermoReactingMultiphaseCloud; - typedef ReactingMultiphaseCloud > + typedef ReactingMultiphaseCloud + < + BasicReactingMultiphaseParcel + < + icoPoly8ThermoPhysics + > + > icoPoly8ThermoReactingMultiphaseCloud; } diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C index 5c3d3d0685..c1914e5ed5 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcel.C @@ -401,7 +401,7 @@ void Foam::ReactingMultiphaseParcel::calc td.cloud().hcTrans()[cellI] += np0 *dMassGas[i] - *td.cloud().mcCarrierThermo().speciesData()[gid].H(T0); + *td.cloud().mcCarrierThermo().speciesData()[gid].Hc(); } forAll(YLiquid_, i) { @@ -410,7 +410,7 @@ void Foam::ReactingMultiphaseParcel::calc td.cloud().hcTrans()[cellI] += np0 *dMassLiquid[i] - *td.cloud().mcCarrierThermo().speciesData()[gid].H(T0); + *td.cloud().mcCarrierThermo().speciesData()[gid].Hc(); } /* // No mapping between solid components and carrier phase @@ -421,7 +421,7 @@ void Foam::ReactingMultiphaseParcel::calc td.cloud().hcTrans()[cellI] += np0 *dMassSolid[i] - *td.cloud().mcCarrierThermo().speciesData()[gid].H(T0); + *td.cloud().mcCarrierThermo().speciesData()[gid].Hc(); } */ forAll(dMassSRCarrier, i) @@ -430,7 +430,7 @@ void Foam::ReactingMultiphaseParcel::calc td.cloud().hcTrans()[cellI] += np0 *dMassSRCarrier[i] - *td.cloud().mcCarrierThermo().speciesData()[i].H(T0); + *td.cloud().mcCarrierThermo().speciesData()[i].Hc(); } // Update momentum transfer diff --git a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C index 524b4d00d5..92525ea2b1 100644 --- a/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C +++ b/src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C @@ -122,7 +122,7 @@ void Foam::ReactingParcel::correctSurfaceValues } // Far field carrier molar fractions - scalarField Xinf(Y_.size()); + scalarField Xinf(td.cloud().mcCarrierThermo().speciesData().size()); forAll(Xinf, i) { diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index cc49e48939..acd1673392 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -2490,13 +2490,18 @@ void Foam::autoLayerDriver::mergePatchFacesUndo << "---------------------------" << nl << " - which are on the same patch" << nl << " - which make an angle < " << layerParams.featureAngle() + << "- which are on the same patch" << nl + << "- which make an angle < " << layerParams.featureAngle() << " degrees" << nl << " (cos:" << minCos << ')' << nl << " - as long as the resulting face doesn't become concave" + << " (cos:" << minCos << ')' << nl + << "- as long as the resulting face doesn't become concave" << " by more than " << layerParams.concaveAngle() << " degrees" << nl << " (0=straight, 180=fully concave)" << nl + << " (0=straight, 180=fully concave)" << nl << endl; label nChanged = mergePatchFacesUndo(minCos, concaveCos, motionDict); @@ -2546,11 +2551,6 @@ void Foam::autoLayerDriver::addLayers ); - // Undistorted edge length - const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); - const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); - - // Point-wise extrusion data // ~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -2612,6 +2612,10 @@ void Foam::autoLayerDriver::addLayers // Disable extrusion on warped faces // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // Undistorted edge length + const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); + const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); + handleWarpedFaces ( pp, @@ -2651,6 +2655,10 @@ void Foam::autoLayerDriver::addLayers } + // Undistorted edge length + const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); + const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); + // Determine (wanted) point-wise layer thickness and expansion ratio scalarField thickness(pp().nPoints()); scalarField minThickness(pp().nPoints()); diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C index 909cf3e683..c4ea71f434 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C @@ -54,6 +54,7 @@ License #include "Random.H" #include "searchableSurfaces.H" #include "treeBoundBox.H" +#include "zeroGradientFvPatchFields.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -2128,7 +2129,8 @@ void Foam::meshRefinement::dumpRefinementLevel() const false ), mesh_, - dimensionedScalar("zero", dimless, 0) + dimensionedScalar("zero", dimless, 0), + zeroGradientFvPatchScalarField::typeName ); const labelList& cellLevel = meshCutter_.cellLevel(); diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C index c910b4c1b5..0f81dff2f5 100644 --- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C @@ -496,23 +496,63 @@ Foam::labelList Foam::refinementSurfaces::getClosedNamedSurfaces() const } -// Count number of triangles per surface region -Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s) -{ - const geometricSurfacePatchList& regions = s.patches(); +// // Count number of triangles per surface region +// Foam::labelList Foam::refinementSurfaces::countRegions(const triSurface& s) +// { +// const geometricSurfacePatchList& regions = s.patches(); +// +// labelList nTris(regions.size(), 0); +// +// forAll(s, triI) +// { +// nTris[s[triI].region()]++; +// } +// return nTris; +// } - labelList nTris(regions.size(), 0); - forAll(s, triI) - { - nTris[s[triI].region()]++; - } - return nTris; -} +// // Pre-calculate the refinement level for every element +// void Foam::refinementSurfaces::wantedRefinementLevel +// ( +// const shellSurfaces& shells, +// const label surfI, +// const List& info, // Indices +// const pointField& ctrs, // Representative coordinate +// labelList& minLevelField +// ) const +// { +// const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; +// +// // Get per element the region +// labelList region; +// geom.getRegion(info, region); +// +// // Initialise fields to region wise minLevel +// minLevelField.setSize(ctrs.size()); +// minLevelField = -1; +// +// forAll(minLevelField, i) +// { +// if (info[i].hit()) +// { +// minLevelField[i] = minLevel(surfI, region[i]); +// } +// } +// +// // Find out if triangle inside shell with higher level +// // What level does shell want to refine fc to? +// labelList shellLevel; +// shells.findHigherLevel(ctrs, minLevelField, shellLevel); +// +// forAll(minLevelField, i) +// { +// minLevelField[i] = max(minLevelField[i], shellLevel[i]); +// } +// } // Precalculate the refinement level for every element of the searchable -// surface. This is currently hardcoded for triSurfaceMesh only. +// surface. void Foam::refinementSurfaces::setMinLevelFields ( const shellSurfaces& shells @@ -522,58 +562,55 @@ void Foam::refinementSurfaces::setMinLevelFields { const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; - if (isA(geom)) + // Precalculation only makes sense if there are different regions + // (so different refinement levels possible) and there are some + // elements. Possibly should have 'enough' elements to have fine + // enough resolution but for now just make sure we don't catch e.g. + // searchableBox (size=6) + if (geom.regions().size() > 1 && geom.globalSize() > 10) { - const triSurfaceMesh& triMesh = refCast(geom); + // Representative local coordinates + const pointField ctrs = geom.coordinates(); - autoPtr minLevelFieldPtr - ( - new triSurfaceLabelField - ( - IOobject - ( - "minLevel", - triMesh.objectRegistry::time().timeName(), // instance - "triSurface", // local - triMesh, - IOobject::NO_READ, - IOobject::AUTO_WRITE - ), - triMesh, - dimless - ) - ); - triSurfaceLabelField& minLevelField = minLevelFieldPtr(); - - const triSurface& s = static_cast(triMesh); - - // Initialise fields to region wise minLevel - forAll(s, triI) + labelList minLevelField(ctrs.size(), -1); { - minLevelField[triI] = minLevel(surfI, s[triI].region()); + // Get the element index in a roundabout way. Problem is e.g. + // distributed surface where local indices differ from global + // ones (needed for getRegion call) + List info; + geom.findNearest + ( + ctrs, + scalarField(ctrs.size(), sqr(GREAT)), + info + ); + + // Get per element the region + labelList region; + geom.getRegion(info, region); + + // From the region get the surface-wise refinement level + forAll(minLevelField, i) + { + if (info[i].hit()) + { + minLevelField[i] = minLevel(surfI, region[i]); + } + } } // Find out if triangle inside shell with higher level - pointField fc(s.size()); - forAll(s, triI) - { - fc[triI] = s[triI].centre(s.points()); - } // What level does shell want to refine fc to? labelList shellLevel; - shells.findHigherLevel(fc, minLevelField, shellLevel); + shells.findHigherLevel(ctrs, minLevelField, shellLevel); - forAll(minLevelField, triI) + forAll(minLevelField, i) { - minLevelField[triI] = max - ( - minLevelField[triI], - shellLevel[triI] - ); + minLevelField[i] = max(minLevelField[i], shellLevel[i]); } - // Store field on triMesh - minLevelFieldPtr.ptr()->store(); + // Store minLevelField on surface + const_cast(geom).setField(minLevelField); } } } @@ -601,44 +638,89 @@ void Foam::refinementSurfaces::findHigherIntersection return; } + if (surfaces_.size() == 1) + { + // Optimisation: single segmented surface. No need to duplicate + // point storage. + + label surfI = 0; + + const searchableSurface& geom = allGeometry_[surfaces_[surfI]]; + + // Do intersection test + List intersectionInfo(start.size()); + geom.findLineAny(start, end, intersectionInfo); + + // See if a cached level field available + labelList minLevelField; + geom.getField(intersectionInfo, minLevelField); + bool haveLevelField = + ( + returnReduce(minLevelField.size(), sumOp