# Extract stiffness tensor C

Hello guys !

So i wanna calculate the components of a 2D stiffness tensor and i already have all the values of the stress (sigma11,sigma22,sigma12) and strains (epsilon11,epsilon22,epsilon12). How can i do that ? Any ideas about how can we solve that ? Hi @CHAIkah,

It would help if you provide a little more detail about your problem. For example:

• Are you trying to solve a homogenization problem? (i.e. to find the effective stiffness tensor)
• What format do you have the stresses and strains in? For example, do you know the strains pointwise throughout some domain, e.g. as a FEniCS Function, or did you compute the volume average so that the strains are stored in an array?
• Are you studying plane stress or plane strain?
1 Like

Hey @conpierce8 !
Yes it’s a homogenization problem so I’m trying to extract the effective stiffness tensor.
It’s a 2D square divided into 10 phases each phase contains a different material under compression on the top and bottom boundaries.
I didn’t specify in the code if it’s plane strain/stress
It’s basically like what we have here Periodic Homogenization except that my model doesn’t have any periodic BC.

Hi @CHAIkah,

I’ll be able to help you more effectively if you provide a code example. In the meantime, note the following:

• You need to study three different load cases (uniaxial strain in x, uniaxial strain in y, and shear strain) to get the full 2D stiffness tensor. (This is the approach taken in the periodic homogenization demo.)
• If you prescribe unit strain on the boundary, the 2D stiffness tensor components are simply equal to the volume-averaged stress. You can compute the volume averaged stress using assemble as in the demo that you mentioned.

As an aside, you implicitly choose either plane strain or plane stress when you create the function sigma(v) that defines your constitutive model, as mentioned in the 2D linear elasticity demo.

1 Like

Hi @CHAIkah,

Since this isn’t a working code, I can’t provide a full solution. But, you should be able to use

sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
sigma_avg = assemble(sum(inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases))) / area


to compute the average stress sigma_avg. You should repeat this analysis three times with three different affine boundary conditions u = \bar{\varepsilon}^{(i)} x\quad \forall x \in \partial \Omega, where:

\bar{\varepsilon} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{12} \end{bmatrix} \\ \bar{\varepsilon}^{(1)} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix},\quad \bar{\varepsilon}^{(2)} = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}, \quad \bar{\varepsilon}^{(3)} = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}

Then the average stress obtained for each case will give you the corresponding column of the 2D stiffness matrix.

you mean i need to study three different load cases (uniaxial strain in x , uniaxial strain in y , and shear strain and each time i’ll get a column of the tensor ?

That’s correct. The “effective stiffness” tensor C^{eff} is defined by the relation

\left< \sigma \right> = C^{eff} \left< \varepsilon \right>

where \left<\cdot\right> denotes the volume average. When you apply an affine boundary condition, the volume average of the strain \left< \varepsilon \right> is equal to the “boundary” strain \bar{\varepsilon}. Thus,

\left< \sigma \right> = C^{eff} \bar{\varepsilon}

If we take our 2D constitutive relation to be

C^{eff} = \begin{bmatrix} C^{eff}_{11} & C^{eff}_{12} & C^{eff}_{16} \\ C^{eff}_{12} & C^{eff}_{22} & C^{eff}_{26} \\ C^{eff}_{16} & C^{eff}_{26} & C^{eff}_{66} \end{bmatrix}

and apply a unit strain in the x-direction, for example, the volume average of the stress will be:

\left< \sigma \right> = \begin{bmatrix} C^{eff}_{11} & C^{eff}_{12} & C^{eff}_{16} \\ C^{eff}_{12} & C^{eff}_{22} & C^{eff}_{26} \\ C^{eff}_{16} & C^{eff}_{26} & C^{eff}_{66} \end{bmatrix} \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} C^{eff}_{11} \\ C^{eff}_{12} \\ C^{eff}_{16} \end{bmatrix}

Repeating for \bar{\varepsilon} = [0\,1\,0]^T and \bar{\varepsilon} = [0\,0\,1]^T gives the second and third columns of the effective stiffness tensor, respectively.

i understand now. Thank you it works !
but @conpierce8 this is what i got and i don’t understand why isn’t symmetrical.

