6. Surface problems

6.1. Line in plane

6.1.1. Cell orientation for line in plane

def set_cell_orientations(mesh):
    from ufl.classes import ReferenceGrad
    import firedrake as fd

    V0 = fd.FunctionSpace(mesh, 'DG', 0)
    X = fd.SpatialCoordinate(mesh)
    flag = fd.Function(V0)
    flag.interpolate(fd.dot(X, fd.as_vector((-ReferenceGrad(X)[1, 0], ReferenceGrad(X)[0, 0]))))
    cell_orientations = fd.Function(V0, dtype=np.int32)
    cell_orientations.dat.data[:] = (flag.dat.data_ro < 0)
    mesh.topology._cell_orientations = cell_orientations


def plot_orientations_1d(mesh):
    import matplotlib.pyplot as plt

    plt.figure(figsize=[4, 4])
    Vc = mesh.coordinates.function_space()
    cell_orientations = mesh.cell_orientations()
    for i, index in enumerate(Vc.cell_node_list):
        coord = mesh.coordinates.dat.data_ro_with_halos[index].real
        o = cell_orientations.dat.data_ro_with_halos[i]

        _x = coord[:, 0]
        _y = coord[:, 1]
        if o > 1/2:
            plt.arrow(_x[0], _y[0], (_x[1]-_x[0])/2, (_y[1]-_y[0])/2, head_width=0.05, head_length=0.05, fc='k', ec='k')
        else:
            plt.arrow(_x[1], _y[1], (_x[0]-_x[1])/2, (_y[0]-_y[1])/2, head_width=0.05, head_length=0.05, fc='k', ec='k')
    bbox = plt.axis('equal')

from firedrake import *
import matplotlib.pyplot as plt

# mesh = Mesh("gmsh/circle_1d.msh", dim=2)
# set_cell_orientations(mesh) 

mesh = CircleManifoldMesh(16)
x = SpatialCoordinate(mesh)
# mesh.init_cell_orientations(x)

set_cell_orientations(mesh)

V = VectorFunctionSpace(mesh, 'CG', degree=1)
n_h = Function(V, name='n_h')
n_h.project(as_vector([-x[1], x[0]]))

plt.figure(figsize=[4, 4])
for coord, vector in zip(mesh.coordinates.dat.data_ro.real, n_h.dat.data_ro.real):
    plt.arrow(coord[0], coord[1], 0.2*vector[0], 0.2*vector[1], head_width=0.05, head_length=0.05, fc='k', ec='k')
bbox = plt.axis('equal')

firedrake:WARNING OMP_NUM_THREADS is not set or is set to a value greater than 1, we suggest setting OMP_NUM_THREADS=1 to improve performance

def test_cell_orientation_curve():
    m = CircleManifoldMesh(3)
    x = SpatialCoordinate(m)
    # m.init_cell_orientations(x)
    set_cell_orientations(m)
    
    U = VectorFunctionSpace(m, 'CG', degree=1)
    V = VectorFunctionSpace(m, 'CG', degree=2)
    f = project(CellNormal(m), U)

    g = interpolate(f, V)
    h = project(f, V)

    assert abs(g.dat.data - h.dat.data).max() < 1e-2
    
    print(g.dat.data - h.dat.data)

test_cell_orientation_curve()

/Users/zzyang/opt/firedrake/firedrake-real-int32-debug/src/firedrake/firedrake/interpolation.py:385: FutureWarning: The use of `interpolate` to perform the numerical interpolation is deprecated.
This feature will be removed very shortly.

Instead, import `interpolate` from the `firedrake.__future__` module to update
the interpolation's behaviour to return the symbolic `ufl.Interpolate` object associated
with this interpolation.

You can then assemble the resulting object to get the interpolated quantity
of interest. For example,

```
from firedrake.__future__ import interpolate
...

assemble(interpolate(expr, V))
```

Alternatively, you can also perform other symbolic operations on the interpolation operator, such as taking
the derivative, and then assemble the resulting form.

  warnings.warn("""The use of `interpolate` to perform the numerical interpolation is deprecated.

[[ 1.38777878e-16  1.11022302e-16]
 [ 5.55111512e-16 -4.79555877e-17]
 [-2.77555756e-16  4.44089210e-16]
 [ 0.00000000e+00  4.62120667e-17]
 [-3.33066907e-16 -4.44089210e-16]
 [ 8.32667268e-17  0.00000000e+00]]

6.2. Surface in 3D space