From 82847ebbaf952a7ae420c8793a700bdce28bc0d5 Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 3 Nov 2022 11:32:32 +0000 Subject: [PATCH 1/2] plane: Raise errors only during dictionary construction This change prevents fatal errors occurring during programmatic construction of a plane object. If an invalid plane is constructed then this can be tested for using a new plane::valid() method. Errors are still generated from dictionary construction as before, and they have been improved to better identify where in the file the erroneous specification is. This change fixes some issues associated with meshToMesh mapping. The cell overlap calculation now detects and skips over degenerate tetrahedra. Previously, it was generating errors as it tried to construct planes from the faces of these degenerate tetrahedra. --- .../meshes/primitiveShapes/plane/plane.C | 205 +++++++----------- .../meshes/primitiveShapes/plane/plane.H | 53 +++-- .../tetOverlapVolume/tetOverlapVolume.C | 16 ++ 3 files changed, 131 insertions(+), 143 deletions(-) diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C index fdba29836a..6b62a8fb06 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -28,29 +28,37 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::plane::calcPntAndVec(const scalarList& C) +void Foam::plane::calcPntAndVec +( + const scalar a, + const scalar b, + const scalar c, + const scalar d +) { - normal_ = vector(C[0], C[1], C[2]); + normal_ = vector(a, b, c); const scalar magNormal = mag(normal_); - if (magNormal == 0) + // Normalise the normal if possible. Set to invalid if not. + if (magNormal > 0) { - FatalErrorInFunction - << "Plane normal has zero length" - << abort(FatalError); + normal_ /= magNormal; + } + else + { + normal_ = vector::zero; } - normal_ /= magNormal; - - if (magNormal < mag(C[3])*vSmall) + // Construct the point if possible. Set to far away if not. + if (magNormal > mag(d)*vSmall) { - FatalErrorInFunction - << "Plane is too far from the origin" - << abort(FatalError); + point_ = - d/magNormal*normal_; + } + else + { + point_ = point::max; } - - point_ = - C[3]/magNormal*normal_; } @@ -61,34 +69,9 @@ void Foam::plane::calcPntAndVec const point& point3 ) { + normal_ = normalised((point1 - point2) ^ (point2 - point3)); + point_ = (point1 + point2 + point3)/3; - vector line12 = point1 - point2; - vector line23 = point2 - point3; - - if - ( - mag(line12) < vSmall - || mag(line23) < vSmall - || mag(point3-point1) < vSmall - ) - { - FatalErrorInFunction - << "Bad points:" << point1 << ' ' << point2 << ' ' << point3 - << abort(FatalError); - } - - normal_ = line12 ^ line23; - scalar magUnitVector(mag(normal_)); - - if (magUnitVector < vSmall) - { - FatalErrorInFunction - << "Plane normal defined with zero length" << nl - << "Bad points:" << point1 << ' ' << point2 << ' ' << point3 - << abort(FatalError); - } - - normal_ /= magUnitVector; } @@ -96,58 +79,38 @@ void Foam::plane::calcPntAndVec Foam::plane::plane(const vector& normalVector) : - normal_(normalVector), + normal_(normalised(normalVector)), point_(Zero) -{ - scalar magUnitVector(mag(normal_)); - - if (magUnitVector > vSmall) - { - normal_ /= magUnitVector; - } - else - { - FatalErrorInFunction - << "plane normal has zero length. basePoint:" << point_ - << abort(FatalError); - } -} +{} Foam::plane::plane(const point& basePoint, const vector& normalVector) : - normal_(normalVector), + normal_(normalised(normalVector)), point_(basePoint) +{} + + +Foam::plane::plane +( + const scalar a, + const scalar b, + const scalar c, + const scalar d +) { - scalar magUnitVector(mag(normal_)); - - if (magUnitVector > vSmall) - { - normal_ /= magUnitVector; - } - else - { - FatalErrorInFunction - << "plane normal has zero length. basePoint:" << point_ - << abort(FatalError); - } -} - - -Foam::plane::plane(const scalarList& C) -{ - calcPntAndVec(C); + calcPntAndVec(a, b, c, d); } Foam::plane::plane ( - const point& a, - const point& b, - const point& c + const point& point1, + const point& point2, + const point& point3 ) { - calcPntAndVec(a, b, c); + calcPntAndVec(point1, point2, point3); } @@ -158,38 +121,35 @@ Foam::plane::plane(const dictionary& dict) { const word planeType(dict.lookup("planeType")); + const dictionary& subDict = dict.optionalSubDict(planeType + "Dict"); + if (planeType == "planeEquation") { - const dictionary& subDict = dict.optionalSubDict("planeEquationDict"); - scalarList C(4); - - C[0] = subDict.lookup("a"); - C[1] = subDict.lookup("b"); - C[2] = subDict.lookup("c"); - C[3] = subDict.lookup("d"); - - calcPntAndVec(C); - + calcPntAndVec + ( + subDict.lookup("a"), + subDict.lookup("b"), + subDict.lookup("c"), + subDict.lookup("d") + ); } else if (planeType == "embeddedPoints") { - const dictionary& subDict = dict.optionalSubDict("embeddedPointsDict"); - - point point1(subDict.lookup("point1")); - point point2(subDict.lookup("point2")); - point point3(subDict.lookup("point3")); - - calcPntAndVec(point1, point2, point3); + calcPntAndVec + ( + subDict.lookup("point1"), + subDict.lookup("point2"), + subDict.lookup("point3") + ); } else if (planeType == "pointAndNormal") { - const dictionary& subDict = dict.optionalSubDict("pointAndNormalDict"); - point_ = subDict.lookupBackwardsCompatible ( {"point", "basePoint"} ); + normal_ = normalised ( @@ -203,47 +163,36 @@ Foam::plane::plane(const dictionary& dict) { FatalIOErrorInFunction(dict) << "Invalid plane type: " << planeType << nl - << "Valid options include: planeEquation, embeddedPoints and " - << "pointAndNormal" + << "Valid options include: " + << "planeEquation, embeddedPoints and pointAndNormal" << abort(FatalIOError); } + + if (normal_ == vector::zero) + { + FatalIOErrorInFunction(subDict) + << "Plane normal has zero length" + << exit(FatalIOError); + } + + if (point_ == point::max) + { + FatalIOErrorInFunction(subDict) + << "Plane is too far from the origin" + << exit(FatalIOError); + } } Foam::plane::plane(Istream& is) : - normal_(is), + normal_(normalised(vector(is))), point_(is) -{ - scalar magUnitVector(mag(normal_)); - - if (magUnitVector > vSmall) - { - normal_ /= magUnitVector; - } - else - { - FatalErrorInFunction - << "plane normal has zero length. basePoint:" << point_ - << abort(FatalError); - } -} +{} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::vector& Foam::plane::normal() const -{ - return normal_; -} - - -const Foam::point& Foam::plane::refPoint() const -{ - return point_; -} - - Foam::FixedList Foam::plane::planeCoeffs() const { FixedList C(4); diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H index 9b1f9d4ab2..f02161dd1c 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -73,6 +73,7 @@ public: class ray { point refPoint_; + vector dir_; public: @@ -109,10 +110,15 @@ private: // Private Member Functions //- Calculates basePoint and normal vector given plane coefficients - void calcPntAndVec(const scalarList& C); + void calcPntAndVec + ( + const scalar C1, + const scalar C2, + const scalar C3, + const scalar C4 + ); //- Calculates basePoint and normal vector given three points - //- Normal vector determined using right hand rule void calcPntAndVec ( const point& point1, @@ -139,27 +145,44 @@ public: const point& point3 ); - //- Construct from coefficients for the - // plane equation: ax + by + cz + d = 0 - explicit plane(const scalarList& C); + //- Construct from coefficients for the plane equation, + // ax + by + cz + d = 0 + explicit plane + ( + const scalar a, + const scalar b, + const scalar c, + const scalar d + ); //- Construct from dictionary explicit plane(const dictionary& planeDict); - //- Construct from Istream. Assumes the base + normal notation. + //- Construct from Istream explicit plane(Istream& is); // Member Functions - //- Return plane normal - const vector& normal() const; + //- Return whether or not this plane is valid + inline bool valid() const + { + return normal_ != vector::zero && point_ != point::max; + } - //- Return or return plane base point - const point& refPoint() const; + //- Return the plane normal + inline const vector& normal() const + { + return normal_; + } - //- Return coefficients for the - // plane equation: ax + by + cz + d = 0 + //- Return the plane base point + inline const point& refPoint() const + { + return point_; + } + + //- Return coefficients for the plane equation, ax + by + cz + d = 0 FixedList planeCoeffs() const; //- Return a point on the plane @@ -176,7 +199,7 @@ public: scalar normalIntersect(const point& pnt0, const vector& dir) const; //- Return cut coefficient for plane and ray - scalar normalIntersect(const ray& r) const + inline scalar normalIntersect(const ray& r) const { return normalIntersect(r.refPoint(), r.dir()); } @@ -184,7 +207,7 @@ public: //- Return the cutting point between the plane and // a line passing through the supplied points template - scalar lineIntersect(const line& l) const + inline scalar lineIntersect(const line& l) const { return normalIntersect(l.start(), l.vec()); } diff --git a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C index 87dcfaa03e..6c80c96133 100644 --- a/src/meshTools/tetOverlapVolume/tetOverlapVolume.C +++ b/src/meshTools/tetOverlapVolume/tetOverlapVolume.C @@ -61,6 +61,10 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol // face 0 const plane pl0(tetB.b(), tetB.d(), tetB.c()); + if (!pl0.valid()) + { + return 0; + } const FixedList t({tetA.a(), tetA.b(), tetA.c(), tetA.d()}); cutTetList1.clear(); tetCut(t, pl0, cut::appendOp(cutTetList1), cut::noOp()); @@ -71,6 +75,10 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol // face 1 const plane pl1(tetB.a(), tetB.c(), tetB.d()); + if (!pl1.valid()) + { + return 0; + } cutTetList2.clear(); for (label i = 0; i < cutTetList1.size(); i++) { @@ -84,6 +92,10 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol // face 2 const plane pl2(tetB.a(), tetB.d(), tetB.b()); + if (!pl2.valid()) + { + return 0; + } cutTetList1.clear(); for (label i = 0; i < cutTetList2.size(); i++) { @@ -97,6 +109,10 @@ Foam::scalar Foam::tetOverlapVolume::tetTetOverlapVol // face 3 const plane pl3(tetB.a(), tetB.b(), tetB.c()); + if (!pl3.valid()) + { + return 0; + } scalar v = 0; for (label i = 0; i < cutTetList1.size(); i++) { From 1a7da6391ee3ca3db5b50956623522d58159832b Mon Sep 17 00:00:00 2001 From: Will Bainbridge Date: Thu, 3 Nov 2022 14:05:12 +0000 Subject: [PATCH 2/2] fvcDdt: Added missing one-field overloads --- src/finiteVolume/finiteVolume/fvc/fvcDdt.C | 63 +++++++++++++++------- src/finiteVolume/finiteVolume/fvc/fvcDdt.H | 56 ++++++++++--------- src/finiteVolume/finiteVolume/fvm/fvmDdt.C | 24 ++++----- src/finiteVolume/finiteVolume/fvm/fvmDdt.H | 14 ++--- 4 files changed, 96 insertions(+), 61 deletions(-) diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C index b74a843bbe..694428d301 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.C +++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.C @@ -102,6 +102,18 @@ ddt } +template +tmp> +ddt +( + const one&, + const GeometricField& vf +) +{ + return ddt(vf); +} + + template tmp> ddt @@ -125,25 +137,11 @@ ddt } -template -tmp> -ddt -( - const GeometricField& sf -) -{ - return fv::ddtScheme::New - ( - sf.mesh(), - sf.mesh().schemes().ddt("ddt(" + sf.name() + ')') - ).ref().fvcDdt(sf); -} - - template tmp> ddt ( + const one&, const one&, const GeometricField& vf ) @@ -156,11 +154,40 @@ template tmp> ddt ( - const GeometricField& vf, - const one& + const one&, + const volScalarField& rho, + const GeometricField& vf ) { - return ddt(vf); + return ddt(rho, vf); +} + + +template +tmp> +ddt +( + const volScalarField& alpha, + const one&, + const GeometricField& vf +) +{ + return ddt(alpha, vf); +} + + +template +tmp> +ddt +( + const GeometricField& sf +) +{ + return fv::ddtScheme::New + ( + sf.mesh(), + sf.mesh().schemes().ddt("ddt(" + sf.name() + ')') + ).ref().fvcDdt(sf); } diff --git a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H index 2331466120..874e49c9b3 100644 --- a/src/finiteVolume/finiteVolume/fvc/fvcDdt.H +++ b/src/finiteVolume/finiteVolume/fvc/fvcDdt.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | Website: https://openfoam.org - \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -80,6 +80,13 @@ namespace fvc const GeometricField& ); + template + tmp> ddt + ( + const one&, + const GeometricField& + ); + template tmp> ddt ( @@ -88,35 +95,36 @@ namespace fvc const GeometricField& ); + template + tmp> ddt + ( + const one&, + const one&, + const GeometricField& + ); + + template + tmp> ddt + ( + const one&, + const volScalarField&, + const GeometricField& + ); + + template + tmp> ddt + ( + const volScalarField&, + const one&, + const GeometricField& + ); + template tmp> ddt ( const GeometricField& ); - template - tmp> ddt - ( - const one&, - const GeometricField& - ); - - template - tmp> ddt - ( - const GeometricField&, - const one& - ); - - inline geometricZeroField ddt - ( - const one&, - const one& - ) - { - return geometricZeroField(); - } - template tmp < diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C index 18ed56e1fa..e5007162c3 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.C +++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.C @@ -55,18 +55,6 @@ ddt } -template -tmp> -ddt -( - const one&, - const GeometricField& vf -) -{ - return ddt(vf); -} - - template tmp> ddt @@ -99,6 +87,18 @@ ddt } +template +tmp> +ddt +( + const one&, + const GeometricField& vf +) +{ + return ddt(vf); +} + + template tmp> ddt diff --git a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H index bb401726be..0231f3afd3 100644 --- a/src/finiteVolume/finiteVolume/fvm/fvmDdt.H +++ b/src/finiteVolume/finiteVolume/fvm/fvmDdt.H @@ -57,13 +57,6 @@ namespace fvm const GeometricField& ); - template - tmp> ddt - ( - const one&, - const GeometricField& - ); - template tmp> ddt ( @@ -78,6 +71,13 @@ namespace fvm const GeometricField& ); + template + tmp> ddt + ( + const one&, + const GeometricField& + ); + template tmp> ddt (