# How to assign different material properties to elements (Linear Elasticity)

Hi

I want to assign different material properties to different elements in a simple linear elasticity problem. The domain is a cube, fixed from left face and a load is applied on the right face. I have 6 tetrahedron elements in this cube. The index of those elements starts from 0 to 5 that could be printed by:

for cell in cells(mesh):
print (cell.index())

Now I want to assign different Young’s modulus to each element. With that being said, I want to assign : E = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]
to the element_indexes = [0 , 1 , 2 , 3 , 4 , 5] respectively.
I was not able to find a clear example on this. I would appreciate any hint clarifying how to do this in FEniCS. The below code solves this linear elasticity problem when the Young’s modulus is constant over all elements (e.g. 1E6):

from dolfin import *

# Define boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

for cell in cells(mesh):
print (cell.index())

n = FacetNormal(mesh)

V = VectorFunctionSpace(mesh, 'CG', degree=1)

left = Left()
right = Right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

E = 1E6  # I want to use 6 different E for 6 elements : E = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]
nu =0.3
force = 1000.
mu = 0.4e6

lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

def sigma(v):
return lmbda * tr(eps(v)) * Identity(3) + 2.0 * mu * eps(v)

du = TrialFunction(V)
u_ = TestFunction(V)

a = inner(sigma(du), eps(u_)) * dx
l =  inner(force * n, u_) * ds(2)

u = Function(V)
solve(a == l, u, DirichletBC(V, Constant((0., 0.,0.)), boundaries, 1))

Hi,

How about using a DG0 function to assign the values to each dof (through a dof-to-cell map) ?

# E =   # I want to use 6 different E for 6 elements : E = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]
DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
dofmap = DG.dofmap()
coords = DG.tabulate_dof_coordinates()
dof_to_cell = [cell for cell in range(mesh.num_cells())
for dof in dofmap.cell_dofs(cell)]
for i in range(coords.shape[0]):
cell_index = dof_to_cell[i]
E.vector()[i] = 1.e6 * (cell_index + 1)

There could be multiple ways of doing this I feel.

Edit 1: I had come accross MiroK’s answer on the old forum sometime back for something in my work.

Edit 2: Also this answer by @dokken

1 Like

The output of the for loop in your code, produces the list of the Young’s modulus values that I could print by: print (E.vector()[i])
But in this case the values of the Young’s modulus are already known as I put them in my code: E = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]
So I guess this part of your snippets code is not necessary. I agree to use DG0 space. The issue is how to write the variational form based on this. I updated the code for DG0 space. Of course it does not work yet as I am not sure how to handle the stress term in variational form based on different Young’s modulus for each element. What changes should I do to get it to work? Here is the code :

from dolfin import *

# Define boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

n = FacetNormal(mesh)

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
dofmap = DG.dofmap()
coords = DG.tabulate_dof_coordinates()
dof_to_cell = [cell for cell in range(mesh.num_cells())
for dof in dofmap.cell_dofs(cell)]

left = Left()
right = Right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

E = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]

nu =0.3
force = 1000.
mu = 0.4e6

# This is dependent to E which is an array. Of course it does not make sense to define it this way!
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

def sigma(v):
return lmbda * tr(eps(v)) * Identity(3) + 2.0 * mu * eps(v)

du = TrialFunction(DG)
u_ = TestFunction(DG)

a = inner(sigma(du), eps(u_)) * dx
l =  inner(force * n, u_) * ds(2)

u = Function(DG)
solve(a == l, u, DirichletBC(DG, Constant((0., 0.,0.)), boundaries, 1))

Right, I meant to point out cell-wise assignment of material properties for a more general case, if needed. For a simple mesh such as this, doing

E.vector()[:] = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]

should suffice.

I guess I meant using DG0 for your material properties and not the unknown. The displacement field should be piecewise continuous (CG1), so I would change the relevant portions to

V = VectorFunctionSpace(mesh, "CG", 1)
du = TrialFunction(V)
u_ = TestFunction(V)

a = inner(sigma(du), eps(u_)) * dx
l =  inner(force * n, u_) * ds(2)

u = Function(V)
solve(a == l, u, DirichletBC(V, Constant((0., 0.,0.)), boundaries, 1))

Also if your Young’s modulus is changing, and the poisson’s ratio is constant, you should consider changing mu as well:

