search

firedrakeproject / firedrake

Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM)
https://firedrakeproject.org
Other
498 stars 157 forks source link

BUG: unable to save CG 2 field to `pvd` when using hexahedral mesh #3315

Closed francesco-ballarin closed 7 months ago

francesco-ballarin commented 8 months ago

Hi, I am running the following code in the latest firedrakeproject/firedrake docker image:

from firedrake import *

# mesh = firedrake.UnitSquareMesh(2, 2, quadrilateral=True)  # quad mesh works
mesh = firedrake.UnitCubeMesh(2, 2, 2, hexahedral=True)  # hex mesh does not
V = FunctionSpace(mesh, "CG", 2)
f = Function(V)
f.interpolate(sin(SpatialCoordinate(mesh)[0]))
outfile = File("output.pvd")
outfile.write(f)

and getting

Traceback (most recent call last):
  File "/home/firedrake/tmp.py", line 9, in <module>
    outfile.write(f)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 650, in write
    vtu = self._write_vtu(*functions)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 537, in _write_vtu
    self._topology = get_topology(coordinates.function)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 180, in get_topology
    perm = vtk_lagrange_hex_reorder(V.ufl_element())
  File "/home/firedrake/firedrake/src/firedrake/firedrake/paraview_reordering.py", line 222, in vtk_lagrange_hex_reorder
    degree = max(ufl_element.degree())
TypeError: 'int' object is not iterable

Considering that the quad case works, I tried replacing https://github.com/firedrakeproject/firedrake/blob/c5e939ddeda6916433cafa251fb341c2e12aa312/firedrake/paraview_reordering.py#L222-L224 with https://github.com/firedrakeproject/firedrake/blob/c5e939ddeda6916433cafa251fb341c2e12aa312/firedrake/paraview_reordering.py#L201 which seems to me a different way of carry out the same task.

After this small change, the new error is

Traceback (most recent call last):
  File "/home/firedrake/tmp.py", line 9, in <module>
    outfile.write(f)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 650, in write
    vtu = self._write_vtu(*functions)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 537, in _write_vtu
    self._topology = get_topology(coordinates.function)
  File "petsc4py/PETSc/Log.pyx", line 188, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "petsc4py/PETSc/Log.pyx", line 189, in petsc4py.PETSc.Log.EventDecorator.decorator.wrapped_func
  File "/home/firedrake/firedrake/src/firedrake/firedrake/output.py", line 180, in get_topology
    perm = vtk_lagrange_hex_reorder(V.ufl_element())
  File "/home/firedrake/firedrake/src/firedrake/firedrake/paraview_reordering.py", line 223, in vtk_lagrange_hex_reorder
    firedrake_local = firedrake_local_to_cart(ufl_element)
  File "/home/firedrake/firedrake/src/firedrake/firedrake/paraview_reordering.py", line 19, in firedrake_local_to_cart
    fiat_element = create_base_element(element).fiat_equivalent
  File "/home/firedrake/firedrake/src/tsfc/gem/utils.py", line 19, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/home/firedrake/firedrake/src/FInAT/finat/cube.py", line 61, in fiat_equivalent
    return FIAT_FlattenedDimensions(self.product.fiat_equivalent)
  File "/home/firedrake/firedrake/src/tsfc/gem/utils.py", line 19, in __get__
    obj.__dict__[self.__name__] = result = self.fget(obj)
  File "/home/firedrake/firedrake/src/FInAT/finat/tensor_product.py", line 77, in fiat_equivalent
    A, B = self.factors
ValueError: too many values to unpack (expected 2)

If I add a print at line 77 of /home/firedrake/firedrake/src/FInAT/finat/tensor_product.py, indeed self.factors is composed of 3 elements. Yet, there is a comment # FIAT TensorProductElement support only 2 factors just before that line, and I don't feel adventorous enough to go and change FIAT and FInAT myself :)

francesco-ballarin commented 8 months ago

Furthermore, https://github.com/firedrakeproject/firedrake/blob/c5e939ddeda6916433cafa251fb341c2e12aa312/firedrake/paraview_reordering.py#L10 seems suspicious: any recent version of paraview uses vtk 9. For instance, I have paraview 5.11.2 packaged by debian and it uses vtk 9.2.