FEniCS / dolfinx

Next generation FEniCS problem solving environment
https://fenicsproject.org
GNU Lesser General Public License v3.0
758 stars 182 forks source link

Codim 0 mixed assembly: Spatial coordinate without constant/coefficient #3098

Open jorgensd opened 7 months ago

jorgensd commented 7 months ago

Summarize the issue

Seems like ufl is not picking up meshes from spatial coordinate. First form has a constant from parent mesh, and it compiles. Second form only has spatial coordinate from parent and fails in signature computation

How to reproduce the bug

Run example below with main branch

Minimal Example (Python)

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 4, 2, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_entities(fdim)

# Get number of cells on process
cell_map = mesh.topology.index_map(tdim)
num_cells = cell_map.size_local + cell_map.num_ghosts

# Create markers for each size of the interface
cell_values = np.ones(num_cells, dtype=np.int32)
cell_values[dolfinx.mesh.locate_entities(
    mesh, tdim, lambda x: x[0] <= 0.5 + 1e-13)] = 2
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(
    num_cells, dtype=np.int32), cell_values)

submesh, sub_cell_to_parent, sub_vertex_to_parent, _ = dolfinx.mesh.create_submesh(mesh, ct.dim, ct.find(2))

mesh_to_submesh = np.full(num_cells, -1)
mesh_to_submesh[sub_cell_to_parent] = np.arange(len(sub_cell_to_parent), dtype=np.int32)

x = ufl.SpatialCoordinate(mesh)
d = dolfinx.fem.Constant(mesh, 1.0)
entity_maps = {mesh._cpp_object: np.array(sub_cell_to_parent, dtype=np.int32)}
dolfinx.fem.form(d * x[0] * ufl.dx(domain=submesh), entity_maps=entity_maps)
dolfinx.fem.form(x[0] * ufl.dx(domain=submesh), entity_maps=entity_maps)

Output (Python)

Traceback (most recent call last):
  File "/root/shared/mwe_2.py", line 34, in <module>
    dolfinx.fem.form(x[0] * ufl.dx(domain=submesh), entity_maps=entity_maps)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 249, in form
    return _create_form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 244, in _create_form
    return _form(form)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/forms.py", line 186, in _form
    ufcx_form, module, code = jit.ffcx_jit(
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 51, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 201, in ffcx_jit
    r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
  File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 225, in compile_forms
    module_name = "libffcx_forms_" + ffcx.naming.compute_signature(
  File "/usr/local/lib/python3.10/dist-packages/ffcx/naming.py", line 39, in compute_signature
    object_signature += ufl_object.signature()
  File "/usr/local/lib/python3.10/dist-packages/ufl/form.py", line 445, in signature
    self._compute_signature()
  File "/usr/local/lib/python3.10/dist-packages/ufl/form.py", line 698, in _compute_signature
    self._signature = compute_form_signature(self, self._compute_renumbering())
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/signature.py", line 137, in compute_form_signature
    terminal_hashdata = compute_terminal_hashdata(integrands, renumbering)
  File "/usr/local/lib/python3.10/dist-packages/ufl/algorithms/signature.py", line 74, in compute_terminal_hashdata
    data = expr._ufl_signature_data_(renumbering)
  File "/usr/local/lib/python3.10/dist-packages/ufl/geometry.py", line 105, in _ufl_signature_data_
    return (self._ufl_class_.__name__,) + self._domain._ufl_signature_data_(renumbering)
  File "/usr/local/lib/python3.10/dist-packages/ufl/domain.py", line 120, in _ufl_signature_data_
    return ("Mesh", renumbering[self], self._ufl_coordinate_element)
KeyError: Mesh(blocked element (Basix element (P, triangle, 1, gll_warped, unset, False, float64, []), (2,)), 0)

Version

main branch

DOLFINx git commit

1ccffb0aa545634ea6e160236a7b043dc2a7fdd2

Installation

Using docker with ghcr.io/fenics/dolfinx/dolfinx:nightly

Additional information

No response

jorgensd commented 7 months ago

The problem is likely in ufl: form.py::_analyze_domains(), as it doesn't extract domains from constants, arguments, coefficients. (https://github.com/FEniCS/ufl/blob/main/ufl/form.py#L604-L616)

jorgensd commented 7 months ago

Or, a problem of the re-numbering algorithm: https://github.com/FEniCS/ufl/blob/main/ufl/form.py#L658-L692 as it does not set numbering for geometric quantities used in https://github.com/FEniCS/ufl/blob/main/ufl/algorithms/signature.py#L73-L74