Hi there! I am experiencing problem defining BCs on submeshes extracted from a parent mesh.
My geometry consists of two concentric circles. I am solving one PDE for the smaller disk and another for the annulus. I successfully recovered the submeshes but having trouble defining the boundaries to my submeshes. Below is the gmsh script followed by the Python script.
gmsh script:
//+
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 5, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 3, 0, 2*Pi};
//+
Curve Loop(1) = {2};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {1};
//+
Curve Loop(3) = {2};
//+
Plane Surface(2) = {2, 3};
//+
Physical Curve(1) = {2};
//+
Physical Curve(2) = {1};
//+
Physical Surface(5) = {1};
//+
Physical Surface(4) = {2};
Python script:
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
from dolfin import *
import meshio
from dolfin import Mesh, XDMFFile, File, MeshValueCollection, cpp
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
"eliminate_zeros": True, \
"precompute_basis_const": True, \
"precompute_ip_const": True}
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
import meshio
msh = meshio.read("example.msh")
triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh)
from dolfin import *
mesh = Mesh()
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
ds_custom = Measure("ds", domain=mesh, subdomain_data=mf)
s5 = MeshView.create(cf, 5) # Extracting subspace 5
s4 = MeshView.create(cf, 4) # Extracting subspace 4
W1 = FunctionSpace(s5,"Lagrange",1) # Defining function space on s5
W2 = FunctionSpace(s4,"Lagrange",1) # Defining function space on s4
def bs5(x): # Defining boundary for subspace s5, curve 1
return mf == 1
def b2s4(x): # Defining boundary for subspace s4, curve 2
return mf == 2
# PROBLEM FOR DOMAIN 5 WITH BOUNDARY 1
# Define boundary condition
G , mu = 1, 0.1
u_D=Constant(0.0)
bc = DirichletBC(W1, u_D, bs5)
# Define variational problem
u = TrialFunction(W1)
v = TestFunction(W1)
f = Constant(-G/mu) # f=-G/mu
a = dot(grad(u), grad(v))*dx
L = -f*v*dx
# Compute solution
us = Function(W1)
solve(a == L, us, bc)
# Save solution to file in PVD format for Paraview
file = File("examplelumenvelocity.pvd");
file << us;
# PROBLEM FOR SURFACE 4 WITH BOUNDARIES 1 & 2
D, k, t = 0.7, 0.3, 3
f0, f1, beta = 5, 2, 0.3
gamma1 = (-1/D)*Constant(f0 + f1*(1-exp(-beta*(t**2)))) # Neumann condition on inner curve -dot(grad(c),n)=gamma1
gamma2 = Constant(0.0) # Dirichlet condition on outer circle
bc2 = DirichletBC(W2, gamma2, b2s4)
bcs=[bc2]
# Define variational problem
w = TestFunction(W2)
c = TrialFunction(W2)
a = inner(grad(c), grad(w))*dx + (((k/D)*c)*w)*dx
a -= - inner(gamma1, w)*ds_custom(bs5)
# Compute solution
cs= Function(W2)
solve(a == 0, cs, bcs)
# Save solution to file in PVD format for Paraview
file = File("example.pvd");
file << cs;
Error message:
*** Warning: Found no facets matching domain for boundary condition.
Invalid subdomain_id <function bs5 at 0x7fc435e3f1f0>.
How to fix this? I am using Dolfin version 2019.2.0.dev0.
