Hello everyone,
I have been trying to couple PDEs of multiple dimensions(3D acoustic domain and 2D elastic plate domain), and fortunately, I have found a useful demo for doing this.
The problem here is that I got the same error as I was trying to fix with my code, which seems to be that the dimension of the integration measure does not match with the dimension of the integrands.
Now, I want to run this demo without errors since it’s such a good reference.
I believe the demo usually runs without errors, so this error could be due to that I am using FEniCSx 0.8.0 for a compatibility reason-- in that case, I would greatly appreciate if someone could give me a hint to modify the code so that this demo runs on 0.8.0.
Below is the block of the code in the above demo, followed by the error produced.
Thank you so much for your time and I am looking forward to hearing your advice.
v, w = ufl.TestFunctions(W)
u = dolfinx.fem.Function(V, name="Displacement")
psi = dolfinx.fem.Function(Q, name="Latent_variable")
psi_k = dolfinx.fem.Function(Q, name="Previous_Latent_variable")
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
def epsilon(w):
return ufl.sym(ufl.grad(w))
def sigma(w, mu, lmbda):
ew = epsilon(w)
gdim = ew.ufl_shape[0]
return 2.0 * mu * epsilon(w) + lmbda * ufl.tr(ufl.grad(w)) * ufl.Identity(gdim)
F = alpha * ufl.inner(sigma(u, mu, lmbda), epsilon(v)) * dx
F -= alpha * ufl.inner(f, v) * dx
F += -ufl.inner(psi - psi_k, ufl.dot(v, n)) * ds
F += ufl.inner(ufl.dot(u, n_g), w) * ds
F += ufl.inner(ufl.exp(psi), w) * ds - ufl.inner(g, w) * ds
residual = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
Cell In[11], line 23
21 F += ufl.inner(ufl.dot(u, n_g), w) * ds
22 F += ufl.inner(ufl.exp(psi), w) * ds - ufl.inner(g, w) * ds
---> 23 residual = dolfinx.fem.form(ufl.extract_blocks(F), entity_maps=entity_maps)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/fem/forms.py:249, in form(form, dtype, form_compiler_options, jit_options, entity_maps)
246 return list(map(lambda sub_form: _create_form(sub_form), form))
247 return form
--> 249 return _create_form(form)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/fem/forms.py:246, in form.<locals>._create_form(form)
244 return _form(form)
245 elif isinstance(form, collections.abc.Iterable):
--> 246 return list(map(lambda sub_form: _create_form(sub_form), form))
247 return form
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/fem/forms.py:246, in form.<locals>._create_form.<locals>.<lambda>(sub_form)
244 return _form(form)
245 elif isinstance(form, collections.abc.Iterable):
--> 246 return list(map(lambda sub_form: _create_form(sub_form), form))
247 return form
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/fem/forms.py:244, in form.<locals>._create_form(form)
241 """Recursively convert ufl.Forms to dolfinx.fem.Form, otherwise
242 return form argument"""
243 if isinstance(form, ufl.Form):
--> 244 return _form(form)
245 elif isinstance(form, collections.abc.Iterable):
246 return list(map(lambda sub_form: _create_form(sub_form), form))
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/fem/forms.py:186, in form.<locals>._form(form)
184 if mesh is None:
185 raise RuntimeError("Expecting to find a Mesh in the form.")
--> 186 ufcx_form, module, code = jit.ffcx_jit(
187 mesh.comm, form, form_compiler_options=form_compiler_options, jit_options=jit_options
188 )
190 # For each argument in form extract its function space
191 V = [arg.ufl_function_space()._cpp_object for arg in form.arguments()]
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/jit.py:51, in mpi_jit_decorator.<locals>.mpi_jit(comm, *args, **kwargs)
47 @functools.wraps(local_jit)
48 def mpi_jit(comm, *args, **kwargs):
49 # Just call JIT compiler when running in serial
50 if comm.size == 1:
---> 51 return local_jit(*args, **kwargs)
53 # Default status (0 == ok, 1 == fail)
54 status = 0
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/dolfinx/jit.py:201, in ffcx_jit(ufl_object, form_compiler_options, jit_options)
199 # Switch on type and compile, returning cffi object
200 if isinstance(ufl_object, ufl.Form):
--> 201 r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
202 elif isinstance(ufl_object, ufl.AbstractFiniteElement):
203 r = ffcx.codegeneration.jit.compile_elements([ufl_object], options=p_ffcx, **p_jit)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:276, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
274 except Exception:
275 pass
--> 276 raise e
278 obj, module = _load_objects(cache_dir, module_name, form_names)
279 return obj, module, (decl, impl)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:256, in compile_forms(forms, options, cache_dir, timeout, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
253 for name in form_names:
254 decl += form_template.format(name=name)
--> 256 impl = _compile_objects(
257 decl,
258 forms,
259 form_names,
260 module_name,
261 p,
262 cache_dir,
263 cffi_extra_compile_args,
264 cffi_verbose,
265 cffi_debug,
266 cffi_libraries,
267 visualise=visualise,
268 )
269 except Exception as e:
270 try:
271 # remove c file so that it will not timeout next time
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/codegeneration/jit.py:383, in _compile_objects(decl, ufl_objects, object_names, module_name, options, cache_dir, cffi_extra_compile_args, cffi_verbose, cffi_debug, cffi_libraries, visualise)
379 libraries = _libraries + cffi_libraries if cffi_libraries is not None else _libraries
381 # JIT uses module_name as prefix, which is needed to make names of all struct/function
382 # unique across modules
--> 383 _, code_body = ffcx.compiler.compile_ufl_objects(
384 ufl_objects, prefix=module_name, options=options, visualise=visualise
385 )
387 ffibuilder = cffi.FFI()
389 ffibuilder.set_source(
390 module_name,
391 code_body,
(...)
394 libraries=libraries,
395 )
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/compiler.py:113, in compile_ufl_objects(ufl_objects, options, object_names, prefix, visualise)
111 # Stage 2: intermediate representation
112 cpu_time = time()
--> 113 ir = compute_ir(analysis, _object_names, _prefix, options, visualise)
114 _print_timing(2, time() - cpu_time)
116 # Stage 3: code generation
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/ir/representation.py:225, in compute_ir(analysis, object_names, prefix, options, visualise)
214 ir_elements = [
215 _compute_element_ir(e, analysis.element_numbers, finite_element_names)
216 for e in analysis.unique_elements
217 ]
219 ir_dofmaps = [
220 _compute_dofmap_ir(e, analysis.element_numbers, dofmap_names)
221 for e in analysis.unique_elements
222 ]
224 irs = [
--> 225 _compute_integral_ir(
226 fd,
227 i,
228 analysis.element_numbers,
229 integral_names,
230 finite_element_names,
231 options,
232 visualise,
233 )
234 for (i, fd) in enumerate(analysis.form_data)
235 ]
236 ir_integrals = list(itertools.chain(*irs))
238 ir_forms = [
239 _compute_form_ir(
240 fd,
(...)
250 for (i, fd) in enumerate(analysis.form_data)
251 ]
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/ir/representation.py:588, in _compute_integral_ir(form_data, form_index, element_numbers, integral_names, finite_element_names, options, visualise)
583 integrand_map: dict[QuadratureRule, ufl.core.expr.Expr] = {
584 rule: integral.integrand() for rule, integral in sorted_integrals.items()
585 }
587 # Build more specific intermediate representation
--> 588 integral_ir = compute_integral_ir(
589 itg_data.domain.ufl_cell(),
590 itg_data.integral_type,
591 ir["entitytype"],
592 integrand_map,
593 ir["tensor_shape"],
594 options,
595 visualise,
596 )
598 ir.update(integral_ir)
600 # Fetch name
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/ir/integral.py:85, in compute_integral_ir(cell, integral_type, entitytype, integrands, argument_shape, p, visualise)
73 # Build terminal_data from V here before factorization. Then we
74 # can use it to derive table properties for all modified
75 # terminals, and then use that to rebuild the scalar graph more
76 # efficiently before argument factorization. We can build
77 # terminal_data again after factorization if that's necessary.
79 initial_terminals = {
80 i: analyse_modified_terminal(v["expression"])
81 for i, v in S.nodes.items()
82 if is_modified_terminal(v["expression"])
83 }
---> 85 mt_table_reference = build_optimized_tables(
86 quadrature_rule,
87 cell,
88 integral_type,
89 entitytype,
90 initial_terminals.values(),
91 ir["unique_tables"],
92 use_sum_factorization=p["sum_factorization"],
93 rtol=p["table_rtol"],
94 atol=p["table_atol"],
95 )
97 # Fetch unique tables for this quadrature rule
98 table_types = {v.name: v.ttype for v in mt_table_reference.values()}
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/ir/elementtables.py:415, in build_optimized_tables(quadrature_rule, cell, integral_type, entitytype, modified_terminals, existing_tables, use_sum_factorization, rtol, atol)
413 t["array"] = np.vstack([td["array"] for td in new_table])
414 else:
--> 415 t = get_ffcx_table_values(
416 quadrature_rule.points,
417 cell,
418 integral_type,
419 element,
420 avg,
421 entitytype,
422 local_derivatives,
423 flat_component,
424 )
425 # Clean up table
426 tbl = clamp_table_small_numbers(t["array"], rtol=rtol, atol=atol)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/ffcx/ir/elementtables.py:139, in get_ffcx_table_values(points, cell, integral_type, element, avg, entitytype, derivative_counts, flat_component)
137 for entity in range(num_entities):
138 entity_points = map_integral_points(points, integral_type, cell, entity)
--> 139 tbl = component_element.tabulate(deriv_order, entity_points)
140 tbl = tbl[basix_index(derivative_counts)]
141 component_tables.append(tbl)
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/basix/ufl.py:660, in _ComponentElement.tabulate(self, nderivs, points)
650 def tabulate(self, nderivs: int, points: _npt.NDArray[np.float64]) -> _npt.NDArray[np.float64]:
651 """Tabulate the basis functions of the element.
652
653 Args:
(...)
658 Tabulated basis functions.
659 """
--> 660 tables = self._element.tabulate(nderivs, points)
661 output = []
662 for tbl in tables:
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/basix/ufl.py:415, in _BasixElement.tabulate(self, nderivs, points)
404 def tabulate(self, nderivs: int, points: _npt.NDArray[np.float64]) -> _npt.NDArray[np.float64]:
405 """Tabulate the basis functions of the element.
406
407 Args:
(...)
413
414 """
--> 415 tab = self._element.tabulate(nderivs, points)
416 # TODO: update FFCx to remove the need for transposing here
417 return tab.transpose((0, 1, 3, 2)).reshape((tab.shape[0], tab.shape[1], -1))
File /opt/homebrew/Caskroom/miniconda/base/envs/try10/lib/python3.12/site-packages/basix/finite_element.py:137, in FiniteElement.tabulate(self, n, x)
105 def tabulate(self, n: int, x: npt.NDArray) -> npt.NDArray[np.floating]:
106 """Compute basis values and derivatives at set of points.
107
108 Note:
(...)
135 basis functions.
136 """
--> 137 return self._e.tabulate(n, x)
RuntimeError: Point dim (2) does not match element dim (1).