# How to determine the angle made by the basis vectors of a tangent plane to the vertices of a mesh?

Hi! In the below MWE, I have a unit square mesh. At each vertex of the mesh, I am defining a tangent plane at each vertex of the mesh whose basis vectors are \mathbf{e}_{1} and \mathbf{e}_{2}. The normal vector to the tangent plane is \mathbf{e}_{3} as shown in Linear shell model — Numerical tours of continuum mechanics using FEniCS master documentation (comet-fenics.readthedocs.io)

I need to determine the angle made by \mathbf{e}_1 with the X axis for each vertex of the mesh. For a square mesh, the answer is trivial. But my actual problem deals with a very complicated 2D geometry. My other question is how do I calculate the angles made by \mathbf{e}_1 at each vertex with the X-axis and construct a rotation matrix using the angle at each vertex?

MWE:

from dolfin import *
import meshio
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from ufl import cofac, sqrt
from scipy import linalg as la
from sympy import symbols
from sympy import atan2,Abs,cos,sin
import math
from numpy import array, zeros
from ufl import *
from ufl.classes import *
from ufl import Jacobian, replace
from ufl.classes import FacetJacobian, ReferenceGrad, ReferenceValue

mesh = UnitSquareMesh(8,8)

def local_frame(mesh):
t = Jacobian(mesh)
if mesh.geometric_dimension() == 2:
t1 = as_vector([t[0, 0], t[1, 0], 0])
t2 = as_vector([t[0, 1], t[1, 1], 0])
else:
t1 = as_vector([t[0, 0], t[1, 0], 0])
t2 = as_vector([t[0, 1], t[1, 1], 0])
e3 = cross(t1, t2)
e3 /= sqrt(dot(e3, e3))
ey = as_vector([0, 1, 0])
ez = as_vector([0, 0, 1])
e1 = cross(ey, e3)
norm_e1 = sqrt(dot(e1, e1))
e1 = conditional(lt(norm_e1, 0.5), ez, e1 / norm_e1)

e2 = cross(e3, e1)
e2 /= sqrt(dot(e2, e2))
return e1, e2, e3

frame = local_frame(mesh)
e1, e2, e3 = frame

VT = VectorFunctionSpace(mesh, "DG", 0, dim=3)
FT = VectorFunctionSpace(mesh, "P", 1)

ffile = XDMFFile("output.xdmf")
ffile.parameters["functions_share_mesh"] = True
for (i, ei) in enumerate(frame):
ei = Function(VT, name="e{}".format(i + 1))
ei.assign(project(frame[i], VT))
#e1n = project(Expression(('1.0','0.0', '0.0'), degree=0), FT)
ffile.write(ei, 0)


In Paraview, \mathbf{e}_{1} looks like (I set “no orientation array” for the glyphs. However, if I set it orientation to \mathbf{e}_{1}, the vectors no longer point all in the same direction. I don’t know the reason yet)

I would like to calculate the angle made by the vectors with respect to the X-axis at each vertex (in this case, it is obvious that the answer is 0. But my actual geometry is a very complex 2D one.) and define a rotation matrix using that angle at each vertex. I am not sure how to do that by not accessing the components of \mathbf{e}_{1}.

Any help would be much appreciated. Let me know if you have any questions.

Hi, if you need a rotation matrix from the global (\mathbf{X}, \mathbf{Y}) frame to (\mathbf{e}_1, \mathbf{e}_2), you can probably use something like:

ex = ufl.as_vector([1, 0])
ey = ufl.as_vector([0, 1])
R = ufl.outer(e1, ex) + ufl.outer(e2, ey)


Otherwise, if you really need the angle, then if \mathbf{e}_1 = \cos\theta\mathbf{e}_X+\sin\theta\mathbf{e}_Y, you can have: \theta=\arctan((\mathbf{e}_1\cdot\mathbf{e}_Y)/(\mathbf{e}_1\cdot\mathbf{e}_X)) which in UFL would be:

theta = ufl.atan_2(e1[1],e1[0])

2 Likes

