Error code while applying Dirichlet BC

Hi there,

Consider a unit disk. I need to apply a constant displacement vector gamma2 = Constant((2.0,1.0)) on the nodes on the boundary of the disk. I tagged the boundary with 1 and applied the boundary condition as
c_out = DirichletBC(V, gamma2, boundary_markers, 1)
This produces the following warning:

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

I have gone through some of the posts related to this issue but none of the suggestions there are working for my case.

Here is an MWE:

C1 = Circle(Point(0, 0),1)


class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(x[0]**2 + x[1]**2, 1, eps=1e-14) and on_boundary
      


class s1(SubDomain):
    def inside(self, x, on_boundary):
        return x[0]**2 + x[1]**2 < 1 

mesh = generate_mesh(C1, 5)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
b1().mark(boundary_markers, 1)
s1().mark(surface_markers, 8)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
dx = Measure('dx', domain=mesh, subdomain_data=surface_markers)


AB = assemble(Constant(1)*dx(8))
print('area before =', AB)

# Normal vector to the intima
n = FacetNormal(mesh)

# Create mesh and define function space

V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
# Define functions
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration
u0  = Function(V) 
P = Constant(12.0)

d = mesh.geometry().dim()
print(d) 
I = Identity(d)             # Identity tensor
F = I + grad(u) 
J_e = det(F)
C_lumen = F.T*F
Ic_lumen = tr(C_lumen)

gamma2 = Constant((2.0,1.0))

alpha, mu_lumen = 11.1, 0.00000000000001
psi_lumen = (mu_lumen)*(Ic_lumen - 3) + (alpha)*(J_e-1)**2 + mu_lumen*ln(J_e) 
Pi_lumen = (psi_lumen*dx(8)) 

F1 = derivative(Pi_lumen, u, v)

# Compute Jacobian of F
Ju = derivative(F1, u, du)



c_out = DirichletBC(V, gamma2, boundary_markers, 1)

bcs=[c_out]

# Solve variational problem
solve(F1 == 0, u, bcs,
form_compiler_parameters=ffc_options)
    
ALE.move(mesh,u)

plot(mesh)
plt.show()

I am not sure where I went wrong with my defining the boundary. Any suggestions would be extremely helpful.

I haven’t managed to look at the entire code, but I think it is this line. Try a more reasonable tolerance on the boundary. Something like

class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(sqrt(x[0]**2 + x[1]**2), 1, eps=1e-1) and on_boundary

You should be able to verify if indeed you have marked the boundary facets correctly by printing boundary_markers.array() after doing b1().mark(boundary_markers, 1). If it is all zeros then you know you haven’t marked it properly.

1 Like

Hi @bhaveshshrimali ,

You were spot on. The markers are not 1 for the boundary, It’s producing an array of 0s although I marked the boundary as 1 through the code

b1().mark(boundary_markers, 1)

Your suggested code

class b1(SubDomain):
  def inside(self, x, on_boundary):
      return near(sqrt(x[0]**2 + x[1]**2), 1, eps=1e-1) and on_boundary

produces

boundary markers = [1 0 0 1 1 0 1 0 1 0 0 1 0 0 1 0 0 1 0 0 0 1 0 0 1 0 1 0 1 0 0 0 0 1 1 0 0
 1 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]

I am confused as to why this is happening. I believe

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
surface_markers = MeshFunction("size_t", mesh, mesh.topology().dim())
b1().mark(boundary_markers, 1)
s1().mark(surface_markers, 8)

is the usual way of tagging the boundaries. Any suggestion as to why this is happening?

I think your eps=1.e-14 is too tight in the definition of b1. A quick way would be to visualize the corresponding mesh function (write to paraview and see if you have marked all the facets properly there).

1 Like

I can do that ight now. But before that, I have a question.

If I am marking the boundary as 1, shouldn’t the array be producing all 1s instead of a mixture of 0s and 1s?

No, you are only marking the facets which correspond to your subdomain b1, which in this case happens to be the boundary of the circle.

Note that there are mesh.num_entities(dim) entities coressponding to dimension dim. And if I remember correctly this should also be equal to mesh.num_facets() for dim=1. Furthermore, this should also be equal to the size of boundary_markers.array(). And out of these only the ones that correspond to b1 (or are inside b1) are being marked. I can’t test this but should be a quick check.

1 Like

I did as you suggested and it seems that the boundary is correctly tagged as 1. Here is the Paraview image: