exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Implementation of constraints, part 2

Browse source
This commit is contained in:
Vuko Vukcevic 2017-03-07 10:28:43 +01:00 committed by Hrvoje Jasak
parent ef9664c191
commit 488ccffca9
11 changed files with 588 additions and 37 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
src/ODE/Make/files
View file

@ -16,6 +16,11 @@ sixDOF = sixDOF
$(sixDOF)/finiteRotation/finiteRotation.C
$(sixDOF)/sixDOFqODE/sixDOFqODE.C
$(sixDOF)/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C
$(sixDOF)/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C
$(sixDOF)/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C
$(sixDOF)/sixDOFODE/sixDOFODE.C
$(sixDOF)/sixDOFODE/newSixDOFODE.C
$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C

29
src/ODE/sixDOF/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C
View file

@ -63,15 +63,36 @@ Foam::constantAngularAcceleration::~constantAngularAcceleration()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::vector
Foam::constantAngularAcceleration::matrixContribution() const
Foam::vector Foam::constantAngularAcceleration::matrixContribution
(
const scalar,
const tensor& toRelative
) const
{
vector mc;
if (inGlobal_)
{
// Constraint is given in global (inertial) coordinate system, transform
// it to local
mc = toRelative & dir_;
}
else
{
// Constraint already in local (body) coordinate system
mc = dir_;
}
return mc;
}
const Foam::vector
Foam::constantAngularAcceleration::sourceContribution() const
Foam::scalar Foam::constantAngularAcceleration::sourceContribution
(
const scalar
) const
{
return alpha_;
}

8
src/ODE/sixDOF/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.H
View file

@ -105,10 +105,14 @@ public:
// Constraint specific functions
//- Return matrix contribution defined by constraint, f(t)
const vector matrixContribution() const = 0;
virtual vector matrixContribution
(
const scalar,
const tensor& toRelative
) const;
//- Return source contribution defined by constraint, a(t)
const vector sourceContribution() const = 0;
virtual scalar sourceContribution(const scalar) const;
};

2
src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationConstraint.C → src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C
View file

@ -65,7 +65,7 @@ Foam::autoPtr<Foam::rotationalConstraint> Foam::rotationalConstraint::New
const word constraintType(dict.lookup("type"));
wordConstructorTable::iterator cstrIter =
wordConstructorTablePtr_->find(constraintTypeType);
wordConstructorTablePtr_->find(constraintType);
if (cstrIter == wordConstructorTablePtr_->end())
{

8
src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationConstraint.H → src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.H
View file

@ -178,10 +178,14 @@ public:
// Constraint specific functions
//- Return matrix contribution defined by constraint, f(t)
const vector matrixContribution() const = 0;
virtual vector matrixContribution
(
const scalar t,
const tensor& toRelative
) const = 0;
//- Return source contribution defined by constraint, a(t)
const vector sourceContribution() const = 0;
virtual scalar sourceContribution(const scalar t) const = 0;
};

84
src/ODE/sixDOF/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C Normal file
View file

@ -0,0 +1,84 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "constantTranslationalAcceleration.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(constantTranslationalAcceleration, 0);
addToRunTimeSelectionTable
(
translationalConstraint,
constantTranslationalAcceleration,
word
);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::constantTranslationalAcceleration::constantTranslationalAcceleration
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
)
:
translationalConstraint(name, sixDOFODE, dict),
dir_(dict.lookup("constraintDirection")),
a_(readScalar(dict.lookup("translationalAcceleration")))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::constantTranslationalAcceleration::~constantTranslationalAcceleration()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::vector Foam::constantTranslationalAcceleration::matrixContribution
(
const scalar,
const tensor&
) const
{
return dir_;
}
Foam::scalar Foam::constantTranslationalAcceleration::sourceContribution
(
const scalar
) const
{
return a_;
}
// ************************************************************************* //

126
src/ODE/sixDOF/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.H Normal file
View file

