Hello everyone,
I am attempting to convert the 2D orthotropic elastic tutorials into 3D, and I receive a very long error message at the point that I try to generate the problem (I get the same error message for both LinearProblem and NonLinearProblem). I tried to trim it down to a MWE, but it’s still 100 lines b/c, unfortunately, everything ultimately feeds into the LinearProblem definition. Any advice is appreciated - thank you.
dolfinx version 0.7.3
HPC system (operating on home node single-processor for now)
# attempted 3D implementation of the 2D tutorials:
# https://bleyerj.github.io/comet-fenicsx/tours/linear_problems/isotropic_orthotropic_elasticity/isotropic_orthotropic_elasticity.html
# https://comet-fenics.readthedocs.io/en/latest/demo/elasticity/orthotropic_elasticity.py.html
import ufl
import numpy as np
from petsc4py import PETSc
default_scalar_type = PETSc.ScalarType
from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc # for some reason, "import dolfinx" doesn't include this
import dolfinx.nls.petsc # or this
domain = dolfinx.mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([1, 1, 1])],
[10, 10, 10], cell_type=dolfinx.mesh.CellType.hexahedron)
Ex, Ey, nuxy, Gxy = 10.0, 10.0, 0.3, 5.0
S0 = np.array([ [ (1/Ex) , (-nuxy / Ex) , (-nuxy / Ex) , 0 , 0 , 0 ],
[ (-nuxy / Ex) , (1/Ex) , (-nuxy / Ex) , 0 , 0 , 0 ],
[ (-nuxy / Ex) , (-nuxy / Ex) , (1/Ex) , 0 , 0 , 0 ],
[ 0 , 0 , 0 , 1/Gxy , 0 , 0 ],
[ 0 , 0 , 0 , 0 , 1/Gxy , 0 ],
[ 0 , 0 , 0 , 0 , 0 , 1/Gxy ] ])
S = ufl.as_matrix(S0)
C = ufl.inv(S)
def epsilon(v):
return ufl.sym(ufl.grad(v))
def strain2voigt(e): # converts stdrain to a vector
return ufl.as_vector([ e[0,0] , e[1,1] , e[2,2] , 2*e[0,1] , 2*e[0,2] , 2*e[1,2] ])
def voigt2stress(s): # converts vector stress to tensor stress
return ufl.as_tensor([[ s[0] , s[3] , s[4] ],
[ s[3] , s[1] , s[5] ],
[ s[4] , s[5] , s[2] ] ])
def sigma(v):
P = (ufl.dot(C, strain2voigt(epsilon(v))))
return P
# Define function space
displacement_element = ufl.VectorElement ("CG", domain.ufl_cell() , 2)
V = dolfinx.fem.functionspace(domain, displacement_element)
# Define variational problem
du = ufl.TrialFunction(V) # I tried w/ and w/out TrialFunction (LinearProblem and NonLinearProblem respectively)
u_ = ufl.TestFunction(V)
u = dolfinx.fem.Function(V, name="Displacement")
a_form = ufl.inner( sigma(du), strain2voigt(epsilon(u_)) ) * ufl.dx
# a_form = ufl.inner( voigt2stress(sigma(du)), epsilon(u_) ) * ufl.dx # this is also logical - but gives the same error message
# uniform traction on top boundary
T = dolfinx.fem.Constant(domain, default_scalar_type((10)))
n = ufl.FacetNormal(domain)
coordinates = domain.geometry.x
xmin, ymin , zmin = coordinates.min(axis=0)
xmax, ymax , zmax = coordinates.max(axis=0)
def top(x):
return np.isclose(x[1], ymax)
def bottom(x):
return np.isclose(x[1], ymin)
def left(x):
return np.isclose(x[0], xmin)
def back(x):
return np.isclose(x[2], zmax)
fdim = domain.topology.dim - 1
top_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, top) # used to apply load
bottom_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, bottom) # used in bcs to fix structure
left_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, left) # used in bcs to fix structure
back_facets = dolfinx.mesh.locate_entities_boundary(domain, fdim, back) # used in bcs to fix structure
marked_facets = np.hstack([top_facets])
marked_values = np.hstack([ np.full_like(top_facets, 1)])
sorted_facets = np.argsort(marked_facets)
facet_tag = dolfinx.mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
ds = ufl.Measure ("ds", domain = domain , subdomain_data = facet_tag ) # this is used to define the traction surface
L_form = ufl.dot(T * n, u_) * (ds(1))
# for fixed BCs: ---------------------------------------
bottom_dofs = dolfinx.fem.locate_dofs_topological((V.sub(1)), fdim, bottom_facets) # my modifications to be more like Dokken's elasticity tutorial
left_dofs = dolfinx.fem.locate_dofs_topological((V.sub(0)), fdim, left_facets)
back_dofs = dolfinx.fem.locate_dofs_topological((V.sub(2)), fdim, back_facets)
bcs = [
dolfinx.fem.dirichletbc(default_scalar_type(0), left_dofs, V.sub(0)),
dolfinx.fem.dirichletbc(default_scalar_type(0), bottom_dofs, V.sub(1)),
dolfinx.fem.dirichletbc(default_scalar_type(0), back_dofs, V.sub(2)),
]
# # for linear:
problem = dolfinx.fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
problem.solve()
# for non-linear:
# problem = dolfinx.fem.petsc.NonlinearProblem ( (a_form+L_form) , u , bcs = bcs ) # this gave the same error message as linear
# solver = dolfinx.nls.petsc.NewtonSolver( domain.comm , problem )
# solver.convergence_criterion = "incremental"
# num_its , converged = solver.solve( u )
# I am happy to add the rest of the code, but it stops with the error when I define "problem", so I end here
The error message:
Traceback (most recent call last):
File "/home/axt/fenicsx-poro1/orthotropic_stuff/posting.py", line 100, in <module>
problem = dolfinx.fem.petsc.LinearProblem(a_form, L_form, u=u, bcs=bcs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/petsc.py", line 588, in __init__
self._a = _create_form(a, form_compiler_options=form_compiler_options, jit_options=jit_options)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 188, in form
return _create_form(form)
^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 183, in _create_form
return _form(form)
^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/forms.py", line 141, in _form
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/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 "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ffcx/codegeneration/jit.py", line 199, in compile_forms
raise e
File "/home/axt/miniconda3/envs/fenicsx-env/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 "/home/axt/miniconda3/envs/fenicsx-env/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 "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 88, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, options) for form in forms)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 88, in <genexpr>
form_data = tuple(_analyze_form(form, options) for form in forms)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ffcx/analysis.py", line 163, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 272, in compute_form_data
form = preprocess_form(form, complex_mode)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/compute_form_data.py", line 231, in preprocess_form
form = apply_algebra_lowering(form)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/apply_algebra_lowering.py", line 148, in apply_algebra_lowering
return map_integrand_dags(LowerCompoundAlgebra(), expr)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py", line 73, in map_integrand_dags
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py", line 25, in map_integrands
mapped_integrals = [map_integrands(function, itg, only_integral_type)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py", line 25, in <listcomp>
mapped_integrals = [map_integrands(function, itg, only_integral_type)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py", line 33, in map_integrands
return itg.reconstruct(function(itg.integrand()))
^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/map_integrands.py", line 73, in <lambda>
return map_integrands(lambda expr: map_expr_dag(function, expr, compress),
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 34, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/corealg/map_dag.py", line 100, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/algorithms/apply_algebra_lowering.py", line 104, in inverse
return inverse_expr(A)
^^^^^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/compound_expressions.py", line 156, in inverse_expr
return adj_expr(A) / determinant_expr(A)
^^^^^^^^^^^
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/ufl/compound_expressions.py", line 174, in adj_expr
raise ValueError(f"adj_expr not implemented for dimension {sh[0]}.")
ValueError: adj_expr not implemented for dimension 6.
Exception ignored in: <function LinearProblem.__del__ at 0x2ac6f33c20c0>
Traceback (most recent call last):
File "/home/axt/miniconda3/envs/fenicsx-env/lib/python3.11/site-packages/dolfinx/fem/petsc.py", line 627, in __del__
self._solver.destroy()
^^^^^^^^^^^^
AttributeError: 'LinearProblem' object has no attribute '_solver'