From c8f03c3bef220ea80809dac3b729ce6e7ec4ae3c Mon Sep 17 00:00:00 2001 From: henry Date: Wed, 2 Dec 2009 22:42:10 +0000 Subject: [PATCH 1/3] Corrected and simplified solver selection for implicit MULES. --- src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C b/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C index 89e0260456..f938c2cf9f 100644 --- a/src/finiteVolume/fvMatrices/solvers/MULES/MULESTemplates.C +++ b/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; From 6b7dd8160d10df27a5b471e697ca7b8fbb6dbc92 Mon Sep 17 00:00:00 2001 From: graham Date: Fri, 4 Dec 2009 16:33:46 +0000 Subject: [PATCH 2/3] Adding a function to calculate the moment of inertia of a triangle to triangle, from: http://en.wikipedia.org/wiki/Inertia_tensor_of_triangle Adding the ability to calculate the inertia tensor of a polygon face by summing the triangle inertias. Added a test application to draw a test face with its principal axes and moments of inertia. --- applications/test/momentOfInertia/Make/files | 3 + .../test/momentOfInertia/Make/options | 5 + .../momentOfInertia/momentOfInertiaTest.C | 112 ++++++++++++++++++ src/OpenFOAM/meshes/meshShapes/face/face.C | 39 +++++- src/OpenFOAM/meshes/meshShapes/face/face.H | 9 ++ .../primitiveShapes/triangle/triangle.H | 11 +- .../primitiveShapes/triangle/triangleI.H | 36 ++++++ 7 files changed, 213 insertions(+), 2 deletions(-) create mode 100644 applications/test/momentOfInertia/Make/files create mode 100644 applications/test/momentOfInertia/Make/options create mode 100644 applications/test/momentOfInertia/momentOfInertiaTest.C diff --git a/applications/test/momentOfInertia/Make/files b/applications/test/momentOfInertia/Make/files new file mode 100644 index 0000000000..89848e34fb --- /dev/null +++ b/applications/test/momentOfInertia/Make/files @@ -0,0 +1,3 @@ +momentOfInertiaTest.C + +EXE = $(FOAM_USER_APPBIN)/momentOfInertiaTest diff --git a/applications/test/momentOfInertia/Make/options b/applications/test/momentOfInertia/Make/options new file mode 100644 index 0000000000..b89c1bfe86 --- /dev/null +++ b/applications/test/momentOfInertia/Make/options @@ -0,0 +1,5 @@ +EXE_INC = \ + -I$(LIB_SRC)/meshTools/lnInclude + +EXE_LIBS = \ + -lfiniteVolume diff --git a/applications/test/momentOfInertia/momentOfInertiaTest.C b/applications/test/momentOfInertia/momentOfInertiaTest.C new file mode 100644 index 0000000000..b4ea06bba2 --- /dev/null +++ b/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; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C index 6e3d88979e..7935c4a9b6 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.C +++ b/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 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // ************************************************************************* // - diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index 25a4f8fbce..c317b3eef9 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/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 diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 3bc3ea38a8..ae24678c4a 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/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. diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H index 40fa909b16..be4c3d95d0 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H @@ -331,6 +331,42 @@ inline scalar triangle::sweptVol(const triangle& t) const } +template +inline tensor triangle::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 inline pointHit triangle::ray ( From f23e3eb3ea98902ed45b7e60af28e2022df8f3f2 Mon Sep 17 00:00:00 2001 From: graham Date: Fri, 4 Dec 2009 17:11:30 +0000 Subject: [PATCH 3/3] Changing linked library from finiteVolume to meshTools. --- applications/test/momentOfInertia/Make/options | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/test/momentOfInertia/Make/options b/applications/test/momentOfInertia/Make/options index b89c1bfe86..54c035b8f5 100644 --- a/applications/test/momentOfInertia/Make/options +++ b/applications/test/momentOfInertia/Make/options @@ -2,4 +2,4 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude EXE_LIBS = \ - -lfiniteVolume + -lmeshTools