Different material properties in subdomains, tensor instead of scalar

General
2

How do you expect a tensor to be assigned to a scalar function space?
The space ("CG", 1) contains one degree of freedom at any point. It does not make sense to try to add a scalar to it.
If you want to do so, you need to use a TensorElement or a TensorFunctionSpace as shown for a dummy rho below:

import dolfinx
from mpi4py import MPI
import ufl
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

el = ufl.TensorElement("Lagrange", mesh.ufl_cell(), 1, shape=(2, 2))

V = dolfinx.fem.FunctionSpace(mesh, el)
u = dolfinx.fem.Function(V)


def rho_B(x):
    tensor = np.array([[1, 2], [3, 4.]])

    values = np.repeat(tensor, x.shape[1])
    return values.reshape(tensor.shape[0]*tensor.shape[1], x.shape[1])


u.interpolate(rho_B)

for i in range(2):
    for j in range(2):
        print(i, j, dolfinx.fem.assemble_scalar(
            dolfinx.fem.form(u[i, j]*ufl.dx)))

if you want the rho to be different in different subdomains, you should use the cells kwarg of interpolate, as it will then only do interpolation on dofs associated to those entities. In your case that would be
u.interpolate(rho_B, cells=ct.find(11))

1 Like