diff --git a/applications/test/momentOfInertia/Make/options b/applications/test/momentOfInertia/Make/options
index 54c035b8f55d183c1ad02bc372398feceaf31718..0b32f3355b3b4bf0fda452da0efb5e5c9e2201fa 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 9a8ae40a5aaa31441ab819518df022157d246182..22327f1bddfa6b18092485d7f87d8fbd589c67af 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
-)
+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;
diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
index f72bcb6e5ebc2288038a6670e98a605ab554c401..2b8183295f77bdbe6811b612db7c8e7bc06e17cd 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 3e5ee425d6ed9a527dbb067732d07e7b809a9174..3a7513770d756d4b7e9817ec1352d6945c3c7b82 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 6b8a7871a04c6fbadb9ee3d61bf47215e8b5764c..9a4287b5c4226e20f6f3e155dd4b2cbc306029f9 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 2c9f81195d082b75e3b08a62ede15f2a8ec37d85..47310b27e3ae0f4470c8c73848f2eb1ed47e1734 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -127,6 +127,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 af3265a5ecf9f95f53e3ac404fddea1a9602b3df..8dc1a53d61f4d0d26bcae4615686bc857327691e 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 <http://www.gnu.org/licenses/>.
-Class
- momentOfInertia
+\*---------------------------------------------------------------------------*/
-Description
- Reimplementation of volInt.c by Brian Mirtich.
- * mirtich@cs.berkeley.edu *
- * http://www.cs.berkeley.edu/~mirtich *
+#include "momentOfInertia.H"
--------------------------------------------------------------------------------
-*/
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
-#include "momentOfInertia.H"
-//#include "pyramidPointFaceRef.H"
-
-// * * * * * * * * * * * * * Private 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
+ const pointField& pts,
+ const triFaceList& triFaces,
+ scalar density,
+ scalar& mass,
+ vector& cM,
+ tensor& J
)
{
- P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;
-
- forAll(f, i)
+ // 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)
{
- 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;
+ 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 Foam::compFaceIntegrals
+void Foam::momentOfInertia::massPropertiesShell
(
- 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
+ const pointField& pts,
+ const triFaceList& triFaces,
+ scalar density,
+ scalar& mass,
+ vector& cM,
+ tensor& J
)
{
- scalar P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb;
+ // Reset properties for accumulation
- compProjectionIntegrals
- (
- points,
- f,
- A,
- B,
-
- P1,
- Pa,
- Pb,
- Paa,
- Pab,
- Pbb,
- Paaa,
- Paab,
- Pabb,
- Pbbb
- );
+ mass = 0.0;
+ cM = vector::zero;
+ J = tensor::zero;
- 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));
-}
+ // Find centre of mass
+ forAll(triFaces, i)
+ {
+ const triFace& tri(triFaces[i]);
-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;
+ triPointRef t
+ (
+ pts[tri[0]],
+ pts[tri[1]],
+ pts[tri[2]]
+ );
- forAll(cFaces, i)
- {
- const vector& n = fNorm[i];
+ scalar triMag = t.mag();
- scalar nx = mag(n[0]);
- scalar ny = mag(n[1]);
- scalar nz = mag(n[2]);
+ cM += triMag*t.centre();
- direction A, B, C;
+ mass += triMag;
+ }
- if (nx > ny && nx > nz)
- {
- C = 0;
- }
- else
- {
- C = (ny > nz) ? 1 : 2;
- }
+ cM /= mass;
- A = (C + 1) % 3;
- B = (A + 1) % 3;
+ mass *= density;
- 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
- );
+ // Find inertia around centre of mass
- T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc));
+ forAll(triFaces, i)
+ {
+ const triFace& tri(triFaces[i]);
+
+ J += triPointRef
+ (
+ pts[tri[0]],
+ pts[tri[1]],
+ pts[tri[2]]
+ ).inertia(cM, density);
+ }
+}
- 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;
+void Foam::momentOfInertia::massPropertiesSolid
+(
+ const triSurface& surf,
+ scalar density,
+ scalar& mass,
+ vector& cM,
+ tensor& J
+)
+{
+ triFaceList faces(surf.size());
- TP[A] += n[A] * Faab;
- TP[B] += n[B] * Fbbc;
- TP[C] += n[C] * Fcca;
+ forAll(surf, i)
+ {
+ faces[i] = triFace(surf[i]);
}
- T1 /= 2;
- T2 /= 3;
- TP /= 2;
+ massPropertiesSolid(surf.points(), faces, density, mass, cM, J);
}
-// Calculate
-// - r: centre of mass
-// - J: inertia around origin (point 0,0,0)
-void Foam::momentOfIntertia
+void Foam::momentOfInertia::massPropertiesShell
(
- const pointField& points,
- const faceList& faces,
- const cell& cFaces,
- point& r,
+ const triSurface& surf,
+ scalar density,
+ scalar& mass,
+ vector& cM,
tensor& J
)
{
- // Face normals
- vectorField fNorm(cFaces.size());
- scalarField fW(cFaces.size());
+ triFaceList faces(surf.size());
- forAll(cFaces, i)
+ forAll(surf, i)
{
- label faceI = cFaces[i];
+ faces[i] = triFace(surf[i]);
+ }
+
+ massPropertiesShell(surf.points(), faces, density, mass, cM, J);
+}
- const face& f = faces[faceI];
- fNorm[i] = f.normal(points);
- fNorm[i] /= mag(fNorm[i]) + VSMALL;
+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);
+}
+
- fW[i] = - (fNorm[i] & points[f[0]]);
+Foam::tmp<Foam::tensorField> Foam::momentOfInertia::meshInertia
+(
+ const polyMesh& mesh
+)
+{
+ tmp<tensorField> tTf = tmp<tensorField>(new tensorField(mesh.nCells()));
+
+ tensorField& tf = tTf();
+
+ forAll(tf, cI)
+ {
+ tf[cI] = meshInertia(mesh, cI);
}
+ return tTf;
+}
- scalar T0;
- vector T1, T2, TP;
- compVolumeIntegrals
+Foam::tensor Foam::momentOfInertia::meshInertia
+(
+ const polyMesh& mesh,
+ label cellI
+)
+{
+ List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
(
- points,
- faces,
- cFaces,
- fNorm,
- fW,
-
- T0,
- T1,
- T2,
- TP
+ mesh,
+ cellI
);
- const scalar density = 1.0; /* assume unit density */
+ triFaceList faces(cellTets.size());
- scalar mass = density * T0;
+ forAll(cellTets, cTI)
+ {
+ faces[cTI] = cellTets[cTI].faceTriIs(mesh);
+ }
- /* compute center of mass */
- r = T1 / T0;
+ scalar m = 0.0;
+ vector cM = vector::zero;
+ tensor J = tensor::zero;
- /* 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];
+ massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J);
- ///* 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];
+ return J;
}
-
// ************************************************************************* //
diff --git a/src/meshTools/momentOfInertia/momentOfInertia.H b/src/meshTools/momentOfInertia/momentOfInertia.H
index 8ff6209f0b18adc821b0657108b0c86717c6e675..bb5ffef6a8fb0968ae8ef954eba515a66bf57d4f 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
-//);
-
-// Calculate
-// - centre of mass
-// - inertia tensor around (0,0,0)
-void momentOfIntertia
-(
- const pointField&,
- const faceList&,
- const cell&,
- point& r,
- tensor& Jorigin
-);
+/*---------------------------------------------------------------------------*\
+ 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<tensorField> meshInertia
+ (
+ const polyMesh& mesh
+ );
+
+ // Calculate the inertia tensor the given cell
+ static tensor meshInertia
+ (
+ const polyMesh& mesh,
+ label cellI
+ );
+};
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //