{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "5514008f-fcce-47bc-87aa-6733d377d745", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import akantu as aka\n", "import subprocess\n", "import pyvista as pv" ] }, { "cell_type": "code", "execution_count": 44, "id": "6b3751d1-1c0b-4875-8014-e5aa30683ca1", "metadata": { "tags": [] }, "outputs": [], "source": [ "class wall_with_rips:\n", " \n", " def __init__(self, name):\n", " \n", " self.name = name\n", " \n", " def initialize(self, length, height, t_w=0.03, t_r=0.03, l_r=0.8, n=10):\n", " \n", " self.l = length\n", " self.h = height\n", " self.t_w = t_w\n", " self.t_r = t_r\n", " self.l_r = l_r\n", " self.n = n\n", " \n", " self.dist = length / (n-1)\n", " \n", " def write_gmsh(self):\n", " \n", " # Geo file\n", " geo_file = \"\"\"\n", " SetFactory(\"OpenCASCADE\");\n", "\n", " length = {}; height = {}; t_w = {};\n", " t_r = {}; l_r = {}; dist = {};\"\"\".format(self.l, self.h, self.t_w, self.t_r, self.l_r, self.dist)\n", "\n", " # Add wall and side rips\n", " geo_file += \"\"\"\n", "\n", " Box(1) = {0, 0, 0, length, t_w, height};\n", "\n", " \"\"\"\n", "\n", " # Add rips according to the set distance\n", " for i in range(self.n):\n", " pos = i * self.dist\n", " # print(pos)\n", " geo_file += \"Box({}) = {{{}-t_r/2, 0, 0, t_r, l_r, height}}; \\n\".format(i+2, pos)\n", "\n", " arr = [\"Volume{{{}}}; \".format(j) for j in range(2, 2+self.n)]\n", " geo_file += \"\\nvol[] = BooleanUnion { Volume{1}; Delete; } { \" + \"\".join(arr) + \" Delete; };\" \n", "\n", " geo_file += \"\"\"\n", " \n", " Printf('%g', vol[0]);\n", "\n", " Physical Volume(\"Volume\") = {vol[0]};\n", " Physical Surface(\"front\") = {1};\n", " Physical Surface(\"bottom\") = {2};\n", " Physical Point(\"fixpoint\") = {1};\n", "\n", " Mesh.MeshSizeMin = 0.02;\n", " Mesh.MeshSizeMax = 0.2;\n", " Mesh.SaveAll = 1;\n", " Mesh.ElementOrder = 2;\n", "\n", " Mesh 3;\n", " \"\"\"\n", " \n", " geo_file += \"Save \\\"{}.msh\\\";\".format(self.name);\n", " \n", " open('{}.geo'.format(self.name), 'w').write(geo_file)\n", "\n", " \n", " def run_gmsh(self):\n", " \n", " ret = subprocess.run('gmsh '+self.name+'.geo -', shell=True)\n", " \n", " if ret.returncode:\n", " print(\"Beware, gmsh could not run: mesh is not regenerated\")\n", " else:\n", " print(\"Mesh generated\")\n", " \n", " \n", " \n", " \n", " # ------------------------------------------------------------------------------------- \n", " # run FEM calculation on wall with hydrostatic load\n", "\n", " def run_akantu(self, material_file):\n", "\n", " class hydrostatic_stress(aka.NeumannFunctor):\n", " def __init__(self, stress):\n", " super().__init__()\n", " self.stress = stress\n", "\n", " def __call__(self, quad_point, dual, coord, normals):\n", " z_coordinate = coord[aka._z]\n", " dual[:] = np.dot(self.stress, normals)*(height - z_coordinate)\n", " \n", " def get_mesh():\n", " \n", " mesh = aka.Mesh(3)\n", " mesh.read(self.name+'.msh')\n", " return mesh\n", "\n", " def applyBC():\n", " for dir in [aka._x, aka._y, aka._z]:\n", " model.applyBC(aka.FixedValue(0, dir), \"bottom\")\n", "\n", " stress = np.eye(3)*-1\n", " model.applyBC(hydrostatic_stress(stress), \"front\")\n", "\n", " def solve():\n", " solver = model.getNonLinearSolver()\n", " solver.set(\"max_iterations\", 10)\n", " solver.set(\"threshold\", 1e-4)\n", " solver.set(\"convergence_type\", aka.SolveConvergenceCriteria.residual)\n", " model.solveStep()\n", "\n", " def dump():\n", " model.setBaseName(self.name)\n", " model.addDumpFieldVector(\"displacement\")\n", " model.addDumpFieldVector(\"external_force\")\n", " model.addDumpField(\"strain\")\n", " model.addDumpField(\"stress\")\n", " model.addDumpField(\"blocked_dofs\")\n", " model.dump()\n", "\n", " aka.parseInput(material_file)\n", " mesh = get_mesh()\n", "\n", " model = aka.SolidMechanicsModel(mesh)\n", " model.initFull(_analysis_method=aka._static)\n", "\n", " applyBC()\n", " solve()\n", " dump()\n", " \n", " print('Model calculated!')\n", " \n", " \n", " \n", " def load_paraview(self):\n", " \n", " para = pv.read(f'paraview/'+self.name+f'_{0:04d}.pvtu')\n", " return para\n", " \n", " \n", "def plotResult(name, field='stress', warped=True):\n", "\n", " p = pv.Plotter(window_size=(600, 600))\n", " p.background_color = 'w'\n", "\n", " wall = pv.read(f'paraview/'+name+f'_{0:04d}.pvtu')\n", " wall.set_active_scalars('stress')\n", " wall_warped = wall.warp_by_vector('displacement')\n", " wall_warped.set_active_scalars('stress')\n", " \n", " if warped:\n", " p.add_mesh(wall_warped, color='grey', scalars=field, show_edges=False)\n", " else:\n", " p.add_mesh(wall, color='grey', scalars=field, show_edges=False)\n", " \n", " p.show(jupyter_backend='panel')\n", " \n", "\n", "def get_max_disp_z(data):\n", " disp_z = data['displacement'][:,2]\n", " return np.max(disp_z)\n", " \n", "def get_volume(data):\n", " return data.volume\n", "\n", "# -----------------------------------------------------------------------------------" ] }, { "cell_type": "code", "execution_count": 41, "id": "998d4512-8823-4a49-8ce0-589f1f42dc79", "metadata": {}, "outputs": [], "source": [ "# GEOMETRY\n", "# ----------------\n", "\n", "length = 3; \n", "height = 1.5;\n", "t_w = 0.02; \n", "t_r = 0.02;\n", "l_r = 0.3; \n", "n = 5; # minimum 2\n", "\n", "# MATERIAL\n", "# ----------------\n", "material_file = \"\"\"\n", "material elastic [\n", " name = polymer\n", " rho = 1000 # density\n", " E = 1e5 # young's modulus\n", " nu = 0.3 # poisson's ratio\n", "]\"\"\"\n", "\n", "# writing the material file\n", "open('material.dat', 'w').write(material_file)\n", "\n", "#reading the material file\n", "material_file = 'material.dat'" ] }, { "cell_type": "markdown", "id": "3bb8dca9-8b90-487b-b338-36c9d543f5ee", "metadata": {}, "source": [ "### Run the model" ] }, { "cell_type": "code", "execution_count": 42, "id": "3790d4eb-542b-4f0f-89e8-a4275dc8bd73", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Mesh generated\n" ] } ], "source": [ "wall = wall_with_rips(name='wall')\n", "\n", "wall.initialize(length=length, height=height, t_w=t_w, t_r=t_r, l_r=l_r, n = n)\n", "wall.write_gmsh()\n", "wall.run_gmsh()" ] }, { "cell_type": "code", "execution_count": 43, "id": "2361d9d2-a87a-460d-bbe8-5551847b07d5", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Model calculated!\n" ] }, { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "
\n", "
\n", "