Hi,
I am trying to solve a 3 point bending problem by solving both a LinearVariational problem for u - displacement, and a NonLinearVariational problem for phi - the phase field crack variable.
When getting into the convergence loop, it doesn’t seem to update the coefficients inside of the residual that depends on another problem for both problems. So when I initialize pnew, pold to 1, it will always keep that value and E_phi will always be zero, like it is not updating.
The output would always look like:
0 SNES Function norm 0.000000000000e+00
Iterations: 1 , time=0.15 L2 error phi 0.0
Surface energy = -1.07439e-14
Is there anything I am doing wrong/ missing?
*** DOLFIN version: 2019.1.0
*** Git changeset: 15b823f2c45d3036b6af931284d0f8e3c77b6845
Please see attached. I tried including the mesh file but the xml is too big, please let me know of a way to post it so the code can be tried.
Thank you so much!
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
import dolfin as d
set_log_active(False)
# set_log_level(ERROR)
def lameConsts(Emodulus, vpoisson):
return [Emodulus / (2 * (1 + vpoisson)),
Emodulus * vpoisson /(1 - (vpoisson**2))]
# Create mesh and function spaces.
filename = './mesh-crack-short.xml'
mesh = Mesh(filename)
# Define Space
V = FunctionSpace(mesh, 'CG', 1)
W = VectorFunctionSpace(mesh, 'CG', 1)
p = Function(V)
dp = TrialFunction(V)
q = TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
up = Function(V)
# Parameters
ATCw = 2. / 3. # 0.716575
Emodulus = 2.857895 # N / mm2
vpoisson = 0.35
[mu, lmbda] = lameConsts(Emodulus, vpoisson)
Gc = 1. # 0.00181133
l = 5
print(mu, lmbda)
# Boundary regions
crackc = [-25, -1.85]
crack_coords1 = [-25.3, 4.5]
crack_coords2 = [-24.7, 4.5]
# Short mesh
rcorner = [0.4, -8.2] # Real location of pin
lcorner = [-50.4000015, -8.19999981]
# Variable definitions
def g_KKL(p):
return (p**3) * (4. - 3. * p)
def w_KKL(p):
return 1. - g_KKL(p)
def gp_KKL(p): # Derivative of g
return 12. * (p**2) * (1. - p)
def wp_KKL(p): # Derivative of w
return -gp_KKL(p)
def gpp_KKL(p): # Derivative of g
return 24. * p - 36. * (p**2)
def wpp_KKL(p): # Derivative of w
return -gpp_KKL(p)
def g_AT2(p):
return (p**2)
def w_AT2(p):
return (1. - p)**2
def gp_AT2(p): # Derivative of g
return 2. * p
def wp_AT2(p): # Derivative of w
return 2. * (p - 1.)
def gpp_AT2(p): # Derivative of g
return 2.
def wpp_AT2(p): # Derivative of w
return 2.
def g_AT1(p):
return (p**2)
def w_AT1(p):
return (1. - p)
def gp_AT1(p): # Derivative of g
return 2.*p
def wp_AT1(p): # Derivative of w
return -1.
def gpp_AT1(p): # Derivative of g
return 2.
def wpp_AT1(p): # Derivative of w
return 0
def epsilon(u):
return sym(grad(u))
def sigma(u):
return 2.0 * mu * epsilon(u) + lmbda * tr(epsilon(u)) * Identity(len(u))
def psi(u):
return 0.5 * inner(sigma(u), epsilon(u))
# Boundary conditions
class toppoint(SubDomain): # push domain
def inside(self, x, on_boundary):
tol = 1e-3
return abs(x[1] - crack_coords1[1]) < tol and \
crack_coords1[0] <= x[0] <= crack_coords2[0]
class crackTip(SubDomain): # push domain
def inside(self, x, on_boundary):
tol = 1e-3
return abs(x[1] - crackc[1]) < tol and abs(x[0] - crackc[0]) < tol
class cornerpoint1(SubDomain):
def inside(self, x, on_boundary):
tol = 1e-3
return abs(x[1] - lcorner[1]) < tol and abs(x[0] - lcorner[0]) < tol
class cornerpoint2(SubDomain):
def inside(self, x, on_boundary):
tol = 1e-3
return abs(x[1] - rcorner[1]) < tol and abs(x[0] - rcorner[0]) < tol
Toppoint = toppoint()
Cornerpoint1 = cornerpoint1()
Cornerpoint2 = cornerpoint2()
CrackTip = crackTip()
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
Toppoint.mark(boundaries, 1)
Cornerpoint1.mark(boundaries, 2)
Cornerpoint2.mark(boundaries, 3)
CrackTip.mark(boundaries, 4)
# Initial value, lower and upper bounds for damage
lb = interpolate(Constant(0), V)
ub = interpolate(Constant(1), V)
el = V.ufl_element()
# Extract top boundary nodes
topbcoords = []
meshitr = SubsetIterator(boundaries, 1)
for itr in meshitr:
for ver in vertices(itr):
topbcoords.append([ver.midpoint().x(), ver.midpoint().y()])
# print(len(topbcoords))
u_Lx = Expression('-t', t=0.0, degree=2) # load
bclx1 = DirichletBC(W, Constant((0.0, 0.0)), Cornerpoint1, method='pointwise')
bclx2 = DirichletBC(W, Constant((0.0, 0.0)), Cornerpoint2, method='pointwise')
bcty = DirichletBC(W.sub(1), u_Lx, Toppoint)
bc_u = [bclx1, bclx2, bcty]
# Crack Boundary condition
bc_phi = [] #DirichletBC(V, Constant(0.0), CrackTip, method='pointwise')
# Define variational form
unew, uold = Function(W), Function(W)
pnew, pold = Function(V), Function(V)
pold.assign(Constant(1.0))
pnew.assign(Constant(1.0))
# Define problems
# Linear problem u
E_du = g_AT1(pnew) * inner(grad(v), sigma(u)) * dx
prob_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
solver_disp = LinearVariationalSolver(prob_disp)
# Nonlinear phi
Gc_fourCw = Gc / (4. * ATCw)
E_phi = (Gc_fourCw * (wp_AT1(p) * inner(p, q) / l +
2 * l * inner(grad(p), grad(q))) +
gp_AT1(p) * psi(unew) * inner(p, q)) * dx
J_phi = (Gc_fourCw * (wpp_AT1(p) * inner(dp, q) / l +
2 * l * inner(grad(dp), grad(q))) +
gpp_AT1(p) * psi(unew) * inner(dp, q)) * dx
prob_phi = NonlinearVariationalProblem(E_phi, pnew, bc_phi, J=J_phi)
prob_phi.set_bounds(lb.vector(), ub.vector())
solver_phi = NonlinearVariationalSolver(prob_phi)
solver_phi.parameters['nonlinear_solver'] = 'snes'
solver_phi.parameters["snes_solver"]["maximum_iterations"] = 50
solver_phi.parameters["snes_solver"]["report"] = True
# Optimization variables
t = 0.
deltaT = 5e-2
ndeltaT = 1e-5
tol = 1e-6
ut = 1e-1 # Magnitude of the remote displacement (mm)
max_t = 100
limit = 1
counter = 0
tolsol = 1e-6
# Outputs
pnew.rename("phi", "phase field")
unew.rename("u", "displacement field")
file_1 = open('LoadDisp.txt', 'w')
output_phi = File("./bendingnewResultsDir/phi.pvd")
output_u = File("./bendingnewResultsDir/u.pvd")
# Staggered scheme
while t <= max_t:
u_Lx.t = 1000 * t
iter = 0
err_phi = 1
err = 1
while err > tol:
iter += 1
solver_disp.solve()
# test error
alpha_error = pnew.vector() - pold.vector()
err_phi = alpha_error.norm('linf')
F_phi = Gc_fourCw * (w_AT1(pnew) / l + l * inner(grad(pnew), grad(pnew))) * dx
solver_phi.solve()
err_u = errornorm(unew, uold, norm_type='l2', mesh=mesh)
err = max(err_u, err_phi)
uold.assign(unew)
pold.assign(pnew)
print('Iterations:', iter, ', time=%2.2g' % t, "L2 error phi",
err_phi, "L2 error u", err_u)
F_u = .5 * g_AT1(pnew) * inner(epsilon(unew), sigma(unew)) * dx
F_phi = Gc_fourCw * (w_AT1(pnew) / l + l * inner(grad(pnew), grad(pnew))) * dx
print(np.array(assemble(E_phi)))
print("Bulk energy= %g" % assemble(F_u),
"Surface energy = %g" % assemble(F_phi))
# Update upper bounds (Only cracks grow)
prob_phi.set_bounds(lb.vector(), pnew.vector())
if counter % 1 == 0:
# Extract reaction force
fyt = 0
F_reaction = assemble(E_du) * unew.vector()
Ft_func = Function(W, F_reaction)
for coord in topbcoords:
fy = Ft_func(coord[0], coord[1])
fyt += abs(fy[1]) # N
# Outputs
print(pnew(crackc[0], crackc[1]))
plt.scatter(t * ut, fyt)
plt.pause(0.05)
output_phi << pnew, t
output_u << unew, t
file_1.write(str(t * ut) + "\t")
file_1.write(str(fyt) + "\n")
counter += 1
t += deltaT
file_1.close()
print('Simulation Done with no error')