Dear Fenicsx users,
I am trying to implement 3D elasto-plasticity using fenicsx.
However, It seems that one of my variable, “c2”, is wrongly calculated.
For example, the variable “c2_check” is 139809.6487373808, which is calculated based on recorded values, but the “test_c2” which is recorded in the loop is zero.
Is there something wrong with my codes?
Very thanks for your assistance.
#-----------------------Import basic library-----------------------
import pyvista
import dolfinx
import pandas as pd
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile, VTKFile
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import numpy as np
import math
import pandas
import basix
import matplotlib.pyplot as plt
i=ufl.Index()
j=ufl.Index()
k=ufl.Index()
l=ufl.Index()
p=ufl.Index()
q=ufl.Index()
r=ufl.Index()
s=ufl.Index()
#-----------------------Definition of quadrature parameters -----------------------
basis_parameter = {"quadrature_degree" : 2, # quadrature_degree
"basix_element_type": basix.CellType.tetrahedron,
"mesh_element_type" : mesh.CellType.tetrahedron,
}
quadrature_degree = basis_parameter["quadrature_degree"] #Set the qudrature degree
quadrature_points, weight = basix.make_quadrature(basis_parameter["basix_element_type"], quadrature_degree)
#-----------------------Geometry-----------------------
x_length = 500 # length on x-direction (mm)
y_length = 100 # length on x-direction (mm)
z_length =100 # length on x-direction (mm)
#----------------------Selection of mesh size -----------------------
estimated_mesh_size=math.gcd(math.gcd(int(x_length*1000), int(y_length*1000)), int(z_length*1000))/1000 #Finding Greatest Common Factor
print("estimated mesh size: "+str(estimated_mesh_size))
mesh_size=50
#-----------------------Definition of the number of nodes-----------------------
number_of_x_element=int(round(x_length/mesh_size))
number_of_y_element=int(round(y_length/mesh_size))
number_of_z_element=int(round(z_length/mesh_size))
domain_1 = mesh.create_box(MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([x_length, y_length, z_length])],
[number_of_x_element, number_of_y_element, number_of_z_element],
cell_type=basis_parameter["mesh_element_type"])
#-----------------------Checking the number of element-----------------------
number_of_nodes_local = domain_1.topology.index_map(0).size_local
number_of_nodes_global= domain_1.topology.index_map(0).size_global
number_of_nodes_ghost= domain_1.topology.index_map(0).num_ghosts
number_of_cells_local = domain_1.topology.index_map(3).size_local
number_of_cells_global= domain_1.topology.index_map(3).size_global
number_of_cells_ghost=domain_1.topology.index_map(3).num_ghosts
#-----------------------Definition of element -----------------------
element_node_Scalar = basix.ufl.element("P", domain_1.topology.cell_name(), degree=quadrature_degree, shape=())
element_node_Tensor_3 = basix.ufl.element("P", domain_1.topology.cell_name(), degree=quadrature_degree, shape=(3,))
element_quad_Scalar = basix.ufl.quadrature_element(domain_1.topology.cell_name(), degree=quadrature_degree)
element_quad_block_Tensor_33 = basix.ufl.blocked_element(element_quad_Scalar, shape=(3,3), symmetry=False)
element_quad_block_Tensor_333 = basix.ufl.blocked_element(element_quad_Scalar, shape=(3,3,3))
element_quad_block_Tensor_3333 = basix.ufl.blocked_element(element_quad_Scalar, shape=(3,3,3,3))
element_quad_block_Tensor_333333 = basix.ufl.blocked_element(element_quad_Scalar, shape=(3,3,3,3,3,3))
#-----------------------Definition of finite functionspace -----------------------
function_space_1_cell_Scalar = fem.functionspace(domain_1,("DG",0,()))
function_space_1_cell_Tensor_3 = fem.functionspace(domain_1,("DG",0,(3,)))
function_space_1_cell_Tensor_33 = fem.functionspace(domain_1,("DG",0,(3,3)))
function_space_1_node_Scalar = fem.functionspace(domain_1,element_node_Scalar)
function_space_1_node_Tensor_3 = fem.functionspace(domain_1,element_node_Tensor_3)
function_space_1_quad_Scalar = fem.functionspace(domain_1, element_quad_Scalar)
function_space_1_quad_block_Tensor_33 = fem.functionspace(domain_1, element_quad_block_Tensor_33,)
function_space_1_quad_block_Tensor_333 = fem.functionspace(domain_1, element_quad_block_Tensor_333,)
function_space_1_quad_block_Tensor_3333 = fem.functionspace(domain_1, element_quad_block_Tensor_3333,)
function_space_1_quad_block_Tensor_333333 = fem.functionspace(domain_1, element_quad_block_Tensor_333333,)
#Check dof
number_of_dofs_local = function_space_1_node_Tensor_3.dofmap.index_map.size_local * function_space_1_node_Tensor_3.dofmap.index_map_bs
number_of_dofs_global = function_space_1_node_Tensor_3.dofmap.index_map.size_global * function_space_1_node_Tensor_3.dofmap.index_map_bs
number_of_dofs_Tensor_33_local = function_space_1_quad_block_Tensor_33.dofmap.index_map.size_local * function_space_1_quad_block_Tensor_33.dofmap.index_map_bs
number_of_dofs_Tensor_33_global = function_space_1_quad_block_Tensor_33.dofmap.index_map.size_global * function_space_1_quad_block_Tensor_33.dofmap.index_map_bs
#-----------------------Default setting-----------------------
boundary_condition_parameter={"displacement" : -10, #(mm)
"point_load":-0, #(N)
"distributed_load":-0/y_length/x_length, # (N/mm^2)
"density": 0, # density (tonne/mm^3) steel:7.85E-9, concrete 2.45E-9
"gravity": 0 # gravity (mm/s^2) -9.81E+3
}
facet_dimension = domain_1.topology.dim - 1 #dimension of facet
edge_dimension = domain_1.topology.dim - 2 #dimension of line
point_dimension = domain_1.topology.dim - 3 #dimension of point
if (boundary_condition_parameter["displacement"]!=0) and ((boundary_condition_parameter["point_load"]==0) or (boundary_condition_parameter["distributed_load"]==0)):
control_target="displacement"
elif (boundary_condition_parameter["displacement"]==0) and ((boundary_condition_parameter["point_load"]!=0) or (boundary_condition_parameter["distributed_load"]!=0)):
control_target="load"
else:
print("You have to select one among the displacement or load boundary condition.")
#-----------------------Definition of fix boundary (support)-----------------------
def fixed_node(x): #Definition of location of fixed_boundary
return np.isclose(x[0], 0)
fixed_boundary_node_1 = mesh.locate_entities_boundary(domain_1,
point_dimension,
fixed_node) # Function for identifying specific condition
fixed_boundary_displacement_1 = np.array([0, 0, 0], dtype=default_scalar_type) # Definition of the magnitute of displacement
fixed_boundary_node_dofs_1=fem.locate_dofs_topological(function_space_1_node_Tensor_3, point_dimension, fixed_boundary_node_1)
fixed_boundary_condition_1 = fem.dirichletbc(fixed_boundary_displacement_1,
fixed_boundary_node_dofs_1,
function_space_1_node_Tensor_3)
#-----------------------Definition of displacement boundary-----------------------
def displacement_node(x): #Definition of location of displacement_boundary
return np.logical_and(np.isclose(x[0], x_length,atol=5),np.isclose(x[2], z_length,atol=5))
displacement_boundary_node_1 = mesh.locate_entities_boundary(domain_1,
point_dimension,
displacement_node) # Function for identifying specific condition
displacement_boundary_node_dofs_1 = fem.locate_dofs_topological(function_space_1_node_Tensor_3.sub(2),
point_dimension,
displacement_boundary_node_1)
displacement_boundary_displacement_1 = np.array(boundary_condition_parameter["displacement"], dtype=default_scalar_type) # Definition of the magnitute of displacement
#----------------------- Definition of Body force ----------------------
body_force = fem.Constant(domain_1, default_scalar_type((0, 0, -boundary_condition_parameter["density"] * boundary_condition_parameter["gravity"])))
#----------------------- Definition of distributed load ----------------------
def distributed_load_facet(x): #Definition of location of top_boundary
return np.isclose(x[2], z_length,atol=5)
distributed_load_facet_1= mesh.locate_entities_boundary(domain_1, facet_dimension ,distributed_load_facet) #select indices of top facets
distributed_load_facet_1_marked = np.full_like(distributed_load_facet_1, 1) # marking using 1
distributed_load_facet_1_tag = mesh.meshtags(domain_1, facet_dimension, distributed_load_facet_1, distributed_load_facet_1_marked) # making tag
# distributed_load_1 = fem.Constant(domain_1, default_scalar_type((0, 0, boundary_condition_parameter["distributed_load"]))) #input the magnitude of load
#function for interpolation
class distributed_load_class():
def __init__(self, load):
self.step = 0
self.load=load
def eval(self, x): # the shape of output of eval function must be matched with the dimension of fem.Function (data dimension1, data dimension2 ,...., the number of points)
load_tensor=np.full((x.shape[0],x.shape[1]), self.load*self.step)
load_tensor[0:2,:]=0
return load_tensor
distributed_load_1_class=distributed_load_class(load=boundary_condition_parameter["distributed_load"])
distributed_load_1=fem.Function(function_space_1_node_Tensor_3, name = "distributed_load_1")
#----------------------- Definition of measures for boundary condition ----------------------
ds = ufl.Measure("ds", domain=domain_1,subdomain_data=distributed_load_facet_1_tag) #for surface integral
# ds_distributed_load = ufl.Measure("ds", domain=domain_1,subdomain_data=distributed_load_facet_1_tag) #for surface integral
# ds_point_load = ufl.Measure("ds", domain=domain_1,subdomain_data=point_load_facet_1_tag) #for surface integral
#-----------------------Definition of basic function -----------------------
def kron_delta_calculation(i, j):
return 1 if i == j else 0
def fourth_order_unit_symmetic_tensor_calculation(type_=None):
I = np.zeros((3, 3, 3, 3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
I[i, j, k, l] = (kron_delta_calculation(i, k) * kron_delta_calculation(j, l) + kron_delta_calculation(i, l) * kron_delta_calculation(j, k)) / 2.0
if type_=="ufl":
return ufl.as_tensor(I.tolist())
else:
return I
def fourth_order_unit_deviatoric_tensor_calculation(type_=None):
I=fourth_order_unit_symmetic_tensor_calculation()
one = np.eye(3)
one_outer_one = np.tensordot(one, one,0)
I_dev = I - (1/3) * one_outer_one
if type_=="ufl":
return ufl.as_tensor(I_dev.tolist())
else:
return I_dev
def lame_parameters(E_0,nu):
lambda_=(E_0 * nu)/((1. + nu) * (1. - 2. * nu))
mu= E_0 / (2. * (1. + nu))
return lambda_, mu
def elastic_stiffness_Tensor_calculation(lambda_,mu,type_=None):
D = np.zeros((3, 3, 3, 3))
for i in range(3):
for j in range(3):
for k in range(3):
for l in range(3):
D[i, j, k, l] = lambda_ * (i == j) * (k == l) + \
2*mu * ((i == k) * (j == l) + (i == l) * (j == k))/2
if type_=="ufl":
return ufl.as_tensor(D.tolist())
else:
return D
def varepsilon_calculation(v_):
return ufl.sym(ufl.grad(v_))
def deviatoric_stress_calculation(sigma_Tensor,type_=None):
if type_=="ufl":
return sigma_Tensor-1/3*ufl.tr(sigma_Tensor)*ufl.Identity(3)
else:
return sigma_Tensor-1/3*np.trace(sigma_Tensor)*np.Identity(3)
def norm(Tensor,type_=None):
if type_=="ufl":
return ufl.sqrt(ufl.inner(Tensor, Tensor))
else:
return np.sqrt(np.tensordot(Tensor, Tensor,2))
def von_mises_yield_criteria_calculation(gamma_hat, mu, eta_Tensor, sigma_y_0, beta, H, n_e_p,type_=None):
if type_=="ufl":
return norm(eta_Tensor,"ufl")-(2*mu+2/3*beta*H)*gamma_hat-ufl.sqrt(2/3)*(sigma_y_0+(1-beta)*H*n_e_p)-2/3*(1-beta)*H*gamma_hat
else:
return norm(eta_Tensor)-(2*mu++2/3*beta*H)*gamma_hat-np.sqrt(2/3)*(sigma_y_0+(1-beta)*H*n_e_p)-2/3*(1-beta)*H*gamma_hat
def elastoplastic_stiffness_Tensor_calculation(D, c1,c2,N_Tensor, type_=None):
if type_=="ufl":
# Compute the outer product N_i N_j N_k N_l
outer_product = ufl.outer(N_Tensor, N_Tensor)
I_dev=fourth_order_unit_deviatoric_tensor_calculation(type_="ufl")
# Compute the tensor D_ijkl^alg
D_alg = D - (c1 - c2) * outer_product - c2 * I_dev
return D_alg
else:
# Compute the outer product N_i N_j N_k N_l
outer_product = np.tensordot(N_Tensor, N_Tensor,0)
I_dev=fourth_order_unit_deviatoric_tensor_calculation()
# Compute the tensor D_ijkl^alg
D_alg = D - (c1 - c2) * outer_product - c2 * I_dev
return D_alg
def return_mapping(delta_u, n_sigma_Tensor, n_e_p, mu, D, sigma_y_0, H, beta, n_alpha_Tensor, n_gamma_hat) :
#trial stress calculation
tr_sigma_Tensor = n_sigma_Tensor+ ufl.as_tensor(D[i,j,k,l]*varepsilon_calculation(delta_u)[k,l],(i,j))
#remove back stress
tr_eta_Tensor=deviatoric_stress_calculation(tr_sigma_Tensor,type_="ufl")-n_alpha_Tensor
#calculate flow direction
norm_tr_eta_Tensor=norm(tr_eta_Tensor,type_="ufl")
N_Tensor=1/norm_tr_eta_Tensor*tr_eta_Tensor
#evaluate yield state
n_gamma_hat.interpolate(lambda x: np.full((x.shape[1]), 0.))
yield_criterion = von_mises_yield_criteria_calculation(n_gamma_hat, mu, tr_eta_Tensor, sigma_y_0, beta, H, n_e_p,"ufl")
#evaluate yield state
n1_gamma_hat = ufl.conditional(yield_criterion < 0., 0., yield_criterion / (2.* mu + 2/3*H))
n1_e_p=ufl.conditional(yield_criterion < 0., n_e_p, n_e_p+ufl.sqrt(2./3.)*n1_gamma_hat)
n1_sigma_Tensor=ufl.conditional(yield_criterion < 0., tr_sigma_Tensor, tr_sigma_Tensor - 2.* mu * n1_gamma_hat * N_Tensor)
n1_alpha_Tensor=ufl.conditional(yield_criterion < 0., n_alpha_Tensor, n_alpha_Tensor+2./3.*beta*H*n1_gamma_hat*N_Tensor)
c1=ufl.conditional(yield_criterion < 0.,0.,4.*mu**2./(2.*mu+2./3.*H))
c2=ufl.conditional(yield_criterion < 0.,0.,4.*mu**2.*n1_gamma_hat/norm_tr_eta_Tensor)
n1_D_alg=ufl.conditional(yield_criterion < 0., D, elastoplastic_stiffness_Tensor_calculation(D, c1,c2,N_Tensor,"ufl"))
return n1_sigma_Tensor, n1_e_p, n1_alpha_Tensor, n1_D_alg, tr_sigma_Tensor, tr_eta_Tensor, norm_tr_eta_Tensor, N_Tensor, yield_criterion, n1_gamma_hat, c1,c2
def quadature_interpolation(target_function, quadrature_points, save_function):
'''
See https://github.com/FEniCS/dolfinx/issues/2243
'''
mesh = save_function.ufl_function_space().mesh
map_c = mesh.topology.index_map(mesh.topology.dim)
num_cells = map_c.size_local + map_c.num_ghosts
cells = np.arange(0, num_cells, dtype=np.int32)
expr_expr = fem.Expression(target_function, quadrature_points)
expr_eval = expr_expr.eval(mesh, cells)
return expr_eval.reshape(-1)
# -----------------------Definition of measures and trial and test functions -----------------------
v_ = ufl.TrialFunction(function_space_1_node_Tensor_3)
u_ = ufl.TestFunction(function_space_1_node_Tensor_3)
#Check the shape of trial and test function
v_.ufl_shape
u_.ufl_shape
dx = ufl.Measure("dx",domain=domain_1, metadata={"quadrature_degree": quadrature_degree, "quadrature_scheme": "default"} )
n = ufl.FacetNormal(domain_1)
#-----------------------Definition of material model -----------------------
elasto_plasticity_3D_parameter = {"E_0" : 200000., # young's modulus (MPa)
"nu" : 0.3, #poisson's ratio
"sigma_y_0" : 400., # yield stress (MPa)
"H" : 0., # Isotropic hardening coefficient (MPa)
"beta" : 1., # Kinematic hardening coefficient (MPa)
}
#Set the invariant constant over the mesh
lambda_,mu=lame_parameters(elasto_plasticity_3D_parameter["E_0"],elasto_plasticity_3D_parameter["nu"])
D=elastic_stiffness_Tensor_calculation(lambda_,mu,type_="ufl")
lambda_= fem.Constant(domain_1, default_scalar_type(lambda_))
mu= fem.Constant(domain_1, default_scalar_type(mu))
E_0 = fem.Constant(domain_1, default_scalar_type(elasto_plasticity_3D_parameter["E_0"]))
nu = fem.Constant(domain_1, default_scalar_type(elasto_plasticity_3D_parameter["nu"]))
sigma_y_0 = fem.Constant(domain_1, default_scalar_type(elasto_plasticity_3D_parameter["sigma_y_0"]))
H = fem.Constant(domain_1, default_scalar_type(elasto_plasticity_3D_parameter["H"]))
beta = fem.Constant(domain_1, default_scalar_type(elasto_plasticity_3D_parameter["beta"])) #Set the constant over the mesh
#Set the variables over the mesh
n_alpha_Tensor = fem.Function(function_space_1_quad_block_Tensor_33, name = "n_alpha_Tensor")
n_sigma_Tensor = fem.Function(function_space_1_quad_block_Tensor_33, name = "n_sigma_Tensor")
n_e_p = fem.Function(function_space_1_quad_Scalar, name = "n_e_p")
n_D_alg = fem.Function(function_space_1_quad_block_Tensor_3333, name = "n_D_alg")
tr_sigma_Tensor = fem.Function(function_space_1_quad_block_Tensor_33, name = "tr_sigma_Tensor")
tr_eta_Tensor = fem.Function(function_space_1_quad_block_Tensor_33, name = "tr_eta_Tensor")
norm_tr_eta_Tensor = fem.Function(function_space_1_quad_Scalar, name = "norm_tr_eta_Tensor")
N_Tensor = fem.Function(function_space_1_quad_block_Tensor_33, name = "N_Tensor")
yield_criterion= fem.Function(function_space_1_quad_Scalar, name = "yield_criterion")
gamma_hat= fem.Function(function_space_1_quad_Scalar, name = "gamma_hat")
c1 = fem.Function(function_space_1_quad_Scalar, name = "c1")
c2 = fem.Function(function_space_1_quad_Scalar, name = "c2")
#initial values
#Set initial values, the dimension of array is (the number of point, data dimension 1, data dimension 2, ...)
n_D_alg.x.array[:]=np.tile(elastic_stiffness_Tensor_calculation(lambda_,mu).reshape(-1),
len(quadrature_points)*number_of_cells_local)
n_u = fem.Function(function_space_1_node_Tensor_3, name = "n_u")
delta_u = fem.Function(function_space_1_node_Tensor_3, name = "delta_u")
Delta_u = fem.Function(function_space_1_node_Tensor_3, name = "Delta_u")
#-----------------------Defenition of record -----------------------
Record_n_u_frame=pd.DataFrame()
Record_delta_u_frame=pd.DataFrame()
Record_n_sigma_Tensor_frame=pd.DataFrame()
Record_n_alpha_Tensor_frame=pd.DataFrame()
Record_n_D_alg_frame=pd.DataFrame()
Record_tr_sigma_Tensor_frame=pd.DataFrame()
Record_tr_eta_Tensor_frame=pd.DataFrame()
Record_norm_tr_eta_Tensor_frame=pd.DataFrame()
Record_N_Tensor_frame=pd.DataFrame()
Record_yield_criterion_frame=pd.DataFrame()
Record_gamma_hat_frame=pd.DataFrame()
Record_c1_frame=pd.DataFrame()
Record_c2_frame=pd.DataFrame()
Record_residual_frame=pd.DataFrame()
#-----------------------Initial record -----------------------
Record_n_u_frame["initial"]=n_u.x.array[:]
Record_delta_u_frame["initial"]=delta_u.x.array[:]
Record_n_sigma_Tensor_frame["initial"]=n_sigma_Tensor.x.array[:]
Record_n_alpha_Tensor_frame["initial"]=n_alpha_Tensor.x.array[:]
Record_n_D_alg_frame["initial"]=n_D_alg.x.array[:]
Record_tr_sigma_Tensor_frame["initial"]=tr_sigma_Tensor.x.array[:]
Record_tr_eta_Tensor_frame["initial"]=tr_eta_Tensor.x.array[:]
Record_norm_tr_eta_Tensor_frame["initial"]=norm_tr_eta_Tensor.x.array[:]
Record_N_Tensor_frame["initial"]=N_Tensor.x.array[:]
Record_yield_criterion_frame["initial"]=yield_criterion.x.array[:]
Record_gamma_hat_frame["initial"]=gamma_hat.x.array[:]
Record_c1_frame["initial"]=c1.x.array[:]
Record_c2_frame["initial"]=c2.x.array[:]
#-----------------------Defenition of solver parameters -----------------------
solver_parameter={"type":"Newton-raphson",
"tolerance":1e-36, # parameters of the Newton-Raphson procedure
"max_step" : 5,
"max_iteration": 20}
step_list = np.linspace(0, 1, solver_parameter["max_step"]+1)[1:]
results = np.zeros((solver_parameter["max_step"]+1, 2))
#For arc-length method
Delta_lambda=0
delta_lambda=0
beta_arc_length=1
for (step_index, step) in enumerate(step_list):
#Define bilinear form for tangential siffness matrix
a = ufl.inner(
ufl.as_tensor(varepsilon_calculation(v_)[k,l]*n_D_alg[i,j,k,l],(i,j)),
varepsilon_calculation(u_))*dx
#Define external force
distributed_load_1_class.step = step
distributed_load_1.interpolate(distributed_load_1_class.eval)
"""
interpolate(u0: Callable | Expression | Function, cells0: ndarray | None = None, cells1: ndarray | None = None)→ None[source]
Interpolate an expression on to the function.
Parameters :
u0 – Callable function, Expression or Function to interpolate.
cells0 – Cells in mesh associated with u0 to interpolate over. If None then all cells are interpolated over.
cells1 – Cells in the mesh associated with self to interpolate over. If None, then taken to be the same cells as cells0. If cells1 is not None, then it must have the same length as cells0.
"""
#Define linear form for internal energy
L = (ufl.dot(body_force,v_)*dx+ufl.dot(distributed_load_1,v_)*ds(1) #External forces
-ufl.inner(varepsilon_calculation(v_),n_sigma_Tensor)*dx) #Internal energy
L_f_ext = (ufl.dot(body_force,v_)*dx+ufl.dot(distributed_load_1,v_)*ds(1)) #External forces
#Initial setting
Delta_u.interpolate(lambda x: np.full((x.shape[0],x.shape[1]), 0))
number_of_iteration= 0
convergence_criteria=10000
# Report
print("step:", str(step_index+1))
while convergence_criteria > solver_parameter["tolerance"] and number_of_iteration < solver_parameter["max_iteration"]:
#Define displacement
if control_target=="load":
delta_displacement_boundary_displacement_1=displacement_boundary_displacement_1
if control_target=="displacement": #If we perform displacement control, the input displacement shoud be delta_displacement, and the delta_displacement have to be zero after first iteration.
if number_of_iteration==0:
if step_index==0:
delta_displacement_boundary_displacement_1=displacement_boundary_displacement_1*(step)
else:
delta_displacement_boundary_displacement_1=displacement_boundary_displacement_1*(step-step_list[step_index-1])
else:
delta_displacement_boundary_displacement_1=displacement_boundary_displacement_1*0
displacement_boundary_condition_1 = fem.dirichletbc(delta_displacement_boundary_displacement_1,
displacement_boundary_node_dofs_1,
function_space_1_node_Tensor_3.sub(2))
if control_target=="displacement":
boundary_conditions=[fixed_boundary_condition_1, displacement_boundary_condition_1]
else:
boundary_conditions=[fixed_boundary_condition_1]
#Solve
if solver_parameter["type"]=="Newton-raphson":
if control_target=="load":
Problem=LinearProblem(a, L, bcs=boundary_conditions, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
delta_u=Problem.solve()
if control_target=="displacement":
Problem=LinearProblem(a, L, bcs=boundary_conditions, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
delta_u=Problem.solve()
elif solver_parameter["type"]=="Arc-length":
if control_target=="load":
raise Exception('Arc-length for load control was not implemented')
if control_target=="displacement":
Problem1=LinearProblem(a, L_f_ext, bcs=boundary_conditions, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
delta_u_I=Problem1.solve()
Problem2=LinearProblem(a, L, bcs=boundary_conditions, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
delta_u_II=Problem2.solve()
L_f_ext_array=fem.assemble_vector(fem.form(L_f_ext)).array
delta_lambda=(-(np.dot(Delta_u.x.array,delta_u_I.x.array)+beta_arc_length**2*Delta_lambda**2*np.dot(L_f_ext_array,L_f_ext_array))
/(np.dot(Delta_u.x.array,delta_u_II.x.array)+beta_arc_length**2*Delta_lambda*np.dot(L_f_ext_array,L_f_ext_array)))
if np.isnan(delta_lambda):
delta_lambda=0
delta_u.x.array[:]=delta_lambda*delta_u_I.x.array+delta_u_II.x.array
#Delta_u and n_u update
Delta_u.x.array[:]=(Delta_u.x.array[:]+delta_u.x.array[:])
n_u.x.array[:]=(n_u.x.array[:]+delta_u.x.array[:])
#Return mapping
n1_sigma_Tensor, n1_e_p, n1_alpha_Tensor, n1_D_alg, n1_tr_sigma_Tensor, n1_tr_eta_Tensor, n1_norm_tr_eta_Tensor, n1_N_Tensor, n1_yield_criterion, n1_gamma_hat, n1_c1, n1_c2= return_mapping(delta_u, n_sigma_Tensor, n_e_p, mu, D, sigma_y_0, H, beta, n_alpha_Tensor, gamma_hat)
#Update
n_sigma_Tensor.x.array[:]=quadature_interpolation(n1_sigma_Tensor, quadrature_points, save_function=n_sigma_Tensor)
n_alpha_Tensor.x.array[:]=quadature_interpolation(n1_alpha_Tensor, quadrature_points, save_function=n_alpha_Tensor)
n_e_p.x.array[:]=quadature_interpolation(n1_e_p, quadrature_points, save_function=n_e_p)
n_D_alg.x.array[:]=quadature_interpolation(n1_D_alg, quadrature_points, save_function=n_D_alg)
tr_sigma_Tensor.x.array[:]=quadature_interpolation(n1_tr_sigma_Tensor, quadrature_points, save_function=tr_sigma_Tensor)
tr_eta_Tensor.x.array[:]=quadature_interpolation(n1_tr_eta_Tensor, quadrature_points, save_function=tr_eta_Tensor)
norm_tr_eta_Tensor.x.array[:]=quadature_interpolation(n1_norm_tr_eta_Tensor, quadrature_points, save_function=norm_tr_eta_Tensor)
N_Tensor.x.array[:]=quadature_interpolation(n1_N_Tensor, quadrature_points, save_function=N_Tensor)
yield_criterion.x.array[:]=quadature_interpolation(n1_yield_criterion, quadrature_points, save_function=yield_criterion)
gamma_hat.x.array[:]=quadature_interpolation(n1_gamma_hat, quadrature_points, save_function=gamma_hat)
c1.x.array[:]=quadature_interpolation(n1_c1, quadrature_points, save_function=c1)
c2.x.array[:]=quadature_interpolation(n1_c2, quadrature_points, save_function=c2)
#Residual calculation
residual = fem.assemble_vector(fem.form(L))
n_residual = sum(residual.array)**2
#converence calculation
if number_of_iteration==0:
delta_u_0_sum=sum(delta_u.x.array**2)
n_residual_0=n_residual
delta_u_sum=sum(delta_u.x.array**2)
convergence_criteria=delta_u_sum/delta_u_0_sum
print(" number of iteration:", number_of_iteration)
print(" convergence_criteria:", convergence_criteria)
print(" residual:", n_residual/n_residual_0)
#Record
Record_n_u_frame[str(step_index+1)+"_"+str(number_of_iteration)]=n_u.x.array[:]
Record_n_D_alg_frame[str(step_index+1)+"_"+str(number_of_iteration)]=n_D_alg.x.array[:]
Record_n_sigma_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=n_sigma_Tensor.x.array[:]
Record_n_alpha_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=n_alpha_Tensor.x.array[:]
Record_tr_sigma_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=tr_sigma_Tensor.x.array[:]
Record_tr_eta_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=tr_eta_Tensor.x.array[:]
Record_norm_tr_eta_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=norm_tr_eta_Tensor.x.array[:]
Record_N_Tensor_frame[str(step_index+1)+"_"+str(number_of_iteration)]=N_Tensor.x.array[:]
Record_yield_criterion_frame[str(step_index+1)+"_"+str(number_of_iteration)]=yield_criterion.x.array[:]
Record_gamma_hat_frame[str(step_index+1)+"_"+str(number_of_iteration)]=gamma_hat.x.array[:]
Record_c1_frame[str(step_index+1)+"_"+str(number_of_iteration)]=c1.x.array[:]
Record_c2_frame[str(step_index+1)+"_"+str(number_of_iteration)]=c2.x.array[:]
Record_delta_u_frame[str(step_index+1)+"_"+str(number_of_iteration)]=delta_u.x.array[:]
Record_residual_frame[str(step_index+1)+"_"+str(number_of_iteration)]=residual.array[:]
number_of_iteration += 1
#-----------------------For checking the material models-----------------------
target_cell=0
test_n_D_alg=Record_n_D_alg_frame.iloc[81*(target_cell):81*(target_cell+1),50].values.reshape(3,3,3,3)
test_n_sigma_Tensor=Record_n_sigma_Tensor_frame.iloc[9*(target_cell):9*(target_cell+1),50].values.reshape(3,3)
test_n_alpha_Tensor=Record_n_alpha_Tensor_frame.iloc[9*(target_cell):9*(target_cell+1),50].values.reshape(3,3)
test_gamma_hat=Record_gamma_hat_frame.iloc[target_cell,50]
test_N_Tensor=Record_N_Tensor_frame.iloc[9*(target_cell):9*(target_cell+1),50].values.reshape(3,3)
test_yield_criterion=Record_yield_criterion_frame.iloc[target_cell,50]
test_tr_sigma_Tensor=Record_tr_sigma_Tensor_frame.iloc[9*(target_cell):9*(target_cell+1),50].values.reshape(3,3)
test_tr_eta_Tensor=Record_tr_eta_Tensor_frame.iloc[9*(target_cell):9*(target_cell+1),50].values.reshape(3,3)
test_norm_tr_eta_Tensor=Record_norm_tr_eta_Tensor_frame.iloc[target_cell,50]
test_c1=Record_c1_frame.iloc[target_cell,50]
test_c2=Record_c2_frame.iloc[target_cell,50]
lambda_check,mu_check=lame_parameters(elasto_plasticity_3D_parameter["E_0"],elasto_plasticity_3D_parameter["nu"])
H_check=elasto_plasticity_3D_parameter["H"]
norm_tr_eta_Tensor_check=norm(test_tr_eta_Tensor)
gamma_hat_check=test_yield_criterion / (2* mu_check + 2/3*H_check)
c1_check=4*mu_check**2/(2*mu_check+2/3*H_check)
c2_check=4*mu_check**2*test_gamma_hat/test_norm_tr_eta_Tensor