Hello everyone,
I have a problem with changing the pressure on a surface through time steps. The pressure rises and then falls to null afterwards. I suppose the object must return to its original configuration, right? But it just keeps getting more and more distorted. And when I decrease the time step (dt) and increasing the number of steps accordingly, the distortion gets more severe. It’s like it is adding up the pressure of the new step to that of the previous one. Here is a code demonstrating what I mean:
import dolfinx as fe
from dolfinx.fem.petsc import NonlinearProblem, LinearProblem, assemble_vector, assemble_matrix, create_vector, apply_lifting, set_bc
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import os
from mpi4py import MPI
import ufl
from petsc4py import PETSc
import pyvista
def main(linear = True, time_dependent = True):
L = 1
W = 0.2
mesh = fe.mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
[20, 6, 6], cell_type=fe.mesh.CellType.hexahedron)
element = ufl.VectorElement("Lagrange", mesh.ufl_cell(), 1)
V = fe.fem.functionspace(mesh, element)
dt = fe.default_scalar_type(0.01)
num_steps = int(1/float(dt))
maximum_pressure = 15
minimum_pressure = 0
pressure = np.append(np.linspace(minimum_pressure, maximum_pressure, num=int(num_steps/2)),\
np.linspace(maximum_pressure, minimum_pressure, num=int(num_steps/2)))
print(pressure)
p = fe.fem.Constant(mesh, pressure[0])
def clamped_boundary(x):
return np.isclose(x[0], 0)
fdim = mesh.topology.dim - 1
clamped_boundary_facets = fe.mesh.locate_entities_boundary(mesh, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=fe.default_scalar_type)
bc = fe.fem.dirichletbc(u_D, fe.fem.locate_dofs_topological(V, fdim, clamped_boundary_facets), V)
boundaries = [(1, lambda x: np.isclose(x[0], L))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = fe.mesh.locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = fe.mesh.meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
######################################################################################################
#########################################Formulation##################################################
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
old_u = fe.fem.Function(V)
old_velocity = fe.fem.Function(V)
old_acceleration = fe.fem.Function(V)
d = len(u)
I = ufl.variable(ufl.Identity(d)) # Identity tensor
F = ufl.variable(I + ufl.grad(u)) # Deformation gradient
C = ufl.variable(F.T*F) # Right Cauchy-Green tensor
J = ufl.variable(ufl.det(F))
####Check out for the metadata
metadata = {"quadrature_degree": 4}
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=mesh, metadata=metadata)
# Generalized-alpha method parameters
alpha_m = fe.fem.Constant(mesh, 0.2)
alpha_f = fe.fem.Constant(mesh, 0.4)
gamma = fe.default_scalar_type(0.5+alpha_f-alpha_m)
beta = fe.default_scalar_type((gamma+0.5)**2/4.)
# Update formula for acceleration
# a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(u_new, u_old, v_old, a_old):
"""Update fields at the end of each time step."""
# Get vectors (references)
u_vec, u0_vec = u_new.x.array[:], u_old.x.array[:]
v0_vec, a0_vec = v_old.x.array[:], a_old.x.array[:]
# use update functions using vector arguments
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)
# Update (u_old <- u)
v_old.x.array[:], a_old.x.array[:] = v_vec, a_vec
u_old.x.array[:] = u_new.x.array[:]
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
normal = -ufl.FacetNormal(mesh)
rho = 1
E = fe.default_scalar_type(1.0e4)
nu = fe.default_scalar_type(0.3)
mu = fe.fem.Constant(mesh, E / (2 * (1 + nu)))
lmbda = fe.fem.Constant(mesh, E * nu / ((1 + nu) * (1 - 2 * nu)))
def S(u):
return 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
acceleration = update_a(u, old_u, old_velocity, old_acceleration, ufl=False)
velocity = update_v(acceleration, old_u, old_velocity, old_acceleration, ufl=True)
formulation = rho * ufl.inner(alpha_m*old_acceleration + (1-alpha_m)* acceleration, v) * dx \
+ ufl.inner(ufl.grad(v), S(avg(old_u, u, alpha_f))) * dx \
- ufl.inner(v, p * normal) * ds(1)
bilinear_form = fe.fem.form(ufl.lhs(formulation))
linear_form = fe.fem.form(ufl.rhs(formulation))
######################################################################################################
######################################################################################################
######################################################################################################
###############################################Solving################################################
A = assemble_matrix(bilinear_form, bcs=[bc])
A.assemble()
b = create_vector(linear_form)
solver = PETSc.KSP().create(mesh.comm)
solver.setInitialGuessNonzero(True)
solver.setOperators(A)
solver.getPC().setType(PETSc.PC.Type.SOR)
solver.view()
displacement = fe.fem.Function(V)
##############Initializing the plot###################
pyvista.start_xvfb()
plotter = pyvista.Plotter()
plotter.open_gif(os.path.join("deformation.gif"), fps=3)
topology, cell_types, geometry = fe.plot.vtk_mesh(mesh, 3)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = displacement.x.array.astype(float).reshape(geometry.shape[0], len(displacement))
grid.point_data["u"] = values
grid.set_active_vectors("u")
# Warp mesh by deformation
warped = grid.warp_by_vector("u", factor=10)
warped.set_active_vectors("u")
# viridis = mpl.colormaps.get_cmap("viridis").resampled(25)
# sargs = dict(title_font_size=25, label_font_size=20, fmt="%.2e", color="black",
# position_x=0.1, position_y=0.8, width=0.8, height=0.1)
# print(max(displacement.x.array))
renderer = plotter.add_mesh(warped, show_edges=True, lighting=False,
# cmap=viridis, scalar_bar_args=sargs,
clim=[0, 1])
Vs = fe.fem.FunctionSpace(mesh, ("Lagrange", 1))
magnitude = fe.fem.Function(Vs)
us = fe.fem.Expression(ufl.sqrt(sum([displacement[i]**2 for i in range(len(displacement))])), Vs.element.interpolation_points())
magnitude.interpolate(us)
warped["mag"] = magnitude.x.array
plotter.write_frame()
for i in range(num_steps):
p.value = pressure[i]
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, linear_form)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [bilinear_form], [[bc]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [bc])
# Solve linear problem
solver.solve(b, displacement.vector)
displacement.x.scatter_forward()
# Update old fields with new quantities
update_fields(displacement, old_u, old_velocity, old_acceleration)
warped["u"][:, :len(displacement)] = displacement.x.array.reshape(geometry.shape[0], len(displacement))
magnitude.interpolate(us)
warped.set_active_scalars("mag")
warped_n = warped.warp_by_vector(factor=10)
plotter.update_coordinates(warped_n.points.copy(), render=False)
plotter.update_scalar_bar_range([0, 10])
plotter.update_scalars(magnitude.x.array)
plotter.write_frame()
plotter.close()
if __name__ == "__main__":
main()