Can anyone help me understand this warning? (I am trying to use quadrature function space)
Ignoring precision in integral metadata compiled using quadrature representation. Not implemented
Without any code to reproduce the error and a version specification of what dolfin you are using and how you installed it, it is hard to Give Amy guidance
I apologize for the delay. ( I am trying to convert the fenics tutorial on plasticity to displacement control problem) This is my code, in addition to this warning my code is diverging if the time step is not too small.
from dolfin import *
from matplotlib import pyplot
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(180e3)
nu = Constant(0.28)
lmbda = Constant(89633)
mu = Constant(70300)
sig0 = Constant(443) # yield strength
H = Constant(300) # hardening modulus
mesh = UnitSquareMesh(200,200)
top = CompiledSubDomain("near(x[1], 1.0) && on_boundary")
bot = CompiledSubDomain("near(x[1], 0.0) && 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 = Expression("t", t = 0.0, degree=1)
bcbot= DirichletBC(V, Constant((0.0,0.0)), bot)
bctop = DirichletBC(V.sub(1), load, top)
bc_u = [bcbot,bctop]
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
top.mark(boundaries, 1)
ds = Measure("ds")(subdomain_data= boundaries)
n = FacetNormal(mesh)
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)
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)
# --------------------------------------------
# Global problem and Newton-Raphson procedure
# --------------------------------------------
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
u_r = 0.03
print("\n Loading: u_r ", float(u_r))
Nitermax, tol = 100, 1e-8 # parameters of the Newton-Raphson procedure
Nincr = 20
load_steps = np.linspace(0, 1.1, Nincr+1)[1:]**0.5
print("time step", Nincr)
for (i, t) in enumerate(load_steps):
load.t = t*u_r
A, Res = assemble_system(a_Newton, res, bc_u)
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_u)
nRes = Res.norm("l2")
print(" Residual:", nRes)
for bc in bc_u:
bc.homogenize()
niter += 1
u.assign(u+Du)
p.assign(p+local_project(dp_, W0))
sig_old.assign(sig)
The code is using old compilation functionality with
which is missing some features.
I wouldn’t worry about the warning if you want to use DOLFIN with quadrature elements.
In general I would strongly suggest to use DOLFINx, see: The new DOLFINx solver is now recommended over DOLFIN
See for instance @bleyerj tutorial on: Welcome — Computational Mechanics Numerical Tours with FEniCSx for various topics in solid mechanics
1 Like
Thank you for the suggestion.