Skip to content

Commit

Permalink
style: apply isort and black formaters
Browse files Browse the repository at this point in the history
  • Loading branch information
@carlos-adir
carlos-adir committed Feb 14, 2025
1 parent d5fcb62 commit 8a6fe0c
Show file tree
Hide file tree
Showing 14 changed files with 453 additions and 304 deletions.
15 changes: 10 additions & 5 deletions src/hyper/elasticity.py
View file
Original file line number Diff line number Diff line change
@@ -1,20 +1,23 @@
# -*- coding: utf-8 -*-

from abc import ABC, abstractmethod # For abstract classes
from . import tensor

import numpy as np

from . import tensor


class Elasticity(ABC):
"""
Data structure for (isotropic) hyperelaticity models
1ST_lamdaE_CONSTANT is the frist lamé parameter lambda
2ST_lamdaE_CONSTANT is the second lamé parameter mu
"""

prop = dict()

def __init__(self, E, nu):
G = E / (2 * (1. + nu))
G = E / (2 * (1.0 + nu))
K = E / (3 * (1 - 2 * nu))
lamda = K - (2 / 3) * G
mu = G
Expand Down Expand Up @@ -88,7 +91,9 @@ def potential(self, C):
lamda = self.get1LAME()
mu = self.get2LAME()
EL = 0.5 * (C - tensor.I(len(C))) # Lagrangian strain E
phi = lamda / 2. * (tensor.trace(EL))**2 + mu * np.tensordot(EL, EL, 2)
phi = lamda / 2.0 * (tensor.trace(EL)) ** 2 + mu * np.tensordot(
EL, EL, 2
)
return phi

def stress(self, C):
Expand Down Expand Up @@ -139,9 +144,9 @@ def potential(self, C):
C = np.zeros((3, 3))
C[:2, :2] = K[:, :]
J = np.sqrt(tensor.det(C)) # J = det(F) and det(C) = J^2
part1 = (mu / 2) * (tensor.trace(C) - 3.)
part1 = (mu / 2) * (tensor.trace(C) - 3.0)
part2 = mu * np.log(J)
part3 = (lamda / 2) * (np.log(J))**2
part3 = (lamda / 2) * (np.log(J)) ** 2
phi = part1 - part2 + part3
return phi

Expand Down
29 changes: 15 additions & 14 deletions src/hyper/fem.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,15 +5,15 @@


import numpy as np
from . import tensor
from . import geometry
from scipy.linalg import polar

from . import geometry, tensor

########################################################################
# Definition of elements for Finite Elements Method #
########################################################################


