I am a beginner in FEnics and trying to set up a time-dependant nonlinear problem similar to the Cahn-Hilliard example and am running into some issues involving Dirichlet boundary conditions and possibly convergence. I am able to set up and change Neumann boundary conditions.
For some reason it appears that my displacement boundary conditions are applied relatively, i.e. even though I define them outside the time loop the boundary moves at every timestep.
from dolfin import *
from ufl import replace
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import brentq
class TriangularLattice(NonlinearProblem):
def __init__(self,L,a,bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
for bc in bcs:
bc.apply(b)
def J(self, A, x):
assemble(self.a, tensor = A)
for bc in bcs:
bc.apply(A)
''' -------------------------
Constants
------------------------- '''
set_log_level(20)
length, H = 1, .1
meshlength = 100*length
meshheight = int(200*H)
sigc = 0# imposed Traction
beta = 1e2 #dissapation constant (includes dependence on time step)
Tend = .2
# Nsteps=int(Tend*50)
Nsteps = 3
dt=Tend/(Nsteps+1)
#Constants for defining energy
eps_s=0.25
phi_s=0.1
b0=60
b1=80
R1=Constant((1,0))
R2=Constant((0.5,sqrt(3)/2))
R3=Constant((-0.5,sqrt(3)/2))
mesh = RectangleMesh(Point(0., 0.), Point(length, H), meshlength, meshheight)
''' -------------------------
BOUNDARY CONDITIONS
------------------------- '''
def left(x, on_boundary):
return near(x[0],0.)
def right(x, on_boundary):
return near(x[0],length)
''' -------------------------
SET UP VARIATIONAL PROBLEM
------------------------- '''
#Interpolation for displacement
W = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
# Initial displacements:
u_0 = Expression(("0.25*x[0]", "0.25*x[1]") ,degree=1)
# u_0 = Expression(("0.1*x[0]", "0.1*x[1]") ,degree=1)
u = interpolate(u_0,W)
u_old = interpolate(u_0,W)
#w, w_old, dw, and w_ are all UNKNOWN functions within the function space W
u_ = TestFunction(W)
du = TrialFunction(W)
Here I define the right side boundary and when I run the whole script the boundary moves by 1e-3 at every time step. I would like to be able to define the boundaries as absolute, i.e. in the reference coordinate system. If I increase the displacement (here, 1e-3), the solver no longer converges – I think this is a related issue but I am not sure. Any help would be much appreciated.
''' -------------------------
DIRICHLET BOUNDARY CONDITIONS
------------------------- '''
#Left and right boundary conditions
u_right = Expression(("0.001*x[0]", "0*x[1]") ,degree=1)
bcl = DirichletBC(W, Constant((0,0)), left)
bcr = DirichletBC(W, Constant((1e-3,0)), right)
bcs = [bcl,bcr]
''' -------------------------
FUNCTIONS TO DEFINE MATERIAL MODEL
------------------------- '''
#strain as a function of displacement
def eps(disp):
return sym(grad(disp))
#strain rate as a function of displacement
def eps_rate(disp,disp_old):
return (eps(disp)-eps(disp_old))/dt
def eps_dev(disp):
return eps(disp)-(tr(eps(disp))/2)*Identity(2)
# #stress as a function of displacement
def sigma(disp,disp_old):
sig_inter = dbistable_energy(dot(R1,dot(eps(disp),R1)))*outer(R1,R1)+dbistable_energy(dot(R2,dot(eps(disp),R2)))*outer(R2,R2)+dbistable_energy(dot(R3,dot(eps(disp),R3)))*outer(R3,R3)
sig_visc = 0.75*beta*eps_rate(disp,disp_old)+0.375*beta*tr(eps_rate(disp,disp_old))*Identity(2)
return sig_inter + sig_visc
# #strain energy as a function of volumetric strain
def psi(eps_vol):
def find0(a):
return 0.5*b1*(a-eps_s)**2+phi_s - 0.5*b0*a**2
root = Constant(brentq(find0, 0, eps_s))
eval = conditional(le(eps_vol, root), 0.5*b0*eps_vol**2, 0.5*b1*(eps_vol-eps_s)**2+phi_s)
return eval
def dpsi(eps_vol):
def find0(a):
return 0.5*b1*(a-eps_s)**2+phi_s - 0.5*b0*a**2
root = Constant(brentq(find0, 0, eps_s))
eval = conditional(le(eps_vol, root), b0*eps_vol, b1*(eps_vol-eps_s))
return eval
def ddpsi(eps_vol):
def find0(a):
return 0.5*b1*(a-eps_s)**2+phi_s - 0.5*b0*a**2
root = Constant(brentq(find0, 0, eps_s))
eval = conditional(le(eps_vol, root), b0, b1)
return eval
def bistable_energy(eps_vol):
return 3*psi(eps_vol)
def dbistable_energy(eps_vol):
return 3*dpsi(eps_vol)
def tr_cross(a,b):
return a[0,0]*b[0,0]+a[0,1]*b[1,0]+a[1,0]*b[0,1]+a[1,1]*b[1,1]
def deviatoric_energy(eps_vol,eps_dev):
1.5*ddpsi(eps_vol)*tr(eps_dev)+0.375*ddpsi(eps_vol)*tr_cross(eps_dev,eps_dev)
def incremental_energy(epsdot):
3/16 *(2*tr_cross(epsdot,epsdot)+pow(tr(epsdot),2))
def total_energy(eps_vol,eps_dev,epsdot):
dev_energy = 1.5*ddpsi(eps_vol)*tr(eps_dev)+0.375*ddpsi(eps_vol)*tr_cross(eps_dev,eps_dev)
incr_energy = 3/16 *(2*tr_cross(epsdot,epsdot)+pow(tr(epsdot),2))
return bistable_energy(eps_vol)+dev_energy+incr_energy
#dissapation potential as a function of viscous strain
def dissipation_potential(u, u_old):
return 0.5*beta*inner(eps(u)-eps(u_old),eps(u)-eps(u_old))
#Total and viscous strain are minimization to the following incremental potential
incremental_potential = total_energy(tr(eps(u))/2,eps_dev(u),eps_rate(u,u_old))*dx + dt*dissipation_potential(u,u_old)*dx
#Find varational solution by setting derivative == 0
L = derivative(incremental_potential, u, u_)
#Calculate Jacobian
a = derivative(L, u, du)
problem = TriangularLattice(L,a,bcs)
# Time-stepping loop
time = np.linspace(0, Tend, Nsteps+1)
''' -------------------------
TIME STEPPING LOOP
------------------------- '''
for (i, ti) in enumerate(time): #For each time step
print("Increment {}/{}".format(i, Nsteps))
u_old.vector()[:] = u.vector()
''' NEWTON SOLVER '''
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 8e-2
solver.parameters["absolute_tolerance"] = 8e-2
solver.parameters["maximum_iterations"] = 100
solver.solve(problem,u.vector())