Implementation of constraints, part 1

src/ODE/sixDOF/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "constantAngularAcceleration.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(constantAngularAcceleration, 0);
+    addToRunTimeSelectionTable
+    (
+        rotationalConstraint,
+        constantAngularAcceleration,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::constantAngularAcceleration::constantAngularAcceleration
+(
+    const word& name,
+    const sixDOFODE& sixDOFODE,
+    const dictionary& dict
+)
+:
+    rotationalConstraint(name, sixDOFODE, dict),
+    dir_(dict.lookup("constraintDirection")),
+    alpha_(readScalar(dict.lookup("angularAcceleration"))),
+    inGlobal_(dict.lookup("inGlobalCoordinateSystem"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantAngularAcceleration::~constantAngularAcceleration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::vector
+Foam::constantAngularAcceleration::matrixContribution() const
+{
+}
+
+
+const Foam::vector
+Foam::constantAngularAcceleration::sourceContribution() const
+{
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.H

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::constantAngularAcceleration
+
+Description
+    Rotational constraint defined by constant angular acceleration:
+
+    g(omegaDot) = omegaDot - a = 0.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    constantAngularAcceleration.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantAngularAcceleration_H
+#define constantAngularAcceleration_H
+
+#include "rotationalConstraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                 Class constantAngularAcceleration Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantAngularAcceleration
+:
+    public rotationalConstraint
+{
+    // Private data
+
+        //- Direction of the constraint (unit vector)
+        const vector dir_;
+
+        //- Constant value of angular acceleration
+        const scalar alpha_;
+
+        //- Switch whether the constraint should be applied in local or global
+        //  coordinate system
+        const Switch inGlobal_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        constantAngularAcceleration(const constantAngularAcceleration&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const constantAngularAcceleration&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantAngularAcceleration");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantAngularAcceleration
+        (
+            const word& name,
+            const sixDOFODE& sixDOFODE,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~constantAngularAcceleration();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            const vector matrixContribution() const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            const vector sourceContribution() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationConstraint.C

@@ -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 "rotationalConstraint.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rotationalConstraint, 0);
+    defineRunTimeSelectionTable(rotationalConstraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::rotationalConstraint
+(
+    const word& name,
+    const sixDOFODE& sixDOFODE,
+    const dictionary& dict
+)
+:
+    name_(name),
+    sixDOFODE_(sixDOFODE)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::~rotationalConstraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::rotationalConstraint> Foam::rotationalConstraint::New
+(
+    const word& name,
+    const sixDOFODE& sixDOFODE,
+    const dictionary& dict
+)
+{
+    const word constraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(constraintTypeType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "rotationalConstraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const sixDOFODE& sixDOFODE,"
+            "\n    const dictionary& dict,"
+            "\n)"
+        )   << "Unknown rotation constraint type: " << constraintType
+            << endl << endl
+            << "Valid rotation constraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<rotationalConstraint>
+    (
+        cstrIter()
+        (
+            name,
+            sixDOFODE,
+            dict
+        )
+    );
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/constraints/rotationalConstraints/rotationalConstraint/rotationConstraint.H

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::rotationalConstraint
+
+Description
+    Abstract base class containing interface for rotation constraints used
+    within sixDOFODE classes.
+
+    The constraint is implicitly defined as:
+
+    g(omegaDot, t) = f(t)*omegaDot + a(t) = 0,
+
+    where omegaDot is angular 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
+                              omegaDot), 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
+    rotationalConstraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rotationalConstraint_H
+#define rotationalConstraint_H
+
+#include "sixDOFODE.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class rotationalConstraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotationalConstraint
+{
+    // Private data
+
+        //- Name of the constraint
+        const word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOFODE_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        rotationalConstraint(const rotationalConstraint&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const rotationalConstraint&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("rotationalConstraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            rotationalConstraint,
+            word,
+            (
+                const word& name,
+                const sixDOFODE& sixDOFODE,
+                const dictionary& dict
+            ),
+            (name, sixDOFODE, dict)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of rotationalConstraint
+        class iNew
+        {
+            const sixDOFODE& sixDOFODE_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOFODE)
+            :
+                sixDOFODE_(sixDOFODE)
+            {}
+
+            autoPtr<rotationalConstraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return rotationalConstraint::New(name, sixDOFODE_, dict);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        rotationalConstraint
+        (
+            const word& name,
+            const sixDOFODE& sixDOFODE,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected rotationalConstraint
+        static autoPtr<rotationalConstraint> New
+        (
+            const word& name,
+            const sixDOFODE& sixDOFODE,
+            const dictionary& dict
+        );
+
+
+    // Destructor
+
+        virtual ~rotationalConstraint();
+
+
+    // 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)
+            const vector matrixContribution() const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            const vector sourceContribution() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -27,16 +27,17 @@ Class
 Description
     6-DOF solver using a geometric method for integration of rotations.
 
+    Run-time selectable constraints are handled via Lagrangian multipliers using
+    the interface from sixDOFConstraint class.
+
 Author
     Viktor Pandza, FSB Zagreb.  All rights reserved.
     Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
-SourceFiles
-    geometricSixDOF.C
-
 \*---------------------------------------------------------------------------*/
 
 #include "geometricSixDOF.H"
+#include "scalarSquareMatrix.H"
 #include "OutputControlDictionary.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -74,19 +75,65 @@ Foam::dimensionedVector Foam::geometricSixDOF::A
     const tensor& R
 ) const
 {
-    // Fix the total force in global coordinate system
-    dimensionedVector fAbs =
-        // External force
-        force()
-        // Spring force in global coordinate system
-      - (linSpringCoeffs() & xR)
-        // Damping force in global coordinate system
-      - (linDampingCoeffs() & uR);
+    // Create a scalar square matrix representing Newton equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(translationalConstraints_.size() + 3, 0.0);
+    scalarSquareMatrix M(rhs.size(), 0.0);
 
-    // Constrain translation simply by setting the total force to zero
-    constrainTranslation(fAbs.value());
+    // Insert mass and explicit forcing into the system. Note: translations are
+    // solved in the global coordinate system
+    const dimensionedVector explicitForcing
+    (
+        force() // External force
+      - (linSpringCoeffs() & xR) // Spring force
+      - (linDampingCoeffs() & uR); // Damping force
+    );
+    const vector& efVal = explicitForcing.value();
+    const scalar& m = mass().value();
 
-    return fAbs/mass();
+    forAll(efVal, dirI)
+    {
+        M[dirI][dirI] = m;
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(translationalConstraints_, tcI)
+    {
+        // Get reference to current constraint
+        const translationalConstraint& curTc = translationalConstraint_[tcI];
+
+        // Get matrix contribution from constraint
+        const vector mc = curTc.matrixContribution();
+
+        // Get matrix index
+        const label index = tcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(lower, dirI)
+        {
+            M[dirI][index] = mc[dirI];
+            M[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] = curTc.sourceContribution();
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains accelerations in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(M, rhs);
+
+    return
+        dimensionedVector
+        (
+            "A",
+            force().dimensions()/mass().dimensions(),
+            vector(rhs.x(), rhs.y(), rhs.z())
+        );
 }
 
 
@@ -96,19 +143,62 @@ Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
     const dimensionedVector& omega
 ) const
 {
-    // External moment (torque) in local coordinate system
-    dimensionedVector mRel =
-        // External moment
-        (dimensionedTensor("R", dimless, R) & moment());
+    // Create a scalar square matrix representing Euler equations with
+    // constraints and the corresponding source (right hand side vector).
+    // Note: size of the matrix is 3 + number of constraints
+    scalarField rhs(rotationalConstraints_.size() + 3, 0.0);
+    scalarSquareMatrix J(rhs.size(), 0.0);
 
-    // Note: constraints not implemented at the moment. They shall be
-    // implemented in terms of Lagrange multipliers.
+    // Insert moment of inertia and explicit forcing into the system
+    const dimensionedVector explicitForcing
+    (
+        E(omega) // Euler part
+      + (dimensionedTensor("R", dimless, R) & moment()) // External torque
+    );
+    const vector& efVal = explicitForcing.value();
+    const diagTensor& I = momentOfInertia().value();
+
+    forAll(efVal, dirI)
+    {
+        J[dirI][dirI] = I[dirI];
+
+        rhs[dirI] = efVal[dirI];
+    }
+
+    // Insert contributions from the constraints
+    forAll(rotationalConstraints_, rcI)
+    {
+        // Get reference to current constraint
+        const rotationalConstraint& curRc = rotationalConstraints_[rcI];
+
+        // Get matrix contribution from the constraint
+        const vector mc = curRc.matrixContribution();
+
+        // Get matrix index
+        const label index = rcI + 3;
+
+        // Insert contributions into the matrix
+        forAll(upper, dirI)
+        {
+            J[dirI][index] = mc[dirI];
+            J[index][dirI] = mc[dirI];
+        }
+
+        // Insert source contribution (remainder of the constraint function)
+        rhs[index] = curRc.sourceContribution();
+    }
+
+    // Solve the matrix using LU decomposition. Note: solution is in the rhs and
+    // it contains OmegaDot's in the first three entries and corresponding
+    // Lagrangian multipliers in other entries.
+    scalarSquareMatrix::LUsolve(J, rhs);
 
     return
-        inv(momentOfInertia())
-      & (
-            E(omega)
-          + mRel
+        dimensionedVector
+        (
+            "OmegaDot",
+            moment().dimensions()/momentOfInertia().dimensions(),
+            vector(rhs.x(), rhs.y(), rhs.z())
         );
 }
 
@@ -122,46 +212,6 @@ Foam::dimensionedVector Foam::geometricSixDOF::E
 }
 
 
-void Foam::geometricSixDOF::constrainTranslation(vector& vec) const
-{
-    // Constrain the vector with respect to referent or global coordinate system
-    if (referentMotionConstraints_)
-    {
-        vector consVec(referentRotation_.R() & vec);
-
-        if (fixedSurge_)
-        {
-            consVec.x() = 0;
-        }
-        if (fixedSway_)
-        {
-            consVec.y() = 0;
-        }
-        if (fixedHeave_)
-        {
-            consVec.z() = 0;
-        }
-
-        vec = referentRotation_.invR() & consVec;
-    }
-    else
-    {
-        if (fixedSurge_)
-        {
-            vec.x() = 0;
-        }
-        if (fixedSway_)
-        {
-            vec.y() = 0;
-        }
-        if (fixedHeave_)
-        {
-            vec.z() = 0;
-        }
-    }
-}
-
-
 Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
 {
     tensor R;
@@ -257,14 +307,8 @@ void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
     // HJ, 23/Mar/2015
     coeffs_ = gsd.coeffs_;
 
-    fixedSurge_ = gsd.fixedSurge_;
-    fixedSway_ = gsd.fixedSway_;
-    fixedHeave_ = gsd.fixedHeave_;
-    fixedRoll_ = gsd.fixedRoll_;
-    fixedPitch_ = gsd.fixedPitch_;
-    fixedYaw_ = gsd.fixedYaw_;
-    referentMotionConstraints_ = gsd.referentMotionConstraints_;
-    referentRotation_ = gsd.referentRotation_;
+    translationalConstraints = gsd.translationalConstraints_;
+    rotationalConstraints = gsd.rotationalConstraints_;
 }
 
 
@@ -285,33 +329,11 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
     omega_(dict().lookup("omega")),
     omegaAverage_("omegaAverage", omega_),
 
-    nEqns_(),
-    coeffs_(),
+    coeffs_(12, 0.0),
 
-    fixedSurge_(dict().lookup("fixedSurge")),
-    fixedSway_(dict().lookup("fixedSway")),
-    fixedHeave_(dict().lookup("fixedHeave")),
-    fixedRoll_(dict().lookup("fixedRoll")),
-    fixedPitch_(dict().lookup("fixedPitch")),
-    fixedYaw_(dict().lookup("fixedYaw")),
-    referentMotionConstraints_
-    (
-        dict().lookupOrDefault<Switch>
-        (
-            "referentMotionConstraints",
-            false
-        )
-    ),
-    referentRotation_(vector(1, 0, 0), 0)
+    translationalConstraints_(),
+    rotationalConstraints_()
 {
-    // Missing constraints. Count how many rotational constraints we have in
-    // order to count the number of equations.
-
-    nEqns_ = 12;
-
-    // Set size for ODE coefficients depending on number of equations
-    coeffs_.setSize(nEqns_);
-
     // Set ODE coefficients from position and rotation
 
     // Linear displacement relative to spring equilibrium
@@ -336,6 +358,31 @@ Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
     coeffs_[9] = 0;
     coeffs_[10] = 0;
     coeffs_[11] = 0;
+
+
+    // Read and construct constraints
+
+    // Read translation constraints if they are present
+    if (dict().found("translationalConstraints"))
+    {
+        PtrList<translationalConstraints> tcList
+        (
+            dict().lookup("translationalConstraints"),
+            translationalConstraints::iNew(*this)
+        );
+        translationalConstraints_.transfer(tcList);
+    }
+
+    // Read rotation constraints if they are present
+    if (dict().found("rotationalConstraints"))
+    {
+        PtrList<rotationalConstraints> tcList
+        (
+            dict().lookup("rotationalConstraints"),
+            rotationalConstraints::iNew(*this)
+        );
+        rotationalConstraints_.transfer(tcList);
+    }
 }
 
 
@@ -354,17 +401,10 @@ Foam::geometricSixDOF::geometricSixDOF
     omega_(gsd.omega_.name(), gsd.omega_),
     omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
 
-    nEqns_(gsd.nEqns_),
     coeffs_(gsd.coeffs_),
 
-    fixedSurge_(gsd.fixedSurge_),
-    fixedSway_(gsd.fixedSway_),
-    fixedHeave_(gsd.fixedHeave_),
-    fixedRoll_(gsd.fixedRoll_),
-    fixedPitch_(gsd.fixedPitch_),
-    fixedYaw_(gsd.fixedYaw_),
-    referentMotionConstraints_(gsd.referentMotionConstraints_),
-    referentRotation_(gsd.referentRotation_)
+    translationalConstraints_(gsd.translationalConstraints_),
+    rotationalConstraints_(gsd.rotationalConstraints_),
 {}
 
 
@@ -581,18 +621,10 @@ bool Foam::geometricSixDOF::writeData(Ostream& os) const
     os.writeKeyword("omega") << tab << omega_
         << token::END_STATEMENT << nl << nl;
 
-    os.writeKeyword("fixedSurge") << tab << fixedSurge_ <<
-        token::END_STATEMENT << nl;
-    os.writeKeyword("fixedSway") << tab << fixedSway_ <<
-        token::END_STATEMENT << nl;
-    os.writeKeyword("fixedHeave") << tab << fixedHeave_ <<
-        token::END_STATEMENT << nl;
-    os.writeKeyword("fixedRoll") << tab << fixedRoll_ <<
-        token::END_STATEMENT << nl;
-    os.writeKeyword("fixedPitch") << tab << fixedPitch_ <<
-        token::END_STATEMENT << nl;
-    os.writeKeyword("fixedYaw") << tab << fixedYaw_ <<
-        token::END_STATEMENT << nl << endl;
+    os.writeKeyword("translationalConstraints") << tab
+        << translationalConstraints_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationalConstraints") << tab
+        << rotationalConstraints_ << token::END_STATEMENT << nl << endl;
 
     return os.good();
 }

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H

@@ -27,6 +27,9 @@ Class
 Description
     6-DOF solver using a geometric method for integration of rotations.
 
+    Run-time selectable constraints are handled via Lagrangian multipliers using
+    the interface from sixDOFConstraint class.
+
     Reference (bibtex):
 
     @article {mullerTerze2016,
@@ -58,6 +61,8 @@ SourceFiles
 #include "sixDOFODE.H"
 #include "finiteRotation.H"
 #include "tolerancesSwitch.H"
+#include "translationalConstraint.H"
+#include "rotationalConstraint.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -101,40 +106,17 @@ class geometricSixDOF
 
         // ODE controls
 
-            //- Number of equations (depending on rotational constraints)
-            scalar nEqns_;
-
             //- ODE coefficients
             scalarField coeffs_;
 
 
-        // Motion constraints (given as fixed motion components)
+        // Motion constraints
 
-            //- Fixed surge (x-translation)
-            Switch fixedSurge_;
+            //- List of translational constraints
+            PtrList<translationalConstraint> translationalConstraints_;
 
-            //- Fixed sway (y-translation)
-            Switch fixedSway_;
-
-            //- Fixed heave (z-translation)
-            Switch fixedHeave_;
-
-            //- Fixed roll (rotation around x)
-            Switch fixedRoll_;
-
-            //- Fixed pitch (rotation around y)
-            Switch fixedPitch_;
-
-            //- Fixed yaw (rotation around z)
-            Switch fixedYaw_;
-
-            //- Constraints in referent coordinate system
-            Switch referentMotionConstraints_;
-
-            //- Rotation of referent coordinate system. Currently defined using
-            //  quaternions for comatibility. Need to reformulate.
-            //  VV, 28/Feb/2017.
-            HamiltonRodriguezRot referentRotation_;
+            //- List of rotational constraints
+            Ptrlist<rotationalConstraint> rotationalConstraints_;
 
 
     // Private Member Functions
@@ -172,10 +154,6 @@ class geometricSixDOF
                 const dimensionedVector& omega
             ) const;
 
-            //- Constrain translation vector in referent or global coordinate
-            //  system
-            void constrainTranslation(vector& vec) const;
-
             //- Exponential map used to calculate increment of the rotation
             //  tensor
             tensor expMap(const vector& rotInc) const;
@@ -235,6 +213,23 @@ public:
 
     // Member Functions
 
+        // Access member functions
+
+            //- Return const reference to translation constraints
+            const PtrList<translationalConstraint>&
+            translationalConstraints() const
+            {
+                return translationalConstraints_;
+            }
+
+            //- Return const reference to translation constraints
+            const PtrList<rotationalConstraint>&
+            rotationalConstraints() const
+            {
+                return rotationalConstraints_;
+            }
+
+
         // Virtual interface for 6DOF motion state
 
             // Variables in relative coordinate system
@@ -284,7 +279,7 @@ public:
                 //- Return number of equations
                 virtual label nEqns() const
                 {
-                    return nEqns_;
+                    return 12;
                 }
 
                 //- Return access to coefficients

src/foam/matrices/Matrix/Matrix.H

@@ -112,7 +112,7 @@ public:
         Matrix(const label n, const label m, const Type&);
 
         //- Copy constructor.
-    Matrix(const Matrix<Form, Type>&);
+        Matrix(const Matrix<Form, Type>&);
 
         //- Construct from Istream.
         Matrix(Istream&);