Homec4science
Diffusion akantu (features/68-automatic-differentiation) Browse

rAKA/examples/python/automatic_differentiation/plate_hole.py
a1ea07d740e6features/68-automatic…

Last Change
Raw File

plate_hole.py
View Options

#!/usr/bin/env python
import subprocess
import akantu as aka
import numpy as np
# import the pyplot submodule to draw figures
import matplotlib.pyplot as plt
# import triangluation routine to plot meshes
import matplotlib.tri as tri
def solve_plate_hole(material_file, field='stress', F=10):
# geometric parameters
w = 10. # width (x-axis)
l = 5. # length (y-axis)
h1 = 0.05 # characteristic mesh size at the hole
h2 = 0.3 # characteristic mesh size in corners
R = 2. # radius of the hole
mesh_file = f"""
Point(1) = {{0, 0, 0, {h2} }};
Point(2) = {{ {R}, 0, 0, {h1} }};
Point(3) = {{ {w}, 0, 0, {h2} }};
Point(4) = {{ {w}, {l}, 0, {h2} }};
Point(5) = {{ 0, {l}, 0, {h2} }};
Point(6) = {{0, {R}, 0, {h1} }};
"""
mesh_file += """
Circle(1) = {6, 1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line Loop(6) = {1, 2, 3, 4, 5};
Plane Surface(7) = {6};
Physical Surface(8) = {7};
Physical Line("XBlocked") = {5};
Physical Line("YBlocked") = {2};
Physical Line("Traction") = {3};
"""
open('plate.geo', 'w').write(mesh_file)
ret = subprocess.run(
"gmsh -2 -order 1 -o plate.msh plate.geo", shell=True, check=True)
if ret.returncode:
print("Beware, gmsh could not run: mesh is not regenerated")
else:
print("Mesh generated")
spatial_dimension = 2
aka.parseInput(material_file)
mesh = aka.Mesh(spatial_dimension)
mesh_file = 'plate.msh'
mesh.read(mesh_file)
# parse input file
aka.parseInput(material_file)
model = aka.SolidMechanicsModel(mesh)
# initialize a static solver
model.initFull(_analysis_method=aka._static)
# configure the linear algebra solver
solver = model.getNonLinearSolver()
solver.set("max_iterations", 100)
solver.set("threshold", 1e-10)
solver.set("convergence_type", aka.SolveConvergenceCriteria.residual)
# set the displacement/Dirichlet boundary conditions
model.applyBC(aka.FixedValue(0.0, aka._x), "XBlocked")
model.applyBC(aka.FixedValue(0.0, aka._y), "YBlocked")
# set the force/Neumann boundary conditions
model.getExternalForce()[:] = 0
trac = [F, 0] # Newtons/m^2
model.applyBC(aka.FromTraction(trac), "Traction")
# compute the solution
try:
model.solveStep()
# n = solver.get('nb_iterations')
# print(f'Converged in {n} iterations')
except Exception as e:
print(e)
# setting a default image size large enough
plt.rcParams['figure.figsize'] = [10, 10]
# extract the mesh
mesh = model.getMesh()
conn = mesh.getConnectivity(aka._triangle_3)
nodes = mesh.getNodes()
triangles = tri.Triangulation(nodes[:, 0], nodes[:, 1], conn)
# plot stress field
plt.axes().set_aspect('equal')
ifield = model.getMaterial(0).getInternalReal(field)(aka._triangle_3)
disp = plt.tripcolor(
triangles, np.linalg.norm(ifield, axis=1))
cbar = plt.colorbar(disp, location="bottom")
cbar.set_label(f'{field} magnitude')
c4science · Help