Commit 2f2ec880 authored by florez's avatar florez
Browse files

wall with rips

parent 973adb6c
This diff is collapsed.
This diff is collapsed.
......@@ -3,9 +3,10 @@ SetFactory("OpenCASCADE");
//+
Cylinder(1) = {0, 0, 0, 0, 0, 1, 0.5, 2*Pi};
//+
Surface Loop(2) = {1, 2, 3};
//Surface Loop(2) = {1, 2, 3};
//+
Volume(2) = {2};
//Volume(2) = {2};
Physical Surface("bottom") = {3};
Physical Surface("top") = {2};
Physical Volume("steel") = {1};
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 1, 1, 1};
//+
Physical Surface("top") = {6};
//+
Physical Surface("bottom") = {5};
Physical Surface("top") = {6};
Physical Volume("cube") = {1};
box1[] = Point{1};
Printf("this is text %g", box1[2]);
Mesh.MeshSizeMin = 1;
Mesh.MeshSizeMax = 1;
Mesh.SaveAll = 1;
Mesh.ElementOrder = 1;
Mesh 3;
Save "cube.msh";
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
3
2 1 "bottom"
2 2 "top"
3 3 "cube"
$EndPhysicalNames
$Entities
8 12 6 1
1 0 0 1 0
2 0 0 0 0
3 0 1 1 0
4 0 1 0 0
5 1 0 1 0
6 1 0 0 0
7 1 1 1 0
8 1 1 0 0
1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1
2 -1e-07 -9.999999994736442e-08 0.9999999000000001 1e-07 1.0000001 1.0000001 0 2 1 -3
3 -1e-07 0.9999999000000001 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 2 4 -3
4 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 1.0000001 1e-07 0 2 2 -4
5 0.9999999000000001 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 2 6 -5
6 0.9999999000000001 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 5 -7
7 0.9999999000000001 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 2 8 -7
8 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 0 2 6 -8
9 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 2 -6
10 -9.999999994736442e-08 -1e-07 0.9999999000000001 1.0000001 1e-07 1.0000001 0 2 1 -5
11 -9.999999994736442e-08 0.9999999000000001 -1e-07 1.0000001 1.0000001 1e-07 0 2 4 -8
12 -9.999999994736442e-08 0.9999999000000001 0.9999999000000001 1.0000001 1.0000001 1.0000001 0 2 3 -7
1 -1e-07 -9.999999994736442e-08 -9.999999994736442e-08 1e-07 1.0000001 1.0000001 0 4 1 2 -3 -4
2 0.9999999000000001 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 5 6 -7 -8
3 -9.999999994736442e-08 -1e-07 -9.999999994736442e-08 1.0000001 1e-07 1.0000001 0 4 9 5 -10 -1
4 -9.999999994736442e-08 0.9999999000000001 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 0 4 11 7 -12 -3
5 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 1.0000001 1e-07 1 1 4 4 11 -8 -9
6 -9.999999994736442e-08 -9.999999994736442e-08 0.9999999000000001 1.0000001 1.0000001 1.0000001 1 2 4 2 12 -6 -10
1 -9.999999994736442e-08 -9.999999994736442e-08 -9.999999994736442e-08 1.0000001 1.0000001 1.0000001 1 3 6 1 2 3 4 5 6
$EndEntities
$Nodes
27 14 1 14
0 1 0 1
1
0 0 1
0 2 0 1
2
0 0 0
0 3 0 1
3
0 1 1
0 4 0 1
4
0 1 0
0 5 0 1
5
1 0 1
0 6 0 1
6
1 0 0
0 7 0 1
7
1 1 1
0 8 0 1
8
1 1 0
1 1 0 0
1 2 0 0
1 3 0 0
1 4 0 0
1 5 0 0
1 6 0 0
1 7 0 0
1 8 0 0
1 9 0 0
1 10 0 0
1 11 0 0
1 12 0 0
2 1 0 1
9
0 0.5 0.5
2 2 0 1
10
1 0.5 0.5
2 3 0 1
11
0.5 0 0.5
2 4 0 1
12
0.5 1 0.5
2 5 0 1
13
0.5 0.5 0
2 6 0 1
14
0.5 0.5 1
3 1 0 0
$EndNodes
$Elements
27 68 1 68
0 1 15 1
1 1
0 2 15 1
2 2
0 3 15 1
3 3
0 4 15 1
4 4
0 5 15 1
5 5
0 6 15 1
6 6
0 7 15 1
7 7
0 8 15 1
8 8
1 1 1 1
9 2 1
1 2 1 1
10 1 3
1 3 1 1
11 4 3
1 4 1 1
12 2 4
1 5 1 1
13 6 5
1 6 1 1
14 5 7
1 7 1 1
15 8 7
1 8 1 1
16 6 8
1 9 1 1
17 2 6
1 10 1 1
18 1 5
1 11 1 1
19 4 8
1 12 1 1
20 3 7
2 1 2 4
21 2 1 9
22 1 3 9
23 4 2 9
24 3 4 9
2 2 2 4
25 6 10 5
26 5 10 7
27 8 10 6
28 7 10 8
2 3 2 4
29 1 2 11
30 5 1 11
31 2 6 11
32 6 5 11
2 4 2 4
33 3 12 4
34 7 12 3
35 4 12 8
36 8 12 7
2 5 2 4
37 2 4 13
38 6 2 13
39 4 8 13
40 8 6 13
2 6 2 4
41 1 14 3
42 5 14 1
43 3 14 7
44 7 14 5
3 1 4 24
45 10 11 12 13
46 9 12 14 11
47 12 14 11 10
48 9 12 11 13
49 2 9 1 11
50 1 9 3 14
51 11 14 1 5
52 4 9 12 3
53 2 4 9 13
54 12 3 14 7
55 5 10 14 7
56 7 10 12 8
57 11 10 5 6
58 12 4 8 13
59 13 8 10 6
60 13 11 2 6
61 14 1 9 11
62 9 12 3 14
63 11 2 9 13
64 12 9 4 13
65 8 10 12 13
66 14 11 10 5
67 12 14 10 7
68 13 10 11 6
$EndElements
material elastic [
name = steel
rho = 7800 # density
E = 210e6 # young's modulus
nu = 0.3 # poisson's ratio
]
\ No newline at end of file
......@@ -3,191 +3,200 @@
{
"cell_type": "code",
"execution_count": 1,
"id": "1dbbab78-6ef1-4e97-8440-f6db326755e9",
"id": "3a4b7e63-10d0-4af2-bc97-5360ea402ff5",
"metadata": {},
"outputs": [],
"source": [
"import akantu as aka\n",
"import numpy as np\n",
"import pygmsh\n",
"\n",
"# import the pyplot submodule to draw figures\n",
"import matplotlib.pyplot as plt\n",
"# import triangluation routine to plot meshes\n",
"import matplotlib.tri as tri\n",
"# setting a default image size large enough\n",
"plt.rcParams['figure.figsize'] = [10, 10]"
"import akantu as aka"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "4f3b63d5-abe9-43ce-9eca-80cd47fb4d49",
"id": "91261137-56d5-49b9-b386-dfb54d40da2a",
"metadata": {},
"outputs": [],
"source": [
"# reading the mesh\n",
"spatial_dimension = 3 \n",
"mesh_file = 'shell_cyl_2.msh'\n",
"mesh = aka.Mesh(spatial_dimension)\n",
"mesh.read(mesh_file)\n"
"# import subprocess\n",
"\n",
"# ret = subprocess.run(\"gmsh -3 -order 1 -o out.msh in.geo\", shell=True)\n",
"# if ret.returncode:\n",
"# print(\"Beware, gmsh could not run: mesh is not regenerated\")\n",
"# else:\n",
"# print(\"Mesh generated\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "ac11734f-f415-43aa-a1a3-92edca236554",
"id": "c1339e28-e6b1-4155-b25e-f3b848025838",
"metadata": {},
"outputs": [],
"source": [
"material_file = \"\"\"\n",
"material elastic [\n",
" name = steel\n",
" rho = 1 # density\n",
" E = 1 # young's modulus\n",
" rho = 7800 # density\n",
" E = 210e6 # 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",
"#reading the material file\n",
"material_file = 'material.dat'\n",
"material_file = 'material.dat'"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "30c1bb68-80b2-41ea-9aa6-a08daf741c9b",
"metadata": {},
"outputs": [],
"source": [
"aka.parseInput(material_file)\n",
"\n",
"# creating the solid mechanics model\n",
"model = aka.SolidMechanicsModel(mesh)\n",
"spatial_dimension = 3\n",
"mesh = aka.Mesh(spatial_dimension)\n",
"mesh.read('cube.msh')\n",
"\n",
"# initialize a static solver\n",
"model = aka.SolidMechanicsModel(mesh)\n",
"model.initFull(_analysis_method=aka._static)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "e2379428-0b40-406f-bfed-160eb297f730",
"execution_count": 5,
"id": "9bb306b5-2e9c-4b3e-aaeb-347e7a1ee5f8",
"metadata": {},
"outputs": [],
"source": [
"# Bottom ring 'bottom', top ring 'top', cylinder surface 'surface'\n",
"\n",
"# set the displacement/Dirichlet boundary conditions\n",
"for dir in [aka._x, aka._y, aka._z]:\n",
" model.applyBC(aka.FixedValue(0, dir), \"bottom\")\n",
" model.applyBC(aka.FixedValue(0.0, dir), \"bottom\")\n",
"\n",
"\n",
"\n",
"# set the force/Neumann boundary conditions\n",
"\n",
"trac = np.eye(3) # Newtons/m^2\n",
"trac = [0,0,1]\n",
"\n",
"model.applyBC(aka.FromTraction(trac), \"top\")\n"
"# trac = np.eye(3)\n",
"trac = [0, 0, 1e8]\n",
"model.applyBC(aka.FromTraction(trac), \"top\")\n",
"np.set_printoptions(threshold=10)\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "6d82b932-ca93-45eb-8e3e-c0ec2c24f381",
"execution_count": 6,
"id": "4d56ac9e-37db-48d3-8366-f8f8cb000bb6",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0., 0., 0.],\n",
" [0., 0., 0.],\n",
" [0., 0., 0.],\n",
" ...,\n",
" [0., 0., 0.],\n",
" [0., 0., 0.],\n",
" [0., 0., 0.]])"
"<py11_akantu.Material at 0x7f2145e75730>"
]
},
"execution_count": 5,
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.set_printoptions(threshold=20)\n",
"steel = model.getMaterial('steel')\n",
"steel"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "845b08f6-cda8-4a72-8fa0-f1ba47fb432f",
"metadata": {},
"outputs": [],
"source": [
"solver = model.getNonLinearSolver()\n",
"solver.set(\"max_iterations\", 6)\n",
"solver.set(\"threshold\", 1e-8)\n",
"solver.set(\"convergence_type\", aka.SolveConvergenceCriteria.residual)\n",
"\n",
"model.getExternalForce()[:]\n",
"# np.size(model.getExternalForce(),0)\n",
"# model.getBlockedDOFs()"
"model.solveStep()"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "09675758-79eb-4314-b23f-f34fe359dc9f",
"execution_count": 10,
"id": "2a96374a-1ade-42cf-8024-69c1dd27f02e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[0., 0., 0.],\n",
" [0., 0., 0.],\n",
" [0., 0., 0.],\n",
" ...,\n",
" [0., 0., 0.],\n",
" [0., 0., 0.],\n",
" [0., 0., 0.]])"
"24"
]
},
"execution_count": 6,
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# configure the linear algebra solver\n",
"solver = model.getNonLinearSolver()\n",
"solver.set(\"max_iterations\", 3)\n",
"solver.set(\"threshold\", 1e-8)\n",
"solver.set(\"convergence_type\", aka.SolveConvergenceCriteria.residual)\n",
"\n",
"# compute the solution\n",
"model.solveStep()\n",
"\n",
"u = model.getDisplacement()\n",
"np.set_printoptions(threshold=120)\n",
"u"
"mesh.getConnectivity(aka._tetrahedron_4)\n",
"stress_field = model.getMaterial(0).getStress(aka._tetrahedron_4)\n",
"np.set_printoptions(threshold=100)\n",
"stress_field\n",
"mesh.getNbElement()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a81b557e-7adb-452a-9982-efeef3070962",
"execution_count": null,
"id": "23e23201-534c-4b3b-8121-6f8504b57243",
"metadata": {},
"outputs": [
{
"ename": "Exception",
"evalue": "akantu::debug::Exception : No element of type (not_ghost:_cohesive_2d_4) in this const ElementTypeMapArray<unsigned int> class(\"mesh:connectivities\") [/tmp/app/spack-stage/spack-stage-akantu-master-t5tv2fws2j5an4dr4sqytng27g7b3fyr/spack-src/src/mesh/element_type_map_tmpl.hh:293]",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mException\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-7-254e5de967ab>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;31m# generate paraview files\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mmodel\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdump\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mException\u001b[0m: akantu::debug::Exception : No element of type (not_ghost:_cohesive_2d_4) in this const ElementTypeMapArray<unsigned int> class(\"mesh:connectivities\") [/tmp/app/spack-stage/spack-stage-akantu-master-t5tv2fws2j5an4dr4sqytng27g7b3fyr/spack-src/src/mesh/element_type_map_tmpl.hh:293]"
"outputs": [],
"source": [
"np.argmax(stress_field,0)\n",
"np.size(stress_field,0)\n",
"stress_normal = stress_field[:,(0,4,8)]\n",
"np.size(stress_normal,0)"
]
}
],
},
{
"cell_type": "code",
"execution_count": null,
"id": "60a70c37-31bf-48bf-a4ac-8333abd62e93",
"metadata": {},
"outputs": [],
"source": [
"# specify what field to output into paraview files\n",
"model.setBaseName(\"shell_cyl\")\n",
"model.setBaseName(\"cube\")\n",
"model.addDumpFieldVector(\"displacement\")\n",
"model.addDumpFieldVector(\"external_force\")\n",
"model.addDumpField(\"strain\")\n",
"model.addDumpField(\"stress\")\n",
"model.addDumpField(\"blocked_dofs\")\n",
"\n",
"# generate paraview files\n",
"model.dump()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b9caba6a-e863-495b-989c-fd3bc3cf44ab",
"id": "ac9c3e33-7117-4809-a433-3e7fef3a3344",
"metadata": {},
"outputs": [],
"source": [
"import pyvista as pv\n",
"\n",
"p = pv.Plotter(off_screen=False, notebook=False)\n",
"p.background_color = 'white'\n",
"\n",
"cyl_msh = pv.read(f'paraview/cube_{0:04d}.pvtu')\n",
"cyl_msh.set_active_scalars('displacement')\n",
"cyl_warped = cyl_msh.warp_by_vector('displacement')\n",
"cyl_warped.set_active_scalars('displacement')\n",
"\n",
"cyl_warped.plot()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "738e4f52-bea5-4ae3-9454-2456af57c208",
"metadata": {},
"outputs": [],
"source": []
......
%% Cell type:code id:1dbbab78-6ef1-4e97-8440-f6db326755e9 tags:
%% Cell type:code id:3a4b7e63-10d0-4af2-bc97-5360ea402ff5 tags:
``` python
import akantu as aka
import numpy as np
import pygmsh
# import the pyplot submodule to draw figures
import matplotlib.pyplot as plt
# import triangluation routine to plot meshes
import matplotlib.tri as tri
# setting a default image size large enough
plt.rcParams['figure.figsize'] = [10, 10]
import akantu as aka
```
%% Cell type:code id:4f3b63d5-abe9-43ce-9eca-80cd47fb4d49 tags:
%% Cell type:code id:91261137-56d5-49b9-b386-dfb54d40da2a tags:
``` python
# reading the mesh
spatial_dimension = 3
mesh_file = 'shell_cyl_2.msh'
mesh = aka.Mesh(spatial_dimension)
mesh.read(mesh_file)
# import subprocess
# ret = subprocess.run("gmsh -3 -order 1 -o out.msh in.geo", shell=True)
# if ret.returncode:
# print("Beware, gmsh could not run: mesh is not regenerated")
# else:
# print("Mesh generated")
```
%% Cell type:code id:ac11734f-f415-43aa-a1a3-92edca236554 tags:
%% Cell type:code id:c1339e28-e6b1-4155-b25e-f3b848025838 tags:
``` python
material_file = """
material elastic [
name = steel
rho = 1 # density
E = 1 # young's modulus
rho = 7800 # density
E = 210e6 # young's modulus
nu = 0.3 # poisson's ratio
]"""
# writing the material file
open('material.dat', 'w').write(material_file)
#reading the material file
material_file = 'material.dat'
```
%% Cell type:code id:30c1bb68-80b2-41ea-9aa6-a08daf741c9b tags:
``` python
aka.parseInput(material_file)
# creating the solid mechanics model
model = aka.SolidMechanicsModel(mesh)
spatial_dimension = 3
mesh = aka.Mesh(spatial_dimension)
mesh.read('cube.msh')
# initialize a static solver
model = aka.SolidMechanicsModel(mesh)
model.initFull(_analysis_method=aka._static)
```
%% Cell type:code id:e2379428-0b40-406f-bfed-160eb297f730 tags:
%% Cell type:code id:9bb306b5-2e9c-4b3e-aaeb-347e7a1ee5f8 tags:
``` python
# Bottom ring 'bottom', top ring 'top', cylinder surface 'surface'
# set the displacement/Dirichlet boundary conditions
for dir in [aka._x, aka._y, aka._z]:
model.applyBC(aka.FixedValue(0, dir), "bottom")
# set the force/Neumann boundary conditions
trac = np.eye(3) # Newtons/m^2
trac = [0,0,1]
model.applyBC(aka.FixedValue(0.0, dir), "bottom")
# trac = np.eye(3)
trac = [0, 0, 1e8]
model.applyBC(aka.FromTraction(trac), "top")
np.set_printoptions(threshold=10)
```
%% Cell type:code id:6d82b932-ca93-45eb-8e3e-c0ec2c24f381 tags:
%% Cell type:code id:4d56ac9e-37db-48d3-8366-f8f8cb000bb6 tags:
``` python
np.set_printoptions(threshold=20)
model.getExternalForce()[:]
# np.size(model.getExternalForce(),0)
# model.getBlockedDOFs()
steel = model.getMaterial('steel')
steel
```
%%%% Output: execute_result
array([[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
...,
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]])
<py11_akantu.Material at 0x7f2145e75730>
%% Cell type:code id:09675758-79eb-4314-b23f-f34fe359dc9f tags:
%% Cell type:code id:845b08f6-cda8-4a72-8fa0-f1ba47fb432f tags:
``` python
# configure the linear algebra solver
solver = model.getNonLinearSolver()
solver.set("max_iterations", 3)
solver.set("max_iterations", 6