Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/test/single_scale/PARADIS/test_barnett/barnett.py
908e42023cdemaster

Last Change
Raw File

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 imp
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)
foo = imp.load_source('myscript', fname)
def computeDisloDisplacements(A, x0):
return foo.compute(A=A, mdpos0=x0)
return computeDisloDisplacements
c4science · Help