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 ?

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,

- Your affine Dirichlet boundary conditions should be applied on all boundaries.
- 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):

u = [x\quad 0]^T\\
u = [0\quad y]^T\\
u = [y\quad x]^T

- Be careful with the ordering of your stresses. FEniCS gives you the four components of \left< \sigma \right> as a vector:

\left< \sigma \right> = [\sigma_{11} \quad \sigma_{12} \quad \sigma_{21} \quad \sigma_{22}]^T

Therefore, after you

```
AA = np.delete(AA, 2)
```

the ordering of the stresses in `AA`

is:

```
AA[0] = sigma_xx
AA[1] = sigma_xy
AA[2] = 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):
return sym(grad(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[1]"), 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[1]", "x[0]"), 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]","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[1]/2", "x[0]/2"), degree=2)
```

rather than

```
bc_xy = Expression(("x[1]", "x[0]"), 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[0]*dx(i))
uy = assemble(w[1]*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):
return sym(grad(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[1]"), 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[1]/2", "x[0]/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]","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

Thanks so much for your answer! @conpierce8

These really help me out. I took some time to generate the corresponding `.geo`

successfully. However, I just went over the codes again and got two new questions:

(1) `area = 100`

which I suppose is pre-defined by scenario, but will it be related to the mesh size? E.g., for this solution, if taking a 4 by 4 mesh, shall we also change the area?

(2) In the file `grid_physical_region.xml`

, four physical regions (1, 2, 3, 4) have been defined. But I am just confused about the line `sum ([inner(sigma(v, i), grad(w))*dx(i) for i in range(nphases)])`

. `range(nphases)`

will iterate as 0, 1, 2, 3 for sure, which sort of conflicts to the sequences of physical region. My feeling is to fix the `dx`

but cannot make sure as a beginner of FEniCS.

Please pardon me if there are any confusing descriptions or format mistakes.

- You can use
`area = assemble(1 * dx, domain=mesh)`

to compute the area of the mesh. - You are correct, it should be

since the domains are marked [1, 2, 3, 4], not [0, 1, 2, 3].`Wint = sum ([inner(sigma(v, i), grad(w))*dx(i + 1) for i in range(nphases)])`