Hello all,
Thanks in advance for considering my question. I am currently trying to set-up a basic thermo-elastic analysis with Dolphinx. As you might have understood, I am a real beginner, so I apologize in advance if my problem is trivial.
I am using the following code to simulate a holed plate submitted to a constant Temperature field with two boundaries clamped. I have built it from several examples.
import gmsh
from mpi4py import MPI
import sys
from dolfinx.io.gmshio import model_to_mesh
import ufl
from dolfinx.fem.petsc import LinearProblem
import pyvista as pv
# Geometry and Mesh
gmsh.initialize(sys.argv)
gmsh.model.add("boolean")
model_rank = 0
mesh_comm = MPI.COMM_WORLD
gdim = 2
L = 100.
l = 25.
R = l/3./2.
plate = gmsh.model.occ.addRectangle(0, 0, 0, L, l)
hole = gmsh.model.occ.addDisk(L/2, l/2, 0, R, R)
holed_plate = gmsh.model.occ.cut([(2, plate)],[(2, hole)])
gmsh.model.occ.synchronize()
surf = gmsh.model.getEntities(2)
gmsh.model.addPhysicalGroup(surf[0][0], [surf[0][1]], 3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 3)
gmsh.model.mesh.generate(gdim)
if(False):
gmsh.fltk.run()
domain, cell_markers, facet_markers = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()
from dolfinx import fem, plot, default_scalar_type, mesh
import numpy as np
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
#V = (domain, ("CG", 1, (domain.geometry.dim, )))
def clamped_boundary(x):
return np.isclose(x[0], 0)
def imposed_disp_boundary(x):
return np.isclose(x[0], L)
fdim = domain.topology.dim - 1
boundary_facets_I = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
boundary_facets_C = mesh.locate_entities_boundary(domain, fdim, imposed_disp_boundary)
u_C = np.array([0, 0], dtype=default_scalar_type)
u_I = np.array([0, 0], dtype=default_scalar_type)
bc_I = fem.dirichletbc(u_C, fem.locate_dofs_topological(V, fdim, boundary_facets_I), V)
bc_C = fem.dirichletbc(u_I, fem.locate_dofs_topological(V, fdim, boundary_facets_C), V)
E_young = 1000
nu = 0.3
lambda_lame = E_young * nu / ((1 + nu) * (1 - 2 *nu))
mu_lame = E_young /(2 * (1 + nu))
alpha = 1e-5
Temp = 100.
T = fem.Constant(domain, default_scalar_type((0, 0)))
ds = ufl.Measure("ds", domain=domain)
def epsilon(u):
return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(v, dT):
eps_th = alpha * dT
eps_el = epsilon(v) - eps_th * ufl.Identity(2)
return lambda_lame * ufl.tr(eps_el) * ufl.Identity(2) + 2.0 * mu_lame * eps_el
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0)))
a = ufl.inner(sigma(u, Temp), epsilon(v)) * ufl.dx
Lin = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
problem = LinearProblem(a, Lin, bcs=[bc_C, bc_I], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Uf = problem.solve()
I get the following error and I cannot figure what is going on… Can somebody help me?
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[29], line 27
24 a = ufl.inner(sigma(u, Temp), epsilon(v)) * ufl.dx
25 Lin = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
---> 27 problem = LinearProblem(a, Lin, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
28 Uf = problem.solve()
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/fem/petsc.py:782, in LinearProblem.__init__(self, a, L, bcs, u, petsc_options, form_compiler_options, jit_options)
744 def __init__(
745 self,
746 a: ufl.Form,
(...)
752 jit_options: typing.Optional[dict] = None,
753 ):
754 """Initialize solver for a linear variational problem.
755
756 Args:
(...)
780 "mumps"})
781 """
--> 782 self._a = _create_form(
783 a,
784 dtype=PETSc.ScalarType,
785 form_compiler_options=form_compiler_options,
786 jit_options=jit_options,
787 )
788 self._A = create_matrix(self._a)
789 self._L = _create_form(
790 L,
791 dtype=PETSc.ScalarType,
792 form_compiler_options=form_compiler_options,
793 jit_options=jit_options,
794 )
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/fem/forms.py:337, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
334 else:
335 return form
--> 337 return _create_form(form)
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/fem/forms.py:331, in form.<locals>._create_form(form)
329 return None
330 else:
--> 331 return _form(form)
332 elif isinstance(form, collections.abc.Iterable):
333 return list(map(lambda sub_form: _create_form(sub_form), form))
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/fem/forms.py:254, in form.<locals>._form(form)
252 if mesh is None:
253 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 254 ufcx_form, module, code = jit.ffcx_jit(
255 mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
256 )
258 # For each argument in form extract its function space
259 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/jit.py:62, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
58 @functools.wraps(local_jit)
59 def mpi_jit(comm, *args, **kwargs):
60 # Just call JIT compiler when running in serial
61 if comm.size == 1:
---> 62 return local_jit(*args, **kwargs)
64 # Default status (0 == ok, 1 == fail)
65 status = 0
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/dolfinx/jit.py:212, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
210 # Switch on type and compile, returning cffi object
211 if isinstance(ufl_object, ufl.Form):
--> 212 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
213 elif isinstance(ufl_object, ufl.Mesh):
214 r = ffcx.codegeneration.jit.compile_coordinate_maps([ufl_object], options=p_ffcx, **p_jit)
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:225, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
223 except Exception:
224 pass
--> 225 raise e
227 obj, module = _load_objects(cache_dir, module_name, form_names)
228 return obj, module, (decl, impl)
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:205, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
202 for name in form_names:
203 decl += form_template.format(name=name)
--> 205 impl = _compile_objects(
206 decl,
207 forms,
208 form_names,
209 module_name,
210 p,
211 cache_dir,
212 cffi_extra_compile_args,
213 cffi_verbose,
214 cffi_debug,
215 cffi_libraries,
216 visualise=visualise,
217 )
218 except Exception as e:
219 try:
220 # remove c file so that it will not timeout next time
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/codegeneration/jit.py:330, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
326 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
328 # JIT uses module_name as prefix, which is needed to make names of all struct/function
329 # unique across modules
--> 330 _, code_body = ffcx.compiler.compile_ufl_objects(
331 ufl_objects, prefix=module_name, options=options, visualise=visualise
332 )
334 # Raise error immediately prior to compilation if no support for C99
335 # _Complex. Doing this here allows FFCx to be used for complex codegen on
336 # Windows.
337 if sys.platform.startswith("win32"):
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/compiler.py:108, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
106 # Stage 1: analysis
107 cpu_time = time()
--> 108 analysis = analyze_ufl_objects(ufl_objects, options["scalar_type"]) # type: ignore
109 _print_timing(1, time() - cpu_time)
111 # Stage 2: intermediate representation
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/analysis.py:94, in analyze_ufl_objects(ufl_objects, scalar_type)
91 else:
92 raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
95 for data in form_data:
96 elements += data.unique_sub_elements
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/analysis.py:94, in <genexpr>(.0)
91 else:
92 raise TypeError("UFL objects not recognised.")
---> 94 form_data = tuple(_analyze_form(form, scalar_type) for form in forms)
95 for data in form_data:
96 elements += data.unique_sub_elements
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ffcx/analysis.py:180, in _analyze_form(form, scalar_type)
177 complex_mode = np.issubdtype(scalar_type, np.complexfloating)
179 # Compute form metadata
--> 180 form_data = ufl.algorithms.compute_form_data(
181 form,
182 do_apply_function_pullbacks=True,
183 do_apply_integral_scaling=True,
184 do_apply_geometry_lowering=True,
185 preserve_geometry_types=(ufl.classes.Jacobian,),
186 do_apply_restrictions=True,
187 do_append_everywhere_integrals=False, # do not add dx integrals to dx(i) in UFL
188 complex_mode=complex_mode,
189 )
191 # Determine unique quadrature degree and quadrature scheme
192 # per each integral data
193 for id, integral_data in enumerate(form_data.integral_data):
194 # Iterate through groups of integral data. There is one integral
195 # data for all integrals with same domain, itype, subdomain_id
(...)
199 # all integrals in this integral data group, i.e. must be the
200 # same for for the same (domain, itype, subdomain_id)
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/algorithms/compute_form_data.py:427, in compute_form_data(form, do_apply_function_pullbacks, do_apply_integral_scaling, do_apply_geometry_lowering, preserve_geometry_types, do_apply_default_restrictions, do_apply_restrictions, do_estimate_degrees, do_append_everywhere_integrals, complex_mode)
424 preprocessed_form = reconstruct_form_from_integral_data(self.integral_data)
426 # TODO: Test how fast this is
--> 427 check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)
429 # TODO: This member is used by unit tests, change the tests to
430 # remove this!
431 self.preprocessed_form = preprocessed_form
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/algorithms/check_arities.py:213, in check_form_arity(form, arguments, complex_mode)
211 """Check the arity of a form."""
212 for itg in form.integrals():
--> 213 check_integrand_arity(itg.integrand(), arguments, complex_mode)
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/algorithms/check_arities.py:194, in check_integrand_arity(expr, arguments, complex_mode)
192 arguments = tuple(sorted(set(arguments), key=lambda x: (x.number(), x.part())))
193 rules = ArityChecker(arguments)
--> 194 arg_tuples = map_expr_dag(rules, expr, compress=False)
195 args = tuple(a[0] for a in arg_tuples)
196 if args != arguments:
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/corealg/map_dag.py:35, in map_expr_dag(function, expression, compress, vcache, rcache)
15 def map_expr_dag(function, expression, compress=True, vcache=None, rcache=None):
16 """Apply a function to each subexpression node in an expression DAG.
17
18 If the same function is called multiple times in a transformation
(...)
33 The result of the final function call
34 """
---> 35 (result,) = map_expr_dags(
36 function, [expression], compress=compress, vcache=vcache, rcache=rcache
37 )
38 return result
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/corealg/map_dag.py:103, in map_expr_dags(function, expressions, compress, vcache, rcache)
101 r = handlers[v._ufl_typecode_](v)
102 else:
--> 103 r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
105 # Optionally check if r is in rcache, a memory optimization
106 # to be able to keep representation of result compact
107 if compress:
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/algorithms/check_arities.py:60, in ArityChecker.sum(self, o, a, b)
57 """Apply to sum."""
58 if a != b:
59 raise ArityMismatch(
---> 60 f"Adding expressions with non-matching form arguments {_afmt(a)} vs {_afmt(b)}."
61 )
62 return a
File ~/miniconda3/envs/fenicsx_env/lib/python3.13/site-packages/ufl/algorithms/check_arities.py:20, in _afmt(atuple)
18 def _afmt(atuple: Tuple[Argument, bool]) -> str:
19 """Return a string representation of an arity tuple."""
---> 20 arg, conj = atuple
21 return f"conj({arg})" if conj else str(arg)
ValueError: not enough values to unpack (expected 2, got 0)