nu = 0.3
force = 1000.
mu = project(E/2./(1.+nu),DG)
# This is dependent to E which is an array. Of course it does not make sense to define it this way!
lmbda = E*nu/(1+nu)/(1-2*nu)
lmbda = project(lmbda, DG)

Having defined the material properties the way you want them to be, I think the variational formulation should be unchanged from thereon.

2 Likes

Perfect. So the below code seems to work:

from dolfin import *

# Define boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 0.0) < DOLFIN_EPS and on_boundary

class Right(SubDomain):
def inside(self, x, on_boundary):
return abs(x[0] - 1.0) < DOLFIN_EPS and on_boundary

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

n = FacetNormal(mesh)

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)

E.vector()[:] = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]

left = Left()
right = Right()

boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())

boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

nu = 0.3
force = 1000.
mu = project(E/2./(1.+nu),DG)
lmbda = E*nu/(1+nu)/(1-2*nu)
lmbda = project(lmbda, DG)

def eps(v):

def sigma(v):
return lmbda * tr(eps(v)) * Identity(3) + 2.0 * mu * eps(v)

V = VectorFunctionSpace(mesh, "CG", 1)
du = TrialFunction(V)
u_ = TestFunction(V)

a = inner(sigma(du), eps(u_)) * dx
l =  inner(force * n, u_) * ds(2)

u = Function(V)
solve(a == l, u, DirichletBC(V, Constant((0., 0.,0.)), boundaries, 1))

In the above code, when we define :

E.vector()[:] = [1E6 , 2E6 , 3E6, 4E6 , 5E6 , 6E6]

How do we know which Young’s modulus is assigned to which element? Does it mean that for example the first one (1E6) goes to element.index = 0 and the last one (6E6) goes to element.index = 5?

So the order of dofs of E is the same as DG.tabulate_dof_coordinates(). This was one of the reasons I had a dof_to_cell array to map the dofs to the corresponding cells and assign E the value you want in that cell.

coords = DG.tabulate_dof_coordinates()
dof_to_cell = [cell for cell in range(mesh.num_cells())
for dof in dofmap.cell_dofs(cell)]
for i in range(coords.shape[0]):
cell_index = dof_to_cell[i] # the i-th dof corresponds to the cell_index-th cell
E.vector()[i] = 1.e6 * (cell_index + 1)
2 Likes

Cool. So if I have a list of Young’s modulus with random values like :

E_VALUES = [7.5E6, 12E6, 83E6, 42E6,33E6,98E6]

Does the below code assign the numbers in the above list to the element index (from 0 to 5) respectively?

from dolfin import *

E_VALUES = [7.5E6, 12E6, 83E6, 42E6,33E6,98E6]
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
dofmap = DG.dofmap()
coords = DG.tabulate_dof_coordinates()
dof_to_cell = [cell for cell in range(mesh.num_cells())
for dof in dofmap.cell_dofs(cell)]
for i in range(coords.shape[0]):
cell_index = dof_to_cell[i]
print (cell_index)
E.vector()[i] = E_VALUES[i] * (cell_index + 1)/(cell_index + 1)

Just doing this should suffice. In fact now that I think, for a piecewise constant function (DG0) the order of dofs (i) and cells (cell_index) should be the same. So I think for this case, you can also do:

from dolfin import *
E_VALUES = [7.5E6, 12E6, 83E6, 42E6,33E6,98E6]
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
E.vector()[:] = E_VALUES

So the question is, if we ignore the cell_index, what is the point in doing: dof_to_cell?
I just double checked that with and without dof_to_cell block and both returned the same results. I mean instead of :

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
E_VALUES = [1E6, 2E6, 3E6, 4E6,5E6,6E6]
dofmap = DG.dofmap()
coords = DG.tabulate_dof_coordinates()
dof_to_cell = [cell for cell in range(mesh.num_cells())
for dof in dofmap.cell_dofs(cell)]
for i in range(coords.shape[0]):
cell_index = dof_to_cell[i]
E.vector()[i] = E_VALUES[i]

We can simply replace:

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
E.vector()[:] = [1E6, 2E6, 3E6, 4E6,5E6,6E6]

They both give the same results. What are your thoughts about it?

