surfaceInertia.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 surfaceInertia
29
30Group
31 grpSurfaceUtilities
32
33Description
34 Calculates the inertia tensor, principal axes and moments of a
35 command line specified triSurface.
36
37 Inertia can either be of the solid body or of a thin shell.
38
39\*---------------------------------------------------------------------------*/
40
41#include "argList.H"
42#include "ListOps.H"
43#include "triSurface.H"
44#include "OFstream.H"
45#include "meshTools.H"
46#include "Random.H"
47#include "transform.H"
48#include "IOmanip.H"
49#include "Pair.H"
50#include "momentOfInertia.H"
51
52// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53
54using namespace Foam;
55
56int main(int argc, char *argv[])
57{
58 argList::addNote
59 (
60 "Calculates the inertia tensor and principal axes and moments"
61 " of the specified surface.\n"
62 "Inertia can either be of the solid body or of a thin shell."
63 );
64
65 argList::noParallel();
66 argList::addArgument("input", "The input surface file");
67
68 argList::addBoolOption
69 (
70 "shellProperties",
71 "Inertia of a thin shell"
72 );
73
74 argList::addOption
75 (
76 "density",
77 "scalar",
78 "Specify density, "
79 "kg/m3 for solid properties, kg/m2 for shell properties"
80 );
81
82 argList::addOption
83 (
84 "referencePoint",
85 "vector",
86 "Inertia relative to this point, not the centre of mass"
87 );
88
89 argList args(argc, argv);
90
91 const auto surfFileName = args.get<fileName>(1);
92 const scalar density = args.getOrDefault<scalar>("density", 1);
93
94 vector refPt = Zero;
95 bool calcAroundRefPt = args.readIfPresent("referencePoint", refPt);
96
97 const triSurface surf(surfFileName);
98
99 scalar m = 0.0;
100 vector cM = Zero;
101 tensor J = Zero;
102
103 if (args.found("shellProperties"))
104 {
105 momentOfInertia::massPropertiesShell(surf, density, m, cM, J);
106 }
107 else
108 {
109 momentOfInertia::massPropertiesSolid(surf, density, m, cM, J);
110 }
111
112 if (m < 0)
113 {
115 << "Negative mass detected, the surface may be inside-out." << endl;
116 }
117
118 vector eVal = eigenValues(symm(J));
119
120 tensor eVec = eigenVectors(symm(J));
121
122 label pertI = 0;
123
124 Random rand(57373);
125
126 while ((magSqr(eVal) < VSMALL) && pertI < 10)
127 {
129 << "No eigenValues found, shape may have symmetry, "
130 << "perturbing inertia tensor diagonal" << endl;
131
132 J.xx() *= 1.0 + SMALL*rand.sample01<scalar>();
133 J.yy() *= 1.0 + SMALL*rand.sample01<scalar>();
134 J.zz() *= 1.0 + SMALL*rand.sample01<scalar>();
135
136 eVal = eigenValues(symm(J));
137
138 eVec = eigenVectors(symm(J));
139
140 pertI++;
141 }
142
143 bool showTransform = true;
144
145 if
146 (
147 (mag(eVec.x() ^ eVec.y()) > (1.0 - SMALL))
148 && (mag(eVec.y() ^ eVec.z()) > (1.0 - SMALL))
149 && (mag(eVec.z() ^ eVec.x()) > (1.0 - SMALL))
150 )
151 {
152 // Make the eigenvectors a right handed orthogonal triplet
153 eVec = tensor
154 (
155 eVec.x(),
156 eVec.y(),
157 eVec.z() * sign((eVec.x() ^ eVec.y()) & eVec.z())
158 );
159
160 // Finding the most natural transformation. Using Lists
161 // rather than tensors to allow indexed permutation.
162
163 // Cartesian basis vectors - right handed orthogonal triplet
164 List<vector> cartesian(3);
165
166 cartesian[0] = vector(1, 0, 0);
167 cartesian[1] = vector(0, 1, 0);
168 cartesian[2] = vector(0, 0, 1);
169
170 // Principal axis basis vectors - right handed orthogonal
171 // triplet
172 List<vector> principal(3);
173
174 principal[0] = eVec.x();
175 principal[1] = eVec.y();
176 principal[2] = eVec.z();
177
178 scalar maxMagDotProduct = -GREAT;
179
180 // Matching axis indices, first: cartesian, second:principal
181
182 Pair<label> match(-1, -1);
183
184 forAll(cartesian, cI)
185 {
186 forAll(principal, pI)
187 {
188 scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
189
190 if (magDotProduct > maxMagDotProduct)
191 {
192 maxMagDotProduct = magDotProduct;
193
194 match.first() = cI;
195
196 match.second() = pI;
197 }
198 }
199 }
200
201 scalar sense = sign
202 (
203 cartesian[match.first()] & principal[match.second()]
204 );
205
206 if (sense < 0)
207 {
208 // Invert the best match direction and swap the order of
209 // the other two vectors
210
211 List<vector> tPrincipal = principal;
212
213 tPrincipal[match.second()] *= -1;
214
215 tPrincipal[(match.second() + 1) % 3] =
216 principal[(match.second() + 2) % 3];
217
218 tPrincipal[(match.second() + 2) % 3] =
219 principal[(match.second() + 1) % 3];
220
221 principal = tPrincipal;
222
223 vector tEVal = eVal;
224
225 tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
226
227 tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
228
229 eVal = tEVal;
230 }
231
232 label permutationDelta = match.second() - match.first();
233
234 if (permutationDelta != 0)
235 {
236 // Add 3 to the permutationDelta to avoid negative indices
237
238 permutationDelta += 3;
239
240 List<vector> tPrincipal = principal;
241
242 vector tEVal = eVal;
243
244 for (label i = 0; i < 3; i++)
245 {
246 tPrincipal[i] = principal[(i + permutationDelta) % 3];
247
248 tEVal[i] = eVal[(i + permutationDelta) % 3];
249 }
250
251 principal = tPrincipal;
252
253 eVal = tEVal;
254 }
255
256 label matchedAlready = match.first();
257
258 match =Pair<label>(-1, -1);
259
260 maxMagDotProduct = -GREAT;
261
262 forAll(cartesian, cI)
263 {
264 if (cI == matchedAlready)
265 {
266 continue;
267 }
268
269 forAll(principal, pI)
270 {
271 if (pI == matchedAlready)
272 {
273 continue;
274 }
275
276 scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
277
278 if (magDotProduct > maxMagDotProduct)
279 {
280 maxMagDotProduct = magDotProduct;
281
282 match.first() = cI;
283
284 match.second() = pI;
285 }
286 }
287 }
288
289 sense = sign
290 (
291 cartesian[match.first()] & principal[match.second()]
292 );
293
294 if (sense < 0 || (match.second() - match.first()) != 0)
295 {
296 principal[match.second()] *= -1;
297
298 List<vector> tPrincipal = principal;
299
300 tPrincipal[(matchedAlready + 1) % 3] =
301 principal[(matchedAlready + 2) % 3]*-sense;
302
303 tPrincipal[(matchedAlready + 2) % 3] =
304 principal[(matchedAlready + 1) % 3]*-sense;
305
306 principal = tPrincipal;
307
308 vector tEVal = eVal;
309
310 tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
311
312 tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
313
314 eVal = tEVal;
315 }
316
317 eVec = tensor(principal[0], principal[1], principal[2]);
318
319 // {
320 // tensor R = rotationTensor(vector(1, 0, 0), eVec.x());
321
322 // R = rotationTensor(R & vector(0, 1, 0), eVec.y()) & R;
323
324 // Info<< "R = " << nl << R << endl;
325
326 // Info<< "R - eVec.T() " << R - eVec.T() << endl;
327 // }
328 }
329 else
330 {
332 << "Non-unique eigenvectors, cannot compute transformation "
333 << "from Cartesian axes" << endl;
334
335 showTransform = false;
336 }
337
338 // calculate the total surface area
339
340 scalar surfaceArea = 0;
341
342 forAll(surf, facei)
343 {
344 const labelledTri& f = surf[facei];
345
346 if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
347 {
349 << "Illegal triangle " << facei << " vertices " << f
350 << " coords " << f.points(surf.points()) << endl;
351 }
352 else
353 {
354 surfaceArea += triPointRef
355 (
356 surf.points()[f[0]],
357 surf.points()[f[1]],
358 surf.points()[f[2]]
359 ).mag();
360 }
361 }
362
363 Info<< nl << setprecision(12)
364 << "Density: " << density << nl
365 << "Mass: " << m << nl
366 << "Centre of mass: " << cM << nl
367 << "Surface area: " << surfaceArea << nl
368 << "Inertia tensor around centre of mass: " << nl << J << nl
369 << "eigenValues (principal moments): " << eVal << nl
370 << "eigenVectors (principal axes): " << nl
371 << eVec.x() << nl << eVec.y() << nl << eVec.z() << endl;
372
373 if (showTransform)
374 {
375 Info<< "Transform tensor from reference state (orientation):" << nl
376 << eVec.T() << nl
377 << "Rotation tensor required to transform "
378 "from the body reference frame to the global "
379 "reference frame, i.e.:" << nl
380 << "globalVector = orientation & bodyLocalVector"
381 << nl;
382
383 Info<< nl
384 << "Entries for sixDoFRigidBodyDisplacement boundary condition:"
385 << nl;
386
387 Info()
388 .writeEntry("mass", m)
389 .writeEntry("centreOfMass", cM)
390 .writeEntry("momentOfInertia", eVal)
391 .writeEntry("orientation", eVec.T());
392 }
393
394 if (calcAroundRefPt)
395 {
396 Info<< nl << "Inertia tensor relative to " << refPt << ": " << nl
397 << momentOfInertia::applyParallelAxisTheorem(m, cM, J, refPt)
398 << endl;
399 }
400
401 OFstream str("axes.obj");
402
403 Info<< nl << "Writing scaled principal axes at centre of mass of "
404 << surfFileName << " to " << str.name() << endl;
405
406 scalar scale = mag(cM - surf.points()[0])/eVal.component(findMin(eVal));
407
408 meshTools::writeOBJ(str, cM);
409 meshTools::writeOBJ(str, cM + scale*eVal.x()*eVec.x());
410 meshTools::writeOBJ(str, cM + scale*eVal.y()*eVec.y());
411 meshTools::writeOBJ(str, cM + scale*eVal.z()*eVec.z());
412
413 for (label i = 1; i < 4; i++)
414 {
415 str << "l " << 1 << ' ' << i + 1 << endl;
416 }
417
418 Info<< "\nEnd\n" << endl;
419
420 return 0;
421}
422
423
424// ************************************************************************* //
Istream and Ostream manipulators taking arguments.
Various functions to operate on Lists.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
Output to file stream, using an OSstream.
Definition: OFstream.H:57
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: Pair.H:69
Random number generator.
Definition: Random.H:60
Extract command arguments and options from the supplied argc and argv parameters.
Definition: argList.H:124
T get(const label index) const
Get a value from the argument at index.
Definition: argListI.H:278
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
bool readIfPresent(const word &optName, T &val) const
Read a value from the named option if present.
Definition: argListI.H:323
T getOrDefault(const word &optName, const T &deflt) const
Get a value from the named option if present, or return default.
Definition: argListI.H:307
A class for handling file names.
Definition: fileName.H:76
A triFace with additional (region) index.
Definition: labelledTri.H:60
Tensor of scalars, i.e. Tensor<scalar>.
Triangulated surface description with patch information.
Definition: triSurface.H:79
A triangle primitive used to calculate face normals and swept volumes.
Definition: triangle.H:80
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:128
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define WarningInFunction
Report a warning using Foam::Warning.
bool match(const UList< wordRe > &patterns, const std::string &text)
Return true if text matches one of the regular expressions.
Definition: stringOps.H:76
Namespace for OpenFOAM.
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
label findMin(const ListType &input, label start=0)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setprecision(const int i)
Definition: IOmanip.H:205
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
Foam::argList args(argc, argv)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
3D tensor transformation operations.