Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-06 04:21:33 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
master
JuliaFEM.jl/docs/tutorials/2015-05-23-one-element-solution.ipynb
T

921 lines
67 KiB
Plaintext
Raw Permalink Normal View History

Added python conf files
2015-06-26 19:53:42 +03:00
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#One element solution with and without nonlinear geometrics.\n",
"\n",
"Author: Jukka Aho\n",
"\n",
"Email: <jukka.aho@kapsi.fi>"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0x00023c5cd89a9c57"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"tic()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We consider here linear quadrangle element and evaluate solution for it. It has been verified using Elmer that the correct solution for minimum displacement in tip is $u_\\mathrm{min,nonlinear}$ = -2.22224475 for geometrically nonlinear case and $u_\\mathrm{min,linear} = -2.18056991$ for linear case.\n",
"\n",
"Solution domain is $\\Omega = [0, 10] \\times [0, 1] \\in \\mathbb{R}^2$ and we are having $F=-2 \\mathrm{N}$ nodal force in $(10, 1)$. $E=90$, $\\nu=0.25$. Left side is fixed, $u|_{x=0}=0$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Geometry, material definitions and shape functions:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(36.0,24.0)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X = [[0.0 10.0 10.0 0.0],\n",
" [0.0 0.0 1.0 1.0]]\n",
"\n",
"# Partial derivatives of bilinear Lagrange polynomials\n",
"dNdξ(ξ) = [[-(1-ξ[2])/4.0 -(1-ξ[1])/4.0],\n",
" [ (1-ξ[2])/4.0 -(1+ξ[1])/4.0],\n",
" [ (1+ξ[2])/4.0 (1+ξ[1])/4.0],\n",
" [-(1+ξ[2])/4.0 (1-ξ[1])/4.0]] \n",
"\n",
"E = 90\n",
"ν = 0.25\n",
"μ = E/(2*(1+ν))\n",
"λ = E*ν/((1+ν)*(1-2*ν))\n",
"λ = 2*λ*μ/(λ + 2*μ)\n",
"μ, λ"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Internal force vector:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"f_int (generic function with 2 methods)"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f_int(X, u, dNdξ, λ, μ, dim=2)\n",
" T = zeros(size(X))\n",
" I = eye(dim)\n",
" \n",
" function J(ξ)\n",
" Jᵀ = X*dNdξ(ξ)\n",
" ∇N = inv(Jᵀ)*dNdξ(ξ)'\n",
" ∇u = u*∇N'\n",
" F = I + ∇u # Deformation gradient\n",
" E = 1/2*(∇u' + ∇u + ∇u'*∇u) # Green-Lagrange strain tensor\n",
" P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor\n",
" S = F*P # PK2 stress tensor\n",
" return S*∇N*det(Jᵀ)\n",
" end\n",
"\n",
" a = 1/sqrt(3)\n",
" ipoints = [[-a -a], [a -a], [a a], [-a a]]\n",
" iweights = [1 1 1 1]\n",
"\n",
" for m = 1:length(iweights)\n",
" w = iweights[m]\n",
" ξ = ipoints[m, :]\n",
" T += w*J(ξ)\n",
" end\n",
"\n",
" return T\n",
"\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Solution"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter f(x) inf-norm Step 2-norm \n",
"------ -------------- --------------\n",
" 0 2.000000e+00 NaN\n",
" 1 1.000142e+00 2.257867e-03\n",
" 2 9.933192e-01 4.422932e-03\n",
" 3 9.961075e-01 8.847632e-03\n",
" 4 9.855335e-01 1.769324e-02\n",
" 5 9.738205e-01 3.538698e-02\n",
" 6 9.612786e-01 7.077382e-02\n",
" 7 9.404424e-01 1.415511e-01\n",
" 8 9.493268e-01 2.831322e-01\n",
" 9 1.157757e+00 5.665408e-01\n",
" 10 6.914588e-01 2.840319e-01\n",
" 11 6.914588e-01 0.000000e+00\n",
" 12 5.805109e-01 2.845511e-01\n",
" 13 5.805109e-01 0.000000e+00\n",
" 14 4.993552e-01 2.851202e-01\n",
" 15 4.993552e-01 0.000000e+00\n",
" 16 4.100475e-01 2.858343e-01\n",
" 17 4.100475e-01 0.000000e+00\n",
" 18 3.195632e-01 2.866934e-01\n",
" 19 3.195632e-01 0.000000e+00\n",
" 20 2.284246e-01 2.877431e-01\n",
" 21 1.444984e-01 3.012556e-01\n",
" 22 8.015667e-05 4.216354e-03\n",
" 23 5.099410e-10 1.611875e-05\n",
"Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Trust-region with dogleg and autoscaling\n",
" * Starting Point: [0.0,0.0,0.0,0.0]\n",
" * Zero: [-0.39914506093725427,-2.1779892316773526,-0.0722858269436904,-2.2222447543690937]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 23\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 24\n",
" * Jacobian Calls (df/dx): 19\n"
]
},
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" 0.0 -0.399145 -0.0722858 0.0\n",
" 0.0 -2.17799 -2.22224 0.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using NLsolve\n",
"\n",
"function solve()\n",
"\n",
" free_dofs = [3, 4, 5, 6]\n",
" u = zeros(2, 4)\n",
" F = zeros(2, 4)\n",
" F[2, 3] = -2.0\n",
"\n",
" function f(u_)\n",
" u[free_dofs] = u_\n",
" T = f_int(X, u, dNdξ, λ, μ)\n",
" R = T-F\n",
" return R[free_dofs]\n",
" end\n",
"\n",
" # Go!\n",
" sol = nlsolve(not_in_place(f), zeros(4), store_trace=true, show_trace=true)\n",
" println(sol)\n",
" u[free_dofs] = sol.zero\n",
" return u\n",
"end\n",
"\n",
"u = solve()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Difference to Elmer solution is"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4.3690935314089074e-9"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abs(u[2, 3] - -2.22224475)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Linear solution is achieved by assuming small deformation, i.e. :\n",
" \n",
" F = I + ∇u ≈ I\n",
" E = 1/2*(∇u' + ∇u + ∇u'*∇u) ≈ 1/2*(∇u' + ∇u)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**OFFTOPIC:** Here we should put something like\n",
"\n",
" using JuliaFEM\n",
" test_allclose(u[2, 3], -2.22224457, tol = ...)\n",
"\n",
"and get a nice green \"**PASSED**\" if everything is ok at this point, otherwise assertionerror"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Notes about linearization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For more robust solution, we shoud get analytical Jacobian for nonlinear system of equations. If not given it's probably some sort of finite difference approximation, i.e."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 119.461 13.4859 61.8587 … -26.0246 -113.693 2.80636\n",
" 13.486 309.575 15.6443 -148.518 -3.10971 -316.139 \n",
" 61.8587 15.6443 130.663 -53.2775 -60.0524 1.5741 \n",
" 9.73232 155.082 36.0592 -303.876 1.51268 -163.427 \n",
" -67.6267 -26.0207 -132.47 63.0313 56.5391 10.2937 \n",
" -26.0246 -148.518 -53.2777 … 299.44 16.2711 152.956 \n",
" -113.693 -3.10985 -60.0524 16.2711 117.207 -14.6739 \n",
" 2.80624 -316.14 1.57406 152.956 -14.674 326.611 "
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function Lin(f, h=1.0e-6)\n",
"\n",
" function D(x)\n",
" J_approx = zeros((8, 8))\n",
" for i=1:8\n",
" Δx = zeros(2, 4)\n",
" Δx[i] += h\n",
" Δf = f(x+Δx) - f(x)\n",
" J_approx[i, :] = Δf[:]/h\n",
" end\n",
" return J_approx\n",
" end\n",
"\n",
" return D\n",
"\n",
"end\n",
"\n",
"J_approx = Lin(u -> f_int(X, u, dNdξ, λ, μ))\n",
"J_approx(u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A more elegant solution is to use automatic differentiation to get analytical tangent stiffness matrix"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 119.461 13.486 61.8587 … -26.0246 -113.693 2.80625\n",
" 13.486 309.576 15.6443 -148.518 -3.10973 -316.14 \n",
" 61.8587 15.6443 130.663 -53.2776 -60.0524 1.57406\n",
" 9.73233 155.082 36.0593 -303.877 1.51267 -163.427 \n",
" -67.6267 -26.0206 -132.47 63.0312 56.539 10.2937 \n",
" -26.0246 -148.518 -53.2776 … 299.439 16.271 152.956 \n",
" -113.693 -3.10973 -60.0524 16.271 117.207 -14.674 \n",
" 2.80625 -316.14 1.57406 152.956 -14.674 326.611 "
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using ForwardDiff\n",
"\n",
"function T!(u, T)\n",
" T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ)\n",
"end\n",
"\n",
"J_analytical = forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)\n",
"J_analytical(reshape(u, 8))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Difference is in smooth situations quite small, but analytical is of course faster."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.000998738079715485"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(J_approx(u) - J_analytical(reshape(u, 8)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we use this analytical Jacobian we should get quadratic convergence when using Newton-Rhapson method:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iter f(x) inf-norm Step 2-norm \n",
"------ -------------- --------------\n",
" 0 2.000000e+00 NaN\n",
" 1 1.985767e+00 3.756746e-01\n",
" 2 1.933311e+00 3.652956e-01\n",
" 3 1.855451e+00 3.896874e-01\n",
" 4 1.749547e+00 4.460338e-01\n",
" 5 1.581001e+00 5.270850e-01\n",
" 6 1.331440e+00 8.104083e-01\n",
" 7 8.751161e-02 2.235404e-01\n",
" 8 2.570554e-04 1.112294e-02\n",
" 9 3.202487e-09 3.961907e-05\n",
"Results of Nonlinear Solver Algorithm\n",
" * Algorithm: Newton with line-search\n",
" * Starting Point: [0.0,0.0,0.0,0.0]\n",
" * Zero: [-0.3991450608223195,-2.1779892314495743,-0.072285826868806,-2.22224475412425]\n",
" * Inf-norm of residuals: 0.000000\n",
" * Iterations: 9\n",
" * Convergence: true\n",
" * |x - x'| < 0.0e+00: false\n",
" * |f(x)| < 1.0e-08: true\n",
" * Function Calls (f): 110\n",
" * Jacobian Calls (df/dx): 18\n"
]
},
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" 0.0 -0.399145 -0.0722858 0.0\n",
" 0.0 -2.17799 -2.22224 0.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function solve2()\n",
"\n",
" free_dofs = [3, 4, 5, 6]\n",
" u = zeros(2, 4)\n",
" F = zeros(2, 4)\n",
" F[2, 3] = -2.0\n",
"\n",
" function f(u_)\n",
" u[free_dofs] = u_\n",
" T = f_int(X, u, dNdξ, λ, μ)\n",
" R = T-F\n",
" return R[free_dofs]\n",
" end\n",
"\n",
" function g(u_)\n",
" u[free_dofs] = u_\n",
" J = J_analytical(reshape(u, 8))\n",
" return J[free_dofs, free_dofs]\n",
" end\n",
"\n",
" # Go!\n",
" sol = nlsolve(not_in_place(f, g), zeros(4), method = :newton, store_trace=true, show_trace=true)\n",
" println(sol)\n",
" u[free_dofs] = sol.zero\n",
" return u\n",
"end\n",
"\n",
"u = solve2()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clearly the last iterations are converging quadratic speeds. Own naive implementation of solver would be:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"norm = 3.0900221367289094\n",
"norm = 0.3212131602153428\n",
"norm = 0.04043178194002758\n",
"norm = 0.0009291101052123099\n",
"norm = 1.5638899074064188e-7\n",
"norm = 1.036219989655546e-14\n",
"norm = 3.130358566071589e-16\n"
]
},
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" 0.0 -0.399145 -0.0722858 0.0\n",
" 0.0 -2.17799 -2.22224 0.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function solve3()\n",
"\n",
" free_dofs = [3, 4, 5, 6]\n",
" u = zeros(2, 4)\n",
" F = zeros(2, 4)\n",
" F[2, 3] = -2.0\n",
"\n",
" function f(u_)\n",
" u[free_dofs] = u_\n",
" T = f_int(X, u, dNdξ, λ, μ)\n",
" R = T-F\n",
" return R[free_dofs]\n",
" end\n",
"\n",
" function g(u_)\n",
" u[free_dofs] = u_\n",
" J = J_analytical(reshape(u, 8))\n",
" return J[free_dofs, free_dofs]\n",
" end\n",
"\n",
" u2 = zeros(4)\n",
" for i=1:7\n",
" du = g(u2) \\ -f(u2)\n",
" u2 += du\n",
" println(\"norm = \",norm(du))\n",
" end\n",
" u[free_dofs] = u2\n",
" return u\n",
"end\n",
"\n",
"u = solve3()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I don't know what NLsolve does because it takes a bit longer to converge?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparison of results"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At end we do compare results between linear and nonlinear systems and put all together."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Float64,2}:\n",
" 0.0 -10.0735 -8.83486 0.0\n",
" 0.0 -10.4167 -11.8171 0.0"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function f_int(X, u, dNdξ, λ, μ, dim=2; nlgeom=true)\n",
" T = zeros(size(X))\n",
" I = eye(dim)\n",
"\n",
" function J(ξ)\n",
" Jᵀ = X*dNdξ(ξ)\n",
" ∇N = inv(Jᵀ)*dNdξ(ξ)'\n",
" ∇u = u*∇N'\n",
" F = I # Deformation gradient for small deformations\n",
" E = 1/2*(∇u' + ∇u) # Small strain tensor\n",
" if nlgeom # .. if nonlinear geometry\n",
" F += ∇u # .. add nonlinear part to deformation gradient, and\n",
" E += 1/2*(∇u'*∇u) # add nonlinear part to strain tensor => Green-Lagrange strain tensor\n",
" end\n",
" P = λ*trace(E)*I + 2*μ*E # PK1 stress tensor\n",
" S = F*P # PK2 stress tensor\n",
" return S*∇N*det(Jᵀ)\n",
" end\n",
"\n",
" a = 1/sqrt(3)\n",
" ipoints = [[-a -a], [a -a], [a a], [-a a]]\n",
" iweights = [1 1 1 1]\n",
"\n",
" for m = 1:length(iweights)\n",
" w = iweights[m]\n",
" ξ = ipoints[m, :]\n",
" T += w*J(ξ)\n",
" end\n",
"\n",
" return T\n",
"\n",
"end\n",
"\n",
"function Lin(f, h=1.0e-6)\n",
"\n",
" function D(x)\n",
" J_approx = zeros((8, 8))\n",
" for i=1:8\n",
" Δx = zeros(2, 4)\n",
" Δx[i] += h\n",
" Δf = f(x+Δx) - f(x)\n",
" J_approx[i, :] = Δf[:]/h\n",
" end\n",
" return J_approx\n",
" end\n",
"\n",
" return D\n",
"\n",
"end\n",
"\n",
"function solve4!(F, u; nlgeom=true, analytical_jacobian=false)\n",
"\n",
" free_dofs = [3, 4, 5, 6]\n",
"\n",
" # Tangent stiffness matrix\n",
" function T!(u, T)\n",
" T[:] = f_int(X, reshape(u, 2, 4), dNdξ, λ, μ, nlgeom=nlgeom)\n",
" end\n",
" if analytical_jacobian\n",
" Jacobian = ForwardDiff.forwarddiff_jacobian(T!, Float64, fadtype=:dual, n=8, m=8)\n",
" else\n",
" Jacobian = Lin(u -> f_int(X, u, dNdξ, λ, μ))\n",
" end\n",
"\n",
" # Assembly of residual vector R(u) = T(u) - F(u)\n",
" function f(u_)\n",
" u[free_dofs] = u_\n",
" T = f_int(X, u, dNdξ, λ, μ, nlgeom=nlgeom)\n",
" return (T-F)[free_dofs]\n",
" end\n",
"\n",
" # Derivative of residual vector w.r.t displacement field u, \"tangent stiffness\"\n",
" function g(u_)\n",
" u[free_dofs] = u_\n",
" if analytical_jacobian\n",
" J = Jacobian(reshape(u, 8))\n",
" else\n",
" J = Jacobian(u)\n",
" end\n",
" return J[free_dofs, free_dofs]\n",
" end\n",
"\n",
" # Go!\n",
" #sol = nlsolve(not_in_place(f, g), zeros(4), method=:newton) # not working as expected\n",
" #sol = nlsolve(not_in_place(f), zeros(4)) # working, but not quadratic convergence\n",
"\n",
" # working, very well :)\n",
" u2 = zeros(4)\n",
" for i=1:100\n",
" du = g(u2) \\ -f(u2)\n",
" u2 += du\n",
" if norm(du) < 1.0e-6\n",
" #println(\"norm = \",norm(du), \" converged in \", i, \" iterations\")\n",
" break\n",
" end\n",
" end\n",
"\n",
" #u0[free_dofs] = sol.zero\n",
" u0[free_dofs] = u2\n",
" return u0\n",
"end\n",
"\n",
"u0 = zeros(2, 4)\n",
"F = zeros(2, 4)\n",
"F[2, 3] = -20.0\n",
"solve4!(F, u0; nlgeom=true, analytical_jacobian=true)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Solving 1 F = -0.000 0.000 0.000\n",
"Solving 2 F = -0.500 -0.545 -0.549\n",
"Solving 3 F = -1.000 -1.090 -1.103\n",
"Solving 4 F = -1.500 -1.635 -1.662\n",
"Solving 5 F = -2.000 -2.181 -2.222\n",
"Solving 6 F = -2.500 -2.727 -2.783\n",
"Solving 7 F = -3.000 -3.291 -3.343\n",
"Solving 8 F = -3.500 -3.822 -3.899\n",
"Solving 9 F = -4.000 -4.361 -4.449\n",
"Solving 10 F = -4.500 -4.906 -4.992\n",
"Solving 11 F = -5.000 -5.451 -5.526\n",
"Solving 12 F = -5.500 -5.997 -6.047\n",
"Solving 13 F = -6.000 -6.542 -6.554\n",
"Solving 14 F = -6.500 -7.087 -7.045\n",
"Solving 15 F = -7.000 -7.632 -7.516\n",
"Solving 16 F = -7.500 -8.177 -7.965\n",
"Solving 17 F = -8.000 -8.722 -8.388\n",
"Solving 18 F = -8.500 -9.267 -8.783\n",
"Solving 19 F = -9.000 -9.813 -9.145\n",
"Solving 20 F = -9.500 -10.358 -9.474\n",
"Solving 21 F = -10.000 -10.903 -9.767\n",
"Solving 22 F = -10.500 -11.448 -10.025\n",
"Solving 23 F = -11.000 -11.993 -10.250\n",
"Solving 24 F = -11.500 -12.538 -10.443\n",
"Solving 25 F = -12.000 -13.083 -10.610\n",
"Solving 26 F = -12.500 -13.629 -10.753\n",
"Solving 27 F = -13.000 -14.174 -10.878\n",
"Solving 28 F = -13.500 -14.719 -10.987\n",
"Solving 29 F = -14.000 -15.264 -11.084\n",
"Solving 30 F = -14.500 -15.809 -11.171\n",
"Solving 31 F = -15.000 -16.354 -11.250\n",
"Solving 32 F = -15.500 -16.899 -11.323\n",
"Solving 33 F = -16.000 -17.445 -11.390\n",
"Solving 34 F = -16.500 -17.990 -11.452\n",
"Solving 35 F = -17.000 -18.535 -11.512\n",
"Solving 36 F = -17.500 -19.080 -11.568\n",
"Solving 37 F = -18.000 -19.625 -11.621\n",
"Solving 38 F = -18.500 -20.170 -11.673\n",
"Solving 39 F = -19.000 -20.715 -11.723\n",
"Solving 40 F = -19.500 -21.261 -11.771\n",
"Solving 41 F = -20.000 -21.806 -11.817\n"
]
}
],
"source": [
"F = zeros(2, 4)\n",
"N = 41\n",
"F[2, 3] = -20.0\n",
"uall = zeros(2, 4, 2, N)\n",
"u0 = zeros(2, 4)\n",
"t = linspace(0, 1, N)\n",
"for i = 1:N\n",
" @printf \"Solving %2i F = %8.3f \" i t[i]*F[2, 3]\n",
" u0 = solve4!(t[i]*F, u0; nlgeom=false, analytical_jacobian=false)\n",
" uall[:, :, 1, i] = u0\n",
" @printf \"%8.3f\" u0[2, 3]\n",
" u0 = solve4!(t[i]*F, u0; nlgeom=true, analytical_jacobian=true)\n",
" uall[:, :, 2, i] = u0\n",
" @printf \"%8.3f\" u0[2, 3]\n",
" println()\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Loading help data...\n"
]
}
],
"source": [
"import PyPlot"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAroAAAItCAYAAADWntniAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcjXX/x/HXmTtr9qXSjahuUXclUlEpsoyloVIiN0ZptSWUiigR0t0y2hUSI2QP2RnUrbG0MLLMkO3O2MJYZ87vj+s2P2NmmLnmnOt7neu8n4/HPKbOueY679On6/i45nt9Lp/f7/cjIiIiIuIxEaYDiIiIiIgEgxpdEREREfEkNboiIiIi4klqdEVERETEk9ToioiIiIgnqdEVEREREU9SoysiIiIinqRGV0REREQ8SY2uiIiIiHiSGl0RERER8SRXNbqrV6+mS5cu3HDDDRQpUoSrrrqK1q1bs3nz5kzbbty4kcjISIoWLUrp0qVp3749ycnJBlKLiIiIiBv5/H6/33SIs1q1asWqVat4+OGHuemmm9izZw8xMTEcPXqUH374gRtuuAGAnTt3csstt1CyZEm6devGkSNHePvtt6lYsSL/+c9/yJcvn+F3IiIiIiKmuarRXbVqFbVq1eKSSy5Jf2zLli3ceOONtGrViq+++gqAZ599lrFjx5KQkED58uUBWLhwIQ0bNuSTTz6hc+fORvKLiIiIiHu4qtHNTs2aNYmIiGD16tUAXH755dSrV4/Y2NgM21WtWpUKFSowf/58EzFFRERExEVctUY3K36/n//+97+UKVMGgF27drFv3z5uvfXWTNvWqlWLtWvXOh1RRERERFzI9Y3u119/ze7du2ndujUAe/bsAaBcuXKZti1XrhwHDhzg9OnTjmYUEREREfdxdaObkJDAc889R506dejQoQMAx48fB6BAgQKZti9YsGCGbUREREQkfF1y8U3M2Lt3L82aNaNkyZJMnjwZn88HQKFChQA4efJkpp85ceJEhm3OlZyczLx586hUqVKWz4uIiIiIWcePHycpKYnGjRunL1vNC1c2uocPH6ZJkyb89ddfLF++nCuuuCL9ubNLFs4uYTjXnj17KF26dJbjxebNm0e7du2CF1pEREREAmLcuHE89thjed6P6xrdEydOcP/997NlyxYWLFhA1apVMzz/97//nbJly6ZPYDjXf/7zH6pXr57lfitVqgRY/+GqVasW8NziPq1atWLy5MmmY4hDQq3eH38MU6fCvHmmk4SmUKu35E1I1Hv1anjhBShfHt5/HwJwNjIcbdy4kXbt2qX3bXnlqkY3NTWV1q1b8+OPPzJ9+nRuv/32LLd76KGHGDNmDDt37swwR3fz5s288MILWf7M2eUK1apVo0aNGsF5A+Iq+fLlU63DSKjV2+eDihUhhCK7SqjVW/LG9fWeMAG6doV69WDyZCha1HSikBeoZaauanRfeOEFZs6cyf33309ycjLjxo3L8PzZpQcvv/wykyZNol69enTv3p0jR44wfPhwbrrpJqKjo01EFxe67rrrTEcQB4VavffsgXNWZUkuhVq9JW9cW2+/H0aMgN69oX17+Pxz0N1ZXcVVje769evx+XzMnDmTmTNnZnjO5/OlN7rly5dn6dKl9OzZk5deeokCBQrQvHlzRowYodv/ikhI2LsXbrrJdAoRsS01FXr2tJYpvPIKvPGG9asacRVXNbqLFy/O8bbXX389c+fODWIaEZHg2bMHGjc2nUJEbDlxAtq1sxbaf/QRPP206USSDVc1uiKB1Lx5c9MRxEGhVO+0NOuMrpYu2BdK9Za8c1W9DxyAli3hp5+sRjcqynQiuQBX3zBCJC9mzZplOoI4KJTqfeAAnDkDWdzgUXIolOoteeeaem/fDnfdBRs2wMKFanJDgM7oimcNGDDAdARxUCjV++wYcDW69oVSvSXvXFHv9euhSRMoWBBWroQqVUwnkhxQo3uelJQUEhISTMeQAFmzZo3pCAFTtWpVChcubDqGa7l69NB59u61vmvpgn2hVG/JO+P1XrAAHnzQam5nz4bLLzebR3JMje55EhISqFmzpukYIpnEx8eb/7CXgNAZXZEQMm4cREdDgwYwaRIUKWI6keSCGt1s6A5q4hZn7xIj3rFnDxQvbv0GVERcyu+HoUOhb1+r0f3kE83IDUFqdLOhO6iJhJZRo0bx+OOPm46RI3v2wJVXmk4R2kKp3pJ3jtc7NRW6d4eRI6FfPxg4UDNyQ5SmLoiIJ4TSeuw9e7RsIa9Cqd6Sd47W+/hxePhhaz7uJ5/A66+ryQ1hOqMrIp4wcuRI0xFybPduqFTJdIrQFkr1lrxzrN7791sjw9auhenTwU3ze8UWndEVEXGYzuiKuFBSEtx5J/z+OyxZoibXI9Toiog4yO+3zuhqja6Ii6xdC7VrW3dyWbUKbrvNdCIJEDW6IiIO+usvawmgzuiKuMT330PdulC+vHUjiGuvNZ1IAkiNrmSSlJREREQE0dHR6Y917NiRiIgIduzYYTCZSPaiQuRWnJqhGxihUm8JjKDVe+xYaNbManSXLIHLLgvO64gxanQlW75zrjL1+XwZ/l3Ebbp06WI6Qo7s3m1919KFvAmVektgBLzefj8MHgwdOkDHjtaFZ5deGtjXEFdQoys5MmTIEDZu3MiV+tNZXKpRo0amI+SIzugGRqjUWwIjoPVOTYVnn4VXXoEBA+DTT+ESDaHyKlVWcuSKK67giiuuMB1DJOTt3g3FiunkkYgRKSnQpg3Mng2ffw666Yjn6Yyu5EhWa3TPXcublJTEo48+SpkyZShUqBC1atVi9uzZ2e5vwoQJ1KtXjxIlSlCoUCGuv/563nzzTU6dOpVp22nTptGuXTuqVKlCkSJFKFKkCLfeeisffPABfr8/26yJiYl88MEH3HTTTRQuXJh69eoF5j+GSB5otJiIIcnJcN99sGABzJihJjdMqNGVHMtuje727du5/fbb2bFjBx06dKB169b8+uuvtGjRgiVLlmTavlOnTjz22GNs27aNhx9+mC5dulCqVCn69etHZGQkqampGbbv27cv69ato3bt2nTr1o0OHTpw9OhRunfvTocOHbLN2717d/r378/NN99Mjx49uOuuu/L0/sXdpk2bZjpCjuzerUY3EEKl3hIYea73tm1Qpw5s3WpddNa0aUByiftp6YLkWFZnTwGWLFnCwIED6devX/pjbdu2JTIykuHDh3PvvfemPz569GhGjx7Ngw8+yNdff02BAgXSnxs4cCADBw5k5MiRdOvWLf3x7777jsqVK2fKEh0dzdixY+nSpQu3ZTHzcO3ataxbt46rrrrK7luWEDJhwgRatmxpOsZF7d4NFSuaThH6QqXeEhh5qnd8vNXYFitmzci95prAhhNXU6MbKCkpkJDg/OtWrQqFCzv/uueoVKkSr776aobHGjVqRIUKFVi9enWGx9977z3y5cvHF198kaHJBXj11VeJiYnh66+/ztDont/kgnV2uVu3bowdO5bvv/8+y0a3T58+anLDyMSJE01HyJHdu+GOO0ynCH2hUm8JDNv1njsXWrWCf/4TZs6EsmUDG0xcT41uoCQkQM2azr9ufDzUqOH8656jevXqWS5rqFChAj/++GP6v6ekpLB+/XrKli3LO++8k+W+8ufPz8aNGzM8tn//foYPH853333Htm3bSElJyfD8rl27stxXVs2viEm6K5qIg778Ejp3ts7mxsYaPykkZqjRDZSqVa2m08TrGlaiRIksH7/kkktIS0tL//eDBw8CsG/fPl5//fVs93du03zo0CFq1apFUlISt99+Ox07dqRUqVJccsklHDx4kPfee4+TJ09muR9NiRC3OXzYuiva3/9uOomIh/n9MGgQ9O8PTz0FMTEaHxbGVPlAKVzY+JlVtytevDgANWrU4KeffsrRz3z++eckJSUxYMAA+vfvn+G5VatW8d5772X7s7rBhbjN2V8+6IyuSJCcOWPNyP3sM6vZffll0J8FYU1TF8QxRYoU4YYbbuDXX39NP7t7MVu2bAHgoYceyvTc0qVLA5pPQtu5t6x2K90VLXBCod4SODmq97Fj8MAD8MUX1rKFV15RkytqdMVZPXv25NSpU3Tq1InDhw9nev7gwYOsXbs2/d/PXoi
],
"text/plain": [
"Figure(PyObject <matplotlib.figure.Figure object at 0x114260a90>)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"PyObject <matplotlib.text.Text object at 0x1158e5c90>"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PyPlot.plot(-reshape(uall[2, 3, 1, :], N), -t*F[2, 3], color=\"red\", linewidth=1.0, label=\"linear\")\n",
"PyPlot.plot(-reshape(uall[2, 3, 2, :], N), -t*F[2, 3], color=\"blue\", linewidth=1.0, label=\"nonlinear\")\n",
"PyPlot.legend(loc=\"best\")\n",
"PyPlot.grid()\n",
"PyPlot.xlabel(\"Displacement [m]\")\n",
"PyPlot.ylabel(\"Force [N]\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Last remarks. For some unknown reason this notebooks takes occasionally a very long time to complete for 40 different solutions. Sometimes NLsolve takes a very long time when using Newton algorithm. Without Newton all passes well. Actually my implementation of Newton method does not have any problems with speed."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"elapsed time: 57.255576466 seconds\n"
]
},
{
"data": {
"text/plain": [
"57.255576466"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"toc()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.8",
"language": "julia",
"name": "julia-0.3"
},
"language_info": {
"name": "julia",
"version": "0.3.8"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Reference in New Issue Copy Permalink