How to implement Dirichlet AND Robin BC at the same facet?

variational formulation
1

Dear community,

I try to solve a transient Poisson equation with Dirichlet BCs at every facet of my 1x1 UnitSquareMesh domain and on top of that, one Robin BC for the left facet.

My MWCE is running, but I expect a different system behavior, when I post process the solution in ParaView, in fact a much higher reduction of u over the domain due to the convective nature of the Robin BC.

So, I´m not sure, if I implemented the Robin AND Dirichlet BC at the left boundary the right way.

I took the following steps to implement the problem:

1. Defined u for t=0 by an expression (Not really important for the implementation of the BC, but probably better for your understanding of the MWCE) 
# Define u for t=0
t = 0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        tol = 1E-14
        
        # u on left boundary
        if x[0] < tol:
            u = Constant(30)
            value[0] = u
            
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(15)
            value[0] = u
            
        # u for the rest of the domain and boundaries
        else:
            u = 40
            value[0] = u
        return value[0]
    def value_shape(self):
        return ()
      
u = MyExpression_t0()
´´´
  1. Implement Dirichlet for every facet
# Define Dirichlet boundary conditions 
dirichlet_bcs = []      

if t>0:
    # Calculate average of u of the left facet 
    # (not that important for demonstration purposes, but part of the project I´m currently working on)
    u_avg_Dirichlet_left = assemble(u*ds(1)) / assemble(Constant(1.0)*ds(1))
    
    # Calculate timedependent left Dirichlet boundary condition
    u_Dirichlet_left = Constant1 * (Constant2 - u_avg_Dirichlet_left)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_left, boundaries, 1))
    
    # Constant top Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_top = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_top, boundaries, 2))
    
    # Constant right Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_right = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_right, boundaries, 3))
        
    # Constant bottom Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_bottom = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_bottom, boundaries, 4))
else:
    # Left Dirichlet boundary condition for t=0
    u_Dirichlet_left = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_left, boundaries, 1))
    
    # Top Dirichlet boundary condition for t=0
    u_Dirichlet_top = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_top, boundaries, 2))
    
    # Right Dirichlet boundary condition for t=0
    u_Dirichlet_right = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_right, boundaries, 3))
    
    # Bottom Dirichlet boundary condition for t=0
    u_Dirichlet_bottom = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_bottom, boundaries, 4))
  1. Define Robin BC
# Define Robin boundary condition -du/dn = p*(u-q)
p = Constant(10000.0)
q = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
  1. Define variational problem and specify the Robin BC by ds(1)
a = u*v*dx + dt*inner(grad(u), grad(v))*dx + dt*p*u*v*ds(1)
L = (u_n + dt*f)*v*dx + dt*p*q*v*ds(1)
  1. Compute solution
solve(a == L, u, bcs=dirichlet_bcs)

I suspect, that either bcs=dirichlet_bcs of step 5 does not consider the Robin BC nor that the variational formulation is wrong.
Or, I underestimate the influence of the Robin BC in the way it is prescribed right now.

Thats´s the full MWCE: 
from fenics import *

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
Constant1 = 10
Constant2 = 4

# Create mesh and define function space
nx = ny = 8
mesh = UnitSquareMesh(nx, ny)
space_dim = mesh.geometry().dim()
V = FunctionSpace(mesh, 'P', 1)

vtkfile = File("20200323_ft03_heat_Robin_AND_Dirichlet/solution.pvd")

# Create classes for defining parts of the boundaries of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)
        
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
    
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, space_dim-1)
#boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define new measures associated with the exterior boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

# Define u for t=0
t = 0
class MyExpression_t0(UserExpression):
    def eval(self, value, x):
        tol = 1E-14
        
        # u on left boundary
        if x[0] < tol:
            u = Constant(30)
            value[0] = u
            
        # u on left domain
        elif x[0] > 0 and x[0] <= 0.1:
            u = Constant(15)
            value[0] = u
            
        # u for the rest of the domain and boundaries
        else:
            u = 40
            value[0] = u
        return value[0]
    def value_shape(self):
        return ()
      
u = MyExpression_t0()

# Define Dirichlet boundary conditions 
dirichlet_bcs = []      

if t>0:
    # Calculate average of u of the left facet 
    # (not that important for demonstration purposes, but part of the project I´m currently working on)
    u_avg_Dirichlet_left = assemble(u*ds(1)) / assemble(Constant(1.0)*ds(1))
    
    # Calculate timedependent left Dirichlet boundary condition
    u_Dirichlet_left = Constant1 * (Constant2 - u_avg_Dirichlet_left)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_left, boundaries, 1))
    
    # Constant top Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_top = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_top, boundaries, 2))
    
    # Constant right Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_right = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_right, boundaries, 3))
        
    # Constant bottom Dirichlet boundary condition to reduce complexity of the MWCE
    u_Dirichlet_bottom = Constant(5)
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_bottom, boundaries, 4))
else:
    # Left Dirichlet boundary condition for t=0
    u_Dirichlet_left = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_left, boundaries, 1))
    
    # Top Dirichlet boundary condition for t=0
    u_Dirichlet_top = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_top, boundaries, 2))
    
    # Right Dirichlet boundary condition for t=0
    u_Dirichlet_right = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_right, boundaries, 3))
    
    # Bottom Dirichlet boundary condition for t=0
    u_Dirichlet_bottom = u
    dirichlet_bcs.append(DirichletBC(V, u_Dirichlet_bottom, boundaries, 4))
    
# Define Robin boundary condition -du/dn = p*(u-q)
p = Constant(10000.0)
q = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

# Define initial value
u_n = interpolate(u, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(40)

a = u*v*dx + dt*inner(grad(u), grad(v))*dx + dt*p*u*v*ds(1)
L = (u_n + dt*f)*v*dx + dt*p*q*v*ds(1)

# Time-stepping
u = Function(V)
for n in range(num_steps):

    # Update current time
    t += dt

    # Compute solution
    solve(a == L, u, bcs=dirichlet_bcs)
    vtkfile << (u, t)
    
    # Update previous solution
    u_n.assign(u)

Thanks a lot and stay healthy everyone! :+1: :+1:

2

Could you please mathematically describe what boundary condition you would like to have on the boundary in question. If you would like to have a the following condition:
\frac{\partial u}{\partial n} + \alpha u = g and u=h on the same boundary, this is equivalent to enforcing a non-zero Neumann condition: \frac{\partial u}{\partial n} = g - \alpha h. In general, you cannot enforce Dirichlet conditions in combination with other conditions such as Neumann conditions, see 1D Cantilever Beam problem.

1 Like
3

I try to model a drug release process.

Equation (1) describes the decreasing amount of drug, that is released from the left boundary
u = f(t)*(u_{average\;of\;left\;facet}-Consant2),\;\;\;\;\;\;\;\;\;(1)
where f(t) is a function, that describes the reducing size of drug particles.
Equation (2) decribes the interaction between the domain of the drug and the surrounding domain by
\frac{du}{dn} = - p*(u-q) , \;\;\;\;\;\;\;\;\;(2)
where p and q (in contrast to the MWCE, where q is a function, because I used the Robin BC implementation from a tutorial) are constants.

Is it legit to describe this boundary behavior as a non-zero Neumann condition by
\frac{du}{dn} = - p*(f(t)*(u_{average\;of\;left\;facet}-Consant2)-q)? \;\;\;\;\;\;\;\;\;(3)