A module for rotation parametrization#

Objectives

This document describes the implementation of generic vectorial parametrization of rotation matrices, essentially following [Bauchau and Trainelli, 2003]. This module is used in the Nonlinear beam model in finite rotations tour.

Implementation aspects#

The module provides the Skew function, mapping a 3D vector to the corresponding skew-symmetric matrix.

An abstract class handles the generic implementation of rotation parametrization based on the corresponding parametrization of the rotation angle \(p(\varphi)\). Considering a generic rotation vector which we denote by \(\boldsymbol{\theta}\), [Bauchau and Trainelli, 2003] works with \(\varphi = \|\boldsymbol{\theta}\|\) and the unit-norm vector \(\boldsymbol{u} = \boldsymbol{\theta}/\varphi\). Note that the involved expressions are usually ill-defined when \(\varphi \to 0\). For this reason, the numerical implementation makes use of a regularized expression for the norm:

(34)#\[\begin{equation} \varphi = \sqrt{\boldsymbol{\theta}^2 + \varepsilon} \end{equation}\]

with \({\varepsilon}\) being a very small number in practice.

The rotation parameter vector \(\boldsymbol{p}\) from [Bauchau and Trainelli, 2003] is given by the rotation_parameter attribute.

The class then implements the following functions:

(35)#\[\begin{align} h_1(\varphi) &= \dfrac{\nu(\varphi)^2}{\epsilon(\varphi)}\\ h_2(\varphi) &= \dfrac{\nu(\varphi)^2}{2}\\ h_3(\varphi) &= \dfrac{\mu(\varphi)-h_1(\varphi)}{p(\varphi)^2}\\ \end{align}\]

where \(\nu(\varphi),\epsilon(\varphi)\) and \(\mu(\varphi)\) are defined in [Bauchau and Trainelli, 2003].

It then provides the expression for the corresponding rotation matrix \(\boldsymbol{R}\):

(36)#\[\begin{equation} \boldsymbol{R} = \boldsymbol{I} + h_1(\varphi)\boldsymbol{P} + h_2(\varphi)\boldsymbol{P}^2 \end{equation}\]

where \(\boldsymbol{P} = \operatorname{skew}(\boldsymbol{p})\), as well as the associated rotation curvature matrix \(\boldsymbol{H}\) involved in the computation of the rotation rate:

(37)#\[\begin{equation} \boldsymbol{H} = \mu(\varphi)\boldsymbol{I} + h_2(\varphi)\boldsymbol{P} + h_3(\varphi)\boldsymbol{P}^2 \end{equation}\]

Available particular cases#

ExponentialMap parametrization#

This parametrization corresponds to the simple choice:

(38)#\[\begin{equation} p(\varphi)=\varphi \end{equation}\]

The corresponding expression for the rotation matrix is the famous Euler-Rodrigues formula.

EulerRodrigues parametrization#

This parametrization corresponds to the simple choice:

(39)#\[\begin{equation} p(\varphi)=2 sin(\varphi/2) \end{equation}\]

SineFamily parametrization#

This generic family for any integer \(m\) corresponds to:

(40)#\[\begin{equation} p(\varphi) = m \sin\left(\frac{\varphi}{m}\right) \end{equation}\]

TangentFamily parametrization#

This generic family for any integer \(m\) corresponds to:

(41)#\[\begin{equation} p(\varphi) = m \tan\left(\frac{\varphi}{m}\right) \end{equation}\]

Code#

The corresponding module is available here rotation_parametrization.py.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Nov 28 14:28:31 2020

@author: bleyerj
"""
from ufl import sqrt, dot, sin, cos, tan, as_matrix, Identity
from dolfinx.fem import Constant

DOLFIN_EPS = 1e-12


def Skew(n):
    """Antisymmetric tensor associated with a vector n"""
    return as_matrix([[0, -n[2], n[1]], [n[2], 0, -n[0]], [-n[1], n[0], 0]])


class RotationParametrization(object):
    """
    A generic class for handling a vectorial parametrization of rotation
    tensors.

        Bauchau, O. A., & Trainelli, L. (2003). The vectorial parameterization of
        rotation. Nonlinear dynamics, 32(1), 71-92. DOI: 10.1023/A:1024265401576
    """

    def __init__(self, p, dp):
        """
        Parameters
        ----------
        p : parametric function (must be such that p(x)/x -> 1 when x->0)

        dp : the derivative of p
        """
        self.p = p
        self.dp = dp

    def ntheta(self, theta):
        """Rotation vector norm"""
        return sqrt(dot(theta, theta) + DOLFIN_EPS)

    def h1(self, theta):
        x = self.ntheta(theta)
        return sin(x) / self.p(x)

    def h2(self, theta):
        x = self.ntheta(theta)
        return 2 * (sin(x / 2) / self.p(x)) ** 2

    def mu(self, theta):
        x = self.ntheta(theta)
        return 1 / self.dp(x)

    def h3(self, theta):
        x = self.ntheta(theta)
        return (self.mu(x) - self.h1(x)) / self.p(x) ** 2

    def rotation_parameter(self, theta):
        """Reparametrized rotation vector"""
        x = self.ntheta(theta)
        return self.p(x) * theta / x

    def rotation_matrix(self, theta):
        """Rotation matrix"""
        P = Skew(self.rotation_parameter(theta))
        return Identity(3) + self.h1(theta) * P + self.h2(theta) * P * P

    def curvature_matrix(self, theta):
        """Curvature matrix involved in the computation of the rotation rate"""
        P = Skew(self.rotation_parameter(theta))
        return (
            self.mu(theta) * Identity(3) + self.h2(theta) * P + self.h3(theta) * P * P
        )


class ExponentialMap(RotationParametrization):
    """Exponential map parametrization (p(x)=x)"""

    def __init__(self):
        RotationParametrization.__init__(self, lambda x: x, lambda x: 1)


class SineFamily(RotationParametrization):
    """Sine family parametrization (p(x)=m*sin(x/m))"""

    def __init__(self, m):
        RotationParametrization.__init__(
            self, lambda x: m * sin(x / m), lambda x: cos(x / m)
        )


class EulerRodrigues(SineFamily):
    """Euler-Rodrigues parametrization (p(x)=2*sin(x/2))"""

    def __init__(self):
        SineFamily.__init__(self, 2)


class TangentFamily(RotationParametrization):
    """Tangent family parametrization (p(x)=m*tan(x/m))"""

    def __init__(self, m):
        RotationParametrization.__init__(
            self, lambda x: m * tan(x / m), lambda x: 1 + tan(x / m) ** 2
        )


class RodriguesParametrization(RotationParametrization):
    def __init__(self):
        pass

    def rotation_parameter(self, theta):
        """Reparametrized rotation vector"""
        return theta

    def rotation_matrix(self, theta):
        P = Skew(self.rotation_parameter(theta))
        h = lambda x: 4 / (4 + x**2)
        return Identity(3) + h(theta) * P + h(theta) / 2 * P * P

    def cruvature_matrix(self, theta):
        """Curvature matrix involved in the computation of the rotation rate"""
        P = Skew(self.rotation_parameter(theta))
        h = lambda x: 4 / (4 + x**2)
        return h(theta) * Identity(3) + h(theta) / 2 * P

References#

[BBT03] (1,2,3,4)

Olivier A Bauchau and Lorenzo Trainelli. The vectorial parameterization of rotation. Nonlinear dynamics, 32:71–92, 2003. doi:10.1023/A:1024265401576.