momentOfInertia.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 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
27\*---------------------------------------------------------------------------*/
28
29#include "momentOfInertia.H"
31
32// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
33
35(
36 const pointField& pts,
37 const triFaceList& triFaces,
38 scalar density,
39 scalar& mass,
40 vector& cM,
41 tensor& J
42)
43{
44 // Reimplemented from: Wm4PolyhedralMassProperties.cpp
45 // File Version: 4.10.0 (2009/11/18)
46
47 // Geometric Tools, LC
48 // Copyright (c) 1998-2010
49 // Distributed under the Boost Software License, Version 1.0.
50 // http://www.boost.org/LICENSE_1_0.txt
51 // http://www.geometrictools.com/License/Boost/LICENSE_1_0.txt
52
53 // Boost Software License - Version 1.0 - August 17th, 2003
54
55 // Permission is hereby granted, free of charge, to any person or
56 // organization obtaining a copy of the software and accompanying
57 // documentation covered by this license (the "Software") to use,
58 // reproduce, display, distribute, execute, and transmit the
59 // Software, and to prepare derivative works of the Software, and
60 // to permit third-parties to whom the Software is furnished to do
61 // so, all subject to the following:
62
63 // The copyright notices in the Software and this entire
64 // statement, including the above license grant, this restriction
65 // and the following disclaimer, must be included in all copies of
66 // the Software, in whole or in part, and all derivative works of
67 // the Software, unless such copies or derivative works are solely
68 // in the form of machine-executable object code generated by a
69 // source language processor.
70
71 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
72 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
73 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
74 // NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR
75 // ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR
76 // OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
77 // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
78 // USE OR OTHER DEALINGS IN THE SOFTWARE.
79
80 const scalar r6 = 1.0/6.0;
81 const scalar r24 = 1.0/24.0;
82 const scalar r60 = 1.0/60.0;
83 const scalar r120 = 1.0/120.0;
84
85 // order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
86 scalarField integrals(10, Zero);
87
88 forAll(triFaces, i)
89 {
90 const triFace& tri(triFaces[i]);
91
92 // vertices of triangle i
93 vector v0 = pts[tri[0]];
94 vector v1 = pts[tri[1]];
95 vector v2 = pts[tri[2]];
96
97 // cross product of edges
98 vector eA = v1 - v0;
99 vector eB = v2 - v0;
100 vector n = eA ^ eB;
101
102 // compute integral terms
103 scalar tmp0, tmp1, tmp2;
104
105 scalar f1x, f2x, f3x, g0x, g1x, g2x;
106
107 tmp0 = v0.x() + v1.x();
108 f1x = tmp0 + v2.x();
109 tmp1 = v0.x()*v0.x();
110 tmp2 = tmp1 + v1.x()*tmp0;
111 f2x = tmp2 + v2.x()*f1x;
112 f3x = v0.x()*tmp1 + v1.x()*tmp2 + v2.x()*f2x;
113 g0x = f2x + v0.x()*(f1x + v0.x());
114 g1x = f2x + v1.x()*(f1x + v1.x());
115 g2x = f2x + v2.x()*(f1x + v2.x());
116
117 scalar f1y, f2y, f3y, g0y, g1y, g2y;
118
119 tmp0 = v0.y() + v1.y();
120 f1y = tmp0 + v2.y();
121 tmp1 = v0.y()*v0.y();
122 tmp2 = tmp1 + v1.y()*tmp0;
123 f2y = tmp2 + v2.y()*f1y;
124 f3y = v0.y()*tmp1 + v1.y()*tmp2 + v2.y()*f2y;
125 g0y = f2y + v0.y()*(f1y + v0.y());
126 g1y = f2y + v1.y()*(f1y + v1.y());
127 g2y = f2y + v2.y()*(f1y + v2.y());
128
129 scalar f1z, f2z, f3z, g0z, g1z, g2z;
130
131 tmp0 = v0.z() + v1.z();
132 f1z = tmp0 + v2.z();
133 tmp1 = v0.z()*v0.z();
134 tmp2 = tmp1 + v1.z()*tmp0;
135 f2z = tmp2 + v2.z()*f1z;
136 f3z = v0.z()*tmp1 + v1.z()*tmp2 + v2.z()*f2z;
137 g0z = f2z + v0.z()*(f1z + v0.z());
138 g1z = f2z + v1.z()*(f1z + v1.z());
139 g2z = f2z + v2.z()*(f1z + v2.z());
140
141 // update integrals
142 integrals[0] += n.x()*f1x;
143 integrals[1] += n.x()*f2x;
144 integrals[2] += n.y()*f2y;
145 integrals[3] += n.z()*f2z;
146 integrals[4] += n.x()*f3x;
147 integrals[5] += n.y()*f3y;
148 integrals[6] += n.z()*f3z;
149 integrals[7] += n.x()*(v0.y()*g0x + v1.y()*g1x + v2.y()*g2x);
150 integrals[8] += n.y()*(v0.z()*g0y + v1.z()*g1y + v2.z()*g2y);
151 integrals[9] += n.z()*(v0.x()*g0z + v1.x()*g1z + v2.x()*g2z);
152 }
153
154 integrals[0] *= r6;
155 integrals[1] *= r24;
156 integrals[2] *= r24;
157 integrals[3] *= r24;
158 integrals[4] *= r60;
159 integrals[5] *= r60;
160 integrals[6] *= r60;
161 integrals[7] *= r120;
162 integrals[8] *= r120;
163 integrals[9] *= r120;
164
165 // mass
166 mass = integrals[0];
167
168 // center of mass
169 cM = vector(integrals[1], integrals[2], integrals[3])/mass;
170
171 // inertia relative to origin
172 J.xx() = integrals[5] + integrals[6];
173 J.xy() = -integrals[7];
174 J.xz() = -integrals[9];
175 J.yx() = J.xy();
176 J.yy() = integrals[4] + integrals[6];
177 J.yz() = -integrals[8];
178 J.zx() = J.xz();
179 J.zy() = J.yz();
180 J.zz() = integrals[4] + integrals[5];
181
182 // inertia relative to center of mass
183 J -= mass*((cM & cM)*I - cM*cM);
184
185 // Apply density
186 mass *= density;
187 J *= density;
188}
189
190
192(
193 const pointField& pts,
194 const triFaceList& triFaces,
195 scalar density,
196 scalar& mass,
197 vector& cM,
198 tensor& J,
199 bool doReduce
200)
201{
202 // Reset properties for accumulation
203
204 mass = 0.0;
205 cM = Zero;
206 J = Zero;
207
208 // Find centre of mass
209
210 forAll(triFaces, i)
211 {
212 const triFace& tri(triFaces[i]);
213
215 (
216 pts[tri[0]],
217 pts[tri[1]],
218 pts[tri[2]]
219 );
220
221 scalar triMag = t.mag();
222
223 cM += triMag*t.centre();
224
225 mass += triMag;
226 }
227
228 if (doReduce)
229 {
230 reduce(cM, sumOp<vector>());
231 reduce(mass, sumOp<scalar>());
232 }
233
234 cM /= mass;
235
236 mass *= density;
237
238 // Find inertia around centre of mass
239
240 forAll(triFaces, i)
241 {
242 const triFace& tri(triFaces[i]);
243
244 J += triPointRef
245 (
246 pts[tri[0]],
247 pts[tri[1]],
248 pts[tri[2]]
249 ).inertia(cM, density);
250 }
251
252 if (doReduce)
253 {
254 reduce(J, sumOp<tensor>());
255 }
256}
257
258
260(
261 const triSurface& surf,
262 scalar density,
263 scalar& mass,
264 vector& cM,
265 tensor& J
266)
267{
268 triFaceList faces(surf.size());
269
270 forAll(surf, i)
271 {
272 faces[i] = triFace(surf[i]);
273 }
274
275 massPropertiesSolid(surf.points(), faces, density, mass, cM, J);
276}
277
278
280(
281 const triSurface& surf,
282 scalar density,
283 scalar& mass,
284 vector& cM,
285 tensor& J,
286 bool doReduce
287)
288{
289 triFaceList faces(surf.size());
290
291 forAll(surf, i)
292 {
293 faces[i] = triFace(surf[i]);
294 }
295
296 massPropertiesShell(surf.points(), faces, density, mass, cM, J, doReduce);
297}
298
299
301(
302 const polyPatch& pp,
303 scalar density,
304 scalar& mass,
305 vector& cM,
306 tensor& J,
307 bool doReduce
308)
309{
310 DynamicList<triFace> faces(3*pp.size());
311
312 // decompose patch faces using triangle fan
313 forAll(pp, faceI)
314 {
315 const face& f = pp[faceI];
316
317 if (f.size() > 2)
318 {
319 const label v0 = 0;
320
321 for (label i = 1; i < f.size() - 1; i++)
322 {
323 faces.append(triFace(f[v0], f[i],f[i + 1]));
324 }
325 }
326 }
327
328 triFaceList triFaces;
329 triFaces.transfer(faces);
330 massPropertiesShell(pp.points(), triFaces, density, mass, cM, J, doReduce);
331}
332
333
335(
336 scalar mass,
337 const vector& cM,
338 const tensor& J,
339 const vector& refPt
340)
341{
342 // The displacement vector (refPt = cM) is the displacement of the
343 // new reference point from the centre of mass of the body that
344 // the inertia tensor applies to.
345
346 vector d = (refPt - cM);
347
348 return J + mass*((d & d)*I - d*d);
349}
350
351
353(
354 const polyMesh& mesh
355)
356{
358
359 tensorField& tf = tTf.ref();
360
361 forAll(tf, cI)
362 {
363 tf[cI] = meshInertia(mesh, cI);
364 }
365
366 return tTf;
367}
368
369
371(
372 const polyMesh& mesh,
373 label celli
374)
375{
377 (
378 mesh,
379 celli
380 );
381
382 triFaceList faces(cellTets.size());
383
384 forAll(cellTets, cTI)
385 {
386 faces[cTI] = cellTets[cTI].faceTriIs(mesh);
387 }
388
389 scalar m = 0.0;
390 vector cM = Zero;
391 tensor J = Zero;
392
393 massPropertiesSolid(mesh.points(), faces, 1.0, m, cM, J);
394
395 return J;
396}
397
398
399// ************************************************************************* //
label n
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
void append(const T &val)
Copy append an element to the end of this list.
Definition: DynamicListI.H:503
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
void transfer(List< T > &list)
Definition: List.C:447
const Field< point_type > & points() const noexcept
Return reference to global points.
const Cmpt & xx() const
Definition: TensorI.H:153
const Cmpt & yx() const
Definition: TensorI.H:174
const Cmpt & yz() const
Definition: TensorI.H:188
const Cmpt & xz() const
Definition: TensorI.H:167
const Cmpt & zz() const
Definition: TensorI.H:209
const Cmpt & xy() const
Definition: TensorI.H:160
const Cmpt & zx() const
Definition: TensorI.H:195
const Cmpt & zy() const
Definition: TensorI.H:202
const Cmpt & yy() const
Definition: TensorI.H:181
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
const Cmpt & z() const
Access to the vector z component.
Definition: VectorI.H:85
const Cmpt & y() const
Access to the vector y component.
Definition: VectorI.H:79
const Cmpt & x() const
Access to the vector x component.
Definition: VectorI.H:73
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
static void massPropertiesShell(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J, bool doReduce=false)
static tensor applyParallelAxisTheorem(scalar mass, const vector &cM, const tensor &J, const vector &refPt)
static tmp< tensorField > meshInertia(const polyMesh &mesh)
static void massPropertiesPatch(const polyPatch &pp, scalar density, scalar &mass, vector &cM, tensor &J, bool doReduce=false)
static void massPropertiesSolid(const pointField &pts, const triFaceList &triFaces, scalar density, scalar &mass, vector &cM, tensor &J)
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
label nCells() const noexcept
Number of mesh cells.
A class for managing temporary objects.
Definition: tmp.H:65
T & ref() const
Definition: tmpI.H:227
A triangular face using a FixedList of labels corresponding to mesh vertices.
Definition: triFace.H:72
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
Point centre() const
Return centre (centroid)
Definition: triangleI.H:105
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:128
dynamicFvMesh & mesh
static const Identity< scalar > I
Definition: Identity.H:94
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
triangle< point, const point & > triPointRef
A triangle using referred points.
Definition: triangle.H:71
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Field< tensor > tensorField
Specialisation of Field<T> for tensor.
Vector< scalar > vector
Definition: vector.H:61
labelList f(nPoints)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333