Yes, you’re right. In fact that’s what I meant above.

However this wouldn’t work straight away for DG1 fields and so on. The dofs will still be arranged in the order of cells, but now with more than one dof per cell.

Makes sense. Anyways, thanks for your help.

Hello
I have 1 more question. What if I have 2 different domains and I only want to apply this DG0 space on 1 of them? For example lets say we have 2 cubes:

There are 12 elements in above figure. 6 for the blue one and 6 for red one.
Now I want to apply the DG0 space on the right cube (e.g. Red cube) and solve the blue one with CG1.
First I found the index of the elements in the red cube as: [1 , 5 ,6 , 8 , 9,11]
Now if I do :

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
E.vector()[:] = [1E6, 0.3E6, 5E6, 3E6,2E6,0.4E6]

The question is:
Does the order of the E values still follow the order of the element indexes?
Does it mean that for example 1E6 is assigned to index 1 , or 5E6 assigned to index 6? I am asking because this time the indexes does not start from 0.
Thanks!

Hi Leo,
I am not sure if I understand your question properly. You want to interpolate E using a CG1 function for one part of the mesh, and DG0 otherwise? In that case, how do you wish to assign the values for the CG1 part? For the snippet that you have provided 1E6 gets assigned to index 0 (the first dof), which corresponds to the first (global) cell, since E is a function on the whole mesh. Also E.vector() should have a shape (12,) instead of 6.

Thanks bhavesh.
If I want to be more specific, for the blue cube, I just want to solve it using a regular linear elasticity (Constant E over all cells). I mean I do not need to use DG0 method. But I want to assign different E to the elements of the red cube. The problem is that lets say the index 0 corresponds to an element in the blue cube. There are 6 elements in each cube and lets say the indexes of the elements in the red cube are [1 , 5 ,6 , 8 , 9,11] . (These are global indexes of the elements in the red cube)
My question is: If I use: E.vector()[:] = [1E6, 0.3E6, 5E6, 3E6,2E6,0.4E6]
The E values will sequentially be assigned to the global indexes : [1 , 5 ,6 , 8 , 9,11] ?

I see.

I don’t think this is true, because the shape of the array on left is 12 whereas for the list on right, it’s 6. Maybe a better way to put it is via the following:

from dolfin import *
import numpy as np
mesh = UnitCubeMesh(2,1,1)
V = FunctionSpace(mesh, "DG", 0)
E = Function(V, name="E")
indices_on_left = np.array([
i for i,cell in enumerate(cells(mesh)) if cell.midpoint().x() < 0.5
], int)
indices_on_right = np.setdiff1d(np.arange(mesh.num_cells()), indices_on_left)
E.vector()[indices_on_left] = 1. #constant
E.vector()[indices_on_right] = np.random.random(indices_on_right.shape) + 1. #chosenby you

This should assign the values on the left and right subdomains (without actually having to create instances of subdomains). Is this what you are looking for?

Cool. Thanks for your suggestion.
I think it is better to forget about parameter E so I can explain what I am looking for in a better way. I want to solve a simple regular linear elasticity problem on the left subdomain (This should be done easily with knowing the values of the E and \nu parameters)

When it comes to the right subdomain,there are 6 elements in there. Lets say the global indexes of the elements in that subdomain are:

right_subdomain_elem_index = [2 , 5 ,6 , 8 , 9, 11]

Please note that the indexes does NOT necessarily start from 0.
In this subdomain, I have new variational form I have a parameter called \alpha and I know the values of \alpha for each element as :

alpha_values = [10 , 6 , 8 , 15 , 3 , 7]

What I am looking for, is to assign the alpha values to the elements in the right subdomain. For example alpha = 10 should be assigned to index = 2, alpha = 6 should be assigned to index = 5, etc. With that being said, we need DG0 on the right subdomain.
Again the \alpha only exists on the right subdomain.
What is the best way to handle such a problem?
Thanks again!

Ok. Could you create a MWE?
If that is the case and your entire mesh comprises of the left and right subdomains, with \alpha only existing on one of them, then one way that I can see is by creating a SubMesh, so effectively doing

