aa847562ce
http://www.geometrictools.com/LibPhysics/RigidBody/RigidBody.html in surfaceInertia utility. Allows calculation of shell or solid properties of a triSurface, with optionally specified volumetric or surface density. Default behaviour is to calculate inertia around centre of mass, with option to specify reference point. Various test surfaces supplied and tested against analytical results.
423 lines
11 KiB
C
423 lines
11 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | Copyright (C) 2009-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 2 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, write to the Free Software Foundation,
|
|
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
|
|
|
|
Application
|
|
momentOfInertiaTest
|
|
|
|
Description
|
|
Calculates the inertia tensor and principal axes and moments of a
|
|
command line specified triSurface. Inertia can either be of the
|
|
solid body or of a thin shell.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "argList.H"
|
|
#include "ListOps.H"
|
|
#include "face.H"
|
|
#include "tetPointRef.H"
|
|
#include "triFaceList.H"
|
|
#include "triSurface.H"
|
|
#include "OFstream.H"
|
|
#include "meshTools.H"
|
|
#include "Random.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
|
|
(
|
|
"Calculates the inertia tensor and principal axes and moments "
|
|
"of a command line specified triSurface. Inertia can either "
|
|
"be of the solid body or of a thin shell."
|
|
);
|
|
|
|
argList::noParallel();
|
|
|
|
argList::validArgs.append("surface file");
|
|
|
|
argList::addBoolOption("shellProperties");
|
|
|
|
argList::addOption
|
|
(
|
|
"density",
|
|
"scalar",
|
|
"Specify density,"
|
|
"kg/m3 for solid properties, kg/m2 for shell properties"
|
|
);
|
|
|
|
argList::addOption
|
|
(
|
|
"referencePoint",
|
|
"vector",
|
|
"Inertia relative to this point, not the centre of mass"
|
|
);
|
|
|
|
argList args(argc, argv);
|
|
|
|
fileName surfFileName(args.additionalArgs()[0]);
|
|
|
|
triSurface surf(surfFileName);
|
|
|
|
scalar density = args.optionLookupOrDefault("density", 1.0);
|
|
|
|
vector refPt = vector::zero;
|
|
|
|
bool calcAroundRefPt = args.optionReadIfPresent("referencePoint", refPt);
|
|
|
|
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
|
|
);
|
|
}
|
|
else
|
|
{
|
|
massPropertiesSolid
|
|
(
|
|
surf.points(),
|
|
faces,
|
|
density,
|
|
m,
|
|
cM,
|
|
J
|
|
);
|
|
}
|
|
|
|
vector eVal = eigenValues(J);
|
|
|
|
tensor eVec = eigenVectors(J);
|
|
|
|
label pertI = 0;
|
|
|
|
Random rand(57373);
|
|
|
|
while ((magSqr(eVal) < VSMALL) && pertI < 10)
|
|
{
|
|
WarningIn(args.executable() + "::main")
|
|
<< "No eigenValues found, shape may have symmetry, "
|
|
<< "perturbing inertia tensor diagonal" << endl;
|
|
|
|
J.xx() *= 1.0 + SMALL*rand.scalar01();
|
|
J.yy() *= 1.0 + SMALL*rand.scalar01();
|
|
J.zz() *= 1.0 + SMALL*rand.scalar01();
|
|
|
|
eVal = eigenValues(J);
|
|
|
|
eVec = eigenVectors(J);
|
|
|
|
pertI++;
|
|
}
|
|
|
|
Info<< nl
|
|
<< "Density = " << density << nl
|
|
<< "Mass = " << m << nl
|
|
<< "Centre of mass = " << cM << nl
|
|
<< "Inertia tensor around centre of mass = " << J << nl
|
|
<< "eigenValues (principal moments) = " << eVal << nl
|
|
<< "eigenVectors (principal axes) = "
|
|
<< eVec.x() << ' ' << eVec.y() << ' ' << eVec.z()
|
|
<< endl;
|
|
|
|
if(calcAroundRefPt)
|
|
{
|
|
Info << "Inertia tensor relative to " << refPt << " = "
|
|
<< applyParallelAxisTheorem(m, cM, J, refPt)
|
|
<< endl;
|
|
}
|
|
|
|
OFstream str("axes.obj");
|
|
|
|
Info<< nl << "Writing scaled principal axes at centre of mass of "
|
|
<< surfFileName << " to " << str.name() << endl;
|
|
|
|
scalar scale = mag(cM - surf.points()[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 = 1; i < 4; i++)
|
|
{
|
|
str << "l " << 1 << ' ' << i + 1 << endl;
|
|
}
|
|
|
|
Info<< nl << "End" << nl << endl;
|
|
|
|
return 0;
|
|
}
|
|
|
|
|
|
// ************************************************************************* //
|