I apologize for the delay in response. I assume my error is coming from displacement control loading. I am trying to convert Elasto-plastic analysis of a 2D von Mises material — Numerical tours of continuum mechanics using FEniCS master documentation to a displacement control loading. I want to plot a load - displacement curve and validate the results. What I am trying to do here is,
Traction = sigma.n
f_y = Traction * ds
This is my code
from dolfin import *
import numpy as np
parameters["form_compiler"]["representation"] = 'quadrature'
import warnings
from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning
warnings.simplefilter("once", QuadratureRepresentationDeprecationWarning)
# material parameters
E = Constant(70e3)
nu = Constant(0.3)
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2./(1+nu)
sig0 = Constant(250.) # yield strength
Et = E/100. # tangent modulus
H = E*Et/(E-Et) # hardening modulus
mesh = Mesh("mesh.xml")
# Boundary conditions
# normal AND
# ds, load:user expression
top = CompiledSubDomain("near(x[1], 0.5) && on_boundary")
# bottom boundary
bot = CompiledSubDomain("near(x[1], -0.5) && on_boundary")
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)
#load, defined using user expression, t
load = Expression("t", t = 0.0, degree=1)
# bc, bottom
bcbot= DirichletBC(V, Constant((0.0,0.0)), bot)
# load applied on top u_y
# boundary condition on top
bctop = DirichletBC(V.sub(1), load, top)
# grouping disp bc together
bc = [bcbot, bctop]
n = FacetNormal(mesh)
sig = Function(W)
sig_old = Function(W)
n_elas = Function(W)
beta = Function(W0)
p = Function(W0, name="Cumulative plastic strain")
f_y = Function(W0, name = "vertical reaction")
u = Function(V, name="Total displacement")
du = Function(V, name="Iteration correction")
Du = Function(V, name="Current increment")
v = TrialFunction(V)
u_ = TestFunction(V)
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_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm
res = -inner(eps(u_), as_3D_tensor(sig))*dxm
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
# XDMF : eXtensible Data Model and Format
# method to exchange scientific data
# XDMF will not handle quadrature space
# we define Functionspace P0
file_results = XDMFFile("plasticity_results.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
P0 = FunctionSpace(mesh, "DG", 0)
p_avg = Function(P0, name="Plastic strain")
load_vertical = Function(P0, name="vertical reaction")
u_r = 0.05
Nitermax, tol = 5, 1e-8 # parameters of the Newton-Raphson procedure
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):
load.t = t*u_r
A, Res = assemble_system(a_Newton, res, bc)
nRes0 = Res.norm("l2")
nRes = nRes0
Du.interpolate(Constant((0, 0)))
print("Increment:", str(i+1))
niter = 0
while nRes/nRes0 > tol and niter < Nitermax:
solve(A, du.vector(), Res, "mumps")
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, Res = assemble_system(a_Newton, res, bc)
nRes = Res.norm("l2")
print(" Residual:", nRes)
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
file_results.write(u, t)
p_avg.assign(project(p, P0))
sig_yy = project(sig_[1],P0)
f_y = sigp[1] *ds(1)
#load_vertical.assign(project(f_y, P0))
file_results.write(str(assemble(f_y)), t)
results[i+1, :] = (t*u_r, t)
The error I am obtaining is
Traceback (most recent call last):
File "/home/anupama/OneDrive/WORK_2023/von_mises_july/vonMises_plasticity/vonMises_plasticity.py", line 166, in <module>
sig_yy = project(sig_[1],P0)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 132, in project
A, b = assemble_system(a, L, bcs=bcs,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 409, in assemble_system
b_dolfin_form = _create_dolfin_form(b_form, form_compiler_parameters)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 60, in _create_dolfin_form
return Form(form,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 100, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
return compile_ufl_objects(forms, "form", object_names,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 190, in compile_ufl_objects
ir = compute_ir(analysis, prefix, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 182, in compute_ir
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 182, in <listcomp>
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 455, in _compute_integral_ir
ir = r.compute_integral_ir(itg_data,
File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 77, in compute_integral_ir
tabulate_basis(sorted_integrals, form_data, itg_data)
File "/usr/lib/python3/dist-packages/ffc/quadrature/tabulate_basis.py", line 224, in tabulate_basis
psi_table = _tabulate_psi_table(integral_type, cellname, tdim,
File "/usr/lib/python3/dist-packages/ffc/quadrature/tabulate_basis.py", line 120, in _tabulate_psi_table
psi_table[entity] = element.tabulate(deriv_order, entity_points)
File "/usr/lib/python3/dist-packages/FIAT/mixed.py", line 76, in tabulate
table = e.tabulate(order, points, entity)
File "/usr/lib/python3/dist-packages/FIAT/quadrature_element.py", line 56, in tabulate
raise AssertionError("Mismatch of quadrature points!")
AssertionError: Mismatch of quadrature points!
Can anyone suggest a way to extract vertical reaction force without this quadrature mismatch coming into the picture?
(I am working with a 2D unit square domain)