This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

0 votes

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:

enter image description here

?

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

asked Oct 25, 2014 by Panda FEniCS Novice (440 points)
edited Oct 26, 2014 by Panda

2 Answers

+2 votes
 
Best answer

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.

answered Oct 27, 2014 by MiroK FEniCS Expert (80,920 points)
selected Oct 27, 2014 by Panda

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!

0 votes
answered Oct 27, 2014 by V_L FEniCS User (4,440 points)
...