This is the code, most of it is the code that was given for an adaptation to the elastodynamics code to Fenicsx, the changes are regarding the mesh, and the time loop, and the final plotting
gmsh.initialize()
Outer_cylinder = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, L, R_outer)
Inner_cylinder = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, L, R_inner)
gmsh.model.occ.synchronize()
tube = gmsh.model.occ.cut([(3, Outer_cylinder)], [(3, Inner_cylinder)])
tube_tag = tube[0][0][1]
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(3, [tube_tag], tag=1)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin",0.01)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax",0.08)
gmsh.model.mesh.generate(3)
gmsh.model.mesh.optimize("Netgen")
domain, ct, _ = model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim=3)
gmsh.finalize()
#### Space definition and function definition
V = fem.functionspace(domain, ("CG", 1, (domain.geometry.dim,)))
Vsig = fem.functionspace(domain, ("DG", 0, (domain.geometry.dim, domain.geometry.dim)))
def boundary_condition(x):
return np.isclose(x[2], 0, atol=1e-1)
fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, boundary_condition)
print("Boundary facets:", boundary_facets)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc_end_0 = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)
bcs = [bc_end_0]
boundary_tag = 1
facet_indices = np.array(boundary_facets, dtype=np.int32)
facet_values = np.full(len(facet_indices), boundary_tag, dtype=np.int32)
boundary_mt = mesh.meshtags(domain, fdim, facet_indices, facet_values)
i=1
# 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.1
eta_m = 0.0
eta_k = 0.0
T = 10
Nsteps = 20
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)
# Force expression
Traction = fem.Constant(domain, default_scalar_type([0,0,0]))
ds = ufl.Measure("ds", domain=domain)
dx = ufl.Measure("dx", domain=domain)
# Stress tensor
def F(u):
return ufl.Identity(len(u)) + ufl.nabla_grad(u)
def F_growth(t, alpha):
return ufl.as_matrix([
[1, 0, 0],
[0, 1, 0],
[0, 0, 1+alpha*t]
])
def epsilon(u, t, alpha):
F_stress = F(u)
F_g = F_growth(t, alpha)
F_stres = ufl.dot(F_stress, ufl.inv(F_g))
return 0.5 * (ufl.dot(F_stress.T,F_stress) - ufl.Identity(len(u)))
def S(E1, E3, nu12, nu13, nu31, G13):
S = np.array([
[1/E1, -nu12/E1, -nu31/E3, 0, 0, 0],
[-nu12/E1, 1/E1, -nu31/E3, 0, 0, 0],
[-nu13/E1, -nu13/E1, 1/E3, 0, 0, 0],
[0, 0, 0, 1/(2*G13), 0, 0],
[0, 0, 0, 0, 1/(2*G13), 0],
[0, 0, 0, 0, 0, (1+nu12)/E1]
]
)
return S
def sigma(u, t, alpha):
E1 = E_list[1][0]
E3 = E_list[1][1]
nu12 = nu_list[1][0]
nu13 = nu_list[1][1]
nu31 = nu_list[1][2]
G13 = G_list[1]
matrix_S = S(E1, E3, nu12, nu13, nu31, G13)
inverse_S = ufl.as_matrix(np.linalg.inv(matrix_S))
eps = epsilon(u, t, alpha)
vector_eps = ufl.as_vector([eps[0, 0], eps[1, 1], eps[2, 2], eps[1, 2], eps[0, 1], eps[0, 2]])
vector_sigma = ufl.dot(inverse_S, vector_eps)
sigma_matrix = ufl.as_matrix([
[vector_sigma[0], vector_sigma[4], vector_sigma[5]],
[vector_sigma[4], vector_sigma[1], vector_sigma[3]],
[vector_sigma[5], vector_sigma[3], vector_sigma[2]]
])
return sigma_matrix
## Below this is the elastodynamics file code (copied from another post on here)
def m(u, u_):
return rho*ufl.inner(u, u_)*dx
def k(u, u_, t, alpha):
return ufl.inner(sigma(u, t, alpha), ufl.grad(u_)) * dx
def c(u, u_, t, alpha):
return eta_m * m(u, u_) + eta_k * k(u, u_, t, alpha)
#Update fields
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_ = dt
beta_ = beta.value
return (u-u_old-float(dt_)*v_old)/float(beta_)/float(dt_)**2 - (1-2*float(beta_))/2/float(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
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)
v_new = update_v(a_new, u_old, v_old, a_old, ufl_=True)
# Quantities to track
times = np.linspace(0, T, Nsteps + 1)
u_tip = np.zeros((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)
plotter = pyvista.Plotter()
plotter.open_gif("growth_trial.gif", fps=3)
topology, cells, geometry = plot.vtk_mesh(u.function_space)
function_grid = pyvista.UnstructuredGrid(topology, cells, geometry)
values = np.zeros((geometry.shape[0], 3))
values[:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
function_grid["u"] = values
function_grid.set_active_vectors("u")
warped = function_grid.warp_by_vector("u", factor=1)
warped.set_active_vectors("u")
# Add mesh to plotter and visualize
actor = plotter.add_mesh(warped, show_edges=True, lighting=False, clim=[0, 10], cmap = "viridis")
solve_time = 0
# Write function at each time step
for i in range(Nsteps):
t = times[i + 1]
print("Time = %s" % t)
# Problem setup
res = m(avg(a_old, a_new, alpha_m), u_) + c(avg(v_old, v_new, alpha_f), u_, t, alpha) + \
k(avg(u_old, u, alpha_f), u_, t, alpha) - rho * \
ufl.inner(Traction, u_)*ds(2) # residual to solve
problem = fem.petsc.NonlinearProblem(res, u, bcs)
# Solver setup
solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
s = time.perf_counter()
num_its, converged = solver.solve(u)
e = time.perf_counter()
solve_time += e-s
assert converged
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)
#normal = nanson_normal(u,normal)
# Update function_grid with new solution
function_grid["u"][:, :len(u)] = u.x.array.reshape(geometry.shape[0], len(u))
warped_n = function_grid.warp_by_vector(factor=1)
warped.points[:, :] = warped_n.points
plotter.update_scalar_bar_range([0, 10])
plotter.show_axes()
plotter.write_frame()
plotter.close()
print("--- %s seconds runtime ---" % (time.time() - start_time))
print(solve_time)