Hello there. Actually, i’m ending to solve my problem, but i’m founding some issues to plot the solutions:
from __future__ import print_function
from petsc4py import PETSc
from fenics import *
import matplotlib.pyplot as plt
from IPython.display import clear_output
import numpy as np
import time
from dolfin import *
combined_file = XDMFFile("combined_solution.xdmf")
combined_file.parameters["flush_output"] = True # Ensure data is written immediately
vtkfile1 = File('Elastic_term.pvd')
vtkfile2 = File('Fracture_term.pvd')
vtkfile3 = File('Total_term.pvd')
T = 25.0 # final time
num_steps = 500 # number of time steps
dt = T / num_steps # time step size
tol = 1E-14; eta = 1*10E-7; eps = 4*10E-2; delta = 0.5; Alternate_tol = 1*10E-4; max_iter = 100
G_c = 1
mesh = UnitSquareMesh(5,5)
V = FunctionSpace(mesh, "Lagrange", 1)
g_t = Expression('t', degree=1, t=0)
def boundary_D(x, on_boundary):
return on_boundary and (near(x[0], 0, tol) or near(x[0], 5, tol))
bc = DirichletBC(V, g_t, boundary_D)
u = Function(V)
v = Function(V)
u1 = TestFunction(V)
v1 = TestFunction(V)
v_previous = Function(V)
# Enforce values between 0 and 1
v_projected = project(v, V) # Project v to ensure it lies in V
v_constrained = Function(V)
v_constrained.assign(v_projected)
apply_constraints = [v_constrained.vector()[i] if v_constrained.vector()[i] <= 1.0 else 1.0 for i in range(len(v_constrained.vector()))]
apply_constraints = [apply_constraints[i] if apply_constraints[i] >= 0.0 else 0.0 for i in range(len(apply_constraints))]
v_constrained.vector()[:] = apply_constraints
# Update the Function object with constrained values
v_previous.assign(v_constrained)
v_previous.interpolate(Constant(1.0)) # Initialize v0 as 1.0
V_scalar = FunctionSpace(mesh, 'CG', 1) # Choose appropriate function space
elastic_energy_func = Function(V_scalar)
fracture_energy_func = Function(V_scalar)
total_energy_func = Function(V_scalar)
norm1 = []
avg_u_values = [] # List to store average u values at each time step
avg_v_values = [] # List to store average u values at each time step
elastic_energy_values = [] # List to store elastic energy at each time step
fracture_energy_values = [] # List to store fracture term energy at each time step
total_energy_values = [] # List to store total energy at each time step
F1 = (v_previous**2 + eta)*u*u1*dx
F2 = ((grad(u)**2)*v*v1 + 2*G_c*(v-1)*v1/(4*eps) + eps*inner(grad(v), grad(v1)))*dx
results = np.zeros((200, 3))
loading = np.linspace(0, 25, num_steps)
for t in range(num_steps):
# Update the projection of u at the beginning of each time step
check_tol = False
i = 0
maxnorm = 1
while maxnorm > Alternate_tol:
solve(F1 == 0, u, bc)
# Solve for v using the previous time step's solution as an upper bound
solve(F2 == 0, v)
u_vec = as_backend_type(u.vector()).vec()
v_vec = as_backend_type(v.vector()).vec()
v_previous_vec = as_backend_type(v_previous.vector()).vec()
diff_vec = v_vec.copy()
diff_vec.axpy(-1.0, v_previous_vec)
maxnorm = diff_vec.norm(PETSc.NormType.INFINITY) # Calculate the maximum norm
print('\n',f"Iteration {i}: MaxNorm(v) = {maxnorm}",'\n')
v_local = v.vector().get_local()
v_previous_local = v_previous.vector().get_local()
# condition_1 = np.any(v_local < v_previous_local)
# condition_2 = np.all((v_local >= 0) & (v_local <= 1))
v_previous.assign(v) # Update v_previous with the values of v
i += 1
if(i > max_iter):
break
# Update for the next time step
print(f"At time step {t}, maxnorm = {maxnorm}")
norm1.append(maxnorm)
avg_u = np.mean(u.vector().get_local()) # Calculate average value of u
avg_u_values.append(avg_u)
avg_v = np.mean(v.vector().get_local()) # Calculate average value of u
avg_v_values.append(avg_v)
elastic_energy = assemble(F1)
fracture_energy = assemble(F2)
elastic_energy_values.append(elastic_energy)
fracture_energy_values.append(fracture_energy)
total_energy = elastic_energy + fracture_energy
total_energy_values.append(total_energy)
# vtkfile1 << (cumulative_elastic_energy, t)
# vtkfile2 << (cumulative_fracture_energy, t)
# vtkfile3 << (cumulative_total_energy, t)
t += dt
g_t.t = t
combined_file.close()
plt.figure()
plt.plot(loading, avg_u_values, label='Average u')
plt.xlabel("Time")
plt.ylabel("Average Values")
plt.legend()
plt.show()
plt.figure()
plt.plot(loading, avg_v_values, label='Average v')
plt.xlabel("Time")
plt.ylabel("Average Values")
plt.legend()
plt.show()
# Plotting individual energy terms
plt.figure()
plt.plot(loading, elastic_energy_values)
plt.xlabel("Time")
plt.ylabel("Elastic Energy")
plt.legend()
plt.show()
plt.figure()
plt.plot(loading, fracture_energy_values)
plt.xlabel("Time")
plt.ylabel("Fracture Term Energy")
plt.legend()
plt.show()
plt.figure()
plt.plot(loading, total_energy_values)
plt.xlabel("Time")
plt.ylabel("Total Energy")
plt.legend()
plt.show()
print(norm1)
I want to plot the graphics depicting the solutions associated with the elastic term energy (F1) and fracture energy term (F2), as well as the total energy (F1 + F2) related to the solutions u and v. However, the current method I’m using isn’t producing the expected results. Could anyone provide tips or guidance on an alternative approach?
