Merge commit '9fbf8f1bfce0b178d61be485fac02007bf9c8b77' into geometricImmersedBoundary

.gitattributes

@@ -0,0 +1,3 @@
+
+*.H gitlab-language=cpp
+*.C gitlab-language=cpp

ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch

@@ -1,6 +1,31 @@
+diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile
+--- ParMGridGen-1.0_orig/Makefile	2017-04-04 15:02:44.020713666 +0200
++++ ParMGridGen-1.0/Makefile	2017-04-04 15:21:48.582647336 +0200
+@@ -1,16 +1,21 @@
+ default:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 
+ serial:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 
+ parallel:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 	(cd ParMGridGen ; make)
+ clean:
++	(mkdir bin)
+ 	(cd MGridGen ; make clean)
+ 	(cd ParMGridGen ; make clean )
+ 
+ realclean:
++	(mkdir bin)
+ 	(cd MGridGen ; make realclean )
+ 	(cd ParMGridGen ; make realclean )
 diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
---- ParMGridGen-1.0_orig/Makefile.in	2001-12-04 16:30:33.000000000 -0800
+--- ParMGridGen-1.0_orig/Makefile.in	2017-04-04 15:02:44.012713543 +0200
-+++ ParMGridGen-1.0/Makefile.in	2013-08-22 20:07:33.491171127 -0700
++++ ParMGridGen-1.0/Makefile.in	2017-04-04 15:06:00.159742074 +0200
 @@ -1,6 +1,6 @@
  #--------------------------------------------------------------------------
  # Which make to use 
@@ -18,6 +43,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
  
  # Which loader to use
  LD = cc
+@@ -22,7 +22,7 @@
+ LDOPTIONS = -O3
+ 
+ # Where to put the executable 
+-BINDIR = ../..
++BINDIR = ../../bin
+ 
+ # Additional libraries
+ DMALLOCDIR = /usr/local
 @@ -33,22 +33,25 @@
  
  # In which directories to look for any additional libraries

ThirdParty/rpmBuild/SOURCES/ParMGridGen-1.0.patch_darwin

@@ -79,8 +79,8 @@ diff -ruN ParMGridGen-1.0_orig/MGridGen/Programs/Makefile ParMGridGen-1.0/MGridG
  ifeq ($(ddmalloc),yes)
    DEBUGFLAGS := $(DEBUGFLAGS) -DDMALLOC -DDEBUG
 diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
---- ParMGridGen-1.0_orig/Makefile.in	2011-12-24 13:54:44.000000000 -0500
+--- ParMGridGen-1.0_orig/Makefile.in	2001-12-05 01:30:33.000000000 +0100
-+++ ParMGridGen-1.0/Makefile.in	2011-12-24 13:49:26.000000000 -0500
++++ ParMGridGen-1.0/Makefile.in	2017-04-04 15:36:04.695980033 +0200
 @@ -1,6 +1,6 @@
  #--------------------------------------------------------------------------
  # Which make to use 
@@ -98,6 +98,15 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
  
  # Which loader to use
  LD = cc
+@@ -22,7 +22,7 @@
+ LDOPTIONS = -O3
+ 
+ # Where to put the executable 
+-BINDIR = ../..
++BINDIR = ../../bin
+ 
+ # Additional libraries
+ DMALLOCDIR = /usr/local
 @@ -33,22 +33,25 @@
  
  # In which directories to look for any additional libraries
@@ -129,3 +138,28 @@ diff -ruN ParMGridGen-1.0_orig/Makefile.in ParMGridGen-1.0/Makefile.in
  
  #--------------------------------------------------------------------------
  # 
+diff -ruN ParMGridGen-1.0_orig/Makefile ParMGridGen-1.0/Makefile
+--- ParMGridGen-1.0_orig/Makefile	2001-11-09 00:41:22.000000000 +0100
++++ ParMGridGen-1.0/Makefile	2017-04-04 14:51:04.033914737 +0200
+@@ -1,16 +1,21 @@
+ default:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 
+ serial:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 
+ parallel:
++	(mkdir bin)
+ 	(cd MGridGen ; make)
+ 	(cd ParMGridGen ; make)
+ clean:
++	(mkdir bin)
+ 	(cd MGridGen ; make clean)
+ 	(cd ParMGridGen ; make clean )
+ 
+ realclean:
++	(mkdir bin)
+ 	(cd MGridGen ; make realclean )
+ 	(cd ParMGridGen ; make realclean )

ThirdParty/rpmBuild/SOURCES/ParaView-4.4.0.patch_darwin

@@ -0,0 +1,12 @@
+diff -ruN ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake
+--- ParaView-4.3.1_orig/Applications/ParaView-4.3.1_extra_install_Darwin.cmake	1969-12-31 19:00:00.000000000 -0500
++++ ParaView-4.3.1/Applications/ParaView-4.3.1_extra_install_Darwin.cmake	2013-10-02 19:00:00.000000000 -0400
+@@ -0,0 +1,8 @@
++#
++# Additional install rules for Mac OS X platforms 
++#
++INSTALL (DIRECTORY buildObj/bin/paraview.app
++           DESTINATION ${CMAKE_INSTALL_PREFIX}/bin
++	   USE_SOURCE_PERMISSIONS
++           COMPONENT Runtime)
++

ThirdParty/rpmBuild/SOURCES/libccmio-2.6.1.patch_0

@@ -1,13 +1,13 @@
 diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/config.gnu.to.star
---- libccmio-2.6.1_orig/config/config.gnu.to.star	2007-12-17 18:32:03.000000000 -0500
+--- libccmio-2.6.1_orig/config/config.gnu.to.star       2007-12-18 00:32:03.000000000 +0100
-+++ libccmio-2.6.1/config/config.gnu.to.star	2010-12-18 14:50:35.000000000 -0500
++++ libccmio-2.6.1/config/config.gnu.to.star    2017-04-04 13:41:35.000000000 +0200
 @@ -1,4 +1,4 @@
 -#!/bin/sh
 +#!/bin/bash
 
  # $Id: config.gnu.to.star,v 1.4 2006/06/05 21:12:16 geoffp Exp $
 
-@@ -34,6 +34,12 @@
+@@ -34,6 +34,13 @@
      x86_64-unknown-linux-gnu-null)     echo linux64_2.4-x86-glibc_2.2.5     ;;
      ppc64-unknown-linux-gnu-null)      echo linux64_2.6-pwr4-glibc_2.3.3    ;;
      i386-apple-darwin8-null)           echo i386-apple-darwin8              ;;
@@ -17,6 +17,7 @@ diff -ruN libccmio-2.6.1_orig/config/config.gnu.to.star libccmio-2.6.1/config/co
 +    i386-apple-darwin13-null)          echo i386-apple-darwin13             ;;
 +    i386-apple-darwin14-null)          echo i386-apple-darwin14             ;;
 +    i386-apple-darwin15-null)          echo i386-apple-darwin15             ;;
++    i386-apple-darwin16-null)          echo i386-apple-darwin16             ;;
      *)                                 echo unknown                         ;;
  esac
 
@@ -30,15 +31,15 @@ diff -ruN libccmio-2.6.1_orig/config/config.guess libccmio-2.6.1/config/config.g
  #   Copyright (C) 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001
  #   Free Software Foundation, Inc.
 diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config.system
---- libccmio-2.6.1_orig/config/config.system	2008-02-25 22:07:16.000000000 -0500
+--- libccmio-2.6.1_orig/config/config.system    2008-02-26 04:07:16.000000000 +0100
-+++ libccmio-2.6.1/config/config.system	2010-12-18 14:51:34.000000000 -0500
++++ libccmio-2.6.1/config/config.system 2017-04-04 13:42:15.000000000 +0200
 @@ -1,4 +1,4 @@
 -#! /bin/sh
 +#! /bin/bash
 
  # $Id: config.system,v 1.2 2005/09/29 22:19:19 geoffp Exp $
 
-@@ -87,6 +87,27 @@
+@@ -87,6 +87,30 @@
      i386-apple-darwin8.11.1)
  	echo i386-apple-darwin8 ;;
 
@@ -62,6 +63,9 @@ diff -ruN libccmio-2.6.1_orig/config/config.system libccmio-2.6.1/config/config.
 +
 +    i386-apple-darwin15.* )
 +       echo i386-apple-darwin15 ;;
++
++    i386-apple-darwin16.* )
++       echo i386-apple-darwin16 ;;
 +
      *)
          echo unknown

ThirdParty/rpmBuild/SPECS/libccmio-2.6.1.spec

