Hello.
The following code runs with version 0.7.0 of DOLFINx and fails with version 0.8.0.
What do I need to make it run in version 0.8.0?. Thanks in advance for your kind answer.
(Continuing the discussion from https://fenicsproject.discourse.group/t/12474)
MWE
I’m sorry that this is long (I did reduce the original code by about 100 lines).
#!/usr/bin/env python
# This file uses the code from
# https://newfrac.github.io/fenicsx-fracture/notebooks/plasticity/plasticity.html
# which holds an MIT license, under the Copyright (c) 2021 Jack S. Hale.
# The listed authors are Gaspard Blondet, Andrey Latyshev, Corrado Maurini
#
# Contents were removed and adapted to version 0.8.0 of DOLFINx from the original
import dolfinx
from dolfinx import (
fem, geometry, default_scalar_type as num_type)
from dolfinx.fem import petsc
import dolfinx.fem.petsc
from dolfinx.io import gmshio
import numpy as np
from petsc4py import PETSc
from mpi4py import MPI
import ufl
import basix
import gmsh
from dolfinx.cpp.log import LogLevel, log
def strain(v):
e = ufl.sym(ufl.grad(v))
return ufl.as_tensor([[e[0, 0], e[0, 1], 0],
[e[0, 1], e[1, 1], 0],
[0, 0, 0]])
def sigma_tr(strain_el, lmbda, mu):
return (1 / 3 * (3 * lmbda + 2 * mu)
* ufl.tr(strain_el) * ufl.Identity(3))
def sigma_dev(strain_el, mu):
return 2 * mu * ufl.dev(strain_el)
def as_3D_tensor(X):
return ufl.as_tensor([[X[0], X[3], 0],
[X[3], X[1], 0],
[0, 0, X[2]]])
def tensor_to_vector(X):
'''
Take a 3x3 tensor and return a vector of size 4 in 2D
'''
return ufl.as_vector([X[0, 0], X[1, 1], X[2, 2], X[0, 1]])
def normVM(sig):
"""Von Mises' equivalent stress"""
s_ = ufl.dev(sig)
return ufl.sqrt(3 / 2 * ufl.inner(s_, s_))
def compute_new_state(du, sig_old, p_old, lmbda, mu, H, sig0):
'''This function returns the updated mechanical state
for a given displacement increment. We separate spheric
and deviatoric parts of the stress to optimize
convergence of the solver
'''
sig_n = as_3D_tensor(sig_old)
sig_el_tr = (1 / 3 * ufl.tr(sig_n) * ufl.Identity(3)
+ sigma_tr(strain(du), lmbda, mu))
sig_el_dev = ufl.dev(sig_n) + sigma_dev(strain(du), mu)
sig_el = sig_el_tr + sig_el_dev
criterion = normVM(sig_el) - sig0 - H * p_old
dp_ = ufl.conditional(
criterion < 0, 0, criterion / (3 * mu + H))
direction = ufl.dev(sig_n)/normVM(sig_n)
new_sig_tr = sig_el_tr
new_sig_dev = ufl.conditional(
criterion < 0, sig_el_dev,
sig_el_dev - 2 * mu * 3 / 2 * dp_ * direction)
return new_sig_tr, new_sig_dev, dp_
# Geometric parameters
# external/internal radius
R_e = 1.3
R_i = 1
# Mesh parameters
lc = 0.03
gdim = 2
verbosity = 0
# Quarter of a ring (pressure vessel)
# Mesh using gmsh ########################################
mesh_comm = MPI.COMM_WORLD
model_rank = 0
facet_tags = {"Lx": 1, "Ly": 2, "inner": 3, "outer": 4}
cell_tags = {"all": 20}
if mesh_comm.rank == model_rank:
gmsh.initialize()
model = gmsh.model()
model.add("Quart_cylinder")
model.setCurrent("Quart_cylinder")
# Create the points
pix = model.occ.addPoint(R_i, 0.0, 0, lc)
pex = model.occ.addPoint(R_e, 0, 0, lc)
piy = model.occ.addPoint(0, R_i, 0, lc)
pey = model.occ.addPoint(0, R_e, 0, lc)
center = model.occ.addPoint(0, 0, 0, lc)
# Create the lines
lx = model.occ.addLine(pix, pex, tag = facet_tags["Lx"])
lout = model.occ.addCircleArc(pex, center, pey, tag = facet_tags["outer"])
ly = model.occ.addLine(pey, piy, tag = facet_tags["Ly"])
lin = model.occ.addCircleArc(piy, center, pix, tag = facet_tags["inner"])
# Create the surface
cloop1 = model.occ.addCurveLoop([lx, lout, ly, lin])
surface_1 = model.occ.addPlaneSurface([cloop1], tag = cell_tags["all"])
model.occ.synchronize()
# Assign mesh and facet tags
surface_entities = [entity[1] for entity in model.getEntities(2)]
model.addPhysicalGroup(2, surface_entities, tag=cell_tags["all"])
model.setPhysicalName(2, 2, "Quart_cylinder surface")
for (key, value) in facet_tags.items():
# 1 : it is the dimension of the object (here a curve)
model.addPhysicalGroup(1, [value], tag=value)
model.setPhysicalName(1, value, key)
# Finalize mesh
model.occ.synchronize()
gmsh.option.setNumber('General.Verbosity', verbosity)
model.mesh.generate(gdim)
# import the mesh in fenicsx with gmshio
msh, cell_tags, facet_tags = gmshio.model_to_mesh(
model, mesh_comm, 0, gdim=gdim)
msh.topology.create_connectivity(msh.topology.dim - 1, msh.topology.dim)
if mesh_comm.rank == model_rank:
gmsh.finalize()
# mesh ends ###########################################
# elasticity modulus MPa
E = fem.Constant(msh, num_type(1))
# Poisson's ratio
nu = fem.Constant(msh, num_type(0.3))
lmbda = fem.Constant(
msh, num_type(E * nu / (1 + nu) / (1 - 2 * nu)))
mu = fem.Constant(msh, num_type(E / (2 * (1 + nu))))
# yield strength MPa
sig0 = fem.Constant(msh, num_type(250 / 70e3))
# hardening modulus MPa
H = fem.Constant(msh, num_type(1 / 99))
# Polynomials
deg_u = 2
deg_stress = 2
# Changed for compatibility #####################
if "0.7.0" in dolfinx.__version__:
Ve = ufl.VectorElement(
'Lagrange', msh.ufl_cell(), degree=deg_u, dim=2)
V = fem.FunctionSpace(msh, Ve)
# Vector 1x4 (not 2x2 tensor; by design: out-of-plane plasticity)
We = ufl.VectorElement(
"Quadrature", msh.ufl_cell(), degree=deg_stress, dim=4)
W = fem.FunctionSpace(msh, We)
W_scal_e = ufl.FiniteElement(
"Quadrature", msh.ufl_cell(), degree=deg_stress)
W_scal = fem.FunctionSpace(msh, W_scal_e)
elif "0.8.0" in dolfinx.__version__:
msh_cell_type = msh.basix_cell()
Ve = basix.ufl.element(
'Lagrange', msh_cell_type, degree=deg_u, shape=(gdim,))
V = fem.functionspace(msh, Ve)
# Vector 1x4 (not 2x2 tensor; by design: out-of-plane plasticity)
We = basix.ufl.quadrature_element(
msh_cell_type, value_shape=(gdim * gdim,),
degree=deg_stress)
W = fem.functionspace(msh, We)
W_scal_e = basix.ufl.quadrature_element(
msh_cell_type, value_shape=(), degree=deg_stress)
W_scal = fem.functionspace(msh, W_scal_e)
sig = fem.Function(W, name="Stress")
p = fem.Function(W_scal, name="Cumulative_plastic_strain")
u = fem.Function(V, name="Total_displacement")
du = fem.Function(V, name="Current_increment")
v = ufl.TrialFunction(V)
u_ = ufl.TestFunction(V)
dx = ufl.Measure(
"dx", domain=msh, metadata={"quadrature_degree": deg_u})
dx_m = ufl.Measure(
"dx", domain=msh, metadata={"quadrature_degree": deg_stress})
new_sig_tr, new_sig_dev, dp_ = compute_new_state(
du, sig, p, lmbda, mu, H, sig0)
residual_u = (ufl.inner(new_sig_tr, strain(v)) * dx_m
+ ufl.inner(new_sig_dev, strain(v)) * dx)
# Works
J_u = ufl.derivative(residual_u, du, u_)
# Does not work with 0.8.0
fem.form(residual_u)
fem.form(J_u)
Error (Not Implemented)
Traceback (most recent call last):
File "<string>", line 19, in __PYTHON_EL_eval
File "/tmp/pyXESJz9", line 203, in <module>
File "/usr/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 188, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
return _form(form)
^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 141, in _form
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
raise e
File "/usr/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 190, in compile_forms
impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 260, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/analysis.py", line 86, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, options) for form in forms)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/analysis.py", line 86, in <genexpr>
form_data = tuple(_analyze_form(form, options) for form in forms)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ffcx/analysis.py", line 160, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 281, in compute_form_data
form = attach_estimated_degrees(form)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 207, in attach_estimated_degrees
degree = estimate_total_polynomial_degree(integral.integrand())
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ufl/algorithms/estimate_degrees.py", line 380, in estimate_total_polynomial_degree
degrees = map_expr_dags(de, [e])
^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 98, in map_expr_dags
r = handlers[v._ufl_typecode_](v)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/ufl/algorithms/estimate_degrees.py", line 87, in coefficient
d = e.embedded_superdegree # FIXME: Use component to improve accuracy for mixed elements
^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/basix/ufl.py", line 132, in embedded_superdegree
return self.highest_degree
^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/basix/ufl.py", line 1204, in highest_degree
return self.sub_element.highest_degree
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/usr/lib/python3.11/site-packages/basix/ufl.py", line 1429, in highest_degree
raise NotImplementedError()
NotImplementedError
Ancillary
I monkey-patched my running ffcx/analysis.py with this[1], but it didn’t work:
--- a/ffcx/ffcx/analysis.py 2023-10-16
+++ a/ffcx/ffcx/analysis.py 2023-10-28
@@ -193,6 +193,16 @@
for i, integral in enumerate(integral_data.integrals):
metadata = integral.metadata()
+ # If form contains a quadrature element, use the custom quadrature scheme
+ custom_q = None
+ for e in ufl.algorithms.extract_elements(integral):
+ if e.has_custom_quadrature:
+ if custom_q is None:
+ custom_q = e.custom_quadrature()
+ else:
+ assert np.allclose(e._points, custom_q[0])
+ assert np.allclose(e._weights, custom_q[1])
+
if custom_q is None:
# Extract quadrature degree
qd_metadata = integral.metadata().get("quadrature_degree", qd_default)
My system
- FEniCSx software
dolfinx: 0.8.0.dev0_r27568.a34b0c9-1
basix: 0.8.0.dev0_r952.eebbc52-1
ufl: 2023.3.0.dev0_r3568.f132188-1
ffcx: 0.8.0.dev0_r7083.b4cb217-1
- Dependencies
python: 3.11.5-2
python-numpy: 1.26.0-1
petsc: 3.19.5.1171.g37df9106526-1
hdf5-openmpi: 1.14.2-1
boost: 1.83.0-2
adios2: 2.8.3-5
scotch: 7.0.4-1
pybind11: 2.11.1-1
python-build: 1.0.1-1
python-cffi: 1.15.1-4
python-cppimport: 22.08.02.r6.g0849d17-1
python-installer: 0.7.0-3
python-mpi4py: 3.1.4-3
python-meshio: 5.3.4-2
python-pytest: 7.4.2-1
python-scikit-build: 0.17.6-1
python-setuptools: 1:68.0.0-1
python-wheel: 0.40.0-3
xtensor: 0.24.0-2
xtensor-blas: 0.20.0-1
- Operating system
Arch Linux: 6.5.4-arch2-1
Footnotes
For reference, ffcx/analysis.py
For reference, ffcx/analysis.py
1c993ce4b1c41d8004b4429cb3b76de42eceeffa
: changes quadrature rule (source)65d96e1b2bcee8c68abbbb6a60ffd7599e316ab2
: uses element without convert ↩︎