How to set diferent materials on difrent parts of a mesh?

Question

So say we have a mesh: UnitCubeMesh(24, 16, 16) and we set its material properties to

# Elasticity parameters
E, nu = 10.0, 0.3

(as in example here). How could we set its top left half to be with diferent material properties say

# Elasticity parameters
E2, nu2 = 4.5, 0.9

and top right say

# Elasticity parameters
E3, nu3 = 15, 0.4

So it would be one mesh composed of 3 materials like so:

?

Could not find such sample in docs/undocumented samples=(

Answer 1

Hi, as an alternative to approach given by @V_L where the material properties are cell wise constant functions you could also consider the following

# Define domain of individual materials
eps = DOLFIN_EPS
domain0 = AutoSubDomain(lambda x: x[2] < 0.5 + eps)
domain1 = AutoSubDomain(lambda x: x[2] > 0.5 - eps and x[0] < 0.5 + eps)
domain2 = AutoSubDomain(lambda x: x[2] > 0.5 - eps and x[0] > 0.5 - eps)

# Have one function with tags of domains
domains = CellFunction('size_t', mesh, 0)
domain0.mark(domains, 0)
domain1.mark(domains, 1)
domain2.mark(domains, 2)

# Material parameters E, nu corresponding to domain 0, 1, 2
materials = [(10, 0.3), (4.5, 0.3), (15, 0.4)]

# Get mu lambda from E, nu
materials = map(lambda (E, nu): (Constant(E/(2*(1+nu))),
                                 Constant(E*nu/((1+nu)*(1-2*nu)))),
                materials)

lmbda0, mu0 = materials[0]
lmbda1, mu1 = materials[1]
lmbda2, mu2 = materials[2]
# Stored strain energy density (compressible neo-Hookean model)
psi0 = (mu0/2)*(Ic - 3) - mu0*ln(J) + (lmbda0/2)*(ln(J))**2
psi1 = (mu1/2)*(Ic - 3) - mu1*ln(J) + (lmbda1/2)*(ln(J))**2
psi2 = (mu2/2)*(Ic - 3) - mu2*ln(J) + (lmbda2/2)*(ln(J))**2

# Total potential energy, each domain with its own strain energy density
dx = Measure('dx')[domains]
Pi = psi0*dx(0) + psi1*dx(1) + psi2*dx(2)
Pi += - dot(B, u)*dx('everywhere') - dot(T, u)*ds  

Btw, should the value of $\nu_2$ really be 0.9? Isn't Poisson ratio always supposed to be less then $1/2$? With your given value the Newton solver actually does not converge.

Comments

Hello Mirok,

Would this dx() approach work in C++ as well? I could not find any demo where this has been attempted in C++. Could you share it if you know of one? I want to create subdomains on a mesh with different properties using the FEniCS Solid Mech. app.

Thanks!

Answer 2

Hi,

Have you looked at this?

http://fenicsproject.org/documentation/tutorial/materials.html