Hi everyone,
I am using docker containr with the image ceciledc/fenics_mixed_dimensional:latest as suggested by user cdaversin for solving a mixed dimensional problem.
I would like to solve a problem consisting of two different Poisson equations on different sinusoidal subdomains coupled by nonlinear boundary condition of Neumann type. Using the mixed_dimensional branch (based on the answer to these 2 questions Coupled Equations Over Multiple Domains , [Internal boundary condition with mixed_dimensional branch] (Internal boundary condition with mixed_dimensional branch), I tried to add the Dirichlet internal boundary condition to the sinusoidal interface, but this gives me the error. I want to set the boundary value = 1 at the top boundary, 0 at the interface and 0 at the bottom boundary.
Once Dirichlet works, I would like to give a non-linear Robin boundary condition like -kappa*gradient(ub) = exp(-ub) at the interface for the top domain and set ua=0 at the interface for the bottom domain.
Error:
Traceback (most recent call last):
File “poisson5.py”, line 135, in
bc_in1 = DirichletBC(W.sub_space(0), a_0, boundary_int, method = “pointwise”) #Bottom domain interface
File “/usr/local/lib/python3.6/dist-packages/fenics_dolfin-2019.2.0.dev0-py3.6-linux-x86_64.egg/dolfin/fem/dirichletbc.py”, line 131, in init
super().init(*args)
TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. dolfin.cpp.fem.DirichletBC(arg0: dolfin.cpp.fem.DirichletBC)
2. dolfin.cpp.fem.DirichletBC(V: dolfin.cpp.function.FunctionSpace, g: dolfin.cpp.function.GenericFunction, sub_domain: dolfin.cpp.mesh.SubDomain, method: str = ‘topological’, check_midpoint: bool = True)
3. dolfin.cpp.fem.DirichletBC(V: dolfin.cpp.function.FunctionSpace, g: dolfin.cpp.function.GenericFunction, sub_domains: dolfin.cpp.mesh.MeshFunctionSizet, sub_domain: int, method: str = ‘topological’)
Invoked with: <dolfin.cpp.function.FunctionSpace object at 0x7fdcbd78da40>, <dolfin.cpp.function.Constant object at 0x7fdcbd7a1500>, <dolfin.cpp.mesh.MeshFunctionSizet object at 0x7fdcc1954b58>, ‘pointwise’
Thank you!
Code:
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
from ufl import nabla_div
from ufl import sinh
import random
# Width of total domain
width = 100e-6
# Height of total domain
height = 200e-6
# Total resolution
resolution = 100
# Forcing function
f = Constant(0)
# Boundary values for the two equations
a_0 = 0 #Bottom boundary Dirichlet
b_0 = 1 #Top Boundary Dirichlet
i0 = 12.6
r = Constant(2*i0)
iApp = -1
F = 96485.33
R = 8.314
T = 300
c = Constant(F/(2*R*T))
kappa = 1
# g = iApp
# funcPhi0 = i0*(exp(C1*phi0) - exp(-C1*phi0))
# dFuncPhi0 = i0*C1*(exp(C1*phi0) + exp(-C1*phi0))
# r = i0*dFuncPhi0
# s = phi0 - funcPhi0/dFuncPhi0
g = iApp/kappa
tol = 1e-14
# Define SubDomain
class TopDomain(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.0001 + 0.00002*sin(5e4*x[0]) + tol
# Define Boundary
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], height)
class BottomBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0)
class Interface(SubDomain):
def inside(self, x, on_boundary):
return near (x[1], 0.0001 + 0.00002*sin(5e4*x[0]))
# # Define SubDomain
# class TopDomain(SubDomain):
# def inside(self, x, on_boundary):
# return x[1] >= 0.0001 + tol
# class BottomDomain(SubDomain):
# def inside(self, x, on_boundary):
# return x[1] <= 0.0001 - tol
# # Define Boundary
# class TopBoundary(SubDomain):
# def inside(self, x, on_boundary):
# return on_boundary and near(x[1], height)
# class BottomBoundary(SubDomain):
# def inside(self, x, on_boundary):
# return on_boundary and near(x[1], 0)
# class Interface(SubDomain):
# def inside(self, x, on_boundary):
# return near (x[1], 0.0001)
# Mesh for the entire domain
mesh_full = RectangleMesh(Point(0, 0),
Point(width, height),
resolution,
resolution)
# Define domains
domains = MeshFunction("size_t", mesh_full, 2)
domains.set_all(0)
top_domain = TopDomain()
top_domain.mark(domains, 1)
# out_sub0 = File("materials.pvd")
# out_sub0 << domains
# Create mesh for the left domain
mesh_top = MeshView.create(domains, 1)
mesh_bottom = MeshView.create(domains, 0)
# Define boundaries
top_boundary = TopBoundary()
bottom_boundary = BottomBoundary()
boundary_top = MeshFunction("size_t", mesh_top, 1)
top_boundary.mark(boundary_top, 1)
boundary_bottom = MeshFunction("size_t", mesh_bottom, 1)
bottom_boundary.mark(boundary_bottom, 1)
interf = Interface()
boundary_int = MeshFunction("size_t", mesh_full, 1)
boundary_int.set_all(0)
interf.mark(boundary_int, 1)
mesh_int = MeshView.create(boundary_int, 1)
# Create the function spaces
W_a = FunctionSpace(mesh_bottom, "Lagrange", 1)
W_b = FunctionSpace(mesh_top, "Lagrange", 1)
W = MixedFunctionSpace(W_a, W_b)
#dS = Measure("dS", domain=mesh_full, subdomain_data=boundary_int)
# ds_int = ds_bottom(1)
ds_int = Measure("dx", domain=mesh_int)
ds_top = Measure("ds", domain = mesh_top, subdomain_data = boundary_top)
# ds_bottom = Measure("ds", domain=mesh_bottom, subdomain_data=boundary_int)
# Define boundary conditions
bc_a = DirichletBC(W.sub_space(0), a_0, boundary_bottom, 1) #Bottom domain
bc_b = DirichletBC(W.sub_space(1), b_0, boundary_top, 1) #Top domain
bc_in1 = DirichletBC(W.sub_space(0), a_0, boundary_int, method = "pointwise") #Bottom domain interface
bc_in2 = DirichletBC(W.sub_space(1), a_0, boundary_int, method = "pointwise") #Top domain interface
# bc_c = DirichletBC(W.sub_space(0), a_0, interf, 1) #Interface for Bottom domain
# bc_d = DirichletBC(W.sub_space(1), a_0, interf, 1) #Interface for Top domain
bcs = [bc_a, bc_b, bc_in1, bc_in2] #for Dirichlet condition
# bcs = [bc_a] #For Neumann condition at top boundary
# Create all of the weak form variables
u = Function(W)
ua, ub = u.sub(0), u.sub(1)
va, vb = TestFunctions(W)
# ub.interpolate(Constant(0.001))
# ua.interpolate(Constant(0.0))
# Need to manually specify these.
dxa = Measure("dx", domain=W.sub_space(0).mesh())
dxb = Measure("dx", domain=W.sub_space(1).mesh())
# g = Expression("1.0", degree=2)
# c = Constant("10.0")
# The weak formulation
F_inner = kappa*( inner(grad(ua), grad(va)) * dxa + kappa*inner(grad(ub), grad(vb)) * dxb
- (f * va * dxa + f * vb * dxb )
)
# F_inner = kappa*( inner(grad(ua), grad(va)) * dxa + kappa*inner(grad(ub), grad(vb)) * dxb
# - (f * va * dxa + f * vb * dxb ) + g*vb*ds_top + r*sinh(C*ub)*vb*ds_int
# )
# F_inner = kappa*( inner(grad(ua), grad(va)) * dxa + kappa*inner(grad(ub), grad(vb)) * dxb
# - (f * va * dxa + f * vb * dxb ) + r*ub*vb*ds_int - r*s*vb*ds_int + g*vb*ds_top
# )
#Neumann boundary condition
# F_inner = ( inner(grad(ua), grad(va)) * dxa + inner(grad(ub), grad(vb)) * dxb
# - (f * va * dxa + f * vb * dxb ) + r*ub*vb*ds_int - r*s*vb*ds_int)
# F_coupling = va*c*ds_int + vb*c*ds_int
F_coupling = 0.0
F = F_inner + F_coupling
# solve(F==0, u, bcs, solver_parameters= {"nonlinear_solver": "newton",
# "newton_solver": {"linear_solver": "gmres"}})
Fi = extract_blocks(F)
dFi = []
for J in Fi:
for wi in u.split():
df = derivative(J, wi)
dFi.append(df)
u_sol = u.split()
problem = MixedNonlinearVariationalProblem(Fi, u_sol, bcs, dFi)
solver = MixedNonlinearVariationalSolver(problem)
solver.solve()
# J = derivative(F, u)
# problem = NonlinearVariationalProblem(F, u, bcs, J)
# solver = NonlinearVariationalSolver(problem)
# solver.solve()
# ## Export last solution
out_sub0 = File("ua.pvd")
out_sub1 = File("ub.pvd")
# out_exact = File("poisson-lm-2D-exact.pvd")
out_sub0 << u.sub(0)
out_sub1 << u.sub(1)
