Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79605189
direct_comparison_small_strain.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Tue, Aug 27, 08:07
Size
8 KB
Mime Type
text/x-python
Expires
Thu, Aug 29, 08:07 (2 d)
Engine
blob
Format
Raw Data
Handle
20210972
Attached To
rMUSPECTRE µSpectre
direct_comparison_small_strain.py
View Options
## Most of the following is copied from https://github.com/tdegeus/GooseFFT (MIT license)
import
numpy
as
np
import
scipy.sparse.linalg
as
sp
import
itertools
import
os
import
sys
sys
.
path
.
append
(
os
.
path
.
join
(
os
.
getcwd
(),
"language_bindings/python"
))
import
muSpectre
as
µ
# ----------------------------------- GRID ------------------------------------
ndim
=
2
# number of dimensions
N
=
31
#31 # number of voxels (assumed equal for all directions)
offset
=
3
#9
ndof
=
ndim
**
2
*
N
**
2
# number of degrees-of-freedom
cell
=
µ
.
SystemFactory
(
µ
.
get_2d_cube
(
N
),
µ
.
get_2d_cube
(
1.
),
µ
.
Formulation
.
small_strain
)
# ---------------------- PROJECTION, TENSORS, OPERATIONS ----------------------
# tensor operations/products: np.einsum enables index notation, avoiding loops
# e.g. ddot42 performs $C_ij = A_ijkl B_lk$ for the entire grid
trans2
=
lambda
A2
:
np
.
einsum
(
'ij... ->ji... '
,
A2
)
ddot42
=
lambda
A4
,
B2
:
np
.
einsum
(
'ijkl...,lk... ->ij... '
,
A4
,
B2
)
ddot44
=
lambda
A4
,
B4
:
np
.
einsum
(
'ijkl...,lkmn...->ijmn...'
,
A4
,
B4
)
dot22
=
lambda
A2
,
B2
:
np
.
einsum
(
'ij... ,jk... ->ik... '
,
A2
,
B2
)
dot24
=
lambda
A2
,
B4
:
np
.
einsum
(
'ij... ,jkmn...->ikmn...'
,
A2
,
B4
)
dot42
=
lambda
A4
,
B2
:
np
.
einsum
(
'ijkl...,lm... ->ijkm...'
,
A4
,
B2
)
dyad22
=
lambda
A2
,
B2
:
np
.
einsum
(
'ij... ,kl... ->ijkl...'
,
A2
,
B2
)
shape
=
tuple
((
N
for
_
in
range
(
ndim
)))
# identity tensor [single tensor]
i
=
np
.
eye
(
ndim
)
def
expand
(
arr
):
new_shape
=
(
np
.
prod
(
arr
.
shape
),
np
.
prod
(
shape
))
ret_arr
=
np
.
zeros
(
new_shape
)
ret_arr
[:]
=
arr
.
reshape
(
-
1
)[:,
np
.
newaxis
]
return
ret_arr
.
reshape
((
*
arr
.
shape
,
*
shape
))
# identity tensors [grid of tensors]
I
=
expand
(
i
)
I4
=
expand
(
np
.
einsum
(
'il,jk'
,
i
,
i
))
I4rt
=
expand
(
np
.
einsum
(
'ik,jl'
,
i
,
i
))
I4s
=
(
I4
+
I4rt
)
/
2.
II
=
dyad22
(
I
,
I
)
# projection operator [grid of tensors]
# NB can be vectorized (faster, less readable), see: "elasto-plasticity.py"
# - support function / look-up list / zero initialize
delta
=
lambda
i
,
j
:
np
.
float
(
i
==
j
)
# Dirac delta function
freq
=
np
.
arange
(
-
(
N
-
1
)
/
2.
,
+
(
N
+
1
)
/
2.
)
# coordinate axis -> freq. axis
Ghat4
=
np
.
zeros
([
ndim
,
ndim
,
ndim
,
ndim
,
N
,
N
])
# zero initialize
# - compute
for
i
,
j
,
l
,
m
in
itertools
.
product
(
range
(
ndim
),
repeat
=
4
):
for
x
,
y
in
itertools
.
product
(
range
(
N
),
repeat
=
ndim
):
q
=
np
.
array
([
freq
[
x
],
freq
[
y
]])
# frequency vector
if
not
q
.
dot
(
q
)
==
0
:
# zero freq. -> mean
Ghat4
[
i
,
j
,
l
,
m
,
x
,
y
]
=
-
(
q
[
i
]
*
q
[
j
]
*
q
[
l
]
*
q
[
m
])
/
(
q
.
dot
(
q
))
**
2
+
\
(
delta
(
j
,
l
)
*
q
[
i
]
*
q
[
m
]
+
delta
(
j
,
m
)
*
q
[
i
]
*
q
[
l
]
+
\
delta
(
i
,
l
)
*
q
[
j
]
*
q
[
m
]
+
delta
(
i
,
m
)
*
q
[
j
]
*
q
[
l
])
/
(
2.
*
q
.
dot
(
q
))
# (inverse) Fourier transform (for each tensor component in each direction)
fft
=
lambda
x
:
np
.
fft
.
fftshift
(
np
.
fft
.
fftn
(
np
.
fft
.
ifftshift
(
x
),
shape
))
ifft
=
lambda
x
:
np
.
fft
.
fftshift
(
np
.
fft
.
ifftn
(
np
.
fft
.
ifftshift
(
x
),
shape
))
# functions for the projection 'G', and the product 'G : K : eps'
G
=
lambda
A2
:
np
.
real
(
ifft
(
ddot42
(
Ghat4
,
fft
(
A2
))
)
)
.
reshape
(
-
1
)
K_deps
=
lambda
depsm
:
ddot42
(
K4
,
depsm
.
reshape
(
ndim
,
ndim
,
N
,
N
))
G_K_deps
=
lambda
depsm
:
G
(
K_deps
(
depsm
))
# ------------------- PROBLEM DEFINITION / CONSTITIVE MODEL -------------------
# phase indicator: cubical inclusion of volume fraction (9**3)/(31**3)
E2
,
E1
=
75e10
,
70e9
poisson
=
.
33
hard
=
µ
.
material
.
MaterialHooke2d
.
make
(
cell
,
"hard"
,
E2
,
poisson
)
soft
=
µ
.
material
.
MaterialHooke2d
.
make
(
cell
,
"soft"
,
E1
,
poisson
)
#for pixel in cell:
# if ((pixel[0] >= N-offset) and
# (pixel[1] < offset)):
# hard.add_pixel(pixel)
# else:
# soft.add_pixel(pixel)
#
phase
=
np
.
zeros
(
shape
);
phase
[
-
offset
:,:
offset
,]
=
1.
phase
=
np
.
zeros
(
shape
);
phase
[
0
,:]
=
1.
phase
=
np
.
zeros
(
shape
);
for
i
,
j
in
itertools
.
product
(
range
(
N
),
repeat
=
ndim
):
c
=
N
//
2
if
(
i
-
c
)
**
2
+
(
j
-
c
)
**
2
<
(
N
//
5
)
**
2
:
#if j<10:
phase
[
i
,
j
]
=
1.
hard
.
add_pixel
([
i
,
j
])
else
:
soft
.
add_pixel
([
i
,
j
])
# material parameters + function to convert to grid of scalars
param
=
lambda
M0
,
M1
:
M0
*
np
.
ones
(
shape
)
*
(
1.
-
phase
)
+
M1
*
np
.
ones
(
shape
)
*
phase
# K = param(0.833,8.33) # bulk modulus [grid of scalars]
# mu = param(0.386,3.86) # shear modulus [grid of scalars]
K2
,
K1
=
(
E
/
(
3
*
(
1
-
2
*
poisson
))
for
E
in
(
E2
,
E1
))
m2
,
m1
=
(
E
/
(
2
*
(
1
+
poisson
))
for
E
in
(
E2
,
E1
))
K
=
param
(
K1
,
K2
)
mu
=
param
(
m1
,
m2
)
# stiffness tensor [grid of tensors]
K4
=
K
*
II
+
2.
*
mu
*
(
I4s
-
1.
/
3.
*
II
)
# ----------------------------- NEWTON ITERATIONS -----------------------------
# initialize stress and strain tensor [grid of tensors]
sig
=
np
.
zeros
([
ndim
,
ndim
,
N
,
N
])
eps
=
np
.
zeros
([
ndim
,
ndim
,
N
,
N
])
# set macroscopic loading
DE
=
np
.
zeros
([
ndim
,
ndim
,
N
,
N
])
DE
[
0
,
1
]
+=
0.01
DE
[
1
,
0
]
+=
0.01
delEps0
=
DE
[:
ndim
,
:
ndim
,
0
,
0
]
µ
DE
=
np
.
zeros
([
ndim
**
2
,
cell
.
size
()])
cell
.
evaluate_stress_tangent
(
µ
DE
);
µ
DE
[:,:]
=
DE
[:,:,
0
,
0
]
.
reshape
([
-
1
,
1
])
# initial residual: distribute "DE" over grid using "K4"
b
=
-
G_K_deps
(
DE
)
G_K_deps2
=
lambda
de
:
cell
.
directional_stiffness
(
de
.
reshape
(
µ
DE
.
shape
))
.
ravel
()
b2
=
-
G_K_deps2
(
µ
DE
)
.
ravel
()
print
(
"b2.shape = {}"
.
format
(
b2
.
shape
))
eps
+=
DE
En
=
np
.
linalg
.
norm
(
eps
)
iiter
=
0
# iterate as long as the iterative update does not vanish
class
accumul
(
object
):
def
__init__
(
self
):
self
.
counter
=
0
def
__call__
(
self
,
x
):
self
.
counter
+=
1
print
(
"at step {}: ||Ax-b||/||b|| = {}"
.
format
(
self
.
counter
,
np
.
linalg
.
norm
(
G_K_deps
(
x
)
-
b
)
/
np
.
linalg
.
norm
(
b
)))
acc
=
accumul
()
cg_tol
=
1e-8
maxiter
=
60
solver
=
µ
.
solvers
.
SolverCGEigen
(
cell
,
cg_tol
,
maxiter
,
verbose
=
True
)
solver
=
µ
.
solvers
.
SolverCG
(
cell
,
cg_tol
,
maxiter
,
verbose
=
True
)
try
:
r
=
µ
.
solvers
.
newton_cg
(
cell
,
delEps0
,
solver
,
1e-5
,
verbose
=
True
)
except
Exception
as
err
:
print
(
err
)
while
True
:
depsm
,
_
=
sp
.
cg
(
tol
=
cg_tol
,
A
=
sp
.
LinearOperator
(
shape
=
(
ndof
,
ndof
),
matvec
=
G_K_deps
,
dtype
=
'float'
),
b
=
b
,
callback
=
acc
)
# solve linear system using CG
#depsm2,_ = sp.cg(tol=1.e-8,
# A = sp.LinearOperator(shape=(ndof,ndof),matvec=G_K_deps2,dtype='float'),
# b = b2,
# callback=acc
#) # solve linear system using CG
eps
+=
depsm
.
reshape
(
ndim
,
ndim
,
N
,
N
)
# update DOFs (array -> tens.grid)
sig
=
ddot42
(
K4
,
eps
)
# new residual stress
b
=
-
G
(
sig
)
# convert residual stress to residual
print
(
'
%10.2e
'
%
(
np
.
max
(
depsm
)
/
En
))
# print residual to the screen
print
(
np
.
linalg
.
norm
(
depsm
)
/
np
.
linalg
.
norm
(
En
))
if
np
.
linalg
.
norm
(
depsm
)
/
En
<
1.e-5
and
iiter
>
0
:
break
# check convergence
iiter
+=
1
print
(
"nb_cg: {0}"
.
format
(
acc
.
counter
))
import
matplotlib.pyplot
as
plt
s11
=
sig
[
0
,
0
,
:,:]
s22
=
sig
[
1
,
1
,
:,:]
s21_2
=
sig
[
0
,
1
,
:,
:]
*
sig
[
1
,
0
,:,
:]
vonM1
=
np
.
sqrt
(
.
5
*
((
s11
-
s22
)
**
2
)
+
s11
**
2
+
s22
**
2
+
6
*
s21_2
)
#vonM1 = np.sqrt(sig[0, 0, :, :]**2 + sig[1, 1, :, :]**2 - sig[0, 0, :, :] * sig[1, 1, :, :] + 3 * sig[0,1]**2)
plt
.
figure
()
img
=
plt
.
pcolormesh
(
vonM1
)
#eps[0,1,:,:])
plt
.
title
(
"goose"
)
plt
.
colorbar
(
img
)
try
:
print
(
r
.
stress
.
shape
)
arr
=
r
.
stress
.
T
.
reshape
(
N
,
N
,
ndim
,
ndim
)
s11
=
arr
[:,:,
0
,
0
]
s22
=
arr
[:,:,
1
,
1
]
s21_2
=
arr
[:,:,
0
,
1
]
*
arr
[:,:,
1
,
0
]
vonM2
=
np
.
sqrt
(
.
5
*
((
s11
-
s22
)
**
2
)
+
s11
**
2
+
s22
**
2
+
6
*
s21_2
)
plt
.
figure
()
plt
.
title
(
"µSpectre"
)
img
=
plt
.
pcolormesh
(
vonM2
)
#eps[0,1,:,:])
plt
.
colorbar
(
img
)
print
(
"err = {}"
.
format
(
np
.
linalg
.
norm
(
vonM1
-
vonM2
)))
print
(
"err2 = {}"
.
format
(
np
.
linalg
.
norm
(
vonM1
-
vonM1
.
T
)))
print
(
"err3 = {}"
.
format
(
np
.
linalg
.
norm
(
vonM2
-
vonM2
.
T
)))
plt
.
figure
()
plt
.
title
(
"diff"
)
img
=
plt
.
pcolormesh
(
vonM1
-
vonM2
.
T
)
#eps[0,1,:,:])
plt
.
colorbar
(
img
)
except
Exception
as
err
:
print
(
err
)
plt
.
show
()
Event Timeline
Log In to Comment