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/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 ( 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++) {