How to plot normal unit vector of faces in a 2D mesh?

4

Not sure what you mean by pre-defining the direction?
Anyhow. You do not use your cell-function in this case, which does not restrict your normal as it should.
Here is an example where i compute the normal between two interface. Changing res from “+” to “-” will change the direction of the normal as described in:

from dolfin import *
import numpy as np

# define the subdomains
class Circle(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2  + x[1]** 2 <= 0.5

# Create meshes on a unit square 2D domain
mesh = UnitSquareMesh(15, 15)
mesh.init(1,2)
facet_to_cells = mesh.topology()(1,2)
m0, m1 = 2,3
# Define markers for Subdomians and boundaries
markers_subdomains = MeshFunction("size_t", mesh, mesh.geometry().dim(), m0)
Circle().mark(markers_subdomains, m1)
File("cf.pvd").write(markers_subdomains)

# Find interface given the two cell markers
marker_array = markers_subdomains.array()
i_facets = []
for facet in facets(mesh):
    f_index = facet.index()
    cells = facet_to_cells(f_index)
    if len(cells) >= 2:
        if np.isin([m0,m1], marker_array[cells]).all():
            print(f_index)
            i_facets.append(f_index)
interface_marker = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
interface_marker.array()[i_facets] = 2

# define the domains for integration
dx = Measure("dx", domain=mesh, subdomain_data=markers_subdomains)
dS = Measure('dS', domain=mesh, subdomain_data=interface_marker, subdomain_id=2)

n = FacetNormal(mesh)
res = "-"
V = VectorFunctionSpace(mesh, "DG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u(res), v(res))*dS + Constant(0)*inner(u,v)*dx
l = inner(n(res), v(res))*dS

A = assemble(a, keep_diagonal=True)
L = assemble(l)

A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
XDMFFile("nh-positive.xdmf").write_checkpoint(nh, "nh", 0, append=False)
1 Like