Dear community,
I am trying to solve the 2-D elasticity equations coupled with the phase field equations to look into crack propagation. While transferring my scripts from dolfin to dolfin-x. a small issue arose for which I haven’t been able to find a satisfactory answer in the demos.
Therefore, I set up a minimal example which works fine in dolfin but does not converge in dolfin-x.
The dolfin-x version:
# First step: Create mesh and mesh tag
import gmsh
import sys
gmsh.initialize()
gmsh.model.add("square with cracks")
a = 0.5 # mm
b = 2*a # mm 2*a
H = a # mm 10*a
gdim = 2 # Geometric dimension of the mesh
ch = 0.00 #
crl = 0.0
lc = 1.0
# Create points
pt1 = gmsh.model.geo.addPoint(0, 0, 0,lc)
pt2 = gmsh.model.geo.addPoint(a, 0, 0,lc)
pt3 = gmsh.model.geo.addPoint(a, H, 0,lc)
pt4 = gmsh.model.geo.addPoint(-a, H, 0,lc)
pt5 = gmsh.model.geo.addPoint(-a, ch, 0,lc)
pt6 = gmsh.model.geo.addPoint(-a, -ch, 0,lc)
pt7 = gmsh.model.geo.addPoint(-a, -H, 0,lc)
pt8 = gmsh.model.geo.addPoint(a, -H, 0,lc)
# #Create line
l1 = gmsh.model.geo.addLine(pt1,pt2) #
l2 = gmsh.model.geo.addLine(pt2,pt3) # right
l3 = gmsh.model.geo.addLine(pt3,pt4) # top
l4 = gmsh.model.geo.addLine(pt4,pt5) # left
l5 = gmsh.model.geo.addLine(pt5,pt1) # crack
l6 = gmsh.model.geo.addLine(pt1,pt6) # crack
l7 = gmsh.model.geo.addLine(pt6,pt7) # left
l8 = gmsh.model.geo.addLine(pt7,pt8) # bottom
l9 = gmsh.model.geo.addLine(pt8,pt2) # right
# Create Lineloop
cl1 = gmsh.model.geo.addCurveLoop([l1,l2,l3,l4,l5])
cl2 = gmsh.model.geo.addCurveLoop([l6,l7,l8,l9,-l1])
# create surface
sur1 = gmsh.model.geo.addPlaneSurface([cl1])
sur2 = gmsh.model.geo.addPlaneSurface([cl2])
# update fisty to find suface and line
gmsh.model.geo.synchronize()
# Add Physical Group
sur_tag, crack_tag, bottom_tag, top_tag, left_tag, right_tag = 100, 101, 102,103,104, 105 # Physical Group tag
gmsh.model.addPhysicalGroup(2, [sur1,sur2], sur_tag)
gmsh.model.setPhysicalName(2, sur_tag, "sur")
gmsh.model.addPhysicalGroup(1, [l5,l6], crack_tag)
gmsh.model.setPhysicalName(1, crack_tag, "crack")
gmsh.model.addPhysicalGroup(1, [l8], bottom_tag)
gmsh.model.setPhysicalName(1, bottom_tag, "bottom")
gmsh.model.addPhysicalGroup(1, [l3], top_tag)
gmsh.model.setPhysicalName(1, top_tag, "top")
gmsh.model.addPhysicalGroup(1, [l4,l7], left_tag)
gmsh.model.setPhysicalName(1, left_tag, "left")
gmsh.model.addPhysicalGroup(1, [l2,l9], right_tag)
gmsh.model.setPhysicalName(1, right_tag, "right")
gmsh.model.geo.synchronize()
### set mesh
n1 = 40 #裂纹表面分为20个单元
num = [n1+1,n1+1]
gmsh.model.geo.mesh.setTransfiniteCurve(l1, num[0]) #线l1 共布num[0]个节点，用于网格划分
gmsh.model.geo.mesh.setTransfiniteCurve(l2, num[1])
gmsh.model.geo.mesh.setTransfiniteCurve(l3, 2*num[0]-1)
gmsh.model.geo.mesh.setTransfiniteCurve(l4, num[1])
gmsh.model.geo.mesh.setTransfiniteCurve(l5, num[0])
gmsh.model.geo.mesh.setTransfiniteCurve(l6, num[0])
gmsh.model.geo.mesh.setTransfiniteCurve(l7, num[1])
gmsh.model.geo.mesh.setTransfiniteCurve(l8, 2*num[0]-1)
gmsh.model.geo.mesh.setTransfiniteCurve(l9, num[1])
gmsh.model.geo.mesh.setTransfiniteSurface(sur1, "Left", [2, 3, 4,5])
gmsh.model.geo.mesh.setTransfiniteSurface(sur2, "Left", [8, 2, 6,7]) #节点排序不一样会导致网格有区别
gmsh.option.setNumber("Mesh.Smoothing", 100)
# update
gmsh.model.geo.synchronize()
# generate mesh
gmsh.model.mesh.generate(2)
gmsh.model.geo.synchronize()
gmsh.write("mesh3D.msh") # save mesh
gmsh.finalize()
# Second step: solve
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
import numpy as np
import ufl
from dolfinx.cpp.mesh import CellType
from dolfinx.fem import FunctionSpace,Function,Constant,LinearProblem,DirichletBC,locate_dofs_topological,assemble_scalar
import gmsh
import numpy as np
from mpi4py import MPI
from dolfinx.io import VTKFile
from gmsh_helpers import read_from_msh
# Create mesh and mesh tags
gmsh.initialize()
mesh, cell_tags, facet_tags = read_from_msh("mesh3D.msh", cell_data=True, facet_data=True, gdim=2)
gmsh.finalize()
# define element and FunctionSpace
e0 = ufl.VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, e0) #VectorFunctionSpace
e1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh,e1) # FunctionSpace(mesh,("CG",1))
WW = FunctionSpace(mesh, ('DG', 0))
# Boundary conditions
bc_u =[] # boundary conditions of displacement
bc_p = [] # boundary conditions of phase field
fdim = mesh.topology.dim - 1
pD = Function(W)
with pD.vector.localForm() as pD_crack_loc:
pD_crack_loc.set(1.0)
dofs_crack = locate_dofs_topological(W, fdim,facet_tags.indices[facet_tags.values == 101]) # x=0,0<y<<a, crack_tag=101
bc_p.append(DirichletBC(pD, dofs_crack))
V0 = V.sub(0).collapse()
u_left = Function(V0)
with u_left.vector.localForm() as u_left_loc:
u_left_loc.set(0)
dofs_left = locate_dofs_topological((V.sub(0), V0), fdim, facet_tags.indices[facet_tags.values == 104]) # x=0 left_tag =104
bc_left = DirichletBC(u_left, dofs_left, V.sub(0))
bc_u.append(bc_left) # add bondary condtion bc_left
V1 = V.sub(1).collapse()
u_bottom = Function(V1)
with u_bottom.vector.localForm() as u_bottom_loc:
u_bottom_loc.set(0)
dofs_bottom = locate_dofs_topological((V.sub(1), V1), fdim, facet_tags.indices[facet_tags.values == 102]) # y=0 bottom_lag=102
bc_bottom = DirichletBC(u_bottom, dofs_bottom, V.sub(1))
bc_u.append(bc_bottom) # add bondary condtion bc_left
class MyExpression:
def __init__(self):
self.t = 0.0
def eval(self, x):
return np.full(x.shape[1], self.t) # Added some spatial variation here. Expression is sin(t)*x
u_ex = MyExpression()
u_top = Function(V1)
u_top.interpolate(u_ex.eval) # interpolate expression into function
u_top.x.scatter_forward()
dofs_top = locate_dofs_topological((V.sub(1), V1), fdim, facet_tags.indices[facet_tags.values == 103]) # y=200 top=102
bc_top =DirichletBC(u_top, dofs_top, V.sub(1))
bc_u.append(bc_top) # add bondary condtion bc_left
# define measure for surface integration
ds = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags)
n = ufl.FacetNormal(mesh) # normal
# Define material properties
young = 210.0e3 # 210.0e3 N/mm/mm 70.0e3
nu = 0.3 # 0.3 0.22
lambda_ = young/((1.0+nu)*(1.0-2.0*nu)) # N/mm/mm
mu = young/(2.0*(1.0+nu)) # N/mm/mm
l, Gc =0.015, 2.7 # N/mm 0.007 2.7
#Variational formulation
uold, unew = Function(V), Function(V)
pnew, pold = Function(W), Function(W)
Hold = Function(WW)
with uold.vector.localForm() as uold_loc:
uold_loc.set(0)
with pold.vector.localForm() as pold_loc:
pold_loc.set(0)
with Hold.vector.localForm() as Hold_loc:
Hold_loc.set(0)
def epsilon(u):
return 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(u.geometric_dimension()) + 2*mu*epsilon(u)
# Define the traction 𝑇 , and body force
T = Constant(mesh, ScalarType((0, 0))) # load force N/mm/mm
f = Constant(mesh, ScalarType((0, 0))) # body force N/mm/mm/mm
u, v = ufl.TrialFunction(V), ufl.TestFunction(V) # the trial and test function for
a = ((1.0-pold)**2)*ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
# Define the linear variational problem
problem_u = LinearProblem(a, L, bcs=bc_u, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
def psi(u):
return 0.5*(lambda_+mu)*(0.5*(ufl.tr(epsilon(u))+ufl.algebra.Abs(ufl.tr(epsilon(u)))))**2+\
mu*ufl.inner(ufl.dev(epsilon(u)),ufl.dev(epsilon(u)))
def H(uold,unew,Hold):
return ufl.conditional(ufl.lt(psi(uold),psi(unew)),psi(unew),Hold)
def fracture(p):
return (p*p/(2.0*l)+0.5*l*ufl.dot(ufl.grad(p),ufl.grad(p)))
p, q = ufl.TrialFunction(W), ufl.TestFunction(W)
f = Constant(mesh, ScalarType(-6))
phi_a = Gc*l*ufl.inner(ufl.grad(p),ufl.grad(q))*ufl.dx + ((Gc/l)+2.0*H(uold,unew,Hold))*ufl.inner(p,q)*ufl.dx
phi_L = 2.0*H(uold,unew,Hold)*q*ufl.dx
# Define the linear variational problem
problem_p = LinearProblem(phi_a, phi_L, bcs=bc_p, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
#project
def project(expression, V):
u_, v_ = ufl.TrialFunction(V), ufl.TestFunction(V)
a_p = ufl.inner(u_, v_) * ufl.dx
L_p = ufl.inner(expression, v_) * ufl.dx
projection = LinearProblem(a_p, L_p)
return projection.solve()
file_p = VTKFile(MPI.COMM_WORLD, "./result/phi.pvd", "w")
t = 0 # Start time
T = 0.007 # End time
num_steps = 20 # Number of time steps
dt = (T-t)/num_steps # Time step size
u_r = 0.007
deltaT = 0.1
tol = 1e-3
load = []
disp = []
ith_tep = 0 # the ith step
#for i in range(num_steps):
while t<=1.0: # real 0.007 test 0.00002
# Update Diriclet boundary condition
t += deltaT
if t >=0.7:
deltaT = 0.0001
u_ex.t+=t*u_r
u_top.interpolate(u_ex.eval)
u_top.x.scatter_forward()
error_L2 = 1
iter = 0
i=-1
while error_L2 > tol:
iter += 1
i=i+1
# solve problem_u
unew = problem_u.solve()
unew.x.scatter_forward()
u_error = ufl.inner(unew - uold, unew - uold) * ufl.dx
error_u = np.sqrt(assemble_scalar(u_error))
# solve problem_p
pnew = problem_p.solve()
pnew.x.scatter_forward()
p_error = ufl.inner(pnew - pold, pnew - pold) * ufl.dx
error_p = np.sqrt(assemble_scalar(p_error))
error_L2 = max(error_u,error_p)
# Update solution at previous time step (u_n)
uold.x.array[:] = unew.x.array
pold.x.array[:] = pnew.x.array
Hnew = project(psi(unew), WW)
Hold.x.array[:] = Hnew.x.array
print ('Iterations:', iter,error_p, error_u)
if error_L2 < tol:
ith_tep =ith_tep+1
print ('ith_tep:', ith_tep, ', Total time', u_ex.t)
tration = ufl.dot(ufl.dot(sigma(unew), n),n)*ds(103)
energy = psi(unew)*ufl.dx
force = assemble_scalar(tration)
all_e = assemble_scalar(energy)
file_p.write_function(pnew, i)
load.append(force)
disp.append(u_ex.t)
#sur = fracture(pnew)*ufl.dx
#print(assemble_scalar(sur),force)
np.savetxt('./result/load.out',load)
np.savetxt('./result/disp.out',disp)
file_p.close()
The dolfin version
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Apr 18 19:18:56 2021
@author: wangzh
"""
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#--------------------------------------
# creare a mesh
#--------------------------------------
set_log_active(False)
nx = 1
ny = 1
num = 1
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh,"triangle", 2, 2)
editor.init_vertices(10)
editor.init_cells(8)
editor.add_vertex(0, Point(0, 0.0))
editor.add_vertex(1, Point(0.5, 0.0))
editor.add_vertex(2, Point(0.5, 0.5))
editor.add_vertex(3, Point(0.0, 0.5))
editor.add_vertex(4, Point(-0.5, 0.5))
#editor.add_vertex(5, Point(-0.5, 0.0001)) #Point(0.5,0.0001)
editor.add_vertex(5, Point(-0.5, 0.0))
editor.add_vertex(6, Point(-0.5, -0.5))
editor.add_vertex(7, Point(0.0, -0.5))
editor.add_vertex(8, Point(0.5, -0.5))
#editor.add_vertex(9, Point(-0.5, -0.0001)) #Point(-0.5, -0.0001)
editor.add_vertex(9, Point(-0.5, 0.0))
editor.add_cell(0, [0, 1, 3])
editor.add_cell(1, [1, 2, 3])
editor.add_cell(2, [0, 3, 4])
editor.add_cell(3, [0, 4, 5])
editor.add_cell(4, [0, 9, 7])
editor.add_cell(5, [6, 7, 9])
#editor.add_cell(4, [0, 6, 7])
#editor.add_cell(5, [0, 9, 6])
editor.add_cell(6, [0, 7, 8])
editor.add_cell(7, [0, 8, 1])
editor.close()
#--------------------------------------
# change num_refine for mesh refinement
#--------------------------------------
k0=-0.47/0.2
k1=-0.48/0.5
num_refine = 6 # real 6 test 4
for i in range(num_refine):
mesh=refine(mesh)
print ("number of vertices",mesh.num_vertices())
print ("number of elements",mesh.num_cells())
print ("mesh.hmin()=",mesh.hmin())
print ("mesh.hmax()=",mesh.hmax())
#plot(mesh)
file_mesh = File("meshwang.pvd")
file_mesh << mesh
# Define Space
V = FunctionSpace(mesh, 'CG', 3)
W = VectorFunctionSpace(mesh, 'CG', 1)
WW = FunctionSpace(mesh, 'DG', 0)
p, q = TrialFunction(V), TestFunction(V)
u, v = TrialFunction(W), TestFunction(W)
# Introduce manually the material parameters
Gc = 2.7
hm = mesh.hmin()
l = 2*hm #l = 0.015
#l = 0.003
lmbda = 121.1538e3
mu = 80.7692e3
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*(lmbda+mu)*(0.5*(tr(epsilon(u))+abs(tr(epsilon(u)))))**2+\
mu*inner(dev(epsilon(u)),dev(epsilon(u)))
def H(uold,unew,Hold):
return conditional(lt(psi(uold),psi(unew)),psi(unew),Hold)
# Boundary conditions
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
def Crack(x):
return abs(x[1]) < 1e-03 and x[0] <= 0.0
load = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(W, Constant((0.0,0.0)), bot)
bctop = DirichletBC(W.sub(0), load, top)
bctop1 = DirichletBC(W.sub(1), Constant(0.0), top)
bc_u = [bcbot, bctop,bctop1]
bc_phi = [DirichletBC(V, Constant(1.0), Crack)]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries,1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
# Variational form
unew, uold = Function(W), Function(W)
pnew, pold, Hold = Function(V), Function(V), Function(V)
E_du = ((1.0-pold)**2+1e-6)*inner(grad(v),sigma(u))*dx
# form linear
E_phi = (Gc*l*inner(grad(p),grad(q))+((Gc/l)+2.0*H(uold,unew,Hold))\
*inner(p,q)-2.0*H(uold,unew,Hold)*q)*dx
p_disp = LinearVariationalProblem(lhs(E_du), rhs(E_du), unew, bc_u)
p_phi = LinearVariationalProblem(lhs(E_phi), rhs(E_phi), pnew, bc_phi)
solver_disp = LinearVariationalSolver(p_disp)
solver_phi = LinearVariationalSolver(p_phi)
# Initialization of the iterative procedure and output requests
# t = 0
u_r = 1.000
# deltaT = 0.00001 # real 0.00001
t = 0 # Start time
T = 0.007 # End time
num_steps = 20 # Number of time steps
deltaT = (T-t)/num_steps # Time step size
tol = 1e-3
#conc_f = XDMFFile("./ResultsDir/phi.xdmf")
conc_f = File ("./ResultsDir/phi.pvd")
fname = open('ForcevsDisp.txt', 'w')
load_dis =[] # load displacement
load_force =[] # load force
forcex = []
forcey = []
print("l=",l,"Gc=",Gc,"deltaT=",deltaT)
x_dofs = W.sub(0).dofmap().dofs()
y_dofs = W.sub(1).dofmap().dofs()
# Staggered scheme
#while t<=0.021: # real 0.007 test 0.00002
for i in range(num_steps):
t += deltaT
# if t >=0.005:
# deltaT = 0.000001 # 0.000001
load.t=t
iter = 0
err = 1
while err > tol:
iter += 1
solver_disp.solve()
solver_phi.solve()
err_u = errornorm(unew,uold,norm_type = 'l2',mesh = None)
err_phi = errornorm(pnew,pold,norm_type = 'l2',mesh = None)
err = max(err_u,err_phi)
uold.assign(unew)
pold.assign(pnew)
Hold.assign(project(psi(unew), WW))
if err < tol:
print ('Iterations:', iter, ', Total time', t)
if round(t*1e6) % 10 == 0:
conc_f<<pnew
Traction = dot(sigma(unew),n)
fy = Traction[0]*ds(1)
fname.write(str(t*u_r) + "\t")
fname.write(str(assemble(fy)) + "\n")
load_dis.append(t*u_r)
load_force.append(assemble(fy)/1000.0)
f_int = assemble(((1.0-pold)**2)*inner(grad(v),sigma(unew))*dx)
bctop.apply(f_int)
Fx = 0
for i in x_dofs:
Fx += f_int[i]
Fy = 0
for i in y_dofs:
Fy += f_int[i]
#print("Fx=",Fx)
#print("Fy=",Fy)
forcex.append(-Fx/1000.0)
forcey.append(-Fy/1000.0)
fname.close()
np.savetxt('./curvedata/load_dis.out',load_dis) #save load_dis
np.savetxt('./curvedata/load_force.out',load_force) #save load_force
np.savetxt('./curvedata/forcex.out',forcex)
np.savetxt('./curvedata/forcey.out',forcey)
print ('Simulation completed')
Any assistance will be greatly appreciated!
Zhihai Wang
Best regards