from dolfin import *
mesh = UnitCubeMesh(2,1,1)
rightMesh = SubMesh(mesh, CompiledSubDomain("x[0] >= 0.5"))
Vr = FunctionSpace(rightMesh, "DG", 0)
alph = Function(Vr)
alpha_values = [10 , 6 , 8 , 15 , 3 , 7]
alph.vector()[:] = alpha_values #this assigns the value in the order of rightmesh.cells()

This should assign the values in the order of rightMesh.cells() which should map to global indices rightMesh.data().array("parent_cell_indices", mesh.geometric_dimension()).

Alright. I just created a MWE. The same geometry. There are 2 cubes. The left cube should be solved using linear elasticity and the right one should be solved using hyperelasticity. The left face of the left cube is fixed (Marked as 1) and a load is applied the right face of the right cube (Marked as 2). The mesh is structured (6 elements on each cube). I created the mesh in GMSH and imported into the FEniCS. Here is the GMSH file:

Point(1) = {0, 0, 0,1};
Point(2) = {1, 0, 0,1};
Point(3) = {2, 0, 0,1};
Point(4) = {0, 1, 0,1};
Point(5) = {1, 1, 0,1};
Point(6) = {2, 1, 0,1};
Point(7) = {0, 0, 1,1};
Point(8) = {1, 0, 1,1};
Point(9) = {2, 0, 1,1};
Point(10) = {0, 1, 1,1};
Point(11) = {1, 1, 1,1};
Point(12) = {2, 1, 1,1};

Line(1) = {10, 11};
Line(2) = {12, 11};
Line(3) = {10, 4};
Line(4) = {4, 1};
Line(5) = {1, 7};
Line(6) = {7, 10};
Line(7) = {5, 11};
Line(8) = {6, 12};
Line(9) = {6, 5};
Line(10) = {5, 4};
Line(11) = {1, 2};
Line(12) = {2, 3};
Line(13) = {3, 9};
Line(14) = {9, 12};
Line(15) = {9, 8};
Line(16) = {8, 7};
Line(17) = {8, 2};
Line(18) = {8, 11};
Line(19) = {6, 3};
Line(20) = {2, 5};

Line Loop(21) = {3, -10, 7, -1};
Plane Surface(22) = {21};
Line Loop(23) = {3, 4, 5, 6};
Plane Surface(24) = {23};
Line Loop(25) = {1, -18, 16, 6};
Plane Surface(26) = {25};
Line Loop(27) = {4, 11, 20, 10};
Plane Surface(28) = {27};
Line Loop(29) = {11, -17, 16, -5};
Plane Surface(30) = {29};
Line Loop(31) = {18, -7, -20, -17};
Plane Surface(32) = {31};
Line Loop(33) = {18, -2, -14, 15};
Plane Surface(34) = {33};
Line Loop(35) = {7, -2, -8, 9};
Plane Surface(36) = {35};
Line Loop(37) = {20, -9, 19, -12};
Plane Surface(38) = {37};
Line Loop(39) = {15, 17, 12, 13};
Plane Surface(40) = {39};
Line Loop(41) = {14, -8, 19, 13};
Plane Surface(42) = {41};
Surface Loop(43) = {22, 24, 28, 30, 26, 32};
Volume(44) = {43};
Surface Loop(45) = {36, 34, 42, 38, 40, 32};
Volume(46) = {45};
Transfinite Surface {22,24,26,28,30,32,34,36,38,40,42};
Transfinite Line {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20} = 2;
Physical Volume(1) = {44};
Physical Volume(2) = {46};
Physical Surface(1) = {24};
Physical Surface(2) = {42};

Now the parameter \alpha ONLY exists on the variational form on the right subdomain and I want to assign 6 different values of \alpha to these elements. The global indexes of the elements on the right cube are: [6,7,8,9,10,11]
Here is the implementation based on previous suggestions:

from dolfin import *

mesh = Mesh("MESH.xml")

boundaries = MeshFunction("size_t", mesh, "MESH_facet_region.xml")
subdomains = MeshFunction("size_t", mesh, "MESH_physical_region.xml")

n = FacetNormal(mesh)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains, metadata={'quadrature_degree': 2})
ds = Measure('ds', domain=mesh, subdomain_data=boundaries, metadata={'quadrature_degree': 2})

V = VectorFunctionSpace(mesh,"Lagrange",degree=2)
u = Function(V)
v = TestFunction(V)

nu = 0.3
force = 1000.
E = 1E6
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)

