mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Merge commit 'OpenCFD/master' into olesenm
This commit is contained in:
Mark Olesen
@ -27,12 +27,14 @@ Application
|
||||
|
||||
Description
|
||||
Calculates the inertia tensor and principal axes and moments of a
|
||||
test face.
|
||||
test face and tetrahedron.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "ListOps.H"
|
||||
#include "face.H"
|
||||
#include "tetPointRef.H"
|
||||
#include "triFaceList.H"
|
||||
#include "OFstream.H"
|
||||
#include "meshTools.H"
|
||||
|
||||
@ -40,67 +42,310 @@ Description
|
||||
|
||||
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[])
|
||||
{
|
||||
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);
|
||||
|
||||
scalar density = 1.0;
|
||||
|
||||
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("momentOfInertiaTest.obj");
|
||||
|
||||
Info<< nl << "Writing test face and scaled principal axes to "
|
||||
<< str.name() << endl;
|
||||
|
||||
forAll(pts, ptI)
|
||||
{
|
||||
meshTools::writeOBJ(str, pts[ptI]);
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
str << "l";
|
||||
|
||||
forAll(f, fI)
|
||||
{
|
||||
str << ' ' << fI + 1;
|
||||
}
|
||||
label nPts = 4;
|
||||
|
||||
str << " 1" << endl;
|
||||
pointField pts(nPts);
|
||||
|
||||
scalar scale = mag(Cf - pts[f[0]])/eVal.component(findMin(eVal));
|
||||
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);
|
||||
|
||||
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());
|
||||
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;
|
||||
}
|
||||
|
||||
for (label i = nPts + 1; i < nPts + 4; i++)
|
||||
{
|
||||
str << "l " << nPts + 1 << ' ' << i + 1 << endl;
|
||||
}
|
||||
|
||||
Info<< nl << "End" << nl << endl;
|
||||
|
||||
3
applications/utilities/surface/surfaceInertia/Make/files
Normal file
3
applications/utilities/surface/surfaceInertia/Make/files
Normal file
@ -0,0 +1,3 @@
|
||||
surfaceInertia.C
|
||||
|
||||
EXE = $(FOAM_APPBIN)/surfaceInertia
|
||||
@ -0,0 +1,7 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/meshTools/lnInclude \
|
||||
-I$(LIB_SRC)/triSurface/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lmeshTools \
|
||||
-ltriSurface
|
||||
296
applications/utilities/surface/surfaceInertia/block.obj
Normal file
296
applications/utilities/surface/surfaceInertia/block.obj
Normal file
@ -0,0 +1,296 @@
|
||||
# Wavefront OBJ file
|
||||
# Regions:
|
||||
# 0 movingWall
|
||||
# 1 fixedWalls
|
||||
# 2 frontAndBack
|
||||
#
|
||||
# points : 96
|
||||
# triangles : 188
|
||||
#
|
||||
v 0.1 0.1 0.5
|
||||
v 0.1 0.1 0.5625
|
||||
v 0.16 0.1 0.5625
|
||||
v 0.16 0.1 0.5
|
||||
v 0.1 0.1 0.625
|
||||
v 0.16 0.1 0.625
|
||||
v 0.1 0.1 0.6875
|
||||
v 0.16 0.1 0.6875
|
||||
v 0.1 0.1 0.75
|
||||
v 0.16 0.1 0.75
|
||||
v 0.22 0.1 0.5625
|
||||
v 0.22 0.1 0.5
|
||||
v 0.22 0.1 0.625
|
||||
v 0.22 0.1 0.6875
|
||||
v 0.22 0.1 0.75
|
||||
v 0.28 0.1 0.5625
|
||||
v 0.28 0.1 0.5
|
||||
v 0.28 0.1 0.625
|
||||
v 0.28 0.1 0.6875
|
||||
v 0.28 0.1 0.75
|
||||
v 0.34 0.1 0.5625
|
||||
v 0.34 0.1 0.5
|
||||
v 0.34 0.1 0.625
|
||||
v 0.34 0.1 0.6875
|
||||
v 0.34 0.1 0.75
|
||||
v 0.4 0.1 0.5625
|
||||
v 0.4 0.1 0.5
|
||||
v 0.4 0.1 0.625
|
||||
v 0.4 0.1 0.6875
|
||||
v 0.4 0.1 0.75
|
||||
v 0.1 -0.3 0.5
|
||||
v 0.1 -0.3 0.5625
|
||||
v 0.1 -0.166667 0.5625
|
||||
v 0.1 -0.166667 0.5
|
||||
v 0.1 -0.0333333 0.5625
|
||||
v 0.1 -0.0333333 0.5
|
||||
v 0.1 -0.3 0.625
|
||||
v 0.1 -0.166667 0.625
|
||||
v 0.1 -0.0333333 0.625
|
||||
v 0.1 -0.3 0.6875
|
||||
v 0.1 -0.166667 0.6875
|
||||
v 0.1 -0.0333333 0.6875
|
||||
v 0.1 -0.3 0.75
|
||||
v 0.1 -0.166667 0.75
|
||||
v 0.1 -0.0333333 0.75
|
||||
v 0.4 -0.3 0.5
|
||||
v 0.4 -0.166667 0.5
|
||||
v 0.4 -0.166667 0.5625
|
||||
v 0.4 -0.3 0.5625
|
||||
v 0.4 -0.0333333 0.5
|
||||
v 0.4 -0.0333333 0.5625
|
||||
v 0.4 -0.166667 0.625
|
||||
v 0.4 -0.3 0.625
|
||||
v 0.4 -0.0333333 0.625
|
||||
v 0.4 -0.166667 0.6875
|
||||
v 0.4 -0.3 0.6875
|
||||
v 0.4 -0.0333333 0.6875
|
||||
v 0.4 -0.166667 0.75
|
||||
v 0.4 -0.3 0.75
|
||||
v 0.4 -0.0333333 0.75
|
||||
v 0.16 -0.3 0.5
|
||||
v 0.16 -0.3 0.5625
|
||||
v 0.16 -0.3 0.625
|
||||
v 0.16 -0.3 0.6875
|
||||
v 0.16 -0.3 0.75
|
||||
v 0.22 -0.3 0.5
|
||||
v 0.22 -0.3 0.5625
|
||||
v 0.22 -0.3 0.625
|
||||
v 0.22 -0.3 0.6875
|
||||
v 0.22 -0.3 0.75
|
||||
v 0.28 -0.3 0.5
|
||||
v 0.28 -0.3 0.5625
|
||||
v 0.28 -0.3 0.625
|
||||
v 0.28 -0.3 0.6875
|
||||
v 0.28 -0.3 0.75
|
||||
v 0.34 -0.3 0.5
|
||||
v 0.34 -0.3 0.5625
|
||||
v 0.34 -0.3 0.625
|
||||
v 0.34 -0.3 0.6875
|
||||
v 0.34 -0.3 0.75
|
||||
v 0.16 -0.166667 0.5
|
||||
v 0.16 -0.0333333 0.5
|
||||
v 0.22 -0.166667 0.5
|
||||
v 0.22 -0.0333333 0.5
|
||||
v 0.28 -0.166667 0.5
|
||||
v 0.28 -0.0333333 0.5
|
||||
v 0.34 -0.166667 0.5
|
||||
v 0.34 -0.0333333 0.5
|
||||
v 0.16 -0.166667 0.75
|
||||
v 0.16 -0.0333333 0.75
|
||||
v 0.22 -0.166667 0.75
|
||||
v 0.22 -0.0333333 0.75
|
||||
v 0.28 -0.166667 0.75
|
||||
v 0.28 -0.0333333 0.75
|
||||
v 0.34 -0.166667 0.75
|
||||
v 0.34 -0.0333333 0.75
|
||||
g movingWall
|
||||
f 1 2 3
|
||||
f 3 4 1
|
||||
f 2 5 6
|
||||
f 6 3 2
|
||||
f 5 7 8
|
||||
f 8 6 5
|
||||
f 7 9 10
|
||||
f 10 8 7
|
||||
f 4 3 11
|
||||
f 11 12 4
|
||||
f 3 6 13
|
||||
f 13 11 3
|
||||
f 6 8 14
|
||||
f 14 13 6
|
||||
f 8 10 15
|
||||
f 15 14 8
|
||||
f 12 11 16
|
||||
f 16 17 12
|
||||
f 11 13 18
|
||||
f 18 16 11
|
||||
f 13 14 19
|
||||
f 19 18 13
|
||||
f 14 15 20
|
||||
f 20 19 14
|
||||
f 17 16 21
|
||||
f 21 22 17
|
||||
f 16 18 23
|
||||
f 23 21 16
|
||||
f 18 19 24
|
||||
f 24 23 18
|
||||
f 19 20 25
|
||||
f 25 24 19
|
||||
f 22 21 26
|
||||
f 26 27 22
|
||||
f 21 23 28
|
||||
f 28 26 21
|
||||
f 23 24 29
|
||||
f 29 28 23
|
||||
f 24 25 30
|
||||
f 30 29 24
|
||||
g fixedWalls
|
||||
f 31 32 33
|
||||
f 33 34 31
|
||||
f 34 33 35
|
||||
f 35 36 34
|
||||
f 36 35 2
|
||||
f 2 1 36
|
||||
f 32 37 38
|
||||
f 38 33 32
|
||||
f 33 38 39
|
||||
f 39 35 33
|
||||
f 35 39 5
|
||||
f 5 2 35
|
||||
f 37 40 41
|
||||
f 41 38 37
|
||||
f 38 41 42
|
||||
f 42 39 38
|
||||
f 39 42 7
|
||||
f 7 5 39
|
||||
f 40 43 44
|
||||
f 44 41 40
|
||||
f 41 44 45
|
||||
f 45 42 41
|
||||
f 42 45 9
|
||||
f 9 7 42
|
||||
f 46 47 48
|
||||
f 48 49 46
|
||||
f 47 50 51
|
||||
f 51 48 47
|
||||
f 50 27 26
|
||||
f 26 51 50
|
||||
f 49 48 52
|
||||
f 52 53 49
|
||||
f 48 51 54
|
||||
f 54 52 48
|
||||
f 51 26 28
|
||||
f 28 54 51
|
||||
f 53 52 55
|
||||
f 55 56 53
|
||||
f 52 54 57
|
||||
f 57 55 52
|
||||
f 54 28 29
|
||||
f 29 57 54
|
||||
f 56 55 58
|
||||
f 58 59 56
|
||||
f 55 57 60
|
||||
f 60 58 55
|
||||
f 57 29 30
|
||||
f 30 60 57
|
||||
f 31 61 62
|
||||
f 62 32 31
|
||||
f 32 62 63
|
||||
f 63 37 32
|
||||
f 37 63 64
|
||||
f 64 40 37
|
||||
f 40 64 65
|
||||
f 65 43 40
|
||||
f 61 66 67
|
||||
f 67 62 61
|
||||
f 62 67 68
|
||||
f 68 63 62
|
||||
f 63 68 69
|
||||
f 69 64 63
|
||||
f 64 69 70
|
||||
f 70 65 64
|
||||
f 66 71 72
|
||||
f 72 67 66
|
||||
f 67 72 73
|
||||
f 73 68 67
|
||||
f 68 73 74
|
||||
f 74 69 68
|
||||
f 69 74 75
|
||||
f 75 70 69
|
||||
f 71 76 77
|
||||
f 77 72 71
|
||||
f 72 77 78
|
||||
f 78 73 72
|
||||
f 73 78 79
|
||||
f 79 74 73
|
||||
f 74 79 80
|
||||
f 80 75 74
|
||||
f 76 46 49
|
||||
f 49 77 76
|
||||
f 77 49 53
|
||||
f 53 78 77
|
||||
f 78 53 56
|
||||
f 56 79 78
|
||||
f 79 56 59
|
||||
f 59 80 79
|
||||
g frontAndBack
|
||||
f 31 34 81
|
||||
f 81 61 31
|
||||
f 34 36 82
|
||||
f 82 81 34
|
||||
f 36 1 4
|
||||
f 4 82 36
|
||||
f 61 81 83
|
||||
f 83 66 61
|
||||
f 81 82 84
|
||||
f 84 83 81
|
||||
f 82 4 12
|
||||
f 12 84 82
|
||||
f 66 83 85
|
||||
f 85 71 66
|
||||
f 83 84 86
|
||||
f 86 85 83
|
||||
f 84 12 17
|
||||
f 17 86 84
|
||||
f 71 85 87
|
||||
f 87 76 71
|
||||
f 85 86 88
|
||||
f 88 87 85
|
||||
f 86 17 22
|
||||
f 22 88 86
|
||||
f 76 87 47
|
||||
f 47 46 76
|
||||
f 87 88 50
|
||||
f 50 47 87
|
||||
f 88 22 27
|
||||
f 27 50 88
|
||||
f 43 65 89
|
||||
f 89 44 43
|
||||
f 44 89 90
|
||||
f 90 45 44
|
||||
f 45 90 10
|
||||
f 10 9 45
|
||||
f 65 70 91
|
||||
f 91 89 65
|
||||
f 89 91 92
|
||||
f 92 90 89
|
||||
f 90 92 15
|
||||
f 15 10 90
|
||||
f 70 75 93
|
||||
f 93 91 70
|
||||
f 91 93 94
|
||||
f 94 92 91
|
||||
f 92 94 20
|
||||
f 20 15 92
|
||||
f 75 80 95
|
||||
f 95 93 75
|
||||
f 93 95 96
|
||||
f 96 94 93
|
||||
f 94 96 25
|
||||
f 25 20 94
|
||||
f 80 59 58
|
||||
f 58 95 80
|
||||
f 95 58 60
|
||||
f 60 96 95
|
||||
f 96 60 30
|
||||
f 30 25 96
|
||||
12608
applications/utilities/surface/surfaceInertia/cylinder.obj
Normal file
12608
applications/utilities/surface/surfaceInertia/cylinder.obj
Normal file
File diff suppressed because it is too large
Load Diff
13174
applications/utilities/surface/surfaceInertia/dangermouse.obj
Normal file
13174
applications/utilities/surface/surfaceInertia/dangermouse.obj
Normal file
File diff suppressed because it is too large
Load Diff
261982
applications/utilities/surface/surfaceInertia/sphere.obj
Normal file
261982
applications/utilities/surface/surfaceInertia/sphere.obj
Normal file
File diff suppressed because it is too large
Load Diff
422
applications/utilities/surface/surfaceInertia/surfaceInertia.C
Normal file
422
applications/utilities/surface/surfaceInertia/surfaceInertia.C
Normal file
@ -0,0 +1,422 @@
|
||||
/*---------------------------------------------------------------------------*\
|
||||
========= |
|
||||
\\ / 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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
Reference in New Issue
Block a user