Skip to content
This repository has been archived by the owner on Aug 8, 2024. It is now read-only.

Commit

Permalink
Merge pull request #6 from eisenforschung/add_doc
Browse files Browse the repository at this point in the history
Add doc
  • Loading branch information
samwaseda authored Aug 4, 2024
2 parents d887c57 + a256b44 commit 57531b8
Show file tree
Hide file tree
Showing 4 changed files with 304 additions and 161 deletions.
49 changes: 33 additions & 16 deletions elaston/linear_elasticity/eshelby.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,10 @@
import numpy as np

__author__ = "Sam Waseda"
__copyright__ = "Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH " \
"- Computational Materials Design (CM) Department"
__copyright__ = (
"Copyright 2021, Max-Planck-Institut für Eisenforschung GmbH "
"- Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Sam Waseda"
__email__ = "[email protected]"
Expand All @@ -21,6 +23,7 @@ class Eshelby:
All notations follow the original paper.
"""

def __init__(self, elastic_tensor, burgers_vector):
self.elastic_tensor = elastic_tensor
self.burgers_vector = burgers_vector
Expand All @@ -33,17 +36,21 @@ def _get_pmat(self, x):
return (
self.elastic_tensor[:, 0, :, 0]
+ np.einsum(
'...,ij->...ij', x, self.elastic_tensor[:, 0, :, 1]+self.elastic_tensor[:, 1, :, 0]
"...,ij->...ij",
x,
self.elastic_tensor[:, 0, :, 1] + self.elastic_tensor[:, 1, :, 0],
)
+ np.einsum('...,ij->...ij', x**2, self.elastic_tensor[:, 1, :, 1])
+ np.einsum("...,ij->...ij", x**2, self.elastic_tensor[:, 1, :, 1])
)

@property
def p(self):
if self._p is None:
coeff = np.polyfit(self.fit_range, np.linalg.det(self._get_pmat(self.fit_range)), 6)
coeff = np.polyfit(
self.fit_range, np.linalg.det(self._get_pmat(self.fit_range)), 6
)
self._p = np.roots(coeff)
self._p = self._p[np.imag(self._p)>0]
self._p = self._p[np.imag(self._p) > 0]
return self._p

@property
Expand All @@ -59,13 +66,15 @@ def Ak(self):
@property
def D(self):
if self._D is None:
F = np.einsum('n,ij->nij', self.p, self.elastic_tensor[:, 1, :, 1])
F = np.einsum("n,ij->nij", self.p, self.elastic_tensor[:, 1, :, 1])
F += self.elastic_tensor[:, 1, :, 0]
F = np.einsum('nik,nk->ni', F, self.Ak)
F = np.einsum("nik,nk->ni", F, self.Ak)
F = np.concatenate((F.T, self.Ak.T), axis=0)
F = np.concatenate((np.real(F), -np.imag(F)), axis=-1)
self._D = np.linalg.solve(F, np.concatenate((np.zeros(3), self.burgers_vector)))
self._D = self._D[:3]+1j*self._D[3:]
self._D = np.linalg.solve(
F, np.concatenate((np.zeros(3), self.burgers_vector))
)
self._D = self._D[:3] + 1j * self._D[3:]
return self._D

@property
Expand All @@ -74,7 +83,7 @@ def dzdx(self):

def _get_z(self, positions):
z = np.stack((np.ones_like(self.p), self.p), axis=-1)
return np.einsum('nk,...k->...n', z, np.asarray(positions)[..., :2])
return np.einsum("nk,...k->...n", z, np.asarray(positions)[..., :2])

def get_displacement(self, positions):
"""
Expand All @@ -87,8 +96,10 @@ def get_displacement(self, positions):
((n,3)-array): Displacement vectors
"""
return np.imag(
np.einsum('nk,n,...n->...k', self.Ak, self.D, np.log(self._get_z(positions)))
)/(2*np.pi)
np.einsum(
"nk,n,...n->...k", self.Ak, self.D, np.log(self._get_z(positions))
)
) / (2 * np.pi)

def get_strain(self, positions):
"""
Expand All @@ -101,7 +112,13 @@ def get_strain(self, positions):
((n,3,3)-array): Strain tensors
"""
strain = np.imag(
np.einsum('ni,n,...n,nj->...ij', self.Ak, self.D, 1/self._get_z(positions), self.dzdx)
np.einsum(
"ni,n,...n,nj->...ij",
self.Ak,
self.D,
1 / self._get_z(positions),
self.dzdx,
)
)
strain = strain+np.einsum('...ij->...ji', strain)
return strain/4/np.pi
strain = strain + np.einsum("...ij->...ji", strain)
return strain / 4 / np.pi
Loading

0 comments on commit 57531b8

Please sign in to comment.