# Define different materials via the sub domain

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
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)

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

My understanding is that you are trying to define a circle/sphere around (x,y) = (0.1, 0.1). You may need to cut the visualization in paraview with a clip to see that.

I’ve already told you in the past to make sure to include all imports in your snippet, so that people who answer you can easily do a copy and paste. Last time I fixed that for you, this time I won’t, hence I won’t be trying to run your code (because it is not runnable). Please try to remember to do this: it’s a simple task, but you must put the others in the best position possible to help you.

added the imports, sorry for that

i tried the cute and as you can see still no change

It seems to me that you are trying to identify a circle in the cross section of the beam, but the mesh of each cross section is so coarse that you’ll never identify the circle.

The following code

``````from dolfin import *
import numpy as np

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, num_elements, num_elements)

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

subdomain_0 = Omega_0()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(2)
subdomain_0.mark(materials, 1)

# rest of the code
File("Materials.pvd") << materials
``````

does produce two subdomains. Note that:

• I had to increase the number of elements in the cross-section
• there is no need to have a `SubDomain` for \Omega_0, and then negate that expression to identify \Omega_1. Simply set the default value to 2 to say that all cells initially belong to \Omega_1, and then identify \Omega_0.

hi thanks for you reply and time checking it
but when run it this is the error i got

``````UMFPACK V5.7.9 (Oct 20, 2019): ERROR: out of memory

Traceback (most recent call last):
File "/home/mirialex/PycharmProjects/fibers/13-06-24-take2.py", line 390, in <module>
solver.solve()
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to successfully call PETSc function 'KSPSolve'.
*** Reason:  PETSc error code is: 76 (Error in external library).
*** Where:   This error was encountered inside ./dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

Process finished with exit code 1

``````

im using a solver once without time dependency and once with
without the time dependency im able to get resolution but once i getting to the loop of the time dependency it fails with the error above
i tried to lower the load steps in the loop or the number of the elements but nothing helped.

if you have any other piece of advice, would be more then happy to read

``````import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np

# rest of the code

du = TrialFunction(V)  # Incremental displacement
v = TestFunction(V)  # Test function
u = Function(V)  # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0))  # Body force per unit volume
T = Constant((0.0, 0.0, 0.0))  # Traction force on the boundary

# Project or interpolate the initial guess onto u
initial_guess = Expression(("x[0]", "x[1]", "x[2]"), degree=1)
# Project or interpolate the initial guess onto u
u.interpolate(initial_guess)
# u.interpolate(Constant((0,0,0)))

# Kinematics
d = len(u)
I = Identity(d)  # Identity tensor
F = variable(I + grad(u))  # Deformation gradient
C = F.T * F  # Right Cauchy-Green tensor

# Invariants of deformation tensors
Ic = tr(C)
J = det(F)

# Stored strain energy density (compressible neo-Hookean model)

W = MU / 2 * (tr(C) - 3 - 2 * ln(J)) + LAMBDA / 2 * pow((J - 1), 2)

n = FacetNormal(mesh)
P = diff(W, F)
# Total potential energy
L = inner(P, grad(v)) * dx - dot(B, v) * dx - dot(T, v) * ds

Je = derivative(L, u, du)

problem = NonlinearVariationalProblem(L, u, bcs=bcs,J=Je)

solver = NonlinearVariationalSolver(problem)
solver.parameters['newton_solver']['relative_tolerance'] = 1e-6
solver.parameters['newton_solver']['linear_solver'] = 'lu'

#rest of the code
dt = 0.0001
t, T = 0.0, 1 * dt
TT = Traction()

while t <= T:
print('time: ', t)

# Increase traction
TT.t = t
tDirBC.time_ = t
tDirBC4.time_ = t
#rest of the code
solver.solve()
t += float(dt)
``````

when i lower it back to number of elemnts of the number 4
its working but i still dont see the 2 materials

Try using `mumps` instead of `umfpack`. If you don’t know how to do that, it would be appreciated that you made a search (useful keywords could be `NonlinearVariationalSolver` and `mumps`) before asking us how to do that.

its working but i still dont see the 2 materials

As I told you, 4 elements in the cross section is too coarse to obtain a decent representation of the circle.
If you want to use a very coarse mesh, the suggestion would be to have it generated by gmsh (you know how to do that, since we discussed it last week), so that the subdomains are build during the mesh generation process, rather than afterwards.