Thanks Emek for the idea of posting mwe and also your on going support.
So far in current simple code the analysis resets on each loop iteration when the domain evolves and subdomain is updated.
So my objective is that the fsm should remember what happened in the current solution subdomain e.g. history variables of the plastic strains and deformations etc. so that when next solution domain is activated e.g. quarter subdomain then the material behaves as it is already loaded when only half domain was active. Here it is very simple example but in my actual application even the loading on previously activated subdomain will also vary in future loading steps and certain parts of domain will have to behave in a certain way depending if they were already plasticaly deformed or not.
Following is the simple testing code from fsm application where I am trying to sequentially evolve the domain of unit cube in x direction from half, quarter to whole cubic domain and solving with prescribed displacement boundary condition on top (y=1):
from dolfin import *
import fsm
#set_log_level(10)
class DirichletBoundaryX(SubDomain):
def inside(self, x, on_boundary):
return (x[0] < DOLFIN_EPS)
class DirichletBoundaryY(SubDomain):
def inside(self, x, on_boundary):
return (x[1] < DOLFIN_EPS)
class DirichletBoundaryZ(SubDomain):
def inside(self, x, on_boundary):
return (x[2] < DOLFIN_EPS)
class PrescribedDisplacementY(SubDomain):
def inside(self, x, on_boundary):
return between(x[1], (1-DOLFIN_EPS, 1) )
class HalfDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (0, 0.5))
class QuarterDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (0, 0.75))
class WholeDomain(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (0, 1))
# making output file
file1 = XDMFFile('output.xdmf')
# Making complete mesh
mesh_Complete = UnitCubeMesh(8,8,8);
mark_sub_domain = MeshFunction("size_t", mesh_Complete, mesh_Complete.topology().dim(), 0)
number_of_submesh = 3
for submesh_number in range(number_of_submesh):
if (submesh_number == 0) :
mark_domain = HalfDomain()
print('HalfDomain')
elif (submesh_number == 1) :
mark_domain = QuarterDomain()
print('QuarterDomain')
elif (submesh_number == 2) :
mark_domain = WholeDomain()
print('WholeDomain')
#mark_domain = WholeDomain()
# marking the subdomain
mark_sub_domain.set_all(0)
mark_domain.mark(mark_sub_domain, 1)
mesh = SubMesh(mesh_Complete, mark_sub_domain, 1)
E = 20000.0;
nu = 0.3;
scheme = "default"
degree = 3
dx = Measure("dx")
dx = dx(degree=degree, scheme=scheme)
V = VectorFunctionSpace(mesh, "Lagrange", 2)
element_t = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=36, quad_scheme=scheme)
Vt = FunctionSpace(mesh, element_t)
element_s = VectorElement("Quadrature", mesh.ufl_cell(), degree=3, dim=6, quad_scheme=scheme)
Vs = FunctionSpace(mesh, element_s)
zero = Constant(0.0)
prescribed_displacement = Constant(1e-3*(submesh_number+1))
bc0 = DirichletBC(V.sub(0), zero, DirichletBoundaryY(), method="pointwise")
bc1 = DirichletBC(V.sub(1), zero, DirichletBoundaryY(), method="pointwise")
bc2 = DirichletBC(V.sub(2), zero, DirichletBoundaryY(), method="pointwise")
bc3 = DirichletBC(V.sub(1), prescribed_displacement, PrescribedDisplacementY(), method="pointwise")
bcs = [bc0, bc1, bc2, bc3]
E_t = 0.3*E
hardening_parameter = E_t/(1.0 - E_t/E)
yield_stress = 9.0
u = Function(V, name="u")
def eps(u):
return as_vector([u[i].dx(i) for i in range(3)] + [u[i].dx(j) + u[j].dx(i) for i, j in [(0, 1), (0, 2), (1, 2)]])
def sigma(s):
#s = ss.function_space()
return as_matrix([[s[0], s[3], s[4]], [s[3], s[1], s[5]], [s[4], s[5], s[2]]])
def tangent(t):
#t = tt.function_space()
return as_matrix([[t[i*6 + j] for j in range(6)] for i in range(6)])
J2 = fsm.python.cpp.plasticity_model.VonMises(E, nu, yield_stress, hardening_parameter)
Qdef = fsm.UFLQuadratureFunction(eps(u), element_s, mesh)
fsm_constitutive_update = fsm.ConstitutiveUpdate(Qdef, J2)
#fsm_tangent = QuadratureFunction(mesh, Vt.element(), fsm_constitutive_update, fsm_constitutive_update.w_tangent())
#fsm_stress = QuadratureFunction(mesh, Vs.element(), fsm_constitutive_update.w_stress())
fsm_tangent = fsm.QuadratureFunction(Vt, fsm_constitutive_update.w_tangent(), fsm_constitutive_update)
fsm_stress = fsm.QuadratureFunction(Vs, fsm_constitutive_update.w_stress())
v = TestFunction(V)
uTrial = TrialFunction(V)
a = inner(eps(v), dot(tangent(fsm_tangent), eps(uTrial)) )*dx
L = inner(grad(v), sigma(fsm_stress))*dx
nonlinear_problem = fsm.PlasticityProblem(a, L, u, fsm_tangent, fsm_stress, bcs)
nonlinear_solver = NewtonSolver()
nonlinear_solver.parameters["convergence_criterion"] = "incremental";
nonlinear_solver.parameters["maximum_iterations"] = 50;
nonlinear_solver.parameters["relative_tolerance"] = 1.0e-6;
nonlinear_solver.parameters["absolute_tolerance"] = 1.0e-15;
# File names for output
eps_p_eq = fsm_constitutive_update.eps_p_eq()
#fsm_constitutive_update.eps_p_eq().compute_mean(eps_eq);
element_eps_p_eq_project = FiniteElement("CG", mesh.ufl_cell(), degree=1)
V_eps_p_eq_project = FunctionSpace(mesh, element_eps_p_eq_project)
eps_p_eq_project = Function(V_eps_p_eq_project, name="eps_p_eq")
# Solve non-linear problem
nonlinear_solver.solve(nonlinear_problem.cpp_object(), u.vector());
# Update variables
fsm_constitutive_update.update_history();
# Write output to files
file1.write(u, t=float(submesh_number));
eps_p_eq_project.vector()[:] = project(eps_p_eq, V_eps_p_eq_project, solver_type='gmres',
form_compiler_parameters={
"representation": parameters["form_compiler"]["representation"],
"quadrature_scheme": scheme,
"quadrature_degree": degree
}
).vector()
file1.write(eps_p_eq_project, t=float(submesh_number));
file1.close()