-1.572410583195420555e-03 2.216522418232252676e+03 1.108260422906548456e+03
4.414543445680684272e+03 -6.237213042085443552e-04 2.207271410980084511e+03
-6.237146402310145103e-04 1.941621083991628257e+04 9.708105108099545760e+03

@CHAIkah can you share your complete code? I suspect there may be an error with your boundary conditions, since the third column is a linear combination of the first two columns (within discretization error).

Hi @CHAIkah,

1. Your affine Dirichlet boundary conditions should be applied on all boundaries.
2. In order to directly extract the stiffness tensor components from the average stress, you need to apply a unit strain. Otherwise you will need to divide the stress by the strain to find the stiffness tensor components. To achieve a unit strain, you can apply the following affine BCs (corresponding to uniaxial strain in x, uniaxial strain in y, and pure shear strain, respectively):
1. Be careful with the ordering of your stresses. FEniCS gives you the four components of \left< \sigma \right> as a vector:

Therefore, after you

AA = np.delete(AA, 2)


the ordering of the stresses in AA is:

AA = sigma_xx
AA = sigma_xy
AA = sigma_yy


To achieve the familiar symmetric stiffness matrix, you must re-order AA (and BB and CC):

AA = np.delete(AA, 2)[[0, 2, 1]]


See the code below.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
fname = "lamin"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
#define materials
Em, num = 50e3 , 0.2
Er, nur = 210e3 , 0.3

material_parameters = [(Em, num), (Er, nur),(Em, num), (Er, nur),(Em, num),\
(Er, nur), (Em, num),(Er, nur), (Em, num),(Er, nur)]

nphases = len(material_parameters)
#define eps and sigma
def eps(v):

def sigma(v,i):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v))*Identity(2) + 2*mu*(eps(v))
# Define function space
W  = VectorFunctionSpace(mesh,"Lagrange", 2)
w  = Function(W)
v = TestFunction(W)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)

Wint = sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)])

#uniaxial strain in y solution
bc_yy = Expression(("0", "x"), degree=2)
bcs1 = [DirichletBC(W, bc_yy, facets, j) for j in [1,2,3,4]]

solve(Wint == 0, w, bcs1)
#compute eff stiffness tensor
area = 100
sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
A = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
AA =np.array(A)
AA = np.delete(AA, 2)[[0,2,1]]

#Shear strain solution
bc_xy = Expression(("x", "x"), degree=2)
bcs2 = [DirichletBC(W, bc_xy, facets, j) for j in [1,2,3,4]]
solve(Wint == 0, w, bcs2)
B = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
BB =np.array(B)
BB = np.delete(BB, 2)[[0,2,1]]

#uniaxial strain in x solution
bc_xx = Expression(("x","0"), degree=2)
bcs3 = [DirichletBC(W, bc_xx, facets, j) for j in [1,2,3,4]]

solve(Wint == 0, w, bcs3)
C = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
CC =np.array(C)
CC = np.delete(CC, 2)[[0,2,1]]

#assemble tensor columns
dd = [[m, n, l]for m,n,l in zip(CC,AA,BB)]
np.savetxt("stiffnessmatrix.txt", dd)


On my machine, this produces the following output in stiffnessmatrix.txt:

  1.526894178315140598e+05  3.272079349476299831e+04 -2.101344925113153399e-03
3.272079349479604207e+04  9.543267626480710169e+04 -4.449651732026041126e-03
-1.050665393418484094e-03 -2.224827810501039949e-03  6.996483562031936890e+04


which is symmetric to within discretization error.

2 Likes

After thinking about this problem some more, I have realized I should make one minor correction to my answer(s). The 2D strain vector is generally formulated as

\varepsilon = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \end{bmatrix} = \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \end{bmatrix}

Then the 2D elasticity relation is written as

\begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{bmatrix} = \begin{bmatrix} C_{11} & C_{12} & C_{16} \\ C_{12} & C_{22} & C_{26} \\ C_{16} & C_{26} & C_{66} \end{bmatrix} \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \gamma_{12} \end{bmatrix}

Thus, the DirichletBC for a unit shear strain should be

bc_xy = Expression(("x/2", "x/2"), degree=2)


rather than

bc_xy = Expression(("x", "x"), degree=2)


as I originally posted.

Yeah i rectified the mistake !
Thank you very much !