@ -0,0 +1,126 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::constantTranslationalAcceleration
Description
Translational constraint defined by constant translational acceleration
g(vDot) = vDot - a = 0.
Author
Viktor Pandza, FSB Zagreb. All rights reserved.
Vuko Vukcevic, FSB Zagreb. All rights reserved.
SourceFiles
constantTranslationalAcceleration.C
\*---------------------------------------------------------------------------*/
#ifndef constantTranslationalAcceleration_H
#define constantTranslationalAcceleration_H
#include "translationalConstraint.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class constantTranslationalAcceleration Declaration
\*---------------------------------------------------------------------------*/
class constantTranslationalAcceleration
:
public translationalConstraint
{
// Private data
//- Direction of the constraint (unit vector)
const vector dir_;
//- Constant value of the translational acceleration
const scalar a_;
// Private Member Functions
//- Disallow default bitwise copy construct
constantTranslationalAcceleration
(
const constantTranslationalAcceleration&
);
//- Disallow default bitwise assignment
void operator=(const constantTranslationalAcceleration&);
public:
//- Runtime type information
TypeName("constantTranslationalAcceleration");
// Constructors
//- Construct from dictionary
constantTranslationalAcceleration
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
);
// Destructor
virtual ~constantTranslationalAcceleration();
// Member Functions
// Constraint specific functions
//- Return matrix contribution defined by constraint, f(t)
virtual vector matrixContribution
(
const scalar,
const tensor&
) const;
//- Return source contribution defined by constraint, a(t)
virtual scalar sourceContribution(const scalar) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

99
src/ODE/sixDOF/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C Normal file
View file

@ -0,0 +1,99 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "translationalConstraint.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(translationalConstraint, 0);
defineRunTimeSelectionTable(translationalConstraint, word);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::translationalConstraint::translationalConstraint
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
)
:
name_(name),
sixDOFODE_(sixDOFODE)
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::translationalConstraint::~translationalConstraint()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::autoPtr<Foam::translationalConstraint> Foam::translationalConstraint::New
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
)
{
const word constraintType(dict.lookup("type"));
wordConstructorTable::iterator cstrIter =
wordConstructorTablePtr_->find(constraintType);
if (cstrIter == wordConstructorTablePtr_->end())
{
FatalErrorIn
(
"translationalConstraint::New"
"\n("
"\n const word& name,"
"\n const sixDOFODE& sixDOFODE,"
"\n const dictionary& dict,"
"\n)"
) << "Unknown translation constraint type: " << constraintType
<< endl << endl
<< "Valid translation constraint types are: " << endl
<< wordConstructorTablePtr_->sortedToc()
<< exit(FatalError);
}
return autoPtr<translationalConstraint>
(
cstrIter()
(
name,
sixDOFODE,
dict
)
);
}
// ************************************************************************* //

200
src/ODE/sixDOF/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H Normal file
View file

@ -0,0 +1,200 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration | Version: 4.0
\\ / A nd | Web: http://www.foam-extend.org
\\/ M anipulation | For copyright notice see file Copyright
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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 3 of the License, or (at your
option) any later version.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::translationalConstraint
Description
Abstract base class containing interface for translational constraints used
within sixDOFODE classes.
The constraint is implicitly defined as:
g(vDot, t) = f(t)*vDot + a(t) = 0,
where vDot is translational acceleration.
Interface provides all the necessary data for inserting the constraint into
the resulting linear system via Lagrangian multipliers:
1. matrixContribution() - corresponding to f(t) (prefactor multiplying
vDot), to be inserted into the matrix.
2. sourceContribution() - corresponding to a(t), to be inserted into right
hand side vector.
Author
Viktor Pandza, FSB Zagreb. All rights reserved.
Vuko Vukcevic, FSB Zagreb. All rights reserved.
SourceFiles
translationalConstraint.C
\*---------------------------------------------------------------------------*/
#ifndef translationalConstraint_H
#define translationalConstraint_H
#include "sixDOFODE.H"
#include "typeInfo.H"
#include "runTimeSelectionTables.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class translationalConstraint Declaration
\*---------------------------------------------------------------------------*/
class translationalConstraint
{
// Private data
//- Name of the constraint
const word name_;
//- Reference to underlying sixDOFODE
const sixDOFODE& sixDOFODE_;
// Private Member Functions
//- Disallow default bitwise copy construct
translationalConstraint(const translationalConstraint&);
//- Disallow default bitwise assignment
void operator=(const translationalConstraint&);
public:
//- Runtime type information
TypeName("translationalConstraint");
// Declare run-time constructor selection table
declareRunTimeSelectionTable
(
autoPtr,
translationalConstraint,
word,
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
),
(name, sixDOFODE, dict)
);
//- Class used for the read-construction of
// PtrLists of translationalConstraint
class iNew
{
const sixDOFODE& sixDOFODE_;
public:
iNew(const sixDOFODE& sixDOFODE)
:
sixDOFODE_(sixDOFODE)
{}
autoPtr<translationalConstraint> operator()(Istream& is) const
{
word name(is);
dictionary dict(is);
return translationalConstraint::New(name, sixDOFODE_, dict);
}
};
// Constructors
//- Construct from dictionary
translationalConstraint
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
);
// Selectors
//- Return a reference to the selected translationalConstraint
static autoPtr<translationalConstraint> New
(
const word& name,
const sixDOFODE& sixDOFODE,
const dictionary& dict
);
// Destructor
virtual ~translationalConstraint();
// Member Functions
// Access functions
//- Return const reference to name of the constraint
const word& name() const
{
return name_;
}
//- Return const reference to underlying sixDOFODE object
const sixDOFODE& sixDOF() const
{
return sixDOFODE_;
}
// Constraint specific functions
//- Return matrix contribution defined by constraint, f(t)
virtual vector matrixContribution
(
const scalar t,
const tensor& toRelative
) const = 0;
//- Return source contribution defined by constraint, a(t)
virtual scalar sourceContribution(const scalar t) const = 0;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

