Define 3d vector field on a 2d mesh

General
1

Hello,
I looked up online for topics similar to this, but nothing helped.

I have a two-dimensional mesh, on which I can successfully define scalar, vectors and tensors. This mesh represents a surface in three-dimensional space, on which a 3d vector field is defined. I would like to define and plot this vector field.

Here is a minimal (non) working example

from __future__ import print_function

from fenics import *
from mshr import *
import numpy as np
import meshio
import ufl as ufl

R = 1.0
r = 0.25
c_R = [0.0, 0.0]
c_r = [0.0, -0.1]

xdmffile_geometry = XDMFFile("geometry.xdmf")

mesh = UnitSquareMesh(1, 1)

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 2)

def norm(x):
    return (np.sqrt(np.dot(x, x)))

class MyVectorFunctionExpression(UserExpression):
    def eval(self, values, x):
        values[0] = x[0]
        values[1] = -x[1]
    def value_shape(self):
        return (2,)
class MyScalarFunctionExpression(UserExpression):
    def eval(self, values, x):
        values[0] = 4*x[0]*x[1]*sin(8*(norm(np.subtract(x, c_r)) - r))*sin(8*(norm(np.subtract(x, c_R)) - R))
    def value_shape(self):
        return (1,)

def X(z):
    return as_tensor([x[0], x[1], z(x)])

z_  = Function(Q)
u_ = interpolate(MyVectorFunctionExpression(element=V.ufl_element()) ,V)
z_ = interpolate(MyScalarFunctionExpression(element=Q.ufl_element()), Q)

xdmffile_geometry.parameters.update(
    {
        "functions_share_mesh": True,
        "rewrite_function_mesh": False
    })
xdmffile_geometry.write(project(z_, Q), 0)
# xdmffile_geometry.write(project(X(z_), V), 0)

If I comment out the last line, this code runs with no errors. However, if I uncomment it, I get

$ python3 mwe.py 
Traceback (most recent call last):
  File "mwe.py", line 58, in <module>
    xdmffile_geometry.write(project(X(z_), V), 0)
  File "mwe.py", line 45, in X
    return as_tensor([x[0], x[1], z(x)])
NameError: name 'x' is not defined

Do you know what I am doing wrong ? I would like to define a 3D vectorial field X that depends on the 2d planar coordinates x,y and on a function z(x,y), where X = [x[0], x[1], z(x)], and plot this vector field just like I did for the scalar field z_ in xdmffile_geometry.write(project(z_, Q), 0). Thanks

2

Please note that your example is not reproducible, as x isn’t defined.
I would suggest providing a minimal reproducible example (using say a unit square mesh from Dolfin).

4

Please find a minimal working example in my revised post.

5

Still not a minimal working example (you are missing many definitions and code cannot be executed as we do not have the mesh)…
See illustrated code below

from fenics import *
import numpy as np
try:
    from ufl import indices
except ImportError:
    from ufl_legacy import indices


mesh = UnitSquareMesh(10, 10)
V3d = VectorFunctionSpace(mesh, 'P', 2, dim=3)
Q = FunctionSpace(mesh, 'P', 2)

def norm(x):
    return (np.sqrt(np.dot(x, x)))
def ufl_norm(x):
    return(sqrt(dot(x, x)))
c_r = 0.1
c_R = 0.3
R = 1
r=0.1
class MyScalarFunctionExpression(UserExpression):
    def eval(self, values, x):
        values[0] = x[0]*x[1]*sin(4*(norm(np.subtract(x, c_r)) - r))*sin(4*(norm(np.subtract(x, c_R)) - R))
        # values[0] = sin(norm(np.subtract(x, c_r)) - r) * sin(norm(np.subtract(x, c_R)) - R)
    def value_shape(self):
        return (1,)

i, j, k, l = indices(4)

def X(z):
    return as_tensor([x[0], x[1], z(x)])

def e(z):
    return as_tensor([[1, 0, z.dx(0)], [0, 1, z.dx(1)]])

#path for mac

z_  = Function(Q)
z_ = interpolate(MyScalarFunctionExpression(element=Q.ufl_element()), Q)

yields

Traceback (most recent call last):
  File "/root/shared/mwe_234.py", line 46, in <module>
    xdmffile_geometry.write(project(X(z_), V3d), 0)
  File "/root/shared/mwe_234.py", line 31, in X
    return as_tensor([x[0], x[1], z(x)])
NameError: name 'x' is not defined. Did you mean: 'X'?

Note that I had to make up values for r, R, c_r, c_R and change to a built-in mesh.

6

I see, I revised my original post again with a minimal working example, which runs.

7

Consider the following example:


from fenics import *
import numpy as np
import ufl as ufl

R = 1.0
r = 0.25
c_R = [0.0, 0.0]
c_r = [0.0, -0.1]


mesh = UnitSquareMesh(1, 1)

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 2)

def norm(x):
    return (np.sqrt(np.dot(x, x)))

class MyVectorFunctionExpression(UserExpression):
    def eval(self, values, x):
        values[0] = x[0]
        values[1] = -x[1]
    def value_shape(self):
        return (2,)
class MyScalarFunctionExpression(UserExpression):
    def eval(self, values, x):
        values[0] = 4*x[0]*x[1]*sin(8*(norm(np.subtract(x, c_r)) - r))*sin(8*(norm(np.subtract(x, c_R)) - R))
    def value_shape(self):
        return (1,)

def X(z):
    x = ufl.SpatialCoordinate(mesh)

    return as_tensor([x[0], x[1], z])

z_  = Function(Q)
u_ = interpolate(MyVectorFunctionExpression(element=V.ufl_element()) ,V)
z_ = interpolate(MyScalarFunctionExpression(element=Q.ufl_element()), Q)

xdmffile_geometry = XDMFFile("geometry.xdmf")

xdmffile_geometry.parameters.update(
    {
        "functions_share_mesh": True,
        "rewrite_function_mesh": False
    })
xdmffile_geometry.write(project(z_, Q), 0)

W = VectorFunctionSpace(mesh, 'P', 2, dim=3)
xdmffile_geometry.write(project(X(z_), W), 0)
8

That works like a charm, thank you. Do you know how to define a 3d vector field such that the origin of the arrow (the point where the arrow is applied) lies on the surface of the 2d mesh (embedded in 3d) and the arrow points in three-dimensional space ?

9

Do you want a vector in a single point on the manifold? What do you want to do with that vector?

10

Consider the figure attached, and say that the x, y axes define the horizontal plane, z the vertical axis. Consider only the half of the torus, and say that this is the two-dimensional manifold defining my mesh.

I want to define a vector field tangent (normal) to such two-dimensional manifold, where the two-dimensional manifold is given by my mesh, as in the figure attached. Then I would like to be able to plot this vector field as in the figure attached, an do tensorial operations with it.

Thanks

Screenshot 2024-04-08 at 20.46.54