I am interested in computing point-wise values of the normal derivative of the solution on the points of an internal boundary.
For the moment I am able to compute the integral of the normal derivative of the solution along the internal boundary, but I would like to have nodal values.
It seems to me that the problem lies in the fact that "v" (the solution) is a TestFunction(V) rather than a Function(V), for which the attribute "vector" exist.
I've added a piece of code. In this problem the domain is a square with a circle at the middle and the circumference is internal boundary.
Thank you for your help!
from dolfin import *
from mshr import *
from math import pi
# -------------------- #
# Create mesh #
# -------------------- #
rectangle = Rectangle(Point(-1., -1.), Point(1., 1.))
R = 0.5
origin = Point(0.,0.)
circle = Circle(origin, R, segments=32)
domain = rectangle
domain.set_subdomain(1, circle)
domain.set_subdomain(2, rectangle - circle)
mesh = generate_mesh(domain, 15)
# Create subdomains
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())
# -------------------- #
# Create boundaries #
# -------------------- #
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] + 1.) < DOLFIN_EPS
boundaries = FacetFunction("size_t", mesh, 0 )
Bottom().mark(boundaries, 1)
# -------------------- #
# Add interface marker #
# -------------------- #
for f in facets(mesh): # For all edges in the mesh
p0 = Vertex(mesh, f.entities(0)[0]) # save the two ending points p0 and p1
p1 = Vertex(mesh, f.entities(0)[1])
if near(p0.x(0)*p0.x(0) + p0.x(1)*p0.x(1), R*R, 0.05) and near(p1.x(0)*p1.x(0) + p1.x(1)*p1.x(1), R*R, 0.05): # check if the points lie on the circle - if yes, put a label on this edge
boundaries[f] = 4
# -------------------------------------------- #
# Mesh / Function Space / Measure / Integrals #
# -------------------------------------------- #
# 1. Save and load mesh/subdomains/boundaries
File("physical_region.xml") << subdomains
File("facet_region.xml") << boundaries
mesh = Mesh(mesh)
#subd = MeshFunction("size_t", mesh, 2)
#bound = MeshFunction("size_t", mesh, boundaries) # why is this not working?
subd = MeshFunction("size_t", mesh, "physical_region.xml")
bound = MeshFunction("size_t", mesh, "facet_region.xml")
# 2. Create Finite Element space (Lagrange P1) + Dirichlet BC
V = FunctionSpace(mesh, "Lagrange", 1)
# Define measures:
dx = Measure("dx")(subdomain_data=subd)
ds = Measure("ds")(subdomain_data=bound)
dS = Measure("dS")(subdomain_data=bound)
# Define trial/test functions
u = TrialFunction(V)
v = TestFunction(V)
# Integral on the external boundary
f0 = v*ds(1) + 1e-15*v*dx
F0 = assemble(f0)
# Integral on the internal boundary (interface)
s0 = inner(grad(v('+')),n('+'))*dS(4) + 1e-15*v*dx
s0 = assemble(s0)
