Hi, I tried to adapt and modify the elasto-plastic code from this amazing tutorial from @bleyerj. But I observed unexpected correlation between iterative solver and input image to extract color data from it and store it in .dat file and use that dat file to store the material information using the mesh cell midpoints.
The code works on mspaint made image like this -
But the iterator doesn’t converge for the image which I generated using python which is below (in this case) -
Here is the minimal working example to reproduce the error -
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
parameters["form_compiler"]["representation"] = "tsfc"
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
from dolfin import *
import matplotlib.pyplot as plt
if not has_linear_algebra_backend("PETSc"):
print("DOLFIN has not been configured with PETSc. Exiting.")
exit()
parameters["linear_algebra_backend"] = "PETSc"
def build_nullspace(V, x):
gdim = V.mesh().geometry().dim()
assert gdim == 2 or gdim == 3
dim = 3 if gdim == 2 else 6
nullspace_basis = [x.copy() for i in range(dim)]
for i in range(gdim):
V.sub(i).dofmap().set(nullspace_basis[i], 1.0);
if gdim == 2:
V.sub(0).set_x(nullspace_basis[2], -1.0, 1);
V.sub(1).set_x(nullspace_basis[2], 1.0, 0);
elif gdim == 3:
V.sub(0).set_x(nullspace_basis[3], -1.0, 1);
V.sub(1).set_x(nullspace_basis[3], 1.0, 0);
V.sub(0).set_x(nullspace_basis[4], 1.0, 2);
V.sub(2).set_x(nullspace_basis[4], -1.0, 0);
V.sub(2).set_x(nullspace_basis[5], 1.0, 1);
V.sub(1).set_x(nullspace_basis[5], -1.0, 2);
for x in nullspace_basis:
x.apply("insert")
return VectorSpaceBasis(nullspace_basis)
import meshio
mesh_from_file = meshio.read("mesh.msh")
import numpy
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
points = mesh.points[:,:2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells})
return out_mesh
triangle_mesh = create_mesh(mesh_from_file, "triangle", prune_z=True)
meshio.write("mesh.xdmf", triangle_mesh)
mesh = Mesh()
with XDMFFile("mesh.xdmf") as infile:
infile.read(mesh)
x_scale = 1.0 / 1023.0
y_scale = 1.0 / 1023.0
for vertex in mesh.coordinates():
vertex[0] *= x_scale
vertex[1] *= y_scale
def extract_fourth_column(dat_file):
fourth_column = []
with open(dat_file, "r") as file:
for line in file:
elements = line.strip().split()
fourth_column.append(int(elements[3]))
return fourth_column
dat_file = "mesh.dat"
fourth_column_array = extract_fourth_column(dat_file)
fourth_column_array
import dolfin
materials = MeshFunction('double', mesh, 2)
G = VectorFunctionSpace(mesh, "Lagrange", 2)
g = Function(G)
local_values_material = np.zeros_like(g.vector().get_local())
for cell in cells(mesh):
midpoint = cell.midpoint().array()
i = (midpoint[0])
j = (midpoint[1])
k = (midpoint[2])
local_values_material[cell.index()] = fourth_column_array[cell.index()]
materials[cell] = int(local_values_material[cell.index()])
g.vector().set_local(local_values_material)
dolfin.XDMFFile(dolfin.MPI.comm_world, "mesh.xdmf").write_checkpoint(g,"g",0)
import dolfin as dol
class al(dol.UserExpression):
def __init__(self, materials, al0, al1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = al0
self.k_1 = al1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
E1 = 21e3
nu1 = 0.3
E2 = 2100e3
nu2 = 0.25
E = al(materials, E1, E2, degree = 0)
nu = al(materials, nu1, nu2, degree = 0)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = Constant(250.)
Et = E/100.
H = E*Et/(E-Et)
deg_u = 2
deg_stress = 2
V = VectorFunctionSpace(mesh, "CG", deg_u)
We = VectorElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, dim=4, quad_scheme='default')
W = FunctionSpace(mesh, We)
W0e = FiniteElement("Quadrature", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')
W0 = FunctionSpace(mesh, W0e)
sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
class bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0., 1e-8)
class top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1., 1e-8)
boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
bottom().mark(boundaries, 1)
top().mark(boundaries, 2)
q_lim = float(sig0)
loading = Expression("q*t", q=q_lim, t=0, degree=2)
bc_bottom = DirichletBC(V, Constant((0., 0.)), boundaries, 1)
ds = Measure("ds")(subdomain_data=boundaries)
n = FacetNormal(mesh)
bc_u = [bc_bottom]
def F_ext(v):
return loading*dot(n, v)*ds(2)
def eps(v):
e = sym(grad(v))
return as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def sigma(eps_el):
return lmbda*tr(eps_el)*Identity(3) + 2*mu*eps_el
def as_3D_tensor(X):
return as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
ppos = lambda x: (x+abs(x))/2.
def proj_sig(deps, old_sig, old_p):
sig_n = as_3D_tensor(old_sig)
sig_elas = sig_n + sigma(deps)
s = dev(sig_elas)
sig_eq = sqrt(3/2.*inner(s, s))
f_elas = sig_eq - sig0 - H*old_p
dp = ppos(f_elas)/(3*mu+H)
n_elas = s/sig_eq*ppos(f_elas)/f_elas
beta = 3*mu*dp/sig_eq
new_sig = sig_elas-beta*s
return as_vector([new_sig[0, 0], new_sig[1, 1], new_sig[2, 2], new_sig[0, 1]]), \
as_vector([n_elas[0, 0], n_elas[1, 1], n_elas[2, 2], n_elas[0, 1]]), \
beta, dp
def sigma_tang(e):
N_elas = as_3D_tensor(n_elas)
return sigma(e) - 3*mu*(3*mu/(3*mu+H)-beta)*inner(N_elas, e)*N_elas-2*mu*beta*dev(e)
metadata = {"quadrature_degree": deg_stress, "quadrature_scheme": "default"}
dxm = dx(metadata=metadata)
a = inner(eps(v), sigma_tang(eps(u_)))*dxm
L = -inner(eps(u_), as_3D_tensor(sig))*dxm + F_ext(u_)
def local_project(v, V, u=None):
dv = TrialFunction(V)
v_ = TestFunction(V)
a_proj = inner(dv, v_)*dxm
b_proj = inner(v, v_)*dxm
solver = LocalSolver(a_proj, b_proj)
solver.factorize()
if u is None:
u = Function(V)
solver.solve_local_rhs(u)
return u
else:
solver.solve_local_rhs(u)
return
Nitermax, tol = 200, 1e-8
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
results = np.zeros((Nincr+1, 2))
for (i, t) in enumerate(load_steps):
loading.t = t
A, b = assemble_system(a, L, bc_u)
b0 = b.norm("l2")
b1 = b0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while b1/b0 > tol and niter < Nitermax:
null_space = build_nullspace(V, du.vector())
as_backend_type(A).set_near_nullspace(null_space)
pc = PETScPreconditioner("petsc_amg")
PETScOptions.set("mg_levels_ksp_type", "chebyshev")
PETScOptions.set("mg_levels_pc_type", "jacobi")
PETScOptions.set("mg_levels_esteig_ksp_type", "cg")
PETScOptions.set("mg_levels_ksp_chebyshev_esteig_steps", 50)
solver = PETScKrylovSolver("cg", pc)
solver.parameters["monitor_convergence"] = False
solver.set_operator(A);
solver.solve(du.vector(), b);
Du.assign(Du+du)
deps = eps(Du)
sig_, n_elas_, beta_, dp_ = proj_sig(deps, sig_old, p)
local_project(sig_, W, sig)
local_project(n_elas_, W, n_elas)
local_project(beta_, W0, beta)
A, b = assemble_system(a, L, bc_u)
b1 = b.norm("l2")
print(" Residual:", b1)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
p_avg.assign(project(p, P0))
results[i+1, :] = (u(0.5, 1.)[1], t)
I apologize if minimal working example is too long, but it is to the best of my knowledge to cut it short.
I have two questions -
- What is the correlation between the image and the iterative solver that it works for one image (mspaint_generated_image) and doesn’t works for another image (python_generated_image) despite of using the same RGB values for both images which is
(237, 28, 36) => red and (63, 72, 204) => blue? - How to solve this issue ?
I have attached msh and dat files for both the images shown above here.
Kindly assist me on how to tackle this error.
Thanks.