### Unable to create Dirichlet boundary condition.

118
views
0
3 months ago by
Dear everyone,

I am new to FEniCS and I am trying to solve a simple PDE system with mixed Dirichlet-Robin condition.
I only would like to ask how to fixed the problem I am encountering when running the following program:

******************************************
from __future__ import print_function
from fenics import *

# Create mesh and define function space
from mshr import *
domain = Circle(Point(0, 0), 1.0) - Circle(Point(0, 0), 0.5)
mesh = generate_mesh(domain, 64)
V = FunctionSpace(mesh, 'P', 1)

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim())

# Define boundaries
class BoundaryGammaM(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and x[0]**2. + x[1]**2. < (1.0 + tol)**2. and x[1] >= 0.0
upper_boundary = BoundaryGammaM()
upper_boundary.mark(boundary_markers, 0)

class BoundaryGammaP(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and x[0]**2. + x[1]**2. < (1.0 + tol)**2. and x[1] <= 0.0
lower_boundary = BoundaryGammaP()
lower_boundary.mark(boundary_markers, 1)

class BoundaryGammaC(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and x[0]**2. + x[1]**2. < (0.5 + tol)**2.
inner_boundary = BoundaryGammaC()
inner_boundary.mark(boundary_markers, 2)

# Define parameters
k = Constant(1.0)
u_p = Constant(40.0)
h_m = Constant(1.1)
h_c = Constant(0.7)
u_c = Constant(22.0)
u_inf = Constant(27.0)
Q_m = Constant(0.0)
Q_c = Constant(0.0)

# Define boundary conditions
boundary_conditions = {0: {'Robin': (h_m, u_inf)},
1: {'Dirichlet': u_p},
2: {'Robin': (h_c, u_c)}}

# Redefine infinitesimals on the boundary
ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

# Collect Dirichlet conditions
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
bcs.append(bc)

# Define state variable and test function
u = TrialFunction(V)
v = TestFunction(V)

# Define bilinear and linear form
L = (h_m + Q_m)*u_inf*ds(0) + (h_c + Q_c)*u_c*ds(2)

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

******************************************
Here's what I get after running the program:

*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  User MeshFunction is not a facet MeshFunction (dimension is wrong).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
Community: FEniCS Project

1
3 months ago by

It can be fixed with a lower topological dimension of the MeshFunction

boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
Thanks for your suggestion! I tried it but now I am getting this error:

*** Error:   Unable to define linear variational problem a(u, v) = L(v) for all v.
*** Reason:  Expecting the right-hand side to be a linear form (not rank 0).
*** Where:   This error was encountered inside LinearVariationalProblem.cpp.
written 3 months ago by jf
nvm, found my mistake!
written 3 months ago by jf