Does anyone have an example of plotting or showing the unit normal vector on an interface between 2 subdomains of a 2D mesh?
I would like to check the direction of unit vector on an interface because I have calculated a sensitivity of an objective function with respect to the location of all the nodes in a domain consisting of 2 subdomains (using adjoint-dolfin), but I only interested in the value of the sensitivity along the interface between these two subdomains.
I compared the sensitivity calculated from the adjoint-dolfin with the sensitivity calculated from finite difference. Both exhibit the same behavior (linear) but with opposite signs.
The facet normal of an interior interface to either side as discussed here: Integrating over an interior surface
You can the project this restriction to a CG1-space and visualize it with Paraview.
Code for projecting a normal to a exterior interface has been done at: How to compute curvature of a boundary
and looks like. This can easily be modified for an interior interface.
mesh = Mesh(UnitDiscMesh.create(MPI.comm_world, 100, 1, 2))
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u,v)*ds
l = inner(n, v)*ds
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
File("nh.pvd") << nh
Thank you for your help. I could get the projected value of the unit normal with the following MWE,
from dolfin import *
# define the edges of the domain
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.5)
# define the subdomains
class Top(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5
# Create meshes on a unit square 2D domain
mesh = UnitSquareMesh(100, 100)
# Define markers for Subdomians and boundaries
markers_subdomains = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
Top().mark(markers_subdomains, 1)
Bottom().mark(markers_subdomains, 2)
markers_boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(markers_boundaries, 4) # Interface edge
# define the domains for integration
dx = Measure("dx", domain=mesh, subdomain_data=markers_subdomains)
ds = Measure("ds", domain=mesh, subdomain_data=markers_boundaries)
dS = Measure('dS', domain=mesh, subdomain_data=markers_boundaries, subdomain_id=4)
n = FacetNormal(mesh)
V = VectorFunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)
a = inner(u, v)('+')*dS
l = inner(n, v)('+')*dS
A = assemble(a, keep_diagonal=True)
L = assemble(l)
A.ident_zeros()
nh = Function(V)
solve(A, nh.vector(), L)
File("nh-positive.pvd") << nh
However, with this method, the direction of the unit normal need to be pre-defined. Is there anyway to get the default projected value without pre-defining the direction of the normal vector?
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)
Sorry for late reply. By pre-defining the direction, I mean that the direction of the unit normal is assigned first (I am not sure if I understand it correctly, but assigning “res” to " + " or " - " to me is like assigning the direction of the unit normal - please correct me if I am wrong). However, I would like to know if I do not assign (or restrict) the direction of the unit normal, by default, which direction it is pointing to?
Also, by cell function, do you mean the following part of the script?
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)
Thank you for your reply.
However, I have deleted the term + Constant(0)*inner(u,v)*dx in your code and I notice that the results with and without this term are the same. Is it supposed to be the case?