/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2009-2010 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 . Application momentOfInertiaTest Description Calculates the inertia tensor and principal axes and moments of a test face and tetrahedron. \*---------------------------------------------------------------------------*/ #include "ListOps.H" #include "face.H" #include "tetPointRef.H" #include "triFaceList.H" #include "OFstream.H" #include "meshTools.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[]) { scalar density = 1.0; { label nPts = 6; pointField pts(nPts); pts[0] = point(4.495, 3.717, -4.112); pts[1] = point(4.421, 3.932, -4.112); pts[2] = point(4.379, 4.053, -4.112); pts[3] = point(4.301, 4.026, -4.300); pts[4] = point(4.294, 4.024, -4.317); pts[5] = point(4.409, 3.687, -4.317); face f(identity(nPts)); point Cf = f.centre(pts); tensor J = tensor::zero; J = f.inertia(pts, Cf, density); vector eVal = eigenValues(J); tensor eVec = eigenVectors(J); Info<< nl << "Inertia tensor of test face " << J << nl << "eigenValues (principal moments) " << eVal << nl << "eigenVectors (principal axes) " << eVec << endl; OFstream str("momentOfInertiaTestFace.obj"); Info<< nl << "Writing test face and scaled principal axes to " << str.name() << endl; forAll(pts, ptI) { meshTools::writeOBJ(str, pts[ptI]); } str << "l"; forAll(f, fI) { str << ' ' << fI + 1; } str << " 1" << endl; scalar scale = mag(Cf - pts[f[0]])/eVal.component(findMin(eVal)); meshTools::writeOBJ(str, Cf); meshTools::writeOBJ(str, Cf + scale*eVal.x()*eVec.x()); meshTools::writeOBJ(str, Cf + scale*eVal.y()*eVec.y()); meshTools::writeOBJ(str, Cf + scale*eVal.z()*eVec.z()); for (label i = nPts + 1; i < nPts + 4; i++) { str << "l " << nPts + 1 << ' ' << i + 1 << endl; } } { label nPts = 4; pointField pts(nPts); pts[0] = point(0, 0, 0); pts[1] = point(1, 0, 0); pts[2] = point(0.5, 1, 0); pts[3] = point(0.5, 0.5, 1); tetPointRef tet(pts[0], pts[1], pts[2], pts[3]); triFaceList tetFaces(4); tetFaces[0] = triFace(0, 2, 1); tetFaces[1] = triFace(1, 2, 3); tetFaces[2] = triFace(0, 3, 2); tetFaces[3] = triFace(0, 1, 3); scalar m = 0.0; vector cM = vector::zero; tensor J = tensor::zero; massPropertiesSolid ( pts, tetFaces, density, m, cM, J ); vector eVal = eigenValues(J); tensor eVec = eigenVectors(J); Info<< nl << "Mass of tetrahedron " << m << nl << "Centre of mass of tetrahedron " << cM << nl << "Inertia tensor of tetrahedron " << J << nl << "eigenValues (principal moments) " << eVal << nl << "eigenVectors (principal axes) " << eVec << endl; OFstream str("momentOfInertiaTestTet.obj"); Info<< nl << "Writing test tetrahedron and scaled principal axes to " << str.name() << endl; forAll(pts, ptI) { meshTools::writeOBJ(str, pts[ptI]); } forAll(tetFaces, tFI) { const triFace& f = tetFaces[tFI]; str << "l"; forAll(f, fI) { str << ' ' << f[fI] + 1; } str << ' ' << f[0] + 1 << endl; } scalar scale = mag(cM - pts[0])/eVal.component(findMin(eVal)); meshTools::writeOBJ(str, cM); meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x()); meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y()); meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z()); for (label i = nPts + 1; i < nPts + 4; i++) { str << "l " << nPts + 1 << ' ' << i + 1 << endl; } } Info<< nl << "End" << nl << endl; return 0; } // ************************************************************************* //