@@ -105,6 +105,7 @@ Patch0:                 libccmio-2.6.1.patch_0
     [ ! -d config/i386-apple-darwin13 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin13
     [ ! -d config/i386-apple-darwin14 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin14
     [ ! -d config/i386-apple-darwin15 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin15
+    [ ! -d config/i386-apple-darwin16 ] && cp -r config/i386-apple-darwin8 config/i386-apple-darwin16
 %endif
     # Warning:
     #  1: The name of the ADF library will be renamed to libadf_ccmio since this

ThirdParty/rpmBuild/SPECS/mesquite-2.1.2.spec

@@ -80,6 +80,7 @@ Source: 		%url/%{name}-%{version}.tar.gz
 Prefix: 		%{_prefix}
 Group: 			Development/Tools
 Patch0:                 mesquite-2.1.2_patch0
+Patch1:                 mesquite-2.1.2_patch1
 
 %define _installPrefix  %{_prefix}/packages/%{name}-%{version}/platforms/%{_WM_OPTIONS}
 
@@ -90,6 +91,7 @@ Patch0:                 mesquite-2.1.2_patch0
 %setup -q
 
 %patch0 -p1
+%patch1 -p1
 
 %build
     # export WM settings in a form that GNU configure recognizes

applications/solvers/basic/sixDOFSolver/sixDOFSolver.C

@@ -36,7 +36,7 @@ Description
 #include "objectRegistry.H"
 #include "foamTime.H"
 #include "ODESolver.H"
-#include "sixDOFbodies.H"
+#include "sixDOFBodies.H"
 #include "OFstream.H"
 
 using namespace Foam;
@@ -48,9 +48,12 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 
-    sixDOFbodies structure(runTime);
+    sixDOFBodies structure(runTime);
     OFstream of(runTime.path()/"motion.dat");
 
+    // Write header for output file
+    of << "# Time, CoG_x, CoG_y, CoG_z, omega_x, omega_y, omega_z" << endl;
+
     Info<< "\nStarting time loop\n" << endl;
 
     for (runTime++; !runTime.end(); runTime++)
@@ -73,11 +76,14 @@ int main(int argc, char *argv[])
                 << structure()[bodyI].omegaAverage().value() << nl
                 << "Current omega in time instant = "
                 << structure()[bodyI].omega().value()  << nl
-                << "Average omegaAbs in time step = "
-                << structure()[bodyI].omegaAverageAbsolute().value() << nl
                 << endl;
 
-            of  << structure()[bodyI].X().value().x() << tab;
+            of  << structure()[bodyI].X().value().x() << tab
+                << structure()[bodyI].X().value().y() << tab
+                << structure()[bodyI].X().value().z() << tab
+                << structure()[bodyI].omega().value().x() << tab
+                << structure()[bodyI].omega().value().y() << tab
+                << structure()[bodyI].omega().value().z() << tab;
         }
 
         of << endl;

applications/solvers/incompressible/RichardsFoam/Make/files

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

applications/solvers/incompressible/RichardsFoam/Make/options

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

applications/solvers/incompressible/RichardsFoam/RichardsFoam.C

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    RichardsFoam
+
+Description
+    Transient solver for flow in unsaturated porous media
+    With chord slope formulation of the Richards equation.
+    van Genuchten laws for unsaturated hydraulic properties parametrisation
+    Global computation of the convergence criterium
+    Adaptative time stepping with a stabilisation procedure
+    NB 1: use backward scheme for time discretisation
+    NB 2: use only mesh with constant cell volumes
+
+References
+   version 0.0 (develloped with OpenFOAM 2.0.1)
+   Details may be found in:
+   Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S.K., Labat, D.,
+    Pokrovsky, O.S., Sekhar, M., Ababou, R., Quintard, M., Submitted.
+   Mechanistic modelling of water fluxes at the watershed scale: An open source
+    massively parallel solver for Richards equation.
+   Submitted to Computer Physics Communications.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pimpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+//    pimpleControl pimple(mesh);
+
+#   include "readPicardControls.H"
+#   include "createFields.H"
+#   include "initContinuityErrs.H"
+#   include "createTimeControls.H"
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    // starting of the time loop.
+    while (runTime.loop())
+    {
+        // time step control operations.
+#       include "readTimeControls.H"
+#       include "setDeltaT.H"
+
+//        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // Beginning of the stabilisation loop for the stabilised adaptive time
+        // step procedure.
+        for (int cyc = 0; cyc < nMaxCycle; cyc++)
+        {
+            // Beginning of the Picard loop.
+            for (int pic = 0; pic < nIterPicard; pic++)
+            {
+#               include "psiEqn.H"
+            }
+
+            // Exit test for the loop associated with the stabilisation cycles
+            // for the adaptive time step procedure.
+            if (crit < precPicard)
+            {
+                break;
+            }
+            else
+            {
+                Info << "Criterion not reached, restart time loop iteration"
+                    << "with a smaller time step / Error = " << crit
+                    << nl << endl;
+
+                runTime.setDeltaT((1/tFact)*runTime.deltaTValue());
+
+                Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+            }
+            // End of the stabilisation cycles loop.
+        }
+
+        // Warning test in case of convergence failure of the Picard loop.
+        if (crit >= precPicard)
+        {
+            Info<< "Convergence failure / Error = " << crit << nl << endl;
+            currentPicard = nIterPicard;
+        }
+
+        // Final updating of the result fields before going to the next time
+        // iteration.
+        psi_tmp = psi;
+
+        thtil_tmp = 0.5*
+        (
+            (1 + sign(psi_tmp)) + (1 - sign(psi_tmp))*
+            pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
+        );
+
+        theta = (thetas - thetar)*thtil + thetar;
+
+        U = - Krel*((fvc::grad(psi)) + vuz);
+
+        // Writting  of the result.
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        // end of the time loop.
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/incompressible/RichardsFoam/createFields.H

@@ -0,0 +1,209 @@
+Info<< "Reading transportProperties\n" << endl;
+
+// reading of the physical constant associated with the considered problem.
+// Localisation of the data within the considered case:
+// constant/transportProperties
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+dimensionedScalar K
+(
+    transportProperties.lookup("K")
+);
+
+dimensionedScalar alpha
+(
+    transportProperties.lookup("alpha")
+);
+
+dimensionedScalar thetas
+(
+    transportProperties.lookup("thetas")
+);
+
+dimensionedScalar thetar
+(
+    transportProperties.lookup("thetar")
+);
+
+dimensionedScalar n
+(
+    transportProperties.lookup("n")
+);
+
+dimensionedScalar C
+(
+    transportProperties.lookup("C")
+);
+
+dimensionedScalar S
+(
+    transportProperties.lookup("S")
+);
+
+// declaration of the variable and results fields
+
+// Water velocity field [m/s]
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Water saturation field [-]
+Info<< "Reading field theta\n" << endl;
+volScalarField theta
+(
+    IOobject
+    (
+        "theta",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Water pressure field [m] - field of resolution.
+Info<< "Reading field psi\n" << endl;
+volScalarField psi
+(
+    IOobject
+    (
+        "psi",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Field of residuals for the Picard loop [m].
+Info<< "Reading field err\n" << endl;
+volScalarField err
+(
+    IOobject
+    (
+        "err",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Dimensionned unit scalar field [m].
+Info<< "Reading field vuz\n" << endl;
+volVectorField vuz
+(
+    IOobject
+    (
+        "vuz",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+// Dimensionless unit vertical upward vector field.
+Info<< "Reading field usf\n" << endl;
+volScalarField usf
+(
+    IOobject
+    (
+        "usf",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#   include "createPhi.H"
+
+// Initialisation of the scalar containing the number of mesh cells. Note the
+// use of gSum instead of sum.
+double nbMesh;
+nbMesh = gSum(usf);
+
+// Initialisation of the scalar containing the residual for the exit test of the
+// Picard loop.
+double crit;
+crit=0.;
+
+// Initialisation of the token which counts the number of Picard iteraion for
+// the adaptive time step procedure.
+int currentPicard;
+currentPicard = nIterPicard-3;
+
+// Initialisation of the token which counts the number of Stabilisation cycle
+// for the stabilisation of the adaptive time step procedure.
+int sc;
+sc = 0;
+
+// Initialisation of the field of altitudes.
+volVectorField positionVector = mesh.C();
+volScalarField z = positionVector.component(vector::Z);
+
+// Initialisation of the intermediate fields for the Picard loop.
+volScalarField psi_tmp = psi;
+volScalarField psim1 = psi;
+
+// Initialisation of the varying transport properties for the Picard loop.
+volScalarField thtil =
+    0.5*
+    (
+        (1 + sign(psi)) + (1 - sign(psi))*
+        pow((1 + pow(mag(alpha*psi),n)), - (1 - (1/n)))
+    );
+
+volScalarField thtil_tmp =
+    0.5*
+    (
+        (1 + sign(psi_tmp)) + (1-sign(psi_tmp))*
+        pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
+    );
+
+volScalarField Krel =
+    0.5*
+    (
+        (1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil,0.5)*
+        pow((1 - pow((1 - pow(thtil,(n/(n - 1)))),(1 - (1/n)))),2)
+    );
+
+volScalarField Crel =
+    S + 0.5*
+    (
+        (1 - sign(psi))*((thetas - thetar)*(thtil - thtil_tmp)*
+        (1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi)) + psi - psi_tmp)))
+    );
+
+// Initialisation of the gravity term.
+volVectorField gradk = fvc::grad(Krel);
+volScalarField gradkz = gradk.component(vector::Z);
+
+// Initialisation of the velocity field.
+U = - Krel*((fvc::grad(psi)) + vuz);

applications/solvers/incompressible/RichardsFoam/psiEqn.H

@@ -0,0 +1,60 @@
+psim1 = psi;
+psi = psi_tmp;
+
+// Resolution of the linear system.
+fvScalarMatrix psiEqn
+(
+    Crel*fvm::ddt(psi)
+ ==
+    fvm::laplacian(Krel, psi, "laplacian(Krel,psi)")
+  + gradkz
+);
+psiEqn.relax();
+psiEqn.solve();
+
+// Update of the varying transport properties.
+thtil = 0.5*
+(
+    (1 + sign(psi)) + (1 - sign(psi))*
+    pow((1+pow(mag(alpha*psi),n)),-(1-(1/n)))
+);
+Krel = 0.5*
+(
+    (1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil, 0.5)*
+    pow((1 - pow((1 - pow(thtil, (n/(n - 1)))),(1 - (1/n)))), 2)
+);
+Crel= S + 0.5*
+(
+    (1 - sign(psi))*
+    (
+        (thetas - thetar)*(thtil - thtil_tmp)*
+        (
+            1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi))
+            + psi - psi_tmp)
+        )
+    )
+);
+
+// Update of the gravity term.
+gradk = fvc::grad(Krel);
+gradkz = gradk.component(2);
+
+// Computation of the field of residuals for the exit test of
+// the Picard loop.
+err = psi - psim1;
+
+// Computation of the residual for the exit test of the Picard
+// loop. Note the use of gSum instead of sum.
+crit = gSumMag(err)/nbMesh;
+
+// exit test for the Picard loop.
+if (crit < precPicard)
+{
+    Info << " Erreur = " << crit
+        << " Picard = " << pic
+        << nl << endl;
+
+    currentPicard=pic;
+
+    break;
+}

applications/solvers/incompressible/RichardsFoam/readPicardControls.H

@@ -0,0 +1,25 @@
+// Reading of the controling parameters of the Picard loop and of the
+// Stabilisation cycles loop for the adaptive time step procedure.
+// Localisation of the data within the considered case: system/fvSolution
+const dictionary& picardDict = mesh.solutionDict().subDict("Picard");
+
+// Maximum number of Picard iterations.
+const int nIterPicard =
+    picardDict.lookupOrDefault<int>("nIterPicard", 1);
+
+// Maximum number of stabilisation cycles.
+const int nMaxCycle =
+    picardDict.lookupOrDefault<int>("nMaxCycle", 1);
+
+// Number of time iterations with low number of Picard iteration beyond which
+// the time step is increased.
+const int stabilisationThreshold =
+    picardDict.lookupOrDefault<int>("stabilisationThreshold", 1);
+
+// Exit criterion for the Picard loop.
+const double precPicard =
+    picardDict.lookupOrDefault<double>("precPicard", 1.);
+
+// Factor of increase/decrease of the time step.
+const double tFact =
+    picardDict.lookupOrDefault<double>("tFact", 1.);

applications/solvers/incompressible/RichardsFoam/setDeltaT.H

@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a low number of Picard iteration and a good
+    convergence.  Reduction of time-step is immediate, but increase is damped to
+    avoid unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+scalar deltaTFact = 1.;
+
+if (adjustTimeStep)
+{
+    deltaTFact = 1.;
+
+    // If There is a too high number of Picard iteration, decrease of the time
+    // step.
+    if (currentPicard > nIterPicard-2)
+    {
+        deltaTFact = 1./tFact;
+    }
+
+    if (currentPicard >= 3)
+    {
+        sc = 0;
+    }
+
+    if ((currentPicard < 3) && (sc < stabilisationThreshold))
+    {
+        sc = sc+1;
+    }
+
+    // If there is more than 'stabilisationThreshold' iteration that the number
+    // of Picard iteration is low, increase of the time step.
+    if ((currentPicard < 3) && (sc == stabilisationThreshold))
+    {
+        deltaTFact = tFact;
+        sc = 0;
+    }
+
+    // Reset of the time step if needed (time step is always lower than
+    // 'maxDeltaT').
+    if (deltaTFact != 1.)
+    {
+        runTime.setDeltaT
+        (
+            min
+            (
+                deltaTFact*runTime.deltaTValue(),
+                maxDeltaT
+            )
+        );
+    }
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //

applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/files

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

applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/options

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

applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/spatialMeanValueRichardsonFoam.C

@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    spatialMeanValue
+
+Description
+    Calculates the spatial mean of the specified volScalarField over the
+    whole domain (for regular grid).
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addRegionOption.H"
+
+    timeSelector::addOptions();
+    argList::validArgs.append("fieldName");
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    instantList timeDirs = timeSelector::select0(runTime, args);
+
+#   include "createNamedMesh.H"
+
+    word fieldName(args.additionalArgs()[0]);
+
+    Info << "Calculating avarage pressure heads:" << endl;
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+
+        IOobject fieldHeader
+        (
+            fieldName,
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ
+        );
+
+        // Check field exists
+        if (fieldHeader.headerOk())
+        {
+            mesh.readUpdate();
+
+            // Read field and calc mean, for regular grid
+            if (fieldHeader.headerClassName() == volScalarField::typeName)
+            {
+                volScalarField field(fieldHeader, mesh);
+
+                int nbMesh;
+	            nbMesh = 0;
+
+	            forAll(field, cellI)
+	            {
+	                nbMesh++;
+	            }
+
+                Info<< runTime.timeName()<< " "
+                    << sum(field).value()/nbMesh<< " "
+                    << endl;
+            }
+            else
+            {
+                FatalError
+                    << "Only possible for volScalarField "<< "s "
+                    << nl << exit(FatalError);
+            }
+        }
+        else
+        {
+            Info<< "    No field " << fieldName << endl;
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/advanced/runDynamicMesh/Make/files

@@ -0,0 +1,3 @@
+runDynamicMesh.C
+
+EXE = $(FOAM_APPBIN)/runDynamicMesh

applications/utilities/mesh/advanced/runDynamicMesh/Make/options

@@ -0,0 +1,12 @@
+EXE_INC = \
+    -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -ldynamicFvMesh \
+    -ldynamicMesh \
+    -lmeshTools \
+    -lfiniteVolume \
+    -L$(MESQUITE_LIB_DIR) -lmesquite

applications/utilities/mesh/advanced/runDynamicMesh/createFields.H

@@ -0,0 +1,13 @@
+    Info<< "Reading field alpha\n" << endl;
+    volScalarField alpha
+    (
+        IOobject
+        (
+            "alpha",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );

applications/utilities/mesh/advanced/runDynamicMesh/runDynamicMesh.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Application
+    runDynamicMesh
+
+Description
+    Creates dynamicFvMesh and updates it in a time loop. Reads indicator field
+    alpha used for dynamic mesh refinement to trigger refinement in certain
+    areas.
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "simpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createDynamicFvMesh.H"
+
+    simpleControl simple(mesh);
+
+#   include "createFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        bool meshChanged = mesh.update();
+
+        // Write only if mesh has changed
+        if (meshChanged)
+        {
+            runTime.write();
+        }
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return(0);
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/star4ToFoam/cellZoneToCellTableId/cellZoneToCellTableId.C

@@ -27,6 +27,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
+#include "objectRegistry.H"
 #include "foamTime.H"
 #include "volFields.H"
 #include "OFstream.H"

applications/utilities/mesh/conversion/star4ToFoam/foamMeshToStar/foamMeshToStar.C

@@ -35,6 +35,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
+#include "objectRegistry.H"
 #include "foamTime.H"
 #include "volFields.H"
 #include "cellModeller.H"

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/files

@@ -4,6 +4,7 @@ linearNormal/linearNormal.C
 linearRadial/linearRadial.C
 sigmaRadial/sigmaRadial.C
 wedge/wedge.C
+gradedNormal/gradedNormal.C
 
 LIB = $(FOAM_LIBBIN)/libextrudeModel
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/Make/options

@@ -1,8 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ldynamicMesh
+    -ldynamicMesh \
-
+    -lODE

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.C

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "gradedNormal.H"
+#include "addToRunTimeSelectionTable.H"
+#include "BisectionRoot.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(gradedNormal, 0);
+
+addToRunTimeSelectionTable(extrudeModel, gradedNormal, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+gradedNormal::gradedNormal(const dictionary& dict)
+:
+    extrudeModel(typeName, dict),
+    thickness_(readScalar(coeffDict_.lookup("thickness"))),
+    delta0_(readScalar(coeffDict_.lookup("initialCellLength"))),
+    expansionFactor_(1.0)
+{
+    // Sanity checks
+    if (thickness_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "thickness should be positive: " << thickness_
+            << exit(FatalError);
+    }
+
+    if (delta0_ <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "initialCellLength should be positive: " << delta0_
+            << exit(FatalError);
+    }
+
+    const scalar maxExpFactor =
+        coeffDict_.lookupOrDefault<scalar>("maxExpansionFactor", 3.0);
+
+    if (maxExpFactor <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "maxExpansionFactor should be positive: " << maxExpFactor
+            << exit(FatalError);
+    }
+
+    const scalar bisectionTol =
+        coeffDict_.lookupOrDefault<scalar>("bisectionTol", 1e-5);
+
+    if (bisectionTol <= SMALL)
+    {
+        FatalErrorIn("gradedNormal(const dictionary&)")
+            << "bisectionTolerance should be positive: " << bisectionTol
+            << exit(FatalError);
+    }
+
+    // Create expansion factor equation represented as a function object
+    expansionFactorEqn eqn(*this);
+
+    // Calculate the expansionFactor using the bisection algorithm with the
+    // default tolerance of 1e-5
+    BisectionRoot<expansionFactorEqn> rootFinder
+    (
+        eqn,
+        bisectionTol
+    );
+
+    // Search for the root in [0, 3], where upper bound 3 is default
+    expansionFactor_ = rootFinder.root
+    (
+        0.0,
+        maxExpFactor
+    );
+
+    // Report the result
+    Info<< "Calculated expansion factor: " << expansionFactor_ << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gradedNormal::~gradedNormal()
+{}
+
+
+// * * * * * * * * * * * * * * * * Operators * * * * * * * * * * * * * * * * //
+
+point gradedNormal::operator()
+(
+    const point& surfacePoint,
+    const vector& surfaceNormal,
+    const label layer
+) const
+{
+    scalar d = 0.0;
+
+    for (label i = 0; i < layer; ++i)
+    {
+        d += delta0_*pow(expansionFactor_, i);
+    }
+
+    return surfacePoint + d*surfaceNormal;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// ************************************************************************* //
+

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/gradedNormal/gradedNormal.H

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::extrudeModels::gradedNormal
+
+Description
+    Extrudes by transforming points normal to the surface. Input parameters are:
+    1. Extrusion thickness (in meters),
+    2. Initial delta (in meters),
+    3. Number of layers.
+
+    Expansion factor is calculated numerically using bisection algorithm.
+
+Author
+    Vuko Vukcevic. FMENA Zagreb.  All rights reserved.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gradedNormal_H
+#define gradedNormal_H
+
+#include "point.H"
+#include "extrudeModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace extrudeModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class gradedNormal Declaration
+\*---------------------------------------------------------------------------*/
+
+class gradedNormal
+:
+    public extrudeModel
+{
+    // Private data
+
+        //- Layer thickness
+        scalar thickness_;
+
+        //- Initial delta (cell size at the extruded patch)
+        scalar delta0_;
+
+        //- Expansion factor
+        scalar expansionFactor_;
+
+
+    // Expansion factor equation (functor class used by BisectionRoot)
+    class expansionFactorEqn
+    {
+        // Private data
+
+            //- Const reference to underlying gradedNormal model
+            const gradedNormal& gnm_;
+
+
+    public:
+
+        //- Construct given gradedNormal model
+        explicit expansionFactorEqn(const gradedNormal& gnm)
+        :
+            gnm_(gnm)
+        {}
+
+
+        //- Function call operator
+        scalar operator()(const scalar& x) const
+        {
+            scalar result = gnm_.thickness();
+
+            for (label i = 0; i <= gnm_.nLayers(); ++i)
+            {
+                result -= gnm_.delta0()*pow(x, i);
+            }
+
+            return result;
+        }
+    };
+
+
+public:
+
+    //- Runtime type information
+    TypeName("gradedNormal");
+
+    // Constructors
+
+        //- Construct from components
+        gradedNormal(const dictionary& dict);
+
+
+    //- Destructor
+    ~gradedNormal();
+
+
+    // Access functions
+
+        //- Return const reference to thickness
+        const scalar& thickness() const
+        {
+            return thickness_;
+        }
+
+
+        //- Return const reference to initial delta
+        const scalar& delta0() const
+        {
+            return delta0_;
+        }
+
+
+    // Member Operators
+
+        point operator()
+        (
+            const point& surfacePoint,
+            const vector& surfaceNormal,
+            const label layer
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace extrudeModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -894,7 +894,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nEnd\n" << endl;
 
-    return 0;
+    return !ok;
 }
 
 

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -660,7 +660,8 @@ autoPtr<mapPolyMesh> createRegionMesh
     const EdgeMap<label>& interfaceToPatch,
     const fvMesh& mesh,
     const label regionI,
-    const word& regionName
+    const word& regionName,
+    const word& newMeshInstance
 )
 {
     // Neighbour cellRegion.
@@ -758,7 +759,7 @@ autoPtr<mapPolyMesh> createRegionMesh
         IOobject
         (
             regionName,
-            mesh.time().timeName(),
+            newMeshInstance,
             mesh.time(),
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
@@ -977,7 +978,8 @@ void createAndWriteRegion
         interfaceToPatch,
         mesh,
         regionI,
-        regionNames[regionI]
+        regionNames[regionI],
+        newMeshInstance
     );
 
     // Read in mesh as fvMesh for mapping.  HJ, 19/May/2011

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -284,6 +284,7 @@ int main(int argc, char *argv[])
                 (
                     "cellDist",
                     runTime.timeName(),
+                    mesh.dbDir(),
                     mesh,
                     IOobject::NO_READ,
                     IOobject::NO_WRITE
@@ -914,7 +915,7 @@ int main(int argc, char *argv[])
         "faBoundary",
         mesh.time().findInstance
         (
-            mesh.dbDir()/polyMesh::meshSubDir,
+            mesh.meshDir(),
             "boundary"
         ),
         faMesh::meshSubDir,
@@ -989,12 +990,20 @@ int main(int argc, char *argv[])
 
             faMesh procMesh(procFvMesh);
 
+            word faceLabelsInstance =
+                procMesh.time().findInstance
+                (
+                    procMesh.meshDir(),
+                    "faceLabels"
+                );
+
+
             labelIOList faceProcAddressing
             (
                 IOobject
                 (
                     "faceProcAddressing",
-                    "constant",
+                    faceLabelsInstance,
                     procMesh.meshSubDir,
                     procFvMesh,
                     IOobject::MUST_READ,
@@ -1007,7 +1016,7 @@ int main(int argc, char *argv[])
                 IOobject
                 (
                     "boundaryProcAddressing",
-                    "constant",
+                    faceLabelsInstance,
                     procMesh.meshSubDir,
                     procFvMesh,
                     IOobject::MUST_READ,
@@ -1031,7 +1040,7 @@ int main(int argc, char *argv[])
                     IOobject
                     (
                         "edgeProcAddressing",
-                        "constant",
+                        faceLabelsInstance,
                         procMesh.meshSubDir,
                         procFvMesh,
                         IOobject::MUST_READ,

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -103,6 +103,37 @@ bool domainDecomposition::writeDecomposition()
     label maxProcFaces = 0;
 
 
+    // Note: get cellLevel and pointLevel. Avoid checking whether they exist or
+    // not by hand. If they don't exist, simpy assume that the level is 0
+    const labelIOList globalCellLevel
+    (
+        IOobject
+        (
+            "cellLevel",
+            this->facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        labelList(nCells(), 0)
+    );
+
+    const labelIOList globalPointLevel
+    (
+        IOobject
+        (
+            "pointLevel",
+            this->facesInstance(),
+            polyMesh::meshSubDir,
+            *this,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        labelList(nPoints(), 0)
+    );
+
+
     // Write out the meshes
     for (label procI = 0; procI < nProcs_; procI++)
     {
@@ -615,6 +646,39 @@ bool domainDecomposition::writeDecomposition()
             procBoundaryAddressing_[procI]
         );
         boundaryProcAddressing.write();
+
+        // Create and write cellLevel and pointLevel information
+        const unallocLabelList& cellMap = cellProcAddressing;
+        labelIOList procCellLevel
+        (
+            IOobject
+            (
+                "cellLevel",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(globalCellLevel, cellMap)
+        );
+        procCellLevel.write();
+
+        const unallocLabelList& pointMap = pointProcAddressing;
+        labelIOList procPointLevel
+        (
+            IOobject
+            (
+                "pointLevel",
+                procMesh.facesInstance(),
+                procMesh.meshSubDir,
+                procMesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            labelList(globalPointLevel, pointMap)
+        );
+        procPointLevel.write();
     }
 
     Info<< nl

applications/utilities/parallelProcessing/decomposePar/faMeshDecomposition.C

@@ -78,7 +78,7 @@ void faMeshDecomposition::distributeFaces()
                     IOobject
                     (
                         "faceProcAddressing",
-                        "constant",
+                        procMesh.facesInstance(),
                         procMesh.meshSubDir,
                         procMesh,
                         IOobject::MUST_READ,
@@ -94,7 +94,7 @@ void faMeshDecomposition::distributeFaces()
                     IOobject
                     (
                         "owner",
-                        "constant",
+                        procMesh.facesInstance(),
                         procMesh.meshSubDir,
                         procMesh,
                         IOobject::MUST_READ,
@@ -127,7 +127,7 @@ void faMeshDecomposition::distributeFaces()
                     IOobject
                     (
                         "faceProcAddressing",
-                        "constant",
+                        procMesh.facesInstance(),
                         procMesh.meshSubDir,
                         procMesh,
                         IOobject::MUST_READ,
@@ -263,7 +263,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
             IOobject
             (
                 "pointProcAddressing",
-                "constant",
+                procFvMesh.facesInstance(),
                 procFvMesh.meshSubDir,
                 procFvMesh,
                 IOobject::MUST_READ,
@@ -279,7 +279,7 @@ void faMeshDecomposition::decomposeMesh(const bool filterEmptyPatches)
                 IOobject
                 (
                     "faceProcAddressing",
-                    "constant",
+                    procFvMesh.facesInstance(),
                     procFvMesh.meshSubDir,
                     procFvMesh,
                     IOobject::MUST_READ,
@@ -1202,7 +1202,7 @@ bool faMeshDecomposition::writeDecomposition()
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                procFvMesh.facesInstance(),
                 procFvMesh.meshSubDir,
                 procFvMesh,
                 IOobject::MUST_READ,
@@ -1340,13 +1340,20 @@ bool faMeshDecomposition::writeDecomposition()
         maxProcPatches = max(maxProcPatches, nProcPatches);
         maxProcEdges = max(maxProcEdges, nProcEdges);
 
+        word faceLabelsInstance =
+            procMesh.time().findInstance
+            (
+                procMesh.meshDir(),
+                "faceLabels"
+            );
+
         // create and write the addressing information
         labelIOList pointProcAddressing
         (
             IOobject
             (
                 "pointProcAddressing",
-                "constant",
+                faceLabelsInstance,
                 procMesh.meshSubDir,
                 procFvMesh,
                 IOobject::NO_READ,
@@ -1361,7 +1368,7 @@ bool faMeshDecomposition::writeDecomposition()
             IOobject
             (
                 "edgeProcAddressing",
-                "constant",
+                faceLabelsInstance,
                 procMesh.meshSubDir,
                 procFvMesh,
                 IOobject::NO_READ,
@@ -1376,7 +1383,7 @@ bool faMeshDecomposition::writeDecomposition()
             IOobject
             (
                 "faceProcAddressing",
-                "constant",
+                faceLabelsInstance,
                 procMesh.meshSubDir,
                 procFvMesh,
                 IOobject::NO_READ,
@@ -1391,7 +1398,7 @@ bool faMeshDecomposition::writeDecomposition()
             IOobject
             (
                 "boundaryProcAddressing",
-                "constant",
+                faceLabelsInstance,
                 procMesh.meshSubDir,
                 procFvMesh,
                 IOobject::NO_READ,

applications/utilities/parallelProcessing/decomposeSets/processorMeshes.C

@@ -157,6 +157,16 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
         // Check if any new meshes need to be read.
         polyMesh::readUpdateState procStat = meshes_[procI].readUpdate();
 
+        // Force changed mesh if the processor mesh holds points. This is
+        // essential for changed region meshes and if the current time is
+        // equal to the initial time during construction.
+        // Pascal Beckstein, 14/Jul/2016
+        fileName timePath = meshes_[procI].time().timePath();
+        if (isFile(timePath/meshes_[procI].meshDir()/"points"))
+        {
+            procStat = polyMesh::POINTS_MOVED;
+        }
+
         // Combine into overall mesh change status
         if (stat == polyMesh::UNCHANGED)
         {
@@ -196,8 +206,36 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
 
 void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
 {
+    // We may not use
+    //   mesh.movePoints(newPoints);
+    //   mesh.write();
+    // here, as this will fail for mesh regions containing
+    //   directMappedPolyPatches or
+    //   directMappedWallPolyPatches
+    // and probably other region-coupled patches. E.g. for both above
+    // mentioned directMapped patch types we have calls to
+    //   directMappedPatchBase::map()
+    // in the corresponding initMovePoints(const pointField& p) member
+    // function. This however involves parallel communication and it would
+    // require all other connected regions to be already reconstructed and
+    // registered during reconstruction of the current region!
+    // Pascal Beckstein, 14/Jul/2016
+
     // Create the new points
-    vectorField newPoints(mesh.nPoints());
+    pointIOField newPoints
+    (
+        IOobject
+        (
+            "points",
+            mesh.time().timeName(),
+            mesh.meshSubDir,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh.allPoints()
+    );
 
     forAll (meshes_, procI)
     {
@@ -225,8 +263,8 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
         }
     }
 
-    mesh.movePoints(newPoints);
+    newPoints.write();
-    mesh.write();
+    mesh.readUpdate();
 }
 
 

applications/utilities/parallelProcessing/reconstructPar/processorMeshes.C

@@ -157,6 +157,16 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
         // Check if any new meshes need to be read.
         polyMesh::readUpdateState procStat = meshes_[procI].readUpdate();
 
+        // Force changed mesh if the processor mesh holds points. This is
+        // essential for changed region meshes and if the current time is
+        // equal to the initial time during construction.
+        // Pascal Beckstein, 14/Jul/2016
+        fileName timePath = meshes_[procI].time().timePath();
+        if (isFile(timePath/meshes_[procI].meshDir()/"points"))
+        {
+            procStat = polyMesh::POINTS_MOVED;
+        }
+
         // Combine into overall mesh change status
         if (stat == polyMesh::UNCHANGED)
         {
@@ -196,8 +206,36 @@ Foam::polyMesh::readUpdateState Foam::processorMeshes::readUpdate()
 
 void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
 {
+    // We may not use
+    //   mesh.movePoints(newPoints);
+    //   mesh.write();
+    // here, as this will fail for mesh regions containing
+    //   directMappedPolyPatches or
+    //   directMappedWallPolyPatches
+    // and probably other region-coupled patches. E.g. for both above
+    // mentioned directMapped patch types we have calls to
+    //   directMappedPatchBase::map()
+    // in the corresponding initMovePoints(const pointField& p) member
+    // function. This however involves parallel communication and it would
+    // require all other connected regions to be already reconstructed and
+    // registered during reconstruction of the current region!
+    // Pascal Beckstein, 14/Jul/2016
+
     // Create the new points
-    vectorField newPoints(mesh.nPoints());
+    pointIOField newPoints
+    (
+        IOobject
+        (
+            "points",
+            mesh.time().timeName(),
+            mesh.meshSubDir,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        mesh.allPoints()
+    );
 
     forAll (meshes_, procI)
     {
@@ -225,8 +263,8 @@ void Foam::processorMeshes::reconstructPoints(fvMesh& mesh)
         }
     }
 
-    mesh.movePoints(newPoints);
+    newPoints.write();
-    mesh.write();
+    mesh.readUpdate();
 }
 
 

applications/utilities/postProcessing/dataConversion/foamToTecplot360/Allwmake

@@ -5,10 +5,16 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
+# foamToTecplot360 will not compile on MS Windows
+if [ "$WM_OSTYPE" = "MSWindows" ]
+then
+    echo MSWindows detected. Skipping compilation of foamToTecplot360.
+else
     # build tecio
     wmake libso tecio/tecsrc
 
     # build converter
     wmake
+fi
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/postProcessing/graphics/PV4FoamReader/Allwclean

@@ -3,6 +3,6 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 rm -rf PV4FoamReader/Make
-wclean libso vtkPV4Foam
+wclean vtkPV4Foam
 
 # ----------------------------------------------------------------- end-of-file

applications/utilities/surface/surfaceCoarsen/Make/files

@@ -1,5 +1,5 @@
 bunnylod/progmesh.C
-bunnylod/vector.C
+bunnylod/vectorb.C
 surfaceCoarsen.C
 
 EXE = $(FOAM_APPBIN)/surfaceCoarsen

applications/utilities/surface/surfaceCoarsen/bunnylod/bunnygut.C

@@ -19,7 +19,7 @@
 #include <GL/gl.h>
 #pragma warning(disable : 4244)
 
-#include "vector.h"
+#include "vectorb.h"
 #include "font.h"
 #include "progmesh.h"
 #include "rabdata.h"

applications/utilities/surface/surfaceCoarsen/bunnylod/progmesh.C

@@ -13,8 +13,8 @@
 #include <assert.h>
 //#include <windows.h>
 
-#include "vector.h"
+#include "vectorb.h"
-#include "list.h"
+#include "listb.h"
 #include "progmesh.h"
 
 #define min(x,y) (((x) <= (y)) ? (x) : (y))

applications/utilities/surface/surfaceCoarsen/bunnylod/progmesh.h

@@ -18,8 +18,8 @@
 #ifndef PROGRESSIVE_MESH_H
 #define PROGRESSIVE_MESH_H
 
-#include "vector.h"
+#include "vectorb.h"
-#include "list.h"
+#include "listb.h"
 
 class tridata {
   public:

applications/utilities/surface/surfaceCoarsen/bunnylod/vectorb.C

@@ -3,7 +3,7 @@
 #include <math.h>
 #include <assert.h>
 
-#include "vector.h"
+#include "vectorb.h"
 
 float  sqr(float a) {return a*a;}
 

bin/tools/RunFunctions

@@ -73,6 +73,11 @@ runParallel ()
         LOG_NAME=log.$APP_NAME
     fi
 
+    if [ "$WM_OSTYPE" = "MSWindows" ]
+    then
+        APP_RUN="${APP_RUN}.exe"
+    fi
+
     if [ -f $LOG_NAME ] ; then
         echo "$APP_NAME already run on $PWD: remove log file to run"
     else
@@ -139,6 +144,17 @@ cloneCase ()
 
 makeFsiCaseLinks ()
 {
+    if [ "$WM_OSTYPE" = "MSWindows" ]
+    then
+        cd $1
+        cd system
+        cp -r ../../$2/system $2
+        cd ../constant
+        cp -r ../../$2/constant $2
+        cd ../0
+        cp -r ../../$2/0 $2
+        cd ../..
+    else
         cd $1
         cd system
         ln -s ../../$2/system $2
@@ -147,10 +163,23 @@ makeFsiCaseLinks ()
         cd ../0
         ln -s ../../$2/0 $2
         cd ../..
+    fi
 }
 
 makeFsiResultsLinks ()
 {
+    if [ "$WM_OSTYPE" = "MSWindows" ]
+    then
+        cd $1
+        TIME_DIRS=`foamInfoExec -times | sed '1,/constant/d'`
+        echo "makeFsiResultsLinks for" $TIME_DIRS
+        cd ../$2
+        for T in $TIME_DIRS
+        do
+            cp -r ../$1/${T}/solid ${T}
+        done
+        cd ..
+    else
         cd $1
         TIME_DIRS=`foamInfoExec -times | sed '1,/constant/d'`
         echo "makeFsiResultsLinks for" $TIME_DIRS
@@ -160,6 +189,7 @@ makeFsiResultsLinks ()
             ln -s ../$1/${T}/solid ${T}
         done
         cd ..
+    fi
 }
 
 #------------------------------------------------------------------------------

etc/bashrc

@@ -160,10 +160,11 @@ done
 : ${WM_OSTYPE:=POSIX}; export WM_OSTYPE
 
 
-# Compiler: set to Gcc or Icc (for Intel's icc)
+# Compiler: set to Gcc, Icc (for Intel's icc) or Clang
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 : ${WM_COMPILER:=Gcc}; export WM_COMPILER
 #: ${WM_COMPILER:=Icc}; export WM_COMPILER
+#: ${WM_COMPILER:=Clang}; export WM_COMPILER
 
 export WM_COMPILER_ARCH=
 export WM_COMPILER_LIB_ARCH=
@@ -316,6 +317,19 @@ Linux)
         export WM_LDFLAGS='-mfloat-abi=hard'
         ;;
 
+    aarch64)
+        # https://gcc.gnu.org/onlinedocs/gcc/AArch64-Options.html
+        # This command will provide more details for your specific processor
+        # gcc -c -Q -march=native -mtune=native --help=target
+        WM_ARCH=linuxARM8
+        export WM_COMPILER_LIB_ARCH=64
+        export WM_CC='gcc'
+        export WM_CXX='g++'
+        export WM_CFLAGS='-fPIC -march=native -mtune=native'
+        export WM_CXXFLAGS='-fPIC -march=native -mtune=native'
+        export WM_LDFLAGS=''
+        ;;
+
     *)
         echo Unknown processor type `uname -m` for Linux
         ;;
@@ -366,7 +380,7 @@ Darwin)
         echo "Using Macports binaries"
     fi
 
-    export WM_USE_MACPORT=1
+    export WM_USE_MACPORT=0
     export WM_BASE_COMPILER=`echo $WM_COMPILER | tr -d "[:digit:]"`
     export WM_MACPORT_MPI_VERSION=`echo $WM_COMPILER | tr "[:upper:]" "[:lower:]"`
     export WM_MACPORT_VERSION=`echo $WM_MACPORT_MPI_VERSION | tr -d "[:alpha:]" | sed -e "s/\(.\)\(.\)/\1\.\2/"`
@@ -434,9 +448,15 @@ Darwin)
             export WM_CXX="g++-mp-$WM_MACPORT_VERSION"
             export WM_FC="gfortran-mp-$WM_MACPORT_VERSION"
         elif [ "$WM_BASE_COMPILER" == "Clang" ]
+        then
+            if [ "$WM_USE_MACPORT" == "1" ]
             then
                 export WM_CC="clang-mp-$WM_MACPORT_VERSION"
                 export WM_CXX="clang++-mp-$WM_MACPORT_VERSION"
+            else
+                export WM_CC="clang"
+                export WM_CXX="clang++"
+            fi
         # Seems like there is no Fortran-frontend for LLVM at thistime
         elif [ "$WM_BASE_COMPILER" == "Dragonegg" ]
         then

etc/prefs.sh-EXAMPLE

@@ -160,8 +160,10 @@ export FOAM_VERBOSE=1
 # System installed bison
 #export BISON_SYSTEM=1
 
-# System installed flex
+# System installed flex. FLEX_DIR should point to the directory where
+# $FLEX_DIR/bin/flex is located
 #export FLEX_SYSTEM=1
+#export FLEX_DIR=/usr
 
 # System installed m4
 #export M4_SYSTEM=1

src/ODE/Make/files

@@ -15,6 +15,23 @@ $(translationODE)/translationODE.C
 sixDOF = sixDOF
 $(sixDOF)/finiteRotation/finiteRotation.C
 $(sixDOF)/sixDOFqODE/sixDOFqODE.C
-$(sixDOF)/sixDOFbodies/sixDOFbodies.C
+
+$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/rotationalConstraint/rotationalConstraint.C
+$(sixDOF)/sixDOFODE/constraints/rotationalConstraints/constantAngularAcceleration/constantAngularAcceleration.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C
+$(sixDOF)/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.C
+
+$(sixDOF)/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.C
+$(sixDOF)/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.C
+$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.C
+$(sixDOF)/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.C
+
+$(sixDOF)/sixDOFODE/sixDOFODE.C
+$(sixDOF)/sixDOFODE/newSixDOFODE.C
+$(sixDOF)/quaternionSixDOF/quaternionSixDOF.C
+$(sixDOF)/geometricSixDOF/geometricSixDOF.C
+
+$(sixDOF)/sixDOFBodies/sixDOFBodies.C
 
 LIB = $(FOAM_LIBBIN)/libODE

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -53,7 +53,7 @@ class ODESolver
 
 protected:
 
-    // Private data
+    // Protected data
 
         //- Reference to ODE
         ODE& ode_;
@@ -62,6 +62,8 @@ protected:
         mutable scalarField dydx_;
 
 
+private:
+
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
@@ -126,7 +128,7 @@ public:
         ) const = 0;
 
 
-        virtual void solve
+        void solve
         (
             const scalar xStart,
             const scalar xEnd,

src/ODE/sixDOF/finiteRotation/finiteRotation.C

@@ -26,22 +26,21 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Update by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 \*---------------------------------------------------------------------------*/
 
 #include "objectRegistry.H"
 #include "finiteRotation.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Static Functions  * * * * * * * * * * * * * //
 
 Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
 {
     vector ur = - *( inv(I + rotT) & (I - rotT) );
 
-    // Scaling to a unit vector.  HJ, problems with round-off
+    // Scaling to a unit vector. Problems with round-off. HJ, 4/Aug/2008
-    // HJ, 4/Aug/2008
-
     if (mag(ur) > SMALL)
     {
         return ur/(mag(ur) + SMALL);
@@ -49,8 +48,7 @@ Foam::vector Foam::finiteRotation::rotVector(const tensor& rotT)
     else
     {
         // Rotation vector is undertermined at zero rotation
-        // Returning arbitrary unit vector
+        // Returning arbitrary unit vector. HJ, 4/Mar/2015
-        // HJ, 4/Mar/2015
         return vector(0, 0, 1);
     }
 }
@@ -79,22 +77,15 @@ Foam::vector Foam::finiteRotation::eulerAngles(const tensor& rotT)
     scalar& pitchAngle = eulerAngles.y();
     scalar& yawAngle = eulerAngles.z();
 
-    // Calculate roll angle
-    rollAngle = atan2(rotT.yz(), rotT.zz());
-
-    // Use mag to avoid negative value due to round-off
-    // HJ, 24/Feb/2016
-    // Bugfix: sqr. SS, 18/Apr/2016
-    const scalar c2 = sqrt(sqr(rotT.xx()) + sqr(rotT.xy()));
-
     // Calculate pitch angle
-    pitchAngle = atan2(-rotT.xz(), c2);
+    pitchAngle = asin(rotT.xz());
 
-    const scalar s1 = sin(rollAngle);
+    // Calculate roll angle
-    const scalar c1 = cos(rollAngle);
+    const scalar cosPitch = cos(pitchAngle);
+    rollAngle = asin(-rotT.yz()/cosPitch);
 
     // Calculate yaw angle
-    yawAngle = atan2(s1*rotT.zx() - c1*rotT.yx(), c1*rotT.yy() - s1*rotT.zy());
+    yawAngle = asin(-rotT.xy()/cosPitch);
 
     return eulerAngles;
 }
@@ -122,6 +113,14 @@ Foam::finiteRotation::finiteRotation
 {}
 
 
+Foam::finiteRotation::finiteRotation(const tensor& R)
+:
+    eInitial_(R),
+    rotTensor_(R),
+    rotIncrementTensor_(tensor::zero)
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::finiteRotation::~finiteRotation()
@@ -130,12 +129,23 @@ Foam::finiteRotation::~finiteRotation()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::finiteRotation::updateRotation(const tensor& R)
+{
+    rotIncrementTensor_ = (R & rotTensor_.T());
+    rotTensor_ = R;
+}
+
+
 void Foam::finiteRotation::updateRotation(const HamiltonRodriguezRot& rot)
 {
-    tensor rotR = rot.R();
+    updateRotation(rot.R());
+}
 
-    rotIncrementTensor_ = (rotR & (rotTensor_.T()));
+
-    rotTensor_ = rotR;
+void Foam::finiteRotation::updateRotationWithIncrement(const tensor& incR)
+{
+    rotIncrementTensor_ = incR;
+    rotTensor_ = (incR & rotTensor_);
 }
 
 

src/ODE/sixDOF/finiteRotation/finiteRotation.H

@@ -26,7 +26,8 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Updated by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 SourceFiles
     finiteRotation.C
@@ -61,20 +62,6 @@ class finiteRotation
         tensor rotIncrementTensor_;
 
 
-    // Private Member Functions
-
-        //- Calculate unit rotation vector from given rotation tensor
-        static vector rotVector(const tensor& rotT);
-
-        //- Calculate rotation angle from given rotation tensor
-        static scalar rotAngle(const tensor& rotT);
-
-        //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention) from
-        //  given rotation tensor. Reference: Mike Day, Insomniac Games,
-        //  Extracting Euler Angles from a Rotation Matrix.
-        static vector eulerAngles(const tensor& rotT);
-
-
 public:
 
     // Constructors
@@ -89,17 +76,40 @@ public:
             const scalar& angle
         );
 
+        //- Construct from rotation tensor
+        explicit finiteRotation(const tensor& R);
+
 
     // Destructor
 
         ~finiteRotation();
 
 
+    // Static Functions
+
+        //- Calculate unit rotation vector from given rotation tensor
+        static vector rotVector(const tensor& rotT);
+
+        //- Calculate rotation angle from given rotation tensor
+        static scalar rotAngle(const tensor& rotT);
+
+        //- Calculate Euler angles (x-y-z (roll-pitch-yaw) convention by Bryan)
+        //  given the rotation tensor (global-to-local transformation).
+        //  Reference: Nikravesh: Computer-Aided Analysis of Mechanical Systems
+        static vector eulerAngles(const tensor& rotT);
+
+
     // Member Functions
 
-        //- Update rotation
+        //- Update rotation given rotation tensor
+        void updateRotation(const tensor& R);
+
+        //- Update rotation given HamiltonRodriguezRot (quaternions)
         void updateRotation(const HamiltonRodriguezRot& rot);
 
+        //- Update rotation given increment rotation tensor
+        void updateRotationWithIncrement(const tensor& incR);
+
         //- Return initial quaternions
         const HamiltonRodriguezRot& eInitial() const;
 

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.C

@@ -0,0 +1,646 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    geometricSixDOF
+
+\*---------------------------------------------------------------------------*/
+
+#include "geometricSixDOF.H"
+#include "scalarSquareMatrix.H"
+#include "OutputControlDictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(geometricSixDOF, 0);
+addToRunTimeSelectionTable(sixDOFODE, geometricSixDOF, dictionary);
+
+}
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncTensorTol_
+(
+    "geometricSixDOFRotIncTensorTol",
+    1e-10
+);
+
+
+const Foam::debug::tolerancesSwitch Foam::geometricSixDOF::rotIncRateTol_
+(
+    "geometricSixDOFRotIncRateTol",
+    1e-6
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::dimensionedVector Foam::geometricSixDOF::A
+(
+    const dimensionedVector& xR,
+    const dimensionedVector& uR,
+    const tensor& R,
+    const scalar t
+) const
+{
+    // 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);
+
+    // Insert mass and explicit forcing into the system. Note: translations are
+    // solved in the global coordinate system and the explicit forcing contains
+    // restraining forces
+    const dimensionedVector explicitForcing =
+        force
+        (
+            t,
+            R.T(),
+            xR.value(),
+            uR.value()
+        );
+
+    const vector& efVal = explicitForcing.value();
+    const scalar& m = mass().value();
+
+    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 = translationalConstraints()[tcI];
+
+        // Get matrix contribution from constraint
+        const vector mc =
+            curTc.matrixContribution
+            (
+                t,
+                R.T(),
+                xR.value(),
+                uR.value()
+            );
+
+        // Get matrix index
+        const label index = tcI + 3;
+
+        // Insert contributions into the matrix
+        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(t, R.T(), xR.value(), uR.value());
+    }
+
+    // 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[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::OmegaDot
+(
+    const tensor& R,
+    const dimensionedVector& omega,
+    const scalar t
+) const
+{
+    // 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);
+
+    // 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
+      + (
+          RT
+        & moment
+          (
+              t,
+              RT.value(),
+              omega.value()
+          )
+        ) // External torque with restraints
+    );
+    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
+            (
+                t,
+                RT.value(),
+                omega.value()
+            );
+
+        // Get matrix index
+        const label index = rcI + 3;
+
+        // Insert contributions into the matrix
+        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(t, RT.value(), omega.value());
+    }
+
+    // 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
+        dimensionedVector
+        (
+            "OmegaDot",
+            moment().dimensions()/momentOfInertia().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::E
+(
+    const dimensionedVector& omega
+) const
+{
+    return (*(momentOfInertia() & omega) & omega);
+}
+
+
+Foam::tensor Foam::geometricSixDOF::expMap(const vector& rotInc) const
+{
+    tensor R;
+
+    // Calculate the magnitude of the rotational increment vector
+    const scalar magRotInc = mag(rotInc);
+
+    if (magRotInc < rotIncTensorTol_())
+    {
+        // No rotational increment
+        R = I;
+    }
+    else
+    {
+        // Calculate rotational increment tensor using the exponential map
+
+        // Skew-symmetric tensor corresponding to the rotation increment
+        const tensor skewRotInc(*rotInc);
+
+        R = I
+          + skewRotInc*sin(magRotInc)/magRotInc
+          + (skewRotInc & skewRotInc)*(1.0 - cos(magRotInc))/sqr(magRotInc);
+    }
+
+    return R;
+}
+
+
+Foam::vector Foam::geometricSixDOF::dexpMap
+(
+    const vector& rotInc,
+    const vector& omega
+) const
+{
+    vector rotIncDot;
+
+    // Calculate the magnitude of the rotational increment vector
+    const scalar magRotInc = mag(rotInc);
+
+    if (magRotInc < rotIncRateTol_())
+    {
+        // Stabilised calculation of rotation increment to avoid small
+        // denominators
+        const tensor lbRotIncOmega(lieBracket(rotInc, omega));
+
+        rotIncDot =
+        (
+            I
+          + lbRotIncOmega/2.0
+          + lieBracket(rotInc, *lbRotIncOmega)/12.0
+        ) & omega;
+    }
+    else
+    {
+        // Calculate rate of the rotational increment vector using the
+        // differential of the exponential map
+
+        // Skew-symmetric tensor corresponding to the rotation increment
+        const tensor skewRotInc(*rotInc);
+
+        rotIncDot =
+        (
+            I
+          + 0.5*skewRotInc
+          - (skewRotInc & skewRotInc)*
+            (magRotInc/tan(magRotInc/2.0) - 2.0)/(2.0*sqr(magRotInc))
+        ) & omega;
+    }
+
+    return rotIncDot;
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::geometricSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to geometricSixDOF&
+    const geometricSixDOF& gsd = refCast<const geometricSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = gsd.Xrel_;
+    U_ = gsd.U_;
+    Uaverage_ = gsd.Uaverage_;
+    rotation_ = gsd.rotation_;
+    omega_ = gsd.omega_;
+    omegaAverage_ = gsd.omegaAverage_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = gsd.coeffs_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::geometricSixDOF(const IOobject& io)
+:
+    sixDOFODE(io),
+
+    Xrel_(dict().lookup("Xrel")),
+    U_(dict().lookup("U")),
+    Uaverage_("Uaverage", U_),
+    rotation_(tensor(dict().lookup("rotationTensor"))),
+    rotIncrement_
+    (
+        dict().lookupOrDefault<tensor>("rotationIncrementTensor", tensor::zero)
+    ),
+    omega_(dict().lookup("omega")),
+    omegaAverage_("omegaAverage", omega_),
+
+    coeffs_(12, 0.0)
+{
+    // Set ODE coefficients from position and rotation
+
+    // Linear displacement relative to spring equilibrium
+    const vector& Xval = Xrel_.value();
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Linear velocity
+    const vector& Uval = U_.value();
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Rotational velocity in non - inertial coordinate system
+    const vector& omegaVal = omega_.value();
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Increment of the rotation vector (zero for initial condition)
+    coeffs_[9] = 0;
+    coeffs_[10] = 0;
+    coeffs_[11] = 0;
+}
+
+
+Foam::geometricSixDOF::geometricSixDOF
+(
+    const word& name,
+    const geometricSixDOF& gsd
+)
+:
+    sixDOFODE(name, gsd),
+
+    Xrel_(gsd.Xrel_.name(), gsd.Xrel_),
+    U_(gsd.U_.name(), gsd.U_),
+    Uaverage_(gsd.Uaverage_.name(), gsd.Uaverage_),
+    rotation_(gsd.rotation_),
+    omega_(gsd.omega_.name(), gsd.omega_),
+    omegaAverage_(gsd.omegaAverage_.name(), gsd.omegaAverage_),
+
+    coeffs_(gsd.coeffs_)
+{}
+
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::geometricSixDOF::clone
+(
+    const word& name
+) const
+{
+    // Create and return an autoPtr to the new geometricSixDOF object with a
+    // new name
+    return autoPtr<sixDOFODE>
+    (
+        new geometricSixDOF
+        (
+            name,
+            *this
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::geometricSixDOF::~geometricSixDOF()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Xrel() const
+{
+    return Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omega() const
+{
+    return omega_;
+}
+
+
+Foam::dimensionedVector Foam::geometricSixDOF::X() const
+{
+    return Xequilibrium() + Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::U() const
+{
+    return U_;
+}
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::Uaverage() const
+{
+    return Uaverage_;
+}
+
+
+const Foam::dimensionedVector& Foam::geometricSixDOF::omegaAverage() const
+{
+    return omegaAverage_;
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toRelative() const
+{
+    return rotation_.T();
+}
+
+
+Foam::tensor Foam::geometricSixDOF::toAbsolute() const
+{
+    return rotation_;
+}
+
+
+const Foam::tensor& Foam::geometricSixDOF::rotIncrementTensor() const
+{
+    return rotIncrement_;
+}
+
+
+void Foam::geometricSixDOF::derivatives
+(
+    const scalar x,
+    const scalarField& y,
+    scalarField& dydx
+) const
+{
+    // Translation
+
+    // Set the derivatives for displacement
+    dydx[0] = y[3];
+    dydx[1] = y[4];
+    dydx[2] = y[5];
+
+    dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
+    dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
+
+    // Get rotational increment vector (u)
+    const vector rotIncrementVector(y[9], y[10], y[11]);
+
+    // Calculate current rotation tensor obtained with exponential map
+    const tensor curRot = (rotation_ & expMap(rotIncrementVector));
+
+    // Calculate translational acceleration using current rotation
+    const vector accel = A(curX, curU, curRot, x).value();
+
+    // Set the derivatives for velocity
+    dydx[3] = accel.x();
+    dydx[4] = accel.y();
+    dydx[5] = accel.z();
+
+    // Rotation
+
+    dimensionedVector curOmega
+    (
+        "curOmega",
+        dimless/dimTime,
+        vector(y[6], y[7], y[8])
+    );
+
+    // Calculate rotational acceleration using current rotation
+    const vector omegaDot = OmegaDot(curRot, curOmega, x).value();
+
+    dydx[6] = omegaDot.x();
+    dydx[7] = omegaDot.y();
+    dydx[8] = omegaDot.z();
+
+    // Calculate derivative of rotIncrementVector using current rotation
+    // information
+    const vector rotIncrementVectorDot =
+        dexpMap(rotIncrementVector, curOmega.value());
+
+    // Set the derivatives for rotation
+    dydx[9] = rotIncrementVectorDot.x();
+    dydx[10] = rotIncrementVectorDot.y();
+    dydx[11] = rotIncrementVectorDot.z();
+}
+
+
+void Foam::geometricSixDOF::update(const scalar delta)
+{
+    // Translation
+
+    // Get displacement
+    const vector Xold = Xrel_.value();
+
+    vector& Xval = Xrel_.value();
+    Xval.x() = coeffs_[0];
+    Xval.y() = coeffs_[1];
+    Xval.z() = coeffs_[2];
+
+    // Get velocity
+    vector& Uval = U_.value();
+    Uval.x() = coeffs_[3];
+    Uval.y() = coeffs_[4];
+    Uval.z() = coeffs_[5];
+
+    // Stabilise translational constraints if necessary
+    forAll(translationalConstraints(), tcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        translationalConstraints()[tcI].stabilise(t, Xval, Uval);
+    }
+
+    // Update (possibly constrained) displacement
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Update (possibly constrained) velocity
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Update average velocity
+    Uaverage_.value() = (Xval - Xold)/delta;
+
+
+    // Rotation
+
+    // Get angular velocity
+    const vector omegaOld = omega_.value();
+
+    vector& omegaVal = omega_.value();
+    omegaVal.x() = coeffs_[6];
+    omegaVal.y() = coeffs_[7];
+    omegaVal.z() = coeffs_[8];
+
+    // Update rotational increment tensor
+    rotIncrement_ = expMap(vector(coeffs_[9], coeffs_[10], coeffs_[11]));
+
+    // Update rotational tensor
+    rotation_ = (rotation_ & rotIncrement_);
+
+    // Stabilise rotational constraints if necessary
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        rotationalConstraints()[rcI].stabilise(t, omegaVal);
+    }
+
+    // Update (possibly constrained) omega
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Reset increment vector in coefficients for the next step
+    coeffs_[9] = 0;
+    coeffs_[10] = 0;
+    coeffs_[11] = 0;
+
+    // Update average omega
+    omegaAverage_.value() = 0.5*(omegaVal + omegaOld);
+}
+
+
+bool Foam::geometricSixDOF::writeData(Ostream& os) const
+{
+    // First write the part related to base class subobject
+    sixDOFODE::writeData(os);
+
+    // Write type name
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    // Write data
+    os.writeKeyword("Xrel") << tab << Xrel_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationTensor") << tab << rotation_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationIncrementTensor") << tab << rotIncrement_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << omega_
+        << token::END_STATEMENT << nl << endl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/geometricSixDOF/geometricSixDOF.H

@@ -0,0 +1,319 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    geometricSixDOF
+
+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 translational/rotationalConstraint classes.
+
+    Reference (bibtex):
+
+    @article {terzeEtAl2016,
+        Author = {Terze, Z. and M\"{u}ller, A. and Zlatar, D.},
+        title = {Singularity-free time integration of rotational quaternions
+        using non-redundant ordinary differential equations},
+        Journal = {Multibody System Dynamics},
+        Year = {2016},
+        Volume = {38},
+        Number = {3},
+        Pages = {201--225}
+    }
+
+    Note on convention: rotation tensor (R, or rotation_) defines
+    body-to-inertial coordinate system transformation
+    (local-to-global). Opposite as in quaternionSixDOF.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    geometricSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef geometricSixDOF_H
+#define geometricSixDOF_H
+
+#include "sixDOFODE.H"
+#include "tolerancesSwitch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class geometricSixDOF Declaration
+\*---------------------------------------------------------------------------*/
+
+class geometricSixDOF
+:
+    public sixDOFODE
+{
+    // Private data
+
+        // Initial body state variables
+
+            //- Displacement relative to spring equilibrium
+            dimensionedVector Xrel_;
+
+            //- Velocity of mass centroid
+            dimensionedVector U_;
+
+            //- Average velocity of mass centroid (evaluated at midstep)
+            dimensionedVector Uaverage_;
+
+            //- Rotation tensor
+            tensor rotation_;
+
+            //- Rotational increment tensor
+            tensor rotIncrement_;
+
+            //- Rotational velocity about mass centroid
+            dimensionedVector omega_;
+
+            //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
+            dimensionedVector omegaAverage_;
+
+
+        //- ODE coefficients
+        scalarField coeffs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        geometricSixDOF(const geometricSixDOF&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const geometricSixDOF&);
+
+
+        // Variables in relative coordinate system (solved for)
+
+            //- Return acceleration in relative coordinate system
+            //  given current values of relative displacement, velocity,
+            //  orientation (rotation tensor) and time
+            dimensionedVector A
+            (
+                const dimensionedVector& xR,
+                const dimensionedVector& uR,
+                const tensor& R,
+                const scalar t
+            ) const;
+
+
+            //- Return rotational acceleration in relative coordinate system
+            //  given current values for relative rotational velocity,
+            //  orientation (rotation tensor) and time
+            dimensionedVector OmegaDot
+            (
+                const tensor& R,
+                const dimensionedVector& omega,
+                const scalar t
+            ) const;
+
+            //- Return the Euler part of moment equation
+            dimensionedVector E
+            (
+                const dimensionedVector& omega
+            ) const;
+
+            //- Exponential map used to calculate increment of the rotation
+            //  tensor
+            tensor expMap(const vector& rotInc) const;
+
+            //- Differential of the expontential map used to calculate the time
+            //  derivative of rotation increment vector
+            vector dexpMap(const vector& rotInc, const vector& omega) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+public:
+
+    // Run-time type information
+    TypeName("geometricSixDOF");
+
+
+    // Static data members
+
+        //- Rotational increment tensor tolerance. Used in expMap member
+        //  function in case the rotation is negligibly small
+        static const debug::tolerancesSwitch rotIncTensorTol_;
+
+        //- Rotational increment rate of change tolerance. Used in dexpMap
+        //  member function in case the rotation rate is negligibly small
+        static const debug::tolerancesSwitch rotIncRateTol_;
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        geometricSixDOF(const IOobject& io);
+
+        //- Construct geometricSixDOF object, changing name
+        geometricSixDOF
+        (
+            const word& name,
+            const geometricSixDOF& gsd
+        );
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const;
+
+
+    // Destructor
+
+        virtual ~geometricSixDOF();
+
+
+    // Member Functions
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const
+                {
+                    return 12;
+                }
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs()
+                {
+                    return coeffs_;
+                }
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const
+                {
+                    return coeffs_;
+                }
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const
+                {
+                    notImplemented
+                    (
+                        "geometricSixDOF::jacobian\n"
+                        "(\n"
+                        "    const scalar x,\n"
+                        "    const scalarField& y,\n"
+                        "    scalarField& dfdx,\n"
+                        "    scalarSquareMatrix& dfdy,\n"
+                        ") const"
+                    );
+                }
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta);
+
+
+        // Write controls
+
+            //- WriteData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.C

@@ -0,0 +1,557 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    quaternionSixDOF
+
+\*---------------------------------------------------------------------------*/
+
+#include "quaternionSixDOF.H"
+#include "OutputControlDictionary.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(quaternionSixDOF, 0);
+addToRunTimeSelectionTable(sixDOFODE, quaternionSixDOF, dictionary);
+
+}
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::dimensionedVector Foam::quaternionSixDOF::A
+(
+    const dimensionedVector& xR,
+    const dimensionedVector& uR,
+    const HamiltonRodriguezRot& rotation,
+    const scalar t
+) const
+{
+    // 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);
+
+    // Insert mass and explicit forcing into the system. Note: translations are
+    // solved in the global coordinate system and the explicit forcing contains
+    // restraining forces
+    const dimensionedVector explicitForcing =
+        force
+        (
+            t,
+            rotation.R(),
+            xR.value(),
+            uR.value()
+        );
+
+    const vector& efVal = explicitForcing.value();
+    const scalar& m = mass().value();
+
+    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 = translationalConstraints()[tcI];
+
+        // Get matrix contribution from constraint
+        const vector mc =
+            curTc.matrixContribution
+            (
+                t,
+                rotation.R(),
+                xR.value(),
+                uR.value()
+            );
+
+        // Get matrix index
+        const label index = tcI + 3;
+
+        // Insert contributions into the matrix
+        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
+           (
+               t,
+               rotation.R(),
+               xR.value(),
+               uR.value()
+           );
+    }
+
+    // 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[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::OmegaDot
+(
+    const HamiltonRodriguezRot& rotation,
+    const dimensionedVector& omega,
+    const scalar t
+) const
+{
+    // 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);
+
+    // Get current inertial-to-local transformation. Note: different convention
+    // (R represents coordinate transformation from global to local)
+    const dimensionedTensor R("R", dimless, rotation.R());
+
+    // Insert moment of inertia and explicit forcing into the system
+    const dimensionedVector explicitForcing
+    (
+        E(omega) // Euler part
+      + (
+          R.value()
+        & moment
+          (
+              t,
+              R.value(),
+              omega.value()
+          )
+        ) // External torque with restraints
+    );
+    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
+            (
+                t,
+                R.value(),
+                omega.value()
+            );
+
+        // Get matrix index
+        const label index = rcI + 3;
+
+        // Insert contributions into the matrix
+        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(t, R.value(), omega.value());
+    }
+
+    // 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
+        dimensionedVector
+        (
+            "OmegaDot",
+            moment().dimensions()/momentOfInertia().dimensions(),
+            vector(rhs[0], rhs[1], rhs[2])
+        );
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::E
+(
+    const dimensionedVector& omega
+) const
+{
+    return (*(momentOfInertia() & omega) & omega);
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::quaternionSixDOF::setState(const sixDOFODE& sd)
+{
+    // First set the state in base class subobject
+    sixDOFODE::setState(sd);
+
+    // Cast sixDOFODE& to quaternionSixDOF&
+    const quaternionSixDOF& qsd = refCast<const quaternionSixDOF>(sd);
+
+    // Set state variables for this class
+    Xrel_ = qsd.Xrel_;
+    U_ = qsd.U_;
+    Uaverage_ = qsd.Uaverage_;
+    rotation_ = qsd.rotation_;
+    omega_ = qsd.omega_;
+    omegaAverage_ = qsd.omegaAverage_;
+
+    // Copy ODE coefficients: this carries actual ODE state
+    // HJ, 23/Mar/2015
+    coeffs_ = qsd.coeffs_;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::quaternionSixDOF::quaternionSixDOF(const IOobject& io)
+:
+    sixDOFODE(io),
+
+    Xrel_(dict().lookup("Xrel")),
+    U_(dict().lookup("U")),
+    Uaverage_("Uaverage", U_),
+    rotation_
+    (
+        vector(dict().lookup("rotationVector")),
+        dimensionedScalar(dict().lookup("rotationAngle")).value()
+    ),
+    omega_(dict().lookup("omega")),
+    omegaAverage_("omegaAverage", omega_),
+
+    coeffs_(13, 0.0)
+{
+    // Set ODE coefficients from position and rotation
+
+    // Linear displacement relative to spring equilibrium
+    const vector& Xval = Xrel_.value();
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Linear velocity
+    const vector& Uval = U_.value();
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Rotational velocity in non - inertial coordinate system
+    const vector& omegaVal = omega_.value();
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Quaternions
+    coeffs_[9] = rotation_.eInitial().e0();
+    coeffs_[10] = rotation_.eInitial().e1();
+    coeffs_[11] = rotation_.eInitial().e2();
+    coeffs_[12] = rotation_.eInitial().e3();
+}
+
+
+Foam::quaternionSixDOF::quaternionSixDOF
+(
+    const word& name,
+    const quaternionSixDOF& qsd
+)
+:
+    sixDOFODE(name, qsd),
+
+    Xrel_(qsd.Xrel_.name(), qsd.Xrel_),
+    U_(qsd.U_.name(), qsd.U_),
+    Uaverage_(qsd.Uaverage_.name(), qsd.Uaverage_),
+    rotation_(qsd.rotation_),
+    omega_(qsd.omega_.name(), qsd.omega_),
+    omegaAverage_(qsd.omegaAverage_.name(), qsd.omegaAverage_),
+
+    coeffs_(qsd.coeffs_)
+{}
+
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::quaternionSixDOF::clone
+(
+    const word& name
+) const
+{
+    // Create and return an autoPtr to the new quaternionSixDOF object with a
+    // new name
+    return autoPtr<sixDOFODE>
+    (
+        new quaternionSixDOF
+        (
+            name,
+            *this
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::quaternionSixDOF::~quaternionSixDOF()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::Xrel() const
+{
+    return Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::omega() const
+{
+    return omega_;
+}
+
+
+Foam::dimensionedVector Foam::quaternionSixDOF::X() const
+{
+    return Xequilibrium() + Xrel_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::U() const
+{
+    return U_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::Uaverage() const
+{
+    return Uaverage_;
+}
+
+
+const Foam::dimensionedVector& Foam::quaternionSixDOF::omegaAverage() const
+{
+    return omegaAverage_;
+}
+
+
+Foam::tensor Foam::quaternionSixDOF::toRelative() const
+{
+    return rotation_.eCurrent().R();
+}
+
+
+Foam::tensor Foam::quaternionSixDOF::toAbsolute() const
+{
+    return rotation_.eCurrent().invR();
+}
+
+
+const Foam::tensor& Foam::quaternionSixDOF::rotIncrementTensor() const
+{
+    return rotation_.rotIncrementTensor();
+}
+
+
+void Foam::quaternionSixDOF::derivatives
+(
+    const scalar x,
+    const scalarField& y,
+    scalarField& dydx
+) const
+{
+    // Set the derivatives for displacement
+    dydx[0] = y[3];
+    dydx[1] = y[4];
+    dydx[2] = y[5];
+
+    dimensionedVector curX("curX", dimLength, vector(y[0], y[1], y[2]));
+    dimensionedVector curU("curU", dimVelocity, vector(y[3], y[4], y[5]));
+    const HamiltonRodriguezRot curRotation
+    (
+        y[9],
+        y[10],
+        y[11],
+        y[12]
+    );
+
+    const vector accel = A(curX, curU, curRotation, x).value();
+
+    dydx[3] = accel.x();
+    dydx[4] = accel.y();
+    dydx[5] = accel.z();
+
+    // Set the derivatives for rotation
+    dimensionedVector curOmega
+    (
+        "curOmega",
+        dimless/dimTime,
+        vector(y[6], y[7], y[8])
+    );
+
+    const vector omegaDot = OmegaDot(curRotation, curOmega, x).value();
+
+    dydx[6] = omegaDot.x();
+    dydx[7] = omegaDot.y();
+    dydx[8] = omegaDot.z();
+
+    dydx[9] = curRotation.eDot(curOmega.value(), 0);
+    dydx[10] = curRotation.eDot(curOmega.value(), 1);
+    dydx[11] = curRotation.eDot(curOmega.value(), 2);
+    dydx[12] = curRotation.eDot(curOmega.value(), 3);
+}
+
+
+void Foam::quaternionSixDOF::update(const scalar delta)
+{
+    // Translation
+
+    // Get displacement
+    const vector Xold = Xrel_.value();
+
+    vector& Xval = Xrel_.value();
+    Xval.x() = coeffs_[0];
+    Xval.y() = coeffs_[1];
+    Xval.z() = coeffs_[2];
+
+    // Get velocity
+    vector& Uval = U_.value();
+    Uval.x() = coeffs_[3];
+    Uval.y() = coeffs_[4];
+    Uval.z() = coeffs_[5];
+
+    // Stabilise translational constraints if necessary
+    forAll(translationalConstraints(), tcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        translationalConstraints()[tcI].stabilise(t, Xval, Uval);
+    }
+
+    // Update (possibly constrained) displacement
+    coeffs_[0] = Xval.x();
+    coeffs_[1] = Xval.y();
+    coeffs_[2] = Xval.z();
+
+    // Update (possibly constrained) velocity
+    coeffs_[3] = Uval.x();
+    coeffs_[4] = Uval.y();
+    coeffs_[5] = Uval.z();
+
+    // Update average velocity
+    Uaverage_.value() = (Xval - Xold)/delta;
+
+
+    // Rotation
+
+    // Get angular velocity
+    vector& omegaVal = omega_.value();
+    omegaVal.x() = coeffs_[6];
+    omegaVal.y() = coeffs_[7];
+    omegaVal.z() = coeffs_[8];
+
+
+    rotation_.updateRotation
+    (
+        HamiltonRodriguezRot
+        (
+            coeffs_[9],
+            coeffs_[10],
+            coeffs_[11],
+            coeffs_[12]
+        )
+    );
+
+    // Stabilise rotational constraints if necessary
+    forAll(rotationalConstraints(), rcI)
+    {
+        // Note: get end time for this ODE step from mesh, assuming that the top
+        // level solves this ode from [t - deltaT, t], yielding solution at
+        // t. Done this way to preserve the interface of ODE class.
+        // VV, 10/Mar/2017.
+        const scalar t = dict().time().value();
+
+        rotationalConstraints()[rcI].stabilise(t, omegaVal);
+    }
+
+    // Update (possibly constrained) omega
+    coeffs_[6] = omegaVal.x();
+    coeffs_[7] = omegaVal.y();
+    coeffs_[8] = omegaVal.z();
+
+    // Update average omega
+    omegaAverage_.value() = rotation_.omegaAverage(delta);
+}
+
+
+bool Foam::quaternionSixDOF::writeData(Ostream& os) const
+{
+    // First write the part related to base class subobject
+    sixDOFODE::writeData(os);
+
+    // Write type name
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    // Write data
+    os.writeKeyword("Xrel") << tab << Xrel_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("U") << tab << U_ << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationVector") << tab << rotation_.rotVector()
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("rotationAngle") << tab
+        << dimensionedScalar("rotAngle", dimless, rotation_.rotAngle())
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("omega") << tab << omega_
+        << token::END_STATEMENT << nl << nl;
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternionSixDOF/quaternionSixDOF.H

@@ -0,0 +1,286 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    quaternionSixDOF
+
+Description
+    6-DOF solver using quaternions
+
+    Run-time selectable constraints are handled via Lagrangian multipliers using
+    the interface from translationa/rotationalConstraint classes.
+
+    Note on convention: rotation tensor R obtained from finiteRotation
+    (rotation_) defines inertial-to-body coordinate system transformation
+    (global-to-local). Opposite as in geometrixSixDOF.
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+
+Reorganization by
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+    Inno Gatin, FSB Zagreb.  All rights reserved.
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    quaternionSixDOF.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef quaternionSixDOF_H
+#define quaternionSixDOF_H
+
+#include "sixDOFODE.H"
+#include "finiteRotation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class quaternionSixDOF Declaration
+\*---------------------------------------------------------------------------*/
+
+class quaternionSixDOF
+:
+    public sixDOFODE
+{
+    // Private data
+
+        // Body state variables
+
+            //- Displacement relative to spring equilibrium
+            dimensionedVector Xrel_;
+
+            //- Velocity of mass centroid
+            dimensionedVector U_;
+
+            //- Average velocity of mass centroid (evaluated at midstep)
+            dimensionedVector Uaverage_;
+
+            //- Finite rotation
+            finiteRotation rotation_;
+
+            //- Rotational velocity about mass centroid
+            dimensionedVector omega_;
+
+            //- Average rotational velocity in relative coordinate system
+            //  (evaluated at midstep)
+            dimensionedVector omegaAverage_;
+
+
+        //- ODE coefficients
+        scalarField coeffs_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        quaternionSixDOF(const quaternionSixDOF&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const quaternionSixDOF&);
+
+
+        // Variables in relative coordinate system (solved for)
+
+            //- Return acceleration in relative coordinate system
+            //  given current values of relative displacement and velocity
+            dimensionedVector A
+            (
+                const dimensionedVector& xR,
+                const dimensionedVector& uR,
+                const HamiltonRodriguezRot& rotation,
+                const scalar t
+            ) const;
+
+
+            //- Return rotational acceleration in relative coordinate system
+            //  given current values for relative rotational velocity
+            dimensionedVector OmegaDot
+            (
+                const HamiltonRodriguezRot& rotation,
+                const dimensionedVector& omega,
+                const scalar t
+            ) const;
+
+            //- Return the Euler part of moment equation
+            dimensionedVector E
+            (
+                const dimensionedVector& omega
+            ) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+public:
+
+    // Run-time type information
+    TypeName("quaternionSixDOF");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        quaternionSixDOF(const IOobject& io);
+
+        //- Construct quaternionSixDOF object, changing name
+        quaternionSixDOF
+        (
+            const word& name,
+            const quaternionSixDOF& qsd
+        );
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const;
+
+
+    // Destructor
+
+        virtual ~quaternionSixDOF();
+
+
+    // Member Functions
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const
+                {
+                    return 13;
+                }
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs()
+                {
+                    return coeffs_;
+                }
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const
+                {
+                    return coeffs_;
+                }
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const
+                {
+                    notImplemented
+                    (
+                        "quaternionSixDOF::jacobian\n"
+                        "(\n"
+                        "    const scalar x,\n"
+                        "    const scalarField& y,\n"
+                        "    scalarField& dfdx,\n"
+                        "    scalarSquareMatrix& dfdy,\n"
+                        ") const"
+                    );
+                }
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta);
+
+
+        // Write controls
+
+            //- WriteData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/quaternions/HamiltonRodriguezRot.H

@@ -26,7 +26,8 @@ Class
 
 Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
-    Updated by Hrvoje Jasak
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
 
 Description
     Rotation defined with 4 parameters: quaternions.
@@ -66,13 +67,13 @@ class HamiltonRodriguezRot
         scalarRectangularMatrix Gt_;
 
         //- Inertial to rotated coordinate system transformation
-        mutable tensor R_;
+        tensor R_;
 
 
     // Private member functions
 
         //- Calculate R_ - inertial to body coord. sys. rotation
-        inline void calculateR() const
+        inline void calculateR()
         {
             R_.xx() = 2*(e1_*e1_ + e0_*e0_ - 0.5);
             R_.xy() = 2*(e1_*e2_ + e0_*e3_);
@@ -133,11 +134,7 @@ public:
         }
 
         //- Construct from rotation vector and angle
-        explicit HamiltonRodriguezRot
+        HamiltonRodriguezRot(const vector& rVect, const scalar& rAngle)
-        (
-            const vector& rVect,
-            const scalar& rAngle
-        )
         :
             Gt_(4, 3),
             R_(tensor::zero)
@@ -164,6 +161,29 @@ public:
             calculateGt();
         }
 
+        //- Construct from rotation tensor
+        explicit HamiltonRodriguezRot(const tensor& R)
+        :
+            Gt_(4, 3),
+            R_(R)
+        {
+            // Calculate Hamilton - Rodriguez (Euler) parameters from rotation
+            // matrix
+
+            // Note: sign of e0_ assumed positive
+            e0_ = Foam::sqrt((tr(R) + 1.0)/4.0);
+
+            // Helper variable
+            const scalar oneByFourEo = 1.0/(4.0*e0_);
+
+            e1_ = oneByFourEo*(R.zy() - R.yz());
+            e2_ = oneByFourEo*(R.xz() - R.zx());
+            e3_ = oneByFourEo*(R.yx() - R.xy());
+
+            // Calculate Gt
+            calculateGt();
+        }
+
 
     // Destructor
 

src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.C

@@ -22,7 +22,7 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    sixDOFbodies
+    sixDOFBodies
 
 Description
     6-DOF solver for multiple bodies
@@ -33,11 +33,11 @@ Author
 \*---------------------------------------------------------------------------*/
 
 #include "objectRegistry.H"
-#include "sixDOFbodies.H"
+#include "sixDOFBodies.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::sixDOFbodies::setBodies()
+void Foam::sixDOFBodies::setBodies()
 {
     // Find if duplicate name existes
     forAll (names_, bodyI)
@@ -51,8 +51,9 @@ void Foam::sixDOFbodies::setBodies()
         {
             if (names_[bodyI] == names_[otherBody])
             {
-                FatalErrorIn("sixDOFbodies::setBodies()")
+                FatalErrorIn("sixDOFBodies::setBodies()")
-                    << "Duplicate names of bodies: this is not allowed"
+                    << "Found duplicate name: " << names_[bodyI]
+                    << " for bodies. This is not allowed."
                     << exit(FatalError);
             }
         }
@@ -66,7 +67,7 @@ void Foam::sixDOFbodies::setBodies()
         odes_.set
         (
             bodyI,
-            new sixDOFqODE
+            sixDOFODE::New
             (
                 IOobject
                 (
@@ -82,19 +83,18 @@ void Foam::sixDOFbodies::setBodies()
         solvers_.set
         (
             bodyI,
-            ODESolver::New
+            ODESolver::New(lookup("solver"), odes_[bodyI])
-            (
-                lookup("solver"),
-                odes_[bodyI]
-            )
         );
+
+        Info<< "Finished creating " << odes_[bodyI].type()
+            << " object for body " << names_[bodyI] << endl;
     }
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::sixDOFbodies::sixDOFbodies
+Foam::sixDOFBodies::sixDOFBodies
 (
     const Time& runTime
 )
@@ -121,7 +121,7 @@ Foam::sixDOFbodies::sixDOFbodies
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::sixDOFbodies::solve()
+void Foam::sixDOFBodies::solve()
 {
     const scalar tol = readScalar(lookup("eps"));
 
@@ -130,10 +130,19 @@ void Foam::sixDOFbodies::solve()
         Info << "Solving 6-DOF for " << names_[bodyI] << " in time"
          << tab << "T = " << runTime_.value() << " s" << endl;
 
+        // Note: set external force and moment needed to initialize the state
+        // of the sixDOFODE to correctly take into account multiple calls per
+        // time step. Using constant force and moment throughout simulation.
+        odes_[bodyI].setExternalForceAndMoment
+        (
+            dimensionedVector(odes_[bodyI].force()),
+            dimensionedVector(odes_[bodyI].moment())
+        );
+
         solvers_[bodyI].solve
         (
+            runTime_.value() - runTime_.deltaT().value(),
             runTime_.value(),
-            runTime_.value() + runTime_.deltaT().value(),
             tol,
             runTime_.deltaT().value()
         );
@@ -141,13 +150,13 @@ void Foam::sixDOFbodies::solve()
 }
 
 
-const Foam::wordList& Foam::sixDOFbodies::names() const
+const Foam::wordList& Foam::sixDOFBodies::names() const
 {
     return names_;
 }
 
 
-const Foam::PtrList<Foam::sixDOFqODE>& Foam::sixDOFbodies::operator()() const
+const Foam::PtrList<Foam::sixDOFODE>& Foam::sixDOFBodies::operator()() const
 {
     return odes_;
 }

src/ODE/sixDOF/sixDOFBodies/sixDOFBodies.H

@@ -22,7 +22,7 @@ License
     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    sixDOFbodies
+    sixDOFBodies
 
 Description
     6-DOF solver for multiple bodies
@@ -31,18 +31,17 @@ Author
     Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
 
 SourceFiles
-    sixDOFbodies.C
+    sixDOFBodies.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef sixDOFbodies_H
+#ifndef sixDOFBodies_H
-#define sixDOFbodies_H
+#define sixDOFBodies_H
 
 #include "foamTime.H"
 #include "IOdictionary.H"
-#include "sixDOFqODE.H"
+#include "sixDOFODE.H"
 #include "ODESolver.H"
-#include "finiteRotation.H"
 #include "primitiveFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -51,10 +50,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class sixDOFbodies Declaration
+                        Class sixDOFBodies Declaration
 \*---------------------------------------------------------------------------*/
 
-class sixDOFbodies
+class sixDOFBodies
 :
     public IOdictionary
 {
@@ -63,10 +62,10 @@ class sixDOFbodies
         //- Reference to time
         const Time& runTime_;
 
-        // Pointer list of Foam::sixDOFqODE objects
+        // Pointer list of Foam::sixDOFODE objects
-        PtrList<sixDOFqODE> odes_;
+        PtrList<sixDOFODE> odes_;
 
-        // Pointer list of Foam::sixDOFqODE solvers
+        // Pointer list of Foam::sixDOFODE solvers
         PtrList<ODESolver> solvers_;
 
         // Name list of solved bodies
@@ -76,10 +75,10 @@ class sixDOFbodies
     // Private Member Functions
 
         //- Disallow default bitwise copy construct
-        sixDOFbodies(const sixDOFbodies&);
+        sixDOFBodies(const sixDOFBodies&);
 
         //- Disallow default bitwise assignment
-        void operator=(const sixDOFbodies&);
+        void operator=(const sixDOFBodies&);
 
         //- Set bodies
         void setBodies();
@@ -90,12 +89,12 @@ public:
     // Constructors
 
         //- Construct from Time
-        sixDOFbodies(const Time& runTime);
+        sixDOFBodies(const Time& runTime);
 
 
     // Destructor
 
-        virtual ~sixDOFbodies()
+        virtual ~sixDOFBodies()
         {}
 
 
@@ -105,7 +104,7 @@ public:
         const wordList& names() const;
 
         //- Return list of bodies
-        const PtrList<sixDOFqODE>& operator()() const;
+        const PtrList<sixDOFODE>& operator()() const;
 
         //- Solve ODEs and update position
         void solve();

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

@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    rotationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("constraintDirection")),
+    alpha_(readScalar(dict.lookup("angularAcceleration"))),
+    inGlobal_(dict.lookup("inGlobalCoordinateSystem"))
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::constantTranslationalAcceleration::"
+            "constantTranslationalAcceleration"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::rotationalConstraint>
+Foam::constantAngularAcceleration::clone() const
+{
+    return autoPtr<rotationalConstraint>
+    (
+        new constantAngularAcceleration(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantAngularAcceleration::~constantAngularAcceleration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::constantAngularAcceleration::matrixContribution
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector&
+) 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;
+}
+
+
+Foam::scalar Foam::constantAngularAcceleration::sourceContribution
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector&
+) const
+{
+    return alpha_;
+}
+
+
+void Foam::constantAngularAcceleration::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+       << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("constraintDirection") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("angularAcceleration") << tab << alpha_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("inGlobalCoordinateSystem") << tab << inGlobal_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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)
+        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_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantAngularAcceleration");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantAngularAcceleration
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~constantAngularAcceleration();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor& toRelative,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar,
+                vector&
+            ) const
+            {
+                // Does nothing: no need to stabilise this constraint
+            }
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rotationalConstraint, 0);
+    defineRunTimeSelectionTable(rotationalConstraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::rotationalConstraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotationalConstraint::~rotationalConstraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::rotationalConstraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::rotationalConstraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::rotationalConstraint> Foam::rotationalConstraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word constraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(constraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "rotationalConstraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\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,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const rotationalConstraint& rc
+)
+{
+    os << rc.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    rc.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const rotationalConstraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,225 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    Notes:
+    1. Constraints are usually used alongside appropriate initial
+       conditions (rotational rate in a given direction),
+    2. According to DAE theory, a stabilisation on orientation/angular velocity
+       is necessary based on "differentiation index" used to obtain the final
+       form of the constraint.
+
+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 "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class rotationalConstraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotationalConstraint
+{
+    // Private Data
+
+        //- Name of the constraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("rotationalConstraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            rotationalConstraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of rotationalConstraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<rotationalConstraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return rotationalConstraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        rotationalConstraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalConstraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected rotationalConstraint
+        static autoPtr<rotationalConstraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~rotationalConstraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+            //- Stabilise the constraint (necessary for constraints with
+            //  differentiation index higher than zero)
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& omega
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const rotationalConstraint& rc
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.C

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("constraintDirection")),
+    a_(readScalar(dict.lookup("translationalAcceleration")))
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::constantTranslationalAcceleration::"
+            "constantTranslationalAcceleration"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint>
+Foam::constantTranslationalAcceleration::clone() const
+{
+    return autoPtr<translationalConstraint>
+    (
+        new constantTranslationalAcceleration(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::constantTranslationalAcceleration::~constantTranslationalAcceleration()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::constantTranslationalAcceleration::matrixContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return dir_;
+}
+
+
+Foam::scalar Foam::constantTranslationalAcceleration::sourceContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return a_;
+}
+
+
+void Foam::constantTranslationalAcceleration::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("constraintDirection") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("translationalAcceleration") << tab << a_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/constantTranslationalAcceleration/constantTranslationalAcceleration.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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)
+        vector dir_;
+
+        //- Constant value of the translational acceleration
+        const scalar a_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantTranslationalAcceleration");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        constantTranslationalAcceleration
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~constantTranslationalAcceleration();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector& 
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar,
+                vector&,
+                vector&
+            ) const
+            {
+                // Does nothing: no need to stabilise this constraint
+            };
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.C

@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "periodicOscillation.H"
+#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(periodicOscillation, 0);
+    addToRunTimeSelectionTable
+    (
+        translationalConstraint,
+        periodicOscillation,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::periodicOscillation::periodicOscillation
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalConstraint(name, dict, sixDOF),
+    dir_(dict.lookup("direction")),
+    a_(readScalar(dict.lookup("motionAmplitude"))),
+    period_(readScalar(dict.lookup("period"))),
+    omega_(mathematicalConstant::twoPi/period_),
+    phi_(readScalar(dict.lookup("phaseShift"))*mathematicalConstant::pi/180.0)
+{
+    // Rescale direction
+    if (mag(dir_) < SMALL)
+    {
+        FatalErrorIn
+        (
+            "Foam::periodicOscillation::periodicOscillation"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Zero direction specified. This is not allowed."
+            << exit(FatalError);
+    }
+    else
+    {
+        dir_ /= mag(dir_);
+    }
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint>
+Foam::periodicOscillation::clone() const
+{
+    return autoPtr<translationalConstraint>
+    (
+        new periodicOscillation(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::periodicOscillation::~periodicOscillation()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::periodicOscillation::matrixContribution
+(
+    const scalar,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return dir_;
+}
+
+
+Foam::scalar Foam::periodicOscillation::sourceContribution
+(
+    const scalar t,
+    const tensor&,
+    const vector&,
+    const vector&
+) const
+{
+    return -a_*sqr(omega_)*(sin(omega_*t + phi_));
+}
+
+
+void Foam::periodicOscillation::stabilise
+(
+    const scalar t,
+    vector& x,
+    vector& u
+) const
+{
+    // Set the displacement according to periodic oscillation
+
+    // First subtract calculated displacement...
+    x -= (x & dir_)*dir_;
+
+    // ... then add the correct displacement
+    x += dir_*sin(omega_*t + phi_);
+
+
+    // Set the velocity according to periodic oscillation
+
+    // First subract calculated velocity...
+    u -= (u & dir_)*dir_;
+
+    // ... then add the correct velocity
+    u += dir_*cos(omega_*t + phi_);
+}
+
+
+void Foam::periodicOscillation::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("direction") << tab << dir_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("motionAmplitude") << tab << a_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("period") << tab << period_
+      << token::END_STATEMENT << nl;
+    os.writeKeyword("phaseShift") << tab
+      << phi_*180.0/mathematicalConstant::pi << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/periodicOscillation/periodicOscillation.H

@@ -0,0 +1,156 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::periodicOscillation
+
+Description
+    Periodic translational motion given by:
+
+    g(vDot, t) = vDot + a*omega^2*sin(omega*t + phi) = 0,
+
+    where a is the amplitude of oscillation, omega radial frequency and phi is
+    the phase shift.
+
+    Note: since this constraint is basically algebraic (dependent on
+    displacement), the differentiation index is two so the stabilisation on
+    displacement and velocity is necessary.
+
+Author
+    Viktor Pandza, FSB Zagreb.  All rights reserved.
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    periodicOscillation.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef periodicOscillation_H
+#define periodicOscillation_H
+
+#include "translationalConstraint.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class periodicOscillation Declaration
+\*---------------------------------------------------------------------------*/
+
+class periodicOscillation
+:
+    public translationalConstraint
+{
+    // Private data
+
+        //- Direction
+        vector dir_;
+    
+        //- Amplitude of the motion
+        const scalar a_;
+
+        //- Period of the motion
+        const scalar period_;
+    
+        //- Radian frequency of the motion
+        const scalar omega_;
+
+        //- Phase shift (in degrees)
+        const scalar phi_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("periodicOscillation");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        periodicOscillation
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~periodicOscillation();
+
+
+    // Member Functions
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar,
+                const tensor&,
+                const vector&,
+                const vector&
+            ) const;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor&,
+                const vector&,
+                const vector& 
+            ) const;
+
+            //- Stabilise the constraint
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& x,
+                vector& u
+            ) const;
+
+
+    // I-O Functions
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(translationalConstraint, 0);
+    defineRunTimeSelectionTable(translationalConstraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::translationalConstraint::translationalConstraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::translationalConstraint::~translationalConstraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::translationalConstraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::translationalConstraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::translationalConstraint> Foam::translationalConstraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    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 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,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const translationalConstraint& tc
+)
+{
+    os << tc.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    tc.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const translationalConstraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraints/translationalConstraints/translationalConstraint/translationalConstraint.H

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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.
+
+    Notes:
+    1. Constraints are usually used alongside appropriate initial
+       conditions (velocity in a given direction),
+    2. According to DAE theory, a stabilisation on displacement/velocity is
+       necessary based on "differentiation index" used to obtain the final form
+       of the constraint.
+
+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 "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class translationalConstraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class translationalConstraint
+{
+    // Private Data
+
+        //- Name of the constraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("translationalConstraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            translationalConstraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of translationalConstraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<translationalConstraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return translationalConstraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        translationalConstraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalConstraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected translationalConstraint
+        static autoPtr<translationalConstraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~translationalConstraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Constraint specific functions
+
+            //- Return matrix contribution defined by constraint, f(t)
+            virtual vector matrixContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+            //- Return source contribution defined by constraint, a(t)
+            virtual scalar sourceContribution
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+            //- Stabilise the constraint (necessary for constraints with
+            //  differentiation index higher than zero)
+            virtual void stabilise
+            (
+                const scalar t,
+                vector& x,
+                vector& u
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const translationalConstraint& tc
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/constraintsAndRestraints.H

@@ -0,0 +1,10 @@
+#ifndef constraintsAndRestraints_H
+#define constraintsAndRestraints_H
+
+#include "translationalConstraint.H"
+#include "rotationalConstraint.H"
+
+#include "translationalRestraint.H"
+#include "rotationalRestraint.H"
+
+#endif

src/ODE/sixDOF/sixDOFODE/newSixDOFODE.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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
+
+Class
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "objectRegistry.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::sixDOFODE> Foam::sixDOFODE::New(const IOobject& io)
+{
+    word sixDOFODETypeName;
+
+    // Get object registry
+    const objectRegistry& database = io.db();
+
+    // Check whether the dictionary is in the database
+    if (database.foundObject<IOdictionary>(io.name()))
+    {
+        sixDOFODETypeName =
+            word
+            (
+                database.lookupObject<IOdictionary>(io.name()).lookup("type")
+            );
+    }
+    else
+    {
+        sixDOFODETypeName = word(IOdictionary(io).lookup("type"));
+    }
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(sixDOFODETypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "sixDOFODE::New(const IOobject& io)"
+        )   << "Unknown sixDOFODE " << sixDOFODETypeName
+            << endl << endl
+            << "Valid sixDOFODE types are:" << endl
+            << dictionaryConstructorTablePtr_->toc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<sixDOFODE>(cstrIter()(io));
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "angularDamper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(angularDamper, 0);
+    addToRunTimeSelectionTable
+    (
+        rotationalRestraint,
+        angularDamper,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::angularDamper::angularDamper
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    rotationalRestraint(name, dict, sixDOF),
+    angDampingCoeffs_(dict.lookup("angularDamping")),
+    inGlobal_(dict.lookup("inGlobalCoordinateSystem"))
+{}
+
+
+Foam::autoPtr<Foam::rotationalRestraint>
+Foam::angularDamper::clone() const
+{
+    return autoPtr<rotationalRestraint>
+    (
+        new angularDamper(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::angularDamper::~angularDamper()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::angularDamper::restrainingMoment
+(
+    const scalar,
+    const tensor& toRelative,
+    const vector& omega
+) const
+{
+    vector rm;
+
+    if (inGlobal_)
+    {
+        // Restraint given in global (inertial) coordinate system, transform it
+        // to local
+        rm = (toRelative & angDampingCoeffs_) & omega;
+    }
+    else
+    {
+        // Restraint already in local (body) coordinate system
+        rm = angDampingCoeffs_ & omega;
+    }
+
+    return -rm;
+}
+
+
+void Foam::angularDamper::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("angularDamping") << tab << angDampingCoeffs_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("inGlobalCoordinateSystem") << tab << inGlobal_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/angularDamper/angularDamper.H

@@ -0,0 +1,119 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::angularDamper
+
+Description
+    Rotational restraint corresponding to angular damper defined by angular
+    damping coefficients.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    angularDamper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef angularDamper_H
+#define angularDamper_H
+
+#include "rotationalRestraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class angularDamper Declaration
+\*---------------------------------------------------------------------------*/
+
+class angularDamper
+:
+    public rotationalRestraint
+{
+    // Private Data
+
+        //- Angular damping coefficients
+        diagTensor angDampingCoeffs_;
+
+        //- Whether the damper is applied in global or local c. s.
+        Switch inGlobal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("angularDamper");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        angularDamper
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalRestraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~angularDamper();
+
+
+    // Member Functions
+
+        // Restraint specific functions
+
+            //- Return restraining moment (in the global coordinate system)
+            virtual vector restrainingMoment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "rotationalRestraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(rotationalRestraint, 0);
+    defineRunTimeSelectionTable(rotationalRestraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::rotationalRestraint::rotationalRestraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::rotationalRestraint::~rotationalRestraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::rotationalRestraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::rotationalRestraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::rotationalRestraint> Foam::rotationalRestraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word restraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(restraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "rotationalRestraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Unknown translation restraint type: " << restraintType
+            << endl << endl
+            << "Valid translation restraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<rotationalRestraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const rotationalRestraint& rr
+)
+{
+    os << rr.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    rr.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const rotationalRestraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/rotationalRestraints/rotationalRestraint/rotationalRestraint.H

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::rotationalRestraint
+
+Description
+    Abstract base class containing interface for rotational restraints used
+    within sixDOFODE classes.
+
+    Interface provides restraining moment to the rotational equations of
+    motion via moment() member function.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    rotationalRestraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef rotationalRestraint_H
+#define rotationalRestraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class rotationalRestraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class rotationalRestraint
+{
+    // Private Data
+
+        //- Name of the restraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("rotationalRestraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            rotationalRestraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of rotationalRestraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<rotationalRestraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return rotationalRestraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        rotationalRestraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<rotationalRestraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected rotationalRestraint
+        static autoPtr<rotationalRestraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~rotationalRestraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Restraint specific functions
+
+            //- Return restraining moment (in the local coordinate system)
+            virtual vector restrainingMoment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const rotationalRestraint& rr
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.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 "linearSpringDamper.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(linearSpringDamper, 0);
+    addToRunTimeSelectionTable
+    (
+        translationalRestraint,
+        linearSpringDamper,
+        word
+    );
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::linearSpringDamper::linearSpringDamper
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    translationalRestraint(name, dict, sixDOF),
+    linSpringCoeffs_(dict.lookup("linearSpring")),
+    linDampingCoeffs_(dict.lookup("linearDamping"))
+{}
+
+
+Foam::autoPtr<Foam::translationalRestraint>
+Foam::linearSpringDamper::clone() const
+{
+    return autoPtr<translationalRestraint>
+    (
+        new linearSpringDamper(*this)
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::linearSpringDamper::~linearSpringDamper()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::linearSpringDamper::restrainingForce
+(
+    const scalar,
+    const tensor&,
+    const vector& x,
+    const vector& u
+) const
+{
+    return - (linSpringCoeffs_ & x) - (linDampingCoeffs_ & u);
+}
+
+
+void Foam::linearSpringDamper::write(Ostream& os) const
+{
+    os.writeKeyword("type") << tab << type()
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("linearSpring") << tab << linSpringCoeffs_
+       << token::END_STATEMENT << nl;
+    os.writeKeyword("linearDamping") << tab << linDampingCoeffs_
+       << token::END_STATEMENT << endl;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/linearSpringDamper/linearSpringDamper.H

@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::linearSpringDamper
+
+Description
+    Translational restraint corresponding to linear spring and linear damper
+    defined by linear spring coefficients and linear damping coefficients.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    linearSpringDamper.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linearSpringDamper_H
+#define linearSpringDamper_H
+
+#include "translationalRestraint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class linearSpringDamper Declaration
+\*---------------------------------------------------------------------------*/
+
+class linearSpringDamper
+:
+    public translationalRestraint
+{
+    // Private Data
+
+        //- Linear spring coefficients
+        diagTensor linSpringCoeffs_;
+    
+        //- Linear damping coefficients
+        diagTensor linDampingCoeffs_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("linearSpringDamper");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        linearSpringDamper
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalRestraint> clone() const;
+
+
+    // Destructor
+
+        virtual ~linearSpringDamper();
+
+
+    // Member Functions
+
+        // Restraint specific functions
+
+            //- Return restraining force (in the global coordinate system)
+            virtual vector restrainingForce
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "translationalRestraint.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(translationalRestraint, 0);
+    defineRunTimeSelectionTable(translationalRestraint, word);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::translationalRestraint::translationalRestraint
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+:
+    name_(name),
+    sixDOF_(sixDOF)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::translationalRestraint::~translationalRestraint()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::word& Foam::translationalRestraint::name() const
+{
+    return name_;
+}
+
+
+const Foam::sixDOFODE& Foam::translationalRestraint::sixDOF() const
+{
+    return sixDOF_;
+}
+
+
+Foam::autoPtr<Foam::translationalRestraint> Foam::translationalRestraint::New
+(
+    const word& name,
+    const dictionary& dict,
+    const sixDOFODE& sixDOF
+)
+{
+    const word restraintType(dict.lookup("type"));
+
+    wordConstructorTable::iterator cstrIter =
+        wordConstructorTablePtr_->find(restraintType);
+
+    if (cstrIter == wordConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "translationalRestraint::New"
+            "\n("
+            "\n    const word& name,"
+            "\n    const dictionary& dict"
+            "\n)"
+        )   << "Unknown translation restraint type: " << restraintType
+            << endl << endl
+            << "Valid translation restraint types are: " << endl
+            << wordConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return autoPtr<translationalRestraint>
+    (
+        cstrIter()
+        (
+            name,
+            dict,
+            sixDOF
+        )
+    );
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const translationalRestraint& tr
+)
+{
+    os << tr.name_ << nl << token::BEGIN_BLOCK << nl;
+
+    tr.write(os);
+
+    os << token::END_BLOCK << endl;
+
+    os.check("Ostream& operator<<(Ostream&, const translationalRestraint&");
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/restraints/translationalRestraints/translationalRestraint/translationalRestraint.H

@@ -0,0 +1,192 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::translationalRestraint
+
+Description
+    Abstract base class containing interface for translational restraints used
+    within sixDOFODE classes.
+
+    Interface provides restraining forces to the translational equations of
+    motion via force() member function.
+
+Author
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    translationalRestraint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef translationalRestraint_H
+#define translationalRestraint_H
+
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "autoPtr.H"
+#include "dimensionedTypes.H"
+#include "dictionary.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class sixDOFODE;
+
+/*---------------------------------------------------------------------------*\
+                   Class translationalRestraint Declaration
+\*---------------------------------------------------------------------------*/
+
+class translationalRestraint
+{
+    // Private Data
+
+        //- Name of the restraint
+        word name_;
+
+        //- Reference to underlying sixDOFODE
+        const sixDOFODE& sixDOF_;
+
+public:
+
+    //- Runtime type information
+    TypeName("translationalRestraint");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            translationalRestraint,
+            word,
+            (
+                const word& name,
+                const dictionary& dict,
+                const sixDOFODE& sixDOF
+            ),
+            (name, dict, sixDOF)
+        );
+
+
+        //- Class used for the read-construction of
+        //  PtrLists of translationalRestraint
+        class iNew
+        {
+            const sixDOFODE& sixDOF_;
+
+        public:
+
+            iNew(const sixDOFODE& sixDOF)
+            :
+                sixDOF_(sixDOF)
+            {}
+
+            autoPtr<translationalRestraint> operator()(Istream& is) const
+            {
+                word name(is);
+                dictionary dict(is);
+                return translationalRestraint::New(name, dict, sixDOF_);
+            }
+        };
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        translationalRestraint
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<translationalRestraint> clone() const = 0;
+
+
+    // Selectors
+
+        //- Return a reference to the selected translationalRestraint
+        static autoPtr<translationalRestraint> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const sixDOFODE& sixDOF
+        );
+
+
+    // Destructor
+
+        virtual ~translationalRestraint();
+
+
+    // Member Functions
+
+        // Access functions
+
+            //- Return name
+            const word& name() const;
+    
+            //- Return underlying sixDOFODE object
+            const sixDOFODE& sixDOF() const;
+
+
+        // Restraint specific functions
+
+            //- Return restraining force (in the global coordinate system)
+            virtual vector restrainingForce
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u 
+            ) const = 0;
+
+
+    // I-O Functions and Operators
+
+        //- Virtual write function
+        virtual void write(Ostream& os) const = 0;
+
+        //- Ostream operator implemented in terms of write operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const translationalRestraint& tr
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODE.C

@@ -0,0 +1,543 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+#include "objectRegistry.H"
+#include "sixDOFODE.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(sixDOFODE, 0);
+defineRunTimeSelectionTable(sixDOFODE, dictionary);
+
+}
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+void Foam::sixDOFODE::updateRelaxFactors()
+{
+    // Calculate translational relax factor
+    const scalar saveOldRelFacT = oldRelaxFactorT_;
+    oldRelaxFactorT_ = relaxFactorT_;
+
+    if(magSqr(A_[0] - An_[1] - A_[1] - An_[2]) > SMALL)
+    {
+        relaxFactorT_ = saveOldRelFacT + (saveOldRelFacT - 1)*
+        (
+            (A_[1] - An_[2])
+          & (A_[0] - An_[1] - A_[1] - An_[2])
+        )/
+        magSqr(A_[0] - An_[1] - A_[1] - An_[2]);
+    }
+    else
+    {
+        relaxFactorT_ = minRelaxFactor_;
+    }
+
+    // Bound the relaxation factor for stability
+    if (relaxFactorT_ > maxRelaxFactor_)
+    {
+        relaxFactorT_ = maxRelaxFactor_;
+    }
+    else if (relaxFactorT_ < minRelaxFactor_)
+    {
+        relaxFactorT_ = minRelaxFactor_;
+    }
+
+    // Calculate rotational relax factor
+    const scalar saveOldRelFacR = oldRelaxFactorR_;
+    oldRelaxFactorR_ = relaxFactorR_;
+
+    if
+    (
+        magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2])
+      > SMALL
+    )
+    {
+        relaxFactorR_ = saveOldRelFacR
+          + (saveOldRelFacR - 1)*
+            (
+                (OmegaDot_[1] - OmegaDotn_[2])
+              & (OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2])
+            )/
+            magSqr(OmegaDot_[0] - OmegaDotn_[1] - OmegaDot_[1] - OmegaDotn_[2]);
+    }
+    else
+    {
+        relaxFactorR_ = minRelaxFactor_;
+    }
+
+    // Bound the relaxation factor for stability
+    if (relaxFactorR_ > maxRelaxFactor_)
+    {
+        relaxFactorR_ = maxRelaxFactor_;
+    }
+    else if (relaxFactorR_ < minRelaxFactor_)
+    {
+        relaxFactorR_ = minRelaxFactor_;
+    }
+}
+
+
+void Foam::sixDOFODE::relaxAcceleration()
+{
+    if (mag(minRelaxFactor_ - maxRelaxFactor_) < SMALL)
+    {
+       // Fixed relaxation
+       relaxFactorT_ = minRelaxFactor_;
+       relaxFactorR_ = minRelaxFactor_;
+    }
+    else
+    {
+        // Use Aitkens relaxation
+
+        // Update Aitkens relaxation factor
+        updateRelaxFactors();
+
+        // Update non relaxed accelerations
+        An_[1] = An_[2];
+        An_[2] = (force_/mass_).value();
+        OmegaDotn_[1] = OmegaDotn_[2];
+        OmegaDotn_[2] = (inv(momentOfInertia_) & moment_).value();
+    }
+
+    const vector Aold = A_[2];
+    const vector OmegaDotold = OmegaDot_[2];
+
+    force_.value() =
+        Aold*mass_.value()
+      + relaxFactorT_*(force_.value() - Aold*mass_.value());
+
+    moment_.value() =
+        (momentOfInertia_.value() & OmegaDotold)
+      + relaxFactorR_*
+        (
+            moment_.value()
+          - (momentOfInertia_.value() & OmegaDotold)
+        );
+
+    // Update relaxed old accelerations
+    A_[0] = A_[1];
+    A_[1] = A_[2];
+    A_[2] = (force_/mass_).value();
+    OmegaDot_[0] = OmegaDot_[1];
+    OmegaDot_[1] = OmegaDot_[2];
+    OmegaDot_[2] = (inv(momentOfInertia_) & moment_).value();
+}
+
+
+void Foam::sixDOFODE::initState()
+{
+    // Get time index
+    const label timeIndex = dict().time().timeIndex();
+
+    if (curTimeIndex_ < timeIndex)
+    {
+        // First call in this time index, store data
+        oldStatePtr_ = this->clone(dict().name() + "_0");
+
+        Info<< "First 6DOF solution within a time step, storing old data..."
+            << endl;
+    }
+    else
+    {
+        // Multiple calls in this time index, reset this data
+        this->setState(oldStatePtr_());
+
+        Info<< "Repeated 6DOF solution within a time step, restoring data..."
+            << endl;
+    }
+
+    // Update local time index
+    curTimeIndex_ = timeIndex;
+}
+
+
+// * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * * //
+
+void Foam::sixDOFODE::setState(const sixDOFODE& sd)
+{
+    // Set state does not copy AList_, AOld_, relaxFactor_ and relaxFactorOld_.
+    // In case of multiple updates, overwriting Aitkens relaxation parameters
+    // would invalidate the underrelaxation.  IG, 5/May/2016
+    mass_ = sd.mass_;
+    momentOfInertia_ = sd.momentOfInertia_;
+
+    Xequilibrium_ = sd.Xequilibrium_;
+
+    force_ = sd.force_;
+    moment_ = sd.moment_;
+
+
+    // Copy constraints
+    
+    // Translational constraints
+    translationalConstraints_.setSize(sd.translationalConstraints_.size());
+    forAll(sd.translationalConstraints_, trI)
+    {
+        translationalConstraints_.set
+        (
+            trI,
+            sd.translationalConstraints_[trI].clone()
+        );
+    }
+
+    // Rotational constraints
+    rotationalConstraints_.setSize(sd.rotationalConstraints_.size());
+    forAll(sd.rotationalConstraints_, rrI)
+    {
+        rotationalConstraints_.set
+        (
+            rrI,
+            sd.rotationalConstraints_[rrI].clone()
+        );
+    }
+
+
+    // Copy restraints
+
+    // Translational restraints
+    translationalRestraints_.setSize(sd.translationalRestraints_.size());
+    forAll(sd.translationalRestraints_, trI)
+    {
+        translationalRestraints_.set
+        (
+            trI,
+            sd.translationalRestraints_[trI].clone()
+        );
+    }
+
+    // Rotational restraints
+    rotationalRestraints_.setSize(sd.rotationalRestraints_.size());
+    forAll(sd.rotationalRestraints_, rrI)
+    {
+        rotationalRestraints_.set
+        (
+            rrI,
+            sd.rotationalRestraints_[rrI].clone()
+        );
+    }
+}
+
+
+Foam::dimensionedVector Foam::sixDOFODE::force
+(
+    const scalar t,
+    const tensor& toRelative,
+    const vector& x,
+    const vector& u
+) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Calculate restraining force
+    dimensionedVector rForce("zero", dimForce, vector::zero);
+
+    forAll(translationalRestraints_, trI)
+    {
+        rForce.value() += translationalRestraints_[trI].restrainingForce
+        (
+            t, // Time
+            toRelative, // Transformation tensor
+            x, // Position in the global c.s.
+            u // Velocity in the global c.s.
+        );
+    }
+
+    // Return linearly interpolated external force with restraining force
+    return (alpha*oldStatePtr_->force() + (1 - alpha)*force()) + rForce;
+}
+
+
+Foam::dimensionedVector Foam::sixDOFODE::moment
+(
+    const scalar t,
+    const tensor& toRelative,
+    const vector& omega
+) const
+{
+    // Get ODE step fraction
+    const scalar alpha = odeStepFraction(t);
+
+    // Calculate restraining moment
+    dimensionedVector rMoment("zero", dimForce*dimLength, vector::zero);
+
+    forAll(rotationalRestraints_, rrI)
+    {
+        rMoment.value() += rotationalRestraints_[rrI].restrainingMoment
+        (
+            t, // Time
+            toRelative, // Transformation tensor
+            omega // Angular velocity in local c.s.
+        );
+    }
+
+    // Return linearly interpolated external moment with restraining moment
+    return (alpha*oldStatePtr_->moment() + (1 - alpha)*moment() + rMoment);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::sixDOFODE(const IOobject& io)
+:
+    ODE(),
+    dict_(io, *this),
+
+    mass_(dict_.lookup("mass")),
+    momentOfInertia_(dict_.lookup("momentOfInertia")),
+
+    Xequilibrium_(dict_.lookup("equilibriumPosition")),
+
+    aitkensRelaxation_
+    (
+        dict_.lookupOrDefault<Switch>("useAitkensRelaxation", false)
+    ),
+    minRelaxFactor_(dict_.lookupOrDefault<scalar>("minRelaxFactor", 0.1)),
+    maxRelaxFactor_(dict_.lookupOrDefault<scalar>("maxRelaxFactor", 0.5)),
+
+    relaxFactorT_(1.0),
+    relaxFactorR_(1.0),
+    oldRelaxFactorT_(1.0),
+    oldRelaxFactorR_(1.0),
+
+    A_(3, vector::zero),
+    OmegaDot_(3, vector::zero),
+    An_(3, vector::zero),
+    OmegaDotn_(3, vector::zero),
+
+    force_(dict_.lookup("force")),
+    moment_(dict_.lookup("moment")),
+
+    curTimeIndex_(-1),
+    oldStatePtr_(),
+
+    translationalConstraints_(),
+    rotationalConstraints_(),
+
+    translationalRestraints_(),
+    rotationalRestraints_()
+{
+    // Sanity checks
+    if (mass_.value() < SMALL)
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Zero or negative mass detected: " << mass_.value()
+            << nl << "Please check " << dict_.name() << "dictionary."
+            << exit(FatalIOError);
+    }
+
+    if (cmptMin(momentOfInertia_.value()) < SMALL)
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Zero or negative moment of inertia detected: "
+            << momentOfInertia_.value()
+            << nl << "Please check " << dict_.name() << "dictionary."
+            << exit(FatalIOError);
+    }
+
+    if
+    (
+        (minRelaxFactor_ < SMALL)
+     || (maxRelaxFactor_ > 1.0)
+     || ((maxRelaxFactor_ - minRelaxFactor_) < 0)
+    )
+    {
+        FatalIOErrorIn
+        (
+            "sixDOFODE::sixDOFODE(const IOobject& io)",
+            dict_
+        )   << "Invalid minRelaxFactor and maxRelaxFactor specified."
+            << nl << "Please use values within 0 and 1."
+            << exit(FatalIOError);
+    }
+
+    // Read and construct constraints and restraints
+
+    // Read translation constraints if they are present
+    if (dict().found("translationalConstraints"))
+    {
+        PtrList<translationalConstraint> tcList
+        (
+            dict().lookup("translationalConstraints"),
+            translationalConstraint::iNew(*this)
+        );
+        translationalConstraints_.transfer(tcList);
+    }
+
+    // Read rotation constraints if they are present
+    if (dict().found("rotationalConstraints"))
+    {
+        PtrList<rotationalConstraint> rcList
+        (
+            dict().lookup("rotationalConstraints"),
+            rotationalConstraint::iNew(*this)
+        );
+        rotationalConstraints_.transfer(rcList);
+    }
+
+    // Read translation restraints if they are present
+    if (dict().found("translationalRestraints"))
+    {
+        PtrList<translationalRestraint> tcList
+        (
+            dict().lookup("translationalRestraints"),
+            translationalRestraint::iNew(*this)
+        );
+        translationalRestraints_.transfer(tcList);
+    }
+
+    // Read rotation restraints if they are present
+    if (dict().found("rotationalRestraints"))
+    {
+        PtrList<rotationalRestraint> rcList
+        (
+            dict().lookup("rotationalRestraints"),
+            rotationalRestraint::iNew(*this)
+        );
+        rotationalRestraints_.transfer(rcList);
+    }
+}
+
+
+Foam::sixDOFODE::sixDOFODE(const word& name, const sixDOFODE& sd)
+:
+    ODE(),
+    dict_(sd.dict_),
+
+    mass_(sd.mass_),
+    momentOfInertia_(sd.momentOfInertia_),
+
+    Xequilibrium_(sd.Xequilibrium_),
+
+    aitkensRelaxation_(sd.aitkensRelaxation_),
+    minRelaxFactor_(sd.minRelaxFactor_),
+    maxRelaxFactor_(sd.maxRelaxFactor_),
+
+    relaxFactorT_(1.0),
+    relaxFactorR_(1.0),
+    oldRelaxFactorT_(1.0),
+    oldRelaxFactorR_(1.0),
+
+    A_(3, vector::zero),
+    OmegaDot_(3, vector::zero),
+    An_(3, vector::zero),
+    OmegaDotn_(3, vector::zero),
+
+    force_(sd.force_),
+    moment_(sd.moment_),
+
+    curTimeIndex_(sd.curTimeIndex_),
+    oldStatePtr_(),
+
+    translationalConstraints_(sd.translationalConstraints_),
+    rotationalConstraints_(sd.rotationalConstraints_),
+
+    translationalRestraints_(sd.translationalRestraints_),
+    rotationalRestraints_(sd.rotationalRestraints_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sixDOFODE::~sixDOFODE()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::OutputControlDictionary<Foam::sixDOFODE>&
+Foam::sixDOFODE::dict() const
+{
+    return dict_;
+}
+
+
+bool Foam::sixDOFODE::writeData(Ostream& os) const
+{
+    os.writeKeyword("mass") << tab << mass_ << token::END_STATEMENT << nl;
+    os.writeKeyword("momentOfInertia") << tab << momentOfInertia_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("equilibriumPosition") << tab << Xequilibrium_
+        << token::END_STATEMENT << nl;
+
+    os.writeKeyword("useAitkensRelaxation") << tab << aitkensRelaxation_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("minRelaxFactor") << tab << minRelaxFactor_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("maxRelaxFactor") << tab << maxRelaxFactor_
+        << token::END_STATEMENT << nl << nl;
+
+    os.writeKeyword("force") << tab << force_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("moment") << tab << moment_
+        << token::END_STATEMENT << nl << nl;
+
+    if (!translationalConstraints_.empty())
+    {
+        os.writeKeyword("translationalConstraints")
+            << translationalConstraints_
+            << token::END_STATEMENT << nl << endl;
+    }
+
+    if (!rotationalConstraints_.empty())
+    {
+        os.writeKeyword("rotationalConstraints")
+            << rotationalConstraints_
+            << token::END_STATEMENT << endl;
+    }
+
+    if (!translationalRestraints_.empty())
+    {
+        os.writeKeyword("translationalRestraints")
+            << translationalRestraints_
+            << token::END_STATEMENT << nl << endl;
+    }
+
+    if (!rotationalRestraints_.empty())
+    {
+        os.writeKeyword("rotationalRestraints")
+            << rotationalRestraints_
+            << token::END_STATEMENT << endl;
+    }
+
+    return os.good();
+}
+
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODE.H

@@ -0,0 +1,444 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    sixDOFODE
+
+Description
+    Abstract base class for six-degrees-of-freedom (6DOF) ordinary differential
+    equations with necessary interface for two-way coupling with CFD solvers.
+
+    Holds list of translational and rotational constraints and restraings to be
+    used in derived classes.
+
+Author
+    Dubravko Matijasevic, FSB Zagreb.  All rights reserved.
+
+Reorganisation by
+    Vuko Vukcevic, FSB Zagreb.  All rights reserved.
+    Inno Gatin, FSB Zagreb. All rights reserved.
+    Hrvoje Jasak, FSB Zagreb.  All rights reserved.
+
+SourceFiles
+    sixDOFODEI.H
+    sixDOFODE.C
+    newSixDOFODE.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sixDOFODE_H
+#define sixDOFODE_H
+
+#include "ODE.H"
+#include "dimensionedTypes.H"
+#include "autoPtr.H"
+#include "OutputControlDictionary.H"
+#include "runTimeSelectionTables.H"
+#include "Switch.H"
+#include "foamTime.H"
+#include "constraintsAndRestraints.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class sixDOFODE Declaration
+\*---------------------------------------------------------------------------*/
+
+class sixDOFODE
+:
+    public ODE
+{
+    // Private data
+
+        //- Dictionary object controlling I/O for sixDOFODE
+        const OutputControlDictionary<sixDOFODE> dict_;
+
+
+        // Body data
+
+            //- Mass
+            dimensionedScalar mass_;
+
+            //- Rotational moment of inertia around centre of mass
+            //  in body (relative) coordinates - aligned with main axes
+            dimensionedDiagTensor momentOfInertia_;
+
+
+        // Equilibrium position
+
+            //- Spring equilibrium position for translation
+            dimensionedVector Xequilibrium_;
+
+
+        // Aitkens relaxation data members
+
+            //- Switch to control whether to use Aitkens relaxation
+            Switch aitkensRelaxation_;
+
+            //- Minimum acceleration relaxation factor (default 0.1)
+            const scalar minRelaxFactor_;
+
+            //- Maximum acceleration relaxation factor (default 0.5)
+            const scalar maxRelaxFactor_;
+
+            //- Translational relaxation factor
+            scalar relaxFactorT_;
+
+            //- Rotational relaxation factor
+            scalar relaxFactorR_;
+
+            //- Old translational relaxation factor
+            scalar oldRelaxFactorT_;
+
+            //- Old rotational relaxation factor
+            scalar oldRelaxFactorR_;
+
+            //- Accelerations from previous iterations
+            //  A_[2] is the old value, A_[1] old old, A_[0] old old old
+            List<vector> A_;
+            List<vector> OmegaDot_;
+
+            //- Previos iteration non relaxed accelerations
+            List<vector> An_;
+            List<vector> OmegaDotn_;
+
+
+        // External forces and moments
+
+            //- Force driving the motion in global (inertial) coord. sys.
+            dimensionedVector force_;
+
+            //- Moment driving the motion in global (inertial) coord. sys.
+            dimensionedVector moment_;
+
+
+        // Private data used to control multiple ODE updates per time step
+        // Note: Before solving the ODE from the top level, we will store the
+        // previous state if this is the first update in a given time step. The
+        // state will be kept as the pointer to sixDOFODE, which is not optimal
+        // because we are keeping track of some unnecessary data. We could use
+        // more elegant and efficient code design. VV, 1/Mar/2017.
+
+            //- Local time index
+            label curTimeIndex_;
+
+            //- Pointer to the sixDOFODE object carrying old state
+            autoPtr<sixDOFODE> oldStatePtr_;
+
+
+        // Motion constraints
+
+            //- List of translational constraints
+            PtrList<translationalConstraint> translationalConstraints_;
+
+            //- List of rotational constraints
+            PtrList<rotationalConstraint> rotationalConstraints_;
+
+
+        // Motion restraints
+
+            //- List of translational restraints
+            PtrList<translationalRestraint> translationalRestraints_;
+
+            //- List of rotational restraints
+            PtrList<rotationalRestraint> rotationalRestraints_;
+
+
+    // Private Member Functions
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            sixDOFODE(const sixDOFODE&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const sixDOFODE&);
+
+
+        // Aitkens relaxation helper function
+
+            //- Update Aitkens relaxation parameters
+            void updateRelaxFactors();
+
+            //- Relax the force (acceleration) using Aitkens or fixed relaxation
+            void relaxAcceleration();
+
+
+        // Helper function for controlling multiple ODE solver calls per
+        // time-step
+
+            //- Initialise ODE before setting external forces/moments and
+            //  solving
+            void initState();
+
+            //- Calculate current ODE step fraction given time from ODE
+            inline scalar odeStepFraction(const scalar odeTime) const;
+
+
+protected:
+
+    // Protected Member Functions
+
+        // Non-access control for setting state variables
+
+            //- Set ODE parameters from another ODE
+            virtual void setState(const sixDOFODE&);
+
+
+        // Get external force and moment (used during the solution process)
+
+            //- Return external force with restraints for given ODE state
+            dimensionedVector force
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& x,
+                const vector& u
+            ) const;
+
+            //- Return external moment with restraints given ODE state
+            dimensionedVector moment
+            (
+                const scalar t,
+                const tensor& toRelative,
+                const vector& omega
+            ) const;
+
+
+public:
+
+    //- Run-time type information
+    TypeName("sixDOFODE");
+
+
+    // Declare run-time constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            sixDOFODE,
+            dictionary,
+            (const IOobject& io),
+            (io)
+        );
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        sixDOFODE(const IOobject& io);
+
+        //- Copy construct given new name
+        sixDOFODE(const word& name, const sixDOFODE& sd);
+
+        //- Return a clone, changing the name
+        virtual autoPtr<sixDOFODE> clone(const word& name) const = 0;
+
+
+    // Selectors
+
+        //- Return autoPtr to the selected sixDOFODE
+        static autoPtr<sixDOFODE> New(const IOobject& io);
+
+
+    // Destructor
+
+        virtual ~sixDOFODE();
+
+
+    // Member Functions
+
+        // Access to common data
+
+            //- Return write controlled dictionary
+            const OutputControlDictionary<sixDOFODE>& dict() const;
+
+            //- Return mass
+            inline const dimensionedScalar& mass() const;
+
+            //- Return access to mass
+            inline dimensionedScalar& mass();
+
+            //- Return moment of inertia
+            inline const dimensionedDiagTensor& momentOfInertia() const;
+
+            //- Return access to moment of inertia
+            inline dimensionedDiagTensor& momentOfInertia();
+
+            //- Return equilibrium position of origin
+            inline const dimensionedVector& Xequilibrium() const;
+
+            //- Return access to equilibrium position of origin
+            inline dimensionedVector& Xequilibrium();
+
+            //- Return linear spring coefficient
+            inline const dimensionedDiagTensor& linSpringCoeffs() const;
+
+            //- Return linear damping coefficient
+            inline const dimensionedDiagTensor& linDampingCoeffs() const;
+
+
+        // Access to forces and moments
+
+            //- Return force in global (inertial) coord. sys.
+            inline const dimensionedVector& force() const;
+
+            //- Return moment in global (inertial) coord. sys.
+            inline const dimensionedVector& moment() const;
+
+            //- Set external force and moment
+            inline void setExternalForceAndMoment
+            (
+                const dimensionedVector& externalForce,
+                const dimensionedVector& externalMoment
+            );
+
+            //- Initialize force and moment for the first time step
+            inline void initExternalForceAndMoment
+            (
+                const dimensionedVector& externalForce,
+                const dimensionedVector& externalMoment
+            );
+
+
+        // Access to motion constraints
+
+            //- Return const reference to translational constraints
+            inline const PtrList<translationalConstraint>&
+            translationalConstraints() const;
+
+            //- Return const reference to rotational constraints
+            inline const PtrList<rotationalConstraint>&
+            rotationalConstraints() const;
+
+
+        // Access to motion restraints
+
+            //- Return const reference to translational restraints
+            inline const PtrList<translationalRestraint>&
+            translationalRestraints() const;
+
+            //- Return const reference to rotational restraints
+            inline const PtrList<rotationalRestraint>&
+            rotationalRestraints() const;
+
+
+        // Virtual interface for 6DOF motion state
+
+            // Variables in relative coordinate system
+
+                //- Return displacement in translated coordinate system
+                //  relative to spring equilibrium
+                virtual const dimensionedVector& Xrel() const = 0;
+
+                //- Return rotational velocity in relative coordinate system
+                virtual const dimensionedVector& omega() const = 0;
+
+
+            // Displacement and velocity in the absolute coordinate system
+
+                //- Return position of origin in absolute coordinate system
+                virtual dimensionedVector X() const = 0;
+
+                //- Return velocity of origin
+                virtual const dimensionedVector& U() const = 0;
+
+                //- Return average velocity of origin (evaluated at midstep)
+                virtual const dimensionedVector& Uaverage() const = 0;
+
+
+            // Average motion per time-step
+
+                //- Return average rotational velocity in relative coordinate
+                //  system (evaluated at midstep)
+                virtual const dimensionedVector& omegaAverage() const = 0;
+
+
+            // Rotations
+
+                //- Return rotation tensor to relative coordinate system
+                virtual tensor toRelative() const = 0;
+
+                //- Return rotation tensor to absolute coordinate system
+                virtual tensor toAbsolute() const = 0;
+
+                //- Return transformation tensor between new and previous
+                //  rotation
+                virtual const tensor& rotIncrementTensor() const = 0;
+
+
+            // ODE parameters
+
+                //- Return number of equations
+                virtual label nEqns() const = 0;
+
+                //- Return access to coefficients
+                virtual scalarField& coeffs() = 0;
+
+                //- Return reference to coefficients
+                virtual const scalarField& coeffs() const = 0;
+
+                //- Evaluate derivatives
+                virtual void derivatives
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dydx
+                ) const = 0;
+
+                //- Evaluate Jacobian
+                virtual void jacobian
+                (
+                    const scalar x,
+                    const scalarField& y,
+                    scalarField& dfdx,
+                    scalarSquareMatrix& dfdy
+                ) const = 0;
+
+                //- Update ODE after the solution, advancing by delta
+                virtual void update(const scalar delta) = 0;
+
+
+        // Write control
+
+            //- writeData member function required by regIOobject
+            virtual bool writeData(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sixDOFODEI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/ODE/sixDOF/sixDOFODE/sixDOFODEI.H

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+    sixDOFODE
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
+
+Foam::scalar Foam::sixDOFODE::odeStepFraction(const scalar odeTime) const
+{
+    // Get current global time. Note: assuming that the ODESolver for a given
+    // time step is integrating from t - deltaT to t.
+    const scalar globalTime = dict().time().value();
+    const scalar globalDeltaT = dict().time().deltaT().value();
+
+    // Calculate current ODE step time step size
+    const scalar odeDeltaT = globalDeltaT - (globalTime - odeTime);
+
+    // Calculate and return ODE step size fraction
+    return odeDeltaT/globalDeltaT;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::dimensionedScalar& Foam::sixDOFODE::mass() const
+{
+    return mass_;
+}
+
+
+Foam::dimensionedScalar& Foam::sixDOFODE::mass()
+{
+    return mass_;
+}
+
+
+const Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia() const
+{
+    return momentOfInertia_;
+}
+
+
+Foam::dimensionedDiagTensor& Foam::sixDOFODE::momentOfInertia()
+{
+    return momentOfInertia_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium() const
+{
+    return Xequilibrium_;
+}
+
+
+Foam::dimensionedVector& Foam::sixDOFODE::Xequilibrium()
+{
+    return Xequilibrium_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::force() const
+{
+    return force_;
+}
+
+
+const Foam::dimensionedVector& Foam::sixDOFODE::moment() const
+{
+    return moment_;
+}
+
+
+void Foam::sixDOFODE::setExternalForceAndMoment
+(
+    const dimensionedVector& externalForce,
+    const dimensionedVector& externalMoment
+)
+{
+    // Initialise the state before setting external forces and moments
+    initState();
+
+    // Set forces and moments
+    force_ = externalForce;
+    moment_ = externalMoment;
+
+    // Relax acceleration if the Aitkens relaxation is turned on
+    if (aitkensRelaxation_)
+    {
+        relaxAcceleration();
+    }
+}
+
+
+void Foam::sixDOFODE::initExternalForceAndMoment
+(
+    const dimensionedVector& externalForce,
+    const dimensionedVector& externalMoment
+)
+{
+    // Initialise force and moment only for the first time step
+    if (curTimeIndex_ == -1)
+    {
+        force_ = externalForce;
+        moment_ = externalMoment;
+    }
+}
+
+
+const Foam::PtrList<Foam::translationalConstraint>&
+Foam::sixDOFODE::translationalConstraints() const
+{
+    return translationalConstraints_;
+}
+
+
+const Foam::PtrList<Foam::rotationalConstraint>&
+Foam::sixDOFODE::rotationalConstraints() const
+{
+    return rotationalConstraints_;
+}
+
+
+const Foam::PtrList<Foam::translationalRestraint>&
+Foam::sixDOFODE::translationalRestraints() const
+{
+    return translationalRestraints_;
+}
+
+
+const Foam::PtrList<Foam::rotationalRestraint>&
+Foam::sixDOFODE::rotationalRestraints() const
+{
+    return rotationalRestraints_;
+}
+
+// ************************************************************************* //

src/OSspecific/MSWindows/printStack.C

@@ -126,7 +126,7 @@ void printSourceFileAndLine
         myAddress = nStream.str();
     }
 
-    if (filename[0] == '/')
+    if (filename[0] == '/' || filename[1] == ':')
     {
         string line = pOpen
         (
@@ -266,7 +266,10 @@ void getSymbolForRaw
     const word& address
 )
 {
-    if (filename.size() && filename[0] == '/')
+    if 
+    (
+        filename.size() && (filename[0] == '/' || filename[1] == ':')
+    )
     {
         string fcnt = pOpen
         (
@@ -316,7 +319,10 @@ void error::printStack(Ostream& os)
             string::size_type space = line.rfind(' ') + 1;
             fileName libPath = line.substr(space, line.size()-space);
 
-            if (libPath.size() && libPath[0] == '/')
+            if 
+            (
+                libPath.size() && (libPath[0] == '/' || libPath[1] == ':')
+            )
             {
                 string offsetString(line.substr(0, line.find('-')));
                 IStringStream offsetStr(offsetString);
@@ -359,10 +365,10 @@ void error::printStack(Ostream& os)
                 programFile = msg.substr(0, min(spacePos, bracketPos));
 
                 // not an absolute path
-                if (programFile[0] != '/')
+                if (programFile[0] != '/' && programFile[1] != ':')
                 {
                     string tmp = pOpen("which " + programFile);
-                    if (tmp[0] == '/' || tmp[0] == '~')
+                    if (tmp[0] == '/' || tmp[1] == ':' || tmp[0] == '~')
                     {
                         programFile = tmp;
                     }

src/OSspecific/POSIX/printStack.C

@@ -130,7 +130,7 @@ void printSourceFileAndLine
     }
 
 #ifndef darwin
-    if (filename[0] == '/')
+    if (filename[0] == '/' || filename[1] == ':')
 #else
     if (1)
 #endif
@@ -175,7 +175,10 @@ void getSymbolForRaw
     const word& address
 )
 {
-    if (filename.size() && filename[0] == '/')
+    if
+    (
+        filename.size() && (filename[0] == '/' || filename[1] == ':')
+    )
     {
         string fcnt = pOpen
         (
@@ -224,7 +227,10 @@ void error::printStack(Ostream& os)
             string::size_type space = line.rfind(' ') + 1;
             fileName libPath = line.substr(space, line.size()-space);
 
-            if (libPath.size() && libPath[0] == '/')
+            if
+            (
+                libPath.size() && (libPath[0] == '/' || libPath[1] == ':')
+            )
             {
                 string offsetString(line.substr(0, line.find('-')));
                 IStringStream offsetStr(offsetString);
@@ -268,10 +274,10 @@ void error::printStack(Ostream& os)
                 programFile = msg.substr(0, min(spacePos, bracketPos));
 
                 // not an absolute path
-                if (programFile[0] != '/')
+                if (programFile[0] != '/' && programFile[1] != ':')
                 {
                     string tmp = pOpen("which " + programFile);
-                    if (tmp[0] == '/' || tmp[0] == '~')
+                    if (tmp[0] == '/' || tmp[1] == ':' || tmp[0] == '~')
                     {
                         programFile = tmp;
                     }

src/dynamicMesh/dynamicFvMesh/Make/files

@@ -10,6 +10,7 @@ dynamicMotionSolverFvMesh/dynamicMotionSolverFvMesh.C
 subsetMotionSolverFvMesh/subsetMotionSolverFvMesh.C
 dynamicInkJetFvMesh/dynamicInkJetFvMesh.C
 dynamicRefineFvMesh/dynamicRefineFvMesh.C
+dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.C
 
 mixerGgiFvMesh/mixerGgiFvMesh.C
 turboFvMesh/turboFvMesh.C

src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C

@@ -31,6 +31,7 @@ License
 #include "syncTools.H"
 #include "pointFields.H"
 #include "directTopoChange.H"
+#include "cellSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,7 +60,16 @@ label dynamicRefineFvMesh::count
         {
             n++;
         }
+
+        // debug also serves to get-around Clang compiler trying to optimise
+        // out this forAll loop under O3 optimisation
+
+        if (debug)
+        {
+            Info<< "n=" << n << endl;
         }
+    }
+
     return n;
 }
 
@@ -860,11 +870,78 @@ void dynamicRefineFvMesh::extendMarkedCells(PackedBoolList& markedCell) const
 }
 
 
+void Foam::dynamicRefineFvMesh::checkEightAnchorPoints
+(
+    PackedBoolList& protectedCell,
+    label& nProtected
+) const
+{
+    const labelList& cellLevel = meshCutter_.cellLevel();
+    const labelList& pointLevel = meshCutter_.pointLevel();
+
+    labelList nAnchorPoints(nCells(), 0);
+
+    forAll(pointLevel, pointI)
+    {
+        const labelList& pCells = pointCells(pointI);
+
+        forAll(pCells, pCellI)
+        {
+            label cellI = pCells[pCellI];
+
+            if (pointLevel[pointI] <= cellLevel[cellI])
+            {
+                // Check if cell has already 8 anchor points -> protect cell
+                if (nAnchorPoints[cellI] == 8)
+                {
+                    if (protectedCell.set(cellI, true))
+                    {
+                        nProtected++;
+                    }
+                }
+
+                if (!protectedCell[cellI])
+                {
+                    nAnchorPoints[cellI]++;
+                }
+            }
+        }
+    }
+
+
+    forAll(protectedCell, cellI)
+    {
+        if (!protectedCell[cellI] && nAnchorPoints[cellI] != 8)
+        {
+            protectedCell.set(cellI, true);
+            nProtected++;
+        }
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
 :
     dynamicFvMesh(io),
+    singleMotionUpdate_
+    (
+        IOdictionary
+        (
+            IOobject
+            (
+                "dynamicMeshDict",
+                time().constant(),
+                *this,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        ).subDict(typeName + "Coeffs")
+        .lookupOrDefault<Switch>("singleMotionUpdate", true)
+    ),
+    curTimeIndex_(-1),
     meshCutter_(*this),
     dumpLevel_(false),
     nRefinementIterations_(0),
@@ -984,12 +1061,63 @@ dynamicRefineFvMesh::dynamicRefineFvMesh(const IOobject& io)
                 nProtected++;
             }
         }
+
+        // Also protect any cells that are less than hex
+        forAll(cells(), cellI)
+        {
+            const cell& cFaces = cells()[cellI];
+
+            if (cFaces.size() < 6)
+            {
+                if (protectedCell_.set(cellI, 1))
+                {
+                    nProtected++;
+                }
+            }
+            else
+            {
+                forAll(cFaces, cFaceI)
+                {
+                    if (faces()[cFaces[cFaceI]].size() < 4)
+                    {
+                        if (protectedCell_.set(cellI, 1))
+                        {
+                            nProtected++;
+                        }
+                        break;
+                    }
+                }
+            }
+        }
+
+        // Check cells for 8 corner points
+        checkEightAnchorPoints(protectedCell_, nProtected);
     }
 
     if (returnReduce(nProtected, sumOp<label>()) == 0)
     {
         protectedCell_.clear();
     }
+    else
+    {
+
+        cellSet protectedCells(*this, "protectedCells", nProtected);
+        forAll(protectedCell_, cellI)
+        {
+            if (protectedCell_[cellI])
+            {
+                protectedCells.insert(cellI);
+            }
+        }
+
+        Info<< "Detected " << returnReduce(nProtected, sumOp<label>())
+            << " cells that are protected from refinement."
+            << " Writing these to cellSet "
+            << protectedCells.name()
+            << "." << endl;
+
+        protectedCells.write();
+    }
 }
 
 
@@ -1003,6 +1131,20 @@ dynamicRefineFvMesh::~dynamicRefineFvMesh()
 
 bool dynamicRefineFvMesh::update()
 {
+    // Handling multiple calls in a single time step
+    if
+    (
+        singleMotionUpdate_
+     && curTimeIndex_ == this->time().timeIndex()
+    )
+    {
+        // This is not the first call to update, simply return false
+        return false;
+    }
+
+    // Update local time index
+    curTimeIndex_ = this->time().timeIndex();
+
     // Re-read dictionary. Choosen since usually -small so trivial amount
     // of time compared to actual refinement. Also very useful to be able
     // to modify on-the-fly.
@@ -1071,9 +1213,9 @@ bool dynamicRefineFvMesh::update()
                 << exit(FatalError);
         }
 
-        word field(refineDict.lookup("field"));
+        const word fieldName(refineDict.lookup("field"));
 
-        const volScalarField& vFld = lookupObject<volScalarField>(field);
+        const volScalarField& vFld = lookupObject<volScalarField>(fieldName);
 
         const scalar lowerRefineLevel =
             readScalar(refineDict.lookup("lowerRefineLevel"));

src/dynamicMesh/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.H

@@ -55,6 +55,20 @@ class dynamicRefineFvMesh
 :
     public dynamicFvMesh
 {
+
+    // Private data
+
+        // Helper variables to enable switching between a single and multiple
+        // mesh motion updates within a time step (if update() is called more
+        // than once in a single time step)
+
+            //- Switch for single motion update (true by default)
+            Switch singleMotionUpdate_;
+
+            //- Helper varaible: current time index
+            label curTimeIndex_;
+
+
 protected:
 
         //- Mesh cutting engine
@@ -148,6 +162,12 @@ protected:
             //- Extend markedCell with cell-face-cell.
             void extendMarkedCells(PackedBoolList& markedCell) const;
 
+            //- Check all cells have 8 anchor points
+            void checkEightAnchorPoints
+            (
+                PackedBoolList& protectedCell,
+                label& nProtected
+            ) const;
 
 private:
 

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicMeshDict

@@ -0,0 +1,77 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | foam-extend: Open Source CFD                    |
+|  \\    /   O peration     | Version:     4.0                                |
+|   \\  /    A nd           | Web:         http://www.foam-extend.org         |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//dynamicFvMeshLib   "libtopoChangerFvMesh.so";
+dynamicFvMesh      dynamicRefinePolyFvMesh;
+//staticFvMesh;
+
+mixerFvMeshCoeffs
+{
+    coordinateSystem
+    {
+        type            cylindrical;
+        origin          (0 0 0);
+        axis            (0 0 1);
+        direction       (1 0 0);
+    }
+
+    rpm             10;
+
+    slider
+    {
+        inside      insideSlider;
+        outside     outsideSlider;
+    }
+}
+
+// Refinement
+dynamicRefineFvMeshCoeffs
+{
+    // Refine every refineInterval timesteps
+    refineInterval 3;
+
+    // Maximum refinement level (starts from 0)
+    maxRefinement 2;
+
+    // Maximum cell limit (approximate)
+    maxCells 1000000;
+
+    // volScalarField to base refinement on
+    field gamma;
+
+    // Which cells to un/refine: based on point values (simple averaging).
+    // - refine pointCells of point value inbetween minLevel..maxLevel
+    // - unrefine pointCells that are within nBufferLayers of points marked
+    //   for refinement.
+    minLevel 0.01;
+    maxLevel 0.99;
+    nBufferLayers 1;
+
+    // Newly introduced patch points optionally get projected onto a surface
+    //projectSurfaces ("fixedWalls4.stl");
+    //projectPatches (fixedWalls);
+    // Maximum project distance
+    //projectDistance 1;
+
+    // Fluxes to adapt. For newly created faces or split faces the flux
+    // gets estimated from an interpolated volVectorField ('velocity')
+    // First is name of the flux to adapt, second is velocity that will
+    // be interpolated and inner-producted with the face area vector.
+    correctFluxes ((phi U));
+}
+
+
+// ************************************************************************* //

src/dynamicMesh/dynamicFvMesh/dynamicRefinePolyFvMesh/dynamicRefinePolyFvMesh.H

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::dynamicRefinePolyFvMesh
+
+Description
+    A fvMesh with built-in refinement of arbitrary polyhedral cells.
+
+    Determines which cells to refine/unrefine and does all in update().
+
+SourceFiles
+    dynamicRefinePolyFvMesh.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+Notes
+    Generalisation of dynamicRefineMesh for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef dynamicRefinePolyFvMesh_H
+#define dynamicRefinePolyFvMesh_H
+
+#include "dynamicFvMesh.H"
+#include "polyRef.H"
+#include "PackedBoolList.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class dynamicRefinePolyFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class dynamicRefinePolyFvMesh
+:
+    public dynamicFvMesh
+{
+
+    // Private data
+
+        // Helper variables to enable switching between a single and multiple
+        // mesh motion updates within a time step (if update() is called more
+        // than once in a single time step)
+
+            //- Switch for single motion update (true by default)
+            Switch singleMotionUpdate_;
+
+            //- Helper varaible: current time index
+            label curTimeIndex_;
+
+
+protected:
+
+        //- Mesh cutting engine
+        polyRef meshCutter_;
+
+        //- Dump cellLevel for postprocessing
+        Switch dumpLevel_;
+
+        //- Fluxes to map
+        List<Pair<word> > correctFluxes_;
+
+        //- Number of refinement/unrefinement steps done so far.
+        label nRefinementIterations_;
+
+
+    // Private Member Functions
+
+        //- Count set/unset elements in packedlist.
+        static label count(const PackedBoolList&, const unsigned int);
+
+        //- Read the projection parameters from dictionary
+        void readDict();
+
+
+        //- Refine cells. Update mesh and fields.
+        autoPtr<mapPolyMesh> refine(const labelList&);
+
+        //- Unrefine cells. Gets passed in centre points of cells to combine.
+        autoPtr<mapPolyMesh> unrefine(const labelList&);
+
+
+        // Selection of cells to un/refine
+
+            //- Get per cell max of connected point
+            scalarField maxPointField(const scalarField&) const;
+
+            //- Get point min of connected cell
+            scalarField minCellField(const volScalarField&) const;
+
+            //- Simple (non-parallel) interpolation by averaging
+            scalarField cellToPoint(const scalarField& vFld) const;
+
+            //- Calculate error (= -1 by default or distance from inbetween
+            //  levels
+            scalarField error
+            (
+                const scalarField& fld,
+                const scalar minLevel,
+                const scalar maxLevel
+            ) const;
+
+            //- Select candidate cells for refinement
+            virtual void selectRefineCandidates
+            (
+                const scalar lowerRefineLevel,
+                const scalar upperRefineLevel,
+                const scalarField& vFld,
+                PackedBoolList& candidateCell
+            ) const;
+
+            //- Subset candidate cells for refinement
+            virtual labelList selectRefineCells
+            (
+                const label maxCells,
+                const label maxRefinement,
+                const PackedBoolList& candidateCell
+            ) const;
+
+            //- Select points that can be unrefined
+            virtual labelList selectUnrefinePoints
+            (
+                const scalar unrefineLevel,
+                const PackedBoolList& markedCell,
+                const scalarField& pFld
+            ) const;
+
+            //- Extend markedCell with cell-face-cell
+            void extendMarkedCells(PackedBoolList& markedCell) const;
+
+
+private:
+
+        //- Disallow default bitwise copy construct
+        dynamicRefinePolyFvMesh(const dynamicRefinePolyFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const dynamicRefinePolyFvMesh&);
+
+public:
+
+    //- Runtime type information
+    TypeName("dynamicRefinePolyFvMesh");
+
+
+    // Constructors
+
+        //- Construct from IOobject
+        explicit dynamicRefinePolyFvMesh(const IOobject& io);
+
+
+    // Destructor
+
+        virtual ~dynamicRefinePolyFvMesh();
+
+
+    // Member Functions
+
+        //- Direct access to the refinement engine
+        const polyRef& meshCutter() const
+        {
+            return meshCutter_;
+        }
+
+        //- Update the mesh for both mesh motion and topology change
+        virtual bool update();
+
+
+    // Writing
+
+        //- Write using given format, version and compression
+        virtual bool writeObject
+        (
+            IOstream::streamFormat fmt,
+            IOstream::versionNumber ver,
+            IOstream::compressionType cmp
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/dynamicMesh/Make/files

@@ -52,6 +52,7 @@ motionSolver/motionSolver.C
 refinementData/refinementData.C
 refinementData/refinementDistanceData.C
 refinementData/refinementHistory.C
+refinementData/polyRefinementHistory.C
 
 directTopoChange/directTopoChange/directTopoChange.C
 
@@ -60,6 +61,7 @@ $(directActions)/addPatchCellLayer.C
 $(directActions)/edgeCollapser.C
 $(directActions)/faceCollapser.C
 $(directActions)/hexRef8.C
+$(directActions)/polyRef.C
 $(directActions)/removeCells.C
 $(directActions)/removeFaces.C
 $(directActions)/removePoints.C

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/hexRef8.C

@@ -353,6 +353,16 @@ void Foam::hexRef8::modFace
 // Bit complex way to determine the unrefined edge length.
 Foam::scalar Foam::hexRef8::getLevel0EdgeLength() const
 {
+    if (cellLevel_.size() != mesh_.nCells())
+    {
+        FatalErrorIn("scalar hexRef8::getLevel0EdgeLength() const")
+            << "Number of cells in mesh:" << mesh_.nCells()
+            << " does not equal size of cellLevel:" << cellLevel_.size()
+            << endl
+            << "This might be because of a restart with inconsistent cellLevel."
+            << abort(FatalError);
+    }
+
     // Determine minimum edge length per refinement level
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 

src/dynamicMesh/dynamicMesh/directTopoChange/directTopoChange/directActions/polyRef.H

@@ -0,0 +1,558 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::polyRef
+
+Description
+    Refinement of (split) polyhedra using directTopoChange.
+
+    Each polyhedron is split by the following procedure:
+    1. Adding points at the edge centre, face centre and cell centre,
+    2. Adding cells n cells where n is the number of points of the cell,
+    3. Splitting each face into multiple faces going from:
+       existing corner point -> new edge centre point -> new face centre
+       point -> other new edge centre point (sharing the same corner point)
+    4. Adding internal faces going from:
+       new edge centre point -> new face centre point -> new cell centre
+       point -> other new face centre point (sharing the same edge)
+
+SourceFiles
+    polyRef.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved.
+
+Notes
+    Generalisation of hexRef8 for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyRef_H
+#define polyRef_H
+
+#include "labelIOList.H"
+#include "face.H"
+#include "HashSet.H"
+#include "DynamicList.H"
+#include "primitivePatch.H"
+#include "removeFaces.H"
+#include "polyRefinementHistory.H"
+#include "PackedList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class polyMesh;
+class polyPatch;
+class directTopoChange;
+class mapPolyMesh;
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                           Class polyRef Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyRef
+{
+    // Private data
+
+        //- Reference to underlying mesh
+        const polyMesh& mesh_;
+
+        //- Per cell the refinement level
+        labelIOList cellLevel_;
+
+        //- Per point the refinement level
+        labelIOList pointLevel_;
+
+        //- Typical edge length between unrefined points
+        const scalar level0Edge_;
+
+        //- Refinement history
+        polyRefinementHistory history_;
+
+        //- Face remover engine
+        removeFaces faceRemover_;
+
+        //- Level of saved points
+        Map<label> savedPointLevel_;
+
+        //- Level of saved cells
+        Map<label> savedCellLevel_;
+
+
+    // Private Member Functions
+
+        //- Reorder according to map
+        static void reorder
+        (
+            const labelList& map,
+            const label len,
+            const label null,
+            labelList& elems
+        );
+
+        //- Get patch and zone info
+        void getFaceInfo
+        (
+            const label faceI,
+            label& patchID,
+            label& zoneID,
+            label& zoneFlip
+        ) const;
+
+        //- Adds a face on top of existing faceI. Reverses if nessecary
+        label addFace
+        (
+            directTopoChange& meshMod,
+            const label faceI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        //- Adds internal face from point. No checks on reversal
+        label addInternalFace
+        (
+            directTopoChange& meshMod,
+            const label meshFaceI,
+            const label meshPointI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        //- Modifies existing faceI for either new owner/neighbour or new face
+        //  points. Reverses if nessecary
+        void modFace
+        (
+            directTopoChange& meshMod,
+            const label faceI,
+            const face& newFace,
+            const label own,
+            const label nei
+        ) const;
+
+        scalar getLevel0EdgeLength() const;
+
+        //- Get cell added to point of cellI (if any)
+        label getAnchorCell
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const label cellI,
+            const label faceI,
+            const label pointI
+        ) const;
+
+        //- Get new owner and neighbour (in unspecified order) of pointI
+        //  on faceI
+        void getFaceNeighbours
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const label faceI,
+            const label pointI,
+
+            label& own,
+            label& nei
+        ) const;
+
+
+        //- Get index of point with minimum pointlevel
+        label findMinLevel(const labelList& f) const;
+
+        //- Get index of point with maximum pointlevel
+        label findMaxLevel(const labelList& f) const;
+
+        //- Count number of vertices <= anchorLevel
+        label countAnchors(const labelList&, const label) const;
+
+        //- Find index of point with wantedLevel, starting from fp
+        label findLevel
+        (
+            const face& f,
+            const label startFp,
+            const bool searchForward,
+            const label wantedLevel
+        ) const;
+
+        ////- Print levels of list of points.
+        //void printLevels(Ostream&, const labelList&) const;
+
+        //- debug: check orientation of added internal face
+        static void checkInternalOrientation
+        (
+            directTopoChange& meshMod,
+            const label cellI,
+            const label faceI,
+            const point& ownPt,
+            const point& neiPt,
+            const face& newFace
+        );
+
+        //- debug: check orientation of new boundary face
+        static void checkBoundaryOrientation
+        (
+            directTopoChange& meshMod,
+            const label cellI,
+            const label faceI,
+            const point& ownPt,
+            const point& boundaryPt,
+            const face& newFace
+        );
+
+        //- If p0 and p1 are existing vertices check if edge is split and insert
+        //  splitPoint
+        void insertEdgeSplit
+        (
+            const labelList& edgeMidPoint,
+            const label p0,
+            const label p1,
+            dynamicLabelList& verts
+        ) const;
+
+        //- Store in maps correspondence from midpoint to anchors and faces
+        label storeMidPointInfo
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            const label faceI,
+            const bool faceOrder,
+            const label midPointI,
+            const label anchorPointI,
+            const label faceMidPointI,
+
+            Map<edge>& midPointToAnchors,
+            Map<edge>& midPointToFaceMids,
+            directTopoChange& meshMod
+        ) const;
+
+        //- Create all internal faces from an unsplit face
+        void createInternalFromSplitFace
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& faceMidPoint,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            const label faceI,
+
+            Map<edge>& midPointToAnchors,
+            Map<edge>& midPointToFaceMids,
+            directTopoChange& meshMod,
+            label& nFacesAdded
+        ) const;
+
+        //- Create all internal faces to split cellI into n cells where n is the
+        //  number of cell points
+        void createInternalFaces
+        (
+            const labelListList& cellAnchorPoints,
+            const labelListList& cellAddedCells,
+            const labelList& cellMidPoint,
+            const labelList& faceMidPoint,
+            const labelList& faceAnchorLevel,
+            const labelList& edgeMidPoint,
+            const label cellI,
+            directTopoChange& meshMod
+        ) const;
+
+        //- Store vertices from startFp upto face split point.
+        //  Used when splitting face into n faces where n is the number of
+        //  points in a face (or number of edges)
+        void walkFaceToMid
+        (
+            const labelList& edgeMidPoint,
+            const label cLevel,
+            const label faceI,
+            const label startFp,
+            dynamicLabelList& faceVerts
+        ) const;
+
+        //- Same as walkFaceToMid but now walk back
+        void walkFaceFromMid
+        (
+            const labelList& edgeMidPoint,
+            const label cLevel,
+            const label faceI,
+            const label startFp,
+            dynamicLabelList& faceVerts
+        ) const;
+
+        //- Updates refineCell so consistent 2:1 refinement. Returns local
+        //  number of cells changed
+        label faceConsistentRefinement
+        (
+            const bool maxSet,
+            PackedList<1>& refineCell
+        ) const;
+
+        //- Check wanted refinement for 2:1 consistency
+        void checkWantedRefinementLevels(const labelList&) const;
+
+
+        // Copy control
+
+            //- Disallow default bitwise copy construct
+            polyRef(const polyRef&);
+
+            //- Disallow default bitwise assignment
+            void operator=(const polyRef&);
+
+
+public:
+
+    //- Runtime type information
+    ClassName("polyRef");
+
+
+    // Constructors
+
+        //- Construct from mesh, read_if_present refinement data
+        //  (from write below)
+        polyRef(const polyMesh& mesh);
+
+        //- Construct from mesh and un/refinement data
+        polyRef
+        (
+            const polyMesh& mesh,
+            const labelList& cellLevel,
+            const labelList& pointLevel,
+            const polyRefinementHistory& history
+        );
+
+        //- Construct from mesh and refinement data.
+        polyRef
+        (
+            const polyMesh& mesh,
+            const labelList& cellLevel,
+            const labelList& pointLevel
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            const labelIOList& cellLevel() const
+            {
+                return cellLevel_;
+            }
+
+            const labelIOList& pointLevel() const
+            {
+                return pointLevel_;
+            }
+
+            const polyRefinementHistory& history() const
+            {
+                return history_;
+            }
+
+            //- Typical edge length between unrefined points
+            scalar level0EdgeLength() const
+            {
+                return level0Edge_;
+            }
+
+        // Refinement
+
+            // Get cell level such that the face has the most points that are
+            // <= level (at least three points)
+            label getAnchorLevel(const label faceI) const;
+
+            //- Helper: get points of a cell without using cellPoints addressing
+            labelList cellPoints(const label cellI) const;
+
+            //- Given valid mesh and current cell level and proposed
+            //  cells to refine calculate any clashes (due to 2:1) and return
+            //  ok list of cells to refine.
+            //  Either adds cells to refine to set (maxSet = true) or
+            //  removes cells to refine (maxSet = false)
+            labelList consistentRefinement
+            (
+                const labelList& cellsToRefine,
+                const bool maxSet
+            ) const;
+
+            //- Like consistentRefinement but slower:
+            //  - specify number of cells between consecutive refinement levels
+            //    (consistentRefinement equivalent to 1)
+            //  - specify max level difference between point-connected cells.
+            //    (-1 to disable). Note that with normal 2:1 limitation
+            //    (maxFaceDiff=1) there can be 8:1 size difference across point
+            //    connected cells so maxPointDiff allows you to make that less.
+            //  cellsToRefine : cells we're thinking about refining. It will
+            //                  extend this set. All refinement levels will be
+            //                  at least maxFaceDiff layers thick.
+            //  facesToCheck  : additional faces where to implement the
+            //                  maxFaceDiff thickness (usually only boundary
+            //                  faces)
+            labelList consistentSlowRefinement
+            (
+                const label maxFaceDiff,
+                const labelList& cellsToRefine,
+                const labelList& facesToCheck,
+                const label maxPointDiff,
+                const labelList& pointsToCheck
+            ) const;
+
+            //- Like consistentSlowRefinement but uses different meshWave
+            //  (proper distance instead of toplogical count). No point checks
+            //  yet
+            labelList consistentSlowRefinement2
+            (
+                const label maxFaceDiff,
+                const labelList& cellsToRefine,
+                const labelList& facesToCheck
+            ) const;
+
+            //- Insert refinement. All selected cells will be split into n cells
+            //  where n is the number of points per cell.
+            //  Returns per element in cells the n cells they were split into.
+            //  Guarantees that the 0th element is the original cell label.
+            //  Mapping:
+            //  -split cells: n new ones get added from original
+            //  -split faces: original gets modified; n - 1 new ones get added
+            //                from original
+            //  -added internal faces: added from original cell face (if
+            //   that was internal) or created out-of-nothing (so will not
+            //   get mapped!). Note: could make this inflate from point but
+            //   that will allocate interpolation.
+            //  -points added to split edge: added from edge start()
+            //  -midpoints added: added from cellPoints[0].
+            labelListList setRefinement
+            (
+                const labelList& cells,
+                directTopoChange&
+            );
+
+            //- Update local numbering for changed mesh
+            void updateMesh(const mapPolyMesh&);
+
+
+            // Restoring : is where other processes delete and reinsert data.
+            //             These callbacks allow this to restore the cellLevel
+            //             and pointLevel for reintroduced points.
+            //             Is not related to undoing my refinement
+
+                //- Signal points/face/cells for which to store data
+                void storeData
+                (
+                    const labelList& pointsToStore,
+                    const labelList& facesToStore,
+                    const labelList& cellsToStore
+                );
+
+                //- Update local numbering + undo
+                //  Data to restore given as new pointlabel + stored pointlabel
+                //  (i.e. what was in pointsToStore)
+                void updateMesh
+                (
+                    const mapPolyMesh&,
+                    const Map<label>& pointsToRestore,
+                    const Map<label>& facesToRestore,
+                    const Map<label>& cellsToRestore
+                );
+
+
+            //- Update local numbering for subsetted mesh.
+            //  Gets new-to-old maps. Not compatible with unrefinement.
+            void subset
+            (
+                const labelList& pointMap,
+                const labelList& faceMap,
+                const labelList& cellMap
+            );
+
+            //- Update local numbering for mesh redistribution
+            void distribute(const mapDistributePolyMesh&);
+
+            //- Debug: Check coupled mesh for correctness
+            void checkMesh() const;
+
+            //- Debug: Check 2:1 consistency across faces.
+            //  maxPointDiff==-1 : only check 2:1 across faces
+            //  maxPointDiff!=-1 : check point-connected cells.
+            void checkRefinementLevels
+            (
+                const label maxPointDiff,
+                const labelList& pointsToCheck
+            ) const;
+
+
+        // Unrefinement (undoing refinement, not arbitrary coarsening)
+
+            //- Return the points at the centre of top-level split cells
+            //  that can be unsplit.
+            labelList getSplitPoints() const;
+
+            //- Given proposed splitPoints to unrefine according to calculate
+            //  any clashes (due to 2:1) and return ok list of points to
+            //  unrefine. Either adds points to refine to set (maxSet = true)
+            //  or removes points to refine (maxSet = false)
+            labelList consistentUnrefinement
+            (
+                const labelList& pointsToUnrefine,
+                const bool maxSet
+            ) const;
+
+            //- Remove some refinement. Needs to be supplied output of
+            //  consistentUnrefinement. Only call if undoable set.
+            //  All n pointCells of a split point will be combined into
+            //  the lowest numbered cell of those n.
+            void setUnrefinement
+            (
+                const labelList& splitPointLabels,
+                directTopoChange&
+            );
+
+
+        // Write
+
+            // Set instance for mesh files
+            void setInstance(const fileName& inst);
+
+            //- Force writing refinement+history to polyMesh directory.
+            bool write() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/dynamicMesh/refinementData/polyRefinementHistory.H

@@ -0,0 +1,378 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::polyRefinementHistory
+
+Description
+    All refinement history. Used in unrefinement.
+
+    - visibleCells: valid for the current mesh and contains per cell -1
+      (cell unrefined) or an index into splitCells_.
+    - splitCells: for every split contains the parent (also index into
+      splitCells) and optionally a subsplit as n indices into splitCells, where
+      n is the number of points of a cell. Note that the numbers in splitCells
+      are not cell labels, they are purely indices into splitCells.
+
+    E.g. 2 cells, cell 1 gets refined so end up with 9 cells:
+    @verbatim
+        // splitCells
+        9
+        (
+        -1 (1 2 3 4 5 6 7 8)
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        0 0()
+        )
+
+        // visibleCells
+        9(-1 1 2 3 4 5 6 7 8)
+    @endverbatim
+
+
+    So cell0 (visibleCells=-1) is unrefined.
+    Cells 1-8 have all valid splitCells entries which are:
+      - parent:0
+      - subsplits:0()
+
+    The parent 0 refers back to the splitcell entries.
+
+
+SourceFiles
+    polyRefinementHistory.C
+
+Author
+    Vuko Vukcevic, Wikki Ltd.  All rights reserved
+
+Notes
+    Generalisation of refinementHistory for polyhedral cells
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef polyRefinementHistory_H
+#define polyRefinementHistory_H
+
+#include "DynamicList.H"
+#include "dynamicLabelList.H"
+#include "labelList.H"
+#include "FixedList.H"
+#include "SLList.H"
+#include "autoPtr.H"
+#include "regIOobject.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class mapPolyMesh;
+class mapDistributePolyMesh;
+
+/*---------------------------------------------------------------------------*\
+                     Class polyRefinementHistory Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyRefinementHistory
+:
+    public regIOobject
+{
+public:
+
+    class splitPolyCell
+    {
+    public:
+
+        // Index to original splitCell this cell was refined off from
+        // -1: top level cell
+        // -2: free splitCell (so should also be in freeSplitCells_)
+        label parent_;
+
+        //- cells this cell was refined into
+        autoPtr<labelList> addedCellsPtr_;
+
+        //- Construct null (parent = -1)
+        splitPolyCell();
+
+        //- Construct from parent
+        splitPolyCell(const label parent);
+
+        //- Construct from Istream
+        splitPolyCell(Istream& is);
+
+        //- Construct as deep copy
+        splitPolyCell(const splitPolyCell&);
+
+        //- Copy operator since autoPtr otherwise 'steals' storage.
+        void operator=(const splitPolyCell& s)
+        {
+            // Check for assignment to self
+            if (this == &s)
+            {
+                FatalErrorIn
+                (
+                    "splitPolyCell::operator=(const Foam::splitPolyCell&)"
+                )   << "Attempted assignment to self"
+                    << abort(FatalError);
+            }
+
+            parent_ = s.parent_;
+
+            addedCellsPtr_.reset
+            (
+                s.addedCellsPtr_.valid()
+              ? new labelList(s.addedCellsPtr_())
+              : NULL
+            );
+        }
+
+        bool operator==(const splitPolyCell& s) const
+        {
+            if (addedCellsPtr_.valid() != s.addedCellsPtr_.valid())
+            {
+                return false;
+            }
+            else if (parent_ != s.parent_)
+            {
+                return false;
+            }
+            else if (addedCellsPtr_.valid())
+            {
+                return addedCellsPtr_() == s.addedCellsPtr_();
+            }
+            else
+            {
+                return true;
+            }
+        }
+
+        bool operator!=(const splitPolyCell& s) const
+        {
+            return !operator==(s);
+        }
+
+        friend Istream& operator>>(Istream&, splitPolyCell&);
+        friend Ostream& operator<<(Ostream&, const splitPolyCell&);
+    };
+
+
+private:
+
+    TypeName("polyRefinementHistory");
+
+    // Private data
+
+        //- Storage for splitCells
+        DynamicList<splitPolyCell> splitCells_;
+
+        //- Unused indices in splitCells
+        dynamicLabelList freeSplitCells_;
+
+        //- Currently visible cells. Indices into splitCells.
+        labelList visibleCells_;
+
+
+    // Private member functions
+
+        //- Debug write
+        static void writeEntry
+        (
+            const List<splitPolyCell>&,
+            const splitPolyCell&
+        );
+
+        //- Debug write
+        static void writeDebug
+        (
+            const labelList&,
+            const List<splitPolyCell>&
+        );
+
+        //- Check consistency of structure, i.e. indices into splitCells_.
+        void checkIndices() const;
+
+        //- Allocate a splitCell. Return index in splitCells_.
+        label allocateSplitCell
+        (
+            const label parent,
+            const label i,
+            const label nCells
+        );
+
+        //- Free a splitCell.
+        void freeSplitCell(const label index);
+
+        //- Mark entry in splitCells. Recursively mark its parent and subs.
+        void markSplit
+        (
+            const label,
+            labelList& oldToNew,
+            DynamicList<splitPolyCell>&
+        ) const;
+
+        void countProc
+        (
+            const label index,
+            const label newProcNo,
+            labelList& splitCellProc,
+            labelList& splitCellNum
+        ) const;
+
+public:
+
+    // Constructors
+
+        //- Construct (read) given an IOobject
+        polyRefinementHistory(const IOobject&);
+
+        //- Construct (read) or construct null
+        polyRefinementHistory
+        (
+            const IOobject&,
+            const List<splitPolyCell>& splitCells,
+            const labelList& visibleCells
+        );
+
+        //- Construct (read) or construct from initial number of cells
+        //  (all visible)
+        polyRefinementHistory(const IOobject&, const label nCells);
+
+        //- Construct as copy
+        polyRefinementHistory(const IOobject&, const polyRefinementHistory&);
+
+        //- Construct from Istream
+        polyRefinementHistory(const IOobject&, Istream&);
+
+
+    // Member Functions
+
+        //- Per cell in the current mesh (i.e. visible) either -1 (unrefined)
+        //  or an index into splitCells.
+        const labelList& visibleCells() const
+        {
+            return visibleCells_;
+        }
+
+        //- Storage for splitPolyCells.
+        const DynamicList<splitPolyCell>& splitCells() const
+        {
+            return splitCells_;
+        }
+
+        //- Cache of unused indices in splitCells
+        const dynamicLabelList& freeSplitCells() const
+        {
+            return freeSplitCells_;
+        }
+
+        //- Is there unrefinement history. Note that this will fall over if
+        //  there are 0 cells in the mesh. But this gives problems with
+        //  lots of other programs anyway.
+        bool active() const
+        {
+            return visibleCells_.size() > 0;
+        }
+
+        //- Get parent of cell
+        label parentIndex(const label cellI) const
+        {
+            label index = visibleCells_[cellI];
+
+            if (index < 0)
+            {
+                FatalErrorIn("polyRefinementHistory::parentIndex(const label)")
+                    << "Cell " << cellI << " is not visible"
+                    << abort(FatalError);
+            }
+            return splitCells_[index].parent_;
+        }
+
+        //- Store splitting of cell into n (number of points for a cell)
+        void storeSplit
+        (
+            const label cellI,
+            const labelList& addedCells
+        );
+
+        //- Store combining n cells into master
+        void combineCells
+        (
+            const label masterCellI,
+            const labelList& combinedCells
+        );
+
+        //- Update numbering for mesh changes
+        void updateMesh(const mapPolyMesh&);
+
+        //- Update numbering for subsetting
+        void subset
+        (
+            const labelList& pointMap,
+            const labelList& faceMap,
+            const labelList& cellMap
+        );
+
+        //- Update local numbering for mesh redistribution.
+        //  Can only distribute clusters sent across in one go; cannot
+        //  handle parts recombined in multiple passes.
+        void distribute(const mapDistributePolyMesh&);
+
+        //- Compact splitCells_. Removes all freeSplitCells_ elements.
+        void compact();
+
+        //- Extend/shrink storage. additional visibleCells_ elements get
+        //  set to -1.
+        void resize(const label nCells);
+
+        //- Debug write
+        void writeDebug() const;
+
+
+        //- ReadData function required for regIOobject read operation
+        virtual bool readData(Istream&);
+
+        //- WriteData function required for regIOobject write operation
+        virtual bool writeData(Ostream&) const;
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, polyRefinementHistory&);
+        friend Ostream& operator<<(Ostream&, const polyRefinementHistory&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/meshMotion/solidBodyMotion/solidBodyMotionFunction/solidBodyMotionFunction.H

@@ -64,7 +64,7 @@ protected:
 
     // Protected data
 
-        //- Solid bodi motion coefficients dictionary
+        //- Solid body motion coefficients dictionary
         dictionary SBMFCoeffs_;
 
         //- Time

src/finiteArea/faMatrices/faMatrix/faMatrix.C

@@ -419,12 +419,13 @@ void faMatrix<Type>::setValues
                 }
                 else
                 {
-                    label patchi = mesh.boundary().whichPatch(edgei);
+                    const label& curEdgeIndex = mesh.edgeIndex()[edgei];
+                    label patchi = mesh.boundary().whichPatch(curEdgeIndex);
 
                     if (internalCoeffs_[patchi].size())
                     {
                         label patchEdgei =
-                            mesh.boundary()[patchi].whichEdge(edgei);
+                            mesh.boundary()[patchi].whichEdge(curEdgeIndex);
 
                         internalCoeffs_[patchi][patchEdgei] =
                             pTraits<Type>::zero;

src/finiteArea/faMesh/faMesh.C

@@ -71,6 +71,7 @@ void Foam::faMesh::setPrimitiveMeshData()
 
     // Set faMesh edges
     edges_.setSize(bp.nEdges());
+    edgeIndex_.setSize(bp.nEdges());
 
     label edgeI = -1;
 
@@ -80,6 +81,7 @@ void Foam::faMesh::setPrimitiveMeshData()
     for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
     {
         edges_[++edgeI] = bp.edges()[curEdge];
+        edgeIndex_[curEdge] = edgeI;
     }
 
     forAll (boundary(), patchI)
@@ -89,6 +91,7 @@ void Foam::faMesh::setPrimitiveMeshData()
         forAll (curP, eI)
         {
             edges_[++edgeI] = bp.edges()[curP[eI]];
+            edgeIndex_[curP[eI]] = edgeI;
         }
     }
 
@@ -259,7 +262,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(pMesh.time().timePath()/"S0"))
+    if (isFile(pMesh.time().timePath()/pMesh.dbDir()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (
@@ -808,7 +811,7 @@ Foam::faMesh::faMesh
     // Calculate the geometry for the patches (transformation tensors etc.)
     boundary_.calcGeometry();
 
-    if (isFile(mesh().time().timePath()/"S0"))
+    if (isFile(mesh().time().timePath()/mesh().dbDir()/"S0"))
     {
         S0Ptr_ = new DimensionedField<scalar, areaMesh>
         (

src/finiteArea/faMesh/faMesh.H

@@ -82,6 +82,9 @@ class faMesh
         //- Boundary mesh
         faBoundaryMesh boundary_;
 
+        //- Edge index
+        labelList edgeIndex_;
+
 
         // Primitive mesh data
 
@@ -397,6 +400,11 @@ public:
                 return faceLabels_;
             }
 
+            //- Return faMesh edge indices
+            const labelList& edgeIndex() const
+            {
+                return edgeIndex_;
+            }
 
             //- Return parallel info
             const faGlobalMeshData& globalData() const;