@conpierce8
if i wanna extract the displacement vector and strain and stress tensors in each subdomain, how can i do that from this code ?

this is what i tried to do for now:

    for i in range(nphases) :
sig11 = assemble(sigma(w,i)[0,0]*dx(i))
sig22 = assemble(sigma(w,i)[1,1]*dx(i))
sig12 = assemble(sigma(w,i)[1,0]*dx(i))
eps11 = assemble(eps(w)[0,0]*dx(i))
eps22 = assemble(eps(w)[1,1]*dx(i))
eps12 = assemble(eps(w)[1,0]*dx(i))
ux = assemble(w*dx(i))
uy = assemble(w*dx(i))


and then print them using np.array

The above code should work, but note that you need to divide by the area of the subdomain in order to obtain the correct average stresses, strains, and displacements.

1 Like

Hi @conpierce8,

Thanks so much for your answering, which really helps me out on characterizing composites by Fenics. However, I am just a beginner of Fenics and I need to specify various stiffness matrices for each element. So would you mind providing a little more details on your geometry?

Hi @ShawnZ,

I don’t have access to the original geometry since it appears the original poster has deleted some messages in this thread. Here is a working code and the corresponding meshes for a 2\times 2 square grid. The meshes were created with Gmsh but unfortunately I don’t have the .geo file.

from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
#Define geometry
# fname = "lamin"
fname = "grid"
mesh = Mesh(fname + ".xml")
subdomains = MeshFunction("size_t", mesh, fname + "_physical_region.xml")
facets = MeshFunction("size_t", mesh, fname + "_facet_region.xml")
subdomains_xdmf = XDMFFile("subdomains.xdmf")
subdomains_xdmf.write(subdomains)
subdomains_xdmf.close()
facets_xdmf = XDMFFile("facets.xdmf")
facets_xdmf.write(facets)
facets_xdmf.close()
#define materials
Em, num = 50e3 , 0.2
Er, nur = 210e3 , 0.3

# material_parameters = [(Em, num), (Er, nur),(Em, num), (Er, nur),(Em, num),\
#                        (Er, nur), (Em, num),(Er, nur), (Em, num),(Er, nur)]
n_phases = 4
nu_random = np.random.normal(0.3, 0.04, 4)
material_parameters = [(Em, nu) for nu in nu_random]

nphases = len(material_parameters)
#define eps and sigma
def eps(v):

def sigma(v,i):
E, nu = material_parameters[i]
lmbda = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
return lmbda*tr(eps(v))*Identity(2) + 2*mu*(eps(v))
# Define function space
W  = VectorFunctionSpace(mesh,"Lagrange", 2)
w  = Function(W)
v = TestFunction(W)
dx = Measure("dx", domain=mesh, subdomain_data=subdomains)

Wint = sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)])

#uniaxial strain in y solution
bc_yy = Expression(("0", "x"), degree=2)
w = project(bc_yy, W)
bcs1 = [DirichletBC(W, bc_yy, facets, j) for j in [1,2,3,4]]

solve(Wint == 0, w, bcs1)
#compute eff stiffness tensor
area = 100
sigma_fs = TensorFunctionSpace(mesh, "Real", 0)
sigma_test = TestFunction(sigma_fs)
A = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
AA =np.array(A)
AA = np.delete(AA, 2)[[0,2,1]]

# #Shear strain solution
# bc_xy = Expression(("x/2", "x/2"), degree=2)
# w = project(bc_xy, W)
# bcs2 = [DirichletBC(W, bc_xy, facets, j) for j in [1,2,3,4]]
# solve(Wint == 0, w, bcs2)
# B = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
# BB =np.array(B)
# BB = np.delete(BB, 2)[[0,2,1]]

# #uniaxial strain in x solution
# bc_xx = Expression(("x","0"), degree=2)
# w = project(bc_xx, W)
# bcs3 = [DirichletBC(W, bc_xx, facets, j) for j in [1,2,3,4]]

# solve(Wint == 0, w, bcs3)
# C = assemble(sum([inner(sigma(w, i), sigma_test)*dx(i) for i in range(nphases)]))/ area
# CC =np.array(C)
# CC = np.delete(CC, 2)[[0,2,1]]

