*** Warning: Found no facets matching domain for boundary condition

mesh
1

What does this error mean?

The full warning message is:

*** Warning: Found no facets matching domain for boundary condition.
Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
Newton solver finished in 0 iterations and 0 linear solver iterations.

Does it mean that I need to mark specific subdomains as part of the boundary condition?

Code:

from fenics import *

"""
    Simulating Magnetic field from a square loop of wire.
    Box is 5x5x1. The wire with lower right hand corner on the origin and at
    height z = 0.5.
    Equation:

    curl (curl A/mu) = 0
    curl A = B (what we're looking for)
    At the boundary (on the wire), B is defined.
"""

#create mesh and define function space
# box is 5x5x1
length=width=5.0
height = 1.0
mesh = BoxMesh(Point(0., 0., 0.), Point(length, width, height), 10, 10, 2)
V = VectorFunctionSpace(mesh, 'P', 1)

loop_side = 1.0 #length of square loop
C = 0.1 #constant (~Remanance/4pi) for boundary eqn component
tol = 1E-14 #tolerance for boundary definition
loop_height_z = 0.5 #the sq loop will be at this height
mu = 1.32E-6; #material mu for neodymium magnet

#-----define boundary components-----#
# x component at boundary
b_D_x_str = "2*K*x[2]*(1/pow(pow(l-x[0], 2) + pow(x[2], 2), 2) - 1/pow(pow(x[0], 2) + pow(x[2], 2), 2))"
# y component at boundary
b_D_y_str = "2*K*x[2]*(1/pow(pow(l-x[1], 2) + pow(x[2], 2), 2) - 1/pow(pow(x[1], 2) + pow(x[2], 2), 2))"
# z component at boundary
b_D_z_str = """-2*K*((x[0]-l)/pow(pow(l-x[0],2) + pow(x[2], 2), 2) + x[1]/(pow(pow(x[1], 2) + pow(x[2], 2), 2)) +
x[0]/(pow(pow(x[0], 2) + pow(x[2], 2), 2)) -
(x[1]-l)/pow(pow(l-x[1],2) + pow(x[2], 2), 2))"""
#-----END define boundary components-----#

b_D = Expression((b_D_x_str, b_D_y_str, b_D_z_str), degree=3, l=loop_side, K=C)

""" Checks if vector x is on the boundary: a square loop with lower left hand
    corner on the origin. Then shifted up by loop_height_z in the z direction"""
def on_boundary(x):
    if near(x[0], 0., tol) or near(x[0], loop_side, tol): #check if 0<x<loop_side
        if near(x[1], 0., tol) or near(x[1], loop_side, tol): #check if 0<y<loop_side
            return near(x[2], loop_height_z, tol) #check if z ~= loop_height_z
    return False

#boundary condition
bc = DirichletBC(V, b_D, on_boundary)

#defining variational problem
B = Function(V); #B field
s = TestFunction(V); #test function
f = Constant((0., 0., 0.)) # zero on RHS
a = (1/mu)*dot(B, curl(s))*dx # Left hand side
F = PETScVector() #right hand side
assemble(inner(f,s)*dx, tensor = F)

#compute solution
B = Function(V)
solve(a==F, B, bc)
2

This means that your on_boundary function always returns False.
I tested your if test, and it does not mark anything.
You could test this by creating a facetfunction that can be visualized in Paraview:

from dolfin import *


#create mesh and define function space
length=width=5.0
height = 1.0
mesh = BoxMesh(Point(0., 0., 0.), Point(length, width, height), 10, 10, 2)
V = VectorFunctionSpace(mesh, 'P', 1)

loop_side = 1.0 #length of square loop
C = 0.1 #constant (~Remanance/4pi) for boundary eqn component
tol = 1E-14 #tolerance for boundary definition
loop_height_z = 0.5 #the sq loop will be at this height

class bndry(SubDomain):
    def inside(self, x, on_boundary):
        # if (near(x[0], 0., tol) or near(x[0], loop_side, tol)) and \
        #    (near(x[1], 0., tol) or near(x[1], loop_side, tol)) and \
        #    (near(x[2], loop_height_z, tol) and on_boundary):
        #     return True
        # else:
        #     return False
        return on_boundary and x[0] <= length/2 and x[1] < width/2 and x[2]< 0.5*height
facetfunc = MeshFunction("size_t", mesh, mesh.geometry().dim()-1)
bnd = bndry()
bnd.mark(facetfunc, 3)

File("facetfunc.pvd").write(facetfunc)
1 Like
3

Thank you for the help! I don’t think I would’ve been able to figure out that bndry class. Installing Paraview and seeing the boundaries helped immensely!