I am sorry, this is the complete code. The functions are inspired from the open source code implemented by Tianju Xue.
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from functools import partial
import ufl
import sys
import time
# ==================================
# Util function
# ===================================
def epsilon(u):
strain = as_tensor([[ u[0].dx(0), (1./2.)*(u[0].dx(1) + u[1].dx(0)) ],
[ (1./2.)*(u[0].dx(1) + u[1].dx(0)) , u[1].dx(1) ]])
return strain
# ----------------------------------------------------------------
# Model : Miehe
#
def psi_plus_elasticity(eps, lamda, mu):
delta = tr(eps)**2 - 4 * det(eps)
sqrt_delta = conditional(gt(delta, 0), sqrt(tr(eps)**2 - 4 * det(eps)), 0)
eigen_value_1 = (tr(eps) + sqrt_delta) / 2
eigen_value_2 = (tr(eps) - sqrt_delta) / 2
tr_epsilon_plus = conditional(gt(tr(eps), 0.), tr(eps), 0.)
eigen_value_1_plus = conditional(gt(eigen_value_1, 0.), eigen_value_1, 0.)
eigen_value_2_plus = conditional(gt(eigen_value_2, 0.), eigen_value_2, 0.)
return lamda / 2 * tr_epsilon_plus**2 + mu * (eigen_value_1_plus**2 + eigen_value_2_plus**2)
def psi_minus_elasticity(eps, lamda, mu):
#eps = epsilon(u)
delta = tr(eps)**2 - 4 * det(eps)
sqrt_delta = conditional(gt(delta, 0), sqrt(tr(eps)**2 - 4 * det(eps)), 0)
eigen_value_1 = (tr(eps) + sqrt_delta) / 2
eigen_value_2 = (tr(eps) - sqrt_delta) / 2
tr_epsilon_minus = conditional(lt(tr(eps), 0.), tr(eps), 0.)
eigen_value_1_minus = conditional(lt(eigen_value_1, 0.), eigen_value_1, 0.)
eigen_value_2_minus = conditional(lt(eigen_value_2, 0.), eigen_value_2, 0.)
return lamda / 2 * tr_epsilon_minus**2 + mu * (eigen_value_1_minus**2 + eigen_value_2_minus**2)
# ----------------------------------------------------------------
# Cauchy stress
def cauchy_stress_plus(eps, psi_plus):
eps = variable(eps)
energy_plus = psi_plus(eps)
sigma_plus = diff(energy_plus, eps)
return sigma_plus
def cauchy_stress_minus(eps, psi_minus):
eps = variable(eps)
energy_minus = psi_minus(eps)
sigma_minus = diff(energy_minus, eps)
return sigma_minus
# ==================================
# Linear elastic problem - pure traction test
# ===================================
ERROR = 40
set_log_level(ERROR) # log level
PROGRESS = 16
#set_log_level(PROGRESS) # log level
parameters.parse() # read paramaters from command line
# set some dolfin specific parameters
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs" # quadrature is deprecated
# ------------------
# folder name : save file Create VTK file for saving solution
# ------------------
savedir = "results_CubeLinElas_Miehe_XUE_nx3/"
vtkfileU = File(savedir+'displacement.pvd')
vtkfileSigma = File(savedir+'stress.pvd')
E, nu = 210000, 0.3 #[MPa]
G = E / (2.0*(1.0 + nu)) # [MPa]
lmbda = E * nu / ((1.0 + nu)* (1.0 - 2.0*nu))
K0, Gc = lmbda+2./3.*G, 2.7
T = 1. # final time
num_steps = 10 # number of time steps
time_step = np.linspace(0, T, num_steps)
n = 5
mesh = UnitSquareMesh(n, n,'crossed')
Nnode = len(mesh.coordinates())
coord = mesh.coordinates()
ndim = mesh.topology().dim()
V = VectorFunctionSpace(mesh,'Lagrange', 1)
Tens = TensorFunctionSpace(mesh, 'Lagrange', 1)
Disp = Expression( 't', t = 0.0, degree = 1 )
# boolean boundary function
boundary_bottom = 'near(x[1], 0)'
boundary_up = 'near(x[1], 1)'
boundary_left = 'near(x[0], 0)'
boundary_right = 'near(x[0], 1)'
bottom = CompiledSubDomain(boundary_bottom)
up = CompiledSubDomain(boundary_up)
left = CompiledSubDomain(boundary_left)
right = CompiledSubDomain(boundary_right)
boundaries = MeshFunction("size_t", mesh,1)
boundaries.set_all(0)
bottom.mark(boundaries, 1) # mark left as 1
up.mark(boundaries, 2) # mark right as 2
left.mark(boundaries, 3) # mark right as 3
right.mark(boundaries, 4) # mark right as 4
ds = Measure("ds",subdomain_data=boundaries) # left: ds(1), right: ds(2), crack : ds(4), ds(0) all the edges
bcu_bottom = DirichletBC( V.sub(1) , Constant((0.)), boundary_bottom)
bcu_right = DirichletBC( V.sub(0) , Constant(0.), boundary_right)
bcu_left = DirichletBC( V.sub(0) , Constant(0.), boundary_left)
bcu_up = DirichletBC( V.sub(1) , Disp, boundary_up )
bcu = [bcu_bottom, bcu_up, bcu_right, bcu_left]
du,v,u_ = TrialFunction(V), TestFunction(V), Function(V)
sig_funct = Function(Tens, name="Stress")
u_old = Function(V)
"""
=========================================================
Define Variational problem
=========================================================
"""
epsilon = epsilon(u_)
psi_plus = partial(psi_plus_linear_elasticity, lmbda, G)
psi_minus = partial(psi_minus_linear_elasticity, lmbda, G)
sigma_minus = cauchy_stress_minus(epsilon, psi_minus)
sigma_plus = cauchy_stress_plus(epsilon, psi_plus)
sigma = sigma_minus + sigma_plus
f = inner( sigma, epsilon(v) ) * dx + dot( Constant((0., 0.)), v ) * ds
J = derivative(f, u_, du)
#u_guess = interpolate(Expression( ('0.', 'x[0]*x[0]'), degree = 2 ),V)
#u_.assign(u_guess)
height = 1.
Uref = 2.
#x, y = x[0],x[1]
def nonzero_initial_guess():
x1 = x[0]
x2 = x[1]
u1 = Constant(0.)
u2 = Uref * x2 / height
u_initial = as_vector([u1, u2])
return u_initial
u_initial = nonzero_initial_guess()
u_.assign(project(u_initial, V))
problem = NonlinearVariationalProblem(f, u_, bcu, J)
solver = NonlinearVariationalSolver(problem)
newton_prm = solver.parameters['newton_solver']
newton_prm['relative_tolerance'] = 1e-8
newton_prm['absolute_tolerance'] = 1e-8
newton_prm['maximum_iterations'] = 100
newton_prm['krylov_solver']['nonzero_initial_guess'] = True
it = 0
for t in time_step[1:]:
it +=1
Disp.t = t
print('\n -> Time iteration : %d, Time [s] : %.3g, Ux_given [mm] : %.3g mm ' %(it, t, Disp.t))
#---------------- RESOLUTION
iter_u = solver.solve()
#---------------- PostProcess
sig_funct.assign(project(sigma, Tens))
elastic_energy = 0.5*inner(sigma, epsilon(u_))*dx
elastic_energy_density = project(0.5*inner(sigma, epsilon(u_)), V_alpha)
elastic_energy_value = assemble(elastic_energy)
Fn = assemble(sigma[0,0]*ds(2))
Ft = assemble(sigma[1,1]*ds(2))
vtkfileU << u_
vtkfileSigma << sig_funct
print( "Eel_density_max: %.8g" % (elastic_energy_density.vector().max()))
print("(Fn, Ft): (%.8g,%.8g)"%(Fn,Ft))
print("\nElastic energies : (%g)"%(elastic_energy_value))
print("-------------------------------------------------------------")
