hi all,

im trying to set different material properties on my model and for some reason in the view of paraview it does not seems to be working, i dont see my model colored in 2 different colors

when i tried different subdomain its working as expected

for exp

```
import os
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin
#rest of code
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 0.5 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 0.5 - DOLFIN_EPS
```

but when i try this subdomain config its not working

```
E1, E2, nu = 30e6, 3e6, 0.49
# Create mesh and define function space
length, width, height = 1.0, 0.2, 0.2 # Dimensions of the rod
num_elements = 40 # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 4, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
#creating subtomains and boundaty conditions
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return np.sqrt(((x[1] - 0.1)*(x[1] - 0.1)) + ((x[2] - 0.1)*(x[2] - 0.1))) <= 0.05
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return np.sqrt(((x[1] - 0.1)*(x[1] - 0.1)) + ((x[2] - 0.1)*(x[2] - 0.1))) > 0.05
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 1))
up = AutoSubDomain(lambda x: near(x[2], 0.2))
down = AutoSubDomain(lambda x: near(x[2], 0))
front = AutoSubDomain(lambda x: near(x[1], 0))
back = AutoSubDomain(lambda x: near(x[1], 0.2))
left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
up.mark(boundary_parts, 5)
down.mark(boundary_parts, 6)
front.mark(boundary_parts, 3)
back.mark(boundary_parts, 4)
bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
bc3 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 3)
bc5 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 5)
bc6 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 6)
tDirBC = Expression(('1.0*time_'), time_=0.0, degree=0)
tDirBC4 = Expression(('((1/(time_+1))-1)*0.2'), time_=0.0, degree=0)
bc2 = DirichletBC(V.sub(0), tDirBC, boundary_parts, 2)
bc4 = DirichletBC(V.sub(1), tDirBC4, boundary_parts, 4)
bcs = [bc1, bc2, bc3, bc4, bc5, bc6]
subdomain_ids = [2, 4]
ds = Measure("ds", subdomain_data=boundary_parts, subdomain_id=(2,4))
ds = ds(degree=4)
#creating material constants
class K(fe.UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
E_global = K(materials, E1, E2, degree=0)
DG = FunctionSpace(mesh, "DG", 2)
materials_function = Function(DG)
E_values = [E1, E2]
materials_function = project(E_global, DG)
MU1, LAMBDA1 = Constant(E1 / (2 * (1 + nu))), Constant(E1 * nu / ((1 + nu) * (1 - 2 * nu)))
MU2, LAMBDA2 = Constant(E2 / (2 * (1 + nu))), Constant(E2 * nu / ((1 + nu) * (1 - 2 * nu)))
class MUMU(fe.UserExpression):
def __init__(self, materials, MU_0, MU_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.MU_0 = MU_0
self.MU_1 = MU_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.MU_0
else:
values[0] = self.MU_1
def value_shape(self):
return ()
MU = MUMU(materials, MU1, MU2, degree=0)
class lamabada(fe.UserExpression):
def __init__(self, materials, lam_0, lam_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.lam_0 = lam_0
self.lam_1 = lam_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.lam_0
else:
values[0] = self.lam_1
def value_shape(self):
return ()
LAMBDA = lamabada(materials, LAMBDA1, LAMBDA2, degree=0)
# rest of the code
File("Materials.pvd") << materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials
#writing into files in order to check it in paraview
materials_function.rename("materials", "")
u.rename("Displacement Vector", "")
beam_deflection_file = fe.XDMFFile("beam_deflection.xdmf")
beam_deflection_file.parameters["flush_output"] = True
beam_deflection_file.parameters["functions_share_mesh"] = True
beam_deflection_file.write(u, 0.0)
beam_deflection_file.write(materials_function, 0.0)
```

would appreciate to know if you have any idea why

this is what i see