class FiniteElement:
"""
Data structure for finite element
Expand All @@ -26,13 +26,13 @@ def __init__(self, t, xNod):
#
self.calculate_shapes(t, xNod)

self.F = [] # deformation gradient F = Grad(u) + I
self.hencky = [] # hencky strain ln(V) with F = V.R
self.E_GL = [] # lagrangian strain E = 1/2 (C - I)
self.E_EA = [] # E_EA strain e = 1/2 (I - b^-1)
self.PK1 = [] # piola kirchoff I : P = F*S
self.sigma = [] # cauchy stress : \sigma
self.K = [] # lagrangian tangent operator dP/dF
self.F = [] # deformation gradient F = Grad(u) + I
self.hencky = [] # hencky strain ln(V) with F = V.R
self.E_GL = [] # lagrangian strain E = 1/2 (C - I)
self.E_EA = [] # E_EA strain e = 1/2 (I - b^-1)
self.PK1 = [] # piola kirchoff I : P = F*S
self.sigma = [] # cauchy stress : \sigma
self.K = [] # lagrangian tangent operator dP/dF
d = self.getDim() # space dimension
npg = self.getNIntPts()
for n in range(npg):
Expand Down Expand Up @@ -133,14 +133,14 @@ def update(self, uNod, mater):
self.E_GL[i] = 0.5 * (C - tensor.I(len(C)))
# compute spatial description
# from scipy.linalg.polar() method with u=R and p=V,
_, V = polar(F, side='left')
_, V = polar(F, side="left")
# replace pure zeros by very low values to prevent "nan" in np.log(V)
V[V < 1e-10] = 1e-15
self.hencky[i] = np.log(V) # ln(V) with F = V.R, "true" strain
b = tensor.leftCauchyGreen(F)
self.E_EA[i] = 0.5 * (tensor.I(len(b)) - tensor.inv(b))
###
if (mater == 0): # skip next lines: do not compute stress
if mater == 0: # skip next lines: do not compute stress
# print("Whahat")
continue
# compute PK2 stress tensor and material tangent operator M=2*dS/dC
Expand Down Expand Up @@ -224,7 +224,7 @@ def __init__(self, theMesh, m=0):
nodei = theMesh.getNode(label_nodei)
xnodei = nodei.getX()
# Since we don't need z, we cut it out
xnodei = xnodei[:self.dim]
xnodei = xnodei[: self.dim]
xNod.append(xnodei)
xNod = np.array(xNod)
new_FiniteElement = FiniteElement(elem.type, xNod)
Expand Down Expand Up @@ -280,7 +280,8 @@ def extract(self, U, e):
"""
if not isinstance(e, (int)):
raise Exception(
"The element 'e' in 'extract' function is not a integer")
"The element 'e' in 'extract' function is not a integer"
)
# print("index = " + str(index))
connect_e = self.connect[e]
# print("connect_e = " + str(connect_e))
Expand Down Expand Up @@ -409,7 +410,7 @@ def getStressPK1(self, n):
return avg

def getStrainEulerAlmansi(self, n):
""" Return average of EA strain at all integration points of element n"""
"""Return average of EA strain at all integration points of element n"""
avg = tensor.tensor(self.dim)
for tens in self.elems[n].E_EA:
avg += tens
Expand Down
80 changes: 57 additions & 23 deletions src/hyper/geometry.py
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,5 @@
import numpy as np


########################################################################
# Gauss Quadrature #
########################################################################
Expand Down Expand Up @@ -33,6 +32,7 @@
# | | (-b, -b), (0, -b), (b, -b) | 25/81, 40/81, 25/81 |
# ‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾


class IPTri:
"""
Numerical integration rules for triangles
Expand All @@ -41,23 +41,28 @@ class IPTri:
P2 = (1, 0)
P3 = (0, 1)
"""

@staticmethod
def getX(npg=3):
if npg == 1:
X = [(1. / 3., 1 / 3.)]
X = [(1.0 / 3.0, 1 / 3.0)]
elif npg == 3:
X = [(2 / 3, 1 / 6), (1 / 6, 1 / 6), (1 / 6, 2 / 3)]
elif npg == 4:
X = [(1 / 3, 1 / 3), (3 / 5, 1 / 5),
(1 / 5, 1 / 5), (1 / 5, 3 / 5)]
X = [
(1 / 3, 1 / 3),
(3 / 5, 1 / 5),
(1 / 5, 1 / 5),
(1 / 5, 3 / 5),
]
else:
raise ValueError("Not found npg = " + str(npg))
return X

@staticmethod
def getW(npg=3):
if npg == 1:
W = [1. / 2.]
W = [1.0 / 2.0]
elif npg == 3:
W = [1 / 6, 1 / 6, 1 / 6]
elif npg == 4:
Expand Down Expand Up @@ -86,9 +91,17 @@ def getX(npg=3):
X = [(-a, a), (-a, -a), (a, -a), (a, a)]
elif npg == 9:
b = np.sqrt(3 / 5)
X = [(-b, b), (0, b), (b, b),
(-b, 0), (0, 0), (b, 0),
(-b, -b), (0, -b), (b, -b)]
X = [
(-b, b),
(0, b),
(b, b),
(-b, 0),
(0, 0),
(b, 0),
(-b, -b),
(0, -b),
(b, -b),
]
else:
raise ValueError("Not found npg = " + str(npg))
return X
Expand All @@ -100,9 +113,17 @@ def getW(npg=4):
elif npg == 4:
W = [1, 1, 1, 1]
elif npg == 9:
W = [25 / 81, 40 / 81, 25 / 81,
40 / 81, 64 / 81, 40 / 81,
25 / 81, 40 / 81, 25 / 81]
W = [
25 / 81,
40 / 81,
25 / 81,
40 / 81,
64 / 81,
40 / 81,
25 / 81,
40 / 81,
25 / 81,
]
else:
raise ValueError("Not found npg = " + str(npg))
return W
Expand Down Expand Up @@ -131,6 +152,7 @@ def getW(npg=4):
# (-1, -1)‾‾‾‾‾‾‾‾‾(1, -1)
#


