Hey y’all. Just as a preface, I am using v 0.10 of DOLFINx. I am building a axisymmetric contact mechanics model between and analytical rigid indenter and an elastic body using a penalty method, currently using the linear elasticity constitutive laws. I have defined the strain tensor as below:
def eps_u(u):
return sym(as_tensor([
[u[0].dx(0), 0, u[0].dx(1)],
[0, u[0]/ r_mesh, 0],
[u[1].dx(0), 0, u[1].dx(1)]
]))
The use of the strain tensor above allows for the FEA to run perfectly, but produces an asymptotic artifact in the volumetric strain graph as the radial distance approaches 0. Through online searching, I am of the belief that this issue is caused by the middle hoop strain term. To accommodate this I was planning on using the modified strain tensor below:
def eps_u(u):
hoop = conditional(lt(r_mesh, 1e-8), u[0].dx(0), u[0] / r_mesh)
return sym(as_tensor([
[u[0].dx(0), 0, u[0].dx(1)],
[0, hoop, 0],
[u[1].dx(0), 0, u[1].dx(1)]
]))
But this returns this error:
Traceback (most recent call last):
File "/home/user/project.py", line 198, in <module>
problem = NonlinearProblem(form, u_displacement, bcs=bcs, J=J, petsc_options_prefix="my_problem_")
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/fem/petsc.py", line 1239, in __init__
self._F = _create_form(
~~~~~~~~~~~~^
F,
^^
...<2 lines>...
entity_maps=entity_maps,
^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 449, in form
return _create_form(form)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 441, in _create_form
return _form(form)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/fem/forms.py", line 361, in _form
ufcx_form, module, code = jit.ffcx_jit(
~~~~~~~~~~~~^
msh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/jit.py", line 60, in mpi_jit
return local_jit(*args, **kwargs)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/jit.py", line 215, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 244, in compile_forms
raise e
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 224, in compile_forms
impl = _compile_objects(
decl,
...<9 lines>...
visualise=visualise,
)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/codegeneration/jit.py", line 349, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^
ufl_objects, prefix=module_name, options=options, visualise=visualise
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
)
^
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/compiler.py", line 113, in compile_ufl_objects
ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/ir/representation.py", line 151, in compute_ir
_compute_integral_ir(
~~~~~~~~~~~~~~~~~~~~^
fd,
^^^
...<4 lines>...
visualise,
^^^^^^^^^^
)
^
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/ir/representation.py", line 402, in _compute_integral_ir
integral_ir = compute_integral_ir(
itg_data.domain.ufl_cell(),
...<5 lines>...
visualise,
)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/ir/integral.py", line 227, in compute_integral_ir
F = compute_argument_factorization(S, rank)
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/ffcx/ir/analysis/factorization.py", line 337, in compute_argument_factorization
F.add_edge(i, F.e2i[o])
~~~~~^^^
KeyError: Zero((), (), ())
Exception ignored while calling deallocator <function NonlinearProblem.__del__ at 0x71e355d33060>:
Traceback (most recent call last):
File "/home/user/anaconda3/envs/fenicsx/lib/python3.14/site-packages/dolfinx/fem/petsc.py", line 1361, in __del__
lambda obj: obj is not None, (self._snes, self._A, self._b, self._x, self._P_mat)
AttributeError: 'NonlinearProblem' object has no attribute '_snes'
I have attached how I have currently written my NonlinearProblem, and all of the elements which go into it along with the solver below:
element_displacement = element("Lagrange", domain.basix_cell(), 2, shape = (2,))
functionspace_displacement = functionspace(domain, element_displacement)
u_displacement = Function(functionspace_displacement, name="displacement")
du_displacement = TrialFunction(functionspace_displacement)
v_displacement = TestFunction(functionspace_displacement)
x = SpatialCoordinate(domain)
r_mesh = x[0]
z_mesh = x[1]
def eps_u(u):
hoop = conditional(lt(r_mesh, 1e-8), u[0].dx(0), u[0] / r_mesh)
return sym(as_tensor([
[u[0].dx(0), 0, u[0].dx(1)],
[0, hoop, 0],
[u[1].dx(0), 0, u[1].dx(1)]
]))
def sig_u(u):
term_1 = lame_lambda*(tr(eps_u(u)))*Identity(3)
term_2 = 2*lame_mu*(eps_u(u))
return term_1+term_2
def ppos(g):
return 0.5*(g + sqrt(g*g + eps_c*eps_c))
delta_val = Constant(domain, 0.0)
n = as_vector((0.0, 1.0, 0.0)) #Normal vector in the z direction.
obstacle = -delta_val + r_mesh**2 / (2.0 * R_sphere)
gap = obstacle - u_displacement[1]
gap_ppos = ppos(-gap)
penalty_factor = scaling_factor * youngs_mod / lc_fine
ds = Measure('ds', domain=domain, subdomain_data=facet_tags)
contact_energy = 0.5 * penalty_factor * gap_ppos**2 * r_mesh * ds(3)
elastic_form = inner(sig_u(u_displacement), eps_u(v_displacement)) * r_mesh * dx
penalty_form = derivative(contact_energy, u_displacement, v_displacement)
f = Constant(domain, np.array([0.0, 0.0]))
L = inner(f, v_displacement) * r_mesh * dx
form = elastic_form + penalty_form - L
J = derivative(form, u_displacement, du_displacement)
problem = NonlinearProblem(form, u_displacement, bcs=bcs, J=J, petsc_options_prefix="my_problem_")
snes_solver = problem.solver
I apologize that I couldn’t provide a full MWE, but does anyone have an idea what could be causing this issue where using the conditional prevents initialization of the problem, and without using conditional, the file runs and solves the problem? Thank you for you time and any help!