def eps(v):

def sigma(v):
return lmbda * tr(eps(v)) * Identity(3) + 2.0 * mu * eps(v)

bc = DirichletBC(V, Constant((0.,0.,0)), boundaries , 1)

F_lin = inner(sigma(u), eps(v))*dx(1)  # left subdomain is solved using linear elastcity

# Right subdomain is solved using nonlinear hyperelasticity
I = Identity(3)
F = I + grad(u)
C = F.T*F
Ic = tr(C)
J  = det(F)

DG = FunctionSpace(mesh, "DG", 0)  # The parameter alpha ONLY exists on the right subdomain
alpha = Function(DG)
alpha.vector()[:] = [1000 , 200 , 700, 100, 2000 , 300] # Assigning different values of alpha to 6 elements on the right subdomain
alpha = project(alpha, DG)

lmbda_2 = 5E6

psi = (alpha/2)*(Ic - 3) - alpha*ln(J) + (lmbda_2/2)*(ln(J))**2
Pi = psi*dx(2) - inner(-force*n, u)*ds(2)

F_hyp = derivative(Pi, u, v)
F = F_lin + F_hyp
J = derivative(F, u)
solve(F==0, u, bc, J=J)

This returns this error:

*** Error:   Unable to set local values of PETSc vector.
*** Reason:  Size of values array is not equal to local vector size.

It seems like FEniCS expects the alpha to exist in all elements (Both subdomains) , while it only appears on the variational form on the right subdomain.
I hope I was able to show what I am trying to do. Do you see any solution to handle such problem?
Thanks!

Hi Leo,

What if you let alpha be defined on the entire domain and only set the values for the elements on the red cube by extracting the corresponding indices (which infact you already have) using

alpha.vector()[[6,7,8,9,10,11]] = [1000 , 200 , 700, 100, 2000 , 300] # Assigning different values of alpha to 6 elements on the right subdomain

Infact as long as you have the indices of the dofs (which in this case are also the cell indices), you can set those manually. And since you are assembling the corresponding form (hyperelasticity) only on the right domain, this should take care of the error you currently have. The rest looks good to me upto the level of solver (that might require tweaking for your real problem).

If not, you can also create a SubMesh as I mentioned before. I believe this should give you identical solution. I can test this sometime tomorrow (with the actual mesh). So the only changes that I make to your code, for now,

import numpy as np
from dolfin import *

snes_solver_parameters = {"nonlinear_solver": "snes",
"snes_solver": {"linear_solver": "lu",
"maximum_iterations": 100,
"report": True,
"line_search": 'basic',
"error_on_nonconvergence": False,
"relative_tolerance":1e-7,
"absolute_tolerance":1e-8}}

mesh = UnitCubeMesh(2,1,1)
left = CompiledSubDomain("x[0] <= 0.5")
right = CompiledSubDomain("x[0] >= 0.5")
leftFace = CompiledSubDomain("near(x[0], 0)")
rightFace = CompiledSubDomain("near(x[0], 1)")

subdomainMark = cpp.mesh.MeshFunctionSizet(mesh, mesh.topology().dim())
faceMark = cpp.mesh.MeshFunctionSizet(mesh, mesh.topology().dim()-1)

rightFace.mark(faceMark, 2)
left.mark(subdomainMark, 1)
right.mark(subdomainMark, 2)

leftMesh = SubMesh(mesh, left)
rightMesh = SubMesh(mesh, right)

n = FacetNormal(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=subdomainMark)
ds = Measure("ds", domain=mesh, subdomain_data=faceMark)
# ...
# ...
indices_on_left = np.array([
i for i,cell in enumerate(cells(mesh)) if cell.midpoint().x() < 0.5
], int)
indices_on_right = np.setdiff1d(np.arange(mesh.num_cells()), indices_on_left)
# ....
Vr = FunctionSpace(mesh, "DG", 0)
alpha = Function(Vr, name="alpha")
alpha.vector()[indices_on_right] = [1000 , 200 , 700, 100, 2000 , 300] # Assigning different values of alpha to 6 elements on the right subdomain

# ...
J = derivative(F, u, TrialFunction(V))
problem = NonlinearVariationalProblem(F, u, bcs=bc, J=J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
1 Like

Awesome. Thanks again for your clarification