diff --git a/applications/test/momentOfInertia/Make/options b/applications/test/momentOfInertia/Make/options index 54c035b8f5..0b32f3355b 100644 --- a/applications/test/momentOfInertia/Make/options +++ b/applications/test/momentOfInertia/Make/options @@ -1,5 +1,6 @@ EXE_INC = \ - -I$(LIB_SRC)/meshTools/lnInclude + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/triSurface/lnInclude EXE_LIBS = \ -lmeshTools diff --git a/applications/test/momentOfInertia/Test-momentOfInertia.C b/applications/test/momentOfInertia/Test-momentOfInertia.C index 9a8ae40a5a..22327f1bdd 100644 --- a/applications/test/momentOfInertia/Test-momentOfInertia.C +++ b/applications/test/momentOfInertia/Test-momentOfInertia.C @@ -26,180 +26,38 @@ Application Description Calculates the inertia tensor and principal axes and moments of a - test face and tetrahedron. + test face, tetrahedron and mesh. \*---------------------------------------------------------------------------*/ +#include "argList.H" +#include "Time.H" +#include "polyMesh.H" #include "ListOps.H" #include "face.H" #include "tetPointRef.H" #include "triFaceList.H" #include "OFstream.H" #include "meshTools.H" +#include "momentOfInertia.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // using namespace Foam; -void massPropertiesSolid -( - const pointField& pts, - const triFaceList triFaces, - scalar density, - scalar& mass, - vector& cM, - tensor& J -) -{ - // Reimplemented from: Wm4PolyhedralMassProperties.cpp - // File Version: 4.10.0 (2009/11/18) - - // Geometric Tools, LC - // Copyright (c) 1998-2010 - // Distributed under the Boost Software License, Version 1.0. - // http://www.boost.org/LICENSE_1_0.txt - // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt - - // Boost Software License - Version 1.0 - August 17th, 2003 - - // Permission is hereby granted, free of charge, to any person or - // organization obtaining a copy of the software and accompanying - // documentation covered by this license (the "Software") to use, - // reproduce, display, distribute, execute, and transmit the - // Software, and to prepare derivative works of the Software, and - // to permit third-parties to whom the Software is furnished to do - // so, all subject to the following: - - // The copyright notices in the Software and this entire - // statement, including the above license grant, this restriction - // and the following disclaimer, must be included in all copies of - // the Software, in whole or in part, and all derivative works of - // the Software, unless such copies or derivative works are solely - // in the form of machine-executable object code generated by a - // source language processor. - - // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES - // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND - // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR - // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR - // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, - // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE - // USE OR OTHER DEALINGS IN THE SOFTWARE. - - const scalar r6 = 1.0/6.0; - const scalar r24 = 1.0/24.0; - const scalar r60 = 1.0/60.0; - const scalar r120 = 1.0/120.0; - - // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx - scalarField integrals(10, 0.0); - - forAll(triFaces, i) - { - const triFace& tri(triFaces[i]); - - // vertices of triangle i - vector v0 = pts[tri[0]]; - vector v1 = pts[tri[1]]; - vector v2 = pts[tri[2]]; - - // cross product of edges - vector eA = v1 - v0; - vector eB = v2 - v0; - vector n = eA ^ eB; - - // compute integral terms - scalar tmp0, tmp1, tmp2; - - scalar f1x, f2x, f3x, g0x, g1x, g2x; - - tmp0 = v0.x() + v1.x(); - f1x = tmp0 + v2.x(); - tmp1 = v0.x()*v0.x(); - tmp2 = tmp1 + v1.x()*tmp0; - f2x = tmp2 + v2.x()*f1x; - f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x; - g0x = f2x + v0.x()*(f1x + v0.x()); - g1x = f2x + v1.x()*(f1x + v1.x()); - g2x = f2x + v2.x()*(f1x + v2.x()); - - scalar f1y, f2y, f3y, g0y, g1y, g2y; - - tmp0 = v0.y() + v1.y(); - f1y = tmp0 + v2.y(); - tmp1 = v0.y()*v0.y(); - tmp2 = tmp1 + v1.y()*tmp0; - f2y = tmp2 + v2.y()*f1y; - f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y; - g0y = f2y + v0.y()*(f1y + v0.y()); - g1y = f2y + v1.y()*(f1y + v1.y()); - g2y = f2y + v2.y()*(f1y + v2.y()); - - scalar f1z, f2z, f3z, g0z, g1z, g2z; - - tmp0 = v0.z() + v1.z(); - f1z = tmp0 + v2.z(); - tmp1 = v0.z()*v0.z(); - tmp2 = tmp1 + v1.z()*tmp0; - f2z = tmp2 + v2.z()*f1z; - f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z; - g0z = f2z + v0.z()*(f1z + v0.z()); - g1z = f2z + v1.z()*(f1z + v1.z()); - g2z = f2z + v2.z()*(f1z + v2.z()); - - // update integrals - integrals[0] += n.x()*f1x; - integrals[1] += n.x()*f2x; - integrals[2] += n.y()*f2y; - integrals[3] += n.z()*f2z; - integrals[4] += n.x()*f3x; - integrals[5] += n.y()*f3y; - integrals[6] += n.z()*f3z; - integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x); - integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y); - integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z); - } - - integrals[0] *= r6; - integrals[1] *= r24; - integrals[2] *= r24; - integrals[3] *= r24; - integrals[4] *= r60; - integrals[5] *= r60; - integrals[6] *= r60; - integrals[7] *= r120; - integrals[8] *= r120; - integrals[9] *= r120; - - // mass - mass = integrals[0]; - - // center of mass - cM = vector(integrals[1], integrals[2], integrals[3])/mass; - - // inertia relative to origin - J.xx() = integrals[5] + integrals[6]; - J.xy() = -integrals[7]; - J.xz() = -integrals[9]; - J.yx() = J.xy(); - J.yy() = integrals[4] + integrals[6]; - J.yz() = -integrals[8]; - J.zx() = J.xz(); - J.zy() = J.yz(); - J.zz() = integrals[4] + integrals[5]; - - // inertia relative to center of mass - J -= mass*((cM & cM)*I - cM*cM); - - // Apply density - mass *= density; - J *= density; -} - - int main(int argc, char *argv[]) { + argList::addOption + ( + "cell", + "label", + "cell to use for inertia calculation, defaults to 0" + ); + + #include "setRootCase.H" + #include "createTime.H" + #include "createPolyMesh.H" + scalar density = 1.0; { @@ -286,16 +144,7 @@ int main(int argc, char *argv[]) vector cM = vector::zero; tensor J = tensor::zero; - massPropertiesSolid - ( - - pts, - tetFaces, - density, - m, - cM, - J - ); + momentOfInertia::massPropertiesSolid(pts, tetFaces, density, m, cM, J); vector eVal = eigenValues(J); @@ -344,7 +193,50 @@ int main(int argc, char *argv[]) { str << "l " << nPts + 1 << ' ' << i + 1 << endl; } + } + { + const label cellI = args.optionLookupOrDefault("cell", 0); + + tensorField mI = momentOfInertia::meshInertia(mesh); + + tensor& J = mI[cellI]; + + vector eVal = eigenValues(J); + + Info<< nl + << "Inertia tensor of cell " << cellI << " " << J << nl + << "eigenValues (principal moments) " << eVal << endl; + + J /= cmptMax(eVal); + + tensor eVec = eigenVectors(J); + + Info<< "eigenVectors (principal axes, from normalised inertia) " << eVec + << endl; + + OFstream str("cell_" + name(cellI) + "_inertia.obj"); + + Info<< nl << "Writing scaled principal axes of cell " << cellI << " to " + << str.name() << endl; + + const point& cC = mesh.cellCentres()[cellI]; + + scalar scale = mag + ( + (cC - mesh.faceCentres()[mesh.cells()[cellI][0]]) + /eVal.component(findMin(eVal)) + ); + + meshTools::writeOBJ(str, cC); + meshTools::writeOBJ(str, cC + scale*eVal.x()*eVec.x()); + meshTools::writeOBJ(str, cC + scale*eVal.y()*eVec.y()); + meshTools::writeOBJ(str, cC + scale*eVal.z()*eVec.z()); + + for (label i = 1; i < 4; i++) + { + str << "l " << 1 << ' ' << i + 1 << endl; + } } Info<< nl << "End" << nl << endl; diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C index f72bcb6e5e..2b8183295f 100644 --- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C +++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C @@ -33,9 +33,6 @@ Description #include "argList.H" #include "ListOps.H" -#include "face.H" -#include "tetPointRef.H" -#include "triFaceList.H" #include "triSurface.H" #include "OFstream.H" #include "meshTools.H" @@ -43,242 +40,12 @@ Description #include "transform.H" #include "IOmanip.H" #include "Pair.H" +#include "momentOfInertia.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // using namespace Foam; -void massPropertiesSolid -( - const pointField& pts, - const triFaceList& triFaces, - scalar density, - scalar& mass, - vector& cM, - tensor& J -) -{ - // Reimplemented from: Wm4PolyhedralMassProperties.cpp - // File Version: 4.10.0 (2009/11/18) - - // Geometric Tools, LC - // Copyright (c) 1998-2010 - // Distributed under the Boost Software License, Version 1.0. - // http://www.boost.org/LICENSE_1_0.txt - // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt - - // Boost Software License - Version 1.0 - August 17th, 2003 - - // Permission is hereby granted, free of charge, to any person or - // organization obtaining a copy of the software and accompanying - // documentation covered by this license (the "Software") to use, - // reproduce, display, distribute, execute, and transmit the - // Software, and to prepare derivative works of the Software, and - // to permit third-parties to whom the Software is furnished to do - // so, all subject to the following: - - // The copyright notices in the Software and this entire - // statement, including the above license grant, this restriction - // and the following disclaimer, must be included in all copies of - // the Software, in whole or in part, and all derivative works of - // the Software, unless such copies or derivative works are solely - // in the form of machine-executable object code generated by a - // source language processor. - - // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, - // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES - // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND - // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR - // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR - // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, - // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE - // USE OR OTHER DEALINGS IN THE SOFTWARE. - - const scalar r6 = 1.0/6.0; - const scalar r24 = 1.0/24.0; - const scalar r60 = 1.0/60.0; - const scalar r120 = 1.0/120.0; - - // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx - scalarField integrals(10, 0.0); - - forAll(triFaces, i) - { - const triFace& tri(triFaces[i]); - - // vertices of triangle i - vector v0 = pts[tri[0]]; - vector v1 = pts[tri[1]]; - vector v2 = pts[tri[2]]; - - // cross product of edges - vector eA = v1 - v0; - vector eB = v2 - v0; - vector n = eA ^ eB; - - // compute integral terms - scalar tmp0, tmp1, tmp2; - - scalar f1x, f2x, f3x, g0x, g1x, g2x; - - tmp0 = v0.x() + v1.x(); - f1x = tmp0 + v2.x(); - tmp1 = v0.x()*v0.x(); - tmp2 = tmp1 + v1.x()*tmp0; - f2x = tmp2 + v2.x()*f1x; - f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x; - g0x = f2x + v0.x()*(f1x + v0.x()); - g1x = f2x + v1.x()*(f1x + v1.x()); - g2x = f2x + v2.x()*(f1x + v2.x()); - - scalar f1y, f2y, f3y, g0y, g1y, g2y; - - tmp0 = v0.y() + v1.y(); - f1y = tmp0 + v2.y(); - tmp1 = v0.y()*v0.y(); - tmp2 = tmp1 + v1.y()*tmp0; - f2y = tmp2 + v2.y()*f1y; - f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y; - g0y = f2y + v0.y()*(f1y + v0.y()); - g1y = f2y + v1.y()*(f1y + v1.y()); - g2y = f2y + v2.y()*(f1y + v2.y()); - - scalar f1z, f2z, f3z, g0z, g1z, g2z; - - tmp0 = v0.z() + v1.z(); - f1z = tmp0 + v2.z(); - tmp1 = v0.z()*v0.z(); - tmp2 = tmp1 + v1.z()*tmp0; - f2z = tmp2 + v2.z()*f1z; - f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z; - g0z = f2z + v0.z()*(f1z + v0.z()); - g1z = f2z + v1.z()*(f1z + v1.z()); - g2z = f2z + v2.z()*(f1z + v2.z()); - - // update integrals - integrals[0] += n.x()*f1x; - integrals[1] += n.x()*f2x; - integrals[2] += n.y()*f2y; - integrals[3] += n.z()*f2z; - integrals[4] += n.x()*f3x; - integrals[5] += n.y()*f3y; - integrals[6] += n.z()*f3z; - integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x); - integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y); - integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z); - } - - integrals[0] *= r6; - integrals[1] *= r24; - integrals[2] *= r24; - integrals[3] *= r24; - integrals[4] *= r60; - integrals[5] *= r60; - integrals[6] *= r60; - integrals[7] *= r120; - integrals[8] *= r120; - integrals[9] *= r120; - - // mass - mass = integrals[0]; - - // center of mass - cM = vector(integrals[1], integrals[2], integrals[3])/mass; - - // inertia relative to origin - J.xx() = integrals[5] + integrals[6]; - J.xy() = -integrals[7]; - J.xz() = -integrals[9]; - J.yx() = J.xy(); - J.yy() = integrals[4] + integrals[6]; - J.yz() = -integrals[8]; - J.zx() = J.xz(); - J.zy() = J.yz(); - J.zz() = integrals[4] + integrals[5]; - - // inertia relative to center of mass - J -= mass*((cM & cM)*I - cM*cM); - - // Apply density - mass *= density; - J *= density; -} - - -void massPropertiesShell -( - const pointField& pts, - const triFaceList& triFaces, - scalar density, - scalar& mass, - vector& cM, - tensor& J -) -{ - // Reset properties for accumulation - - mass = 0.0; - cM = vector::zero; - J = tensor::zero; - - // Find centre of mass - - forAll(triFaces, i) - { - const triFace& tri(triFaces[i]); - - triPointRef t - ( - pts[tri[0]], - pts[tri[1]], - pts[tri[2]] - ); - - scalar triMag = t.mag(); - - cM += triMag*t.centre(); - - mass += triMag; - } - - cM /= mass; - - mass *= density; - - // Find inertia around centre of mass - - forAll(triFaces, i) - { - const triFace& tri(triFaces[i]); - - J += triPointRef - ( - pts[tri[0]], - pts[tri[1]], - pts[tri[2]] - ).inertia(cM, density); - } -} - - -tensor applyParallelAxisTheorem -( - scalar m, - const vector& cM, - const tensor& J, - const vector& refPt -) -{ - // The displacement vector (refPt = cM) is the displacement of the - // new reference point from the centre of mass of the body that - // the inertia tensor applies to. - - vector d = (refPt - cM); - - return J + m*((d & d)*I - d*d); -} - - int main(int argc, char *argv[]) { argList::addNote @@ -321,40 +88,17 @@ int main(int argc, char *argv[]) triSurface surf(surfFileName); - triFaceList faces(surf.size()); - - forAll(surf, i) - { - faces[i] = triFace(surf[i]); - } - scalar m = 0.0; vector cM = vector::zero; tensor J = tensor::zero; if (args.optionFound("shellProperties")) { - massPropertiesShell - ( - surf.points(), - faces, - density, - m, - cM, - J - ); + momentOfInertia::massPropertiesShell(surf, density, m, cM, J); } else { - massPropertiesSolid - ( - surf.points(), - faces, - density, - m, - cM, - J - ); + momentOfInertia::massPropertiesSolid(surf, density, m, cM, J); } if (m < 0) @@ -583,7 +327,7 @@ int main(int argc, char *argv[]) showTransform = false; } - Info<< nl << setprecision(10) + Info<< nl << setprecision(12) << "Density: " << density << nl << "Mass: " << m << nl << "Centre of mass: " << cM << nl @@ -615,7 +359,7 @@ int main(int argc, char *argv[]) if (calcAroundRefPt) { Info<< nl << "Inertia tensor relative to " << refPt << ": " << nl - << applyParallelAxisTheorem(m, cM, J, refPt) + << momentOfInertia::applyParallelAxisTheorem(m, cM, J, refPt) << endl; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H index 3e5ee425d6..3a7513770d 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H @@ -41,6 +41,7 @@ SourceFiles #include "tetPointRef.H" #include "triPointRef.H" #include "polyMesh.H" +#include "triFace.H" #include "face.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -146,6 +147,10 @@ public: // mesh face for this tet from the supplied mesh inline triPointRef faceTri(const polyMesh& mesh) const; + //- Return the point indices corresponding to the tri on the mesh + // face for this tet from the supplied mesh + inline triFace faceTriIs(const polyMesh& mesh) const; + //- Return the geometry corresponding to the tri on the // mesh face for this tet from the supplied mesh using // the old position diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H index 6b8a7871a0..9a4287b5c4 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H @@ -122,6 +122,21 @@ Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const } +Foam::triFace Foam::tetIndices::faceTriIs(const polyMesh& mesh) const +{ + const faceList& pFaces = mesh.faces(); + + const Foam::face& f = pFaces[faceI_]; + + return triFace + ( + f[faceBasePtI_], + f[facePtAI_], + f[facePtBI_] + ); +} + + Foam::triPointRef Foam::tetIndices::oldFaceTri(const polyMesh& mesh) const { const pointField& oldPPts = mesh.oldPoints(); diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index 0a31e2b362..3a27ffd29d 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -128,6 +128,7 @@ $(cellZoneSources)/setToCellZone/setToCellZone.C pointZoneSources = sets/pointZoneSources $(pointZoneSources)/setToPointZone/setToPointZone.C +momentOfInertia/momentOfInertia.C surfaceSets/surfaceSets.C diff --git a/src/meshTools/momentOfInertia/momentOfInertia.C b/src/meshTools/momentOfInertia/momentOfInertia.C index af3265a5ec..8dc1a53d61 100644 --- a/src/meshTools/momentOfInertia/momentOfInertia.C +++ b/src/meshTools/momentOfInertia/momentOfInertia.C @@ -21,382 +21,328 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Class - momentOfInertia - -Description - Reimplementation of volInt.c by Brian Mirtich. - * mirtich@cs.berkeley.edu * - * http://www.cs.berkeley.edu/~mirtich * - -------------------------------------------------------------------------------- -*/ +\*---------------------------------------------------------------------------*/ #include "momentOfInertia.H" -//#include "pyramidPointFaceRef.H" -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // -//Foam::tensor Foam::momentOfInertia -//( -// const pointField& points, -// const faceList& faces, -// const cell& cFaces, -// const point& cc -//) -//{ -// tensor t(tensor::zero); -// -// forAll(cFaces, i) -// { -// const face& f = faces[cFaces[i]]; -// -// scalar pyrVol = pyramidPointFaceRef(f, cc).mag(points); -// -// vector pyrCentre = pyramidPointFaceRef(f, cc).centre(points); -// -// vector d = pyrCentre - cc; -// -// t.xx() += pyrVol*(sqr(d.y()) + sqr(d.z())); -// t.yy() += pyrVol*(sqr(d.x()) + sqr(d.z())); -// t.zz() += pyrVol*(sqr(d.x()) + sqr(d.y())); -// -// t.xy() -= pyrVol*d.x()*d.y(); -// t.xz() -= pyrVol*d.x()*d.z(); -// t.yz() -= pyrVol*d.y()*d.z(); -// } -// -// // Symmetric -// t.yx() = t.xy(); -// t.zx() = t.xz(); -// t.zy() = t.yz(); -// -// return t; -//} - - -#define sqr(x) ((x)*(x)) -#define pow3(x) ((x)*(x)*(x)) - -// compute various integrations over projection of face -void Foam::compProjectionIntegrals +void Foam::momentOfInertia::massPropertiesSolid ( - const pointField& points, - const face& f, - const direction A, - const direction B, - - scalar& P1, - scalar& Pa, - scalar& Pb, - scalar& Paa, - scalar& Pab, - scalar& Pbb, - scalar& Paaa, - scalar& Paab, - scalar& Pabb, - scalar& Pbbb -) -{ - P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0; - - forAll(f, i) - { - scalar a0 = points[f[i]][A]; - scalar b0 = points[f[i]][B]; - scalar a1 = points[f[(i+1) % f.size()]][A]; - scalar b1 = points[f[(i+1) % f.size()]][B]; - scalar da = a1 - a0; - scalar db = b1 - b0; - - scalar a0_2 = a0 * a0; - scalar a0_3 = a0_2 * a0; - scalar a0_4 = a0_3 * a0; - - scalar b0_2 = b0 * b0; - scalar b0_3 = b0_2 * b0; - scalar b0_4 = b0_3 * b0; - - scalar a1_2 = a1 * a1; - scalar a1_3 = a1_2 * a1; - - scalar b1_2 = b1 * b1; - scalar b1_3 = b1_2 * b1; - - scalar C1 = a1 + a0; - - scalar Ca = a1*C1 + a0_2; - scalar Caa = a1*Ca + a0_3; - scalar Caaa = a1*Caa + a0_4; - - scalar Cb = b1*(b1 + b0) + b0_2; - scalar Cbb = b1*Cb + b0_3; - scalar Cbbb = b1*Cbb + b0_4; - - scalar Cab = 3*a1_2 + 2*a1*a0 + a0_2; - scalar Kab = a1_2 + 2*a1*a0 + 3*a0_2; - - scalar Caab = a0*Cab + 4*a1_3; - scalar Kaab = a1*Kab + 4*a0_3; - - scalar Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3; - scalar Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3; - - P1 += db*C1; - Pa += db*Ca; - Paa += db*Caa; - Paaa += db*Caaa; - Pb += da*Cb; - Pbb += da*Cbb; - Pbbb += da*Cbbb; - Pab += db*(b1*Cab + b0*Kab); - Paab += db*(b1*Caab + b0*Kaab); - Pabb += da*(a1*Cabb + a0*Kabb); - } - - P1 /= 2.0; - Pa /= 6.0; - Paa /= 12.0; - Paaa /= 20.0; - Pb /= -6.0; - Pbb /= -12.0; - Pbbb /= -20.0; - Pab /= 24.0; - Paab /= 60.0; - Pabb /= -60.0; -} - - -void Foam::compFaceIntegrals -( - const pointField& points, - const face& f, - const vector& n, - const scalar w, - const direction A, - const direction B, - const direction C, - - scalar& Fa, - scalar& Fb, - scalar& Fc, - scalar& Faa, - scalar& Fbb, - scalar& Fcc, - scalar& Faaa, - scalar& Fbbb, - scalar& Fccc, - scalar& Faab, - scalar& Fbbc, - scalar& Fcca -) -{ - scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb; - - compProjectionIntegrals - ( - points, - f, - A, - B, - - P1, - Pa, - Pb, - Paa, - Pab, - Pbb, - Paaa, - Paab, - Pabb, - Pbbb - ); - - scalar k1 = 1 / n[C]; - scalar k2 = k1 * k1; - scalar k3 = k2 * k1; - scalar k4 = k3 * k1; - - Fa = k1 * Pa; - Fb = k1 * Pb; - Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1); - - Faa = k1 * Paa; - Fbb = k1 * Pbb; - Fcc = k3 * (sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb - + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1)); - - Faaa = k1 * Paaa; - Fbbb = k1 * Pbbb; - Fccc = -k4 * (pow3(n[A])*Paaa + 3*sqr(n[A])*n[B]*Paab - + 3*n[A]*sqr(n[B])*Pabb + pow3(n[B])*Pbbb - + 3*w*(sqr(n[A])*Paa + 2*n[A]*n[B]*Pab + sqr(n[B])*Pbb) - + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1)); - - Faab = k1 * Paab; - Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb); - Fcca = k3 * (sqr(n[A])*Paaa + 2*n[A]*n[B]*Paab + sqr(n[B])*Pabb - + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa)); -} - - -void Foam::compVolumeIntegrals -( - const pointField& points, - const faceList& faces, - const cell& cFaces, - const vectorField& fNorm, - const scalarField& fW, - - scalar& T0, - vector& T1, - vector& T2, - vector& TP -) -{ - T0 = 0; - T1 = vector::zero; - T2 = vector::zero; - TP = vector::zero; - - forAll(cFaces, i) - { - const vector& n = fNorm[i]; - - scalar nx = mag(n[0]); - scalar ny = mag(n[1]); - scalar nz = mag(n[2]); - - direction A, B, C; - - if (nx > ny && nx > nz) - { - C = 0; - } - else - { - C = (ny > nz) ? 1 : 2; - } - - A = (C + 1) % 3; - B = (A + 1) % 3; - - scalar Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca; - compFaceIntegrals - ( - points, - faces[cFaces[i]], - n, - fW[i], - A, - B, - C, - - Fa, - Fb, - Fc, - Faa, - Fbb, - Fcc, - Faaa, - Fbbb, - Fccc, - Faab, - Fbbc, - Fcca - ); - - T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc)); - - T1[A] += n[A] * Faa; - T1[B] += n[B] * Fbb; - T1[C] += n[C] * Fcc; - - T2[A] += n[A] * Faaa; - T2[B] += n[B] * Fbbb; - T2[C] += n[C] * Fccc; - - TP[A] += n[A] * Faab; - TP[B] += n[B] * Fbbc; - TP[C] += n[C] * Fcca; - } - - T1 /= 2; - T2 /= 3; - TP /= 2; -} - - -// Calculate -// - r: centre of mass -// - J: inertia around origin (point 0,0,0) -void Foam::momentOfIntertia -( - const pointField& points, - const faceList& faces, - const cell& cFaces, - point& r, + const pointField& pts, + const triFaceList& triFaces, + scalar density, + scalar& mass, + vector& cM, tensor& J ) { - // Face normals - vectorField fNorm(cFaces.size()); - scalarField fW(cFaces.size()); + // Reimplemented from: Wm4PolyhedralMassProperties.cpp + // File Version: 4.10.0 (2009/11/18) - forAll(cFaces, i) + // Geometric Tools, LC + // Copyright (c) 1998-2010 + // Distributed under the Boost Software License, Version 1.0. + // http://www.boost.org/LICENSE_1_0.txt + // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt + + // Boost Software License - Version 1.0 - August 17th, 2003 + + // Permission is hereby granted, free of charge, to any person or + // organization obtaining a copy of the software and accompanying + // documentation covered by this license (the "Software") to use, + // reproduce, display, distribute, execute, and transmit the + // Software, and to prepare derivative works of the Software, and + // to permit third-parties to whom the Software is furnished to do + // so, all subject to the following: + + // The copyright notices in the Software and this entire + // statement, including the above license grant, this restriction + // and the following disclaimer, must be included in all copies of + // the Software, in whole or in part, and all derivative works of + // the Software, unless such copies or derivative works are solely + // in the form of machine-executable object code generated by a + // source language processor. + + // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, + // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES + // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND + // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR + // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR + // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, + // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE + // USE OR OTHER DEALINGS IN THE SOFTWARE. + + const scalar r6 = 1.0/6.0; + const scalar r24 = 1.0/24.0; + const scalar r60 = 1.0/60.0; + const scalar r120 = 1.0/120.0; + + // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx + scalarField integrals(10, 0.0); + + forAll(triFaces, i) { - label faceI = cFaces[i]; + const triFace& tri(triFaces[i]); - const face& f = faces[faceI]; + // vertices of triangle i + vector v0 = pts[tri[0]]; + vector v1 = pts[tri[1]]; + vector v2 = pts[tri[2]]; - fNorm[i] = f.normal(points); - fNorm[i] /= mag(fNorm[i]) + VSMALL; + // cross product of edges + vector eA = v1 - v0; + vector eB = v2 - v0; + vector n = eA ^ eB; - fW[i] = - (fNorm[i] & points[f[0]]); + // compute integral terms + scalar tmp0, tmp1, tmp2; + + scalar f1x, f2x, f3x, g0x, g1x, g2x; + + tmp0 = v0.x() + v1.x(); + f1x = tmp0 + v2.x(); + tmp1 = v0.x()*v0.x(); + tmp2 = tmp1 + v1.x()*tmp0; + f2x = tmp2 + v2.x()*f1x; + f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x; + g0x = f2x + v0.x()*(f1x + v0.x()); + g1x = f2x + v1.x()*(f1x + v1.x()); + g2x = f2x + v2.x()*(f1x + v2.x()); + + scalar f1y, f2y, f3y, g0y, g1y, g2y; + + tmp0 = v0.y() + v1.y(); + f1y = tmp0 + v2.y(); + tmp1 = v0.y()*v0.y(); + tmp2 = tmp1 + v1.y()*tmp0; + f2y = tmp2 + v2.y()*f1y; + f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y; + g0y = f2y + v0.y()*(f1y + v0.y()); + g1y = f2y + v1.y()*(f1y + v1.y()); + g2y = f2y + v2.y()*(f1y + v2.y()); + + scalar f1z, f2z, f3z, g0z, g1z, g2z; + + tmp0 = v0.z() + v1.z(); + f1z = tmp0 + v2.z(); + tmp1 = v0.z()*v0.z(); + tmp2 = tmp1 + v1.z()*tmp0; + f2z = tmp2 + v2.z()*f1z; + f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z; + g0z = f2z + v0.z()*(f1z + v0.z()); + g1z = f2z + v1.z()*(f1z + v1.z()); + g2z = f2z + v2.z()*(f1z + v2.z()); + + // update integrals + integrals[0] += n.x()*f1x; + integrals[1] += n.x()*f2x; + integrals[2] += n.y()*f2y; + integrals[3] += n.z()*f2z; + integrals[4] += n.x()*f3x; + integrals[5] += n.y()*f3y; + integrals[6] += n.z()*f3z; + integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x); + integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y); + integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z); } + integrals[0] *= r6; + integrals[1] *= r24; + integrals[2] *= r24; + integrals[3] *= r24; + integrals[4] *= r60; + integrals[5] *= r60; + integrals[6] *= r60; + integrals[7] *= r120; + integrals[8] *= r120; + integrals[9] *= r120; - scalar T0; - vector T1, T2, TP; + // mass + mass = integrals[0]; - compVolumeIntegrals - ( - points, - faces, - cFaces, - fNorm, - fW, + // center of mass + cM = vector(integrals[1], integrals[2], integrals[3])/mass; - T0, - T1, - T2, - TP - ); + // inertia relative to origin + J.xx() = integrals[5] + integrals[6]; + J.xy() = -integrals[7]; + J.xz() = -integrals[9]; + J.yx() = J.xy(); + J.yy() = integrals[4] + integrals[6]; + J.yz() = -integrals[8]; + J.zx() = J.xz(); + J.zy() = J.yz(); + J.zz() = integrals[4] + integrals[5]; - const scalar density = 1.0; /* assume unit density */ + // inertia relative to center of mass + J -= mass*((cM & cM)*I - cM*cM); - scalar mass = density * T0; - - /* compute center of mass */ - r = T1 / T0; - - /* compute inertia tensor */ - J.xx() = density * (T2[1] + T2[2]); - J.yy() = density * (T2[2] + T2[0]); - J.zz() = density * (T2[0] + T2[1]); - J.xy() = J.yx() = - density * TP[0]; - J.yz() = J.zy() = - density * TP[1]; - J.zx() = J.xz() = - density * TP[2]; - - ///* translate inertia tensor to center of mass */ - //J[XX] -= mass * (r[1]*r[1] + r[2]*r[2]); - //J[YY] -= mass * (r[2]*r[2] + r[0]*r[0]); - //J[ZZ] -= mass * (r[0]*r[0] + r[1]*r[1]); - //J[XY] = J[YX] += mass * r[0] * r[1]; - //J[YZ] = J[ZY] += mass * r[1] * r[2]; - //J[ZX] = J[XZ] += mass * r[2] * r[0]; + // Apply density + mass *= density; + J *= density; } +void Foam::momentOfInertia::massPropertiesShell +( + const pointField& pts, + const triFaceList& triFaces, + scalar density, + scalar& mass, + vector& cM, + tensor& J +) +{ + // Reset properties for accumulation + + mass = 0.0; + cM = vector::zero; + J = tensor::zero; + + // Find centre of mass + + forAll(triFaces, i) + { + const triFace& tri(triFaces[i]); + + triPointRef t + ( + pts[tri[0]], + pts[tri[1]], + pts[tri[2]] + ); + + scalar triMag = t.mag(); + + cM += triMag*t.centre(); + + mass += triMag; + } + + cM /= mass; + + mass *= density; + + // Find inertia around centre of mass + + forAll(triFaces, i) + { + const triFace& tri(triFaces[i]); + + J += triPointRef + ( + pts[tri[0]], + pts[tri[1]], + pts[tri[2]] + ).inertia(cM, density); + } +} + + +void Foam::momentOfInertia::massPropertiesSolid +( + const triSurface& surf, + scalar density, + scalar& mass, + vector& cM, + tensor& J +) +{ + triFaceList faces(surf.size()); + + forAll(surf, i) + { + faces[i] = triFace(surf[i]); + } + + massPropertiesSolid(surf.points(), faces, density, mass, cM, J); +} + + +void Foam::momentOfInertia::massPropertiesShell +( + const triSurface& surf, + scalar density, + scalar& mass, + vector& cM, + tensor& J +) +{ + triFaceList faces(surf.size()); + + forAll(surf, i) + { + faces[i] = triFace(surf[i]); + } + + massPropertiesShell(surf.points(), faces, density, mass, cM, J); +} + + +Foam::tensor Foam::momentOfInertia::applyParallelAxisTheorem +( + scalar mass, + const vector& cM, + const tensor& J, + const vector& refPt +) +{ + // The displacement vector (refPt = cM) is the displacement of the + // new reference point from the centre of mass of the body that + // the inertia tensor applies to. + + vector d = (refPt - cM); + + return J + mass*((d & d)*I - d*d); +} + + +Foam::tmp Foam::momentOfInertia::meshInertia +( + const polyMesh& mesh +) +{ + tmp tTf = tmp(new tensorField(mesh.nCells())); + + tensorField& tf = tTf(); + + forAll(tf, cI) + { + tf[cI] = meshInertia(mesh, cI); + } + + return tTf; +} + + +Foam::tensor Foam::momentOfInertia::meshInertia +( + const polyMesh& mesh, + label cellI +) +{ + List cellTets = polyMeshTetDecomposition::cellTetIndices + ( + mesh, + cellI + ); + + triFaceList faces(cellTets.size()); + + forAll(cellTets, cTI) + { + faces[cTI] = cellTets[cTI].faceTriIs(mesh); + } + + scalar m = 0.0; + vector cM = vector::zero; + tensor J = tensor::zero; + + massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J); + + return J; +} + // ************************************************************************* // diff --git a/src/meshTools/momentOfInertia/momentOfInertia.H b/src/meshTools/momentOfInertia/momentOfInertia.H index 8ff6209f0b..bb5ffef6a8 100644 --- a/src/meshTools/momentOfInertia/momentOfInertia.H +++ b/src/meshTools/momentOfInertia/momentOfInertia.H @@ -25,6 +25,9 @@ Class momentOfInertia Description + Calculates the inertia tensor and principal axes and moments of a + polyhedra/cells/triSurfaces. Inertia can either be of the solid body or + of a thin shell. SourceFiles momentOfInertia.H @@ -34,34 +37,86 @@ SourceFiles #ifndef momentOfInertia_H #define momentOfInertia_H -#include "tensor.H" -#include "primitiveMesh.H" +#include "tetPointRef.H" +#include "triFaceList.H" +#include "triSurface.H" +#include "polyMesh.H" +#include "polyMeshTetDecomposition.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { -////- Moment of inertia around cell centre for single cell. -//tensor momentOfInertia -//( -// const pointField&, -// const faceList&, -// const cell&, -// const point& cc -//); +/*---------------------------------------------------------------------------*\ + Class momentOfInertia Declaration +\*---------------------------------------------------------------------------*/ + +class momentOfInertia +{ + +public: + + static void massPropertiesSolid + ( + const pointField& pts, + const triFaceList& triFaces, + scalar density, + scalar& mass, + vector& cM, + tensor& J + ); + + static void massPropertiesShell + ( + const pointField& pts, + const triFaceList& triFaces, + scalar density, + scalar& mass, + vector& cM, + tensor& J + ); + + static void massPropertiesSolid + ( + const triSurface& surf, + scalar density, + scalar& mass, + vector& cM, + tensor& J + ); + + static void massPropertiesShell + ( + const triSurface& surf, + scalar density, + scalar& mass, + vector& cM, + tensor& J + ); + + static tensor applyParallelAxisTheorem + ( + scalar mass, + const vector& cM, + const tensor& J, + const vector& refPt + ); + + // Calculate the inertia tensor for all cells in the mesh + static tmp meshInertia + ( + const polyMesh& mesh + ); + + // Calculate the inertia tensor the given cell + static tensor meshInertia + ( + const polyMesh& mesh, + label cellI + ); +}; -// Calculate -// - centre of mass -// - inertia tensor around (0,0,0) -void momentOfIntertia -( - const pointField&, - const faceList&, - const cell&, - point& r, - tensor& Jorigin -); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //