I am wondering how I could apply Dirichlet boundary conditions if they come from a complicated function, for example something that interpolates some experimental random data. There is no explicit form for this function (cannot be written as a sin(t)). Instead, when we call it, it yields the desired value.
So, say we have
def complicated_func(t):
do some things like interpolate linearly random data between different points.
return desired_value_at_time_t
Then somewhere in the code we would do:
u_D = Expression('complicated_func(t)', complicated_func(t), t=0, degree = 2)
bc = DirichletBC(V, u_D, Top_region())
And in the time loop we would update the function like so:
u_D.t = t
But the problem is that it doesnât work. I havenât been able to make it work unless I explicitly write down the time varying expression, like sin(t) for example. Is there a way to make it work for a user-defined function like above? I am guessing yes, but I havenât found any working example so far.
Edit: After a discussion below, I add a MWE:
from dolfin import *
import numpy as np
import mshr
from scipy.interpolate import InterpolatedUnivariateSpline
def complicated_func(t):
return t
cylinder = mshr.Cylinder(Point(0, 0, 0), Point(0, 0, -9), 8.0, 8.0)
domain = cylinder
mesh = mshr.generate_mesh(domain, 20)
V = FunctionSpace(mesh, 'CG', 1)
# Define boundary conditions
tol = 1E-14
class Top(SubDomain):
def inside(self, x , on_boundary):
return on_boundary and near(x[2], 0.0, tol)
# Define subdomains.
# top
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[2] >= -6 - tol
# Other region.
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[2] <= -6 + tol
# Define the mesh
region = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
# Use the subdomains to mark the cells belonging to each subdomain
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(region, 0)
subdomain_1.mark(region, 1)
# k values of the layers
k_0 = 75
k_1 = 7
# Define the values of kappa in each subdomain
class K(UserExpression):
def __init__(self, region, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.region = region
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.region[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
kappa = K(region, k_0, k_1)
# Apply boundary conditions
u_D0 = Constant(0)
bc0 = [DirichletBC(V, u_D0, Top())]
# Define initial value
u_0 = Constant(5.0)
u_n = interpolate(u_0, V)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
F = u*v*dx + kappa*dot(grad(u), grad(v))*dx - u_n*v*dx
a, L = lhs(F), rhs(F)
u = Function(V)
num_steps= 10
t=0
dt=0.1
for n in range(num_steps):
# Update boundary condition (because it's time-dependent)
u_D0.assign = complicated_func(t)
#bc0 = DirichletBC(V, u_D0, Top())
#print('this line', complicated_func(t))
print('this line', u_D0.values(), t) # Debug.
# Update current time
t += dt
# Compute solution
solve(a == L, u, bcs = bc0)
#print(t, u.vector()[0,0,0] - u.vector()[0,0,-2.0])
#print('the value at the origin is ', u.vector()[0,0,0])
#print(u.vector()[:])
# Update previous solution
u_n.assign(u)