diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H index 25b0e6a5594dc6edd9af5051fcf86fe072772fb5..28911cf012de791dd9b5f7eb2389681529a38305 100644 --- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/createFields.H @@ -53,7 +53,6 @@ mesh.setFluxRequired(p.name()); suppressDict.add("cellMask", true); suppressDict.add("cellDisplacement", true); suppressDict.add("interpolatedCells", true); - suppressDict.add("cellInterpolationWeight", true); } const_cast<dictionary&> diff --git a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H index 5a3386be7786995ca9e4019c67dc579b0ea6fdf3..00838f5fea4ec9cdb363d547805095a6bf784ba8 100644 --- a/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H +++ b/applications/solvers/incompressible/pimpleFoam/overPimpleDyMFoam/interpolatedFaces.H @@ -239,11 +239,9 @@ forAll(isNeiInterpolatedFace, faceI) n1 = vector(-n.z(), 0, n.x()); } } + n1.normalise(); - n1 /= mag(n1); - - vector n2 = n ^ n1; - n2 /= mag(n2); + const vector n2 = normalised(n ^ n1); tensor rot = tensor diff --git a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C index f4a90a275bd769935a00932bac691d2d290e55fb..762f487ad16937cf14dc75ae9fa0ce87efaef143 100644 --- a/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C +++ b/applications/solvers/multiphase/compressibleInterFoam/twoPhaseMixtureThermo/twoPhaseMixtureThermo.C @@ -26,7 +26,7 @@ License #include "twoPhaseMixtureThermo.H" #include "gradientEnergyFvPatchScalarField.H" #include "mixedEnergyFvPatchScalarField.H" - +#include "collatedFileOperation.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -60,6 +60,10 @@ Foam::twoPhaseMixtureThermo::twoPhaseMixtureThermo T2.write(); } + // Note: we're writing files to be read in immediately afterwards. + // Avoid any thread-writing problems. + fileHandler().flush(); + thermo1_ = rhoThermo::New(U.mesh(), phase1Name()); thermo2_ = rhoThermo::New(U.mesh(), phase2Name()); diff --git a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C index 3ba6e82b720abc2e08810d4c0bbfd14a5965a139..a60ba21ae89f2cd23fb47e4c09a461b6f32148c7 100644 --- a/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C +++ b/applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/laserDTRM.C @@ -165,8 +165,7 @@ void Foam::radiation::laserDTRM::initialise() const scalar t = mesh_.time().value(); const vector lPosition = focalLaserPosition_->value(t); - vector lDir = laserDirection_->value(t); - lDir /= mag(lDir); + const vector lDir = normalised(laserDirection_->value(t)); if (debug) { @@ -175,7 +174,7 @@ void Foam::radiation::laserDTRM::initialise() } // Find a vector on the area plane. Normal to laser direction - vector rArea = vector::zero; + vector rArea = Zero; scalar magr = 0.0; { @@ -188,7 +187,7 @@ void Foam::radiation::laserDTRM::initialise() magr = mag(rArea); } } - rArea /= mag(rArea); + rArea.normalise(); scalar dr = focalLaserRadius_/ndr_; scalar dTheta = mathematical::twoPi/ndTheta_; diff --git a/applications/test/cubicEqn/Test-cubicEqn.C b/applications/test/cubicEqn/Test-cubicEqn.C index d0ab606666275b1dc614a7bed3fd31eab26509f5..aa0579a41940d23aac180a523d7b871a73f0d84f 100644 --- a/applications/test/cubicEqn/Test-cubicEqn.C +++ b/applications/test/cubicEqn/Test-cubicEqn.C @@ -51,7 +51,7 @@ void test(const Type& polynomialEqn, const scalar tol) case roots::real: v[i] = polynomialEqn.value(r[i]); e[i] = polynomialEqn.error(r[i]); - ok = ok && mag(v[i]) < tol*mag(e[i]); + ok = ok && mag(v[i]) <= tol*mag(e[i]); break; case roots::posInf: v[i] = + inf; @@ -79,8 +79,10 @@ void test(const Type& polynomialEqn, const scalar tol) int main() { - const int nTests = 1000000; - for (int t = 0; t < nTests; ++ t) + const scalar tol = 5; + + const label nTests = 1000000; + for (label t = 0; t < nTests; ++ t) { test ( @@ -91,11 +93,26 @@ int main() randomScalar(1e-50, 1e+50), randomScalar(1e-50, 1e+50) ), - 100 + tol ); } + Info << nTests << " random cubics tested" << endl; - Info << nTests << " cubics tested" << endl; + const label coeffMin = -9, coeffMax = 10, nCoeff = coeffMax - coeffMin; + for (label a = coeffMin; a < coeffMax; ++ a) + { + for (label b = coeffMin; b < coeffMax; ++ b) + { + for (label c = coeffMin; c < coeffMax; ++ c) + { + for (label d = coeffMin; d < coeffMax; ++ d) + { + test(cubicEqn(a, b, c, d), tol); + } + } + } + } + Info<< nCoeff*nCoeff*nCoeff*nCoeff << " integer cubics tested" << endl; return 0; } diff --git a/applications/test/faces/Test-faces.C b/applications/test/faces/Test-faces.C index 280cdf70069dcf59656e1fa355dbe3a90884ebd4..f8278272313f44b90e0d80a97715c1a10e0473be 100644 --- a/applications/test/faces/Test-faces.C +++ b/applications/test/faces/Test-faces.C @@ -31,19 +31,81 @@ Description #include "argList.H" #include "labelledTri.H" +#include "pointList.H" +#include "ListOps.H" using namespace Foam; + +template<class Face> +void faceInfo(const Face& f, const UList<point>& points) +{ + Info<< f + << " points:" << f.points(points) + << " normal:" << f.unitNormal(points); +} + + +template<class Face> +void testSign +( + const Face& f, + const UList<point>& points, + const UList<point>& testPoints +) +{ + for (const point& p : testPoints) + { + Info<< " point:" << p << " sign=" << f.sign(p, points) << nl; + } +} + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: int main(int argc, char *argv[]) { + pointList points1 + ({ + { 0, 0, 0}, + { -1, -1, 1}, + { 1, -1, -1}, + { 1, 1, -1}, + { -1, 1, 1} + }); + + pointList points2 = ListOps::create<point> + ( + points1, + [](const point& p){ return point(p.x(), p.y(), -p.z()); } + ); + + pointList testPoints + ({ + { -2, -2, -2}, + { -2, -2, 2}, + { 0, 0, 0}, + { 2, 2, -2}, + { 2, 2, 2} + }); + face f1{ 1, 2, 3, 4 }; - Info<< "face:" << f1 << nl; + Info<< "face:"; faceInfo(f1, points1); Info << nl; + testSign(f1, points1, testPoints); + + Info<< "face:"; faceInfo(f1, points2); Info << nl; + testSign(f1, points2, testPoints); + Info<< nl; triFace t1{ 1, 2, 3 }; - Info<< "triFace:" << t1 << nl; + Info<< "triFace:"; faceInfo(t1, points1); Info << nl; + testSign(t1, points1, testPoints); + + Info<< "triFace:"; faceInfo(t1, points2); Info << nl; + testSign(t1, points2, testPoints); + Info<< nl; + f1 = t1; Info<< "face:" << f1 << nl; diff --git a/applications/utilities/mesh/advanced/splitCells/splitCells.C b/applications/utilities/mesh/advanced/splitCells/splitCells.C index a1e5a8fc0199cc995bd47d94ddb7fef29131b3b5..ea6f01a8ef345bae25419cf7fdddfe9cd13aff50 100644 --- a/applications/utilities/mesh/advanced/splitCells/splitCells.C +++ b/applications/utilities/mesh/advanced/splitCells/splitCells.C @@ -135,8 +135,11 @@ bool largerAngle // Get cos between faceCentre and normal vector to determine in // which quadrant angle is. (Is correct for unwarped faces only!) // Correct for non-outwards pointing normal. - vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]); - c1c0 /= mag(c1c0) + VSMALL; + const vector c1c0 = + normalised + ( + mesh.faceCentres()[f1] - mesh.faceCentres()[f0] + ); scalar fcCosAngle = n0 & c1c0; @@ -319,8 +322,7 @@ bool splitCell { const edge& e = mesh.edges()[edgeI]; - vector eVec = e.vec(mesh.points()); - eVec /= mag(eVec); + const vector eVec = e.unitVec(mesh.points()); vector planeN = eVec ^ halfNorm; @@ -328,8 +330,7 @@ bool splitCell // halfway on fully regular meshes (since we want cuts // to be snapped to vertices) planeN += 0.01*halfNorm; - - planeN /= mag(planeN); + planeN.normalise(); // Define plane through edge plane cutPlane(mesh.points()[e.start()], planeN); @@ -410,11 +411,8 @@ void collectCuts label f0, f1; meshTools::getEdgeFaces(mesh, celli, edgeI, f0, f1); - vector n0 = faceAreas[f0]; - n0 /= mag(n0); - - vector n1 = faceAreas[f1]; - n1 /= mag(n1); + const vector n0 = normalised(faceAreas[f0]); + const vector n1 = normalised(faceAreas[f1]); if ( diff --git a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C index d0dd7c4a4e490fd26fb0b40159579e4347d9281e..51e92215775b25f51df82fe73a02cd3224549e8f 100644 --- a/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C +++ b/applications/utilities/mesh/conversion/gmshToFoam/gmshToFoam.C @@ -188,19 +188,17 @@ label findInternalFace(const primitiveMesh& mesh, const labelList& meshF) bool correctOrientation(const pointField& points, const cellShape& shape) { // Get centre of shape. - point cc(shape.centre(points)); + const point cc(shape.centre(points)); // Get outwards pointing faces. faceList faces(shape.faces()); - forAll(faces, i) + for (const face& f : faces) { - const face& f = faces[i]; - - vector n(f.normal(points)); + const vector areaNorm(f.areaNormal(points)); // Check if vector from any point on face to cc points outwards - if (((points[f[0]] - cc) & n) < 0) + if (((points[f[0]] - cc) & areaNorm) < 0) { // Incorrectly oriented return false; diff --git a/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C b/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C index 809f996afc3b1e1e645c34da30f0a13c76d6a113..48b55bd14b36fec4c8e3c42c6100fd65c192d3a9 100644 --- a/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C +++ b/applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C @@ -225,9 +225,10 @@ int main(int argc, char *argv[]) // Determine orientation of tri v.s. cell centre. point cc(cll.centre(points)); point fc(tri.centre(points)); - vector fn(tri.normal(points)); - if (((fc - cc) & fn) < 0) + const vector areaNorm(tri.areaNormal(points)); + + if (((fc - cc) & areaNorm) < 0) { // Boundary face points inwards. Flip. boundaryFaces[facei].flip(); diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C index 191ddf53a93d8086dcf95cdf90f9b0a2db07e239..0a0b9731229809c3ce3ef0fb2d4fc3244650ca56 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellAspectRatioControl/cellAspectRatioControl.C @@ -45,7 +45,7 @@ Foam::cellAspectRatioControl::cellAspectRatioControl ) { // Normalise the direction - aspectRatioDirection_ /= mag(aspectRatioDirection_) + SMALL; + aspectRatioDirection_.normalise(); Info<< nl << "Cell Aspect Ratio Control" << nl diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C index 25b8507bf53ddafd81cfdcce84c0e9c183426b44..6da2b5e6c0aa989e77584bcece6fd4240cb4d7cd 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C @@ -428,7 +428,7 @@ void Foam::searchableSurfaceControl::initialVertices if (mag(normals[0]) < SMALL) { - normals[0] = vector(1, 1, 1); + normals[0] = vector::one; } pointAlignment.reset(new triad(normals[0])); diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C index f92686b321cbbe8fbc102e03d0056e4d2a376f67..edb36f528c7b2c44589bd1104f7edd333ec683c1 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C @@ -1137,7 +1137,7 @@ void Foam::conformalVoronoiMesh::move() alignmentDirs[aA] = a + sign(dotProduct)*b; - alignmentDirs[aA] /= mag(alignmentDirs[aA]); + alignmentDirs[aA].normalise(); } } } diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C index fb9899438a7c766ac5d604e4e3b46a7850501f4b..00bad93cd888b6c2ae008d252f06b23fe1028e2c 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C @@ -1832,11 +1832,11 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches } vector correctNormal = calcSharedPatchNormal(vc1, vc2); - correctNormal /= mag(correctNormal); + correctNormal.normalise(); Info<< " cN " << correctNormal << endl; - vector fN = f.normal(pts); + vector fN = f.areaNormal(pts); if (mag(fN) < SMALL) { @@ -1844,7 +1844,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches continue; } - fN /= mag(fN); + fN.normalise(); Info<< " fN " << fN << endl; if ((fN & correctNormal) > 0) diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C index a4a4d0a8f09a78f7ce1a0e23f6d244bd082c59c9..9a653df699c59890a155865d38782d92d22c3801 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C @@ -220,8 +220,7 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating // << endl; // Calculate master point - vector masterPtVec(normalDir + nextNormalDir); - masterPtVec /= mag(masterPtVec) + SMALL; + const vector masterPtVec = normalised(normalDir + nextNormalDir); if ( diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C index 1f42b3e429d3fc0ad7473e352055c29e3f61025d..2a37049388999a8024defbb3a211357de420c2ff 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C @@ -481,10 +481,9 @@ void Foam::conformalVoronoiMesh::calcFaceZones norm ); - vector fN = faces[facei].normal(mesh.points()); - fN /= mag(fN) + SMALL; + const vector areaNorm = faces[facei].areaNormal(mesh.points()); - if ((norm[0] & fN) < 0) + if ((norm[0] & areaNorm) < 0) { flipMap[facei] = true; } diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C index 0afaae03f389e225b883b53389a161e519752be6..bc71e6d0f197262f86b0771395344a3f22b3fe3e 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C @@ -89,11 +89,9 @@ void Foam::conformationSurfaces::hasBoundedVolume Info<< " Index = " << surfaces_[s] << endl; Info<< " Offset = " << regionOffset_[s] << endl; - forAll(triSurf, sI) + for (const labelledTri& f : triSurf) { - const label patchID = - triSurf[sI].region() - + regionOffset_[s]; + const label patchID = f.region() + regionOffset_[s]; // Don't include baffle surfaces in the calculation if @@ -102,15 +100,15 @@ void Foam::conformationSurfaces::hasBoundedVolume != extendedFeatureEdgeMesh::BOTH ) { - sum += triSurf[sI].normal(surfPts); + sum += f.areaNormal(surfPts); } else { - nBaffles++; + ++nBaffles; } } Info<< " has " << nBaffles << " baffles out of " - << triSurf.size() << " triangles" << endl; + << triSurf.size() << " triangles" << nl; totalTriangles += triSurf.size(); } @@ -740,8 +738,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide surface.findNearest(sample, nearestDistSqr, info); - vector hitDir = info[0].rawPoint() - samplePts[i]; - hitDir /= mag(hitDir) + SMALL; + const vector hitDir = + normalised + ( + info[0].rawPoint() - samplePts[i] + ); pointIndexHit surfHit; label hitSurface; diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C index 278bc76f92b977d210f29474430d5a1b57645413..e0ce3b1ce0b7b6090490824513c985fc5436b192 100644 --- a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C +++ b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/CV2D.C @@ -558,7 +558,7 @@ void Foam::CV2D::newPoints() alignmentDirs[aA] = a + sign(dotProduct)*b; - alignmentDirs[aA] /= mag(alignmentDirs[aA]); + alignmentDirs[aA].normalise(); } } } @@ -845,7 +845,7 @@ void Foam::CV2D::newPoints() cd0 = vector2D(cd0.x() + cd0.y(), cd0.y() - cd0.x()); // Normalise the primary coordinate direction - cd0 /= mag(cd0); + cd0.normalise(); // Calculate the orthogonal coordinate direction vector2D cd1(-cd0.y(), cd0.x()); diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C index c746d95dc62021238e97951f1c270d1699a547b5..45a51ee384414cc0d66f2e4d1bebc9d99399ee64 100644 --- a/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C +++ b/applications/utilities/mesh/generation/foamyMesh/foamyQuadMesh/insertFeaturePoints.C @@ -123,7 +123,7 @@ void Foam::CV2D::insertFeaturePoints() vector2DField fpn = toPoint2D(feMesh.edgeNormals(edgeI)); vector2D cornerNormal = sum(fpn); - cornerNormal /= mag(cornerNormal); + cornerNormal.normalise(); if (debug) { diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 8d521b254ca5162053f8ec011bf6b2194975a363..cef09c6e94a7c2cf48cdc714c92c8509c4e27019 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -665,6 +665,13 @@ addLayersControls // Default is 0.5*featureAngle. Set to -180 always attempt extrusion //layerTerminationAngle 25; + // Optional: disable shrinking of edges that have one (or two) points + // on an extruded patch. + // Default is false to enable single/two cell thick channels to still + // have layers. In <=1806 this was true by default. On larger gaps it + // should have no effect. + //disableWallEdges true; + // Optional: at non-patched sides allow mesh to slip if extrusion // direction makes angle larger than slipFeatureAngle. Default is // 0.5*featureAngle. diff --git a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C index 40ee2c76c1e9b76f99d59929641b5e908abd81ee..74101860ca84296849a7cfeae8b3d2c88a974a42 100644 --- a/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C +++ b/applications/utilities/mesh/manipulation/mergeMeshes/mergeMeshes.C @@ -112,7 +112,7 @@ int main(int argc, char *argv[]) const word addRegion = args.lookupOrDefault<word> ( - "masterRegion", + "addRegion", polyMesh::defaultRegion ); diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C index cef8d2695b18f85299456b28a35549ad8b795bf4..ede9932cd367460046f8e8d18e386653bc2476bb 100644 --- a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C +++ b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C @@ -271,8 +271,8 @@ Foam::label Foam::meshDualiser::addInternalFace ); //pointField dualPoints(meshMod.points()); - //vector n(newFace.normal(dualPoints)); - //n /= mag(n); + //const vector n(newFace.unitNormal(dualPoints)); + // //Pout<< "Generated internal dualFace:" << dualFacei // << " verts:" << newFace // << " points:" << UIndirectList<point>(meshMod.points(), newFace) @@ -298,8 +298,8 @@ Foam::label Foam::meshDualiser::addInternalFace ); //pointField dualPoints(meshMod.points()); - //vector n(newFace.normal(dualPoints)); - //n /= mag(n); + //const vector n(newFace.unitNormal(dualPoints)); + // //Pout<< "Generated internal dualFace:" << dualFacei // << " verts:" << newFace // << " points:" << UIndirectList<point>(meshMod.points(), newFace) @@ -355,8 +355,8 @@ Foam::label Foam::meshDualiser::addBoundaryFace ); //pointField dualPoints(meshMod.points()); - //vector n(newFace.normal(dualPoints)); - //n /= mag(n); + //const vector n(newFace.unitNormal(dualPoints)); + // //Pout<< "Generated boundary dualFace:" << dualFacei // << " verts:" << newFace // << " points:" << UIndirectList<point>(meshMod.points(), newFace) diff --git a/applications/utilities/mesh/manipulation/setSet/setSet.C b/applications/utilities/mesh/manipulation/setSet/setSet.C index 7735602a13716f7bc5c662fc3bc25513672a6f6e..bc800b920aee19d401b79881caa4a5a4d8de8e14 100644 --- a/applications/utilities/mesh/manipulation/setSet/setSet.C +++ b/applications/utilities/mesh/manipulation/setSet/setSet.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -51,7 +51,6 @@ Description #include "faceZoneSet.H" #include "pointZoneSet.H" #include "timeSelector.H" -#include "collatedFileOperation.H" #include <stdio.h> @@ -283,6 +282,9 @@ void removeZone { WarningInFunction << "Failed writing zone " << setName << endl; } + zones.write(); + // Force flushing so we know it has finished writing + fileHandler().flush(); } } @@ -541,6 +543,8 @@ bool doCommand << "Failed writing set " << currentSet.objectPath() << endl; } + // Make sure writing is finished + fileHandler().flush(); } } } @@ -732,8 +736,6 @@ int main(int argc, char *argv[]) // Specific to topoSet/setSet: quite often we want to block upon writing // a set so we can immediately re-read it. So avoid use of threading // for set writing. - fileOperations::collatedFileOperation::maxThreadFileBufferSize = 0; - timeSelector::addOptions(true, false); #include "addRegionOption.H" argList::addBoolOption("noVTK", "Do not write VTK files"); diff --git a/applications/utilities/mesh/manipulation/topoSet/topoSet.C b/applications/utilities/mesh/manipulation/topoSet/topoSet.C index 225cd611d4aeae75c56ac1e48f71fcf667307e6e..9878b75b3c56394048cabd9d8fba4da17b6880d7 100644 --- a/applications/utilities/mesh/manipulation/topoSet/topoSet.C +++ b/applications/utilities/mesh/manipulation/topoSet/topoSet.C @@ -43,7 +43,6 @@ Description #include "faceZoneSet.H" #include "pointZoneSet.H" #include "IOdictionary.H" -#include "collatedFileOperation.H" using namespace Foam; @@ -91,6 +90,7 @@ void removeZone { WarningInFunction << "Failed writing zone " << setName << endl; } + fileHandler().flush(); } } @@ -199,11 +199,6 @@ polyMesh::readUpdateState meshReadUpdate(polyMesh& mesh) int main(int argc, char *argv[]) { - // Specific to topoSet/setSet: quite often we want to block upon writing - // a set so we can immediately re-read it. So avoid use of threading - // for set writing. - fileOperations::collatedFileOperation::maxThreadFileBufferSize = 0; - timeSelector::addOptions(true, false); #include "addDictOption.H" #include "addRegionOption.H" @@ -307,6 +302,7 @@ int main(int argc, char *argv[]) << "Failed writing set " << currentSet().objectPath() << endl; } + fileHandler().flush(); } break; @@ -347,6 +343,7 @@ int main(int argc, char *argv[]) << "Failed writing set " << currentSet().objectPath() << endl; } + fileHandler().flush(); } break; @@ -359,6 +356,7 @@ int main(int argc, char *argv[]) << "Failed writing set " << currentSet().objectPath() << endl; } + fileHandler().flush(); break; case topoSetSource::INVERT: @@ -370,6 +368,7 @@ int main(int argc, char *argv[]) << "Failed writing set " << currentSet().objectPath() << endl; } + fileHandler().flush(); break; case topoSetSource::REMOVE: diff --git a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C index a7132dae225e135ff405302716e27a9ca67a8ca1..9356e6d1a94462e840f11caab8edb60ceb73d117 100644 --- a/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C +++ b/applications/utilities/mesh/manipulation/transformPoints/transformPoints.C @@ -280,8 +280,8 @@ int main(int argc, char *argv[]) ( args.lookup("rotate")() ); - n1n2[0] /= mag(n1n2[0]); - n1n2[1] /= mag(n1n2[1]); + n1n2[0].normalise(); + n1n2[1].normalise(); const tensor rotT = rotationTensor(n1n2[0], n1n2[1]); diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C index 92cbf96ba0423eefc1360d369a97a90a3c332eb3..f11942fba950439eb6b8c2b89dbf8b5f19da97d4 100644 --- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C +++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C @@ -108,7 +108,6 @@ Usage #include "pointFieldDecomposer.H" #include "lagrangianFieldDecomposer.H" #include "decompositionModel.H" -#include "collatedFileOperation.H" #include "faCFD.H" #include "emptyFaPatch.H" @@ -523,12 +522,6 @@ int main(int argc, char *argv[]) // Decompose the mesh if (!decomposeFieldsOnly) { - // Disable buffering when writing mesh since we need to read - // it later on when decomposing the fields - float bufSz = - fileOperations::collatedFileOperation::maxThreadFileBufferSize; - fileOperations::collatedFileOperation::maxThreadFileBufferSize = 0; - mesh.decomposeMesh(); mesh.writeDecomposition(decomposeSets); @@ -587,8 +580,7 @@ int main(int argc, char *argv[]) << " for use in manual decomposition." << endl; } - fileOperations::collatedFileOperation::maxThreadFileBufferSize = - bufSz; + fileHandler().flush(); } @@ -1486,8 +1478,14 @@ int main(int argc, char *argv[]) fieldDecomposer.decomposeFields(areaScalarFields); fieldDecomposer.decomposeFields(areaVectorFields); - fieldDecomposer.decomposeFields(areaSphericalTensorFields); - fieldDecomposer.decomposeFields(areaSymmTensorFields); + fieldDecomposer.decomposeFields + ( + areaSphericalTensorFields + ); + fieldDecomposer.decomposeFields + ( + areaSymmTensorFields + ); fieldDecomposer.decomposeFields(areaTensorFields); fieldDecomposer.decomposeFields(edgeScalarFields); diff --git a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.C b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.C index 3dc56dc5d0fe2e968cc69b2ca5893c6b8640566e..1cf6e69aaf718c77be7f598239d47d94d0c4a38b 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoam.C @@ -194,7 +194,8 @@ int Foam::vtkPVFoam::setTime(const std::vector<double>& requestTimes) Time& runTime = dbPtr_(); - // Get times list + // Get times list. Flush first to force refresh. + fileHandler().flush(); instantList Times = runTime.times(); int nearestIndex = timeIndex_; @@ -301,11 +302,6 @@ Foam::vtkPVFoam::vtkPVFoam fileName FileName(vtkFileName); - // Make sure not to use the threaded version - it does not like - // being loaded as a shared library - static cleanup order is problematic. - // For now just disable the threaded writer. - fileOperations::collatedFileOperation::maxThreadFileBufferSize = 0; - // avoid argList and get rootPath/caseName directly from the file fileName fullCasePath(FileName.path()); @@ -729,6 +725,8 @@ std::vector<double> Foam::vtkPVFoam::findTimes(const bool skipZero) const if (dbPtr_.valid()) { const Time& runTime = dbPtr_(); + // Get times list. Flush first to force refresh. + fileHandler().flush(); instantList timeLst = runTime.times(); // find the first time for which this mesh appears to exist diff --git a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamUpdateInfo.C b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamUpdateInfo.C index 26ef1d660c2db83fff233abd025d8f09203237c8..5a5b0222eaa11d7d7bdc975a2535d2cfd1389483 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamUpdateInfo.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamUpdateInfo.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. ------------------------------------------------------------------------------- License @@ -225,6 +225,9 @@ void Foam::vtkPVFoam::updateInfoLagrangian // List of lagrangian objects across all times HashSet<fileName> names; + // Get times list. Flush first to force refresh. + fileHandler().flush(); + for (const instant& t : dbPtr_().times()) { names.insert @@ -697,6 +700,7 @@ void Foam::vtkPVFoam::updateInfoLagrangianFields // to some of the clouds. HashTable<wordHashSet> fields; + fileHandler().flush(); for (const instant& t : dbPtr_().times()) { for (const auto& cloudName : cloudNames) diff --git a/applications/utilities/postProcessing/noise/Allwmake b/applications/utilities/postProcessing/noise/Allwmake index 7452621cd23b522eab94d29ac635ec519d777629..2c9d265a5d535fe22c8e4d1fe0c912c7ce7a15ba 100755 --- a/applications/utilities/postProcessing/noise/Allwmake +++ b/applications/utilities/postProcessing/noise/Allwmake @@ -1,12 +1,13 @@ #!/bin/sh cd ${0%/*} || exit 1 # Run from this directory +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments # (for error catching) . $WM_PROJECT_DIR/wmake/scripts/have_fftw #------------------------------------------------------------------------------ if have_fftw then - wmake + wmake $targetType else echo "==> skip noise utility (no FFTW)" fi diff --git a/applications/utilities/preProcessing/boxTurb/Allwmake b/applications/utilities/preProcessing/boxTurb/Allwmake index 21e6762b4ae17ae2857043fe1ddadf0cca45acb0..820f2da5524b50cced8f461d559d4c1d1b927201 100755 --- a/applications/utilities/preProcessing/boxTurb/Allwmake +++ b/applications/utilities/preProcessing/boxTurb/Allwmake @@ -1,12 +1,13 @@ #!/bin/sh cd ${0%/*} || exit 1 # Run from this directory +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments # (for error catching) . $WM_PROJECT_DIR/wmake/scripts/have_fftw #------------------------------------------------------------------------------ if have_fftw then - wmake + wmake $targetType else echo "==> skip boxTurb utility (no FFTW)" fi diff --git a/applications/utilities/preProcessing/engineSwirl/createFields.H b/applications/utilities/preProcessing/engineSwirl/createFields.H index a77933c8f78a660bdbe94bf0606865a21f2f5d3c..bb8b2b0d23a15ba0d7b669315ef38372e1ef61cb 100644 --- a/applications/utilities/preProcessing/engineSwirl/createFields.H +++ b/applications/utilities/preProcessing/engineSwirl/createFields.H @@ -74,6 +74,6 @@ if (mag(yT) < SMALL) } //swirlAxis doesn't have to be of unit length. -xT /= mag(xT); -yT /= mag(yT); -zT /= mag(zT); +xT.normalise(); +yT.normalise(); +zT.normalise(); diff --git a/applications/utilities/surface/surfaceBooleanFeatures/Allwmake b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake index 0aae5a7cfb892026e95f8ea8e2f13deeb8b2a3b8..933ee039b9819872139af740cdb18924e56f1983 100755 --- a/applications/utilities/surface/surfaceBooleanFeatures/Allwmake +++ b/applications/utilities/surface/surfaceBooleanFeatures/Allwmake @@ -1,5 +1,6 @@ #!/bin/sh cd ${0%/*} || exit 1 # Run from this directory +. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments # (for error catching) . $WM_PROJECT_DIR/wmake/scripts/have_cgal #------------------------------------------------------------------------------ @@ -14,6 +15,6 @@ else export COMP_FLAGS="-DNO_CGAL" fi -wmake +wmake $targetType #------------------------------------------------------------------------------ diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C index 814cec9ce9070ba514ee83e603807d82d0d1b74b..d0cdb3a94415c2377a89427c13ab726b7799ff10 100644 --- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C +++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C @@ -317,11 +317,9 @@ label calcNormalDirection const vector& pointOnEdge ) { - vector cross = (normal ^ edgeDir); - cross /= mag(cross); + const vector cross = normalised(normal ^ edgeDir); - vector fC0tofE0 = faceCentre - pointOnEdge; - fC0tofE0 /= mag(fC0tofE0); + const vector fC0tofE0 = normalised(faceCentre - pointOnEdge); label nDir = ((cross & fC0tofE0) > 0.0 ? 1 : -1); diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C index 9e28c823f8ea6d12458a26173220c762fb498ce8..a62ecb2877eae30e71465b1d3d3f52cc1786a9bd 100644 --- a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C +++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C @@ -127,8 +127,7 @@ void greenRefine //{ // const edge& e = surf.edges()[edgeIndex]; -// vector eVec = e.vec(surf.localPoints()); -// eVec /= mag(eVec) + SMALL; +// const vector eVec = e.unitVec(surf.localPoints()); // const labelList& pEdges = surf.pointEdges()[pointIndex]; // @@ -136,8 +135,7 @@ void greenRefine // { // const edge& nearE = surf.edges()[pEdges[eI]]; -// vector nearEVec = nearE.vec(surf.localPoints()); -// nearEVec /= mag(nearEVec) + SMALL; +// const vector nearEVec = nearE.unitVec(surf.localPoints()); // const scalar dot = eVec & nearEVec; // const scalar minCos = degToRad(angle); diff --git a/applications/utilities/surface/surfaceInflate/surfaceInflate.C b/applications/utilities/surface/surfaceInflate/surfaceInflate.C index 83c5d147a0949fb0e3f393286acaa14cd1cf9c84..731d583065c6ccc836592f7111c718345575bcce 100644 --- a/applications/utilities/surface/surfaceInflate/surfaceInflate.C +++ b/applications/utilities/surface/surfaceInflate/surfaceInflate.C @@ -89,11 +89,8 @@ tmp<vectorField> calcVertexNormals(const triSurface& surf) { // Weighted average of normals of faces attached to the vertex // Weight = fA / (mag(e0)^2 * mag(e1)^2); - tmp<vectorField> tpointNormals - ( - new pointField(surf.nPoints(), Zero) - ); - vectorField& pointNormals = tpointNormals.ref(); + auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero); + auto& pointNormals = tpointNormals.ref(); const pointField& points = surf.points(); const labelListList& pointFaces = surf.pointFaces(); @@ -108,20 +105,20 @@ tmp<vectorField> calcVertexNormals(const triSurface& surf) const label faceI = pFaces[fI]; const triFace& f = surf[faceI]; - vector fN = f.normal(points); + vector areaNorm = f.areaNormal(points); scalar weight = calcVertexNormalWeight ( f, meshPoints[pI], - fN, + areaNorm, points ); - pointNormals[pI] += weight*fN; + pointNormals[pI] += weight * areaNorm; } - pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL; + pointNormals[pI].normalise(); } return tpointNormals; @@ -168,9 +165,9 @@ tmp<vectorField> calcPointNormals // Get average edge normal vector n = Zero; - forAll(eFaces, i) + for (const label facei : eFaces) { - n += s.faceNormals()[eFaces[i]]; + n += s.faceNormals()[facei]; } n /= eFaces.size(); @@ -192,7 +189,7 @@ tmp<vectorField> calcPointNormals { if (nNormals[pointI] > 0) { - pointNormals[pointI] /= mag(pointNormals[pointI]); + pointNormals[pointI].normalise(); } } } diff --git a/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C b/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C index 9d590fb9aa5d393d398880316e13963b4e39bd85..1ea1c8a175162eb2bcb6f0c0a30ed3a7400bb37f 100644 --- a/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C +++ b/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C @@ -565,16 +565,24 @@ void calcPointVecs if (face0I != -1) { label v0 = triSurfaceTools::oppositeVertex(surf, face0I, edgeI); - vector e0 = (points[v0] - points[e.start()]) ^ eVec; - e0 /= mag(e0); + const vector e0 = + normalised + ( + (points[v0] - points[e.start()]) ^ eVec + ); + midVec = e0; } if (face1I != -1) { label v1 = triSurfaceTools::oppositeVertex(surf, face1I, edgeI); - vector e1 = (points[e.start()] - points[v1]) ^ eVec; - e1 /= mag(e1); + const vector e1 = + normalised + ( + (points[e.start()] - points[v1]) ^ eVec + ); + midVec += e1; } @@ -895,8 +903,7 @@ int main(int argc, char *argv[]) { scalar minLen = minEdgeLen(surf, pointi); - vector n = borderPointVec[pointi]; - n /= mag(n); + const vector n = normalised(borderPointVec[pointi]); newPoints[newPointi] = newPoints[pointi] + 0.1 * minLen * n; } diff --git a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C index 10d389e85711ae04d1a23b4bf7151470a123c22b..51499c59ec95f2c9e9b47d6fbb82d771763b4791 100644 --- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C +++ b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C @@ -175,8 +175,8 @@ int main(int argc, char *argv[]) ( args.lookup("rotate")() ); - n1n2[0] /= mag(n1n2[0]); - n1n2[1] /= mag(n1n2[1]); + n1n2[0].normalise(); + n1n2[1].normalise(); const tensor rotT = rotationTensor(n1n2[0], n1n2[1]); diff --git a/src/OpenFOAM/containers/Bits/bitSet/bitSet.C b/src/OpenFOAM/containers/Bits/bitSet/bitSet.C index 941b5c4a083b99b2d7962c2dde2e168eadd6e013..3218145d86ae06ba0ac0f3c8c63e0e26cda15e53 100644 --- a/src/OpenFOAM/containers/Bits/bitSet/bitSet.C +++ b/src/OpenFOAM/containers/Bits/bitSet/bitSet.C @@ -159,9 +159,6 @@ Foam::bitSet& Foam::bitSet::orEq(const bitSet& other, const bool strict) << "Perform |= on dissimilar sized bitSets: " << size() << " vs. " << other.size() << nl; } - { - return *this; - } label minpos = -1; // Min trim point diff --git a/src/OpenFOAM/db/IOobjects/decomposedBlockData/decomposedBlockData.C b/src/OpenFOAM/db/IOobjects/decomposedBlockData/decomposedBlockData.C index 3ff5b05c2f9d326e452ab1f3d25ebf4bb784bee9..cca72be73ffa9acacf4b993ab64812d7216f6425 100644 --- a/src/OpenFOAM/db/IOobjects/decomposedBlockData/decomposedBlockData.C +++ b/src/OpenFOAM/db/IOobjects/decomposedBlockData/decomposedBlockData.C @@ -938,7 +938,7 @@ bool Foam::decomposedBlockData::writeBlocks label startProc = 1; label nSendProcs = nProcs-1; - while (nSendProcs > 0) + while (nSendProcs > 0 && startProc < nProcs) { nSendProcs = calcNumProcs ( @@ -952,7 +952,7 @@ bool Foam::decomposedBlockData::writeBlocks startProc ); - if (startProc == nProcs || nSendProcs == 0) + if (nSendProcs == 0) { break; } diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.C b/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.C index 616846413134533dcfab185ff0aaa7eeaf542526..7ce3d250d9b3a493a2c0e04f72bfdbaedc5b9db4 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.C +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.C @@ -265,7 +265,11 @@ void Foam::OFstreamCollator::waitForBufferSpace(const off_t wantedSize) const } } - if (totalSize == 0 || (totalSize+wantedSize) <= maxBufferSize_) + if + ( + totalSize == 0 + || (wantedSize >= 0 && (totalSize+wantedSize) <= maxBufferSize_) + ) { break; } @@ -354,7 +358,8 @@ bool Foam::OFstreamCollator::write IOstream::streamFormat fmt, IOstream::versionNumber ver, IOstream::compressionType cmp, - const bool append + const bool append, + const bool useThread ) { // Determine (on master) sizes to receive. Note: do NOT use thread @@ -374,7 +379,7 @@ bool Foam::OFstreamCollator::write Pstream::scatter(maxLocalSize, Pstream::msgType(), localComm_); } - if (maxBufferSize_ == 0 || maxLocalSize > maxBufferSize_) + if (!useThread || maxBufferSize_ == 0 || maxLocalSize > maxBufferSize_) { if (debug) { @@ -589,4 +594,20 @@ bool Foam::OFstreamCollator::write } +void Foam::OFstreamCollator::waitAll() +{ + // Wait for all buffer space to be available i.e. wait for all jobs + // to finish + if (Pstream::master(localComm_)) + { + if (debug) + { + Pout<< "OFstreamCollator : waiting for thread to have consumed all" + << endl; + } + waitForBufferSpace(-1); + } +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.H b/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.H index 92ef9068bfb5355a2e64d6ab81b58334d1ddc760..59463419921f2276139bceea78e535ab3c885cbb 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.H +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/OFstreamCollator.H @@ -206,8 +206,12 @@ public: IOstream::streamFormat, IOstream::versionNumber, IOstream::compressionType, - const bool append + const bool append, + const bool useThread = true ); + + //- Wait for all thread actions to have finished + void waitAll(); }; diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.C b/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.C index f0bbe6d6496642d274669f03d920a539646a6353..67d1985ec05ccfad6f6057370675ceaf483df212 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.C +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.C @@ -559,14 +559,32 @@ bool Foam::fileOperations::collatedFileOperation::writeObject } else { + // Re-check static maxThreadFileBufferSize variable to see + // if needs to use threading + bool useThread = (maxThreadFileBufferSize > 0); + if (debug) { Pout<< "collatedFileOperation::writeObject :" << " For object : " << io.name() - << " starting collating output to " << pathName << endl; + << " starting collating output to " << pathName + << " useThread:" << useThread << endl; + } + + if (!useThread) + { + writer_.waitAll(); } - threadedCollatedOFstream os(writer_, pathName, fmt, ver, cmp); + threadedCollatedOFstream os + ( + writer_, + pathName, + fmt, + ver, + cmp, + useThread + ); // If any of these fail, return (leave error handling to Ostream // class) @@ -593,6 +611,18 @@ bool Foam::fileOperations::collatedFileOperation::writeObject } } +void Foam::fileOperations::collatedFileOperation::flush() const +{ + if (debug) + { + Pout<< "collatedFileOperation::flush : clearing and waiting for thread" + << endl; + } + masterUncollatedFileOperation::flush(); + // Wait for thread to finish (note: also removes thread) + writer_.waitAll(); +} + Foam::word Foam::fileOperations::collatedFileOperation::processorsDir ( diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.H b/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.H index f03e2ae43a112ea652897728ff25668d4f4e902d..f24dc52c9df501738c9ce869e53ab9c8d327372b 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.H +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/collatedFileOperation.H @@ -155,6 +155,9 @@ public: // Other + //- Forcibly wait until all output done. Flush any cached data + virtual void flush() const; + //- Actual name of processors dir virtual word processorsDir(const IOobject&) const; diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.C b/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.C index df919680fbd92c71163012ec84d5cc0305b59dd4..c35ba0c1fa7919d2d5a064cc9d53d7e9fe7a8be6 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.C +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -35,13 +35,15 @@ Foam::threadedCollatedOFstream::threadedCollatedOFstream const fileName& pathName, streamFormat format, versionNumber version, - compressionType compression + compressionType compression, + const bool useThread ) : OStringStream(format, version), writer_(writer), pathName_(pathName), - compression_(compression) + compression_(compression), + useThread_(useThread) {} @@ -57,7 +59,8 @@ Foam::threadedCollatedOFstream::~threadedCollatedOFstream() IOstream::BINARY, version(), compression_, - false // append + false, // append + useThread_ ); } diff --git a/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.H b/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.H index 7efd8d0daf22508793d8cb5e700c4b5095be5db5..54217aa69ab095a32b51a1b7b30fa7fe4c7ee024 100644 --- a/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.H +++ b/src/OpenFOAM/global/fileOperations/collatedFileOperation/threadedCollatedOFstream.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2017 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2017-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -60,6 +60,8 @@ class threadedCollatedOFstream const IOstream::compressionType compression_; + const bool useThread_; + public: @@ -72,7 +74,8 @@ public: const fileName& pathname, streamFormat format=ASCII, versionNumber version=currentVersion, - compressionType compression=UNCOMPRESSED + compressionType compression=UNCOMPRESSED, + const bool useThread = true ); diff --git a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C index 927ca6f96c5f9103562857c35bcc91245dddcb31..03cf4bc80c5022a1322e95fe78adda2a91f16ac0 100644 --- a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C +++ b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.C @@ -784,7 +784,11 @@ Foam::IOobject Foam::fileOperation::findInstance for (; instanceI >= 0; --instanceI) { // Shortcut: if actual directory is the timeName we've already tested it - if (ts[instanceI].name() == startIO.instance()) + if + ( + ts[instanceI].name() == startIO.instance() + && ts[instanceI].name() != stopInstance + ) { continue; } @@ -991,6 +995,17 @@ Foam::label Foam::fileOperation::nProcs } +void Foam::fileOperation::flush() const +{ + if (debug) + { + Pout<< "fileOperation::flush : clearing processor directories cache" + << endl; + } + procsDirs_.clear(); +} + + Foam::fileName Foam::fileOperation::processorsCasePath ( const IOobject& io, diff --git a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H index b236e9eb45887cbdcaa99b39f4cc904fecafc79d..214ead34db443904b2496033b63afe7160910998 100644 --- a/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H +++ b/src/OpenFOAM/global/fileOperations/fileOperation/fileOperation.H @@ -514,6 +514,9 @@ public: virtual void setTime(const Time&) const {} + //- Forcibly wait until all output done. Flush any cached data + virtual void flush() const; + //- Generate path (like io.path) from root+casename with any // 'processorXXX' replaced by procDir (usually 'processsors') fileName processorsCasePath diff --git a/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C b/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C index b1d1f03130864892203f380b5a3cceca282f43b9..18ff87a6cb66d69710a6a444aa5143377c4d58c9 100644 --- a/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C +++ b/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.C @@ -2561,6 +2561,13 @@ Foam::fileOperations::masterUncollatedFileOperation::NewOFstream } +void Foam::fileOperations::masterUncollatedFileOperation::flush() const +{ + fileOperation::flush(); + times_.clear(); +} + + Foam::label Foam::fileOperations::masterUncollatedFileOperation::addWatch ( const fileName& fName diff --git a/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.H b/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.H index 03c100921d16b95109dddeddc1f462d88759e7e8..545f20041b7975d041400e0fe82770e58a8f4373 100644 --- a/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.H +++ b/src/OpenFOAM/global/fileOperations/masterUncollatedFileOperation/masterUncollatedFileOperation.H @@ -769,6 +769,9 @@ public: //- Callback for time change virtual void setTime(const Time&) const; + //- Forcibly wait until all output done. Flush any cached data + virtual void flush() const; + //- Return cached times const HashPtrTable<instantList>& times() const { diff --git a/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H b/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H index a9f3c395542b99c1dc24b0aa50d4c67d1cc1da53..22efb3d438a0c3f22d28d0b550bc5d566033d6d2 100644 --- a/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H +++ b/src/OpenFOAM/meshes/meshShapes/edge/edgeI.H @@ -421,10 +421,10 @@ inline Foam::vector Foam::edge::unitVec(const UList<point>& pts) const } #endif - Foam::vector v = pts[second()] - pts[first()]; - v /= ::Foam::mag(v) + VSMALL; + const vector v = (pts[second()] - pts[first()]); + const scalar s(Foam::mag(v)); - return v; + return s < ROOTVSMALL ? Zero : v/s; } diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index 15aaa8e20012074c5f623e9078d2ec6dd4ebf70a..986a65fcb897124d49435db8e802458975709ee2 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/src/OpenFOAM/meshes/meshShapes/face/face.H @@ -299,6 +299,20 @@ public: label& nearLabel ) const; + //- The sign for the side of the face plane the point is on, + //- using three evenly distributed face points for the estimated normal. + // Uses the supplied tolerance for rounding around zero. + // \return + // - 0: on plane + // - +1: front-side + // - -1: back-side + int sign + ( + const point& p, + const UList<point>& points, + const scalar tol = SMALL + ) const; + //- Return contact sphere diameter scalar contactSphereDiameter ( diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C index 1b250c401a4c3578266a13cdb09fe20000e10f48..bf5008c03fa80ecb9328d4e5b8a0ee63e76a1013 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C +++ b/src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C @@ -309,4 +309,23 @@ Foam::pointHit Foam::face::nearestPointClassify } +int Foam::face::sign +( + const point& p, + const UList<point>& points, + const scalar tol +) const +{ + // Take three points [0, 1/3, 2/3] from the face + // - assumes the face is not severely warped + + return triPointRef + ( + points[operator[](0)], + points[operator[](size()/3)], + points[operator[]((2*size())/3)] + ).sign(p, tol); +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H index aa6778ef6aca00774afb03a8c71272ac7273bdec..4b599dbd949eff28e6cfedb5d8ce4b7d2e556fdd 100644 --- a/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H +++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFace.H @@ -221,6 +221,19 @@ public: label& nearLabel ) const; + //- The sign for which side of the face plane the point is on. + // Uses the supplied tolerance for rounding around zero. + // \return + // - 0: on plane + // - +1: front-side + // - -1: back-side + inline int sign + ( + const point& p, + const UList<point>& points, + const scalar tol = SMALL + ) const; + //- Return number of edges inline label nEdges() const; diff --git a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H index f8afa79b6ebc36b1c886ecc8625509b16e308f57..0e7f1dae9d5292bd93c593fd0add9359f72bf60e 100644 --- a/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H +++ b/src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H @@ -332,6 +332,17 @@ inline Foam::pointHit Foam::triFace::nearestPointClassify } +inline int Foam::triFace::sign +( + const point& p, + const UList<point>& points, + const scalar tol +) const +{ + return this->tri(points).sign(p, tol); +} + + inline Foam::label Foam::triFace::nEdges() const { return 3; diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C index 68254970189c0ad9d691ef096d47f82731249f81..6b4c68d157f3f7e86ae3842cb70fd60dbb8ed499 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -27,18 +27,27 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::globalIndex::globalIndex +Foam::globalIndex::globalIndex(Istream& is) +{ + is >> offsets_; +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::globalIndex::reset ( const label localSize, const int tag, const label comm, const bool parallel ) -: - offsets_(Pstream::nProcs(comm)+1) { - labelList localSizes(Pstream::nProcs(comm), 0); + offsets_.resize(Pstream::nProcs(comm)+1); + + labelList localSizes(Pstream::nProcs(comm), Zero); localSizes[Pstream::myProcNo(comm)] = localSize; + if (parallel) { Pstream::gatherList(localSizes, tag, comm); @@ -47,9 +56,9 @@ Foam::globalIndex::globalIndex label offset = 0; offsets_[0] = 0; - for (label proci = 0; proci < Pstream::nProcs(comm); proci++) + for (label proci = 0; proci < Pstream::nProcs(comm); ++proci) { - label oldOffset = offset; + const label oldOffset = offset; offset += localSizes[proci]; if (offset < oldOffset) @@ -65,20 +74,21 @@ Foam::globalIndex::globalIndex } -Foam::globalIndex::globalIndex(const label localSize) -: - offsets_(Pstream::nProcs()+1) +void Foam::globalIndex::reset(const label localSize) { - labelList localSizes(Pstream::nProcs(), 0); + offsets_.resize(Pstream::nProcs()+1); + + labelList localSizes(Pstream::nProcs(), Zero); localSizes[Pstream::myProcNo()] = localSize; + Pstream::gatherList(localSizes, Pstream::msgType()); Pstream::scatterList(localSizes, Pstream::msgType()); label offset = 0; offsets_[0] = 0; - for (label proci = 0; proci < Pstream::nProcs(); proci++) + for (label proci = 0; proci < Pstream::nProcs(); ++proci) { - label oldOffset = offset; + const label oldOffset = offset; offset += localSizes[proci]; if (offset < oldOffset) @@ -94,12 +104,6 @@ Foam::globalIndex::globalIndex(const label localSize) } -Foam::globalIndex::globalIndex(Istream& is) -{ - is >> offsets_; -} - - // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // Foam::Istream& Foam::operator>>(Istream& is, globalIndex& gi) diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H index 40f8ad1741fc12744a23425d226a79c19aebcae8..f89d245192e9826700a214b5eae19c2359256982 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -78,11 +78,11 @@ public: //- Construct from local max size. // Does communication with default communicator and message tag. - globalIndex(const label localSize); + inline explicit globalIndex(const label localSize); //- Construct from local max size. // Does communication with given communicator and message tag - globalIndex + inline globalIndex ( const label localSize, const int tag, @@ -91,10 +91,10 @@ public: ); //- Copy construct from list of labels - inline globalIndex(const labelUList& offsets); + inline explicit globalIndex(const labelUList& offsets); //- Move construct from list of labels - inline globalIndex(labelList&& offsets); + inline explicit globalIndex(labelList&& offsets); //- Construct from Istream globalIndex(Istream& is); @@ -102,14 +102,38 @@ public: // Member Functions - // Edit + //- Check for null constructed or global sum == 0 + inline bool empty() const; - //- Change after construction - inline labelList& offsets(); + //- Const-access to the offsets + inline const labelList& offsets() const; + + + // Edit + + //- Write-access to the offsets, for changing after construction + inline labelList& offsets(); + + //- Reset from local size. + // Does communication with default communicator and message tag. + void reset(const label localSize); + + //- Reset from local size. + // Does communication with given communicator and message tag + void reset + ( + const label localSize, + const int tag, + const label comm, + const bool parallel // use parallel comms + ); // Queries relating to my processor (using world communicator) + //- My local start + inline label localStart() const; + //- My local size inline label localSize() const; @@ -132,6 +156,12 @@ public: //- Global sum of localSizes inline label size() const; + //- Start of proci data + inline label offset(const label proci) const; + + //- Start of proci data + inline label localStart(const label proci) const; + //- Size of proci data inline label localSize(const label proci) const; @@ -150,9 +180,6 @@ public: //- Which processor does global come from? Binary search. inline label whichProcID(const label i) const; - //- Start of proci data - inline label offset(const label proci) const; - // Other diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H index 801e69957330a38418f710219a7dc2c32fb6b56c..3adfb13ab09f4d2bedf968e85347f25fa6761392 100644 --- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H +++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H @@ -28,6 +28,28 @@ License // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +inline Foam::globalIndex::globalIndex(const label localSize) +: + globalIndex() +{ + reset(localSize); +} + + +inline Foam::globalIndex::globalIndex +( + const label localSize, + const int tag, + const label comm, + const bool parallel +) +: + globalIndex() +{ + reset(localSize, tag, comm, parallel); +} + + inline Foam::globalIndex::globalIndex(const labelUList& offsets) : offsets_(offsets) @@ -42,6 +64,18 @@ inline Foam::globalIndex::globalIndex(labelList&& offsets) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +inline bool Foam::globalIndex::empty() const +{ + return offsets_.empty() || offsets_.last() == 0; +} + + +inline const Foam::labelList& Foam::globalIndex::offsets() const +{ + return offsets_; +} + + inline Foam::labelList& Foam::globalIndex::offsets() { return offsets_; @@ -54,6 +88,18 @@ inline Foam::label Foam::globalIndex::offset(const label proci) const } +inline Foam::label Foam::globalIndex::localStart(const label proci) const +{ + return offsets_[proci]; +} + + +inline Foam::label Foam::globalIndex::localStart() const +{ + return localStart(Pstream::myProcNo()); +} + + inline Foam::label Foam::globalIndex::localSize(const label proci) const { return offsets_[proci+1] - offsets_[proci]; diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index f66b48555f3dc72127b7cb0ad549e4d6b1952852..6db8ff11fa0a37ec1cfb393b9786fe84a7112555 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2017 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2016-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -94,7 +94,7 @@ void Foam::polyMesh::calcDirections() const { reduce(emptyDirVec, sumOp<vector>()); - emptyDirVec /= mag(emptyDirVec); + emptyDirVec.normalise(); for (direction cmpt=0; cmpt<vector::nComponents; cmpt++) { @@ -118,7 +118,7 @@ void Foam::polyMesh::calcDirections() const { reduce(wedgeDirVec, sumOp<vector>()); - wedgeDirVec /= mag(wedgeDirVec); + wedgeDirVec.normalise(); for (direction cmpt=0; cmpt<vector::nComponents; cmpt++) { @@ -216,7 +216,13 @@ Foam::polyMesh::polyMesh(const IOobject& io) IOobject ( "boundary", - time().findInstance(meshDir(), "boundary"), + time().findInstance // allow 'newer' boundary file + ( + meshDir(), + "boundary", + IOobject::MUST_READ, + faces_.instance() + ), meshSubDir, *this, IOobject::MUST_READ, @@ -235,12 +241,13 @@ Foam::polyMesh::polyMesh(const IOobject& io) IOobject ( "pointZones", - time().findInstance - ( - meshDir(), - "pointZones", - IOobject::READ_IF_PRESENT - ), + //time().findInstance + //( + // meshDir(), + // "pointZones", + // IOobject::READ_IF_PRESENT + //), + faces_.instance(), meshSubDir, *this, IOobject::READ_IF_PRESENT, @@ -253,12 +260,13 @@ Foam::polyMesh::polyMesh(const IOobject& io) IOobject ( "faceZones", - time().findInstance - ( - meshDir(), - "faceZones", - IOobject::READ_IF_PRESENT - ), + //time().findInstance + //( + // meshDir(), + // "faceZones", + // IOobject::READ_IF_PRESENT + //), + faces_.instance(), meshSubDir, *this, IOobject::READ_IF_PRESENT, @@ -271,12 +279,13 @@ Foam::polyMesh::polyMesh(const IOobject& io) IOobject ( "cellZones", - time().findInstance - ( - meshDir(), - "cellZones", - IOobject::READ_IF_PRESENT - ), + //time().findInstance + //( + // meshDir(), + // "cellZones", + // IOobject::READ_IF_PRESENT + //), + faces_.instance(), meshSubDir, *this, IOobject::READ_IF_PRESENT, @@ -1375,7 +1384,7 @@ bool Foam::polyMesh::pointInCell vector proj = p - faceTri.centre(); - if ((faceTri.normal() & proj) > 0) + if ((faceTri.areaNormal() & proj) > 0) { return false; } @@ -1405,7 +1414,7 @@ bool Foam::polyMesh::pointInCell vector proj = p - faceTri.centre(); - if ((faceTri.normal() & proj) > 0) + if ((faceTri.areaNormal() & proj) > 0) { return false; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C index 4257571d72a5668a8585f01f98c093f939154e40..845977c910fd12a804053bfc368659cc9546d7ff 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C @@ -60,9 +60,9 @@ Foam::label Foam::cyclicPolyPatch::findMaxArea forAll(faces, facei) { - scalar areaSqr = magSqr(faces[facei].normal(points)); + scalar areaSqr = magSqr(faces[facei].areaNormal(points)); - if (areaSqr > maxAreaSqr) + if (maxAreaSqr < areaSqr) { maxAreaSqr = areaSqr; maxI = facei; @@ -81,7 +81,7 @@ void Foam::cyclicPolyPatch::calcTransforms() vectorField half0Areas(half0.size()); forAll(half0, facei) { - half0Areas[facei] = half0[facei].normal(half0.points()); + half0Areas[facei] = half0[facei].areaNormal(half0.points()); } // Half1 @@ -89,7 +89,7 @@ void Foam::cyclicPolyPatch::calcTransforms() vectorField half1Areas(half1.size()); forAll(half1, facei) { - half1Areas[facei] = half1[facei].normal(half1.points()); + half1Areas[facei] = half1[facei].areaNormal(half1.points()); } calcTransforms @@ -260,19 +260,17 @@ void Foam::cyclicPolyPatch::calcTransforms // use calculated normals. vector n0 = findFaceMaxRadius(half0Ctrs); vector n1 = -findFaceMaxRadius(half1Ctrs); - n0 /= mag(n0) + VSMALL; - n1 /= mag(n1) + VSMALL; + n0.normalise(); + n1.normalise(); if (debug) { - scalar theta = radToDeg(acos(n0 & n1)); - Pout<< "cyclicPolyPatch::calcTransforms :" << " patch:" << name() << " Specified rotation :" << " n0:" << n0 << " n1:" << n1 - << " swept angle: " << theta << " [deg]" - << endl; + << " swept angle: " << radToDeg(acos(n0 & n1)) + << " [deg]" << endl; } // Extended tensor from two local coordinate systems calculated @@ -420,18 +418,17 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors { vector n0 = findFaceMaxRadius(half0Ctrs); vector n1 = -findFaceMaxRadius(half1Ctrs); - n0 /= mag(n0) + VSMALL; - n1 /= mag(n1) + VSMALL; + n0.normalise(); + n1.normalise(); if (debug) { - scalar theta = radToDeg(acos(n0 & n1)); - Pout<< "cyclicPolyPatch::getCentresAndAnchors :" << " patch:" << name() << " Specified rotation :" << " n0:" << n0 << " n1:" << n1 - << " swept angle: " << theta << " [deg]" + << " swept angle: " + << radToDeg(acos(n0 & n1)) << " [deg]" << endl; } @@ -499,12 +496,10 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors // Determine the face with max area on both halves. These // two faces are used to determine the transformation tensors label max0I = findMaxArea(pp0.points(), pp0); - vector n0 = pp0[max0I].normal(pp0.points()); - n0 /= mag(n0) + VSMALL; + const vector n0 = pp0[max0I].unitNormal(pp0.points()); label max1I = findMaxArea(pp1.points(), pp1); - vector n1 = pp1[max1I].normal(pp1.points()); - n1 /= mag(n1) + VSMALL; + const vector n1 = pp1[max1I].unitNormal(pp1.points()); if (mag(n0 & n1) < 1-matchTolerance()) { diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C index 80714b0bb471e8e7ace2b5a72da25eed90406fce..b63e957266e90b19f54294a1ffbf47a5e4f3ed0d 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/oldCyclic/oldCyclicPolyPatch.C @@ -92,9 +92,9 @@ Foam::label Foam::oldCyclicPolyPatch::findMaxArea forAll(faces, facei) { - scalar areaSqr = magSqr(faces[facei].normal(points)); + scalar areaSqr = magSqr(faces[facei].areaNormal(points)); - if (areaSqr > maxAreaSqr) + if (maxAreaSqr < areaSqr) { maxAreaSqr = areaSqr; maxI = facei; @@ -308,10 +308,17 @@ void Foam::oldCyclicPolyPatch::getCentresAndAnchors label face0 = getConsistentRotationFace(half0Ctrs); label face1 = getConsistentRotationFace(half1Ctrs); - vector n0 = ((half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_); - vector n1 = ((half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_); - n0 /= mag(n0) + VSMALL; - n1 /= mag(n1) + VSMALL; + const vector n0 = + normalised + ( + (half0Ctrs[face0] - rotationCentre_) ^ rotationAxis_ + ); + + const vector n1 = + normalised + ( + (half1Ctrs[face1] - rotationCentre_) ^ -rotationAxis_ + ); if (debug) { @@ -359,12 +366,10 @@ void Foam::oldCyclicPolyPatch::getCentresAndAnchors // Determine the face with max area on both halves. These // two faces are used to determine the transformation tensors label max0I = findMaxArea(pp.points(), half0Faces); - vector n0 = half0Faces[max0I].normal(pp.points()); - n0 /= mag(n0) + VSMALL; + const vector n0 = half0Faces[max0I].unitNormal(pp.points()); label max1I = findMaxArea(pp.points(), half1Faces); - vector n1 = half1Faces[max1I].normal(pp.points()); - n1 /= mag(n1) + VSMALL; + const vector n1 = half1Faces[max1I].unitNormal(pp.points()); if (mag(n0 & n1) < 1-matchTolerance()) { diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C index 9a69a136c9d84315b365c10d646f571066ab927e..d1a6160d624f9238f47bf2a40fe9434c396154b8 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C +++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/wedge/wedgePolyPatch.C @@ -86,7 +86,7 @@ void Foam::wedgePolyPatch::calcGeometry(PstreamBuffers&) sign(n_.y())*(max(mag(n_.y()), 0.5) - 0.5), sign(n_.z())*(max(mag(n_.z()), 0.5) - 0.5) ); - centreNormal_ /= mag(centreNormal_); + centreNormal_.normalise(); cosAngle_ = centreNormal_ & n_; diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C index 4f70b4ca8515f91848acd9944230812b3c6942df..0abf2c14317e132cd2a273f9924bf833af7bb699 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C +++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsSortEdges.C @@ -63,8 +63,7 @@ Foam::PatchTools::sortedEdgeFaces const point& edgePt = localPoints[e.start()]; - vector e2 = e.vec(localPoints); - e2 /= mag(e2) + VSMALL; + const vector e2 = e.unitVec(localPoints); // Get the vertex on 0th face that forms a vector with the first // edge point that has the largest angle with the edge @@ -77,8 +76,11 @@ Foam::PatchTools::sortedEdgeFaces { if (f0[fpI] != e.start()) { - vector faceEdgeDir = localPoints[f0[fpI]] - edgePt; - faceEdgeDir /= mag(faceEdgeDir) + VSMALL; + const vector faceEdgeDir = + normalised + ( + localPoints[f0[fpI]] - edgePt + ); const scalar angle = e2 & faceEdgeDir; @@ -92,11 +94,10 @@ Foam::PatchTools::sortedEdgeFaces // Get vector normal both to e2 and to edge from opposite vertex // to edge (will be x-axis of our coordinate system) - vector e0 = e2 ^ maxAngleEdgeDir; - e0 /= mag(e0) + VSMALL; + const vector e0 = normalised(e2 ^ maxAngleEdgeDir); // Get y-axis of coordinate system - vector e1 = e2 ^ e0; + const vector e1 = e2 ^ e0; SortableList<scalar> faceAngles(faceNbs.size()); @@ -116,8 +117,11 @@ Foam::PatchTools::sortedEdgeFaces { if (f[fpI] != e.start()) { - vector faceEdgeDir = localPoints[f[fpI]] - edgePt; - faceEdgeDir /= mag(faceEdgeDir) + VSMALL; + const vector faceEdgeDir = + normalised + ( + localPoints[f[fpI]] - edgePt + ); const scalar angle = e2 & faceEdgeDir; @@ -129,8 +133,7 @@ Foam::PatchTools::sortedEdgeFaces } } - vector vec = e2 ^ maxAngleEdgeDir; - vec /= mag(vec) + VSMALL; + const vector vec = normalised(e2 ^ maxAngleEdgeDir); faceAngles[nbI] = pseudoAngle ( diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C index 50e4222be0a54fee5b911ae385cacf592d34ce6f..9a76270169c9c7e4de2a84777819fef9cc5d2318 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C +++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C @@ -281,12 +281,12 @@ calcPointNormals() const const labelList& curFaces = pf[pointi]; - forAll(curFaces, facei) + for (const label facei : curFaces) { - curNormal += faceUnitNormals[curFaces[facei]]; + curNormal += faceUnitNormals[facei]; } - curNormal /= mag(curNormal) + VSMALL; + curNormal.normalise(); } if (debug) @@ -425,7 +425,7 @@ calcFaceAreas() const forAll(n, facei) { - n[facei] = this->operator[](facei).normal(points_); + n[facei] = this->operator[](facei).areaNormal(points_); } if (debug) @@ -472,8 +472,7 @@ calcFaceNormals() const forAll(n, facei) { - n[facei] = this->operator[](facei).normal(points_); - n[facei] /= mag(n[facei]) + VSMALL; + n[facei] = this->operator[](facei).unitNormal(points_); } if (debug) diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index 5c860a0c2c4c479b6961570c3716989dc6d999ca..3fb4b324da972096ab85b7888813e9d5eae056e1 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -78,8 +78,7 @@ Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness { vector Cpf = fCtrs[facei] - ownCc; - vector normal = fAreas[facei]; - normal /= mag(normal) + ROOTVSMALL; + vector normal = normalised(fAreas[facei]); vector d = normal*(normal & Cpf); diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H index c20edf59a087a3142517ac88a3ca34488199fdbb..4235bddd607585e6f3c4032cbaaa7553ff9bdff5 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H +++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H @@ -34,9 +34,8 @@ SourceFiles #ifndef cut_H #define cut_H -#include "FixedList.H" -#include "nil.H" #include "plane.H" +#include "FixedList.H" #include "tetPointRef.H" #include "triPointRef.H" #include "zero.H" @@ -55,8 +54,6 @@ namespace cut template<class Type> class uniformOp { -private: - //- Data Type data_; @@ -92,7 +89,7 @@ public: class noOp : - public uniformOp<nil> + public uniformOp<zero::null> { public: @@ -109,9 +106,9 @@ public: {} //- Construct from base - noOp(const uniformOp<nil>& op) + noOp(const uniformOp<zero::null>& op) : - uniformOp<nil>(op) + uniformOp<zero::null>(op) {} @@ -138,7 +135,7 @@ public: class areaOp : - public uniformOp<nil> + public uniformOp<zero::null> { public: @@ -155,9 +152,9 @@ public: {} //- Construct from base - areaOp(const uniformOp<nil>& op) + areaOp(const uniformOp<zero::null>& op) : - uniformOp<nil>(op) + uniformOp<zero::null>(op) {} @@ -183,7 +180,7 @@ public: class volumeOp : - public uniformOp<nil> + public uniformOp<zero::null> { public: @@ -200,9 +197,9 @@ public: {} //- Construct from base - volumeOp(const uniformOp<nil>& op) + volumeOp(const uniformOp<zero::null>& op) : - uniformOp<nil>(op) + uniformOp<zero::null>(op) {} @@ -315,7 +312,7 @@ public: template<unsigned Size> class listOp : - public uniformOp<nil> + public uniformOp<zero::null> { public: @@ -356,9 +353,9 @@ public: {} //- Construct from base - listOp(const uniformOp<nil>& op) + listOp(const uniformOp<zero::null>& op) : - uniformOp<nil>(op) + uniformOp<zero::null>(op) {} @@ -379,8 +376,6 @@ public: typedef listOp<3> listTriOp; - - typedef listOp<4> listTetOp; @@ -437,8 +432,7 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // //- Trait to determine the result of the addition of two operations -template<class AheadOp, class BehindOp> -class opAddResult; +template<class AheadOp, class BehindOp> class opAddResult; template<class Op> class opAddResult<Op, Op> diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/line.H b/src/OpenFOAM/meshes/primitiveShapes/line/line.H index a3ccc2319f24d7254408f18d1ed083dcf747432a..f07522832431971785aa31c84d52dd160cabc614 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/line/line.H +++ b/src/OpenFOAM/meshes/primitiveShapes/line/line.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -98,24 +98,27 @@ public: inline line(Istream&); - // Member functions + // Member Functions - // Access + // Access //- Return first point inline PointRef first() const; - //- Return second point + //- Return second (last) point inline PointRef second() const; + //- Return last (second) point + inline PointRef last() const; + //- Return first point inline PointRef start() const; - //- Return second point + //- Return second (last) point inline PointRef end() const; - // Properties + // Properties //- Return centre (centroid) inline Point centre() const; diff --git a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H index b2ee4d40ea132bfc26517b1b77617546e237aa19..01b018e50fdf86c901d221a89ffc793b0da8f144 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/line/lineI.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -23,6 +23,7 @@ License \*---------------------------------------------------------------------------*/ +#include "zero.H" #include "IOstreams.H" // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -70,6 +71,13 @@ inline PointRef Foam::line<Point, PointRef>::second() const } +template<class Point, class PointRef> +inline PointRef Foam::line<Point, PointRef>::last() const +{ + return b_; +} + + template<class Point, class PointRef> inline PointRef Foam::line<Point, PointRef>::start() const { @@ -107,10 +115,10 @@ inline Point Foam::line<Point, PointRef>::vec() const template<class Point, class PointRef> inline Point Foam::line<Point, PointRef>::unitVec() const { - Point v = b_ - a_; - v /= ::Foam::mag(v) + VSMALL; + const Point v = (b_ - a_); + const scalar s(::Foam::mag(v)); - return v; + return s < ROOTVSMALL ? Zero : v/s; } diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index d8de8e96096ec011914f3d884f47bea04ba81565..c5fbb32c3914926b933ade10dff7f133e6ffd92e 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -323,6 +323,14 @@ public: pointHit& edgePoint ) const; + //- The sign for which side of the face plane the point is on. + // Uses the supplied tolerance for rounding around zero. + // \return + // - 0: on plane + // - +1: front-side + // - -1: back-side + inline int sign(const point& p, const scalar tol = SMALL) const; + //- Decompose triangle into triangles above and below plane template<class AboveOp, class BelowOp> inline void sliceWithPlane diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index ffd9c20cab727505a810a4812408114033b98dec..00ab9eceee6f61f2634d5c90537a945cbea0f1c6 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -815,6 +815,19 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint } +template<class Point, class PointRef> +inline int Foam::triangle<Point, PointRef>::sign +( + const point& p, + const scalar tol +) const +{ + const scalar dist = ((p - a_) & unitNormal()); + + return ((dist < -tol) ? -1 : (dist > tol) ? +1 : 0); +} + + template<class Point, class PointRef> inline void Foam::triangle<Point, PointRef>::dummyOp::operator() ( diff --git a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H b/src/OpenFOAM/primitives/Vector2D/Vector2DI.H index 3c5fcdb7d31b6fa467539af0d7f00f473ad9cb1d..95c0514943e0064f746d5b7ba5a6fcfc76d4c86b 100644 --- a/src/OpenFOAM/primitives/Vector2D/Vector2DI.H +++ b/src/OpenFOAM/primitives/Vector2D/Vector2DI.H @@ -105,6 +105,8 @@ inline Foam::Vector2D<Cmpt>& Foam::Vector2D<Cmpt>::normalise() { *this /= s; } + + return *this; } diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C index c46b794c38ea16d67c8dd3dbc8f6a697d3567fa5..3aad674001853afde2e1b543303909ba736d0f16 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqn.C @@ -81,7 +81,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const // This is assumed not to over- or under-flow. If it does, all bets are off. const scalar p = c*a - b*b/3; - const scalar q = b*b*b*(2.0/27.0) - b*c*a/3 + d*a*a; + const scalar q = b*b*b*scalar(2)/27 - b*c*a/3 + d*a*a; const scalar disc = p*p*p/27 + q*q/4; // How many roots of what types are available? @@ -96,7 +96,7 @@ Foam::Roots<3> Foam::cubicEqn::roots() const if (oneReal) { - const Roots<1> r = linearEqn(- a, b/3).roots(); + const Roots<1> r = linearEqn(a, b/3).roots(); return Roots<3>(r.type(0), r[0]); } else if (twoReal) diff --git a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H index da5fc3cba14ef44e5f226d1bcdaef220455d540a..cb528b6be4a65a8ed36c6fe3b91bebf79a480d01 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H +++ b/src/OpenFOAM/primitives/polynomialEqns/cubicEqn/cubicEqnI.H @@ -114,7 +114,10 @@ inline Foam::scalar Foam::cubicEqn::derivative(const scalar x) const inline Foam::scalar Foam::cubicEqn::error(const scalar x) const { - return mag(SMALL*x*(x*(x*3*a() + 2*b()) + c())); + return + SMALL*magSqr(x)*(mag(x*a()) + mag(b())) + + SMALL*mag(x)*(mag(x*(x*a() + b())) + mag(c())) + + SMALL*(mag(x*(x*(x*a() + b()) + c())) + mag(d())); } diff --git a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H index 17bfdf452da1f4e0a576b199770a0b5fbe20ea86..f5c4f2b7c94cdb859ac550bcb107c74588aa817e 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H +++ b/src/OpenFOAM/primitives/polynomialEqns/linearEqn/linearEqnI.H @@ -82,7 +82,7 @@ inline Foam::scalar Foam::linearEqn::derivative(const scalar x) const inline Foam::scalar Foam::linearEqn::error(const scalar x) const { - return mag(SMALL*x*a()); + return SMALL*(mag(x*a()) + mag(b())); } diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C index d84c2928c25b133d4712829768c11439cf8179a5..d2ac0ea8a1284f29b508f13ff6e4c5066af73f26 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqn.C @@ -71,7 +71,7 @@ Foam::Roots<2> Foam::quadraticEqn::roots() const if (oneReal) { - const Roots<1> r = linearEqn(- a, b/2).roots(); + const Roots<1> r = linearEqn(a, b/2).roots(); return Roots<2>(r, r); } else if (twoReal) diff --git a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H index 7672a74d72458e4b70d5b2912a79c249bb0a8035..15ffe4e4e227a54cdf827a4b108a0d76f5e0f314 100644 --- a/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H +++ b/src/OpenFOAM/primitives/polynomialEqns/quadraticEqn/quadraticEqnI.H @@ -100,7 +100,9 @@ inline Foam::scalar Foam::quadraticEqn::derivative(const scalar x) const inline Foam::scalar Foam::quadraticEqn::error(const scalar x) const { - return mag(SMALL*x*(x*2*a() + b())); + return + SMALL*mag(x)*(mag(x*a()) + mag(b())) + + SMALL*(mag(x*(x*a() + b())) + mag(c())); } diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C index 3fac6847e7d38d092f59f4158f3abbf9875e4aad..f87c8c9b394ea5a2045db945187d94a355d9a96e 100644 --- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C +++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/externalWallHeatFluxTemperature/externalWallHeatFluxTemperatureFvPatchScalarField.C @@ -215,6 +215,7 @@ externalWallHeatFluxTemperatureFvPatchScalarField if (qrName_ != "none") { + qrPrevious_.setSize(mapper.size()); qrPrevious_.map(ptf.qrPrevious_, mapper); } } diff --git a/src/conversion/polyDualMesh/polyDualMesh.C b/src/conversion/polyDualMesh/polyDualMesh.C index b0c1f2eaead38c3ef52c1e89c6c23aac8d0b7fe9..9b325c94157f8c1672cb026b57e1b86a9e0e6d8e 100644 --- a/src/conversion/polyDualMesh/polyDualMesh.C +++ b/src/conversion/polyDualMesh/polyDualMesh.C @@ -1006,8 +1006,9 @@ void Foam::polyDualMesh::calcDual { // Check orientation. const face& f = dynDualFaces.last(); - vector n = f.normal(dualPoints); - if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0) + const vector areaNorm = f.areaNormal(dualPoints); + + if (((mesh.points()[owner] - dualPoints[f[0]]) & areaNorm) > 0) { WarningInFunction << " on boundary edge:" << edgeI @@ -1121,8 +1122,9 @@ void Foam::polyDualMesh::calcDual { // Check orientation. const face& f = dynDualFaces.last(); - vector n = f.normal(dualPoints); - if (((mesh.points()[owner] - dualPoints[f[0]]) & n) > 0) + const vector areaNorm = f.areaNormal(dualPoints); + + if (((mesh.points()[owner] - dualPoints[f[0]]) & areaNorm) > 0) { WarningInFunction << " on internal edge:" << edgeI diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C index 9d8ad5de9c39d8968b28a19cb4e5af38b6f274c5..d6c18d1cb89acc33f307a118e92b005cec03c688 100644 --- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C +++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C @@ -389,8 +389,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::fvMeshDistribute::generateTestField const fvMesh& mesh ) { - vector testNormal(1, 1, 1); - testNormal /= mag(testNormal); + const vector testNormal = normalised(vector::one); tmp<surfaceScalarField> tfld ( @@ -440,8 +439,7 @@ void Foam::fvMeshDistribute::testField(const surfaceScalarField& fld) { const fvMesh& mesh = fld.mesh(); - vector testNormal(1, 1, 1); - testNormal /= mag(testNormal); + const vector testNormal = normalised(vector::one); const surfaceVectorField n(mesh.Sf()/mesh.magSf()); diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C index 1dbbc324efb4cba52e89f8d117b19090efe908c6..edb8209d6f5f4e9a43a1cb7fa98b0f342d040d43 100644 --- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C +++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C @@ -1363,28 +1363,26 @@ bool Foam::cellCuts::loopAnchorConsistent // Create identity face for ease of calculation of normal etc. face f(identity(loopPts.size())); - vector normal = f.normal(loopPts); - point ctr = f.centre(loopPts); + const vector areaNorm = f.areaNormal(loopPts); + const point ctr = f.centre(loopPts); // Get average position of anchor points. vector avg(Zero); - forAll(anchorPoints, ptI) + for (const label pointi : anchorPoints) { - avg += mesh().points()[anchorPoints[ptI]]; + avg += mesh().points()[pointi]; } avg /= anchorPoints.size(); - if (((avg - ctr) & normal) > 0) + if (((avg - ctr) & areaNorm) > 0) { return true; } - else - { - return false; - } + + return false; } diff --git a/src/dynamicMesh/meshCut/cellLooper/geomCellLooper.C b/src/dynamicMesh/meshCut/cellLooper/geomCellLooper.C index 5528d9dad5348e0950abee38367fb3a5fb4b78d0..fbe7c6469f0832f3a18ddc028698edbad58cf99a 100644 --- a/src/dynamicMesh/meshCut/cellLooper/geomCellLooper.C +++ b/src/dynamicMesh/meshCut/cellLooper/geomCellLooper.C @@ -159,9 +159,7 @@ void Foam::geomCellLooper::getBase(const vector& n, vector& e0, vector& e1) // Use component normal to n as base vector. - e0 = base - nComp*n; - - e0 /= mag(e0) + VSMALL; + e0 = normalised(base - nComp * n); e1 = n ^ e0; @@ -404,8 +402,7 @@ bool Foam::geomCellLooper::cut forAll(sortedAngles, i) { - vector toCtr(loopPoints[i] - ctr); - toCtr /= mag(toCtr); + const vector toCtr = normalised(loopPoints[i] - ctr); sortedAngles[i] = pseudoAngle(e0, e1, toCtr); } diff --git a/src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C b/src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C index b7637869d601bffcfa2562cc796679aa2f1b54a5..483960ca24eca69e8288403c806becc89c6c069b 100644 --- a/src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C +++ b/src/dynamicMesh/meshCut/cellLooper/topoCellLooper.C @@ -262,8 +262,7 @@ Foam::label Foam::topoCellLooper::getAlignedNonFeatureEdge vector e0 = mesh().points()[e.start()] - ctr; vector e1 = mesh().points()[e.end()] - ctr; - vector n = e0 ^ e1; - n /= mag(n); + const vector n = normalised(e0 ^ e1); scalar cosAngle = mag(refDir & n); diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C index b5258913447df7894863134abd695c164729c034..68deece2545eee83ff710c61b8ae51616f6b8ff5 100644 --- a/src/dynamicMesh/meshCut/directions/directions.C +++ b/src/dynamicMesh/meshCut/directions/directions.C @@ -315,8 +315,7 @@ Foam::directions::directions vector tan2(globalDict.lookup("tan2")); check2D(correct2DPtr, tan2); - vector normal = tan1 ^ tan2; - normal /= mag(normal); + const vector normal = normalised(tan1 ^ tan2); Info<< "Global Coordinate system:" << endl << " normal : " << normal << endl diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index 39807dc780e23b3f7b334f6775e0078a7e5e71be..2a0b294eba0bb5b39ae38621e5341f04ffc84382 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -1498,8 +1498,7 @@ bool Foam::polyMeshGeometry::checkFaceAngles const face& f = fcs[facei]; - vector faceNormal = faceAreas[facei]; - faceNormal /= mag(faceNormal) + VSMALL; + const vector faceNormal = normalised(faceAreas[facei]); // Get edge from f[0] to f[size-1]; vector ePrev(p[f.first()] - p[f.last()]); @@ -1646,7 +1645,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist // p[f[fpI]], // p[f.nextLabel(fpI)], // fc -// ).normal() +// ).areaNormal() // ); // // scalar magTri = mag(triArea); @@ -1692,20 +1691,30 @@ bool Foam::polyMeshGeometry::checkFaceTwist if (mesh.isInternalFace(facei)) { - nf = cellCentres[nei[facei]] - cellCentres[own[facei]]; - nf /= mag(nf) + VSMALL; + nf = + normalised + ( + cellCentres[nei[facei]] + - cellCentres[own[facei]] + ); } else if (patches[patches.whichPatch(facei)].coupled()) { nf = - neiCc[facei-mesh.nInternalFaces()] - - cellCentres[own[facei]]; - nf /= mag(nf) + VSMALL; + normalised + ( + neiCc[facei-mesh.nInternalFaces()] + - cellCentres[own[facei]] + ); } else { - nf = faceCentres[facei] - cellCentres[own[facei]]; - nf /= mag(nf) + VSMALL; + nf = + normalised + ( + faceCentres[facei] + - cellCentres[own[facei]] + ); } if (nf != vector::zero) @@ -1721,7 +1730,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist p[f[fpI]], p[f.nextLabel(fpI)], fc - ).normal() + ).areaNormal() ); scalar magTri = mag(triArea); @@ -1826,7 +1835,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist p[f[fp]], p[f.nextLabel(fp)], fc - ).normal(); + ).areaNormal(); scalar magTri = mag(prevN); @@ -1853,7 +1862,7 @@ bool Foam::polyMeshGeometry::checkTriangleTwist p[f[fp]], p[f.nextLabel(fp)], fc - ).normal() + ).areaNormal() ); scalar magTri = mag(triN); diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C index 79efcbc1be1223b12a34585ca44020e465bb0355..be4b283afd7a4463020b530ff51759aa95e99572 100644 --- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C +++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C @@ -458,8 +458,11 @@ Foam::label Foam::faceCoupleInfo::mostAlignedCutEdge eVec /= magEVec; - vector eToEndPoint(localPoints[edgeEnd] - localPoints[otherPointi]); - eToEndPoint /= mag(eToEndPoint); + const vector eToEndPoint = + normalised + ( + localPoints[edgeEnd] - localPoints[otherPointi] + ); scalar cosAngle = eVec & eToEndPoint; diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C index 426866008b5094dd8e70de8ae1b7feae598c1669..44c3ca0afce090a18ee6eaa42a242f34622edee4 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C @@ -55,9 +55,8 @@ bool Foam::combineFaces::convexFace const face& f ) { - // Get outwards pointing normal of f. - vector n = f.normal(points); - n /= mag(n); + // Get outwards pointing normal of f, only the sign matters. + const vector areaNorm = f.areaNormal(points); // Get edge from f[0] to f[size-1]; vector ePrev(points[f.first()] - points[f.last()]); @@ -78,7 +77,7 @@ bool Foam::combineFaces::convexFace { vector edgeNormal = ePrev ^ e10; - if ((edgeNormal & n) < 0) + if ((edgeNormal & areaNorm) < 0) { // Concave. Check angle. if ((ePrev & e10) < minConcaveCos) @@ -150,12 +149,10 @@ void Foam::combineFaces::regioniseFaces // - small angle if (p0 != -1 && p0 == p1 && !patches[p0].coupled()) { - vector f0Normal = mesh_.faceAreas()[f0]; - f0Normal /= mag(f0Normal); - vector f1Normal = mesh_.faceAreas()[f1]; - f1Normal /= mag(f1Normal); + vector f0Normal = normalised(mesh_.faceAreas()[f0]); + vector f1Normal = normalised(mesh_.faceAreas()[f1]); - if ((f0Normal&f1Normal) > minCos) + if ((f0Normal & f1Normal) > minCos) { Map<label>::const_iterator f0Fnd = faceRegion.find(f0); diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C index 1df70acf0767d887426dd2eca5be3c492bd33bfb..0608d45f7b3239a5a3aa826ffff5922ff0f869bd 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C @@ -410,13 +410,9 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio if (detJ < 1e-5) { - collapseAxis = f.edges()[f.longestEdge(pts)].vec(pts); - // It is possible that all the points of a face are the same - if (magSqr(collapseAxis) > VSMALL) - { - collapseAxis /= mag(collapseAxis); - } + + collapseAxis = f.edges()[f.longestEdge(pts)].unitVec(pts); // Empirical correlation for high aspect ratio faces @@ -433,9 +429,7 @@ void Foam::edgeCollapser::faceCollapseAxisAndAspectRatio // Cannot necessarily determine linearly independent // eigenvectors, or any at all, use longest edge direction. - collapseAxis = f.edges()[f.longestEdge(pts)].vec(pts); - - collapseAxis /= mag(collapseAxis); + collapseAxis = f.edges()[f.longestEdge(pts)].unitVec(pts); aspectRatio = 1.0; } diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C index e58aa2ae8ea6e69ff5d9104bb3aa5760648efcfb..8d1b96c318499668457b4aa963d92125ced6d1a9 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C @@ -835,11 +835,11 @@ void Foam::hexRef8::checkInternalOrientation face compactFace(identity(newFace.size())); pointField compactPoints(meshMod.points(), newFace); - vector n(compactFace.normal(compactPoints)); + const vector areaNorm(compactFace.areaNormal(compactPoints)); - vector dir(neiPt - ownPt); + const vector dir(neiPt - ownPt); - if ((dir & n) < 0) + if ((dir & areaNorm) < 0) { FatalErrorInFunction << "cell:" << celli << " old face:" << facei @@ -850,9 +850,9 @@ void Foam::hexRef8::checkInternalOrientation << abort(FatalError); } - vector fcToOwn(compactFace.centre(compactPoints) - ownPt); + const vector fcToOwn(compactFace.centre(compactPoints) - ownPt); - scalar s = (fcToOwn&n) / (dir&n); + const scalar s = (fcToOwn & areaNorm) / (dir & areaNorm); if (s < 0.1 || s > 0.9) { @@ -881,11 +881,11 @@ void Foam::hexRef8::checkBoundaryOrientation face compactFace(identity(newFace.size())); pointField compactPoints(meshMod.points(), newFace); - vector n(compactFace.normal(compactPoints)); + const vector areaNorm(compactFace.areaNormal(compactPoints)); - vector dir(boundaryPt - ownPt); + const vector dir(boundaryPt - ownPt); - if ((dir & n) < 0) + if ((dir & areaNorm) < 0) { FatalErrorInFunction << "cell:" << celli << " old face:" << facei @@ -896,9 +896,9 @@ void Foam::hexRef8::checkBoundaryOrientation << abort(FatalError); } - vector fcToOwn(compactFace.centre(compactPoints) - ownPt); + const vector fcToOwn(compactFace.centre(compactPoints) - ownPt); - scalar s = (fcToOwn&dir) / magSqr(dir); + const scalar s = (fcToOwn & dir) / magSqr(dir); if (s < 0.7 || s > 1.3) { diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C index 3e0e48fe3e3ce5bcc12a6a90a8608c8b74511d95..4cd0f0e0ae754f162da251008cfe5d483218a704 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C @@ -213,11 +213,8 @@ Foam::label Foam::removePoints::countPointUsage label vLeft = e0.otherVertex(common); label vRight = e1.otherVertex(common); - vector e0Vec = points[common] - points[vLeft]; - e0Vec /= mag(e0Vec) + VSMALL; - - vector e1Vec = points[vRight] - points[common]; - e1Vec /= mag(e1Vec) + VSMALL; + const vector e0Vec = normalised(points[common] - points[vLeft]); + const vector e1Vec = normalised(points[vRight] - points[common]); if ((e0Vec & e1Vec) > minCos) { diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C index 3dcabff88d44a4eea85502d054c9cfb70e050e15..dea8299a49faca7e2e7c502b4af675f7f22258ad 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C @@ -514,8 +514,8 @@ void Foam::tetDecomposer::setRefinement -1, //edge -1, //face -1, //patchi - zoneI, - zoneFlip + -1, //zone + false ); } // 2b. Within neighbour cell - to cell centre @@ -542,8 +542,8 @@ void Foam::tetDecomposer::setRefinement -1, //edge -1, //face -1, //patchi - zoneI, - zoneFlip + -1, //zone + false ); } } @@ -628,8 +628,8 @@ void Foam::tetDecomposer::setRefinement -1, //edge -1, //face -1, //patchi - zoneI, - zoneFlip + -1, //zone + false ); } @@ -658,8 +658,8 @@ void Foam::tetDecomposer::setRefinement -1, //edge -1, //face -1, //patchi - zoneI, - zoneFlip + -1, //zone + false ); } } @@ -702,14 +702,6 @@ void Foam::tetDecomposer::setRefinement if (decomposeCell[celli]) { - label zoneI = mesh_.faceZones().whichZone(facei); - bool zoneFlip = false; - if (zoneI != -1) - { - const faceZone& fz = mesh_.faceZones()[zoneI]; - zoneFlip = fz.flipMap()[fz.whichFace(facei)]; - } - const face& f = mesh_.faces()[facei]; //const labelList& fEdges = mesh_.faceEdges()[facei]; forAll(f, fp) @@ -843,8 +835,8 @@ void Foam::tetDecomposer::setRefinement -1, //fEdges[fp], //masterEdge facei, //masterFace -1, //patchi - zoneI, - zoneFlip + -1, //zone + false ); } } diff --git a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C index 4484c7f85241d3b795d77ebb57819d91fa1e4a54..2b6354aa8cecffb742b7f079977ba624b09d01e8 100644 --- a/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C +++ b/src/dynamicMesh/slidingInterface/enrichedPatch/enrichedPatchCutFaces.C @@ -150,8 +150,7 @@ void Foam::enrichedPatch::calcCutFaces() const } // Grab face normal - vector normal = curLocalFace.normal(lp); - normal /= mag(normal); + const vector normal = curLocalFace.unitNormal(lp); while (edgeSeeds.size()) { @@ -221,10 +220,9 @@ void Foam::enrichedPatch::calcCutFaces() const // Get the vector along the edge and the right vector vector ahead = curPoint - lp[prevPointLabel]; ahead -= normal*(normal & ahead); - ahead /= mag(ahead); + ahead.normalise(); - vector right = normal ^ ahead; - right /= mag(right); + const vector right = normalised(normal ^ ahead); // Pout<< "normal: " << normal // << " ahead: " << ahead diff --git a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C index ca6d503b0d2170565b916ca1b18d463428fffc39..cbcdb3a68275fccfbd952717eb7f59c161a89aa9 100644 --- a/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C +++ b/src/dynamicMesh/slidingInterface/slidingInterfaceProjectPoints.C @@ -920,15 +920,16 @@ bool Foam::slidingInterface::projectPoints() const // projected-master-to-edge distance is used, to avoid // problems with badly defined master planes. HJ, // 17/Oct/2004 - vector edgeNormalInPlane = - edgeVec - ^ ( - slavePointNormals[curEdge.start()] - + slavePointNormals[curEdge.end()] + const vector edgeNormalInPlane = + normalised + ( + edgeVec + ^ ( + slavePointNormals[curEdge.start()] + + slavePointNormals[curEdge.end()] + ) ); - edgeNormalInPlane /= mag(edgeNormalInPlane); - for (const label cmp : curMasterPoints) { // Skip the current point if the edge start or end has diff --git a/src/fileFormats/stl/STLtriangleI.H b/src/fileFormats/stl/STLtriangleI.H index 541bec22bfe0bc723ae0664fd1ef9ed6d8edb2be..70529371b48d3d35c59acb293b2530f6d16eed91 100644 --- a/src/fileFormats/stl/STLtriangleI.H +++ b/src/fileFormats/stl/STLtriangleI.H @@ -143,9 +143,8 @@ inline void Foam::STLtriangle::write const point& pt2 ) { - // calculate the normal ourselves - vector norm = triPointRef(pt0, pt1, pt2).normal(); - norm /= mag(norm) + VSMALL; + // Calculate the normal ourselves + const vector norm = triPointRef(pt0, pt1, pt2).unitNormal(); write(os, norm, pt0, pt1, pt2); } diff --git a/src/fileFormats/vtk/read/vtkUnstructuredReader.C b/src/fileFormats/vtk/read/vtkUnstructuredReader.C index c20f6276678a8a918d29fc553156d15c038240c0..3094f203e35de4490b82798a3d9446ba73a3aa3a 100644 --- a/src/fileFormats/vtk/read/vtkUnstructuredReader.C +++ b/src/fileFormats/vtk/read/vtkUnstructuredReader.C @@ -979,10 +979,10 @@ void Foam::vtkUnstructuredReader::read(ISstream& inFile) ); const point bottomCc(bottom.centre()); - const vector bottomNormal(bottom.normal()); + const vector bottomNormal(bottom.areaNormal()); const point topCc(top.centre()); - if (((topCc-bottomCc)&bottomNormal) < 0) + if (((topCc - bottomCc) & bottomNormal) < 0) { // Flip top and bottom Swap(shape[0], shape[3]); diff --git a/src/finiteArea/faMesh/faMeshDemandDrivenData.C b/src/finiteArea/faMesh/faMeshDemandDrivenData.C index 82ae59060daa1771d10b47637b7d9044158c3ca1..e4050c7610d95879937bc53601fc79d323210bc5 100644 --- a/src/finiteArea/faMesh/faMeshDemandDrivenData.C +++ b/src/finiteArea/faMesh/faMeshDemandDrivenData.C @@ -425,9 +425,7 @@ void Foam::faMesh::calcFaceAreaNormals() const vectorField& nInternal = faceAreaNormals.ref(); forAll(localFaces, faceI) { - nInternal[faceI] = - localFaces[faceI].normal(localPoints)/ - localFaces[faceI].mag(localPoints); + nInternal[faceI] = localFaces[faceI].unitNormal(localPoints); } forAll(boundary(), patchI) @@ -529,8 +527,7 @@ void Foam::faMesh::calcEdgeAreaNormals() const forAll(edgeAreaNormals.internalField(), edgeI) { - vector e = edges()[edgeI].vec(points()); - e /= mag(e); + const vector e = edges()[edgeI].unitVec(points()); // scalar wStart = // 1.0 - sqr(mag(e^pointNormals[edges()[edgeI].end()])); @@ -580,8 +577,7 @@ void Foam::faMesh::calcEdgeAreaNormals() const pointNormals[patchEdges[edgeI].start()] + pointNormals[patchEdges[edgeI].end()]; - vector e = patchEdges[edgeI].vec(points()); - e /= mag(e); + const vector e = patchEdges[edgeI].unitVec(points()); edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] -= e*(e&edgeAreaNormals.boundaryField()[patchI][edgeI]); @@ -1285,10 +1281,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const // Transform points const vector& origin = points[curPoint]; - vector axis(result[curPoint]/mag(result[curPoint])); + const vector axis = normalised(result[curPoint]); vector dir(allPoints[0] - points[curPoint]); dir -= axis*(axis&dir); - dir /= mag(dir); + dir.normalise(); + coordinateSystem cs("cs", origin, axis, dir); forAll(allPoints, pI) @@ -1567,10 +1564,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const // Transform points const vector& origin = points[curPoint]; - vector axis(result[curPoint]/mag(result[curPoint])); + const vector axis = normalised(result[curPoint]); vector dir(allPoints[0] - points[curPoint]); dir -= axis*(axis&dir); - dir /= mag(dir); + dir.normalise(); + const coordinateSystem cs("cs", origin, axis, dir); scalarField W(allPoints.size(), 1.0); @@ -1738,10 +1736,11 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const // Transform points const vector& origin = points[curPoint]; - vector axis(result[curPoint]/mag(result[curPoint])); + const vector axis = normalised(result[curPoint]); vector dir(allPoints[0] - points[curPoint]); dir -= axis*(axis&dir); - dir /= mag(dir); + dir.normalise(); + coordinateSystem cs("cs", origin, axis, dir); scalarField W(allPoints.size(), 1.0); diff --git a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C index b0c96720a158370161e3aadb142dfe6a6392e979..8c64de4c8a6012a7dffff4e4cf62ffdcab0997c9 100644 --- a/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C +++ b/src/finiteArea/faMesh/faPatches/faPatch/faPatch.C @@ -306,8 +306,8 @@ Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const { - tmp<vectorField> tfN(new vectorField()); - vectorField& fN = tfN.ref(); + auto tfN = tmp<vectorField>::New(); + auto& fN = tfN.ref(); if (ngbPolyPatchIndex() == -1) { @@ -325,8 +325,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const forAll(fN, faceI) { - fN[faceI] = faces[ngbFaces[faceI]].normal(points) - /faces[ngbFaces[faceI]].mag(points); + fN[faceI] = faces[ngbFaces[faceI]].unitNormal(points); } return tfN; diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C index 47292d13774b45940b0c35ff1469493e451ae185..692d6f79f971b0b08fb9602f9b1b49a1522eac56 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C @@ -404,19 +404,16 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const forAll(owner, edgeI) { // Edge normal - area normal - vector edgeNormal = lengths[edgeI]^edges[edgeI].vec(points); - edgeNormal /= mag(edgeNormal); - + vector edgeNormal = + normalised(lengths[edgeI] ^ edges[edgeI].vec(points)); // Unit delta vector vector unitDelta = faceCentres[neighbour[edgeI]] - faceCentres[owner[edgeI]]; - unitDelta -= - edgeNormal*(edgeNormal&unitDelta); - - unitDelta /= mag(unitDelta); + unitDelta -= edgeNormal*(edgeNormal & unitDelta); + unitDelta.normalise(); // Calc PN arc length @@ -447,7 +444,7 @@ void Foam::edgeInterpolation::makeDeltaCoeffs() const // Edge normal - area tangent - edgeNormal = lengths[edgeI]/mag(lengths[edgeI]); + edgeNormal = normalised(lengths[edgeI]); // Delta coefficient as per definition // dc[edgeI] = 1.0/(lPN*(unitDelta & edgeNormal)); @@ -496,7 +493,6 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const const labelUList& neighbour = mesh().neighbour(); const edgeVectorField& lengths = mesh().Le(); - const edgeScalarField& magLengths = mesh().magLe(); const edgeList& edges = mesh().edges(); const pointField& points = mesh().points(); @@ -508,8 +504,8 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const forAll(owner, edgeI) { // Edge normal - area normal - vector edgeNormal = lengths[edgeI] ^ edges[edgeI].vec(points); - edgeNormal /= mag(edgeNormal); + vector edgeNormal = + normalised(lengths[edgeI] ^ edges[edgeI].vec(points)); // Unit delta vector vector unitDelta = @@ -517,11 +513,10 @@ void Foam::edgeInterpolation::makeCorrectionVectors() const - faceCentres[owner[edgeI]]; unitDelta -= edgeNormal*(edgeNormal & unitDelta); - - unitDelta /= mag(unitDelta); + unitDelta.normalise(); // Edge normal - area tangent - edgeNormal = lengths[edgeI]/magLengths[edgeI]; + edgeNormal = normalised(lengths[edgeI]); // Delta coeffs deltaCoeffs[edgeI] = 1.0/(unitDelta & edgeNormal); diff --git a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C index 9d7d2c2516ef1137cb90e1c9c9e69aa292bf9802..d8ef7244b557516dafbaabcd6b4db4ea61674e93 100644 --- a/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C +++ b/src/finiteVolume/cfdTools/general/SRF/SRFModel/SRFModel/SRFModel.C @@ -59,14 +59,11 @@ Foam::SRF::SRFModel::SRFModel ), Urel_(Urel), mesh_(Urel_.mesh()), - origin_("origin", dimLength, lookup("origin")), - axis_(lookup("axis")), + origin_("origin", dimLength, get<vector>("origin")), + axis_(normalised(get<vector>("axis"))), SRFModelCoeffs_(optionalSubDict(type + "Coeffs")), omega_(dimensionedVector("omega", dimless/dimTime, Zero)) -{ - // Normalise the axis - axis_ /= mag(axis_); -} +{} // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // @@ -82,21 +79,19 @@ bool Foam::SRF::SRFModel::read() if (regIOobject::read()) { // Re-read origin - lookup("origin") >> origin_; + readEntry("origin", origin_); // Re-read axis - lookup("axis") >> axis_; - axis_ /= mag(axis_); + readEntry("axis", axis_); + axis_.normalise(); // Re-read sub-model coeffs SRFModelCoeffs_ = optionalSubDict(type() + "Coeffs"); return true; } - else - { - return false; - } + + return false; } diff --git a/src/finiteVolume/cfdTools/general/coupling/externalFileCoupler.H b/src/finiteVolume/cfdTools/general/coupling/externalFileCoupler.H index ec961499fbdd393e71485cea739fad9d64c23171..b0aab9c3474502dd714f3e7080b0767a05e67e91 100644 --- a/src/finiteVolume/cfdTools/general/coupling/externalFileCoupler.H +++ b/src/finiteVolume/cfdTools/general/coupling/externalFileCoupler.H @@ -243,13 +243,13 @@ public: // content which corresponds to particular return values: // - \c status=done // \return Foam::Time::saEndTime + // - \c action=noWriteNow + // \return Foam::Time::saNoWriteNow // - \c action=writeNow // \return Foam::Time::saWriteNow // - \c action=nextWrite // \return Foam::Time::saNextWrite - // - \c action=noNextWrite - // \return Foam::Time::saNoNextWrite - // - Anything else (empty file, no action= or status=, etc) + // - Anything else (empty file, no action= or status=, etc) // \return Foam::Time::saUnknown enum Time::stopAtControls waitForSlave() const; diff --git a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C index 3dff767145ec5af98ee80e81d5a938006cd9f407..74f0903c712d06d2b04a2375616c409ce56014c4 100644 --- a/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C +++ b/src/finiteVolume/fields/fvPatchFields/derived/mappedField/mappedFieldFvPatchField.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -144,6 +144,7 @@ void Foam::mappedFieldFvPatchField<Type>::autoMap const fvPatchFieldMapper& m ) { + fixedValueFvPatchField<Type>::autoMap(m); mappedPatchBase::clearOut(); } @@ -155,6 +156,7 @@ void Foam::mappedFieldFvPatchField<Type>::rmap const labelList& addr ) { + fixedValueFvPatchField<Type>::rmap(ptf, addr); mappedPatchBase::clearOut(); } diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C index 0b3179edd9b94e094d87591cd9cc678c46aadb03..7bd47b8991960bdb7d748b77aea8e0acc4068e27 100644 --- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C +++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutCell/isoCutCell.C @@ -267,12 +267,11 @@ void Foam::isoCutCell::calcIsoFacePointsFromEdges() // Defining local coordinates with zhat along isoface normal and xhat from // isoface centre to first point in isoFaceEdges_ - const vector zhat = isoFaceArea_/mag(isoFaceArea_); + const vector zhat = normalised(isoFaceArea_); vector xhat = isoFaceEdges_[0][0] - isoFaceCentre_; xhat = (xhat - (xhat & zhat)*zhat); - xhat /= mag(xhat); - vector yhat = zhat^xhat; - yhat /= mag(yhat); + xhat.normalise(); + vector yhat = normalised(zhat ^ xhat); DebugPout << "Calculated local coordinates" << endl; diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C index e0e335b55d1b672a71e2c88d607b43abed8df312..62a6ed602a60478b2c39ab310cb340955bbba915 100644 --- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C +++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoCutFace/isoCutFace.C @@ -735,12 +735,12 @@ void Foam::isoCutFace::quadAreaCoeffs if (Bx > 10*SMALL) { // If |AB| > 0 ABCD we use AB to define xhat - xhat = (B - A)/mag(B - A); + xhat = normalised(B - A); } else if (mag(C - D) > 10*SMALL) { // If |AB| ~ 0 ABCD is a triangle ACD and we use CD for xhat - xhat = (C - D)/mag(C - D); + xhat = normalised(C - D); } else { diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H index 07baf8e69e83da79c8298e106b5f3aa80a6d8e5c..11b6377918a63487eac2247639a78292e8962488 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTet.H @@ -81,10 +81,11 @@ bool interpolationCellPointFace<Type>::findTet referencePoint = tetPoints[p1]; faceNormal = - (tetPoints[p3] - tetPoints[p1]) - ^ (tetPoints[p2] - tetPoints[p1]); - - faceNormal /= mag(faceNormal); + normalised + ( + (tetPoints[p3] - tetPoints[p1]) + ^ (tetPoints[p2] - tetPoints[p1]) + ); // correct normal to point into the tet vector v0 = tetPoints[n] - referencePoint; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H index e6b2fea84b9fc06c5a6feb88c042b0a123009525..2c721d1556c804eb6945f55db2240135043b713a 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/findCellPointFaceTriangle.H @@ -57,8 +57,7 @@ bool interpolationCellPointFace<Type>::findTriangle // calculate edge normal (pointing inwards) for (label i=0; i<3; i++) { - normal[i] = triangleFaceNormal ^ edge[i]; - normal[i] /= mag(normal[i]) + VSMALL; + normal[i] = normalised(triangleFaceNormal ^ edge[i]); } // check if position is inside triangle diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C index d7b130a1ca64d37e813054bbfaf709c8497fad52..9b218164303d475c061cebb37febb6e26ccd68a3 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointFace/interpolationCellPointFace.C @@ -88,12 +88,9 @@ Type Foam::interpolationCellPointFace<Type>::interpolate label closestFace = -1; scalar minDistance = GREAT; - forAll(cellFaces, facei) + for (const label nFace : cellFaces) { - label nFace = cellFaces[facei]; - - vector normal = this->pMeshFaceAreas_[nFace]; - normal /= mag(normal); + const vector normal = normalised(this->pMeshFaceAreas_[nFace]); const vector& faceCentreTmp = this->pMeshFaceCentres_[nFace]; diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C index 064bf94aea345e9fe3a4fa89faac176d326a89ef..f9f98c5ba229eb1bc8b45a3bddfdbf1f40c35fb9 100644 --- a/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C +++ b/src/finiteVolume/interpolation/interpolation/interpolationPointMVC/pointMVCWeight.C @@ -156,10 +156,8 @@ void Foam::pointMVCWeight::calcWeights label jPlus1 = f.fcIndex(j); //Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl; - vector n0 = u[j]^v; - n0 /= mag(n0); - vector n1 = u[jPlus1]^v; - n1 /= mag(n1); + const vector n0 = normalised(u[j] ^ v); + const vector n1 = normalised(u[jPlus1] ^ v); scalar l = min(mag(n0 - n1), 2.0); //Pout<< " l:" << l << endl; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C index 6fe277786b6ad4b26b7731190980908c38c775c8..d157e4478c8c8838b12f2d6a21f80367f761b6cf 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C @@ -76,8 +76,7 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::findFaceDirs { const fvMesh& mesh = this->mesh(); - idir = mesh.faceAreas()[facei]; - idir /= mag(idir); + idir = normalised(mesh.faceAreas()[facei]); #ifndef SPHERICAL_GEOMETRY if (mesh.nGeometricD() <= 2) // find the normal direction diff --git a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C index 41b6eb06bfc5a6a2df4870c81e14d07647b18bd2..ed0c304a218900a07c1d4206b0b6aa2abee618ed 100644 --- a/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C +++ b/src/functionObjects/field/wallBoundedStreamLine/wallBoundedParticleTemplates.C @@ -127,7 +127,7 @@ Foam::scalar Foam::wallBoundedParticle::trackToEdge const bool posVol = (nbrTi.tet(mesh()).mag() > 0); const vector path(endPosition - localPosition_); - if (posVol == ((nbrTi.faceTri(mesh()).normal() & path) < 0)) + if (posVol == ((nbrTi.faceTri(mesh()).areaNormal() & path) < 0)) { // Change into nbrCell. No need to change tetFace, tetPt. //Pout<< " crossed from cell:" << celli_ diff --git a/src/functionObjects/utilities/Make/options b/src/functionObjects/utilities/Make/options index a9844b6d0e5f6524940b3d9a2b76073f6b8f6d79..9cccdbc1a446f565c9c6c70be62b2ef71afe289b 100644 --- a/src/functionObjects/utilities/Make/options +++ b/src/functionObjects/utilities/Make/options @@ -4,6 +4,7 @@ EXE_INC = \ -I$(LIB_SRC)/fileFormats/lnInclude \ -I$(LIB_SRC)/conversion/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/ODE/lnInclude \ -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ -I$(LIB_SRC)/transportModels/compressible/lnInclude @@ -13,6 +14,7 @@ LIB_LIBS = \ -lfiniteArea \ -lconversion \ -lsampling \ + -ldynamicMesh \ -lfluidThermophysicalModels \ -lcompressibleTransportModels \ -lODE diff --git a/src/functionObjects/utilities/ensightWrite/ensightWrite.C b/src/functionObjects/utilities/ensightWrite/ensightWrite.C index c228b5964a49144c9b95336b119f8088286d0d1a..518c74c17ecd9f0d508bd7457f3c7352ad5306a8 100644 --- a/src/functionObjects/utilities/ensightWrite/ensightWrite.C +++ b/src/functionObjects/utilities/ensightWrite/ensightWrite.C @@ -47,6 +47,39 @@ namespace functionObjects // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +Foam::bitSet Foam::functionObjects::ensightWrite::cellSelection() const +{ + bitSet cellsToSelect; + + // Could also deal with cellZones here, as required + + if (bounds_.empty()) + { + return cellsToSelect; + } + + + const auto& cellCentres = static_cast<const fvMesh&>(mesh_).C(); + + const label len = mesh_.nCells(); + + cellsToSelect.resize(len); + + for (label celli=0; celli < len; ++celli) + { + const point& cc = cellCentres[celli]; + + if (bounds_.contains(cc)) + { + cellsToSelect.set(celli); + } + } + + return cellsToSelect; +} + + + int Foam::functionObjects::ensightWrite::process(const word& fieldName) { int state = 0; @@ -83,7 +116,8 @@ Foam::functionObjects::ensightWrite::ensightWrite caseOpts_(writeOpts_.format()), selectFields_(), dirName_("ensightWrite"), - consecutive_(false) + consecutive_(false), + bounds_() { if (postProcess) { @@ -104,6 +138,12 @@ bool Foam::functionObjects::ensightWrite::read(const dictionary& dict) { fvMeshFunctionObject::read(dict); + // Ensure consistency + + meshSubset_.clear(); + ensMesh_.clear(); + + // // writer options // @@ -123,6 +163,10 @@ bool Foam::functionObjects::ensightWrite::read(const dictionary& dict) } + bounds_.clear(); + dict.readIfPresent("bounds", bounds_); + + // // case options // @@ -197,7 +241,25 @@ bool Foam::functionObjects::ensightWrite::write() if (!ensMesh_.valid()) { writeGeom = true; - ensMesh_.reset(new ensightMesh(mesh_, writeOpts_)); + + bitSet selection = this->cellSelection(); + + if (returnReduce(!selection.empty(), orOp<bool>())) + { + meshSubset_.clear(); + meshSubset_.reset(new fvMeshSubset(mesh_, selection)); + ensMesh_.reset + ( + new ensightMesh(meshSubset_->subMesh(), writeOpts_) + ); + } + else + { + ensMesh_.reset + ( + new ensightMesh(mesh_, writeOpts_) + ); + } } if (ensMesh_().needsUpdate()) { @@ -219,12 +281,12 @@ bool Foam::functionObjects::ensightWrite::write() DynamicList<word> missing(selectFields_.size()); DynamicList<word> ignored(selectFields_.size()); - // check exact matches first + // Check exact matches first for (const wordRe& select : selectFields_) { if (!select.isPattern()) { - const word& fieldName = static_cast<const word&>(select); + const word& fieldName = select; if (!candidates.erase(fieldName)) { @@ -271,7 +333,14 @@ void Foam::functionObjects::ensightWrite::updateMesh(const mapPolyMesh& mpm) { // fvMeshFunctionObject::updateMesh(mpm); - if (ensMesh_.valid()) + // This is heavy-handed, but with a bounding-box limited sub-mesh, + // we don't readily know if the updates affect the subsetted mesh. + if (!bounds_.empty()) + { + ensMesh_.clear(); + meshSubset_.clear(); + } + else if (ensMesh_.valid()) { ensMesh_->expire(); } @@ -282,7 +351,14 @@ void Foam::functionObjects::ensightWrite::movePoints(const polyMesh& mpm) { // fvMeshFunctionObject::updateMesh(mpm); - if (ensMesh_.valid()) + // This is heavy-handed, but with a bounding-box limited sub-mesh, + // we don't readily know if the updates affect the subsetted mesh. + if (!bounds_.empty()) + { + ensMesh_.clear(); + meshSubset_.clear(); + } + else if (ensMesh_.valid()) { ensMesh_->expire(); } diff --git a/src/functionObjects/utilities/ensightWrite/ensightWrite.H b/src/functionObjects/utilities/ensightWrite/ensightWrite.H index 25954ea89c050565f3d8f9e51282109f2ca0d40b..d549623a20721c46f579686c38dee8105f8fdac0 100644 --- a/src/functionObjects/utilities/ensightWrite/ensightWrite.H +++ b/src/functionObjects/utilities/ensightWrite/ensightWrite.H @@ -67,6 +67,7 @@ Description faceZones | Select faceZones to write | no | consecutive | Consecutive output numbering | no | false nodeValues | Write values at nodes | no | false + bounds | Limit with a bounding box | no | \endtable Note that if the \c patches entry is an empty list, this will select all @@ -93,6 +94,7 @@ SourceFiles #include "interpolation.H" #include "volFields.H" +#include "fvMeshSubset.H" #include "surfaceFields.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -129,6 +131,12 @@ class ensightWrite //- Consecutive output numbering bool consecutive_; + //- Restrict to bounding box + boundBox bounds_; + + //- Mesh subset handler + autoPtr<fvMeshSubset> meshSubset_; + //- Ensight case handler autoPtr<ensightCase> ensCase_; @@ -150,6 +158,9 @@ class ensightWrite return *ensMesh_; } + //- Define cell selection from bounding box. + // An empty set if no bounding box is specified. + bitSet cellSelection() const; //- Apply for the volume field type template<class Type> diff --git a/src/functionObjects/utilities/ensightWrite/ensightWriteTemplates.C b/src/functionObjects/utilities/ensightWrite/ensightWriteTemplates.C index e68858f3fd5d80481c52f646cfa4a99a1e8b9fb1..a1da699b44c75e329d8eeace09918a3a761a1912 100644 --- a/src/functionObjects/utilities/ensightWrite/ensightWriteTemplates.C +++ b/src/functionObjects/utilities/ensightWrite/ensightWriteTemplates.C @@ -38,20 +38,39 @@ int Foam::functionObjects::ensightWrite::writeVolField // State: return 0 (not-processed), -1 (skip), +1 ok typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType; + const VolFieldType* fldPtr; + // Already done, or not available - if (state || !foundObject<VolFieldType>(inputName)) + if (state || !(fldPtr = lookupObjectPtr<VolFieldType>(inputName))) { return state; } autoPtr<ensightFile> os = ensCase().newData<Type>(inputName); - ensightOutput::writeField<Type> - ( - lookupObject<VolFieldType>(inputName), - ensMesh(), - os, - caseOpts_.nodeValues() - ); + + + if (meshSubset_.valid()) + { + tmp<VolFieldType> tfield = meshSubset_->interpolate(*fldPtr); + + ensightOutput::writeField<Type> + ( + tfield(), + ensMesh(), + os, + caseOpts_.nodeValues() + ); + } + else + { + ensightOutput::writeField<Type> + ( + *fldPtr, + ensMesh(), + os, + caseOpts_.nodeValues() + ); + } Log << " " << inputName; diff --git a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C index cdc5ec680c77ef187546eff4fe6949912a3b9e5e..843dfb1a39c523c30d0a4a0412e8906201e6dd29 100644 --- a/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C +++ b/src/lagrangian/DSMC/parcels/Templates/DSMCParcel/DSMCParcel.C @@ -116,8 +116,7 @@ void Foam::DSMCParcel<ParcelType>::hitWallPatch scalar m = constProps.mass(); - vector nw = wpp.faceAreas()[wppLocalFace]; - nw /= mag(nw); + const vector nw = normalised(wpp.faceAreas()[wppLocalFace]); scalar U_dot_nw = U_ & nw; diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C index 154faee32bef5adb2f706b8d1d616322d91a68dc..4d6b0a184c9b2e14cdc13d6f6b79fba3155bd1db 100644 --- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C +++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C @@ -265,18 +265,15 @@ void Foam::FreeStream<CloudType>::inflow() // Normal unit vector *negative* so normal is pointing into the // domain - vector n = patch.faceAreas()[pFI]; - n /= -mag(n); + const vector n = -normalised(patch.faceAreas()[pFI]); // Wall tangential unit vector. Use the direction between the // face centre and the first vertex in the list - vector t1 = fC - (mesh.points()[f[0]]); - t1 /= mag(t1); + const vector t1 = normalised(fC - (mesh.points()[f[0]])); // Other tangential unit vector. Rescaling in case face is not // flat and n and t1 aren't perfectly orthogonal - vector t2 = n^t1; - t2 /= mag(t2); + const vector t2 = normalised(n ^ t1); scalar faceTemperature = boundaryT[patchi][pFI]; diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C index dd145c62e1c2c966ee43f07801aaa9803ed9fdd0..ca2bc100418c05ef57c49fdfd0252b6c65b6ed3d 100644 --- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C +++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C @@ -68,8 +68,7 @@ void Foam::MaxwellianThermal<CloudType>::correct label wppLocalFace = wpp.whichFace(p.face()); - vector nw = p.normal(); - nw /= mag(nw); + const vector nw = p.normal(); // Normal velocity magnitude scalar U_dot_nw = U & nw; diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C index 21c2cd73a8b56b775b3807fe33fa5113c3182c56..4fe4f0649d798a393d9de5c0d84e9b9ff2da5786 100644 --- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C +++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C @@ -66,8 +66,7 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct label wppLocalFace = wpp.whichFace(p.face()); - vector nw = p.normal(); - nw /= mag(nw); + const vector nw = p.normal(); // Normal velocity magnitude scalar U_dot_nw = U & nw; diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C index a3d182d07b91067f9adfd6f4d2f9fe0c8864929f..2d9a00df536d4c5b01e3cd08d4f66bdc15d061a5 100644 --- a/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C +++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C @@ -57,8 +57,7 @@ void Foam::SpecularReflection<CloudType>::correct { vector& U = p.U(); - vector nw = p.normal(); - nw /= mag(nw); + const vector nw = p.normal(); scalar U_dot_nw = U & nw; diff --git a/src/lagrangian/basic/particle/particle.H b/src/lagrangian/basic/particle/particle.H index 4d6928e2b1842dbd8290bbb42c07b28da96e1613..f7e9f445327e2f0e8a73f602a429c099f7acc269 100644 --- a/src/lagrangian/basic/particle/particle.H +++ b/src/lagrangian/basic/particle/particle.H @@ -522,8 +522,7 @@ public: //- Return the current tet transformation tensor inline barycentricTensor currentTetTransform() const; - //- Return the normal of the tri on tetFacei_ for the - // current tet. + //- The (unit) normal of the tri on tetFacei_ for the current tet. inline vector normal() const; //- Is the particle on a face? diff --git a/src/lagrangian/basic/particle/particleI.H b/src/lagrangian/basic/particle/particleI.H index 771eeec0575bcd9da12ef9d1e213ec1ce92890e3..3698b17f2232b08aae8f030ec1b4e76afa388ca7 100644 --- a/src/lagrangian/basic/particle/particleI.H +++ b/src/lagrangian/basic/particle/particleI.H @@ -276,7 +276,7 @@ inline Foam::barycentricTensor Foam::particle::currentTetTransform() const inline Foam::vector Foam::particle::normal() const { - return currentTetIndices().faceTri(mesh_).normal(); + return currentTetIndices().faceTri(mesh_).unitNormal(); } @@ -324,8 +324,7 @@ void Foam::particle::patchData(vector& n, vector& U) const Pair<vector> centre, base, vertex1, vertex2; movingTetGeometry(1, centre, base, vertex1, vertex2); - n = triPointRef(base[0], vertex1[0], vertex2[0]).normal(); - n /= mag(n); + n = triPointRef(base[0], vertex1[0], vertex2[0]).unitNormal(); // Interpolate the motion of the three face vertices to the current // coordinates @@ -343,8 +342,7 @@ void Foam::particle::patchData(vector& n, vector& U) const vector centre, base, vertex1, vertex2; stationaryTetGeometry(centre, base, vertex1, vertex2); - n = triPointRef(base, vertex1, vertex2).normal(); - n /= mag(n); + n = triPointRef(base, vertex1, vertex2).unitNormal(); U = Zero; } diff --git a/src/lagrangian/basic/particle/particleTemplates.C b/src/lagrangian/basic/particle/particleTemplates.C index a43c7763972a286767974e6d8d8c44a02816d5a2..fffeacd3f0bac13864d530481a75aa62db4d071b 100644 --- a/src/lagrangian/basic/particle/particleTemplates.C +++ b/src/lagrangian/basic/particle/particleTemplates.C @@ -250,8 +250,7 @@ void Foam::particle::hitSymmetryPlanePatch template<class TrackCloudType> void Foam::particle::hitSymmetryPatch(TrackCloudType&, trackingData&) { - vector nf = normal(); - nf /= mag(nf); + const vector nf = normal(); transformProperties(I - 2.0*nf*nf); } diff --git a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C index 1a8255ae7d615056a3be9057c3adac92d0b371c6..d8bca67835d45cc6f662099cbf22e5880485b77e 100644 --- a/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C +++ b/src/lagrangian/intermediate/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C @@ -167,7 +167,7 @@ void Foam::ParticleCollector<CloudType>::initConcentricCircles() { refDir = this->coeffDict().lookup("refDir"); refDir -= normal_[0]*(normal_[0] & refDir); - refDir /= mag(refDir); + refDir.normalise(); } else { @@ -315,13 +315,14 @@ void Foam::ParticleCollector<CloudType>::collectParcelPolygon // the face's decomposed triangles does not work due to ambiguity along // the diagonals. const face& f = faces_[facei]; - const vector n = f.normal(points_); + const vector areaNorm = f.areaNormal(points_); bool inside = true; - for (label i = 0; i < f.size(); ++ i) + for (label i = 0; i < f.size(); ++i) { const label j = f.fcIndex(i); const triPointRef t(pIntersect, points_[f[i]], points_[f[j]]); - if ((n & t.normal()) < 0) + + if ((areaNorm & t.areaNormal()) < 0) { inside = false; break; @@ -580,8 +581,7 @@ Foam::ParticleCollector<CloudType>::ParticleCollector forAll(polygons, polyI) { polygons[polyI] = polygonAndNormal[polyI].first(); - normal_[polyI] = polygonAndNormal[polyI].second(); - normal_[polyI] /= mag(normal_[polyI]) + ROOTVSMALL; + normal_[polyI] = normalised(polygonAndNormal[polyI].second()); } initPolygons(polygons); diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C index 7d68828701da8dd440f70f354ac2d6650df8fd50..96407263374762f04cb4ed0c794aebc40055adf6 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairCollision.C @@ -238,15 +238,14 @@ void Foam::PairCollision<CloudType>::wallInteraction() if (nearest.distance() < r) { - vector normal = mesh.faceAreas()[realFacei]; - - normal /= mag(normal); + const vector normal = + normalised(mesh.faceAreas()[realFacei]); const vector& nearPt = nearest.rawPoint(); - vector pW = nearPt - pos; + const vector pW = normalised(nearPt - pos); - scalar normalAlignment = normal & pW/(mag(pW) + SMALL); + const scalar normalAlignment = normal & pW; // Find the patchIndex and wallData for WallSiteData object label patchi = patchID[realFacei - mesh.nInternalFaces()]; @@ -315,15 +314,13 @@ void Foam::PairCollision<CloudType>::wallInteraction() if (nearest.distance() < r) { - vector normal = rwf.normal(pts); - - normal /= mag(normal); + const vector normal = rwf.unitNormal(pts); const vector& nearPt = nearest.rawPoint(); - vector pW = nearPt - pos; + const vector pW = normalised(nearPt - pos); - scalar normalAlignment = normal & pW/mag(pW); + const scalar normalAlignment = normal & pW; // Find the patchIndex and wallData for WallSiteData object diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C index db0296b9147308243770bc074301dbdf3185407d..1e203baca64f63ee2e0a37d4e1fe7d89de83d426 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeInjection/ConeInjection.C @@ -104,8 +104,7 @@ Foam::ConeInjection<CloudType>::ConeInjection forAll(positionAxis_, i) { vector& axis = positionAxis_[i].second(); - - axis /= mag(axis); + axis.normalise(); vector tangent = Zero; scalar magTangent = 0.0; @@ -277,7 +276,7 @@ void Foam::ConeInjection<CloudType>::setProperties vector normal = alpha*(tanVec1_[i]*cos(beta) + tanVec2_[i]*sin(beta)); vector dirVec = dcorr*positionAxis_[i].second(); dirVec += normal; - dirVec /= mag(dirVec); + dirVec.normalise(); parcel.U() = Umag_.value(t)*dirVec; diff --git a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C index 065bf8124840aa067ce9a6763e894a7e1f26a880..eba6578fcbe66d053a928ec14441864c0bdcab2e 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C +++ b/src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/ConeNozzleInjection/ConeNozzleInjection.C @@ -206,7 +206,7 @@ Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection Random& rndGen = this->owner().rndGen(); // Normalise direction vector - direction_ /= mag(direction_); + direction_.normalise(); // Determine direction vectors tangential to direction vector tangent = Zero; @@ -436,7 +436,7 @@ void Foam::ConeNozzleInjection<CloudType>::setProperties vector normal = alpha*normal_; vector dirVec = dcorr*direction_; dirVec += normal; - dirVec /= mag(dirVec); + dirVec.normalise(); switch (flowType_) { diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C index e819c7fc9e9f4191d0bc81acae178eb5f7f809a7..7adb7b5806bd20a6e5caa026a6c6f13b24f71aa0 100644 --- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C +++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C @@ -344,7 +344,7 @@ Foam::vector Foam::PackingModels::Implicit<CloudType>::velocityCorrection const vector U = uCorrect_()[celli]; // face geometry - vector nHat = mesh.faces()[facei].normal(mesh.points()); + vector nHat = mesh.faces()[facei].areaNormal(mesh.points()); const scalar nMag = mag(nHat); nHat /= nMag; diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C index d93a6985c92dba86afc63c268a37f3ebced629f6..52e3a0a8b9ef2b22a8f50813932935b1da1d5d7d 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C +++ b/src/lagrangian/molecularDynamics/molecule/molecule/molecule.C @@ -249,10 +249,9 @@ void Foam::molecule::hitProcessorPatch(moleculeCloud&, trackingData& td) void Foam::molecule::hitWallPatch(moleculeCloud&, trackingData&) { - vector nw = normal(); - nw /= mag(nw); + const vector nw = normal(); - scalar vn = v_ & nw; + const scalar vn = v_ & nw; // Specular reflection if (vn > 0) diff --git a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H index 947998d4c07bd29d9f7951728f806f178804e709..80b5045dfac3fb8c788980bbe5e46a8aa7e21dd4 100644 --- a/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H +++ b/src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H @@ -90,9 +90,11 @@ inline Foam::molecule::constantProperties::constantProperties Info<< nl << "Linear molecule." << endl; - vector dir = siteReferencePositions_[1] - siteReferencePositions_[0]; - - dir /= mag(dir); + const vector dir = + normalised + ( + siteReferencePositions_[1] - siteReferencePositions_[0] + ); tensor Q = rotationTensor(dir, vector(1,0,0)); @@ -334,9 +336,11 @@ inline bool Foam::molecule::constantProperties::linearMoleculeTest() const return true; } - vector refDir = siteReferencePositions_[1] - siteReferencePositions_[0]; - - refDir /= mag(refDir); + const vector refDir = + normalised + ( + siteReferencePositions_[1] - siteReferencePositions_[0] + ); for ( @@ -345,9 +349,11 @@ inline bool Foam::molecule::constantProperties::linearMoleculeTest() const i++ ) { - vector dir = siteReferencePositions_[i] - siteReferencePositions_[i-1]; - - dir /= mag(dir); + const vector dir = + normalised + ( + siteReferencePositions_[i] - siteReferencePositions_[i-1] + ); if (mag(refDir & dir) < 1 - SMALL) { diff --git a/src/lagrangian/solidParticle/solidParticle.C b/src/lagrangian/solidParticle/solidParticle.C index 568c1b4c7d66721048ca45761f81cc39a162a6e3..dfc1426e760655ee56d70359855b3cade3225650 100644 --- a/src/lagrangian/solidParticle/solidParticle.C +++ b/src/lagrangian/solidParticle/solidParticle.C @@ -104,8 +104,7 @@ void Foam::solidParticle::hitProcessorPatch void Foam::solidParticle::hitWallPatch(solidParticleCloud& cloud, trackingData&) { - vector nw = normal(); - nw /= mag(nw); + const vector nw = normal(); scalar Un = U_ & nw; vector Ut = U_ - Un*nw; diff --git a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C index bc228d90fbdfc5616382066a61f6f0c41eff74e7..3a0be877113a67459c55268794b919e9b3a74195 100644 --- a/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C +++ b/src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C @@ -43,7 +43,7 @@ Foam::LISAAtomization<CloudType>::LISAAtomization SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod")) { // Note: Would be good if this could be picked up from the injector - injectorDirection_ /= mag(injectorDirection_); + injectorDirection_.normalise(); if (SMDCalcMethod_ == "method1") { diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C index 7593b3763e16bbe71e1a4663924269d1cd55fc7b..175ecf96888184028224cedf90e0a71afef86f04 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C @@ -59,7 +59,8 @@ void Foam::blockDescriptor::check(const Istream& is) forAll(faces, i) { point faceCentre(faces[i].centre(vertices_)); - vector faceNormal(faces[i].normal(vertices_)); + vector faceNormal(faces[i].areaNormal(vertices_)); + if (mag(faceNormal) > SMALL) { if (((faceCentre - blockCentre) & faceNormal) > 0) diff --git a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C index fc0c08b8623443bfcd277132f852d8755a958a00..25f1ddfdf9f2019e0b6f7ee06f0836232d8fa1a5 100644 --- a/src/mesh/blockMesh/blockMesh/blockMeshCheck.C +++ b/src/mesh/blockMesh/blockMesh/blockMeshCheck.C @@ -187,9 +187,9 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const if ( ( - patchFace.normal(points) - & faces[cellFaces[cellFacei]].normal(points) - ) < 0.0 + patchFace.areaNormal(points) + & faces[cellFaces[cellFacei]].areaNormal(points) + ) < 0 ) { Info<< tab << tab diff --git a/src/mesh/extrudeModel/linearDirection/linearDirection.C b/src/mesh/extrudeModel/linearDirection/linearDirection.C index 36a0094d2eddc13365f4b157246b65ae68d580a7..fba682302216616616aba6d1f845b513ce9e0905 100644 --- a/src/mesh/extrudeModel/linearDirection/linearDirection.C +++ b/src/mesh/extrudeModel/linearDirection/linearDirection.C @@ -45,11 +45,9 @@ addToRunTimeSelectionTable(extrudeModel, linearDirection, dictionary); linearDirection::linearDirection(const dictionary& dict) : extrudeModel(typeName, dict), - direction_(coeffDict_.lookup("direction")), - thickness_(readScalar(coeffDict_.lookup("thickness"))) + direction_(coeffDict_.get<vector>("direction").normalise()), + thickness_(coeffDict_.get<scalar>("thickness")) { - direction_ /= mag(direction_); - if (thickness_ <= 0) { FatalErrorInFunction diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C index 9253faac3bca0014a7d472a2cac44e921ee55006..42560a26393da1466f18a96dc105c048a9664785 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/fieldSmoother/fieldSmoother.C @@ -129,8 +129,7 @@ void Foam::fieldSmoother::smoothNormals { //full smoothing neighbours + point value average[pointI] = 0.5*(normals[pointI]+average[pointI]); - normals[pointI] = average[pointI]; - normals[pointI] /= mag(normals[pointI]) + VSMALL; + normals[pointI] = normalised(average[pointI]); } } } @@ -196,8 +195,7 @@ void Foam::fieldSmoother::smoothPatchNormals { // full smoothing neighbours + point value average[pointI] = 0.5*(normals[pointI]+average[pointI]); - normals[pointI] = average[pointI]; - normals[pointI] /= mag(normals[pointI]) + VSMALL; + normals[pointI] = normalised(average[pointI]); } } } diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C index 094d83af79f4bb05c0b7ffd21ffcb232888f14c9..6bb99d604db052694bf7ee2c7c31808597db1586 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C @@ -58,33 +58,42 @@ bool Foam::medialAxisMeshMover::isMaxEdge ( const List<PointData<vector>>& pointWallDist, const label edgeI, - const scalar minCos + const scalar minCos, + const bool disableWallEdges ) const { const pointField& points = mesh().points(); - - // Do not mark edges with one side on moving wall. - const edge& e = mesh().edges()[edgeI]; - vector v0(points[e[0]] - pointWallDist[e[0]].origin()); - scalar magV0(mag(v0)); - - if (magV0 < SMALL) + if (disableWallEdges) { - return false; - } - - vector v1(points[e[1]] - pointWallDist[e[1]].origin()); - scalar magV1(mag(v1)); + // 1. Do not mark edges with one side on moving wall. + vector v0(points[e[0]] - pointWallDist[e[0]].origin()); + scalar magV0(mag(v0)); + if (magV0 < SMALL) + { + return false; + } - if (magV1 < SMALL) - { - return false; + vector v1(points[e[1]] - pointWallDist[e[1]].origin()); + scalar magV1(mag(v1)); + if (magV1 < SMALL) + { + return false; + } } + //// 2. Do not mark edges with both sides on a moving wall. + //vector v0(points[e[0]] - pointWallDist[e[0]].origin()); + //scalar magV0(mag(v0)); + //vector v1(points[e[1]] - pointWallDist[e[1]].origin()); + //scalar magV1(mag(v1)); + //if (magV0 < SMALL && magV1 < SMALL) + //{ + // return false; + //} - //- Detect based on vector to nearest point differing for both endpoints + //// 3. Detect based on vector to nearest point differing for both endpoints //v0 /= magV0; //v1 /= magV1; // @@ -170,6 +179,13 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) mesh().globalData().nTotalPoints() ); + const bool disableWallEdges = coeffDict.lookupOrDefault<bool> + ( + "disableWallEdges", + false + ); + + // Predetermine mesh edges // ~~~~~~~~~~~~~~~~~~~~~~~ @@ -336,7 +352,16 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) } } - else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos)) + else if + ( + isMaxEdge + ( + pointWallDist, + edgeI, + minMedialAxisAngleCos, + disableWallEdges + ) + ) { // Both end points of edge have very different nearest wall // point. Mark both points as medial axis points. @@ -351,41 +376,56 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) // Calculate distance along edge const point& p0 = points[e[0]]; + const point& origin0 = pointWallDist[e[0]].origin(); const point& p1 = points[e[1]]; - scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec; - scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec; + const point& origin1 = pointWallDist[e[1]].origin(); + scalar dist0 = (p0-origin0) & eVec; + scalar dist1 = (origin1-p1) & eVec; scalar s = 0.5*(dist1+eMag+dist0); - point medialAxisPt; + point medialAxisPt(vector::max); if (s <= dist0) { - medialAxisPt = p0; + // Make sure point is not on wall. Note that this + // check used to be inside isMaxEdge. + if (magSqr((p0-origin0)) > Foam::sqr(SMALL)) + { + medialAxisPt = p0; + } } else if (s >= dist0+eMag) { - medialAxisPt = p1; + // Make sure point is not on wall. Note that this + // check used to be inside isMaxEdge. + if (magSqr((p1-origin1)) > Foam::sqr(SMALL)) + { + medialAxisPt = p1; + } } else { medialAxisPt = p0+(s-dist0)*eVec; } - forAll(e, ep) + if (medialAxisPt != vector::max) { - label pointI = e[ep]; - - if (!pointMedialDist[pointI].valid(dummyTrackData)) + forAll(e, ep) { - maxPoints.append(pointI); - maxInfo.append - ( - pointEdgePoint + label pointI = e[ep]; + + if (!pointMedialDist[pointI].valid(dummyTrackData)) + { + maxPoints.append(pointI); + maxInfo.append ( - medialAxisPt, //points[pointI], - magSqr(points[pointI]-medialAxisPt)//0.0 - ) - ); - pointMedialDist[pointI] = maxInfo.last(); + pointEdgePoint + ( + medialAxisPt, //points[pointI], + magSqr(points[pointI]-medialAxisPt)//0.0 + ) + ); + pointMedialDist[pointI] = maxInfo.last(); + } } } } @@ -1441,12 +1481,14 @@ void Foam::medialAxisMeshMover::calculateDisplacement //- Option 2: Look at component in the direction // of nearest (medial axis or static) point. - vector n = - patchDisp[patchPointI] - / (mag(patchDisp[patchPointI]) + VSMALL); - vector mVec = medialVec_[pointI]-mesh().points()[pointI]; - mVec /= mag(mVec)+VSMALL; - thicknessRatio *= (n&mVec); + const vector n = normalised(patchDisp[patchPointI]); + const vector mVec = + normalised + ( + medialVec_[pointI] - mesh().points()[pointI] + ); + + thicknessRatio *= (n & mVec); if (thicknessRatio > maxThicknessToMedialRatio) { diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H index 6fa9a02405abd3f50de1b05a0b26966bcf1a132c..8f3ea990f4cc8d08febdf74099c5f6b52e8bdaa0 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H @@ -114,7 +114,8 @@ class medialAxisMeshMover ( const List<PointData<vector>>& pointWallDist, const label edgeI, - const scalar minCos + const scalar minCos, + const bool checkWall ) const; //- Initialise medial axis. Uses pointDisplacement field only diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C index d0da764d781740750d278785420de09639a04139..da0ea0bee73d3decf77ae499a7d44a99c3f02933 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C @@ -172,8 +172,7 @@ void Foam::meshRefinement::calcNeighbourData forAll(faceCells, i) { // Extrapolate the face centre. - vector fn = faceAreas[i]; - fn /= mag(fn)+VSMALL; + const vector fn = normalised(faceAreas[i]); label own = faceCells[i]; label ownLevel = cellLevel[own]; @@ -2344,7 +2343,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMeshRegions label nRemove = findRegions ( mesh_, - mergeDistance_*vector(1,1,1), // perturbVec + mergeDistance_ * vector::one, // perturbVec locationsInMesh, locationsOutsideMesh, cellRegion.nRegions(), diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C index 30abfceaa118dcd3c30665b9c6c3d364a4d10580..5ca36811c73a440b3494d41dec8d623a2c102a39 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C @@ -1649,7 +1649,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk ( mesh_, cellRegion, - mergeDistance_*vector(1,1,1), + mergeDistance_ * vector::one, insidePoint ); @@ -1916,7 +1916,7 @@ void Foam::meshRefinement::findCellZoneTopo ( mesh_, cellRegion, - mergeDistance_*vector(1,1,1), + mergeDistance_ * vector::one, keepPoint ); @@ -3875,7 +3875,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh findRegions ( mesh_, - mergeDistance_*vector(1,1,1), // perturbVec + mergeDistance_ * vector::one, // perturbVec locationsInMesh, locationsOutsideMesh, cellRegion.nRegions(), diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C index 714aad7bc93f768f8eb51bab853c37c278bc05ab..8b147255b8259c6e247f9e33e98eabf3a4a8456f 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C @@ -243,8 +243,7 @@ Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells { label facei = candidateFaces[i]; - vector n = mesh_.faceAreas()[facei]; - n /= mag(n); + const vector n = normalised(mesh_.faceAreas()[facei]); label region = surfaces_.globalRegion ( @@ -296,7 +295,7 @@ bool Foam::meshRefinement::isCollapsedFace const scalar severeNonorthogonalityThreshold = ::cos(degToRad(maxNonOrtho)); - vector s = mesh_.faces()[facei].normal(points); + vector s = mesh_.faces()[facei].areaNormal(points); scalar magS = mag(s); // Check face area diff --git a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C index d9cb3ecddfa6c5b3ee3ee55673d611a1ab0db2c6..6a3c996c9adae5a16436595a4217ffd4eaab2f36 100644 --- a/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C +++ b/src/mesh/snappyHexMesh/refinementFeatures/refinementFeatures.C @@ -204,6 +204,14 @@ void Foam::refinementFeatures::read distances_[featI][j] = distLevels[j].first(); levels_[featI][j] = distLevels[j].second(); + if (levels_[featI][j] < 0) + { + FatalErrorInFunction + << "Feature " << featFileName + << " has illegal refinement level " << levels_[featI][j] + << exit(FatalError); + } + // Check in incremental order if (j > 0) { @@ -575,8 +583,8 @@ void Foam::refinementFeatures::findNearestEdge const treeDataEdge& td = tree.shapes(); const edge& e = td.edges()[nearInfo[sampleI].index()]; - nearNormal[sampleI] = e.vec(td.points()); - nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL; + + nearNormal[sampleI] = e.unitVec(td.points()); } } } @@ -638,8 +646,8 @@ void Foam::refinementFeatures::findNearestRegionEdge ); const edge& e = td.edges()[nearInfo[sampleI].index()]; - nearNormal[sampleI] = e.vec(td.points()); - nearNormal[sampleI] /= mag(nearNormal[sampleI])+VSMALL; + + nearNormal[sampleI] = e.unitVec(td.points()); } } } diff --git a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C index f693ca674a80c0e6d626713af444b7cf8aab5914..0726bd0a3a23b96b3144e6401a569f3c4715e4f6 100644 --- a/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C +++ b/src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C @@ -207,6 +207,7 @@ Foam::refinementSurfaces::refinementSurfaces ( globalMinLevel[surfI] < 0 || globalMaxLevel[surfI] < globalMinLevel[surfI] + || globalMaxLevel[surfI] < 0 || globalLevelIncr[surfI] < 0 ) { diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C index ea8c00447a1443db1e231e73f939faf7f1949dff..6713ba2d83253dbaadf32828b9e5b601e1bbeabd 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C @@ -55,6 +55,8 @@ void Foam::shellSurfaces::setAndCheckLevels const List<Tuple2<scalar, label>>& distLevels ) { + const searchableSurface& shell = allGeometry_[shells_[shellI]]; + if (modes_[shellI] != DISTANCE && distLevels.size() != 1) { FatalErrorInFunction @@ -73,6 +75,16 @@ void Foam::shellSurfaces::setAndCheckLevels distances_[shellI][j] = distLevels[j].first(); levels_[shellI][j] = distLevels[j].second(); + if (levels_[shellI][j] < -1) + { + FatalErrorInFunction + << "Shell " << shell.name() + << " has illegal refinement level " + << levels_[shellI][j] + << exit(FatalError); + } + + // Check in incremental order if (j > 0) { @@ -95,8 +107,6 @@ void Foam::shellSurfaces::setAndCheckLevels } } - const searchableSurface& shell = allGeometry_[shells_[shellI]]; - if (modes_[shellI] == DISTANCE) { Info<< "Refinement level according to distance to " diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C index a00e7804014ca33059eaa883e19d22dd20972b4b..a2e7969f90805d33f0bf6856cf6464716e64339e 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C @@ -1101,9 +1101,7 @@ void Foam::snappySnapDriver::detectNearSurfaces //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // //{ - // const scalar cos45 = Foam::cos(degToRad(45.0)); - // vector n(cos45, cos45, cos45); - // n /= mag(n); + // const vector n = normalised(vector::one); // // pointField start(14*pp.nPoints()); // pointField end(start.size()); @@ -2458,10 +2456,11 @@ void Foam::snappySnapDriver::detectWarpedFaces //Info<< "Splitting face:" << f << " into f0:" << f0 // << " f1:" << f1 << endl; - vector n0 = f0.normal(localPoints); - scalar n0Mag = mag(n0); - vector n1 = f1.normal(localPoints); - scalar n1Mag = mag(n1); + const vector n0 = f0.areaNormal(localPoints); + const scalar n0Mag = mag(n0); + + const vector n1 = f1.areaNormal(localPoints); + const scalar n1Mag = mag(n1); if (n0Mag > ROOTVSMALL && n1Mag > ROOTVSMALL) { diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C index 7f7a3ac556a71fd8e2f5cab95f3eed5b0194afca..24e02f5ebb0d2555d9f23fd453ca11b023640bd4 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriverFeature.C @@ -1860,9 +1860,9 @@ Foam::labelPair Foam::snappySnapDriver::findDiagonalAttraction isConcave ( compact0.centre(points0), - compact0.normal(points0), + compact0.areaNormal(points0), compact1.centre(points1), - compact1.normal(points1), + compact1.areaNormal(points1), concaveCos ) ) diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C index 0e503da082850722d4316e8e2b876c11a2f8d0e7..9b1659e4b9d70f4dacfb61d0a2b947ed5ebd3e20 100644 --- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C +++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/directAMI/directAMI.C @@ -87,7 +87,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds // second attempt: match by shooting a ray into the tgt face if (!found) { - const vector srcN = srcF.normal(srcPoints); + const vector srcN = srcF.areaNormal(srcPoints); for (const label tgtI : tgtNbr) { @@ -117,7 +117,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds { Pout<< "source face not found: id=" << srcI << " centre=" << srcCf[srcI] - << " normal=" << srcF.normal(srcPoints) + << " normal=" << srcF.areaNormal(srcPoints) << " points=" << srcF.points(srcPoints) << endl; @@ -128,7 +128,7 @@ void Foam::directAMI<SourcePatch, TargetPatch>::appendToDirectSeeds Pout<< "face id: " << tgtI << " centre=" << tgtF.centre(tgtPoints) - << " normal=" << tgtF.normal(tgtPoints) + << " normal=" << tgtF.areaNormal(tgtPoints) << " points=" << tgtF.points(tgtPoints) << endl; } diff --git a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C index 26ef1811a2f0891b8d8a3b7c8aaac12596de7177..795397fc2910430bbdb43872db34e63537718f37 100644 --- a/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C +++ b/src/meshTools/AMIInterpolation/patches/cyclicAMI/cyclicAMIPolyPatch/cyclicAMIPolyPatch.C @@ -197,8 +197,8 @@ void Foam::cyclicAMIPolyPatch::calcTransforms reduce(n0, maxMagSqrOp<point>()); reduce(n1, maxMagSqrOp<point>()); - n0 /= mag(n0) + VSMALL; - n1 /= mag(n1) + VSMALL; + n0.normalise(); + n1.normalise(); // Extended tensor from two local coordinate systems calculated // using normal and rotation axis @@ -367,14 +367,14 @@ void Foam::cyclicAMIPolyPatch::calcTransforms() vectorField half0Areas(half0.size()); forAll(half0, facei) { - half0Areas[facei] = half0[facei].normal(half0.points()); + half0Areas[facei] = half0[facei].areaNormal(half0.points()); } const cyclicAMIPolyPatch& half1 = neighbPatch(); vectorField half1Areas(half1.size()); forAll(half1, facei) { - half1Areas[facei] = half1[facei].normal(half1.points()); + half1Areas[facei] = half1[facei].areaNormal(half1.points()); } calcTransforms diff --git a/src/meshTools/cellFeatures/cellFeatures.C b/src/meshTools/cellFeatures/cellFeatures.C index 8e0ab780d19fefccbb5226eea69db35ba032e4a7..3f68579b516d219e7b1907f661d9c0bdca4f59ca 100644 --- a/src/meshTools/cellFeatures/cellFeatures.C +++ b/src/meshTools/cellFeatures/cellFeatures.C @@ -124,15 +124,11 @@ bool Foam::cellFeatures::isCellFeatureEdge // Check the angle between them by comparing the face normals. - vector n0 = mesh_.faceAreas()[face0]; - n0 /= mag(n0); - - vector n1 = mesh_.faceAreas()[face1]; - n1 /= mag(n1); + const vector n0 = normalised(mesh_.faceAreas()[face0]); + const vector n1 = normalised(mesh_.faceAreas()[face1]); scalar cosAngle = n0 & n1; - const edge& e = mesh_.edges()[edgeI]; const face& f0 = mesh_.faces()[face0]; @@ -167,10 +163,8 @@ bool Foam::cellFeatures::isCellFeatureEdge { return true; } - else - { - return false; - } + + return false; } @@ -429,13 +423,11 @@ bool Foam::cellFeatures::isFeaturePoint(const label edge0, const label edge1) const edge& e0 = mesh_.edges()[edge0]; - vector e0Vec = e0.vec(mesh_.points()); - e0Vec /= mag(e0Vec); + const vector e0Vec = e0.unitVec(mesh_.points()); const edge& e1 = mesh_.edges()[edge1]; - vector e1Vec = e1.vec(mesh_.points()); - e1Vec /= mag(e1Vec); + const vector e1Vec = e1.unitVec(mesh_.points()); scalar cosAngle; diff --git a/src/meshTools/indexedOctree/treeDataFace.C b/src/meshTools/indexedOctree/treeDataFace.C index 2548717c6914a96ef334209124bf85c70423dc4f..7efb12332b370b553cce333668b5990b5e5bee2e 100644 --- a/src/meshTools/indexedOctree/treeDataFace.C +++ b/src/meshTools/indexedOctree/treeDataFace.C @@ -250,14 +250,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType vector pointNormal(Zero); - forAll(pFaces, i) + for (const label facei : pFaces) { - if (isTreeFace_.test(pFaces[i])) + if (isTreeFace_.test(facei)) { - vector n = mesh_.faceAreas()[pFaces[i]]; - n /= mag(n) + VSMALL; - - pointNormal += n; + pointNormal += normalised(mesh_.faceAreas()[facei]); } } @@ -319,14 +316,11 @@ Foam::volumeType Foam::treeDataFace::getVolumeType vector edgeNormal(Zero); - forAll(eFaces, i) + for (const label facei : eFaces) { - if (isTreeFace_.test(eFaces[i])) + if (isTreeFace_.test(facei)) { - vector n = mesh_.faceAreas()[eFaces[i]]; - n /= mag(n) + VSMALL; - - edgeNormal += n; + edgeNormal += normalised(mesh_.faceAreas()[facei]); } } @@ -369,11 +363,8 @@ Foam::volumeType Foam::treeDataFace::getVolumeType vector ePrev = points[f[f.rcIndex(fp)]] - fc; vector eNext = points[f[f.fcIndex(fp)]] - fc; - vector nLeft = ePrev ^ e; - nLeft /= mag(nLeft) + VSMALL; - - vector nRight = e ^ eNext; - nRight /= mag(nRight) + VSMALL; + vector nLeft = normalised(ePrev ^ e); + vector nRight = normalised(e ^ eNext); if (debug & 2) { diff --git a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C index 04d67382580e86cddce90113e639314eb4dba78e..3fde178bcf1cea5cc3471129753281cd4110d79e 100644 --- a/src/meshTools/indexedOctree/treeDataPrimitivePatch.C +++ b/src/meshTools/indexedOctree/treeDataPrimitivePatch.C @@ -172,7 +172,7 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType // already have point. pointHit curHit = f.nearestPoint(sample, points); - const vector area = f.normal(points); + const vector area = f.areaNormal(points); const point& curPt = curHit.rawPoint(); // @@ -310,11 +310,8 @@ Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType vector ePrev = points[f[f.rcIndex(fp)]] - fc; vector eNext = points[f[f.fcIndex(fp)]] - fc; - vector nLeft = ePrev ^ e; - nLeft /= mag(nLeft) + VSMALL; - - vector nRight = e ^ eNext; - nRight /= mag(nRight) + VSMALL; + vector nLeft = normalised(ePrev ^ e); + vector nRight = normalised(e ^ eNext); if (debug & 2) { diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C index e49b2a76752ea18bcde782d04953ab0f49019a48..9a9c7a2fd101b18a716dc1e91c91250e5f5bf91d 100644 --- a/src/meshTools/meshSearch/meshSearch.C +++ b/src/meshTools/meshSearch/meshSearch.C @@ -810,8 +810,7 @@ Foam::List<Foam::pointIndexHit> Foam::meshSearch::intersections { DynamicList<pointIndexHit> hits; - vector edgeVec = pEnd - pStart; - edgeVec /= mag(edgeVec); + const vector edgeVec = normalised(pEnd - pStart); point pt = pStart; diff --git a/src/meshTools/meshTools/meshTools.C b/src/meshTools/meshTools/meshTools.C index 07fc46bec4784bdd753becb229767aee99a8bc53..96300c8cad31cf52cb3fa5c47b3736f027b10551 100644 --- a/src/meshTools/meshTools/meshTools.C +++ b/src/meshTools/meshTools/meshTools.C @@ -192,11 +192,7 @@ Foam::vector Foam::meshTools::normEdgeVec const label edgeI ) { - vector eVec = mesh.edges()[edgeI].vec(mesh.points()); - - eVec /= mag(eVec); - - return eVec; + return mesh.edges()[edgeI].unitVec(mesh.points()); } @@ -804,7 +800,7 @@ Foam::vector Foam::meshTools::edgeToCutDir edgeI = meshTools::walkFace(mesh, facei, edgeI, vertI, 2); } - avgVec /= mag(avgVec) + VSMALL; + avgVec.normalise(); return avgVec; } diff --git a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C index 05e7f05fa2ec0f733cd2f54fe4d140ce96ffba1a..65180a349999b380dcda8043eaa1a6d46d8f0a9e 100644 --- a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C +++ b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C @@ -579,8 +579,7 @@ bool Foam::primitiveMeshGeometry::checkFaceSkewness // Boundary faces: consider them to have only skewness error. // (i.e. treat as if mirror cell on other side) - vector faceNormal = faceAreas[facei]; - faceNormal /= mag(faceNormal) + VSMALL; + const vector faceNormal = normalised(faceAreas[facei]); vector dOwn = faceCentres[facei] - cellCentres[own[facei]]; @@ -767,8 +766,7 @@ bool Foam::primitiveMeshGeometry::checkFaceAngles const face& f = fcs[facei]; - vector faceNormal = faceAreas[facei]; - faceNormal /= mag(faceNormal) + VSMALL; + const vector faceNormal = normalised(faceAreas[facei]); // Get edge from f[0] to f[size-1]; vector ePrev(p[f.first()] - p[f.last()]); @@ -1023,13 +1021,11 @@ bool Foam::primitiveMeshGeometry::checkFaceTwist label nWarped = 0; - forAll(checkFaces, i) + for (const label facei : checkFaces) { - label facei = checkFaces[i]; - const face& f = fcs[facei]; - scalar magArea = mag(faceAreas[facei]); + const scalar magArea = mag(faceAreas[facei]); if (f.size() > 3 && magArea > VSMALL) { @@ -1039,21 +1035,21 @@ bool Foam::primitiveMeshGeometry::checkFaceTwist forAll(f, fpI) { - vector triArea + const vector triArea ( triPointRef ( p[f[fpI]], p[f.nextLabel(fpI)], fc - ).normal() + ).areaNormal() ); - scalar magTri = mag(triArea); + const scalar magTri = mag(triArea); if (magTri > VSMALL && ((nf & triArea/magTri) < minTwist)) { - nWarped++; + ++nWarped; if (setPtr) { diff --git a/src/meshTools/searchableSurfaces/searchableCone/searchableCone.C b/src/meshTools/searchableSurfaces/searchableCone/searchableCone.C index a28dc3600d857e2a6e207ff39d3b6a538641ef1b..1ffa831b510f1b506e646062779062078689e42f 100644 --- a/src/meshTools/searchableSurfaces/searchableCone/searchableCone.C +++ b/src/meshTools/searchableSurfaces/searchableCone/searchableCone.C @@ -107,15 +107,8 @@ void Foam::searchableCone::findNearestAndNormal // Remove the parallel component and normalise v -= parallel*unitDir_; - scalar magV = mag(v); - if (magV < ROOTVSMALL) - { - v = Zero; - } - else - { - v /= magV; - } + const scalar magV = mag(v); + v.normalise(); // Nearest and normal on disk at point1 point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v); @@ -143,31 +136,30 @@ void Foam::searchableCone::findNearestAndNormal p1 /= mag(p1); // Find vector along the two end of cone - vector b(projPt2 - projPt1); - scalar magS = mag(b); - b /= magS; + const vector b = normalised(projPt2 - projPt1); // Find the vector along sample pt and pt at one end of cone - vector a(sample - projPt1); + vector a = (sample - projPt1); if (mag(a) <= ROOTVSMALL) { // Exception: sample on disk1. Redo with projPt2. - vector a(sample - projPt2); + a = (sample - projPt2); + // Find normal unitvector - nearCone = (a & b)*b+projPt2; + nearCone = (a & b)*b + projPt2; + vector b1 = (p1 & b)*b; - normalCone = p1 - b1; - normalCone /= mag(normalCone); + normalCone = normalised(p1 - b1); } else { - // Find neartest point on cone surface - nearCone = (a & b)*b+projPt1; + // Find nearest point on cone surface + nearCone = (a & b)*b + projPt1; + // Find projection along surface of cone vector b1 = (p1 & b)*b; - normalCone = p1 - b1; - normalCone /= mag(normalCone); + normalCone = normalised(p1 - b1); } if (innerRadius1_ > 0 || innerRadius2_ > 0) @@ -176,13 +168,10 @@ void Foam::searchableCone::findNearestAndNormal point iCprojPt1 = point1_+ innerRadius1_*v; point iCprojPt2 = point2_+ innerRadius2_*v; - vector iCp1 = (iCprojPt1 - point1_); - iCp1 /= mag(iCp1); + const vector iCp1 = normalised(iCprojPt1 - point1_); // Find vector along the two end of cone - vector iCb(iCprojPt2 - iCprojPt1); - magS = mag(iCb); - iCb /= magS; + const vector iCb = normalised(iCprojPt2 - iCprojPt1); // Find the vector along sample pt and pt at one end of conde @@ -190,22 +179,22 @@ void Foam::searchableCone::findNearestAndNormal if (mag(iCa) <= ROOTVSMALL) { - vector iCa(sample - iCprojPt2); + iCa = (sample - iCprojPt2); // Find normal unitvector iCnearCone = (iCa & iCb)*iCb+iCprojPt2; + vector b1 = (iCp1 & iCb)*iCb; - iCnormalCone = iCp1 - b1; - iCnormalCone /= mag(iCnormalCone); + iCnormalCone = normalised(iCp1 - b1); } else { // Find nearest point on cone surface iCnearCone = (iCa & iCb)*iCb+iCprojPt1; + // Find projection along surface of cone vector b1 = (iCp1 & iCb)*iCb; - iCnormalCone = iCp1 - b1; - iCnormalCone /= mag(iCnormalCone); + iCnormalCone = normalised(iCp1 - b1); } } } @@ -234,8 +223,9 @@ void Foam::searchableCone::findNearestAndNormal scalar para = (v1 & unitDir_); // Remove the parallel component and normalise v1 -= para*unitDir_; - scalar magV1 = mag(v1); + const scalar magV1 = mag(v1); v1 = v1/magV1; + if (para < 0.0 && magV1 >= radius1_) { // Near point 1. Set point to intersection of disk and cone. @@ -281,7 +271,8 @@ void Foam::searchableCone::findNearestAndNormal scalar para = (v1 & unitDir_); // Remove the parallel component and normalise v1 -= para*unitDir_; - scalar magV1 = mag(v1); + + const scalar magV1 = mag(v1); v1 = v1/magV1; if (para < 0.0 && magV1 >= innerRadius1_) @@ -382,7 +373,7 @@ void Foam::searchableCone::findLineAll { // Find dot product: mag(s)>VSMALL suggest that it is greater - scalar s = (V&unitDir_); + scalar s = (V & unitDir_); if (mag(s) > VSMALL) { tPoint1 = -s1/s; @@ -471,8 +462,7 @@ void Foam::searchableCone::findLineAll else { vector va = cone.unitDir_; - vector v1(end-start); - v1 = v1/mag(v1); + vector v1 = normalised(end-start); scalar p = (va&v1); vector a1 = (v1-p*va); diff --git a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C index f531578355f402785ef4603e1d3ddc42555f383c..226c782d498a038923f3349d1b59252ca70bb9cc 100644 --- a/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C +++ b/src/meshTools/searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C @@ -89,8 +89,8 @@ Foam::searchableExtrudedCircle::searchableExtrudedCircle vector halfSpan(0.5*bounds().span()); point ctr(bounds().midpoint()); - bounds().min() = ctr - mag(halfSpan)*vector(1, 1, 1); - bounds().max() = ctr + mag(halfSpan)*vector(1, 1, 1); + bounds().min() = ctr - mag(halfSpan) * vector::one; + bounds().max() = ctr + mag(halfSpan) * vector::one; // Calculate bb of all points treeBoundBox bb(bounds()); @@ -182,8 +182,9 @@ void Foam::searchableExtrudedCircle::findNearest if (info[i].hit()) { - vector d(samples[i]-info[i].hitPoint()); - info[i].setPoint(info[i].hitPoint() + d/mag(d)*radius_); + const vector d = normalised(samples[i] - info[i].hitPoint()); + + info[i].setPoint(info[i].hitPoint() + d*radius_); } } } @@ -360,7 +361,8 @@ void Foam::searchableExtrudedCircle::findParametricNearest { radialStart = start-curvePoints[0]; radialStart -= (radialStart&axialVecs[0])*axialVecs[0]; - radialStart /= mag(radialStart); + radialStart.normalise(); + qStart = quaternion(radialStart, 0.0); info[0] = pointIndexHit(true, start, 0); @@ -370,11 +372,12 @@ void Foam::searchableExtrudedCircle::findParametricNearest { vector radialEnd(end-curvePoints.last()); radialEnd -= (radialEnd&axialVecs.last())*axialVecs.last(); - radialEnd /= mag(radialEnd); + radialEnd.normalise(); vector projectedEnd = radialEnd; projectedEnd -= (projectedEnd&axialVecs[0])*axialVecs[0]; - projectedEnd /= mag(projectedEnd); + projectedEnd.normalise(); + qProjectedEnd = quaternion(projectedEnd, 0.0); info.last() = pointIndexHit(true, end, 0); @@ -385,8 +388,8 @@ void Foam::searchableExtrudedCircle::findParametricNearest quaternion q(slerp(qStart, qProjectedEnd, lambdas[i])); vector radialDir(q.transform(radialStart)); - radialDir -= (radialDir&axialVecs[i])*axialVecs.last(); - radialDir /= mag(radialDir); + radialDir -= (radialDir & axialVecs[i]) * axialVecs.last(); + radialDir.normalise(); info[i] = pointIndexHit(true, curvePoints[i]+radius_*radialDir, 0); } @@ -432,10 +435,10 @@ void Foam::searchableExtrudedCircle::getNormal normal[i] = info[i].hitPoint()-curvePt.hitPoint(); // Subtract axial direction - vector axialVec = edges[curvePt.index()].vec(points); - axialVec /= mag(axialVec); - normal[i] -= (normal[i]&axialVec)*axialVec; - normal[i] /= mag(normal[i]); + const vector axialVec = edges[curvePt.index()].unitVec(points); + + normal[i] -= (normal[i] & axialVec) * axialVec; + normal[i].normalise(); } } } diff --git a/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C b/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C index 16b9430ad441c9a7c5104911453e7a2774bee89b..5691a89e31f75fa714003737e96798cdc88632ef 100644 --- a/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C +++ b/src/meshTools/searchableSurfaces/searchableSphere/searchableSphere.C @@ -336,9 +336,7 @@ void Foam::searchableSphere::getNormal { if (info[i].hit()) { - normal[i] = info[i].hitPoint() - centre_; - - normal[i] /= mag(normal[i])+VSMALL; + normal[i] = normalised(info[i].hitPoint() - centre_); } else { diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C index 752f245196606eaddb198c7633b3c6efd5f080a2..110f57feeb9b75ccdf069d64022a7dc4c689b01b 100644 --- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C +++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C @@ -752,8 +752,8 @@ void Foam::triSurfaceMesh::getNormal { if (info[i].hit()) { - label facei = info[i].index(); - normal[i] = s[facei].normal(pts); + const label facei = info[i].index(); + normal[i] = s[facei].unitNormal(pts); scalar qual = s[facei].tri(pts).quality(); @@ -768,12 +768,10 @@ void Foam::triSurfaceMesh::getNormal if (nbrQual > qual) { qual = nbrQual; - normal[i] = s[nbri].normal(pts); + normal[i] = s[nbri].unitNormal(pts); } } } - - normal[i] /= mag(normal[i]) + VSMALL; } else { @@ -789,12 +787,9 @@ void Foam::triSurfaceMesh::getNormal if (info[i].hit()) { const label facei = info[i].index(); - // Cached: - //normal[i] = faceNormals()[facei]; // Uncached - normal[i] = s[facei].normal(pts); - normal[i] /= mag(normal[i]) + VSMALL; + normal[i] = s[facei].unitNormal(pts); } else { diff --git a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C index 8ab9e1ac9494dcbdd73f21a5a8b2eac50796c31b..581214ba5489f8cb1873cb742d21a8aff93f0ddc 100644 --- a/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C +++ b/src/meshTools/sets/cellSources/rotatedBoxToCell/rotatedBoxToCell.C @@ -62,11 +62,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const boxPoints[6] = origin_ + k_ + i_ + j_; boxPoints[7] = origin_ + k_ + j_; - labelList boxVerts(8); - forAll(boxVerts, i) - { - boxVerts[i] = i; - } + labelList boxVerts(identity(8)); const cellModel& hex = cellModel::ref(cellModel::HEX); @@ -77,7 +73,7 @@ void Foam::rotatedBoxToCell::combine(topoSet& set, const bool add) const vectorField boxFaceNormals(boxFaces.size()); forAll(boxFaces, i) { - boxFaceNormals[i] = boxFaces[i].normal(boxPoints); + boxFaceNormals[i] = boxFaces[i].areaNormal(boxPoints); //Pout<< "Face:" << i << " position:" << boxFaces[i].centre(boxPoints) // << " normal:" << boxFaceNormals[i] << endl; diff --git a/src/meshTools/sets/faceSources/normalToFace/normalToFace.C b/src/meshTools/sets/faceSources/normalToFace/normalToFace.C index efc6d136ce7210e805f056d129142bf667ac79eb..d16998d626586f61e7d4f720d437844ebd31cf6d 100644 --- a/src/meshTools/sets/faceSources/normalToFace/normalToFace.C +++ b/src/meshTools/sets/faceSources/normalToFace/normalToFace.C @@ -51,7 +51,7 @@ Foam::topoSetSource::addToUsageTable Foam::normalToFace::usage_ void Foam::normalToFace::setNormal() { - normal_ /= mag(normal_) + VSMALL; + normal_.normalise(); Info<< " normalToFace : Normalized vector to " << normal_ << endl; @@ -116,8 +116,7 @@ void Foam::normalToFace::applyToSet forAll(mesh_.faceAreas(), facei) { - vector n = mesh_.faceAreas()[facei]; - n /= mag(n) + VSMALL; + const vector n = normalised(mesh_.faceAreas()[facei]); if (mag(1 - (n & normal_)) < tol_) { @@ -136,8 +135,7 @@ void Foam::normalToFace::applyToSet { const label facei = iter.key(); - vector n = mesh_.faceAreas()[facei]; - n /= mag(n) + VSMALL; + const vector n = normalised(mesh_.faceAreas()[facei]); if (mag(1 - (n & normal_)) < tol_) { diff --git a/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C b/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C index 391ba268250fbd0d06f9f8f5334c900e70fc48eb..3b4e6ed5b7a7978beca40adcbbcf66936279b8b8 100644 --- a/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C +++ b/src/meshTools/sets/faceZoneSources/setAndNormalToFaceZone/setAndNormalToFaceZone.C @@ -124,7 +124,7 @@ void Foam::setAndNormalToFaceZone::applyToSet { newAddressing.append(facei); - vector n = faces[facei].normal(points); + const vector n = faces[facei].areaNormal(points); if ((n & normal_) > 0) { newFlipMap.append(false); diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C index 34eded1ff313478509b37fc627b996297d358213..0852aa7a8d8f18596d77019544dcfd9772e3f638 100644 --- a/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C +++ b/src/meshTools/tetOverlapVolume/tetOverlapVolumeTemplates.C @@ -52,7 +52,7 @@ void Foam::tetOverlapVolume::tetTetOverlap tetPointRef tetATet = tetA.tet(); for (label facei = 0; facei < 4; ++facei) { - tetAFaceAreas[facei] = -tetATet.tri(facei).normal(); + tetAFaceAreas[facei] = -tetATet.tri(facei).areaNormal(); tetAMag2FaceAreas[facei] = magSqr(tetAFaceAreas[facei]); if (tetAMag2FaceAreas[facei] < ROOTVSMALL) { @@ -65,7 +65,7 @@ void Foam::tetOverlapVolume::tetTetOverlap tetPointRef tetBTet = tetB.tet(); for (label facei = 0; facei < 4; ++facei) { - tetBFaceAreas[facei] = -tetBTet.tri(facei).normal(); + tetBFaceAreas[facei] = -tetBTet.tri(facei).areaNormal(); tetBMag2FaceAreas[facei] = magSqr(tetBFaceAreas[facei]); if (tetBMag2FaceAreas[facei] < ROOTVSMALL) { diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C index a142a37ee6be4e4ad58b7d09d8e6b35c5a275e7e..9cfed82774976f44eba2e859cb5ed6559e257f9b 100644 --- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C +++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C @@ -433,11 +433,14 @@ Foam::label Foam::intersectedSurface::nextEdge const edge& prevE = edges[prevEdgei]; // x-axis of coordinate system - vector e0 = n ^ (points[prevE.otherVertex(prevVerti)] - points[prevVerti]); - e0 /= mag(e0) + VSMALL; + const vector e0 = + normalised + ( + n ^ (points[prevE.otherVertex(prevVerti)] - points[prevVerti]) + ); - // Get y-axis of coordinate system - vector e1 = n ^ e0; + // y-axis of coordinate system + const vector e1 = n ^ e0; if (mag(mag(e1) - 1) > SMALL) { @@ -489,10 +492,12 @@ Foam::label Foam::intersectedSurface::nextEdge ) { // Calculate angle of edge with respect to base e0, e1 - vector vec = - n ^ (points[e.otherVertex(prevVerti)] - points[prevVerti]); - - vec /= mag(vec) + VSMALL; + const vector vec = + normalised + ( + n + ^ (points[e.otherVertex(prevVerti)] - points[prevVerti]) + ); scalar angle = pseudoAngle(e0, e1, vec); @@ -1021,7 +1026,7 @@ Foam::faceList Foam::intersectedSurface::splitFace // See if normal needs flipping. faces.shrink(); - vector n = faces[0].normal(eSurf.points()); + const vector n = faces[0].areaNormal(eSurf.points()); if ((n & surf.faceNormals()[facei]) < 0) { diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C index 2f9101a461d7bb52db7c49b3d64dc3d0bed351e8..612ffb3bf18561771e2e550f23380c55994c7f21 100644 --- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C +++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C @@ -118,8 +118,7 @@ void Foam::edgeIntersections::intersectEdges const point& pStart = points1[meshPoints[e.start()]]; const point& pEnd = points1[meshPoints[e.end()]]; - const vector eVec(pEnd - pStart); - const vector n(eVec/(mag(eVec) + VSMALL)); + const vector n = normalised(pEnd - pStart); // Start tracking somewhat before pStart and up to somewhat after p1. // Note that tolerances here are smaller than those used to classify @@ -252,8 +251,7 @@ bool Foam::edgeIntersections::inlinePerturb label v0 = surf1.meshPoints()[e[0]]; label v1 = surf1.meshPoints()[e[1]]; - vector eVec(points1[v1] - points1[v0]); - vector n = eVec/mag(eVec); + const vector n = normalised(points1[v1] - points1[v0]); if (perturbStart) { @@ -326,9 +324,7 @@ bool Foam::edgeIntersections::rotatePerturb n /= magN; rndVec -= n*(n & rndVec); - - // Normalize - rndVec /= mag(rndVec) + VSMALL; + rndVec.normalise(); // Scale to be moved by tolerance. rndVec *= 0.01*magN; diff --git a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C index c083f60146787e2f8879ec4aa068542444acfdf6..8bf1e6332f0d15768db6e1f6ee03cc4187489474 100644 --- a/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C +++ b/src/meshTools/triSurface/faceTriangulation/faceTriangulation.C @@ -54,8 +54,8 @@ Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges const pointField& points ) { - tmp<vectorField> tedges(new vectorField(f.size())); - vectorField& edges = tedges.ref(); + auto tedges = tmp<vectorField>::New(f.size()); + auto& edges = tedges.ref(); forAll(f, i) { @@ -63,9 +63,8 @@ Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges point nextPt = points[f[f.fcIndex(i)]]; vector vec(nextPt - thisPt); - vec /= mag(vec) + VSMALL; - edges[i] = vec; + edges[i] = vec.normalise(); } return tedges; @@ -159,9 +158,9 @@ bool Foam::faceTriangulation::triangleContainsPoint const point& pt ) { - scalar area01Pt = triPointRef(p0, p1, pt).normal() & n; - scalar area12Pt = triPointRef(p1, p2, pt).normal() & n; - scalar area20Pt = triPointRef(p2, p0, pt).normal() & n; + const scalar area01Pt = triPointRef(p0, p1, pt).areaNormal() & n; + const scalar area12Pt = triPointRef(p1, p2, pt).areaNormal() & n; + const scalar area20Pt = triPointRef(p2, p0, pt).areaNormal() & n; if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0)) { @@ -172,10 +171,8 @@ bool Foam::faceTriangulation::triangleContainsPoint FatalErrorInFunction << abort(FatalError); return false; } - else - { - return false; - } + + return false; } @@ -209,7 +206,7 @@ void Foam::faceTriangulation::findDiagonal ); // rayDir should be normalized already but is not due to rounding errors // so normalize. - rayDir /= mag(rayDir) + VSMALL; + rayDir.normalise(); // @@ -309,8 +306,8 @@ void Foam::faceTriangulation::findDiagonal { // pt inside triangle (so perhaps visible) // Select based on minimal angle (so guaranteed visible). - vector edgePt0 = pt - startPt; - edgePt0 /= mag(edgePt0); + vector edgePt0 = (pt - startPt); + edgePt0.normalise(); scalar cos = rayDir & edgePt0; if (cos > maxCos) @@ -598,14 +595,12 @@ bool Foam::faceTriangulation::split // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -// Null constructor Foam::faceTriangulation::faceTriangulation() : triFaceList() {} -// Construct from components Foam::faceTriangulation::faceTriangulation ( const pointField& points, @@ -615,8 +610,7 @@ Foam::faceTriangulation::faceTriangulation : triFaceList(f.size()-2) { - vector avgNormal = f.normal(points); - avgNormal /= mag(avgNormal) + VSMALL; + const vector avgNormal = f.unitNormal(points); label triI = 0; @@ -629,7 +623,6 @@ Foam::faceTriangulation::faceTriangulation } -// Construct from components Foam::faceTriangulation::faceTriangulation ( const pointField& points, @@ -651,7 +644,6 @@ Foam::faceTriangulation::faceTriangulation } -// Construct from Istream Foam::faceTriangulation::faceTriangulation(Istream& is) : triFaceList(is) diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C index 4ef1974fe50f80f0abf89982b67a7b156f311548..abe113cceb7410c9b7ba4b093b145a3bd23930dd 100644 --- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C +++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C @@ -410,8 +410,7 @@ void Foam::triSurfaceSearch::findLineAll if (inter.hit()) { - vector lineVec = end[pointi] - start[pointi]; - lineVec /= mag(lineVec) + VSMALL; + const vector lineVec = normalised(end[pointi] - start[pointi]); if ( diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C index 5f87ba232edf27a987664b68ef5f1d90af251b7b..ac28865272bc3c0291efe196d58c568d0c657c7d 100644 --- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C +++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C @@ -107,8 +107,7 @@ Foam::pointToPointPlanarInterpolation::calcCoordinateSystem << exit(FatalError); } - vector n = e1^(points[index2]-p0); - n /= mag(n); + const vector n = normalised(e1 ^ (points[index2]-p0)); if (debug) { diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C index f64cb0eea2bf2d992f70eac5a851dd9cb3c7babd..e179edd35b7be16b27e4806e8ee3cb301714c639 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C @@ -68,8 +68,8 @@ Foam::triSurfaceTools::vertexNormals(const triSurface& surf) Info<< "Calculating vertex normals" << endl; - tmp<vectorField> tfld(new vectorField(surf.nPoints(), Zero)); - vectorField& pointNormals = tfld.ref(); + auto tpointNormals = tmp<vectorField>::New(surf.nPoints(), Zero); + auto& pointNormals = tpointNormals.ref(); const pointField& points = surf.points(); const labelListList& pointFaces = surf.pointFaces(); @@ -79,28 +79,27 @@ Foam::triSurfaceTools::vertexNormals(const triSurface& surf) { const labelList& pFaces = pointFaces[pI]; - forAll(pFaces, fI) + for (const label facei : pFaces) { - const label facei = pFaces[fI]; const triFace& f = surf[facei]; - vector fN = f.normal(points); + const vector areaNorm = f.areaNormal(points); const scalar weight = vertexNormalWeight ( f, meshPoints[pI], - fN, + areaNorm, points ); - pointNormals[pI] += weight*fN; + pointNormals[pI] += weight * areaNorm; } - pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL; + pointNormals[pI].normalise(); } - return tfld; + return tpointNormals; } @@ -114,8 +113,8 @@ Foam::triSurfaceTools::vertexTriads const pointField& points = surf.points(); const Map<label>& meshPointMap = surf.meshPointMap(); - tmp<triadField> tfld(new triadField(points.size())); - triadField& pointTriads = tfld.ref(); + auto tpointTriads = tmp<triadField>::New(points.size()); + auto& pointTriads = tpointTriads.ref(); forAll(points, pI) { @@ -131,16 +130,13 @@ Foam::triSurfaceTools::vertexTriads plane p(pt, normal); // Pick arbitrary point in plane - vector dir1 = pt - p.somePointInPlane(1e-3); - dir1 /= mag(dir1); - - vector dir2 = dir1 ^ normal; - dir2 /= mag(dir2); + vector dir1 = normalised(pt - p.somePointInPlane(1e-3)); + vector dir2 = normalised(dir1 ^ normal); pointTriads[meshPointMap[pI]] = triad(dir1, dir2, normal); } - return tfld; + return tpointTriads; } @@ -187,7 +183,7 @@ Foam::triSurfaceTools::curvatures // Set up a local coordinate system for the face const vector& e0 = edgeVectors[0]; - const vector eN = f.normal(points); + const vector eN = f.areaNormal(points); const vector e1 = (e0 ^ eN); if (magSqr(eN) < ROOTVSMALL) @@ -276,7 +272,7 @@ Foam::triSurfaceTools::curvatures ( f, meshPoints[patchPointIndex], - f.normal(points), + f.areaNormal(points), points ); @@ -304,8 +300,8 @@ Foam::triSurfaceTools::curvatures } } - tmp<scalarField> tfld(new scalarField(points.size(), Zero)); - scalarField& curvatureAtPoints = tfld.ref(); + auto tcurvatureAtPoints = tmp<scalarField>::New(points.size(), Zero); + scalarField& curvatureAtPoints = tcurvatureAtPoints.ref(); forAll(curvatureAtPoints, pI) { @@ -325,7 +321,7 @@ Foam::triSurfaceTools::curvatures curvatureAtPoints[meshPoints[pI]] = curvature; } - return tfld; + return tcurvatureAtPoints; } @@ -356,8 +352,8 @@ Foam::triSurfaceTools::writeCurvature { Info<< "Extracting curvature of surface at the points." << endl; - tmp<scalarField> tfld = triSurfaceTools::curvatures(surf); - scalarField& curv = tfld.ref(); + tmp<scalarField> tcurv = triSurfaceTools::curvatures(surf); + scalarField& curv = tcurv.ref(); triSurfacePointScalarField outputField ( @@ -379,7 +375,7 @@ Foam::triSurfaceTools::writeCurvature outputField.write(); outputField.swap(curv); - return tfld; + return tcurv; } diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 75808f63eaf5b127691cddc5051255360f3a259d..46c36f7deb3e144fad3d4dc00db94d3aed062b60 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -315,11 +315,8 @@ Foam::scalar Foam::triSurfaceTools::faceCosAngle const vector base0(pLeft - pStart); const vector base1(pRight - pStart); - vector n0(common ^ base0); - n0 /= Foam::mag(n0); - - vector n1(base1 ^ common); - n1 /= Foam::mag(n1); + const vector n0 = normalised(common ^ base0); + const vector n1 = normalised(base1 ^ common); return n0 & n1; } @@ -2068,11 +2065,11 @@ Foam::vector Foam::triSurfaceTools::surfaceNormal vector edgeNormal(Zero); - forAll(eFaces, i) + for (const label facei : eFaces) { - edgeNormal += surf.faceNormals()[eFaces[i]]; + edgeNormal += surf.faceNormals()[facei]; } - return edgeNormal/(mag(edgeNormal) + VSMALL); + return normalised(edgeNormal); } else { @@ -2626,14 +2623,13 @@ void Foam::triSurfaceTools::calcInterpolationWeights edge[1] = tri.a()-tri.c(); edge[2] = tri.b()-tri.a(); - vector triangleFaceNormal = edge[1] ^ edge[2]; + const vector triangleFaceNormal = edge[1] ^ edge[2]; // calculate edge normal (pointing inwards) FixedList<vector, 3> normal; for (label i=0; i<3; i++) { - normal[i] = triangleFaceNormal ^ edge[i]; - normal[i] /= mag(normal[i]) + VSMALL; + normal[i] = normalised(triangleFaceNormal ^ edge[i]); } weights[0] = ((p-tri.b()) & normal[0]) / max(VSMALL, normal[0] & edge[1]); diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C index 5cafe990b1873874c45451ddc1756fc1d13782b6..9e26175a1966dbcfeececf6ddd2c2e7b28b804f2 100644 --- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C +++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.C @@ -434,6 +434,48 @@ void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles } +bool Foam::cellCellStencils::inverseDistance::betterDonor +( + const label destMesh, + const label currentDonorMesh, + const label newDonorMesh +) const +{ + // This determines for multiple overlapping meshes which one provides + // the best donors. Is very basic and only looks at indices of meshes: + // - 'nearest' mesh index wins, i.e. on mesh 0 it preferentially uses donors + // from mesh 1 over mesh 2 (if applicable) + // - if same 'distance' the highest mesh wins. So on mesh 1 it + // preferentially uses donors from mesh 2 over mesh 0. This particular + // rule helps to avoid some interpolation loops where mesh 1 uses donors + // from mesh 0 (usually the background) but mesh 0 then uses + // donors from 1. + + if (currentDonorMesh == -1) + { + return true; + } + else + { + const label currentDist = mag(currentDonorMesh-destMesh); + const label newDist = mag(newDonorMesh-destMesh); + + if (newDist < currentDist) + { + return true; + } + else if (newDist == currentDist && newDonorMesh > currentDonorMesh) + { + return true; + } + else + { + return false; + } + } +} + + void Foam::cellCellStencils::inverseDistance::markDonors ( const globalIndex& globalCells, @@ -471,11 +513,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors // TBD: check for multiple donors. Maybe better one? For // now check 'nearer' mesh - if - ( - allStencil[celli].empty() - || (mag(srcI-tgtI) < mag(allDonor[celli]-tgtI)) - ) + if (betterDonor(tgtI, allDonor[celli], srcI)) { label globalDonor = globalCells.toGlobal(srcCellMap[srcCelli]); @@ -605,11 +643,7 @@ void Foam::cellCellStencils::inverseDistance::markDonors // TBD: check for multiple donors. Maybe better one? For // now check 'nearer' mesh - if - ( - allStencil[celli].empty() - || (mag(srcI-tgtI) < mag(allDonor[celli]-tgtI)) - ) + if (betterDonor(tgtI, allDonor[celli], srcI)) { allStencil[celli].setSize(1); allStencil[celli][0] = globalDonor; @@ -1642,6 +1676,7 @@ Foam::cellCellStencils::inverseDistance::inverseDistance // Protect local fields from interpolation nonInterpolatedFields_.insert("cellInterpolationWeight"); nonInterpolatedFields_.insert("cellTypes"); + nonInterpolatedFields_.insert("maxMagWeight"); // Read zoneID this->zoneID(); @@ -2066,6 +2101,9 @@ bool Foam::cellCellStencils::inverseDistance::update() if (debug&2) { + // Dump mesh + mesh_.time().write(); + // Dump stencil mkDir(mesh_.time().timePath()); OBJstream str(mesh_.time().timePath()/"injectionStencil.obj"); @@ -2100,29 +2138,79 @@ bool Foam::cellCellStencils::inverseDistance::update() cellInterpolationWeight_.instance() = mesh_.time().timeName(); cellInterpolationWeight_.write(); - // Dump cell types - volScalarField volTypes - ( - IOobject + // Dump max weight + { + volScalarField maxMagWeight ( - "cellTypes", - mesh_.time().timeName(), + IOobject + ( + "maxMagWeight", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar(dimless, Zero), - zeroGradientFvPatchScalarField::typeName - ); + dimensionedScalar(dimless, Zero), + zeroGradientFvPatchScalarField::typeName + ); + forAll(cellStencil_, celli) + { + const scalarList& wghts = cellInterpolationWeights_[celli]; + forAll(wghts, i) + { + if (mag(wghts[i]) > mag(maxMagWeight[celli])) + { + maxMagWeight[celli] = wghts[i]; + } + } + if (mag(maxMagWeight[celli]) > 1) + { + const pointField& cc = mesh_.cellCentres(); + Pout<< "cell:" << celli + << " at:" << cc[celli] + << " zone:" << zoneID[celli] + << " donors:" << cellStencil_[celli] + << " weights:" << wghts + << " coords:" + << UIndirectList<point>(cc, cellStencil_[celli]) + << " donorZone:" + << UIndirectList<label>(zoneID, cellStencil_[celli]) + << endl; + } + } + maxMagWeight.correctBoundaryConditions(); + maxMagWeight.write(); + } - forAll(volTypes.internalField(), cellI) + // Dump cell types { - volTypes[cellI] = cellTypes_[cellI]; + volScalarField volTypes + ( + IOobject + ( + "cellTypes", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh_, + dimensionedScalar(dimless, Zero), + zeroGradientFvPatchScalarField::typeName + ); + + forAll(volTypes.internalField(), cellI) + { + volTypes[cellI] = cellTypes_[cellI]; + } + volTypes.correctBoundaryConditions(); + volTypes.write(); } - volTypes.correctBoundaryConditions(); - volTypes.write(); + + // Dump stencil mkDir(mesh_.time().timePath()); diff --git a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H index 00fddad93f7c60260df73f4c14f63128f69fd2a3..22bf2531c8e3e9f45fcaaf55d982ca5d5bae9def 100644 --- a/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H +++ b/src/overset/cellCellStencil/inverseDistance/inverseDistanceCellCellStencil.H @@ -181,6 +181,14 @@ protected: labelList& allCellTypes ) const; + //- If multiple donors meshes: decide which is best + bool betterDonor + ( + const label destMesh, + const label currentDonorMesh, + const label newDonorMesh + ) const; + //- Determine donors for all tgt cells void markDonors ( diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C index 07e3fb0ed412a1f28229decf434cb08ddaf5e7ef..64766a1df7f018cb315b40c304be0f02ca7b46e7 100644 --- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C +++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C @@ -1891,8 +1891,7 @@ void Foam::distributedTriSurfaceMesh::getNormal forAll(triangleIndex, i) { label trii = triangleIndex[i]; - normal[i] = s[trii].normal(s.points()); - normal[i] /= mag(normal[i]) + VSMALL; + normal[i] = s[trii].unitNormal(s.points()); } diff --git a/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C b/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C index b3ee98eeff5f7fa509dc1a755fc0f1edb80144d4..92fe0e91160c23482fd77b7f6f484c70778a5e5c 100644 --- a/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C +++ b/src/rigidBodyDynamics/rigidBodyModelState/rigidBodyModelState.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2016-2017 OpenFOAM Foundation - \\/ M anipulation | + \\/ M anipulation | Copyright (C) 2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -51,7 +51,21 @@ Foam::RBD::rigidBodyModelState::rigidBodyModelState qDdot_(dict.lookupOrDefault("qDdot", scalarField(model.nDoF(), Zero))), t_(dict.lookupOrDefault<scalar>("t", -1)), deltaT_(dict.lookupOrDefault<scalar>("deltaT", 0)) -{} +{ + if + ( + q_.size() != model.nDoF() + || qDot_.size() != model.nDoF() + || qDdot_.size() != model.nDoF() + ) + { + FatalErrorInFunction << "State parameters 'q', 'qDot', 'qDdot'" + << " do not have the same size as the number of DoF " + << model.nDoF() + << ". Is your \"rigidBodyMotionState\" state file consistent?" + << exit(FatalError); + } +} // ************************************************************************* // diff --git a/src/sampling/sampledSet/circle/circleSet.C b/src/sampling/sampledSet/circle/circleSet.C index b5387e23ca47b821664dd423053efdf53cca849d..062bcdcb068d6986fff1397c393903afabf8150c 100644 --- a/src/sampling/sampledSet/circle/circleSet.C +++ b/src/sampling/sampledSet/circle/circleSet.C @@ -94,8 +94,7 @@ void Foam::circleSet::calcSamples label nPoint = 1; while (theta < 360) { - axis1 = axis1*cosAlpha + (axis1^circleAxis_)*sinAlpha; - axis1 /= mag(axis1); + axis1 = normalised(axis1*cosAlpha + (axis1^circleAxis_)*sinAlpha); point pt = origin_ + radius*axis1; label celli = searchEngine().findCell(pt); diff --git a/src/sampling/sampledSet/sampledSet/sampledSet.C b/src/sampling/sampledSet/sampledSet/sampledSet.C index cc053293305dca8956d93f37726281867f9207a9..856ce6a4016c186e0a74bfaa06e0e1f843a796c9 100644 --- a/src/sampling/sampledSet/sampledSet/sampledSet.C +++ b/src/sampling/sampledSet/sampledSet/sampledSet.C @@ -172,9 +172,7 @@ Foam::scalar Foam::sampledSet::calcSign vec /= magVec; - vector n = mesh().faceAreas()[facei]; - - n /= mag(n) + VSMALL; + const vector n = normalised(mesh().faceAreas()[facei]); return n & vec; } diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C index d6c0090be4898ac567ce3be758cdc9cca8b3a97f..6cecc030c2180a00fc5820411c3e65ec68fefa62 100644 --- a/src/sampling/surface/isoSurface/isoSurfaceCell.C +++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C @@ -295,8 +295,8 @@ Foam::pointIndexHit Foam::isoSurfaceCell::collapseSurface if (shared[0] != -1) { - const vector n0 = tri0.normal(localPoints); - const vector n1 = tri1.normal(localPoints); + const vector n0 = tri0.areaNormal(localPoints); + const vector n1 = tri1.areaNormal(localPoints); // Merge any zero-sized triangles, // or if they point in the same direction. diff --git a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C index fa421aa73dbc4100baaa973158156c6d4ea5714c..a6cf610c7466b1e90fd824b14b21c728696b1c78 100644 --- a/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C +++ b/src/surfMesh/surfaceFormats/stl/STLsurfaceFormat.C @@ -37,9 +37,9 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell const Face& f ) { - // calculate the normal ourselves, for flexibility and speed - vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal(); - norm /= mag(norm) + VSMALL; + // Calculate the normal ourselves, for flexibility and speed + const vector norm = + triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal(); // simple triangulation about f[0]. // better triangulation should have been done before @@ -71,9 +71,9 @@ inline void Foam::fileFormats::STLsurfaceFormat<Face>::writeShell const label zoneI ) { - // calculate the normal ourselves, for flexibility and speed - vector norm = triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).normal(); - norm /= mag(norm) + VSMALL; + // Calculate the normal ourselves, for flexibility and speed + const vector norm = + triPointRef(pts[f[0]], pts[f[1]], pts[f[2]]).unitNormal(); // simple triangulation about f[0]. // better triangulation should have been done before diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C index 85740809aec6adc2e56614990e84d4f12929e067..973482e8ddb603ffd7a577b11cc80b01d11b3943 100644 --- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C +++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/solarLoad.C @@ -294,10 +294,9 @@ void Foam::radiation::solarLoad::updateSkyDiffusiveRadiation void Foam::radiation::solarLoad::initialise(const dictionary& coeffs) { - if (coeffs.found("gridUp")) + if (coeffs.readIfPresent("gridUp", verticalDir_)) { - coeffs.lookup("gridUp") >> verticalDir_; - verticalDir_ /= mag(verticalDir_); + verticalDir_.normalise(); } else if (mesh_.foundObject<uniformDimensionedVectorField>("g")) { diff --git a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C index 7c720f5ab0a5981c915035471f22a2709e1eb73e..f95d517c188f42f99b7c20d8baa318b9127176f6 100644 --- a/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C +++ b/src/thermophysicalModels/radiation/submodels/solarCalculator/solarCalculator.C @@ -123,12 +123,8 @@ void Foam::solarCalculator::calculateBetaTetha() void Foam::solarCalculator::calculateSunDirection() { - - dict_.lookup("gridUp") >> gridUp_; - gridUp_ /= mag(gridUp_); - - dict_.lookup("gridEast") >> eastDir_; - eastDir_ /= mag(eastDir_); + gridUp_ = normalised(dict_.get<vector>("gridUp")); + eastDir_ = normalised(dict_.get<vector>("gridEast")); coord_.reset ( @@ -140,7 +136,7 @@ void Foam::solarCalculator::calculateSunDirection() direction_.y() = cos(beta_)*cos(tetha_); // South axis direction_.x() = cos(beta_)*sin(tetha_); // West axis - direction_ /= mag(direction_); + direction_.normalise(); if (debug) { @@ -163,10 +159,9 @@ void Foam::solarCalculator::init() { case mSunDirConstant: { - if (dict_.found("sunDirection")) + if (dict_.readIfPresent("sunDirection", direction_)) { - dict_.lookup("sunDirection") >> direction_; - direction_ /= mag(direction_); + direction_.normalise(); } else { diff --git a/src/waveModels/waveModel/waveModel.C b/src/waveModels/waveModel/waveModel.C index 324ac5ed5b8f75d37cbe91a8c9ad90c44a69183d..829cd89d59aee7233a77661edc1a56a77d694a96 100644 --- a/src/waveModels/waveModel/waveModel.C +++ b/src/waveModels/waveModel/waveModel.C @@ -57,8 +57,7 @@ void Foam::waveModel::initialiseGeometry() // - X: streamwise: patch normal // - Y: spanwise: Z^X // - Z: up: (negative) gravity direction - vector x(-gAverage(patch_.faceAreas())); - x /= mag(x) + ROOTVSMALL; + vector x = normalised(-gAverage(patch_.faceAreas())); vector z = -g_/mag(g_); vector y = z^x; diff --git a/tutorials/IO/dictionary/fatal-ending1.dict b/tutorials/IO/dictionary/fatal-ending1.dict index 2fa8e588d643af00eee8112e17ddfdd65e434559..76aa004c086b339a19314b68f0a9d28517de4ba6 100644 --- a/tutorials/IO/dictionary/fatal-ending1.dict +++ b/tutorials/IO/dictionary/fatal-ending1.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-ending2.dict b/tutorials/IO/dictionary/fatal-ending2.dict index fb4f85fd944bc5fa0b1a6c698ddc20171be406d6..1b6aa2b8df077c23ae3145e6f121bfae362224b9 100644 --- a/tutorials/IO/dictionary/fatal-ending2.dict +++ b/tutorials/IO/dictionary/fatal-ending2.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-ending3.dict b/tutorials/IO/dictionary/fatal-ending3.dict index 073e489efa63c34c9ece1c14974bc1cf11c73ab1..1359d7d51f407f7e780fbdaa22073eb882915118 100644 --- a/tutorials/IO/dictionary/fatal-ending3.dict +++ b/tutorials/IO/dictionary/fatal-ending3.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-ending4.dict b/tutorials/IO/dictionary/fatal-ending4.dict index 7f97b1578df1c962d7908d9ce2843fb2f8af5c75..b98fe3363fbe271fa2bd0ed6479f6037a422dc22 100644 --- a/tutorials/IO/dictionary/fatal-ending4.dict +++ b/tutorials/IO/dictionary/fatal-ending4.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-premature-ending1.dict b/tutorials/IO/dictionary/fatal-premature-ending1.dict index 40ee0b48db8518843c8be926ce1da5c0642f0bf8..50663197eda75c8e9da5d3e0ad4b1384c0e260fa 100644 --- a/tutorials/IO/dictionary/fatal-premature-ending1.dict +++ b/tutorials/IO/dictionary/fatal-premature-ending1.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-premature-ending2.dict b/tutorials/IO/dictionary/fatal-premature-ending2.dict index 7df892f725532938bd20945c20eff590e9a602ee..8a1bb29663af986a1668687ff20c08966ee2bf7d 100644 --- a/tutorials/IO/dictionary/fatal-premature-ending2.dict +++ b/tutorials/IO/dictionary/fatal-premature-ending2.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-primitive-ending1.dict b/tutorials/IO/dictionary/fatal-primitive-ending1.dict index d097357658e4704d17426dc8e0a7bbe98c7d5c73..f1d8fb07fb14a8927f5e946a2f05bcf596a6e811 100644 --- a/tutorials/IO/dictionary/fatal-primitive-ending1.dict +++ b/tutorials/IO/dictionary/fatal-primitive-ending1.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-primitive-ending2.dict b/tutorials/IO/dictionary/fatal-primitive-ending2.dict index 918718394926066304c082c43caf25a457f2191a..3fcd7cb303af0da525fb0fb91aa28a7022e21ab2 100644 --- a/tutorials/IO/dictionary/fatal-primitive-ending2.dict +++ b/tutorials/IO/dictionary/fatal-primitive-ending2.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/fatal-primitive-ending3.dict b/tutorials/IO/dictionary/fatal-primitive-ending3.dict index 0c326f6569b609c308b07398965a9c02cd687f6d..45fc555eb4828f9c913b7a51507632a6ed2a3141 100644 --- a/tutorials/IO/dictionary/fatal-primitive-ending3.dict +++ b/tutorials/IO/dictionary/fatal-primitive-ending3.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/good-empty2.dict b/tutorials/IO/dictionary/good-empty2.dict index 9c4a3f7a0a77be8f93d6d9154cf4e7b33cd5c57e..69d260404857e354d588063952ae5863e10db1a2 100644 --- a/tutorials/IO/dictionary/good-empty2.dict +++ b/tutorials/IO/dictionary/good-empty2.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/good-ending1.dict b/tutorials/IO/dictionary/good-ending1.dict index 29078dcbe883531b45bd79bd7c48f05b45fd982d..3eeb1c2b7d83689560cb39441cbbc7dfe66e5cf1 100644 --- a/tutorials/IO/dictionary/good-ending1.dict +++ b/tutorials/IO/dictionary/good-ending1.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/good-primitive-ending1.dict b/tutorials/IO/dictionary/good-primitive-ending1.dict index 62c70cee37ad7e4d6c92d3fbb6b9d8e25238d6dc..b46b73ce3b6f3034cada026e703463519ca33cdf 100644 --- a/tutorials/IO/dictionary/good-primitive-ending1.dict +++ b/tutorials/IO/dictionary/good-primitive-ending1.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/IO/dictionary/missed-ending3.dict b/tutorials/IO/dictionary/missed-ending3.dict index ec002dad47f8c9ba0256b7c6ca82f0744bca3fbe..bc53b15ab56b5c6f9b62f102dc90bcdba93ba1b5 100644 --- a/tutorials/IO/dictionary/missed-ending3.dict +++ b/tutorials/IO/dictionary/missed-ending3.dict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict b/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict index d66267b73f66121ca1354096f94588a09bac687f..348bce125d29eb05227b3e7bed563d0062a6ad2d 100644 --- a/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict +++ b/tutorials/basic/overPotentialFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict @@ -25,9 +25,4 @@ displacementLaplacianCoeffs dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - // ************************************************************************* // diff --git a/tutorials/compressible/overRhoSimpleFoam/hotCylinder/cylinderAndBackground/constant/dynamicMeshDict b/tutorials/compressible/overRhoSimpleFoam/hotCylinder/cylinderAndBackground/constant/dynamicMeshDict index d66267b73f66121ca1354096f94588a09bac687f..348bce125d29eb05227b3e7bed563d0062a6ad2d 100644 --- a/tutorials/compressible/overRhoSimpleFoam/hotCylinder/cylinderAndBackground/constant/dynamicMeshDict +++ b/tutorials/compressible/overRhoSimpleFoam/hotCylinder/cylinderAndBackground/constant/dynamicMeshDict @@ -25,9 +25,4 @@ displacementLaplacianCoeffs dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - // ************************************************************************* // diff --git a/tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict b/tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict index d66267b73f66121ca1354096f94588a09bac687f..348bce125d29eb05227b3e7bed563d0062a6ad2d 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict +++ b/tutorials/incompressible/overPimpleDyMFoam/cylinder/cylinderAndBackground/constant/dynamicMeshDict @@ -25,9 +25,4 @@ displacementLaplacianCoeffs dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - // ************************************************************************* // diff --git a/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/constant/dynamicMeshDict b/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/constant/dynamicMeshDict index de44bc5768ecd9631e05029e2da58e069b6822e8..9c10f9c182ee09b691355dc9239defd85325d8db 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/constant/dynamicMeshDict +++ b/tutorials/incompressible/overPimpleDyMFoam/simpleRotor/constant/dynamicMeshDict @@ -16,11 +16,6 @@ FoamFile dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - //motionSolverLibs ( "libfvMotionSolvers.so" ); // //solver displacementLaplacian; diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U index b1d4e2b1b0a9926d0fad92fe922e738dcc687a73..686016a479e433cf611ea29b492d314103d38fa4 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/U @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/epsilon b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/epsilon index 96e651d850e01ba7788cae39fd43064610ffcae6..983cc058b14e80463fb96c18be8f7135e5bc55e3 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/epsilon +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/epsilon @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/k b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/k index 70c5926719dcb361fd0868cb6e960c333aeda8de..d186d01d4fb88d56555b43ba813a0b86650a724d 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/k +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/k @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/nut b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/nut index 8254830002a70ffc3d912c9718c41f5afefb2eed..fa2790c57b7b855317beaaf6ca399d662ae34c2e 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/nut +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/nut @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p index 1818c10f3d491c49ec479edbf831ab74151567c6..c6900c3805dc9f395a1e7f5d03172c71098674ec 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/p @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/pointDisplacement b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/pointDisplacement index 7252992283e152d7ed13e531cf8598a2756dea1a..f20520b4046fc013bb7f626d672731efbac9d662 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/pointDisplacement +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/0.orig/pointDisplacement @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: plus | +| \\ / O peration | Version: v1806 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ diff --git a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/constant/dynamicMeshDict b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/constant/dynamicMeshDict index be7e563482406028d2a28b294a6b373397777d75..e95cf2990064eabe10c04666ac2dd0769f89851d 100644 --- a/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/constant/dynamicMeshDict +++ b/tutorials/incompressible/overPimpleDyMFoam/twoSimpleRotors/constant/dynamicMeshDict @@ -16,11 +16,6 @@ FoamFile dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - solver multiSolidBodyMotionSolver; multiSolidBodyMotionSolverCoeffs diff --git a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/constant/dynamicMeshDict b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/constant/dynamicMeshDict index d66267b73f66121ca1354096f94588a09bac687f..348bce125d29eb05227b3e7bed563d0062a6ad2d 100644 --- a/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/constant/dynamicMeshDict +++ b/tutorials/incompressible/overSimpleFoam/aeroFoil/background_overset/constant/dynamicMeshDict @@ -25,9 +25,4 @@ displacementLaplacianCoeffs dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -// layerRelax 0.3; -} - // ************************************************************************* // diff --git a/tutorials/incompressible/simpleFoam/windAroundBuildings/Allrun-parallel b/tutorials/incompressible/simpleFoam/windAroundBuildings/Allrun-parallel new file mode 100755 index 0000000000000000000000000000000000000000..444dc97cc96e17f0f4f92da7ad4b4befff57c704 --- /dev/null +++ b/tutorials/incompressible/simpleFoam/windAroundBuildings/Allrun-parallel @@ -0,0 +1,16 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory +. $WM_PROJECT_DIR/bin/tools/RunFunctions # Tutorial run functions + +runApplication surfaceFeatureExtract + +runApplication blockMesh +runApplication snappyHexMesh -overwrite + +runApplication decomposePar + +runParallel $(getApplication) + +runApplication reconstructPar + +#------------------------------------------------------------------------------ diff --git a/tutorials/incompressible/simpleFoam/windAroundBuildings/system/decomposeParDict b/tutorials/incompressible/simpleFoam/windAroundBuildings/system/decomposeParDict new file mode 100644 index 0000000000000000000000000000000000000000..c2c0de8e69c9e3fb0a1c60784624541423a31c61 --- /dev/null +++ b/tutorials/incompressible/simpleFoam/windAroundBuildings/system/decomposeParDict @@ -0,0 +1,23 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: v1806 | +| \\ / A nd | Web: www.OpenFOAM.com | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object decomposeParDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +numberOfSubdomains 6; + +// method kahip; +method scotch; + + +// ************************************************************************* // diff --git a/tutorials/incompressible/simpleFoam/windAroundBuildings/system/ensightWrite b/tutorials/incompressible/simpleFoam/windAroundBuildings/system/ensightWrite index 238ba0196394ddc145d3d86988c95aa2e6e3865d..4fb96ac6f44555e28e7a976e230cf5a5f3b2c554 100644 --- a/tutorials/incompressible/simpleFoam/windAroundBuildings/system/ensightWrite +++ b/tutorials/incompressible/simpleFoam/windAroundBuildings/system/ensightWrite @@ -12,6 +12,9 @@ ensightWrite // Fields to output (words or regex) fields (U p "(k|epsilon|omega)"); + // Limit output region + bounds (0 0 0) (245 180 80); + //- Write more frequent than fields writeControl timeStep; writeInterval 5; diff --git a/tutorials/mesh/parallel/cavity/system/decomposeParDict b/tutorials/mesh/parallel/cavity/system/decomposeParDict index 978d27ee1daa736885470f955538307240862852..7318ec1d80c51e66905dea144eb60ef0c506358a 100644 --- a/tutorials/mesh/parallel/cavity/system/decomposeParDict +++ b/tutorials/mesh/parallel/cavity/system/decomposeParDict @@ -107,6 +107,7 @@ structuredCoeffs method scotch; } +/* constraints { //- Keep owner and neighbour on same processor for faces in zones: @@ -148,7 +149,7 @@ constraints enabled false; } } - +*/ //// Is the case distributed? Note: command-line argument -roots takes //// precedence diff --git a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/constant/dynamicMeshDict b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/constant/dynamicMeshDict index 739c75ecc4414eb68145543c7550249c6e43a99d..103e65031d12de027c58045e9f4556192f88e071 100644 --- a/tutorials/multiphase/overInterDyMFoam/floatingBody/background/constant/dynamicMeshDict +++ b/tutorials/multiphase/overInterDyMFoam/floatingBody/background/constant/dynamicMeshDict @@ -18,10 +18,6 @@ motionSolverLibs ("libsixDoFRigidBodyMotion.so"); dynamicFvMesh dynamicOversetFvMesh; -dynamicOversetFvMeshCoeffs -{ -} - solver sixDoFRigidBodyMotion; sixDoFRigidBodyMotionCoeffs