/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd. \\/ M anipulation | ------------------------------------------------------------------------------- License This file is part of OpenFOAM. OpenFOAM is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. OpenFOAM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . \*---------------------------------------------------------------------------*/ #include "plane.H" #include "tensor.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // Calculate base point and unit normal vector from plane equation void Foam::plane::calcPntAndVec(const scalarList& C) { if (mag(C[0]) > VSMALL) { basePoint_ = vector((-C[3]/C[0]), 0, 0); } else { if (mag(C[1]) > VSMALL) { basePoint_ = vector(0, (-C[3]/C[1]), 0); } else { if (mag(C[2]) > VSMALL) { basePoint_ = vector(0, 0, (-C[3]/C[2])); } else { FatalErrorIn("void plane::calcPntAndVec(const scalarList&)") << "At least one plane coefficient must have a value" << abort(FatalError); } } } unitVector_ = vector(C[0], C[1], C[2]); scalar magUnitVector(mag(unitVector_)); if (magUnitVector < VSMALL) { FatalErrorIn("void plane::calcPntAndVec(const scalarList&)") << "Plane normal defined with zero length" << abort(FatalError); } unitVector_ /= magUnitVector; } void Foam::plane::calcPntAndVec ( const point& point1, const point& point2, const point& point3 ) { basePoint_ = (point1 + point2 + point3)/3; vector line12 = point1 - point2; vector line23 = point2 - point3; if ( mag(line12) < VSMALL || mag(line23) < VSMALL || mag(point3-point1) < VSMALL ) { FatalErrorIn ( "void plane::calcPntAndVec\n" "(\n" " const point&,\n" " const point&,\n" " const point&\n" ")\n" ) << "Bad points." << abort(FatalError); } unitVector_ = line12 ^ line23; scalar magUnitVector(mag(unitVector_)); if (magUnitVector < VSMALL) { FatalErrorIn ( "void plane::calcPntAndVec\n" "(\n" " const point&,\n" " const point&,\n" " const point&\n" ")\n" ) << "Plane normal defined with zero length" << abort(FatalError); } unitVector_ /= magUnitVector; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // Construct from normal vector through the origin Foam::plane::plane(const vector& normalVector) : unitVector_(normalVector), basePoint_(vector::zero) { scalar magUnitVector(mag(unitVector_)); if (magUnitVector > VSMALL) { unitVector_ /= magUnitVector; } else { FatalErrorIn("plane::plane(const vector&)") << "plane normal has zero length" << abort(FatalError); } } // Construct from point and normal vector Foam::plane::plane(const point& basePoint, const vector& normalVector) : unitVector_(normalVector), basePoint_(basePoint) { scalar magUnitVector(mag(unitVector_)); if (magUnitVector > VSMALL) { unitVector_ /= magUnitVector; } else { FatalErrorIn("plane::plane(const point&, const vector&)") << "plane normal has zero length" << abort(FatalError); } } // Construct from plane equation Foam::plane::plane(const scalarList& C) { calcPntAndVec(C); } // Construct from three points Foam::plane::plane ( const point& a, const point& b, const point& c ) { calcPntAndVec(a, b, c); } // Construct from dictionary Foam::plane::plane(const dictionary& dict) : unitVector_(vector::zero), basePoint_(point::zero) { word planeType(dict.lookup("planeType")); if (planeType == "planeEquation") { const dictionary& subDict = dict.subDict("planeEquationDict"); scalarList C(4); C[0] = readScalar(subDict.lookup("a")); C[1] = readScalar(subDict.lookup("b")); C[2] = readScalar(subDict.lookup("c")); C[3] = readScalar(subDict.lookup("d")); calcPntAndVec(C); } else if (planeType == "embeddedPoints") { const dictionary& subDict = dict.subDict("embeddedPoints"); point point1(subDict.lookup("point1")); point point2(subDict.lookup("point2")); point point3(subDict.lookup("point3")); calcPntAndVec(point1, point2, point3); } else if (planeType == "pointAndNormal") { const dictionary& subDict = dict.subDict("pointAndNormalDict"); basePoint_ = subDict.lookup("basePoint"); unitVector_ = subDict.lookup("normalVector"); unitVector_ /= mag(unitVector_); } else { FatalIOErrorIn ( "plane::plane(const dictionary&)", dict ) << "Invalid plane type: " << planeType << abort(FatalIOError); } } // Construct from Istream. Assumes point and normal vector. Foam::plane::plane(Istream& is) : unitVector_(is), basePoint_(is) { scalar magUnitVector(mag(unitVector_)); if (magUnitVector > VSMALL) { unitVector_ /= magUnitVector; } else { FatalErrorIn("plane::plane(Istream& is)") << "plane normal has zero length" << abort(FatalError); } } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // // Return plane normal vector const Foam::vector& Foam::plane::normal() const { return unitVector_; } // Return plane base point const Foam::point& Foam::plane::refPoint() const { return basePoint_; } // Return coefficcients for plane equation: ax + by + cz + d = 0 Foam::scalarList Foam::plane::planeCoeffs() const { scalarList C(4); scalar magX = mag(unitVector_.x()); scalar magY = mag(unitVector_.y()); scalar magZ = mag(unitVector_.z()); if (magX > magY) { if (magX > magZ) { C[0] = 1; C[1] = unitVector_.y()/unitVector_.x(); C[2] = unitVector_.z()/unitVector_.x(); } else { C[0] = 0; C[1] = 0; C[2] = 1; } } else { if (magY > magZ) { C[0] = 0; C[1] = 1; C[2] = unitVector_.z()/unitVector_.y(); } else { C[0] = 0; C[1] = 0; C[2] = 1; } } C[3] = - C[0] * basePoint_.x() - C[1] * basePoint_.y() - C[2] * basePoint_.z(); return C; } // Return nearest point in the plane for the given point Foam::point Foam::plane::nearestPoint(const point& p) const { return p - unitVector_*((p - basePoint_) & unitVector_); } // Return distance from the given point to the plane Foam::scalar Foam::plane::distance(const point& p) const { return mag((p - basePoint_) & unitVector_); } // Cutting point for plane and line defined by origin and direction Foam::scalar Foam::plane::normalIntersect ( const point& pnt0, const vector& dir ) const { scalar denom = stabilise((dir & unitVector_), VSMALL); return ((basePoint_ - pnt0) & unitVector_)/denom; } // Cutting line of two planes Foam::plane::ray Foam::plane::planeIntersect(const plane& plane2) const { // Mathworld plane-plane intersection. Assume there is a point on the // intersection line with z=0 and solve the two plane equations // for that (now 2x2 equation in x and y) // Better: use either z=0 or x=0 or y=0. const vector& n1 = normal(); const vector& n2 = plane2.normal(); const point& p1 = refPoint(); const point& p2 = plane2.refPoint(); scalar n1p1 = n1&p1; scalar n2p2 = n2&p2; vector dir = n1 ^ n2; // Determine zeroed out direction (can be x,y or z) by looking at which // has the largest component in dir. scalar magX = mag(dir.x()); scalar magY = mag(dir.y()); scalar magZ = mag(dir.z()); direction iZero, i1, i2; if (magX > magY) { if (magX > magZ) { iZero = 0; i1 = 1; i2 = 2; } else { iZero = 2; i1 = 0; i2 = 1; } } else { if (magY > magZ) { iZero = 1; i1 = 2; i2 = 0; } else { iZero = 2; i1 = 0; i2 = 1; } } vector pt; pt[iZero] = 0; pt[i1] = (n2[i2]*n1p1 - n1[i2]*n2p2) / (n1[i1]*n2[i2] - n2[i1]*n1[i2]); pt[i2] = (n2[i1]*n1p1 - n1[i1]*n2p2) / (n1[i2]*n2[i1] - n1[i1]*n2[i2]); return ray(pt, dir); } // Cutting point of three planes Foam::point Foam::plane::planePlaneIntersect ( const plane& plane2, const plane& plane3 ) const { List pcs(3); pcs[0]= planeCoeffs(); pcs[1]= plane2.planeCoeffs(); pcs[2]= plane3.planeCoeffs(); tensor a ( pcs[0][0],pcs[0][1],pcs[0][2], pcs[1][0],pcs[1][1],pcs[1][2], pcs[2][0],pcs[2][1],pcs[2][2] ); vector b(pcs[0][3],pcs[1][3],pcs[2][3]); return (inv(a) & (-b)); } // * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * // bool Foam::operator==(const plane& a, const plane& b) { if (a.basePoint_ == b.basePoint_ && a.unitVector_ == b.unitVector_) { return true; } else { return false; } } bool Foam::operator!=(const plane& a, const plane& b) { return !(a == b); } // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const plane& a) { os << a.unitVector_ << token::SPACE << a.basePoint_; return os; } // ************************************************************************* //