Surface problems

7. Surface problems#

7.1. Line in plane#

7.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')
../_images/c5bcc9d7a6bac90fea94f0ff780308c565ac0534405cec4f3203d6261ed2e765.png 
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()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[4], line 1
----> 1 test_cell_orientation_curve()

Cell In[3], line 14, in test_cell_orientation_curve()
     11 g = interpolate(f, V)
     12 h = project(f, V)
---> 14 assert abs(g.dat.data - h.dat.data).max() < 1e-2
     16 print(g.dat.data - h.dat.data)

AttributeError: 'Interpolate' object has no attribute 'dat'

7.2. Surface in 3D space#

Contents