Discontiuity of divergence and curl of known function on surface

[Alpha-Legion]

I have a smooth surface S embedded in 3D and a smooth vector field u_3D defined on the 3D space.
I project u_3D onto S and obtained a tangent vector field u_2D of S.
Then I want to compute the divergence and curl of u_2D.
When S is a sphere modeled by the built-in mesh of Firedrake, the results seems OK.
But when S is a surface modeled by a mesh generated by myself, the div(u_2D) and curl(u_2D) are not smooth.

Since the quality of my mesh is fine as I feel, I guess the problem might be related to function spaces, about which I might misunderstand. But I am really not sure about this. Did I make any mistake? Any hint will be helpful. Thank you very much!
My codes are as follows:

from firedrake import *
import numpy as np

data_file = 'S.msh'
# mesh & coordinates
mesh = Mesh(data_file, dim=3)

# function spaces
deg = 1
SRT = FunctionSpace(mesh, 'RT', deg) # 1-forms (closed for degree = 1)
SRTE = FunctionSpace(mesh, 'RTE', deg) # 1-forms (coclosed for 1</= degree </=2)
SCG = FunctionSpace(mesh, 'CG', deg) # scalars
VCG = VectorFunctionSpace(mesh, 'CG', deg) # vectors

# orientation from analytic expression
X = SpatialCoordinate(mesh)
x, y, z = X
r2 = x**2 + y**2
f1 = (x**2 - y**2 - 1)**2 + (2*x*y)**2 - 1
f2 = 30 * sqrt(r2) * (z**2 - 0.2**2)
dFun = as_vector([(f2 + 8*f1*(r2 - 1)) * x,
                  (f2 + 8*f1*(r2 + 1)) * y,
                  2*(10 * r2**1.5 + 1) * z])
dFun_norm = sqrt(inner(dFun, dFun))
orientation = project(dFun/dFun_norm, VCG, name='Orientation')
mesh.init_cell_orientations(orientation)
normal = project(CellNormal(mesh), VCG, name='Normal')

# define vector field
s1 = 2 # scaling factor
u_3D = as_vector([
sin(s1*z) + cos(s1*y),
sin(s1*x) + cos(s1*z),
sin(s1*y) + cos(s1*x)
])

# project to tangent space
u_N = inner(normal, u_3D)
u_2D = u_3D - u_N * normal
u2D = project(u_2D, VCG, name='u_2D')

u_H1 = project(u_2D, VCG, name='H^1 input vector field')
#u_Hdiv = project(u_2D, SRT, name='H(div) input vector field')
#u_Hcurl = project(u_2D, SRTE, name='H(curl) input vector field')

# differentiation
div_u = project(div(u_H1), SCG, name='div_u')
curl_u = project(curl(u_H1), VCG, name='curl_u')
curl2_u = project(inner(normal, curl_u), SCG, name='curl2_u')

# output to file
VTKFile('Double-Torus_HHD_inputs.pvd').write(normal, u2D, u_H1, u_Hdiv, u_Hcurl, div_u, curl2_u, star_u)

I've tried with

SDG = FunctionSpace(mesh, 'DG', deg) # scalars
VDG = VectorFunctionSpace(mesh, 'DG', deg) # vectors

and with degree deg-1 as well, but the visualizations of the results are similar.

[colinjcotter]

The divergence of an Hdiv finite element space is in DG, so if you want to properly visualise it you need to project to the DG space of the appropriate size instead, and similar for the Hcurl space.

[dham]

I think the reason you don't observe the same thing for the built-in sphere mesh is due to the regularity of the mesh structure in that case.

[Alpha-Legion]

The divergence of an Hdiv finite element space is in DG, so if you want to properly visualise it you need to project to the DG space of the appropriate size instead, and similar for the Hcurl space.

Thanks for the help. Sorry, I forgot to mention that I've tried with

SDG = FunctionSpace(mesh, 'DG', deg) # scalars
VDG = VectorFunctionSpace(mesh, 'DG', deg) # vectors

and with degree deg-1 as well, but the visualizations of the results are similar.

[Alpha-Legion]

I think the reason you don't observe the same thing for the built-in sphere mesh is due to the regularity of the mesh structure in that case.

Thanks. It's possible. But if that is the case, then the requirement for the regularity of mesh is really beyond my expectation.