class SFT3:
"""
Shape functions for 3-node triangle
Expand All @@ -139,6 +161,7 @@ class SFT3:
P2 = (1, 0)
P3 = (0, 1)
"""

@staticmethod
def N(x):
N1 = 1 - x[0] - x[1]
Expand Down Expand Up @@ -178,6 +201,7 @@ class SFT6:
P5 = (0, 1)
P6 = (0, 0.5)
"""

@staticmethod
def N(x):
# Ni(Pj) = delta_ij
Expand Down Expand Up @@ -219,6 +243,7 @@ class SFQ4:
P3 = (1, 1)
P4 = (-1, 1)
"""

@staticmethod
def N(x):
N1 = (x[0] - 1) * (x[1] - 1) / 4
Expand Down Expand Up @@ -257,6 +282,7 @@ class SFQ8:
P7 = (-1, 1)
P8 = (-1, 0)
"""

@staticmethod
def N(x):
N1 = -(x[0] - 1) * (x[1] - 1) * (x[0] + x[1] + 1) / 4
Expand All @@ -271,18 +297,26 @@ def N(x):

@staticmethod
def dN(x):
dN1 = [(2 * x[0] + x[1]) * (1 - x[1]) / 4,
(1 - x[0]) * (x[0] + 2 * x[1]) / 4]
dN2 = [x[0] * (x[1] - 1), (x[0]**2 - 1) / 2]
dN3 = [(-2 * x[0] + x[1]) * (x[1] - 1) / 4,
(-x[0] + 2 * x[1]) * (x[0] + 1) / 4]
dN4 = [(1 - x[1]**2) / 2, -x[1] * (x[0] + 1)]
dN5 = [(2 * x[0] + x[1]) * (x[1] + 1) / 4,
(x[0] + 1) * (x[0] + 2 * x[1]) / 4]
dN6 = [-x[0] * (x[1] + 1), (1 - x[0]**2) / 2]
dN7 = [(2 * x[0] - x[1]) * (x[1] + 1) / 4,
(x[0] - 1) * (x[0] - 2 * x[1]) / 4]
dN8 = [(x[1]**2 - 1) / 2, x[1] * (x[0] - 1)]
dN1 = [
(2 * x[0] + x[1]) * (1 - x[1]) / 4,
(1 - x[0]) * (x[0] + 2 * x[1]) / 4,
]
dN2 = [x[0] * (x[1] - 1), (x[0] ** 2 - 1) / 2]
dN3 = [
(-2 * x[0] + x[1]) * (x[1] - 1) / 4,
(-x[0] + 2 * x[1]) * (x[0] + 1) / 4,
]
dN4 = [(1 - x[1] ** 2) / 2, -x[1] * (x[0] + 1)]
dN5 = [
(2 * x[0] + x[1]) * (x[1] + 1) / 4,
(x[0] + 1) * (x[0] + 2 * x[1]) / 4,
]
dN6 = [-x[0] * (x[1] + 1), (1 - x[0] ** 2) / 2]
dN7 = [
(2 * x[0] - x[1]) * (x[1] + 1) / 4,
(x[0] - 1) * (x[0] - 2 * x[1]) / 4,
]
dN8 = [(x[1] ** 2 - 1) / 2, x[1] * (x[0] - 1)]
return [dN1, dN2, dN3, dN4, dN5, dN6, dN7, dN8]

@staticmethod
Expand Down
Loading

0 comments on commit 8a6fe0c

Please sign in to comment.