50
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C
View file

@ -72,7 +72,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
(
const dimensionedVector& xR,
const dimensionedVector& uR,
const tensor& R
const tensor& R,
const scalar t
) const
{
// Create a scalar square matrix representing Newton equations with
@ -87,7 +88,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
(
force() // External force
- (linSpringCoeffs() & xR) // Spring force
- (linDampingCoeffs() & uR); // Damping force
- (linDampingCoeffs() & uR) // Damping force
);
const vector& efVal = explicitForcing.value();
const scalar& m = mass().value();
@ -103,23 +104,23 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
forAll(translationalConstraints_, tcI)
{
// Get reference to current constraint
const translationalConstraint& curTc = translationalConstraint_[tcI];
const translationalConstraint& curTc = translationalConstraints_[tcI];
// Get matrix contribution from constraint
const vector mc = curTc.matrixContribution();
const vector mc = curTc.matrixContribution(t, R);
// Get matrix index
const label index = tcI + 3;
// Insert contributions into the matrix
forAll(lower, dirI)
forAll(mc, dirI)
{
M[dirI][index] = mc[dirI];
M[index][dirI] = mc[dirI];
}
// Insert source contribution (remainder of the constraint function)
rhs[index] = curTc.sourceContribution();
rhs[index] = curTc.sourceContribution(t);
}
// Solve the matrix using LU decomposition. Note: solution is in the rhs and
@ -132,7 +133,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
(
"A",
force().dimensions()/mass().dimensions(),
vector(rhs.x(), rhs.y(), rhs.z())
vector(rhs[0], rhs[1], rhs[2])
);
}
@ -140,7 +141,8 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
(
const tensor& R,
const dimensionedVector& omega
const dimensionedVector& omega,
const scalar t
) const
{
// Create a scalar square matrix representing Euler equations with
@ -149,11 +151,14 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
scalarField rhs(rotationalConstraints_.size() + 3, 0.0);
scalarSquareMatrix J(rhs.size(), 0.0);
// Get current inertial-to-local transformation
const dimensionedTensor RT("RT", dimless, R.T());
// Insert moment of inertia and explicit forcing into the system
const dimensionedVector explicitForcing
(
E(omega) // Euler part
+ (dimensionedTensor("R", dimless, R) & moment()) // External torque
+ (RT & moment()) // External torque
);
const vector& efVal = explicitForcing.value();
const diagTensor& I = momentOfInertia().value();
@ -172,20 +177,20 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
const rotationalConstraint& curRc = rotationalConstraints_[rcI];
// Get matrix contribution from the constraint
const vector mc = curRc.matrixContribution();
const vector mc = curRc.matrixContribution(t, RT.value());
// Get matrix index
const label index = rcI + 3;
// Insert contributions into the matrix
forAll(upper, dirI)
forAll(mc, dirI)
{
J[dirI][index] = mc[dirI];
J[index][dirI] = mc[dirI];
}
// Insert source contribution (remainder of the constraint function)
rhs[index] = curRc.sourceContribution();
rhs[index] = curRc.sourceContribution(t);
}
// Solve the matrix using LU decomposition. Note: solution is in the rhs and
@ -198,7 +203,7 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
(
"OmegaDot",
moment().dimensions()/momentOfInertia().dimensions(),
vector(rhs.x(), rhs.y(), rhs.z())
vector(rhs[0], rhs[1], rhs[2])
);
}
@ -307,8 +312,8 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
// HJ, 23/Mar/2015
coeffs_ = gsd.coeffs_;
translationalConstraints = gsd.translationalConstraints_;
rotationalConstraints = gsd.rotationalConstraints_;
translationalConstraints_ = gsd.translationalConstraints_;
rotationalConstraints_ = gsd.rotationalConstraints_;
}
@ -365,10 +370,10 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
// Read translation constraints if they are present
if (dict().found("translationalConstraints"))
{
PtrList<translationalConstraints> tcList
PtrList<translationalConstraint> tcList
(
dict().lookup("translationalConstraints"),
translationalConstraints::iNew(*this)
translationalConstraint::iNew(*this)
);
translationalConstraints_.transfer(tcList);
}
@ -376,10 +381,10 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
// Read rotation constraints if they are present
if (dict().found("rotationalConstraints"))
{
PtrList<rotationalConstraints> tcList
PtrList<rotationalConstraint> tcList
(
dict().lookup("rotationalConstraints"),
rotationalConstraints::iNew(*this)
rotationalConstraint::iNew(*this)
);
rotationalConstraints_.transfer(tcList);
}
@ -404,7 +409,7 @@ Foam::geometricSixDOF::geometricSixDOF
coeffs_(gsd.coeffs_),
translationalConstraints_(gsd.translationalConstraints_),
rotationalConstraints_(gsd.rotationalConstraints_),
rotationalConstraints_(gsd.rotationalConstraints_)
{}
@ -511,7 +516,7 @@ void Foam::geometricSixDOF::derivatives
const tensor curRot = (rotation_ & expMap(rotIncrementVector));
// Calculate translational acceleration using current rotation
const vector accel = A(curX, curU, curRot).value();
const vector accel = A(curX, curU, curRot, x).value();
// Set the derivatives for velocity
dydx[3] = accel.x();
@ -528,7 +533,7 @@ void Foam::geometricSixDOF::derivatives
);
// Calculate rotational acceleration using current rotation
const vector omegaDot = OmegaDot(curRot, curOmega).value();
const vector omegaDot = OmegaDot(curRot, curOmega, x).value();
dydx[6] = omegaDot.x();
dydx[7] = omegaDot.y();
@ -569,7 +574,6 @@ void Foam::geometricSixDOF::update(const scalar delta)
Uval.z() = coeffs_[5];
// Constrain velocity and re-set coefficients
constrainTranslation(Uval);
coeffs_[3] = Uval.x();
coeffs_[4] = Uval.y();
coeffs_[5] = Uval.z();

14
src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H
View file

@ -116,7 +116,7 @@ class geometricSixDOF
PtrList<translationalConstraint> translationalConstraints_;
//- List of rotational constraints
Ptrlist<rotationalConstraint> rotationalConstraints_;
PtrList<rotationalConstraint> rotationalConstraints_;
// Private Member Functions
@ -131,21 +131,25 @@ class geometricSixDOF
// Variables in relative coordinate system (solved for)
//- Return acceleration in relative coordinate system
// given current values of relative displacement and velocity
// given current values of relative displacement, velocity,
// orientation (rotation tensor) and time
dimensionedVector A
(
const dimensionedVector& xR,
const dimensionedVector& uR,
const tensor& R
const tensor& R,
const scalar t
) const;
//- Return rotational acceleration in relative coordinate system
// given current values for relative rotational velocity
// given current values for relative rotational velocity,
// orientation (rotation tensor) and time
dimensionedVector OmegaDot
(
const tensor& R,
const dimensionedVector& omega
const dimensionedVector& omega,
const scalar t
) const;
//- Return the Euler part of moment equation