#
=========================================================================
#
Python
implementation:
Auxiliary
function
to
rotate
the
cubic
elasticity
#
matrix
around
an
arbitrary
axis
(
after
https://github.com
/
MikulaJakub)
#
J.
Ulloa,
M.P.
Ariza,
J.E.
Andrade
&
M.
Ortiz
#
Homogenized
models
of
mechanical
metamaterials
#
Computer
Methods
in
Applied
Mechanics
and
Engineering,
August
2024
#
=========================================================================
from
dolfin
import
*
import
numpy
as
np
def
rotate_elast
(
C11
,
C12
,
C44
,
o1
,
o2
,
o3
,
ang
):
#
Cubic
elasticity
tensor
in
matrix
form
D
=
np
.
array
([
[C11
,
C12
,
C12
,
0
,
0
,
0
]
,
[C12
,
C11
,
C12
,
0
,
0
,
0
]
,
[C12
,
C12
,
C11
,
0
,
0
,
0
]
,
[
0
,
0
,
0
,
C44
,
0
,
0
]
,
[
0
,
0
,
0
,
0
,
C44
,
0
]
,
[
0
,
0
,
0
,
0
,
0
,
C44]])
#
Conversion
matrix
from
engineering
to
tensorial
strain
R
=
np
.
zeros
((
6
,
6
))
R[
0
,
0
]
=
1
.
0
R[
1
,
1
]
=
1
.
0
R[
2
,
2
]
=
1
.
0
R[
3
,
3
]
=
2
.
0
R[
4
,
4
]
=
2
.
0
R[
5
,
5
]
=
2
.
0
invR
=
R
.
copy
()
invR[
3
,
3
]
=
0
.
5
invR[
4
,
4
]
=
0
.
5
invR[
5
,
5
]
=
0
.
5
#
Rotation
matrix
about
an
arbitrary
axis
R_cg
=
np
.
array
([
[o1
**
2
+
(
1
-
o1
**
2
)
*
np
.
cos
(ang)
,
o1
*
o2
*
(
1
-
np
.
cos
(ang))
-
o3
*
np
.
sin
(ang)
,
o1
*
o3
*
(
1
-
np
.
cos
(ang))
+
o2
*
np
.
sin
(ang)]
,
[o1
*
o2
*
(
1
-
np
.
cos
(ang))
+
o3
*
np
.
sin
(ang)
,
o2
**
2
+
(
1
-
o2
**
2
)
*
np
.
cos
(ang)
,
o2
*
o3
*
(
1
-
np
.
cos
(ang))
-
o1
*
np
.
sin
(ang)]
,
[o1
*
o3
*
(
1
-
np
.
cos
(ang))
-
o2
*
np
.
sin
(ang)
,
o2
*
o3
*
(
1
-
np
.
cos
(ang))
+
o1
*
np
.
sin
(ang)
,
o3
**
2
+
(
1
-
o3
**
2
)
*
np
.
cos
(ang)]])
K1
=
np
.
zeros
((
3
,
3
))
K2
=
np
.
zeros
((
3
,
3
))
K3
=
np
.
zeros
((
3
,
3
))
K4
=
np
.
zeros
((
3
,
3
))
for
i
in
range
(
3
):
for
j
in
range
(
3
):
K1[i
,
j]
=
R_cg[i
,
j]
**
2
K2[i
,
j]
=
R_cg[i
,
(j
+
1
)
%
3
]
*
R_cg[i
,
(j
+
2
)
%
3
]
K3[i
,
j]
=
R_cg[(i
+
1
)
%
3
,
j]
*
R_cg[(i
+
2
)
%
3
,
j]
K4[i
,
j]
=
R_cg[(i
+
1
)
%
3
,
(j
+
1
)
%
3
]
*
R_cg[(i
+
2
)
%
3
,
(j
+
2
)
%
3
]
\
+
R_cg[(i
+
1
)
%
3
,
(j
+
2
)
%
3
]
*
R_cg[(i
+
2
)
%
3
,
(j
+
1
)
%
3
]
#
Populate
T_cg
and
invT_cg
T_cg
=
np
.
zeros
((
6
,
6
))
T_cg[
:
3
,
:
3
]
=
K1
T_cg[
:
3
,
3
:
6
]
=
2
.
0
*
K2
T_cg[
3
:
6
,
:
3
]
=
K3
T_cg[
3
:
6
,
3
:
6
]
=
K4
invT_cg
=
np
.
zeros
((
6
,
6
))
invT_cg[
:
3
,
:
3
]
=
K1
.
T
invT_cg[
:
3
,
3
:
6
]
=
2
.
0
*
K3
.
T
invT_cg[
3
:
6
,
:
3
]
=
K2
.
T
invT_cg[
3
:
6
,
3
:
6
]
=
K4
.
T
#
Rotated
cubic
elasticity
tensor
in
matrix
form
C_rot
=
T_cg
@
D
@
R
@
invT_cg
@
invR
return
C_rot