Define multiple sets of material properties

General
1

Hello, everyone, I wonder if I was correct in defining the material properties (stiffness matrices) of the five physics groups in this way, as well as the weak integral form.

# read the mesh file
mesh_cell, cell_tag, facet_tag = dolfinx.io.gmshio.read_from_msh(msh_path, mpi4py.MPI.COMM_WORLD, 0, gdim=3)

# extract the corresponding cells by physical group tags
tags = cell_tag.values
layer_0 = np.where(tags == 1)[0]
layer_1 = np.where(tags == 2)[0]
layer_2 = np.where(tags == 3)[0]
layer_3 = np.where(tags == 4)[0]
layer_4 = np.where(tags == 5)[0]

# Create a DG0 function space
Q = fem.functionspace(mesh_cell, ("DG", 0))
indicator = fem.Function(Q)
indicator.x.array[layer_0] = 0
indicator.x.array[layer_1] = 1
indicator.x.array[layer_2] = 2
indicator.x.array[layer_3] = 3
indicator.x.array[layer_4] = 4
indicator.x.scatter_forward()

C_layer_0 = fem.Constant(mesh_cell,  C_layers[1].astype(PETSc.ScalarType))
C_layer_1 = fem.Constant(mesh_cell,  C_layers[2].astype(PETSc.ScalarType))
C_layer_2 = fem.Constant(mesh_cell,  C_layers[3].astype(PETSc.ScalarType))
C_layer_3 = fem.Constant(mesh_cell,  C_layers[4].astype(PETSc.ScalarType))
C_layer_4 = fem.Constant(mesh_cell,  C_layers[5].astype(PETSc.ScalarType))


def get_C(indicator):
    C = ufl.conditional(ufl.eq(indicator, 0), C_layer_0,
                        ufl.conditional(ufl.eq(indicator, 1), C_layer_1,
                                        ufl.conditional(ufl.eq(indicator, 2), C_layer_2,
                                                        ufl.conditional(ufl.eq(indicator, 3), C_layer_3,
                                                                        C_layer_4))))
    return C


# 3D strain-stress relationship
def strain_3D(u, repr="vectorial"):
    eps = ufl.sym(ufl.grad(u))
    if repr == "vectorial":
        return ufl.as_vector([
            eps[0, 0], eps[1, 1], eps[2, 2],
            2*eps[1, 2], 2*eps[0, 2], 2*eps[0, 1]
        ])
    elif repr == "tensorial":
        return eps


def stress_3D(u, indicator, repr="vectorial"):
    eps_vec = strain_3D(u, "vectorial")
    C = get_C(indicator)
    sigma_vec = ufl.dot(C, eps_vec)
    if repr == "vectorial":
        return sigma_vec
    elif repr == "tensorial":
        return ufl.as_matrix([
            [sigma_vec[0], sigma_vec[5]/2, sigma_vec[4]/2],
            [sigma_vec[5]/2, sigma_vec[1], sigma_vec[3]/2],
            [sigma_vec[4]/2, sigma_vec[3]/2, sigma_vec[2]]
        ])


V = fem.functionspace(mesh_cell, ("Lagrange", 1, (mesh_cell.geometry.dim,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

dx = ufl.Measure("dx", domain=mesh_cell, metadata={"quadrature_degree": 2})
a = ufl.inner(stress_3D(u, indicator), strain_3D(v)) * dx

I have made an extension based on this example.> Blockquote
Thanks!

2

Hi. This should work, you can also adapt this function to matrix-valued function. This would avoid defining conditionals which is generally feasible only for simple geometric subdomains

3

Thanks a lot! This has helped me a lot.