# Support needed to setup thermoelastic analysis

**URL:** <https://fenicsproject.discourse.group/t/support-needed-to-setup-thermoelastic-analysis/16169>\
**Category:** dolfinx\
**Created:** [October 31, 2024, 3:46pm UTC](https://fenicsproject.discourse.group/t/support-needed-to-setup-thermoelastic-analysis/16169 "2024-10-31T15:46:24Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jean-phi66](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jean-phi66/32/7744_2.png) [@jean-phi66](https://fenicsproject.discourse.group/u/jean-phi66)\
**Post date:** [October 31, 2024, 3:46pm UTC](https://fenicsproject.discourse.group/t/support-needed-to-setup-thermoelastic-analysis/16169/1 "2024-10-31T15:46:24Z")

</div>

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.

```auto
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?

```auto
---------------------------------------------------------------------------
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)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 31, 2024, 5:00pm UTC](https://fenicsproject.discourse.group/t/support-needed-to-setup-thermoelastic-analysis/16169/2 "2024-10-31T17:00:18Z")

</div>

> [@jean-phi66](#):
>
> ```auto
> 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
> 
> ```

since dT is a scalar,

> [@jean-phi66](#):
>
> `Temp = 100.`

> [@jean-phi66](#):
>
> `a = ufl.inner(sigma(u, Temp), epsilon(v)) * ufl.dx`

`a` is not a bilinear form.

I would rather write this on residual form,  
`F= ufl.inner(sigma(u, Temp), epsilon(v)) * ufl.dx- ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds`  
and then use `a, Lin = ufl.system(F)`

---

<div class="post-metadata">

**Author:** ![jean-phi66](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/jean-phi66/32/7744_2.png) [@jean-phi66](https://fenicsproject.discourse.group/u/jean-phi66)\
**Post date:** [November 1, 2024, 12:16pm UTC](https://fenicsproject.discourse.group/t/support-needed-to-setup-thermoelastic-analysis/16169/3 "2024-11-01T12:16:28Z")

</div>

It works! Thank you very much for your prompt and efficient answer 😃  
I am posting the corrected part if it can help somebody else in the future.

```auto
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

F= ufl.inner(sigma(u, Temp), epsilon(v)) * ufl.dx- ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds
a, Lin = ufl.system(F)

problem = LinearProblem(a, Lin, bcs=[bc_C, bc_I], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Uf = problem.solve()

```