Skip to content
GitLab
Commit aa15f447 authored by graham's avatar graham
Browse files

ENH: Transferring momentOfInertia calc from utils to meshTools lib. 

Adding mesh cell inertia calc.
parent cde4e666
Hide whitespace changes
Inline Side-by-side
applications/test/momentOfInertia/Make/options
View file @ aa15f447
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude
EXE_LIBS = \
-lmeshTools
applications/test/momentOfInertia/Test-momentOfInertia.C
View file @ aa15f447
......@@ -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
)
int main(int argc, char *argv[])
{
// 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;
}
argList::addOption
(
"cell",
"label",
"cell to use for inertia calculation, defaults to 0"
);
#include "setRootCase.H"
#include "createTime.H"
#include "createPolyMesh.H"
int main(int argc, char *argv[])
{
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;
......
applications/utilities/surface/surfaceInertia/surfaceInertia.C
View file @ aa15f447
......@@ -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;
}
......
src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
View file @ aa15f447
......@@ -41,6 +41,7 @@ SourceFiles
#include "tetPointRef.H"
#include "triPointRef.H"