# #assemble tensor columns
# dd = [[m, n, l]for m,n,l in zip(CC,AA,BB)]
# np.savetxt("stiffnessmatrix.txt", dd)


grid.xml:

<?xml version="1.0" encoding="UTF-8"?>

<dolfin xmlns:dolfin="http://www.fenicsproject.org">
<mesh celltype="triangle" dim="2">
<vertices size="13">
<vertex index="0" x="0.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="1" x="1.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="2" x="2.0000000000000000e+00" y="0.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="3" x="0.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="4" x="1.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="5" x="2.0000000000000000e+00" y="1.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="6" x="0.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="7" x="1.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="8" x="2.0000000000000000e+00" y="2.0000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="9" x="5.0000000000000000e-01" y="5.0000000000000000e-01" z="0.0000000000000000e+00"/>
<vertex index="10" x="1.5000000000000000e+00" y="5.0000000000000000e-01" z="0.0000000000000000e+00"/>
<vertex index="11" x="1.5000000000000000e+00" y="1.5000000000000000e+00" z="0.0000000000000000e+00"/>
<vertex index="12" x="5.0000000000000000e-01" y="1.5000000000000000e+00" z="0.0000000000000000e+00"/>
</vertices>
<cells size="16">
<triangle index="0" v0="0" v1="1" v2="9"/>
<triangle index="1" v0="3" v1="0" v2="9"/>
<triangle index="2" v0="1" v1="4" v2="9"/>
<triangle index="3" v0="4" v1="3" v2="9"/>
<triangle index="4" v0="1" v1="2" v2="10"/>
<triangle index="5" v0="4" v1="1" v2="10"/>
<triangle index="6" v0="2" v1="5" v2="10"/>
<triangle index="7" v0="5" v1="4" v2="10"/>
<triangle index="8" v0="4" v1="5" v2="11"/>
<triangle index="9" v0="7" v1="4" v2="11"/>
<triangle index="10" v0="5" v1="8" v2="11"/>
<triangle index="11" v0="8" v1="7" v2="11"/>
<triangle index="12" v0="3" v1="4" v2="12"/>
<triangle index="13" v0="6" v1="3" v2="12"/>
<triangle index="14" v0="4" v1="7" v2="12"/>
<triangle index="15" v0="7" v1="6" v2="12"/>
</cells>
</mesh>
</dolfin>


grid_facet_region.xml:

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_function type="uint" dim="1" size="28">
<entity index="0" value="1"/>
<entity index="1" value="4"/>
<entity index="2" value="0"/>
<entity index="3" value="1"/>
<entity index="4" value="0"/>
<entity index="5" value="0"/>
<entity index="6" value="0"/>
<entity index="7" value="2"/>
<entity index="8" value="0"/>
<entity index="9" value="0"/>
<entity index="10" value="4"/>
<entity index="11" value="0"/>
<entity index="12" value="0"/>
<entity index="13" value="0"/>
<entity index="14" value="0"/>
<entity index="15" value="0"/>
<entity index="16" value="0"/>
<entity index="17" value="0"/>
<entity index="18" value="0"/>
<entity index="19" value="2"/>
<entity index="20" value="0"/>
<entity index="21" value="0"/>
<entity index="22" value="3"/>
<entity index="23" value="0"/>
<entity index="24" value="3"/>
<entity index="25" value="0"/>
<entity index="26" value="0"/>
<entity index="27" value="0"/>
</mesh_function>
</dolfin>


grid_physical_region.xml:

<?xml version="1.0" encoding="UTF-8"?>
<dolfin xmlns:dolfin="http://fenicsproject.org">
<mesh_function type="uint" dim="2" size="16">
<entity index="0" value="1"/>
<entity index="1" value="1"/>
<entity index="2" value="1"/>
<entity index="3" value="1"/>
<entity index="4" value="2"/>
<entity index="5" value="2"/>
<entity index="6" value="2"/>
<entity index="7" value="2"/>
<entity index="8" value="3"/>
<entity index="9" value="3"/>
<entity index="10" value="3"/>
<entity index="11" value="3"/>
<entity index="12" value="4"/>
<entity index="13" value="4"/>
<entity index="14" value="4"/>
<entity index="15" value="4"/>
</mesh_function>
</dolfin>

1 Like