Error code while applying Dirichlet BC

1

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.

2

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
3

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?

4

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
5

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?

6

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.

7

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
8

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

boundary marker
boundary marker1472Ă—805 115 KB