1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
11 OpenFOAM is free software: you can redistribute it and/or modify it
12 under the terms of the GNU General Public License as published by
13 the Free Software Foundation, either version 3 of the License, or
14 (at your option) any later version.
16 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
28 Calculates the inertia tensor and principal axes and moments of a
29 test face, tetrahedron and mesh.
31 \*---------------------------------------------------------------------------*/
38 #include "tetPointRef.H"
39 #include "triFaceList.H"
41 #include "meshTools.H"
42 #include "momentOfInertia.H"
44 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
48 int main(int argc, char *argv[])
54 "cell to use for inertia calculation, defaults to 0"
57 #include "setRootCase.H"
58 #include "createTime.H"
59 #include "createPolyMesh.H"
68 pts[0] = point(4.495, 3.717, -4.112);
69 pts[1] = point(4.421, 3.932, -4.112);
70 pts[2] = point(4.379, 4.053, -4.112);
71 pts[3] = point(4.301, 4.026, -4.300);
72 pts[4] = point(4.294, 4.024, -4.317);
73 pts[5] = point(4.409, 3.687, -4.317);
75 face f(identity(nPts));
77 point Cf = f.centre(pts);
79 tensor J = tensor::zero;
81 J = f.inertia(pts, Cf, density);
83 vector eVal = eigenValues(J);
85 tensor eVec = eigenVectors(J);
87 Info<< nl << "Inertia tensor of test face " << J << nl
88 << "eigenValues (principal moments) " << eVal << nl
89 << "eigenVectors (principal axes) " << eVec
92 OFstream str("momentOfInertiaTestFace.obj");
94 Info<< nl << "Writing test face and scaled principal axes to "
95 << str.name() << endl;
99 meshTools::writeOBJ(str, pts[ptI]);
106 str << ' ' << fI + 1;
111 scalar scale = mag(Cf - pts[f[0]])/eVal.component(findMin(eVal));
113 meshTools::writeOBJ(str, Cf);
114 meshTools::writeOBJ(str, Cf + scale*eVal.x()*eVec.x());
115 meshTools::writeOBJ(str, Cf + scale*eVal.y()*eVec.y());
116 meshTools::writeOBJ(str, Cf + scale*eVal.z()*eVec.z());
118 for (label i = nPts + 1; i < nPts + 4; i++)
120 str << "l " << nPts + 1 << ' ' << i + 1 << endl;
127 pointField pts(nPts);
129 pts[0] = point(0, 0, 0);
130 pts[1] = point(1, 0, 0);
131 pts[2] = point(0.5, 1, 0);
132 pts[3] = point(0.5, 0.5, 1);
134 tetPointRef tet(pts[0], pts[1], pts[2], pts[3]);
136 triFaceList tetFaces(4);
138 tetFaces[0] = triFace(0, 2, 1);
139 tetFaces[1] = triFace(1, 2, 3);
140 tetFaces[2] = triFace(0, 3, 2);
141 tetFaces[3] = triFace(0, 1, 3);
144 vector cM = vector::zero;
145 tensor J = tensor::zero;
147 momentOfInertia::massPropertiesSolid(pts, tetFaces, density, m, cM, J);
149 vector eVal = eigenValues(J);
151 tensor eVec = eigenVectors(J);
154 << "Mass of tetrahedron " << m << nl
155 << "Centre of mass of tetrahedron " << cM << nl
156 << "Inertia tensor of tetrahedron " << J << nl
157 << "eigenValues (principal moments) " << eVal << nl
158 << "eigenVectors (principal axes) " << eVec
161 OFstream str("momentOfInertiaTestTet.obj");
163 Info<< nl << "Writing test tetrahedron and scaled principal axes to "
164 << str.name() << endl;
168 meshTools::writeOBJ(str, pts[ptI]);
171 forAll(tetFaces, tFI)
173 const triFace& f = tetFaces[tFI];
179 str << ' ' << f[fI] + 1;
182 str << ' ' << f[0] + 1 << endl;
185 scalar scale = mag(cM - pts[0])/eVal.component(findMin(eVal));
187 meshTools::writeOBJ(str, cM);
188 meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
189 meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
190 meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
192 for (label i = nPts + 1; i < nPts + 4; i++)
194 str << "l " << nPts + 1 << ' ' << i + 1 << endl;
199 const label cellI = args.optionLookupOrDefault("cell", 0);
201 tensorField mI(momentOfInertia::meshInertia(mesh));
203 tensor& J = mI[cellI];
205 vector eVal = eigenValues(J);
208 << "Inertia tensor of cell " << cellI << " " << J << nl
209 << "eigenValues (principal moments) " << eVal << endl;
213 tensor eVec = eigenVectors(J);
215 Info<< "eigenVectors (principal axes, from normalised inertia) " << eVec
218 OFstream str("cell_" + name(cellI) + "_inertia.obj");
220 Info<< nl << "Writing scaled principal axes of cell " << cellI << " to "
221 << str.name() << endl;
223 const point& cC = mesh.cellCentres()[cellI];
227 (cC - mesh.faceCentres()[mesh.cells()[cellI][0]])
228 /eVal.component(findMin(eVal))
231 meshTools::writeOBJ(str, cC);
232 meshTools::writeOBJ(str, cC + scale*eVal.x()*eVec.x());
233 meshTools::writeOBJ(str, cC + scale*eVal.y()*eVec.y());
234 meshTools::writeOBJ(str, cC + scale*eVal.z()*eVec.z());
236 for (label i = 1; i < 4; i++)
238 str << "l " << 1 << ' ' << i + 1 << endl;
242 Info<< nl << "End" << nl << endl;
248 // ************************************************************************* //