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.
