hi all,
not sure if its ok to ask help with this, but i will try
basically im trying to create a graph of the sigma xx (y axis) and the deformation gradient (f11) and im facing some strange results, im receiving a straight line instead of the parabolic line as expected.
i have noticed that i have defined my equation with v the test function and not u function.
L = inner(P, grad(v))*dx + dot(B, v)*dx - dot(T, v)*ds
any chance you can advise if my definition is not correct and this is what is causing the strange behavior ?
my code is a combination of few documentations i have found online
this is my code
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from petsc4py import PETSc
import numpy as np
from matplotlib.ticker import FormatStrFormatter
def build_nullspace(V):
# Initialize the nullspace basis with the required size
x = Function(V).vector()
nullspace_basis = [x.copy() for _ in range(3)] # Adjusted for 3 vectors
# Define the basis for translations in the y and z directions
V.sub(1).dofmap().set(nullspace_basis[0], 1.0) # Translation in y
V.sub(2).dofmap().set(nullspace_basis[1], 1.0) # Translation in z
# Define the basis for rotation about the x-axis
V.sub(1).set_x(nullspace_basis[2], 1.0, 2) # y component affects z direction
V.sub(2).set_x(nullspace_basis[2], -1.0, 1) # z component affects y direction
for x in nullspace_basis:
x.apply("insert")
basis = VectorSpaceBasis(nullspace_basis)
basis.orthonormalize()
return basis
class Problem(NonlinearProblem):
def __init__(self, J, F, bcs):
self.bilinear_form = J
self.linear_form = F
self.bcs = bcs
NonlinearProblem.__init__(self)
def F(self, b, x):
assemble(self.linear_form, tensor=b)
for bc in self.bcs:
bc.apply(b, x)
NS = build_nullspace(V)
NS.orthogonalize(b)
def J(self, A, x):
assemble(self.bilinear_form, tensor=A)
for bc in self.bcs:
bc.apply(A)
class CustomSolver(NewtonSolver):
def __init__(self):
NewtonSolver.__init__(self, mesh.mpi_comm(),
PETScKrylovSolver(), PETScFactory.instance())
# Set solver options directly
self.parameters["linear_solver"] = "gmres"
self.parameters["preconditioner"] = "ilu"
self.parameters["maximum_iterations"] = 100000 # Increase Newton solver iterations
self.parameters["error_on_nonconvergence"] = False
self.parameters["relative_tolerance"] = 1E-6
self.parameters["absolute_tolerance"] = 1E-7
self.parameters["krylov_solver"]["absolute_tolerance"] = 1E-7
self.parameters["krylov_solver"]["relative_tolerance"] = 1E-6
self.parameters["krylov_solver"]["maximum_iterations"] = 100000 # Increase Krylov solver iterations
def solver_setup(self, A, P, problem, iteration):
self.linear_solver().set_operator(A)
NS = build_nullspace(V)
as_backend_type(A).set_nullspace(NS)
# These options are now set in __init__, but you can modify them here if needed
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
PETScOptions.set("pc_factor_mat_solver_type", "mumps")
self.linear_solver().set_from_options()
# The update_solution method within the CustomSolver class is designed to update the solution vector x during
# the Newton-Raphson iteration process for solving a nonlinear problem.
def update_solution(self, x, dx, relaxation_parameter, nonlinear_problem, iteration):
tau = 1.0
theta = min(sqrt(2.0*tau/norm(dx, norm_type="l2", mesh=V.mesh())), 1.0)
info("Newton damping parameter: %.3e" % theta)
x.axpy(-theta, dx)
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
E1, E2, nu = 50e4, 50e4, 0.49
# Create mesh and define function space
length, width, height = 1.0, 0.2, 0.2 # Dimensions of the rod
num_elements = 40 # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 5, 5)
#mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 1/2 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1/2 - DOLFIN_EPS
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
# Boundary conditions
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 1))
left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
bc2 = DirichletBC(V.sub(0), Constant(1), boundary_parts, 2)
bcs = [bc1,bc2]
#neumann = Expression(("9000"), degree=0)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_parts)
class K(fe.UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
E_global = K(materials, E1, E2, degree=0)
DG = FunctionSpace(mesh, "DG", 2)
materials_function = Function(DG)
E_values = [E1, E2]
materials_function = project(E_global, DG)
MU1, LAMBDA1 = Constant(E1/(2*(1 + nu))), Constant(E1*nu/((1 + nu)*(1 - 2*nu)))
MU2, LAMBDA2 = Constant(E2/(2*(1 + nu))), Constant(E2*nu/((1 + nu)*(1 - 2*nu)))
class MUMU(fe.UserExpression):
def __init__(self, materials, MU_0, MU_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.MU_0 = MU_0
self.MU_1 = MU_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.MU_0
else:
values[0] = self.MU_1
def value_shape(self):
return ()
MU = MUMU(materials, MU1, MU2, degree=0)
class lamabada(fe.UserExpression):
def __init__(self, materials, lam_0, lam_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.lam_0 = lam_0
self.lam_1 = lam_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.lam_0
else:
values[0] = self.lam_1
def value_shape(self):
return ()
LAMBDA = lamabada(materials, LAMBDA1, LAMBDA2, degree=0)
points = [Point(0.25, 0.1, 0.1), Point(0.75, 0.1, 0.1), Point(0.5, 0.1, 0.1), Point(0.0, 0.1, 0.1), Point(1.0, 0.1, 0.1)] # Example points
# checking the values of the parameters on each point on the rod
# Create a function space
V_scalar = FunctionSpace(mesh, 'DG', 2)
# Project LAMBDA onto this function space
lambda_projected = project(LAMBDA, V_scalar)
# Evaluate LAMBDA at the specified points
for pt in points:
print(f"LAMBDA at {pt} point: {lambda_projected(pt)}")
# Create a function space
E_scalar = FunctionSpace(mesh, 'DG', 2)
# Project LAMBDA onto this function space
E_projected = project(E_global, E_scalar)
# Evaluate LAMBDA at the specified points
for pt in points:
print(f"E at {pt}: {E_projected(pt)}")
MU_scalar = FunctionSpace(mesh, 'DG', 2)
# Project LAMBDA onto this function space
MU_projected = project(MU, MU_scalar)
# Evaluate LAMBDA at the specified points
for pt in points:
print(f"MU at {pt}: {MU_projected(pt)}")
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0)) # Body force per unit volume
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
# Project or interpolate the initial guess onto u
u.interpolate(Constant((0,0,0)))
# Kinematics
d = len(u)
I = Identity(d) # Identity tensor
F = variable(I + grad(u)) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Stored strain energy density (compressible neo-Hookean model)
W = MU/2*(tr(C) - 3 - 2*ln(J)) + LAMBDA/2*pow((J-1), 2)
n = FacetNormal(mesh)
# Total potential energy
P = diff(W, F)
# Compute first variation of Pi (directional derivative about u in the direction of v)
L = inner(P, grad(v))*dx + dot(B, v)*dx - dot(T, v)*ds(2)
# Solve variational problem
Je=derivative(L,u)
problem=Problem(Je, L, bcs)
# Reset PETSc options to a known state
PETSc.Options().clear()
# Set only the essential solver and preconditioner types
PETSc.Options().setValue("ksp_type", "gmres")
PETSc.Options().setValue("pc_type", "ilu")
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())
CL = F*F.T
# Deviatoric part
dev_stress = MU * (CL - I)
# Volumetric part
vol_stress = LAMBDA * (J-1) * I
#calculating the cauchy stress
sigma = (1/J) * dev_stress + vol_stress
S = FunctionSpace(mesh, "P", 1) # Example function space
sigma_xx = project(sigma[0, 0], S)
points_of_interest = [(0.0, 0.1, 0.1),(0.1, 0.1, 0.1), (0.4, 0.1, 0.1),(0.6, 0.1, 0.1), (0.9, 0.1, 0.1)] # Example points
for point in points_of_interest:
try:
sigma_xx_value = sigma_xx(point)
print(f"sigma_xx at {point}: {sigma_xx_value}")
except Exception as e:
print(f"Could not evaluate sigma_xx at {point}: {e}")
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
File("Materials.pvd")<< materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials
materials_function.rename("materials", "")
u.rename("Displacement Vector", "")
beam_deflection_file = fe.XDMFFile("beam_deflection.xdmf")
beam_deflection_file.parameters["flush_output"] = True
beam_deflection_file.parameters["functions_share_mesh"] = True
beam_deflection_file.write(u, 0.0)
beam_deflection_file.write(materials_function, 0.0)
# Plot solution
plot(u)
plt.show()
# List to store stress and deformation gradient values for all points
stress_values_all = []
deformation_gradient_values_all = []
# Calculate stress and deformation gradient at all points
for x in S.tabulate_dof_coordinates():
#is to return the coordinates of all degrees of freedom (DOFs) associated with the finite element function space "S."
point = Point(x)
try:
sigma_xx_value = sigma_xx(point)
deformation_gradient_xx_value = project(grad(u)[0, 0], S)(
point) # Assuming deformation_xx is the first component of the gradient
stress_values_all.append(sigma_xx_value)
deformation_gradient_values_all.append(deformation_gradient_xx_value)
except Exception as e:
print(f"Could not evaluate at {point}: {e}")
# Convert lists to numpy arrays
# Define a function to normalize data to [0, 1] range
def normalize_data(values):
min_val = np.min(values)
max_val = np.max(values)
return (values - min_val) / (max_val - min_val)
# Normalize deformation gradient and stress values
deformation_gradient_normalized = normalize_data(deformation_gradient_values_all)
stress_normalized = normalize_data(stress_values_all)
# Fit a polynomial to the normalized data
degree = 2 # Keeping the polynomial degree low
coefficients_normalized = np.polyfit(deformation_gradient_normalized, stress_normalized, degree)
# Create a polynomial function from the coefficients
polynomial_function_normalized = np.poly1d(coefficients_normalized)
# Generate points for the polynomial curve on normalized scale
deformation_gradient_fit_normalized = np.linspace(0, 1, 100)
stress_fit_normalized = polynomial_function_normalized(deformation_gradient_fit_normalized)
# Convert normalized fit data back to original scale for plotting
deformation_gradient_fit = deformation_gradient_fit_normalized * (np.max(deformation_gradient_values_all) - np.min(deformation_gradient_values_all)) + np.min(deformation_gradient_values_all)
stress_fit = stress_fit_normalized * (np.max(stress_values_all) - np.min(stress_values_all)) + np.min(stress_values_all)
# Plot the polynomial curve on original scale
plt.plot(deformation_gradient_fit_normalized, stress_fit_normalized, label=f'Normalized Polynomial',color='red',linewidth=2)
'''
plt.xlim([0.999999, 1.0000022]) # Set x-axis limits
plt.ylim([580602.05904136, 580602.05904139]) # Set y-axis limits
plt.xticks(np.linspace(0.999999, 1.0000022, 10)) # Adjust x-axis ticks
plt.yticks(np.linspace(580602.05904136, 580602.05904139, 10))
'''
plt.gca().xaxis.set_major_formatter(FormatStrFormatter('%.5f'))
plt.gca().yaxis.set_major_formatter(FormatStrFormatter('%.5f'))
# Adding labels, title, and legend as before
plt.xlabel('Deformation Gradient')
plt.ylabel('Sigma_xx')
plt.title('Stress vs. Deformation Gradient')
plt.legend()
plt.grid(True)
# Finally, display the plot
plt.show()
# Close the XDMF files
beam_deflection_file.close()