Hi, I am trying to implement time dependent Dirichlet boundary conditions to a single point in the middle of the side of a square bar (in this case a displacement ramping up to some arbitrary amplitude).
I followed similar threads on the topic such as this one, but I couldn’t seem to apply it to my geometry. It also feels like there is a better way to define the BC rather than using the interpolate function, since it is only at a single point. See the code below:
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh, default_scalar_type, log
from dolfinx.io import XDMFFile
import ufl
import dolfinx.fem.petsc
import dolfinx.nls.petsc
import time
start_time = time.time()
# Create mesh and function space
L = 1.0 # Length of the rod
nx, ny, nz = 4, 4, 20 # Number of elements in x, y, and z directions
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array(
[0.1, 0.1, L])], [nx, ny, nz], cell_type=mesh.CellType.hexahedron)
V = fem.VectorFunctionSpace(domain, ("CG", 1))
fdim = domain.topology.dim - 1
#no penetration in z dir boundary conditions
def nopen_boundary(x):
return np.isclose(x[2], 0)
boundary_facets = mesh.locate_entities_boundary(domain, fdim, nopen_boundary)
bc = fem.dirichletbc(default_scalar_type(0), fem.locate_dofs_topological(V.sub(2), fdim, boundary_facets), V.sub(2))
#time dep boundary conditions
Amp=1e-3 #amplitude of perturbation
def point1(x):
return np.allclose(x, [0.05,0.0,0.0])
class time_dependent_dbc1:
def __init__(self):
self.t = 0.0
def eval(self, x):
#return np.stack((np.full(x.shape[0], Amp*self.t), (np.zeros(x.shape[1])),(np.zeros(x.shape[2])) ))
return np.array([Amp*self.t, 0.0,0.0])
V2 = V.sub(2).collapse()[0]
point1_pert = time_dependent_dbc1()
point1_pert.t = 0.
bc_fun = fem.Function(V2)
bc_fun.interpolate(point1_pert.eval)
point1_dofs = fem.locate_dofs_geometrical((V.sub(2),V2), point1)
bc1 = fem.dirichletbc(bc_fun, point1_dofs, V.sub(2))
allbcs=[bc, bc1]
# Elasticity constants
E = 1.0e6
nu = 0.3
mu = E / (2.0 * (1 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
# Time integration parameters
am = 0.2
af = 0.4
g = 0.5 + af - am
alpha_m = fem.Constant(domain, am)
alpha_f = fem.Constant(domain, af)
gamma = fem.Constant(domain, g)
beta = fem.Constant(domain, (g + 0.5)**2 / 4.0)
rho = 1.0
eta_m = 0.0
eta_k = 0.0
T = 1.0
Nsteps = 50
dt = fem.Constant(domain, T / Nsteps)
# Functions and fields
u_ = ufl.TestFunction(V)
u = fem.Function(V, name='Displacement')
u_old = fem.Function(V)
v_old = fem.Function(V)
a_old = fem.Function(V)
# External load
# B = fem.Constant(domain, PETSc.ScalarType((1.0, 0.0, 0.0))
# ) # Applied load in the x-direction
# Measures
dx = ufl.Measure("dx", domain=domain)
# GA method functions
def sigma(r):
return lmbda*ufl.nabla_div(r)*ufl.Identity(len(r)) + 2*mu*ufl.sym(ufl.grad(r))
def m(u, u_):
return rho*ufl.inner(u, u_)*dx
def k(u, u_):
return ufl.inner(sigma(u), ufl.grad(u_))*dx
def c(u, u_):
return eta_m*m(u, u_) + eta_k*k(u, u_)
# update variables
def update_a(u, u_old, v_old, a_old, ufl_=True):
if ufl_:
dt_ = dt
beta_ = beta
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
def update_v(a, u_old, v_old, a_old, ufl_=True):
if ufl_:
return v_old + dt*((1-gamma)*a_old + gamma*a)
else:
return v_old + float(dt) * ((1 - float(gamma)) * a_old + float(gamma)*a)
def update_fields(u, u_old, v_old, a_old):
u_vec, u0_vec = u.x.array[:], u_old.x.array[:]
v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl_=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl_=False)
v_old.x.array[:] = v_vec #numpy.ndarray stored inside dolfinx.fem.function.Function
a_old.x.array[:] = a_vec
u_old.x.array[:] = u_vec
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
a_new = update_a(u, u_old, v_old, a_old, ufl_=True) #ufl.algebra.Sum
v_new = update_v(a_new, u_old, v_old, a_old, ufl_=True)
## Problem setup
res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_) + \
k(avg(u_old, u, alpha_f), u_) #- rho * ufl.inner(B, u_)*dx # residual to solve
problem = fem.petsc.NonlinearProblem(res, u, bcs=allbcs) #petsc.linear problem?
## Solver setup
solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
# solver options
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
# Quantities to track
times = np.linspace(0, T, Nsteps + 1)
# Create file
with XDMFFile(domain.comm, "rod_dbc0.xdmf", "w") as xdmf_file:
xdmf_file.write_mesh(domain)
xdmf_file.write_function(u, 0)
# Solve time to assess time spent in each component
# solve_time = 0
# Write function at each time step
for i in range(Nsteps):
t = times[i + 1]
print("Time = %s" % t)
point1_pert.t=t
bc_fun.interpolate(point1_pert.eval)
problem = fem.petsc.NonlinearProblem(res, u, bcs=allbcs)
solver = dolfinx.nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
# log.set_log_level(log.LogLevel.INFO) # Convergence info
# s = time.perf_counter()
num_its, converged = solver.solve(u)
# e = time.perf_counter()
# solve_time += e-s
# print(solve_time)
assert converged
# print(f"Number of interations: {num_its:d}")
u.x.scatter_forward()
# Update old fields with new quantities
update_fields(u, u_old, v_old, a_old)
# Write the function to the XDMF file
xdmf_file.write_function(u, t)
print("--- %s seconds runtime ---" % (time.time() - start_time))
Note: the code runs fine if the time dependent BCs are taken out.
Thanks!