# Brief doubt with plotting fenics solution

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

results = np.zeros((200, 3))

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.xlabel("Time")
plt.ylabel("Average Values")
plt.legend()
plt.show()

plt.figure()
plt.xlabel("Time")
plt.ylabel("Average Values")
plt.legend()
plt.show()

# Plotting individual energy terms
plt.figure()
plt.xlabel("Time")
plt.ylabel("Elastic Energy")
plt.legend()
plt.show()

plt.figure()
plt.xlabel("Time")
plt.ylabel("Fracture Term Energy")
plt.legend()
plt.show()

plt.figure()
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?

Could you please state what you expect to get, and show what plots you are currently getting?
As far as I can tell from looking at your code F1 and F2 are residuals. These can be associated with the finite element space of the test function u1 and v1, and plotted over space for each time step.

Is that what you want? or do you want a single accumulated energy value for each time step?

Hello, sorry for this. what I’m getting is:

There is a lot of lines and I don’t know why, it should be the single accumulated energy at each time step, That is, the total elastic, fracture and total energy at each time step.

Also, it seems like that is not adding both energy terms the total energy

Start by considering these variables. What do you see when you print them? I would expect these to be a vector, as F1 is a residual formulation.

This means that you have to perform the reduction of this vector to a given number.

Ok, I’ve been thinking for a while and I don’t know if it’s correct, but with your idea, I’ve got:

I summed all values and, as a good news, I’ve got only positive terms:
Cumulative Elastic Energy: 0.011589729018203784
Cumulative Fracture Energy: 2.1286995202790905e-15
Cumulative Total Energy: 0.011589729018205912
What I did:

``````    elastic_energy = assemble(F1)
fracture_energy = assemble(F2)
elastic_energy_values.append(elastic_energy.get_local().sum())
fracture_energy_values.append(fracture_energy.get_local().sum())

total_energy = assemble(F1) + assemble(F2)
total_energy_values.append(elastic_energy.get_local().sum() + fracture_energy.get_local().sum())

cumulative_elastic_energy += elastic_energy.get_local().sum()
cumulative_fracture_energy += fracture_energy.get_local().sum()
cumulative_total_energy = cumulative_elastic_energy + cumulative_fracture_energy
``````

And, two questions, these lines:

``````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
``````

Are really restricting the values of v to be between 0 and 1? Because there is some negative values for v;
Where can I found more about get_local() functions? They are really good to extract numerical values;
Also, thanks.

1 Like