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
512 stars 160 forks source link

Error assembling forms involving vector R0 functions #1456

Open jrmaddison opened 5 years ago

jrmaddison commented 5 years ago

For example:

from firedrake import *

mesh = UnitSquareMesh(10, 10)
vspace = VectorFunctionSpace(mesh, "R", 0)
test, trial = TestFunction(vspace), TrialFunction(vspace)

M = assemble(inner(test, trial) * dx)

leads to the error

UFL:ERROR Replacement expressions must have the same shape as what they replace.
Traceback (most recent call last):
  File "test.py", line 7, in <module>
    M = assemble(inner(test, trial) * dx)
  File "[..]/firedrake/src/firedrake/firedrake/assemble.py", line 105, in assemble
    loops = tuple(loops)
  File "[..]/firedrake/src/firedrake/firedrake/assemble.py", line 228, in _assemble
    kernels = tsfc_interface.compile_form(f, "form", parameters=form_compiler_parameters, inverse=inverse)
  File "[..]/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 206, in compile_form
    f = _real_mangle(f)
  File "[..]/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 233, in _real_mangle
    return ufl.replace(form, replacements)
  File "[..]/firedrake/src/ufl/ufl/algorithms/replace.py", line 69, in replace
    return map_integrand_dags(Replacer(mapping2), e)
  File "[..]/firedrake/src/ufl/ufl/algorithms/replace.py", line 38, in __init__
    error("Replacement expressions must have the same shape as what they replace.")
  File "[..]/firedrake/src/ufl/ufl/log.py", line 172, in error
    raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Replacement expressions must have the same shape as what they replace.
wence- commented 5 years ago

There are two issues here:

  1. _real_mangle assumes scalar RealFunctionSpace (one could fix that by just making a vector of 1 of appropriate shape)

  2. I don't think the rest of the code will work.

jrmaddison commented 4 years ago

Some related examples:

from firedrake import *

mesh = UnitIntervalMesh(10)

vspace = VectorFunctionSpace(mesh, "R", degree=0, dim=2)
v_F = Function(vspace, name="v_F")
v_F.assign(Constant((1.0, 1.0)))

print(assemble(dot(v_F, Constant((1.0, 1.0))) * dx))

and:

from firedrake import *

mesh = UnitIntervalMesh(10)

tspace = TensorFunctionSpace(mesh, "R", degree=0, shape=(2, 2))
t_F = Function(tspace, name="t_F")
t_F.assign(Constant(((1.0, 1.0), (1.0, 1.0))))

print(assemble(inner(t_F, Constant(((1.0, 1.0), (1.0, 1.0)))) * dx)) 

lead to LoopyError s (long error messages). Switching "R" to "DG" yields the correct results.