Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F88854649
barnett.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
Mon, Oct 21, 00:34
Size
10 KB
Mime Type
text/x-python
Expires
Wed, Oct 23, 00:34 (2 d)
Engine
blob
Format
Raw Data
Handle
21829558
Attached To
rLIBMULTISCALE LibMultiScale
barnett.py
View Options
#!/usr/bin/env python
"""
Module in charge of
- providing the tools to compute dislocation
displacements due to the presence of a set of closed loops
from the Barnet/Fivel method.
"""
################################################################
import
os
import
sys
import
importlib.util
import
numpy
as
np
################################################################
def
computeClosurePoint
(
A
,
B
,
b
,
slip_plane
=
None
):
# Marc Fivel & Christophe Depres,
# Philosophical Magazine, Volume 94, Issue 28, 2014
# The closure point is defined as the intersection point
# between arbitrary vector of u to the current slip plane.
# The slip plane can be defined by cross product of burgers vector
# and vector of AB
O
=
np
.
array
([
0.0
,
0.0
,
0.0
])
u
=
np
.
array
([
1.0
,
1.0
,
1.0
])
# arbitrary vector from the origin point O
# http://mathworld.wolfram.com/Line-PlaneIntersection.html
slip_normal
=
np
.
cross
(
B
-
A
,
b
)
# print slip_normal
# print np.linalg.norm(slip_normal)
if
np
.
linalg
.
norm
(
slip_normal
)
<
1e-15
:
if
slip_plane
is
None
:
raise
Exception
(
'slip normal must be specified'
)
slip_direction
=
np
.
cross
(
b
,
slip_plane
)
B
=
A
+
slip_direction
X1
=
A
+
b
X2
=
A
.
copy
()
X3
=
B
.
copy
()
X4
=
O
.
copy
()
X5
=
u
.
copy
()
nom_mat
=
np
.
zeros
((
4
,
4
))
nom_mat
[
0
,
:]
=
1
nom_mat
[
1
:,
0
]
=
X1
nom_mat
[
1
:,
1
]
=
X2
nom_mat
[
1
:,
2
]
=
X3
nom_mat
[
1
:,
3
]
=
X4
nom
=
np
.
linalg
.
det
(
nom_mat
)
denom_mat
=
np
.
zeros
((
4
,
4
))
denom_mat
[
0
,
:
3
]
=
1
denom_mat
[
1
:,
0
]
=
X1
denom_mat
[
1
:,
1
]
=
X2
denom_mat
[
1
:,
2
]
=
X3
denom_mat
[
1
:,
3
]
=
X5
-
X4
denom
=
np
.
linalg
.
det
(
denom_mat
)
# print denom_mat
# print denom
t
=
-
nom
/
denom
C
=
X4
+
(
X5
-
X4
)
*
t
return
C
################################################################
def
SolidAngleSegmentBarnett
(
vec_la
,
vec_lb
,
vec_lc
,
N
):
# print("vec_la:", vec_la)
# print("vec_lb:", vec_lb)
# print("vec_lc:", vec_lc)
val_a
=
np
.
arccos
(
np
.
einsum
(
'ai,ai->a'
,
vec_lb
,
vec_lc
))
val_b
=
np
.
arccos
(
np
.
einsum
(
'ai,ai->a'
,
vec_la
,
vec_lc
))
val_c
=
np
.
arccos
(
np
.
einsum
(
'ai,ai->a'
,
vec_la
,
vec_lb
))
val_s
=
.
5
*
(
val_a
+
val_b
+
val_c
)
# print("val_a:", val_a)
# print("val_b:", val_b)
# print("val_c:", val_c)
# print("val_s:", val_s)
tan2e_4
=
(
np
.
tan
(
0.5
*
val_s
)
*
np
.
tan
(
0.5
*
(
val_s
-
val_a
))
*
np
.
tan
(
0.5
*
(
val_s
-
val_b
))
*
np
.
tan
(
0.5
*
(
val_s
-
val_c
)))
# print("tan2e_4:", tan2e_4)
# print("std::sqrt(tan2e_4):", np.sqrt(tan2e_4))
val_e
=
4.0
*
np
.
arctan
(
np
.
sqrt
(
tan2e_4
))
val_e
[
val_e
>
2.0
*
np
.
pi
]
=
0.
val_e
[
val_e
<
0.0
]
=
0.
# print("val_e:", val_e)
# print("la dot N:", vec_la.dot(N))
res
=
-
np
.
sign
(
vec_la
.
dot
(
N
))
*
val_e
# print 'omega_abc:', res
return
res
################################################################
def
computeProjectedPoint
(
A
,
C
,
burg
):
# Marc Fivel & Christophe Depres,
# Philosophical Magazine, Volume 94, Issue 28, 2014
X1
=
A
+
burg
*
1e+8
X2
=
A
-
burg
*
1e+8
X0
=
C
AC
=
X0
-
X1
AB
=
X2
-
X1
n
=
AB
/
np
.
linalg
.
norm
(
AB
)
dotACn
=
AC
.
dot
(
n
)
P
=
X1
+
dotACn
*
n
return
P
################################################################
def
segmentDisplacementsOfAtom
(
A
,
B
,
C
,
X
,
burg
,
pois
,
**
kwargs
):
""" computes the displacement field at point X of
a triangular loop made of 3 segments AB, BC, CA
burg: the burgers vector
pois: the poisson ratio
"""
# print("A: ", A)
# print("B: ", B)
# print("C: ", C)
vec_ab
=
B
-
A
vec_bc
=
C
-
B
vec_ca
=
C
-
A
# print("AB: ", vec_ab)
# print("BC: ", vec_bc)
# print("CA: ", vec_ca)
ab_norm
=
np
.
linalg
.
norm
(
vec_ab
)
bc_norm
=
np
.
linalg
.
norm
(
vec_bc
)
ca_norm
=
np
.
linalg
.
norm
(
vec_ca
)
# Zero length segment; nothing to do.
if
ab_norm
<
1e-20
or
bc_norm
<
1e-20
or
ca_norm
<
1e-20
:
return
np
.
zeros
(
X
.
shape
)
# compute t_ab,t_bc,t_ca
t_ab
=
vec_ab
/
ab_norm
t_bc
=
vec_bc
/
bc_norm
t_ca
=
vec_ca
/
ca_norm
# compute R_A,R_B,R_C
vec_ra
=
A
-
X
vec_rb
=
B
-
X
vec_rc
=
C
-
X
# print("vec_ra:", vec_ra)
# print("vec_rb:", vec_rb)
# print("vec_rc:", vec_rc)
# compute lambda_A,lambda_B,lambda_C
ra_norm
=
np
.
linalg
.
norm
(
vec_ra
,
axis
=
1
)
rb_norm
=
np
.
linalg
.
norm
(
vec_rb
,
axis
=
1
)
rc_norm
=
np
.
linalg
.
norm
(
vec_rc
,
axis
=
1
)
# print("ra_norm:", ra_norm)
# print("rb_norm:", rb_norm)
# print("rc_norm:", rc_norm)
vec_la
=
np
.
einsum
(
'ai->ia'
,
vec_ra
)
/
ra_norm
vec_lb
=
np
.
einsum
(
'ai->ia'
,
vec_rb
)
/
rb_norm
vec_lc
=
np
.
einsum
(
'ai->ia'
,
vec_rc
)
/
rc_norm
vec_la
=
np
.
einsum
(
'ia->ai'
,
vec_la
)
vec_lb
=
np
.
einsum
(
'ia->ai'
,
vec_lb
)
vec_lc
=
np
.
einsum
(
'ia->ai'
,
vec_lc
)
# print("vec_la:", vec_la)
# print("vec_lb:", vec_lb)
# print("vec_lc:", vec_lc)
# Calculate the solid angle.
# There will be no contribution to solid angle if AB and BC are colinear.
N
=
np
.
cross
(
vec_ab
,
vec_bc
)
n_norm
=
np
.
linalg
.
norm
(
N
)
if
n_norm
<
1e-20
:
omega_abc
=
np
.
zeros
(
X
.
shape
[
0
])
else
:
N
/=
n_norm
omega_abc
=
SolidAngleSegmentBarnett
(
vec_la
,
vec_lb
,
vec_lc
,
N
)
# print("omega_abc:", omega_abc)
# Calculate Barnett
# fAB = (b x tAB) ln[(Rb/Ra) . ((1 + lB . tAB)/(1 + lA . tAB))]
c_ab
=
np
.
log
(
rb_norm
/
ra_norm
*
(
(
1
+
vec_lb
.
dot
(
t_ab
))
/
(
1
+
vec_la
.
dot
(
t_ab
))))
c_bc
=
np
.
log
(
rc_norm
/
rb_norm
*
(
(
1
+
vec_lc
.
dot
(
t_bc
))
/
(
1
+
vec_lb
.
dot
(
t_bc
))))
c_ca
=
np
.
log
(
ra_norm
/
rc_norm
*
(
(
1
+
vec_la
.
dot
(
t_ca
))
/
(
1
+
vec_lc
.
dot
(
t_ca
))))
f_ab
=
np
.
einsum
(
'i,a->ai'
,
np
.
cross
(
burg
,
t_ab
),
c_ab
)
f_bc
=
np
.
einsum
(
'i,a->ai'
,
np
.
cross
(
burg
,
t_bc
),
c_bc
)
f_ca
=
np
.
einsum
(
'i,a->ai'
,
np
.
cross
(
burg
,
t_ca
),
c_ca
)
# Calculate Barnett gAB = [b . (lA x lB)](lA + lB) / (1 + lA . lB)
g_ab
=
np
.
einsum
(
'i,ai->a'
,
burg
,
np
.
cross
(
vec_la
,
vec_lb
))
/
(
1
+
np
.
einsum
(
'ai,ai->a'
,
vec_la
,
vec_lb
))
g_ab
=
np
.
einsum
(
'a,ai->ai'
,
g_ab
,
(
vec_la
+
vec_lb
))
g_bc
=
np
.
einsum
(
'i,ai->a'
,
burg
,
np
.
cross
(
vec_lb
,
vec_lc
))
/
(
1
+
np
.
einsum
(
'ai,ai->a'
,
vec_lb
,
vec_lc
))
g_bc
=
np
.
einsum
(
'a,ai->ai'
,
g_bc
,
(
vec_lb
+
vec_lc
))
g_ca
=
np
.
einsum
(
'i,ai->a'
,
burg
,
np
.
cross
(
vec_lc
,
vec_la
))
/
(
1
+
np
.
einsum
(
'ai,ai->a'
,
vec_lc
,
vec_la
))
g_ca
=
np
.
einsum
(
'a,ai->ai'
,
g_ca
,
(
vec_lc
+
vec_la
))
# Update the displacements.
coeff1
=
omega_abc
/
(
4.0
*
np
.
pi
)
# print(burg)
coeff3
=
1.
/
(
8.
*
np
.
pi
*
(
1
-
pois
))
coeff2
=
(
1.0
-
2.0
*
pois
)
*
coeff3
# U = - c2*burg - c2*(fAB+fBC+fCA) + c3*(gAB+gBC+gCA)
# print("coeff1 = ", coeff1)
# print("coeff2 = ", coeff2)
# print("coeff3 = ", coeff3)
U
=
(
-
np
.
einsum
(
'a,i->ai'
,
coeff1
,
burg
)
-
coeff2
*
(
f_ab
+
f_bc
+
f_ca
)
+
coeff3
*
(
g_ab
+
g_bc
+
g_ca
))
# print('U', U)
return
U
################################################################
def
computeAnalyticDisplacements
(
A
,
B
,
X
,
burg
,
pois
=.
3
,
**
kwargs
):
# print 'A:', A
# print 'B:', B
# Compute closure point for given AB segment
C
=
computeClosurePoint
(
A
,
B
,
burg
,
**
kwargs
)
# print("ClosurePoint=" + str(C) + " " + str(A) + " " + str(B))
# fivel's original method
# Compute point P1, projected point of C on line [A, A + b]
P1
=
computeProjectedPoint
(
A
,
C
,
burg
)
# print("P1=" + str(P1) + " " + str(A) + " " + str(B))
# Compute point P2, projected point of C on line [B, B + b]
P2
=
computeProjectedPoint
(
B
,
C
,
burg
)
# print("P2=" + str(P2) + " " + str(A) + " " + str(B))
#
# Compute displacement from the 2 triangles identified
U
=
segmentDisplacementsOfAtom
(
P2
,
A
,
B
,
X
,
burg
,
pois
)
# print U
U
+=
segmentDisplacementsOfAtom
(
P1
,
A
,
P2
,
X
,
burg
,
pois
)
# barnett's original method
# U = segmentDisplacementsOfAtom(C, A, B, X, burg, pois)
# print 'U:', U[:, 0]
# print("For X=" + str(X) + " U=" + str(U))
return
U
################################################################
def
computeDisloDisplacements
(
nodes
,
edges
,
X
,
burg
,
**
kwargs
):
disp
=
np
.
zeros
(
X
.
shape
)
for
e
in
edges
:
A
=
nodes
[
e
[
0
],
:]
B
=
nodes
[
e
[
1
],
:]
disp
+=
computeAnalyticDisplacements
(
A
,
B
,
X
,
burg
,
**
kwargs
)
.
copy
()
return
disp
################################################################
def
createDDCircularLoop
(
radius
,
npoints
):
nodes
=
np
.
zeros
((
npoints
,
3
))
nodes
[:,
0
]
=
radius
*
np
.
cos
(
np
.
linspace
(
0
,
2
*
np
.
pi
,
npoints
))
nodes
[:,
1
]
=
radius
*
np
.
sin
(
np
.
linspace
(
0
,
2
*
np
.
pi
,
npoints
))
edges
=
np
.
fromfunction
(
lambda
i
,
j
:
i
+
j
,
(
npoints
-
1
,
2
),
dtype
=
int
)
return
nodes
,
edges
################################################################
def
createDDSquareLoop
(
corners
,
npoints
):
npoints
=
np
.
array
(
npoints
)
corners
=
np
.
array
(
corners
)
nodes
=
list
()
for
i
in
[
0
,
1
,
2
,
3
]:
weights
=
np
.
linspace
(
0.
,
1.
,
npoints
[
i
])
# print weights
c1
=
corners
[
i
]
if
i
+
1
<=
3
:
c2
=
corners
[
i
+
1
]
else
:
c2
=
corners
[
0
]
nds
=
np
.
einsum
(
'k,i->ki'
,
(
1
-
weights
),
c1
)
nds
+=
np
.
einsum
(
'k,i->ki'
,
weights
,
c2
)
# print nds
nodes
+=
list
(
nds
)
nodes
=
np
.
array
(
nodes
)
# print nodes
edges
=
np
.
fromfunction
(
lambda
i
,
j
:
i
+
j
,
(
nodes
.
shape
[
0
],
2
),
dtype
=
int
)
edges
[
-
1
][
1
]
=
0
# print edges
return
nodes
,
edges
################################################################
# read paradis data
def
readParadisData
(
filename
):
" Read a paradis node file and build a numpy array of points "
_file
=
open
(
filename
)
flag
=
False
lines
=
_file
.
readlines
()
_file
.
close
()
for
i
,
line
in
enumerate
(
lines
):
line
=
line
.
strip
()
if
flag
is
True
:
break
if
line
==
'nodalData ='
:
flag
=
True
index
=
i
+
1
# print index
lines
=
lines
[
index
:]
lines
=
lines
[::
5
]
points
=
[]
for
line
in
lines
:
entries
=
line
.
split
()
entries
=
entries
[
1
:
4
]
coord
=
float
(
entries
[
0
]),
float
(
entries
[
1
]),
float
(
entries
[
2
])
points
.
append
(
coord
)
points
=
np
.
array
(
points
)
return
points
def
readParadisLoop
(
filename
):
nodes
=
readParadisData
(
filename
)
edges
=
np
.
fromfunction
(
lambda
i
,
j
:
i
+
j
,
(
nodes
.
shape
[
0
],
2
),
dtype
=
int
)
edges
[
-
1
][
1
]
=
0
return
nodes
,
edges
def
readScriptLoop
(
dirname
):
fname
=
os
.
path
.
join
(
dirname
,
'apply_barnett.py'
)
# print fname
sys
.
path
.
append
(
dirname
)
spec
=
importlib
.
util
.
spec_from_file_location
(
fname
,
"myscript"
)
foo
=
importlib
.
util
.
module_from_spec
(
spec
)
spec
.
loader
.
exec_module
(
foo
)
def
computeDisloDisplacements
(
A
,
x0
):
return
foo
.
compute
(
A
=
A
,
mdpos0
=
x0
)
return
computeDisloDisplacements
Event Timeline
Log In to Comment