Commit 4747a0c2 authored by mattijs's avatar mattijs
Browse files

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

parents d7d97f1f b98a01b2
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -41,6 +41,9 @@ Description
#include "OFstream.H"
#include "meshTools.H"
#include "Random.H"
#include "transform.H"
#include "IOmanip.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -355,6 +358,12 @@ int main(int argc, char *argv[])
);
}
if (m < 0)
{
WarningIn(args.executable() + "::main")
<< "Negative mass detected" << endl;
}
vector eVal = eigenValues(J);
tensor eVec = eigenVectors(J);
......@@ -380,19 +389,221 @@ int main(int argc, char *argv[])
pertI++;
}
Info<< nl
<< "Density = " << density << nl
<< "Mass = " << m << nl
<< "Centre of mass = " << cM << nl
<< "Inertia tensor around centre of mass = " << J << nl
<< "eigenValues (principal moments) = " << eVal << nl
<< "eigenVectors (principal axes) = "
<< eVec.x() << ' ' << eVec.y() << ' ' << eVec.z()
<< endl;
bool showTransform = true;
if
(
(mag(eVec.x() ^ eVec.y()) > (1.0 - SMALL))
&& (mag(eVec.y() ^ eVec.z()) > (1.0 - SMALL))
&& (mag(eVec.z() ^ eVec.x()) > (1.0 - SMALL))
)
{
// Make the eigenvectors a right handed orthogonal triplet
eVec.z() *= sign((eVec.x() ^ eVec.y()) & eVec.z());
// Finding the most natural transformation. Using Lists
// rather than tensors to allow indexed permutation.
// Cartesian basis vectors - right handed orthogonal triplet
List<vector> cartesian(3);
cartesian[0] = vector(1, 0, 0);
cartesian[1] = vector(0, 1, 0);
cartesian[2] = vector(0, 0, 1);
// Principal axis basis vectors - right handed orthogonal
// triplet
List<vector> principal(3);
principal[0] = eVec.x();
principal[1] = eVec.y();
principal[2] = eVec.z();
scalar maxMagDotProduct = -GREAT;
// Matching axis indices, first: cartesian, second:principal
Pair<label> match(-1, -1);
forAll(cartesian, cI)
{
forAll(principal, pI)
{
scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
if (magDotProduct > maxMagDotProduct)
{
maxMagDotProduct = magDotProduct;
match.first() = cI;
match.second() = pI;
}
}
}
scalar sense = sign
(
cartesian[match.first()] & principal[match.second()]
);
if (sense < 0)
{
// Invert the best match direction and swap the order of
// the other two vectors
List<vector> tPrincipal = principal;
tPrincipal[match.second()] *= -1;
tPrincipal[(match.second() + 1) % 3] =
principal[(match.second() + 2) % 3];
tPrincipal[(match.second() + 2) % 3] =
principal[(match.second() + 1) % 3];
principal = tPrincipal;
vector tEVal = eVal;
tEVal[(match.second() + 1) % 3] = eVal[(match.second() + 2) % 3];
tEVal[(match.second() + 2) % 3] = eVal[(match.second() + 1) % 3];
eVal = tEVal;
}
label permutationDelta = match.second() - match.first();
if (permutationDelta != 0)
{
// Add 3 to the permutationDelta to avoid negative indices
permutationDelta += 3;
List<vector> tPrincipal = principal;
vector tEVal = eVal;
for (label i = 0; i < 3; i++)
{
tPrincipal[i] = principal[(i + permutationDelta) % 3];
tEVal[i] = eVal[(i + permutationDelta) % 3];
}
principal = tPrincipal;
eVal = tEVal;
}
label matchedAlready = match.first();
match =Pair<label>(-1, -1);
maxMagDotProduct = -GREAT;
forAll(cartesian, cI)
{
if (cI == matchedAlready)
{
continue;
}
forAll(principal, pI)
{
if (pI == matchedAlready)
{
continue;
}
scalar magDotProduct = mag(cartesian[cI] & principal[pI]);
if (magDotProduct > maxMagDotProduct)
{
maxMagDotProduct = magDotProduct;
match.first() = cI;
match.second() = pI;
}
}
}
sense = sign
(
cartesian[match.first()] & principal[match.second()]
);
if (sense < 0 || (match.second() - match.first()) != 0)
{
principal[match.second()] *= -1;
List<vector> tPrincipal = principal;
tPrincipal[(matchedAlready + 1) % 3] =
principal[(matchedAlready + 2) % 3]*-sense;
tPrincipal[(matchedAlready + 2) % 3] =
principal[(matchedAlready + 1) % 3]*-sense;
principal = tPrincipal;
vector tEVal = eVal;
tEVal[(matchedAlready + 1) % 3] = eVal[(matchedAlready + 2) % 3];
tEVal[(matchedAlready + 2) % 3] = eVal[(matchedAlready + 1) % 3];
eVal = tEVal;
}
eVec.x() = principal[0];
eVec.y() = principal[1];
eVec.z() = principal[2];
// {
// tensor R = rotationTensor(vector(1, 0, 0), eVec.x());
// R = rotationTensor(R & vector(0, 1, 0), eVec.y()) & R;
// Info<< "R = " << nl << R << endl;
// Info<< "R - eVec.T() " << R - eVec.T() << endl;
// }
}
else
{
WarningIn(args.executable() + "::main")
<< "Non-unique eigenvectors, cannot compute transformation "
<< "from Cartesian axes" << endl;
showTransform = false;
}
Info<< nl << setprecision(10)
<< "Density: " << density << nl
<< "Mass: " << m << nl
<< "Centre of mass: " << cM << nl
<< "Inertia tensor around centre of mass: " << nl << J << nl
<< "eigenValues (principal moments): " << eVal << nl
<< "eigenVectors (principal axes): " << nl
<< eVec.x() << nl << eVec.y() << nl << eVec.z() << endl;
if (showTransform)
{
Info<< "Transform tensor from reference state (Q). " << nl
<< "Rotation tensor required to transform "
"from the body reference frame to the global "
"reference frame, i.e.:" << nl
<< "globalVector = Q & bodyLocalVector"
<< nl << eVec.T()
<< endl;
}
if (calcAroundRefPt)
{
Info << "Inertia tensor relative to " << refPt << " = "
Info << "Inertia tensor relative to " << refPt << ": "
<< applyParallelAxisTheorem(m, cM, J, refPt)
<< endl;
}
......
......@@ -30,6 +30,11 @@ License
void Foam::sixDoFRigidBodyMotion::applyRestraints()
{
if (restraints_.size() == 0)
{
return;
}
if (Pstream::master())
{
forAll(restraints_, rI)
......@@ -57,6 +62,11 @@ void Foam::sixDoFRigidBodyMotion::applyRestraints()
void Foam::sixDoFRigidBodyMotion::applyConstraints(scalar deltaT)
{
if (constraints_.size() == 0)
{
return;
}
if (Pstream::master())
{
label iter = 0;
......
......@@ -69,6 +69,9 @@ bool Foam::sixDoFRigidBodyMotionConstraint::read
const dictionary& sDoFRBMCDict
)
{
name_ =
fileName(sDoFRBMCDict.name().name()).components(token::COLON).last();
tolerance_ = (readScalar(sDoFRBMCDict.lookup("tolerance")));
relaxationFactor_ = sDoFRBMCDict.lookupOrDefault<scalar>
......
......@@ -64,6 +64,9 @@ bool Foam::sixDoFRigidBodyMotionRestraint::read
const dictionary& sDoFRBMRDict
)
{
name_ =
fileName(sDoFRBMRDict.name().name()).components(token::COLON).last();
sDoFRBMRCoeffs_ = sDoFRBMRDict.subDict(type() + "Coeffs");
return true;
......
......@@ -35,7 +35,7 @@ boundaryField
{
type sixDoFRigidBodyDisplacement;
centreOfMass (0.5 0.5 0.5);
momentOfInertia (0.08622222 0.8622222 0.144);
momentOfInertia (0.08622222 0.08622222 0.144);
mass 9.6;
rhoInf 1; // for forces calculation
// See sixDoFRigidBodyMotionState
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment