interpolate of symmetric tensor is broken

mmorandi commented 4 months ago

Summarize the issue

Consider the attached test case. Interpolate is asking for 36 numbers even if the tensor is symmetric. If one passes a full tensor (mat_modulus_symm1.interpolate) then it silently truncates the tensor; if, on the other hand, one passes the upper triangular portion of the tensor components (mat_modulus_symm2.interpolate ), then it bails out

How to reproduce the bug

Check the reproducer below

Minimal Example (Python)

from mpi4py import MPI
import numpy as np

import basix
import dolfinx

def tensor(x, tens):
    tens_reshaped = tens.reshape((6 * 6, 1))
    return np.broadcast_to(tens_reshaped, (6 * 6, x.shape[1]))

def symm_tensor(x, tens):
    tens_reshaped = tens[np.triu_indices(6)].reshape((21, 1))
    return np.broadcast_to(tens_reshaped, (21, x.shape[1]))

mesh = dolfinx.mesh.create_interval(nx=1, points=((-1),(1)),comm=MPI.COMM_WORLD)
MODULUS_ELEMENT = basix.ufl.element("DG", mesh.basix_cell(), 0, shape=(6,6))
MODULUS_FS = dolfinx.fem.FunctionSpace(mesh, MODULUS_ELEMENT)
mat_modulus = dolfinx.fem.Function(MODULUS_FS, name = "ElasticModulus")

MODULUS_ELEMENT_SYMM = basix.ufl.element("DG", mesh.basix_cell(), 0, shape=(6,6), symmetry=True)
MODULUS_FS_SYMM = dolfinx.fem.FunctionSpace(mesh, MODULUS_ELEMENT_SYMM)
mat_modulus_symm1 = dolfinx.fem.Function(MODULUS_FS_SYMM, name = "ElasticModulusSymm")
mat_modulus_symm2 = dolfinx.fem.Function(MODULUS_FS_SYMM, name = "ElasticModulusSymm")

mat = np.array([
        [11, 12, 13, 14, 15, 16],
        [12, 22, 23, 24, 25, 26],
        [13, 23, 33, 34, 35, 36],
        [14, 24, 34, 44, 45, 46],
        [15, 25, 35, 45, 55, 56],
        [16, 26, 36, 46, 56, 66]])

def Domain_0(x):
    return np.ones(x.shape[1])

domain = dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, Domain_0)

#this works
mat_modulus.interpolate(lambda x: tensor(x, mat), domain)
print(mat_modulus.x.array)
print("==================")

#this works, but is wrong
mat_modulus_symm1.interpolate(lambda x: tensor(x, mat), domain)
print(mat_modulus_symm1.x.array)

#this should work, should be correct, but does not work
mat_modulus_symm2.interpolate(lambda x: symm_tensor(x, mat), domain)

Output (Python)

RuntimeError: Interpolation data has the wrong shape/size.
(fenicsx-env) marco@dark:~/Projects/Chierichetti/AnbaDolfin/ReleaseGPL (master *=)> python3 anbax/TestInterpolate.py 
[11. 12. 13. 14. 15. 16. 12. 22. 23. 24. 25. 26. 13. 23. 33. 34. 35. 36.
 14. 24. 34. 44. 45. 46. 15. 25. 35. 45. 55. 56. 16. 26. 36. 46. 56. 66.]
==================
[11. 12. 13. 14. 15. 16. 12. 22. 23. 24. 25. 26. 13. 23. 33. 34. 35. 36.
 14. 24. 34.]
Traceback (most recent call last):
  File "/home/marco/local/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 397, in interpolate
    _interpolate(u, cells)
  File "/home/marco/local/miniconda3/envs/fenicsx-env/lib/python3.12/functools.py", line 909, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/marco/local/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 373, in _interpolate
    self._cpp_object.interpolate(u, cells, nmm_interpolation_data)
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.fem.Function_float64, f: numpy.ndarray[numpy.float64], cells: numpy.ndarray[numpy.int32]) -> None
    2. (self: dolfinx.cpp.fem.Function_float64, u: dolfinx.cpp.fem.Function_float64, cells: numpy.ndarray[numpy.int32], nmm_interpolation_data: Tuple[List[int], List[int], List[float], List[int]]) -> None
    3. (self: dolfinx.cpp.fem.Function_float64, expr: dolfinx::fem::Expression<double, double>, cells: numpy.ndarray[numpy.int32]) -> None

Invoked with: <dolfinx.cpp.fem.Function_float64 object at 0x7ff9301f47f0>, <function <lambda> at 0x7ff9301d6e80>, array([0], dtype=int32), ((), (), (), ())

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/marco/Projects/Chierichetti/AnbaDolfin/ReleaseGPL/anbax/TestInterpolate.py", line 49, in <module>
    mat_modulus_symm2.interpolate(lambda x: symm_tensor(x, mat), domain)
  File "/home/marco/local/miniconda3/envs/fenicsx-env/lib/python3.12/site-packages/dolfinx/fem/function.py", line 402, in interpolate
    self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
RuntimeError: Interpolation data has the wrong shape/size.

Version

0.7.3

DOLFINx git commit

b488a495d3d74a048e5866f437f2df20b075d907

Installation

conda on a linux box

Additional information

No response

jorgensd commented 4 months ago

@mscroggs The issue here is that we use u.function_space().value_size() within interpolate, which is (6, 6), while the basix element reference-value-shape is (21, ).

Similarly, the value size is 36.

What would be the best thing to do here?

mscroggs commented 4 months ago

Following #3158, the first option for interpolation into a symmetric element will now be correct rather than doing something incorrect, ie:

mat_modulus_symm1.interpolate(lambda x: tensor(x, mat), domain)
print(mat_modulus_symm1.x.array)

This fix will be included in the 0.8 release.

I'd like to get the second option working too. I've opened #3159 for this, with a unit test for this included

bleyerj commented 4 months ago

Just to be clear, what happens if with option 1 you pass a full array which is non-symmetric ?

mscroggs commented 4 months ago

Just to be clear, what happens if with option 1 you pass a full array which is non-symmetric ?

It just looks at the lower part of the matrix and ignores the upper part. Maybe we need an issue to change this at some point

I've opened an issue to discuss how to handle this case: https://github.com/FEniCS/dolfinx/issues/3165

(edit: swapped upper and lower)

FEniCS / dolfinx