Merge branch 'master' of /home/noisy3/OpenFOAM/OpenFOAM-dev

applications/test/momentOfInertia/Make/files

@@ -0,0 +1,3 @@
+momentOfInertiaTest.C
+
+EXE = $(FOAM_USER_APPBIN)/momentOfInertiaTest

applications/test/momentOfInertia/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools

applications/test/momentOfInertia/momentOfInertiaTest.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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
+    test face.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ListOps.H"
+#include "face.H"
+#include "OFstream.H"
+#include "meshTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+using namespace Foam;
+
+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]);
+    }
+
+    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;
+    }
+
+    Info<< nl << "End" << nl << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -691,6 +691,44 @@ Foam::scalar Foam::face::sweptVol
 }
 
 
+Foam::tensor Foam::face::inertia
+(
+    const pointField& p,
+    const point& refPt,
+    scalar density
+) const
+{
+    // If the face is a triangle, do a direct calculation
+    if (size() == 3)
+    {
+        return triPointRef
+        (
+            p[operator[](0)],
+            p[operator[](1)],
+            p[operator[](2)]
+        ).inertia(refPt, density);
+    }
+
+    point c = centre(p);
+
+    tensor J = tensor::zero;
+
+    forAll(*this, i)
+    {
+        triPointRef t
+        (
+            p[operator[](i)],
+            p[operator[](fcIndex(i))],
+            c
+        );
+
+        J += t.inertia(refPt, density);
+    }
+
+    return J;
+}
+
+
 Foam::edgeList Foam::face::edges() const
 {
     const labelList& points = *this;
@@ -808,4 +846,3 @@ Foam::label Foam::face::trianglesQuads
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //
-

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -200,6 +200,15 @@ public:
             const pointField& newPoints
         ) const;
 
+        //- Return the inertia tensor, with optional reference
+        //  point and density specification
+        tensor inertia
+        (
+            const pointField&,
+            const point& refPt = vector::zero,
+            scalar density = 1.0
+        ) const;
+
         //- Return potential intersection with face with a ray starting
         //  at p, direction n (does not need to be normalized)
         //  Does face-center decomposition and returns triangle intersection

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -38,6 +38,7 @@ SourceFiles
 
 #include "intersection.H"
 #include "vector.H"
+#include "tensor.H"
 #include "pointHit.H"
 
 
@@ -100,7 +101,7 @@ public:
     //- Return types for classify
     enum proxType
     {
-        NONE, 
+        NONE,
         POINT,  // Close to point
         EDGE    // Close to edge
     };
@@ -152,6 +153,14 @@ public:
             //- Return swept-volume
             inline scalar sweptVol(const triangle& t) const;
 
+            //- Return the inertia tensor, with optional reference
+            //  point and density specification
+            inline tensor inertia
+            (
+                PointRef refPt = vector::zero,
+                scalar density = 1.0
+            ) const;
+
             //- Return point intersection with a ray.
             //  For a hit, the distance is signed. Positive number
             //  represents the point in front of triangle.

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -331,6 +331,42 @@ inline scalar triangle<Point, PointRef>::sweptVol(const triangle& t) const
 }
 
 
+template<class Point, class PointRef>
+inline tensor triangle<Point, PointRef>::inertia
+(
+    PointRef refPt,
+    scalar density
+) const
+{
+    Point aRel = a_ - refPt;
+    Point bRel = b_ - refPt;
+    Point cRel = c_ - refPt;
+
+    tensor V
+    (
+        aRel.x(), aRel.y(), aRel.z(),
+        bRel.x(), bRel.y(), bRel.z(),
+        cRel.x(), cRel.y(), cRel.z()
+    );
+
+    scalar a = Foam::mag((b_ - a_)^(c_ - a_));
+
+    tensor S = 1/24.0*(tensor::one + I);
+
+    return
+    (
+        a*I/24.0*
+        (
+            (aRel & aRel)
+          + (bRel & bRel)
+          + (cRel & cRel)
+          + ((aRel + bRel + cRel) & (aRel + bRel + cRel))
+        )
+      - a*(V.T() & S & V)
+    )
+   *density;
+}
+
 template<class Point, class PointRef>
 inline pointHit triangle<Point, PointRef>::ray
 (

src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C

@@ -142,8 +142,7 @@ void Foam::MULES::implicitSolve
 {
     const fvMesh& mesh = psi.mesh();
 
-    const dictionary& MULEScontrols = 
-       mesh.solverDict(psi.name()).subDict("MULESImplicit");
+    const dictionary& MULEScontrols = mesh.solverDict(psi.name());
 
     label maxIter
     (
@@ -255,7 +254,7 @@ void Foam::MULES::implicitSolve
         solve
         (
             psiConvectionDiffusion + fvc::div(lambda*phiCorr),
-            MULEScontrols.lookup("solver")
+            MULEScontrols
         );
 
         scalar maxPsiM1 = gMax(psi.internalField()) - 1.0;