@bleyerj Thanks for the solution. I think atan2 needs to be changed to atan_2. I thought it was resolved in this post ( atan2 vs atan_2 - documentaion error · Issue #102 · FEniCS/ufl (github.com)). But ufl.atan2() returns the following error:

AttributeError: module 'ufl' has no attribute 'atan2'. Did you mean: 'atan'?


Please let me correct if I am wrong.

I defined theta and the rotation matrix as follows:

theta=ufl.atan_2(frame[0][1],frame[0][0])
R = as_tensor([[cos(theta),-sin(theta),0.0],[sin(theta),cos(theta),0.0],[0.0,0.0,1.0]])


I think it will give me the same R that you defined.

Right, I edited the solution above.

Hi @bleyerj ,

I just realized that the orientation of the vector \mathbf{e}_{1} isn’t always along the axis of the body for any arbitrary geometry. For example, for the below geometry

I want the orientation of \mathbf{e}_{1} to be along the axis of the body (the red arrow in the below figure)

Instead, the orientation of \mathbf{e}_{1} is as follows (I selected a few points of the mesh for displaying \mathbf{e}_{1})

Is there any way I can modify the below code to get the orientation of \mathbf{e}_{1} along the axis of the beam. In my real problem, the beam might get bent over time. So the axis is not always a straight line and might be hard to get an expression for the axis.

msh = meshio.read("singleslantchain.msh")

######################################################################################################################

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

######################  Defining vector and scalar function spaces  ##################

F_bact = VectorFunctionSpace(mesh, 'P', 1)
V_bact = VectorFunctionSpace(mesh, 'P', 1)
T_bact = TensorFunctionSpace(mesh, 'P', 1)

######################  DEFINING INITIAL ROTATION MATRIX  ##############################################

def local_frame(mesh):
t = Jacobian(mesh)
if mesh.geometric_dimension() == 2:
t1 = as_vector([t[0, 0], t[1, 0], 0])
t2 = as_vector([t[0, 1], t[1, 1], 0])
else:
t1 = as_vector([t[0, 0], t[1, 0], 0])
t2 = as_vector([t[0, 1], t[1, 1], 0])
e3 = cross(t1, t2)
e3 /= sqrt(dot(e3, e3))
ey = as_vector([0, 1, 0])
ez = as_vector([0, 0, 1])
e1 = cross(ey, e3)
norm_e1 = sqrt(dot(e1, e1))
e1 = conditional(lt(norm_e1, 0.5), ez, e1 / norm_e1)

e2 = cross(e3, e1)
e2 /= sqrt(dot(e2, e2))
return e1, e2, e3

frame = local_frame(mesh)
e1, e2, e3 = frame

ex = ufl.as_vector([1, 0, 0])
ey = ufl.as_vector([0, 1, 0])

theta = ufl.atan_2(e1[1],e1[0])


I need \theta to be angle between the X-axis and the axis of the beam (the red arrow in the previous figure; the axis might be a curved line for a curved beam)

In case you need the Gmsh file for the geometry, here it is:

SetFactory("OpenCASCADE");
//+
Point(1) = {0.0, 0.0, 0, 1.0};
//+
Point(2) = {5.0*Cos(Pi/4), 5.0*Sin(Pi/4), 0, 1.0};
//+
Point(3) = {5.0*Cos(Pi/4), 5.0*Sin(Pi/4)+0.5, 0, 1.0};
//+
Point(4) = {0.0, 0.5, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {3, 4, 1, 2};
//+
Plane Surface(1) = {1};
//+
Physical Curve("left", 1) = {4};
//+
Physical Curve("right", 2) = {2};
//+
Physical Curve("top", 3) = {3};
//+
Physical Curve("bottom", 4) = {1};
//+
Physical Surface("body", 5) = {1};


Please let me know if you have any ideas on how to achieve that. Also, please let me know if you have any questions.

P.S. I did not start a new topic because I thought that this would be a nice extension of the existing topic.

You will need to find some way to define it from your geometry for the general case.

@bleyerj Yes, you are right. Finding the axis would have made life easy. However, I was thinking of an interpolation technique illustrated in the schematic diagram below. Suppose we have a curved beam. For every point on the curves 3 and 4, we find the unit tangent vectors and then interpolate in the region in between. So for every point on the domain, we will get unit vectors that are close to the direction of the axis (the mathematical form of the axis is not known). We calculate the angle made by the vectors at each node/vertex with the X-axis and define a rotation matrix for each vertex using the angle. We don’t need to define the axis then. Let me know what you think of this.

Schematic diagram: