In solid mechanics it is usual to define a routine as a material model and usually in non elastic models some history variables should be kept. In order to handle this I defined a function for keeping the history variable in the system. The following code shows what I did:
--- Define variational problem
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = TrialFunction(V) # Displacement at time n+1
u_n= Function(V) # Displacement at time n
Du = Function(V) # Function for u_(n+1)=u_n - Du
alpha = Function(H2) # History variables at time n+1
alpha_n= Function(H2) # History variables at time n
alpha_n = project(Constant( ((0,0),(0,0)) ), H2)
Kinematics and material model
epsilon = sym(grad(u))
Material model routine calculates sigma_(n+1) and alpha_(n+1)
mat = material_model(epsilon,alpha_n,parameters)
sigma = mat['sigma']
alpha = mat['alpha']
Bilinear form definitions
L = inner(Sigma, nabla_grad(v))dx
F = dot(Tr, v)ds(1)
B = L - F
Jacobian is int (B^T j B) but for now I use numerical jacobian !!
J = derivative(B, u, du)
Now having the above form I want to integrate over time the equation and in each time step I should solve a nonlinear equation
--- Compute time dependent solution
T = 0.22 # total simulation time
dt= 0.001
t = dt # initial time step
while t <= T:
Tr.t = t
omega = 1.0 # relaxation parameter
eps = 1.0
tol = 1.0E-5
iter = 0
maxiter = 25
while eps > tol and iter < maxiter:
iter += 1
j, b = assemble_system(J, B, bc)
solve(j, Du.vector(), b)
eps = numpy.linalg.norm(Du.vector().array(), ord=numpy.Inf)
print 'Norm:', eps
u.vector()[:] = u_n.vector() - omega*Du.vector()
u_n.assign(u)
alpha_n.assign(mat['alpha'])
t += dt
The problem is that I do not know how to make sure that the line
alpha_n.assign(mat['alpha'])
works. I do not call material routine directly in the iterations and assemble_system function does this behind the scene, so I do not have a direct access to alpha here.
Any help on this issue is very welcomed.
Regards
