Error defining the boundary conditions: *** Warning: Found no facets matching domain for boundary condition

General
1

Hello, I am currently stuck at defining appropriate boundary conditions for my dogbone specimen geometry (kindly see the image below). I also went through some posts where the similar error had been occured like this where @dokken introduces a facetfunction where we can also visualize domain in paraview and this where @MiroK marks the facets and applies DBC therefore I modified my code accordingly but still the error persists.
My MWE is:

import meshio
from dolfin import *
import mgis.fenics as mf
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import clear_output

msh = meshio.read("/mnt/d/Research Projects/FEniCS/doggy.msh")

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
    return out_mesh

triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
meshio.write("doggy.xdmf", triangle_mesh)

# Read the modified mesh with swapped axes
modified_mesh = Mesh()
with XDMFFile("doggy.xdmf") as infile:
    infile.read(modified_mesh)

# Get the vertex coordinates of the modified mesh
coordinates = modified_mesh.coordinates()

# Swap the x and y coordinates
modified_coordinates = coordinates[:, ::-1]

# Update the vertex coordinates in the modified mesh
modified_mesh.coordinates()[:, :] = modified_coordinates

# Plot the modified mesh
plot(modified_mesh)

Vu = VectorFunctionSpace(modified_mesh, "CG", 1)
u  = Function(Vu, name="Displacement")

class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 31.75)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (x[1] == -31.75)


boundaries = MeshFunction("size_t", modified_mesh, modified_mesh.geometry().dim()-1)

top().mark(boundaries, 1) 
bottom().mark(boundaries, 2) 

Uimp = Expression(("0", "t"), t=1, degree=0)
bcu = [DirichletBC(Vu, Constant((0, 0)), boundaries, 2),
       DirichletBC(Vu, Uimp, boundaries, 1)]

v = Function(Vu)  
bcu[1].apply(v.vector()) #<--!!!This line is giving error specifically!!!

Vd = FunctionSpace(modified_mesh, "CG", 1)
d = Function(Vd, name="Damage")
dold = Function(Vd, name="Previous damage")
bcd = DirichletBC(Vd, Constant(0.), boundaries, 2)

Error:

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

Image of my domain:

image
image257×746 8.99 KB

I humbly request to kindly suggest me on how to go about this error. I appreciate any inputs regarding the same.
Thank You for the help.

2

Using

Without a tolerance in near, can be the cause.
Also

You should use near with a tolerance here as well, as == never works when you use floats. See
https://fenicsproject.org/olddocs/dolfin/latest/python/_autogenerated/dolfin.cpp.math.html?highlight=#dolfin.cpp.math.near
I.e, add a third argument to near that is bigger than 3e-16.

3

@dokken Thanks for the reply, I added the tolerances like this:

class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 31.75, 3e-16)

class bottom(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], -31.75, 3e-16)

but still this line is causing the error:

bcu[1].apply(v.vector())

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

What can possibly be the cause now ?
Thank you

4

As @dokken pointed, you should use tolerance bigger than 3e-16, so you could try something bigger, like 1e-8

1 Like
5

@Kei_Yamamoto Thanks for the suggestion, I tried with 1e-8 also, but no luck, it is giving the same error.

6

As you haven’t provided your msh/geo file, one cannot reproduce the behavior, and thus its tricky to help you any further.

7

Kindly follow this link to access my .msh file:
https://filetransfer.io/data-package/Eu8WMcCU#link

8

Hi @dokken, I tried several other methods (such as making my bottom surface as x[1] == 0, scaling the domain y-axis length to 1 unit) to make boundary conditions work but none of them worked and lead to the same error. I think the error lies in how this v vector (kindly see the code below) works when top boundary condition is applied.

v = Function(Vu)
bcu[1].apply(v.vector())

I took the code from here where my main motive is to import gmsh generated mesh and apply boundary conditions to see the effect of damage field and displacement.
I hope to get guidance regarding the same.

9

The error here is that your .msh file represents a 3D geometry, i.e. the cells are tetrahedra and the facets are triangles. You can’t convert this to a 2D mesh by simply pruning the z coordinate. You need to create a 2D mesh, i.e. one with cells as triangles and facets as lines.

Since your mesh is originally composed of tetrahedra, all the lines are shared by at least two triangles. Therefore, when you load the pruned mesh into FEniCS, on_boundary is false for every facet (line), because the test for a line to be on the boundary is that it is connected to exactly one cell. Hence your subdomains do not find any matching facets.

10

Thank you @conpierce8 for reminding me about importing a 2D mesh instead of 3D. I appreciate it.