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 8cda37b
Show file tree
Hide file tree
Showing 15 changed files with 508 additions and 348 deletions.
99 changes: 55 additions & 44 deletions examples/cylinder.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
"""
import os
import sys

sys.path.append("../src/")
sys.path.append("../geo/")
sys.path.append("tests/msh/")
Expand All @@ -55,7 +56,7 @@


def mkdir_p(path):
""" make directory at 'path' (full or relative) if not existing """
"""make directory at 'path' (full or relative) if not existing"""
try:
os.makedirs(path)
except OSError as exc: # Python >2.5
Expand Down Expand Up @@ -117,11 +118,11 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
raise Exception("U[bcDofs] != Uref[bcDofs]")


print('============================================')
print('= =')
print('= 2D plane strain thick cylinder =')
print('= =')
print('============================================')
print("============================================")
print("= =")
print("= 2D plane strain thick cylinder =")
print("= =")
print("============================================")


VERBOSE = True
Expand All @@ -140,12 +141,12 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):

if material == "steel":
# steel
matmodel = 'StVenant'
matmodel = "StVenant"
Eyoung = 210e5 # Young's modulus (Pa)
nu = 0.3 # Poisson's ratio
else:
# filled rubber (very stiff):
matmodel = 'NeoHookean'
matmodel = "NeoHookean"
Eyoung = 10e6 # Young's modulus (Pa)
nu = 0.45 # Poisson's ratio

Expand All @@ -164,7 +165,7 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# ITERATION PARAMETERS #
#################################
nSteps = 10 # number of steps in the incremental loading loop
tStep = 1. / nSteps # size of the steps
tStep = 1.0 / nSteps # size of the steps
# rdisp = 0.00001 # constant radial displacement on interior boundary (StVenant)
rdisp = 0.0001 # constant radial displacement on interior boundary (StVenant)
# rdisp = 0.01 #constant radial displacement on interior boundary (NeoHookean)
Expand Down Expand Up @@ -209,7 +210,7 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# basemesh = "cylinder-quad_coarse.msh"
# basemesh = "cylinder-quad_refined.msh"

basegeo = 'cylinder.geo'
basegeo = "cylinder.geo"
if VERBOSE:
print("############################")
print("# FILES AND FOLDERS #")
Expand All @@ -227,10 +228,12 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):

with open(geo_folder + basegeo) as geomfile:
lines = geomfile.readlines()
Ri = float(lines[0].strip().strip(";").split(
" = ")[1]) # interior radius (m)
Re = float(lines[1].strip().strip(";").split(
" = ")[1]) # exterior radius (m)
Ri = float(
lines[0].strip().strip(";").split(" = ")[1]
) # interior radius (m)
Re = float(
lines[1].strip().strip(";").split(" = ")[1]
) # exterior radius (m)

if VERBOSE:
print(" Geometry cylinder:")
Expand All @@ -239,7 +242,7 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):


inputfile = msh_folder + basemesh
meshfile = open(inputfile, 'r')
meshfile = open(inputfile, "r")
mesh = gmsh.gmshInput_mesh(meshfile, d=2)
meshfile.close()
dim = mesh.dim # dimension of the problem
Expand All @@ -257,17 +260,17 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# GSMH OUTPUT #
#################################

mesh_name = basemesh.split('-')[0]
mesh_name = basemesh.split("-")[0]
# 'tri_coarse' = element type + mesh refinement
mesh_type = basemesh.split('-')[1][:-4]
mesh_type = basemesh.split("-")[1][:-4]
# subdirectory for output files
new_folder_name = mesh_type

dirout = os.path.join(out_folder + matmodel, new_folder_name)
mkdir_p(dirout) # create directory for output files if not already existing
outputfile = mesh_name + '-results_Nsteps-' + str(int(nSteps)) + ".msh"
outputfile = mesh_name + "-results_Nsteps-" + str(int(nSteps)) + ".msh"
outputfile = os.path.join(dirout, outputfile)
gmsh_out = open(outputfile, 'w')
gmsh_out = open(outputfile, "w")
gmsh.gmshOutput_mesh(gmsh_out, mesh)

plot.plot_mesh(mesh, label="Before", color="b")
Expand All @@ -277,7 +280,7 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# PREPROSSESING #
#################################

#--- boundary conditions
# --- boundary conditions
#
# On the left side of the quarter-cylinder (X=0), the displacements in the x-direction are fixed (u_x =0)
# At the bottom of the quarter-cylinder (Y=0), the displacements in the y-direction are fixed (u_y=0)
Expand Down Expand Up @@ -352,9 +355,9 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
freeDofs = list(np.sort(freeDofs))
# display fixed degrees of freedom
if VERBOSE:
print('############################')
print('# PRE PROCESSING #')
print('############################')
print("############################")
print("# PRE PROCESSING #")
print("############################")
print(" Mesh information:")
print(" Total nodes: %d" % nNodes)
print(" Total elements: %d" % nElements)
Expand All @@ -380,17 +383,17 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# stop()

#
#--- create material model instance
# --- create material model instance
#
if matmodel == "StVenant":
material = elasticity.StVenantKirchhoffElasticity(Eyoung, nu)
elif matmodel == "NeoHookean":
material = elasticity.NeoHookeanElasticity(Eyoung, nu)
else:
raise(KeyError, "matmodel must be one of 'StVenant', 'NeoHookean'")
raise (KeyError, "matmodel must be one of 'StVenant', 'NeoHookean'")

#
#--- create FE model instance
# --- create FE model instance
#
model = fem.FEModel(mesh, material)

Expand Down Expand Up @@ -426,12 +429,12 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
# SIMULATION #
#################################
if VERBOSE:
print('############################')
print('# SIMULATION #')
print('############################')
print("############################")
print("# SIMULATION #")
print("############################")
start_time = datetime.now() # measure time for the simulation
time = 0.0 # final time is nSteps x tStep = nSteps x (1/nSteps) = 1
itermoy = 0. # average number of iterations per time step
itermoy = 0.0 # average number of iterations per time step
# prescribed displacement per step on the right hand side of the membrane
displacements = []
reactions = [] # corresponding reaction per step
Expand Down Expand Up @@ -524,14 +527,16 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):

norm = np.linalg.norm(R)
test = TOLERANCE * norm
if (test < PRECISION):
if test < PRECISION:
test = PRECISION
iteration = 0
if VERBOSE:
print(' *** Iteration %02d: residual norm =%.3e ***' %
(iteration, norm))
print(
" *** Iteration %02d: residual norm =%.3e ***"
% (iteration, norm)
)

while (norm > test): # Newton-Raphson iterations
while norm > test: # Newton-Raphson iterations

# compute correction
# (U1, U2) = cut(U, bcDofs, freeDofs)
Expand Down Expand Up @@ -562,8 +567,10 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
norm = np.linalg.norm(R2)
iteration = iteration + 1
if VERBOSE:
print(' *** Iteration %02d: residual norm =%.3e ***' %
(iteration, norm))
print(
" *** Iteration %02d: residual norm =%.3e ***"
% (iteration, norm)
)

# =====================================================================
if iteration == ITERMAX:
Expand Down Expand Up @@ -591,7 +598,7 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
reactions.append(nodalForce) # (N)
displacements.append(uint) # (m)
#
#--- get strain and stress tensor fields for current time step
# --- get strain and stress tensor fields for current time step
#
E = []
P = []
Expand All @@ -605,17 +612,21 @@ def verify(U, Uref, bcDofs, TOLERANCE=1e-9):
euler.append(model.getStrainEulerAlmansi(n).flatten())
sigma.append(model.getStressCauchy(n).flatten())
#
#--- gmsh output for current time step
# --- gmsh output for current time step
#
gmsh.gmshOutput_nodal(gmsh_out, "Residual",
R.reshape(nNodes, dim), iStep, time)
gmsh.gmshOutput_nodal(gmsh_out, "Displacements",
U.reshape((nNodes, dim)), iStep, time)
gmsh.gmshOutput_nodal(
gmsh_out, "Residual", R.reshape(nNodes, dim), iStep, time
)
gmsh.gmshOutput_nodal(
gmsh_out, "Displacements", U.reshape((nNodes, dim)), iStep, time
)
gmsh.gmshOutput_element(gmsh_out, "Green-Lagrange strain", E, iStep, time)
gmsh.gmshOutput_element(
gmsh_out, "Piola Kirchhoff I stress", P, iStep, time)
gmsh_out, "Piola Kirchhoff I stress", P, iStep, time
)
gmsh.gmshOutput_element(
gmsh_out, "Euler-Almansi strain", euler, iStep, time)
gmsh_out, "Euler-Almansi strain", euler, iStep, time
)
gmsh.gmshOutput_element(gmsh_out, "Hencky strain", hencky, iStep, time)
gmsh.gmshOutput_element(gmsh_out, "Cauchy stress", sigma, iStep, time)

Expand Down
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
Loading

0 comments on commit 8cda37b

Please sign in to comment.