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:
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:
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}\):
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:
Available particular cases#
ExponentialMap
parametrization#
This parametrization corresponds to the simple choice:
The corresponding expression for the rotation matrix is the famous Euler-Rodrigues formula.
EulerRodrigues
parametrization#
This parametrization corresponds to the simple choice:
SineFamily
parametrization#
This generic family for any integer \(m\) corresponds to:
TangentFamily
parametrization#
This generic family for any integer \(m\) corresponds to:
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#
Olivier A Bauchau and Lorenzo Trainelli. The vectorial parameterization of rotation. Nonlinear dynamics, 32:71–92, 2003. doi:10.1023/A:1024265401576.