Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-30 13:29:22 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
0bc9d27ffc794338a09ee9bbc9ed71a74cfa021f
JuliaFEM.jl/notebooks/2015-06-25-elasticity-solver-example.ipynb
T
Jukka Aho 9d199a5b65 performance improvements.
2015-12-02 17:42:48 +02:00

1242 lines
37 KiB
Plaintext
Raw Blame History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Solving elasticity problems using JuliaFEM\n",
"\n",
"Author(s): Jukka Aho\n",
"\n",
"**Abstract**: Solving elasticity equations using JuliaFEM.\n",
"\n",
"### Weak form\n",
"\n",
"Given function spaces\n",
"\\begin{align}\n",
"\\boldsymbol{\\mathcal{U}} & =\\left\\{ \\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\boldsymbol{u}\\left(\\boldsymbol{X},t\\right)=\\hat{\\boldsymbol{u}}\\left(\\boldsymbol{X},t\\right)\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\\\\\n",
"\\boldsymbol{\\mathcal{V}} & =\\left\\{ \\delta\\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\delta\\boldsymbol{u}\\left(\\boldsymbol{X}\\right)=0\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\n",
"\\end{align}\n",
"find $\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{U}}$ such that\n",
"\\begin{equation}\n",
"\\delta\\mathcal{W}:=\\int_{\\Omega_{0}}\\rho_{0}\\ddot{\\boldsymbol{u}}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}+\\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\hat{\\boldsymbol{b}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}-\\int_{\\Gamma_{\\sigma}}\\hat{\\boldsymbol{t}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}A_{0} =0 \\qquad\\forall\\delta\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{V}}\n",
"\\end{equation}\n",
"\n",
"### Some formulas\n",
"\\begin{align}\n",
"J & =\\det\\left(F\\right)\\\\\n",
"I_{c} & =\\mbox{tr}\\left(C\\right)\\\\\n",
"\\mathbf{C} & =\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}\\\\\n",
"\\mathbf{F} & =\\mathbf{I}+\\nabla\\mathbf{u}\\\\\n",
"\\mathbf{E} & =\\frac{1}{2}\\left(\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}-\\mathbf{I}\\right)\n",
"\\end{align}\n",
"\n",
"### Potential energy\n",
"\n",
"\\begin{equation}\n",
"\\underset{u\\in\\boldsymbol{\\mathcal{U}}}{\\min}\\Pi\\left(\\mathbf{u}\\right)\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"\\Pi\\left(\\mathbf{u}\\right)=\\int_{\\Omega}\\psi\\left(\\mathbf{u}\\right)-\\int_{\\Omega}\\hat{\\mathbf{b}}_{0}\\cdot\\mathbf{u}-\\int_{\\Gamma_{\\sigma}}\\hat{\\mathbf{t}}_{0}\\cdot\\mathbf{u}\\,\\mathrm{d}A_{0}\n",
"\\end{equation}\n",
"\n",
"### Material models\n",
"\n",
"https://en.wikipedia.org/wiki/Strain_energy_density_function\n",
"\n",
"Saint-Venant-Kirchhoff model https://en.wikipedia.org/wiki/Hyperelastic_material\n",
"\\begin{equation}\n",
"\\psi\\left(\\mathbf{E}\\right)=\\frac{\\lambda}{2}\\left[\\mbox{tr}\\left(\\mathbf{E}\\right)\\right]^{2}+\\mu\\mbox{tr}\\left(\\mathbf{E}^2\\right)\n",
"\\end{equation}\n",
"\n",
"neo-Hookean material https://en.wikipedia.org/wiki/Neo-Hookean_solid\n",
"\\begin{equation}\n",
"\\psi=\\frac{\\mu}{2}\\left(I_{c}-3\\right)-\\mu\\ln\\left(J\\right)+\\frac{\\lambda}{2}\\ln\\left(J\\right)^{2}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"addprocs(7)'\n",
"@everywhere using JuliaFEM"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM.Core: Element, Tri3, Tet4, Tri6, Tet10\n",
"using JuliaFEM.Core: ElasticityProblem, DirichletProblem\n",
"using JuliaFEM.Core: DirectSolver"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\",\"ELSET\"]\n",
"INFO: Parsing elements\n",
"INFO: 62454 elements found\n",
"INFO: Creating ELSET PISTON\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" "
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Parsing elements\n",
"INFO: 3494 elements found\n",
"INFO: Creating element set BC1\n",
"INFO: Creating element set BC2\n",
"INFO: model loaded.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 2.269207 seconds (14.52 M allocations: 433.635 MB, 13.38% gc time)\n"
]
}
],
"source": [
"@time begin\n",
" # Linear models\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_8789_P1.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_16436_P1.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_27343_P1.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_45510_P1.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_75470_P1.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/wrench/wrench_128903_P1.inp\")\n",
"\n",
" # Quadratic models\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_19611_P2.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_55950_P2.inp\")\n",
" model = open(JuliaFEM.Core.parse_abaqus, \"/Temp/piston_107168_P2.inp\")\n",
" #model = open(JuliaFEM.Core.parse_abaqus, \"../geometry/piston/piston_345757_P2.inp\")\n",
"\n",
" info(\"model loaded.\")\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: 65948 elements created.\n",
"INFO: all ready for solver\n"
]
}
],
"source": [
"perm = Dict{Int64, Int64}()\n",
"for (i, nid) in enumerate(keys(model[\"nodes\"]))\n",
" perm[nid] = i\n",
"end\n",
"\n",
"elements = Dict{Int64, Element}()\n",
"for (elid, node_ids) in model[\"elements\"]\n",
" connectivity = Int64[perm[nid] for nid in node_ids]\n",
" coords = Vector{Float64}[model[\"nodes\"][nid] for nid in node_ids]\n",
" elmap = Dict(3 => Tri3, 4 => Tet4, 6 => Tri6, 10 => Tet10)\n",
" for (dim, eltype) in elmap\n",
" if length(coords) == dim\n",
" element = eltype(connectivity)\n",
" element[\"geometry\"] = coords\n",
" element[\"id\"] = elid\n",
" elements[elid] = element\n",
" end\n",
" end\n",
"end\n",
"\n",
"info(\"$(length(elements)) elements created.\")\n",
"\n",
"problem = ElasticityProblem()\n",
"#for elid in model[\"elsets\"][\"WRENCH\"]\n",
"for elid in model[\"elsets\"][\"PISTON\"]\n",
" elements[elid][\"youngs modulus\"] = 210.0e3\n",
" elements[elid][\"poissons ratio\"] = 0.3\n",
" push!(problem, elements[elid])\n",
"end\n",
"\n",
"traction = Vector{Float64}[[0.0, 0.0, -10.0] for i in 1:6]\n",
"#traction = Vector{Float64}[[100.0, 100.0, 100.0] for i in 1:3]\n",
"#traction = Vector{Float64}[[0.0, -100.0, 20.0] for i in 1:3]\n",
"for elid in model[\"elsets\"][\"BC1\"]\n",
" elements[elid][\"displacement traction force\"] = traction\n",
" push!(problem, elements[elid])\n",
"end\n",
"\n",
"bc = DirichletProblem(\"displacement\", 3)\n",
"for elid in model[\"elsets\"][\"BC2\"]\n",
" elements[elid][\"displacement\"] = 0.0\n",
" push!(bc, elements[elid])\n",
"end\n",
"\n",
"solver = DirectSolver()\n",
"push!(solver, problem)\n",
"push!(solver, bc)\n",
"info(\"all ready for solver\")"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"empty!(model)\n",
"gc()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: # of field problems: 1\n",
"INFO: # of boundary problems: 1\n",
"INFO: Starting iteration 1\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 11.391177449683427\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 3.686000108718872\n",
"INFO: field assembly : 162.32500004768372\n",
"INFO: reduce stiffness matrix : 5.069999933242798\n",
"INFO: create sparse matrices : 0.1399998664855957\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 29.73900008201599\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.621999979019165\n",
"INFO: non-linear iteration : 204.1119999885559\n",
"INFO: Starting iteration 2\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 0.11391526987712108\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 1.9660000801086426\n",
"INFO: field assembly : 160.2829999923706\n",
"INFO: reduce stiffness matrix : 5.057999849319458\n",
"INFO: create sparse matrices : 0.1420001983642578\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.674999952316284\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.6159999370574951\n",
"INFO: non-linear iteration : 231.24799990653992\n",
"INFO: Starting iteration 3\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 0.0004840672223822576\n",
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 2.003999948501587\n",
"INFO: field assembly : 168.06199979782104\n",
"INFO: reduce stiffness matrix : 5.039999961853027\n",
"INFO: create sparse matrices : 0.14100003242492676\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.27800011634827\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.61899995803833\n",
"INFO: non-linear iteration : 238.61800003051758\n",
"INFO: Starting iteration 4\n",
"INFO: Assembling boundary problems...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: # of interface dofs: 8568\n",
"INFO: Assembling field problems...\n",
"INFO: Assembling body 1...\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: # of dofs in problem 1: 321504\n",
"INFO: Solving system\n",
"INFO: Solved, calculating interior dofs...\n",
"INFO: Problem solved. solution norm: 5.318293313386342e-10\n"
]
},
{
"data": {
"text/plain": [
"(4,true)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: timing info for non-linear iteration:\n",
"INFO: boundary assembly : 2.0319998264312744\n",
"INFO: field assembly : 166.31699991226196\n",
"INFO: reduce stiffness matrix : 5.141999959945679\n",
"INFO: create sparse matrices : 0.1419999599456787\n",
"INFO: dump matrices to disk : 0.0\n",
"INFO: solve problem : 60.019999980926514\n",
"INFO: back substitute : 0.0\n",
"INFO: update element data : 1.621999979019165\n",
"INFO: non-linear iteration : 236.77899980545044\n",
"INFO: solver finished in 912.1029999256134 seconds.\n"
]
}
],
"source": [
"solver.reduce_stiffness_matrix = false\n",
"#solver.dump_matrices = true\n",
"call(solver, 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: 1350 boundary dofs\n",
"INFO: Assemble: 10 % done\n",
"INFO: Assemble: 20 % done\n",
"INFO: Assemble: 30 % done\n",
"INFO: Assemble: 40 % done\n",
"INFO: Assemble: 50 % done\n",
"INFO: Assemble: 60 % done\n",
"INFO: Assemble: 70 % done\n",
"INFO: Assemble: 80 % done\n",
"INFO: Assemble: 90 % done\n",
"INFO: Assemble: 100 % done\n",
"INFO: 26367 field dofs\n",
"INFO: 25017 interior dofs\n",
"INFO: Assembly: 17.84435486793518\n",
"INFO: det = Inf\n",
"INFO: Factorization: 1.87113618850708\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 19.775338 seconds (152.31 M allocations: 8.037 GB, 10.16% gc time)\n"
]
}
],
"source": [
"using JuliaFEM.Core: assemble\n",
"\n",
"for elid in model[\"elsets\"][\"PISTON\"]\n",
" elements[elid][\"displacement\"] = (0.0 => Vector{Float64}[zeros(Float64, 3) for i = 1:4])\n",
"end\n",
"\n",
"function do_assembly()\n",
" t0 = time()\n",
" boundary_assembly = sum(map((p)->assemble(p, 0.0), solver.boundary_problems))\n",
" boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)\n",
" info(\"$(length(boundary_dofs)) boundary dofs\")\n",
" field_assembly = sum(map((p)->assemble(p, 0.0), solver.field_problems))\n",
" field_dofs = unique(field_assembly.stiffness_matrix.I)\n",
" info(\"$(length(field_dofs)) field dofs\")\n",
" interior_dofs = setdiff(field_dofs, boundary_dofs)\n",
" info(\"$(length(interior_dofs)) interior dofs\")\n",
" info(\"Assembly: \", time()-t0)\n",
"\n",
" t0 = time()\n",
" K = sparse(field_assembly.stiffness_matrix)\n",
" f = sparse(field_assembly.force_vector)\n",
" @assert maximum(abs(1/2*(K + K') - K)) < 1.0e-6\n",
" K = 1/2*(K + K')\n",
" Kii = K[interior_dofs, interior_dofs]\n",
" info(\"det = \", det(Kii))\n",
" F = cholfact(Kii)\n",
" Kbb = K[boundary_dofs, boundary_dofs]\n",
" Kib = K[interior_dofs, boundary_dofs]\n",
" fi = f[interior_dofs]\n",
" info(\"Factorization: \", time()-t0)\n",
" dim = size(K, 1)\n",
" return dim, F, Kbb, Kib, fi, boundary_dofs, interior_dofs\n",
"end\n",
"\n",
"@time dim, F, Kbb, Kib, fi, boundary_dofs, interior_dofs = do_assembly();"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"58833x1 sparse matrix with 0 Float64 entries:"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Kc = spzeros(dim, dim)\n",
"fc = spzeros(dim, 1)"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: 10 % done.\n",
"INFO: 20 % done.\n",
"INFO: 30 % done.\n",
"INFO: 40 % done.\n",
"INFO: 50 % done.\n",
"INFO: 60 % done.\n",
"INFO: 70 % done.\n",
"INFO: 80 % done.\n",
"INFO: 90 % done.\n",
"INFO: 100 % done.\n",
"INFO: norm(2) = 535881.6985084609\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 20.500665 seconds (13.36 M allocations: 3.063 GB, 3.34% gc time)\n"
]
}
],
"source": [
"@time begin\n",
" nb = length(boundary_dofs)\n",
" p = round(Int, nb/10)\n",
" Kd = zeros(nb, nb)\n",
" Vd = zeros(nb)\n",
" for bi in 1:nb\n",
" mod(bi, p) == 0 && info(round(Int, bi/nb*100), \" % done.\")\n",
" C = full(F \\ Kib[:, bi])\n",
" fill!(Vd, 0.0)\n",
" for bj in bi:nb\n",
" d = Kib[:, bj]\n",
" Kd[bj,bi] = dot(C[rowvals(d)], nonzeros(d))\n",
" Vd[bj] = dot(C[rowvals(d)], nonzeros(d))\n",
" end\n",
" Kd[:, bi] = Vd\n",
" end\n",
" Kd += tril(Kd, -1)'\n",
" Kc = spzeros(dim, dim)\n",
" Kc[boundary_dofs, boundary_dofs] = Kbb - Kd\n",
" info(\"norm(2) = \", norm(Kd))\n",
" @assert isapprox(norm(Kd), 535881.698508461)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"79868.97400906119"
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"maximum(abs(1/2*(Kd + Kd') - Kd))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"30% memory usage before start. 100% CPU usage. 85 % maximum memory usage."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_slices (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function get_slices(nb, n)\n",
"# n = 3\n",
"# nb = 25\n",
" kk = round(Int, collect(linspace(0, nb, n+1)))\n",
" return [kk[j]+1:kk[j+1] for j=1:length(kk)-1]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: 10.0 % done\n",
"INFO: 20.0 % done\n",
"INFO: 30.0 % done\n",
"INFO: 40.0 % done\n",
"INFO: 50.0 % done\n",
"INFO: 60.0 % done\n",
"INFO: 70.0 % done\n",
"INFO: 80.0 % done\n",
"INFO: 90.0 % done\n",
"INFO: 100.0 % done\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" 47"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: 3.525162568997601e6\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
".470587 seconds (251.48 k allocations: 10.678 GB, 5.06% gc time)\n"
]
}
],
"source": [
"@time begin\n",
" chunks = 10\n",
" nb = length(boundary_dofs)\n",
" kk = round(Int, collect(linspace(0, nb, chunks+1)))\n",
" sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]\n",
" Kc = spzeros(dim, dim)\n",
" for (k, sli) in enumerate(sl)\n",
" b1 = boundary_dofs[sli]\n",
" Sc = F \\ Kib[:,sli]\n",
" for slj in sl\n",
" b2 = boundary_dofs[slj]\n",
" Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc\n",
" end\n",
" info(round(k/chunks*100, 0), \" % done\")\n",
" end\n",
"end\n",
"info(norm(Kc, 1))\n",
"@assert isapprox(norm(Kc, 1), 3.525162568997601e6)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Saving results to file"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"<Grid Name=\"Grid\">\n",
" <Time Value=\"0\"/>\n",
"</Grid>\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n",
"temporal_collection = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n",
"grid = JuliaFEM.Postprocess.xdmf_new_grid(temporal_collection; time=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save geometry to xdmf file"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: number of nodes in model: 8789\n",
"INFO: Number of elements in model: 31437\n"
]
}
],
"source": [
"nnodes = length(model[\"nodes\"])\n",
"info(\"number of nodes in model: $nnodes\")\n",
"\n",
"X = zeros(3, nnodes)\n",
"for nid in keys(model[\"nodes\"])\n",
" X[:, perm[nid]] = model[\"nodes\"][nid]\n",
"end\n",
"\n",
"nelements = length(model[\"elsets\"][\"PISTON\"])\n",
"info(\"Number of elements in model: $nelements\")\n",
"elmap = zeros(Int64, 5, nelements)\n",
"#elmap[1,:] = 0x0026\n",
"elmap[1,:] = 0x6\n",
"for (i, elid) in enumerate(model[\"elsets\"][\"PISTON\"])\n",
" elmap[2:end,i] = Int64[perm[nid] for nid in model[\"elements\"][elid]]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.Postprocess.xdmf_new_mesh(grid, X, elmap)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Save nodal data to model"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(3,8789)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM.Core: get_connectivity\n",
"\n",
"u = zeros(3, nnodes)\n",
"\n",
"for element in values(elements)\n",
" isa(element, Element{Tet4}) || continue\n",
" connectivity = get_connectivity(element)\n",
" field = element[\"displacement\"](0.0)\n",
" for (i, nid) in enumerate(connectivity)\n",
" u[:, nid] = field[i]\n",
" end\n",
"end\n",
"size(u)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.Postprocess.xdmf_new_field(grid, \"Displacement\", \"nodes\", u)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1525505"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/piston_8789_P1.xmf\")"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3x10 Array{Float64,2}:\n",
" 0.0582911 0.339694 0.431768 … 0.124351 0.346492 0.000617283\n",
" 0.0900323 0.22174 0.433814 0.140503 0.275036 0.0044649 \n",
" 0.0838928 -0.148224 0.00262483 -0.11171 0.0843743 -0.0175412 "
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u[:, 1:10]"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x3 Array{Float64,2}:\n",
" -0.0471355 -0.468786 -0.222571"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"minimum(u, 2)'"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x3 Array{Float64,2}:\n",
" 0.511168 0.563004 0.66328"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"maximum(u, 2)'"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"38.52232721814439"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"piston_8789_P1_solution_norm = 38.52232721814439"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"@assert isapprox(norm(vec(u)), piston_8789_P1_solution_norm)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.1",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Reference in New Issue View Git Blame Copy Permalink