Explicitly defining the subdomains and boundaries return warning Found no facets matching domain for boundary condition

1

Hello everyone,

I am trying to solve for a fluid flow problem in a cylinder that is in a rectangular box (which representing solid). I define subdomains for the fluid and solid and the boundaries for both subdomains. After running the calculation, I got the warning *** Warning: Found no facets matching domain for boundary condition. and the result is ridiculously high. I suspect that the way I define the subdomians and the boundaries are the cause of the problem. Could anyone help me point out what is wrong in my subdomians and the boundaries definitions?

The following is the MWE that reproduces the warning:

from fenics import *

# classes for subdomains
class fluid(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.0, 2.0)) and (x[1]**2 + x[2]**2) <= 0.5625  

class solid(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.0, 2.0)) and (x[1]**2 + x[2]**2) >= 0.5625

# classes for boundaries
class fluid_in(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) and (x[1]**2 + x[2]**2) <= 0.5625 

class fluid_side(SubDomain):
    def inside(self, x, on_boundary):
        return between(x[0], (0.0, 2.0)) and near(x[1]**2 + x[2]**2, 0.5625)

class fluid_out(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 2.0) and (x[1]**2 + x[2]**2) <= 0.5625    

class solid_in(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) and (x[1]**2 + x[2]**2) >= 0.5625

class solid_sides(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0) and near(x[2], -1.0) and near(x[2], 1.0)

class solid_bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -1.0)

class solid_out(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 2.0) and (x[1]**2 + x[2]**2) >= 0.5625

# Varying cases
Re   = Constant(5E-1, name="Re")   #500
Pr   = Constant(100, name="Pr")    #0.01, 100
kCHT = Constant(20.0, name="kCHT") #1.0, 5.0, 20.0

# Solid properties
rhoS = Constant(1.0, name="rhoS")                  # density
cpS  = Constant(100.0, name="cpS")                 # heat capacity
cdS  = Constant(100.0, name="cdS")                 # heat conductivity
alpS = Constant(cdS/(rhoS*cpS), name="alpS")       # heat diffusivity

# Parameters for the problem
U_avg = Constant(0.1, name="U_avg")                # inlet velocity
L_ch = Constant(3.5, name="L_channel")             # length of the plate
h_ch = Constant(0.25, name="h_channel")            # half of the channel height

# Fluid properties
cdF  = Constant(cdS/kCHT, name="cdF")              # heat conductivity
rhoF = Constant(1.0, name="rhoF")                  # density
muF  = Constant(rhoF*U_avg*2*h_ch/Re, name="muF")  # dynamic viscosity
cpF  = Constant(cdF*Pr/muF, name="cpF")            # heat capacity
alpF = Constant(cdF/(rhoF*cpF), name="alpF")       # heat diffusivity

# Constant for Bcs
P_in  =  Constant(3*muF*L_ch*U_avg/(h_ch*h_ch), name="P_inlet")
P_out =  Constant(0.0, name="P_outlet")

# create mesh
mesh = BoxMesh(Point(0, -1, -1), Point(2, 1, 1), 20, 20, 20)

# Define markers for Subdomians, Dirichlet and Neuman bcs
subdomain_marker = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
fluid().mark(subdomain_marker, 1)
solid().mark(subdomain_marker, 2)

boundary_marker = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
fluid_in().mark(boundary_marker, 1)
fluid_side().mark(boundary_marker, 2)
fluid_out().mark(boundary_marker, 3)
solid_in().mark(boundary_marker, 4)
solid_sides().mark(boundary_marker, 5)
solid_bottom().mark(boundary_marker, 6)
solid_out().mark(boundary_marker, 7)

dx = Measure("dx", domain=mesh, subdomain_data=subdomain_marker)
ds = Measure("ds", domain=mesh, subdomain_data=boundary_marker)

V  = VectorElement("CG",  mesh.ufl_cell(), 2)       # Velocity 
P  = FiniteElement("CG",  mesh.ufl_cell(), 1)       # Pressure 
element = MixedElement([V, P])
W = FunctionSpace(mesh, element)

# Define boundary conditions for fluid problem 
bc1 = DirichletBC(W.sub(0), Constant((0.0, 0.0, 0.0)), boundary_marker, 2) 
bcs = [bc1]

# define trial functions and test functions
(u, p) = TrialFunctions(W) 
(v, q) = TestFunctions(W)

# unit normal vector on the surfaces of the fluid mesh
n = FacetNormal(mesh)

# weak form for the Stroke-flow
a = muF*inner(grad(u), grad(v))*dx(1)\
    - p*div(v)*dx(1)\
    - q*div(u)*dx(1)

l = - P_in*inner(n, v)*ds(1)\
    - P_out*inner(n, v)*ds(3)

# Assembly the lhs and rhs matrices
A = assemble(a, keep_diagonal=True)
L = assemble(l)

# Apply boundary conditions
[bc.apply(A) for bc in bcs]
[bc.apply(L) for bc in bcs]

# solve for the solution
A.ident_zeros()
w = Function(W, name="fluid variables")

# solve(A, w.vector(), L) 
solver_lin = LUSolver("mumps")
solver_lin.solve(A, w.vector(), L)

# Split using deep copy
(u, p) = w.split(True)

print('finish solving linear problem')
# save the fluid solution
out_u = File('workspace/dolfin-adjoint/3-D_application/results/u.pvd')
out_u << u
out_p = File('workspace/dolfin-adjoint/3-D_application/results/p.pvd')
out_p << p

Thank you very much

2

Your domain’s boundary probably doesn’t lie exactly on the circle of radius 0.5625. You should consider a geometry comparison tolerance.

3

Thank you for the hint, I have exported the subdomains and boundaries and have found that some of them are not highlighted properly.

1 Like