Finite rotation library
This file implements the functions that allows to deal with finite rotations tensors in dolfin.
First of all, import dolfin
from dolfin import *
import ufl
# and immediately after that all the constants that allows compute
# the scalar coefficients involved in the closed form of the rotation
# formulae.
# The file is documented here: :doc:`RotCoef.py`
# ::
import RotCoef
At this point we need to define some convenience functions:
Axial vector of a second-order tensor:
given the second order tensor \(\boldsymbol{T}\) returns its axial vector \(\boldsymbol{v}\) such that \(\boldsymbol{v} \times = \frac{1}{2}(\boldsymbol{T} - \boldsymbol{T}^\mathrm{T})\).
Note that, since the input is a second-order tensor, and the output must be vector, the components \(T_{i,j}\) of tensor \(\boldsymbol{T}\) are extracted with
T[i, j], and vector \(\boldsymbol{v}\) is built using the functionas_vector
def axial(Tensor):
return as_vector([
(Tensor[1, 2] - Tensor[2, 1])*0.5,
(Tensor[2, 0] - Tensor[0, 2])*0.5,
(Tensor[0, 1] - Tensor[1, 0])*0.5])
Symmetric part of a second-order tensor:
compute the symmetric part \(\boldsymbol{S}\) of tensor \(\boldsymbol{T}\) : \(\boldsymbol{S} = \frac{1}{2}(\boldsymbol{T} +\boldsymbol{T}^\mathrm{T})\)
def symm(Tensor):
return (Tensor + Tensor.T) * 0.5
Skew-symmetrice tensor:
given vector \(\boldsymbol{v}\) returns the skew-symmetric tensor \(\boldsymbol{v}\times\).
Here, since the argument is a vector and the returned entity is a function, the tensor is built using the function
as_tensor, whith the components \(v_i\) of vector \(\boldsymbol{v}\) extracted asV[i]
def crossTensor(V):
return as_tensor([
[0.0, V[2], -V[1]],
[-V[2], 0.0, V[0]],
[V[1],-V[0], 0.0]])
Norm ov vector \(\boldsymbol{v}\):
given vector \(\boldsymbol{v}\) returns \(||\boldsymbol{v}||=\sqrt{\boldsymbol{v}\boldsymbol{v}}\)
def vnorm(V):
return sqrt(dot(V, V))
Helper function:
computes the Maclaurin expansion of coefficient
idx
def coef_pow_expansion(idx, phi2):
coef = 0.;
for i in range(RotCoef.SerExpTrunc[idx]-1, -1, -1):
coef = coef * phi2 + RotCoef.SerExpCoeffs[idx][i]
return coef
Helper functions:
compute the the coefficient \(a=\frac{\sin(\varphi)}{\varphi}\), \(b=\frac{1-\cos(\varphi)}{\varphi^2}\), \(c=\frac{1-\sin(\varphi)/\varphi}{\varphi}\) where \(\varphi=||\boldsymbol{\varphi}||\);
Since their expression is ill-defined for \(\varphi\rightarrow 0\) the coefficients are computed using their Maclaurin expansion for \(\varphi^2<\overline{\varphi}^2\), where the threshold \(\overline{\varphi}^2\) is stored into
RotCoef.SerExpThrsh.The use of the UFL
conditionalexpression allows to propagate the conditional evaluation down to the form expression, thus selecting at run-time the correct expression, either the Maclaurin expansion or the trigonometric one.
def coefa(phi2):
return conditional(le(phi2, RotCoef.SerExpThrsh[0]),
coef_pow_expansion(0, phi2),
sin(sqrt(phi2))/sqrt(phi2))
def coefb(phi2):
return conditional(le(phi2, RotCoef.SerExpThrsh[1]),
coef_pow_expansion(1, phi2),
(-cos(sqrt(phi2))+1.)/phi2)
def coefc(phi2):
return conditional(le(phi2, RotCoef.SerExpThrsh[2]),
coef_pow_expansion(2, phi2),
(-sin(sqrt(phi2))/sqrt(phi2)+1.)/phi2)
Rotation tensor:
given vector \(\boldsymbol{\varphi}\) returns the othogonal tensor \(\boldsymbol{\varPhi}\in\mathrm{SO}(3)\)
\(\boldsymbol{\varPhi}=\exp(\boldsymbol{\varphi}\times)= \boldsymbol{I} + a(\varphi)\boldsymbol{\varphi}\times + b(\varphi)\boldsymbol{\varphi}\times\boldsymbol{\varphi}\times\).
def RotTensor(phi):
PhiCross = crossTensor(phi)
phi2 = dot(phi, phi)
a = coefa(phi2)
b = coefb(phi2)
return Identity(3) + PhiCross * a + dot(PhiCross, PhiCross) * b
Rotation differential tensor:
Given vector \(\boldsymbol{\varphi}\) returns the othogonal tensor \(\boldsymbol{\varGamma}\) such that \(\boldsymbol{\varphi}_\delta\times= \delta\boldsymbol{\varPhi}\boldsymbol{\varPhi}^{\mathrm{T}}=\boldsymbol{\varGamma}\delta \boldsymbol{\varphi}\)
\(\boldsymbol{\varGamma}=\mathrm{dexp}(\boldsymbol{\varphi}\times) = \boldsymbol{I}+b(\varphi)\boldsymbol{\varphi}\times+c(\varphi)\boldsymbol{\varphi}\times\boldsymbol{\varphi}\times\).
def RotDiffTensor(phi):
PhiCross = crossTensor(phi)
phi2 = dot(phi, phi)
b = coefb(phi2)
c = coefc(phi2)
return Identity(3) + PhiCross * b + dot(PhiCross, PhiCross) * c
A compiled expression:
it changes the rotation vector into its complementary whenever its modulus is greater than \(\pi\).
The expression dimension is equal to 6:
unwrap_phi() : Expression(6):the first three components should hold the displacement vector, the last three components the rotation vector. A
std::shared_ptr<Function> phi_functionpublic variable stores a pointer to the
Functionstoring the unknowns of the problem, with three position components followed by the three rotation components that should be eventually changed into their complementary.Finally, the
evalmethod returns the result computed at positionx.The value of
Functionphi_functionat positionxmust first be computed:the local
Arrayphi, with length equal to 6, is declared; it will be used to store the values computed byphi_function:Array<double> phi(6);the
FunctionSpacepointer of functionphi_functionis retrieved using thefunction_space()method. From it, the mesh pointer retrieved by using themesh()method;a
Cellis built from the mesh knowing the indexufc.cell_indexof the cell where theExpressionneeds to be evaluated:phi_function->eval(phi, x, dolfin_cell_phi, ufc_cell);evaluate (
eval) theFunctionphi_functionat positionx, and store the result into theArray<double> phi(6)variable:phi_function->eval(phi, x, dolfin_cell_phi, ufc_cell);
It is now possible to store into the return vector
valuesthe unchanged three position components, followed by the components of the rotation vector; the latter are left unchanged if the their modulus is less than \(\pi\), otherwise are changed into the complementary rotation vector.The complete expression code follows:
cppcode_unwrap_phi = """
#include <pybind11/pybind11.h>
#include <dolfin/mesh/Cell.h>
#include <dolfin/function/FunctionSpace.h>
#include <dolfin/function/Function.h>
#include <dolfin/function/Expression.h>
namespace dolfin {
class unwrap_phi : public Expression {
public:
unwrap_phi() : Expression(6) {}
std::shared_ptr<Function> phi_function;
void eval(Array<double>& values, const Array<double>& x, const ufc::cell& ufc_cell) const {
Array<double> phi(6);
constexpr double pi = std::acos(-1);
Cell dolfin_cell_phi(*phi_function->function_space()->mesh(), ufc_cell.index);
phi_function->eval(phi, x, dolfin_cell_phi, ufc_cell);
double phi_norm = std::sqrt(phi[3] * phi[3] + phi[4] * phi[4] + phi[5] * phi[5]);
if (phi_norm > pi) {
constexpr double pi2 = 2. * std::acos(-1);
const double c = 1. - pi2 / phi_norm;
values[0] = phi[0];
values[1] = phi[1];
values[2] = phi[2];
values[3] = phi[3] * c;
values[4] = phi[4] * c;
values[5] = phi[5] * c;
} else {
values[0] = phi[0];
values[1] = phi[1];
values[2] = phi[2];
values[3] = phi[3];
values[4] = phi[4];
values[5] = phi[5];
}
}
};
}
PYBIND11_MODULE(SIGNATURE, m) {
pybind11::class_<dolfin::unwrap_phi, std::shared_ptr<dolfin::unwrap_phi>, dolfin::Expression>
(m, "unwrap_phi")
.def(pybind11::init<>())
.def_readwrite("phi_function", &dolfin::unwrap_phi::phi_function)
;
}
"""
Rotation vector:
given the othogonal tensor \(\boldsymbol{\varPhi}\in\mathrm{SO}(3)\) returns the rotation vector \(\boldsymbol{\varphi}\) such that \(\exp(\boldsymbol{\varphi}\times)=\boldsymbol{\varPhi}\) This method is not used in the code, and is provided for completeness. It makes heavy use of conditionals.
def RotVector(PHI):
CosPhi = (tr(PHI)-1.0)*0.5
def CosPhi1():
V = axial(PHI)
SinPhi = vnorm(V)
b = ufl.atan_2(SinPhi, CosPhi)
a = coefa(b*b)
V = V / a
return V
def CosPhi2():
T = (symm(PHI) - Identity(3) * CosPhi) / (1. - CosPhi)
V = conditional( ufl.And(ge(T[0, 0], T[1, 1]), ge(T[0, 0], T[2, 2])),
as_vector(T[:, 0] / sqrt(T[0, 0])),
conditional(ge(T[1, 1], T[2, 2]),
as_vector(T[:, 1] / sqrt(T[1, 1])),
as_vector(T[:, 2] / sqrt(T[2, 2])),
)
)
SinPhi = tr(dot(crossTensor(V), PHI)) * (-0.5)
V = V * ufl.atan_2(SinPhi, CosPhi)
return V
V = conditional(gt(CosPhi, 0.), CosPhi1(), CosPhi2())
return V