Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License

The corresponding files can be obtained from:

A module for rotation parametrization

This document describes the implementation of generic vectorial parametrization of rotation matrices, essentially following [BAU03]. 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}\), [BAU03] 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:

\begin{equation} \varphi = \sqrt{\boldsymbol{\theta}^2 + \varepsilon} \end{equation}

with \({\varepsilon}=\) DOLFIN_EPS in practice.

The rotation parameter vector \(\boldsymbol{p}\) from [BAU03] is given by the rotation_parameter attribute.

The class then implements the following functions:

\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 [BAU03].

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

\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:

\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:

\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:

\begin{equation} p(\varphi)=2 sin(\varphi/2) \end{equation}

SineFamily parametrization

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

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

TangentFamily parametrization

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

\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 dolfin import Constant, DOLFIN_EPS


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):
        m = Constant(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):
        m = Constant(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

[BAU03] Bauchau, O. A., & Trainelli, L. (2003). The vectorial parameterization of rotation. Nonlinear dynamics, 32(1), 71-92.

[ ]: