mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-10-01 22:09:20 +00:00
Jukka Aho
1641 lines
61 KiB
Plaintext
1641 lines
61 KiB
Plaintext
{
|
||
"cells": [
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"# Solving elasticity problems using JuliaFEM\n",
|
||
"\n",
|
||
"Author(s): Jukka Aho\n",
|
||
"\n",
|
||
"**Abstract**: A workflow to solve typical elasticity problem. This document also tries to give some quidelines how to develop JuliaFEM."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Bottom-up design\n",
|
||
"\n",
|
||
"We go piece by piece starting from something simple and going up to more complicated programming model."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"*Design principle 1*: we introduce new ideas using Notebooks."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"*Design principle 2*: we use ``Logging``. Forget ``println``."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 1,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
|
||
]
|
||
},
|
||
"execution_count": 1,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
},
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
|
||
" in int32 at deprecated.jl:49\n",
|
||
" in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
|
||
" in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
|
||
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
|
||
" in anonymous at task.jl:365\n",
|
||
"while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n",
|
||
"WARNING: int32(x) is deprecated, use Int32(x) instead.\n",
|
||
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
|
||
" in int32 at deprecated.jl:49\n",
|
||
" in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
|
||
" in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
|
||
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
|
||
" in anonymous at task.jl:365\n",
|
||
"while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"using Logging\n",
|
||
"using ForwardDiff\n",
|
||
"Logging.configure(level=DEBUG)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"*Design principle 3*: we write docstrings using [numpy style](https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt).\n",
|
||
"\n",
|
||
"*Design principle 4*: we don't use greek characters in code which is implemented to JuliaFEM. In notebooks they are ok.\n",
|
||
"\n",
|
||
"*Design principle 5*: we use 4 space indentation like in Python.\n",
|
||
"\n",
|
||
"First we write some elementary functions to calculate stiffness matrix."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 2,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"calc_local_matrices! (generic function with 1 method)"
|
||
]
|
||
},
|
||
"execution_count": 2,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"\"\"\"\n",
|
||
"Calculate local tangent stiffness matrix and residual force vector\n",
|
||
"R = T - F for elasticity problem.\n",
|
||
"\n",
|
||
"Parameters\n",
|
||
"----------\n",
|
||
"X : Element coordinates\n",
|
||
"u : Displacement field\n",
|
||
"R : Residual force vector\n",
|
||
"K : Tangent stiffness matrix\n",
|
||
"basis : Basis functions\n",
|
||
"dbasis : Derivative of basis functions\n",
|
||
"lambda : Material parameter\n",
|
||
"mu : Material parameter\n",
|
||
"ipoints : integration points\n",
|
||
"iweights : integration weights\n",
|
||
"\n",
|
||
"Returns\n",
|
||
"-------\n",
|
||
"None\n",
|
||
"\n",
|
||
"Notes\n",
|
||
"-----\n",
|
||
"If material parameters are given in list, they are interpolated to gauss\n",
|
||
"points using shape functions.\n",
|
||
"\n",
|
||
"Examples\n",
|
||
"--------\n",
|
||
"\n",
|
||
"\"\"\"\n",
|
||
"function calc_local_matrices2!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)\n",
|
||
" dim, nnodes = size(X)\n",
|
||
" I = eye(dim)\n",
|
||
" R[:,:] = 0.0\n",
|
||
" K[:,:] = 0.0\n",
|
||
"\n",
|
||
" dF = zeros(dim, dim)\n",
|
||
"\n",
|
||
" for m = 1:length(iweights)\n",
|
||
" w = iweights[m]\n",
|
||
" xi = ipoints[m, :]\n",
|
||
" # calculate material parameters\n",
|
||
" lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))\n",
|
||
" mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))\n",
|
||
" Jt = X*dbasis(xi)\n",
|
||
" detJ = det(Jt)\n",
|
||
" dbasisdX = dbasis(xi)*inv(Jt)\n",
|
||
"\n",
|
||
" gradu = u*dbasisdX\n",
|
||
" F = I + gradu # Deformation gradient\n",
|
||
" E = 1/2*(gradu' + gradu + gradu'*gradu) # Green-Lagrange strain tensor\n",
|
||
" S = lambda*trace(E)*I + 2*mu*E # PK2 stress tensor\n",
|
||
" P = F*S # PK1 stress tensor\n",
|
||
"\n",
|
||
" R[:,:] += w*P*dbasisdX'*detJ\n",
|
||
"\n",
|
||
" for p = 1:nnodes\n",
|
||
" for i = 1:dim\n",
|
||
" dF[:,:] = 0.0\n",
|
||
" dF[i,:] = dbasisdX[p,:]\n",
|
||
" dE = 1/2*(F'*dF + dF'*F)\n",
|
||
" dS = lambda*trace(dE)*I + 2*mu*dE\n",
|
||
" dP = dF*S + F*dS\n",
|
||
" for q = 1:nnodes\n",
|
||
" for j = 1:dim\n",
|
||
" K[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*dbasisdX[q,:]')[1]*detJ\n",
|
||
" end\n",
|
||
" end\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" end\n",
|
||
"end\n",
|
||
"\n",
|
||
"\"\"\"\n",
|
||
"Autodiff version.\n",
|
||
"\"\"\"\n",
|
||
"function calc_local_matrices!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)\n",
|
||
" dim, nnodes = size(X)\n",
|
||
" I = eye(dim)\n",
|
||
" R[:,:] = 0.0\n",
|
||
"\n",
|
||
" #dF = zeros(dim, dim)\n",
|
||
"\n",
|
||
" function calc_R!(u, R)\n",
|
||
" for m = 1:length(iweights)\n",
|
||
" w = iweights[m]\n",
|
||
" xi = ipoints[m, :]\n",
|
||
" # calculate material parameters\n",
|
||
" lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))\n",
|
||
" mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))\n",
|
||
" Jt = X*dbasis(xi)\n",
|
||
" detJ = det(Jt)\n",
|
||
" dbasisdX = dbasis(xi)*inv(Jt)\n",
|
||
"\n",
|
||
" gradu = u*dbasisdX\n",
|
||
" F = I + gradu # Deformation gradient\n",
|
||
" E = 1/2*(gradu' + gradu + gradu'*gradu) # Green-Lagrange strain tensor\n",
|
||
" S = lambda*trace(E)*I + 2*mu*E # PK2 stress tensor\n",
|
||
" P = F*S # PK1 stress tensor\n",
|
||
"\n",
|
||
" R[:,:] += w*P*dbasisdX'*detJ\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" # herlper for tangent stiffness matrix\n",
|
||
" function R!(u, R)\n",
|
||
" R[:] = 0\n",
|
||
" calc_R!(reshape(u, dim, nnodes), reshape(R, dim, nnodes))\n",
|
||
" #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))\n",
|
||
" end\n",
|
||
" Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
|
||
"\n",
|
||
" K[:, :] = Jacobian(reshape(u, dim*nnodes))\n",
|
||
" R!(reshape(u, dim*nnodes), reshape(R, dim*nnodes))\n",
|
||
"\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"*Design principle 6*: we test our code. We use FactCheck for testing."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 3,
|
||
"metadata": {
|
||
"collapsed": true
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"using FactCheck"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 4,
|
||
"metadata": {
|
||
"collapsed": false,
|
||
"scrolled": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"test solve one element model\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:26:DEBUG:root:Converged in 6 iterations.\n",
|
||
"13-Aug 00:37:26:DEBUG:root:solution vector: \n",
|
||
" [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n",
|
||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
|
||
"13-Aug 00:37:26:DEBUG:root:norm of u: 3.1292483947150047\n",
|
||
"13-Aug 00:37:27:DEBUG:root:Converged in 6 iterations.\n",
|
||
"13-Aug 00:37:27:DEBUG:root:solution vector: \n",
|
||
" [0.0 0.7433248532717796 1.048521014723486 0.0\n",
|
||
" 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n",
|
||
"13-Aug 00:37:27:DEBUG:root:norm of u: 3.1292483947150056\n",
|
||
"13-Aug 00:37:27:DEBUG:root:Iteration 1\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"Array(Float64,(12,12)) 12x12 Array{Float64,2}"
|
||
]
|
||
},
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:27:DEBUG:root:Norm of du: 0.5992228342549063\n",
|
||
"13-Aug 00:37:27:DEBUG:root:Element displacement: [0.0 -0.02264423092574128 0.022536491965822688 0.0\n",
|
||
" 0.0 -0.12688379170176511 -0.12679760053383046 0.0]\n",
|
||
"13-Aug 00:37:27:DEBUG:root:Element displacement: [-0.02264423092574128 -0.029998807330416523 0.030242156525002267 0.022536491965822688\n",
|
||
" -0.12688379170176511 -0.4041002710283739 -0.40446733271991214 -0.12679760053383046]\n",
|
||
"13-Aug 00:37:27:DEBUG:root:solution vector: \n",
|
||
" [-0.02264423092574128 -0.029998807330416523 0.022536491965822688 0.030242156525002267 0.0 0.0\n",
|
||
" -0.12688379170176511 -0.4041002710283739 -0.12679760053383046 -0.40446733271991214 0.0 0.0]\n",
|
||
"13-Aug 00:37:27:DEBUG:root:norm of u: 0.5992228342549063\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
":\n",
|
||
" 132.8 0.0 23.6 -3.0 -113.6 … 23.6 3.0 -33.2 15.0\n",
|
||
" 0.0 324.8 3.0 77.6 0.0 -3.0 77.6 15.0 -81.2\n",
|
||
" 23.6 3.0 66.4 -15.0 -33.2 0.0 0.0 0.0 0.0\n",
|
||
" -3.0 77.6 -15.0 162.4 15.0 0.0 0.0 0.0 0.0\n",
|
||
" -113.6 0.0 -33.2 15.0 132.8 -33.2 -15.0 23.6 -3.0\n",
|
||
" 0.0 -317.6 15.0 -81.2 0.0 … -15.0 -81.2 3.0 77.6\n",
|
||
" -33.2 -15.0 -56.8 3.0 23.6 0.0 0.0 0.0 0.0\n",
|
||
" -15.0 -81.2 -3.0 -158.8 3.0 0.0 0.0 0.0 0.0\n",
|
||
" 23.6 -3.0 0.0 0.0 -33.2 66.4 15.0 -56.8 3.0\n",
|
||
" 3.0 77.6 0.0 0.0 -15.0 15.0 162.4 -3.0 -158.8\n",
|
||
" -33.2 15.0 0.0 0.0 23.6 … -56.8 -3.0 66.4 -15.0\n",
|
||
" 15.0 -81.2 0.0 0.0 -3.0 3.0 -158.8 -15.0 162.4\n",
|
||
"K norm = 708.0378644377365\n",
|
||
"du = [-0.02264423092574128 -0.029998807330416523 0.022536491965822688 0.030242156525002267 0.0 0.0\n",
|
||
" -0.12688379170176511 -0.4041002710283739 -0.12679760053383046 -0.40446733271991214 0.0 0.0]\n",
|
||
"Out of 3 total facts:"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"delayed_handler (generic function with 4 methods)"
|
||
]
|
||
},
|
||
"execution_count": 4,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"facts(\"test solve one element model\") do\n",
|
||
" X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
|
||
" F = [0 0; 0 0; 0 -2; 0 0]'\n",
|
||
"\n",
|
||
" # Material properties\n",
|
||
" E = 90\n",
|
||
" nu = 0.25\n",
|
||
" mu = E/(2*(1+nu))\n",
|
||
" la = E*nu/((1+nu)*(1-2*nu))\n",
|
||
" la = 2*la*mu/(la + 2*mu)\n",
|
||
"\n",
|
||
" u = zeros(2, 4)\n",
|
||
" du = zeros(2, 4)\n",
|
||
" R = zeros(2, 4)\n",
|
||
" K = zeros(8, 8)\n",
|
||
"\n",
|
||
" basis(xi) = [\n",
|
||
" (1-xi[1])*(1-xi[2])/4\n",
|
||
" (1+xi[1])*(1-xi[2])/4\n",
|
||
" (1+xi[1])*(1+xi[2])/4\n",
|
||
" (1-xi[1])*(1+xi[2])/4]\n",
|
||
"\n",
|
||
" dbasis(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n",
|
||
" (1-xi[2])/4.0 -(1+xi[1])/4.0\n",
|
||
" (1+xi[2])/4.0 (1+xi[1])/4.0\n",
|
||
" -(1+xi[2])/4.0 (1-xi[1])/4.0]\n",
|
||
"\n",
|
||
" ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
|
||
" iweights = [1, 1, 1, 1]\n",
|
||
" free_dofs = [3, 4, 5, 6]\n",
|
||
"\n",
|
||
" for i=1:10\n",
|
||
" calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
|
||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
|
||
" u += du\n",
|
||
" if norm(du) < 1.0e-9\n",
|
||
" Logging.debug(\"Converged in $i iterations.\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Tested against Elmer solution\n",
|
||
" Logging.debug(\"solution vector: \\n $u\")\n",
|
||
" @fact u[2, 3] --> roughly(-2.222244754401764)\n",
|
||
" norm1 = norm(u)\n",
|
||
" Logging.debug(\"norm of u: $(norm(u))\")\n",
|
||
"\n",
|
||
" # We rotate model a bit and make sure that L2 norm is same\n",
|
||
" phi = 30/180*pi\n",
|
||
" rmat = [\n",
|
||
" cos(phi) -sin(phi)\n",
|
||
" sin(phi) cos(phi)]\n",
|
||
" X = rmat*X\n",
|
||
" F = rmat*F\n",
|
||
" u = zeros(2, 4)\n",
|
||
" for i=1:10\n",
|
||
" calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
|
||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
|
||
" u += du\n",
|
||
" if norm(du) < 1.0e-9\n",
|
||
" Logging.debug(\"Converged in $i iterations.\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
" Logging.debug(\"solution vector: \\n $u\")\n",
|
||
" Logging.debug(\"norm of u: $(norm(u))\")\n",
|
||
" @fact norm(u) --> roughly(norm1) \n",
|
||
"\n",
|
||
" # test two element model\n",
|
||
" X = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'\n",
|
||
" u = zeros(2, 6)\n",
|
||
" du = zeros(2, 6)\n",
|
||
" R = zeros(2, 4)\n",
|
||
" K = zeros(8, 8)\n",
|
||
" ass1 = [9, 10, 1, 2, 5, 6, 11, 12]\n",
|
||
" ass2 = [1, 2, 3, 4, 7, 8, 5, 6]\n",
|
||
" free_dofs = collect(1:8)\n",
|
||
" F = [0 0; 0 0; 0 0; 0 -0.1; 0 0; 0 0]'\n",
|
||
"\n",
|
||
" A = zeros(12, 12)\n",
|
||
" b = zeros(2, 6)\n",
|
||
" for i=1:1\n",
|
||
" Logging.debug(\"Iteration $i\")\n",
|
||
" A[:,:] = 0.0\n",
|
||
" b[:] = 0.0\n",
|
||
" #Logging.debug(\"Assembling\")\n",
|
||
" for ass in (ass1, ass2)\n",
|
||
" #Logging.debug(\"ass = $ass, u[ass] = $(u[ass])\")\n",
|
||
" calc_local_matrices!(X, u[ass], R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
|
||
" A[ass,ass] += K\n",
|
||
" b[ass] += R[:]\n",
|
||
" end\n",
|
||
" dump(round(A, 2))\n",
|
||
" println(\"K norm = $(norm(A[free_dofs, free_dofs]))\")\n",
|
||
" du[free_dofs] = A[free_dofs, free_dofs] \\ -(b - F)[free_dofs]\n",
|
||
" println(\"du = $du\")\n",
|
||
" u += du\n",
|
||
" Logging.debug(\"Norm of du: $(norm(du))\")\n",
|
||
" for ass in (ass1, ass2)\n",
|
||
" Logging.debug(\"Element displacement: $(reshape(u[ass], 2, 4))\")\n",
|
||
" end\n",
|
||
" if norm(du) < 1.0e-9\n",
|
||
" Logging.debug(\"Converged in $i iterations.\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
" Logging.debug(\"solution vector: \\n $u\")\n",
|
||
" Logging.debug(\"norm of u: $(norm(u))\")\n",
|
||
" @pending norm(u) --> :something\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"One element solutions are not particularly interesting so next step is to create function that assembles global matrix from local matrices. Some data types:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 5,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"#type Node\n",
|
||
"# id :: Int\n",
|
||
"# #elements :: Array{Int64, 1}\n",
|
||
"#end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 6,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"type Element\n",
|
||
" id :: Int\n",
|
||
" node_ids :: Array{Int64, 1}\n",
|
||
" coordinates :: Array{Float64, 2}\n",
|
||
" attributes :: Dict{ASCIIString, Any}\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 7,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"\n",
|
||
" Verified: 2\n",
|
||
" Pending: 1\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"type Assembly\n",
|
||
" # LHS\n",
|
||
" I :: Array{Int64, 1}\n",
|
||
" J :: Array{Int64, 1}\n",
|
||
" A :: Array{Float64, 1}\n",
|
||
" # RHS\n",
|
||
" i :: Array{Int64, 1}\n",
|
||
" b :: Array{Float64, 1}\n",
|
||
" # global dofs for each element\n",
|
||
" gdofs :: Dict{Int64, Array{Int64, 1}}\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 8,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"get_integration_scheme (generic function with 2 methods)"
|
||
]
|
||
},
|
||
"execution_count": 8,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"\"\"\"\n",
|
||
"Return shape functions and their derivatives for a element.\n",
|
||
"\n",
|
||
"Parameters\n",
|
||
"----------\n",
|
||
"element::Element\n",
|
||
"\n",
|
||
"Returns\n",
|
||
"-------\n",
|
||
"tuple (basis, dbasis)\n",
|
||
"\"\"\"\n",
|
||
"function get_shape_functions(el::Element)\n",
|
||
" ndim, nnodes = size(el.coordinates)\n",
|
||
" #Logging.debug(\"ndim = $ndim, nnodes=$nnodes\")\n",
|
||
" if (nnodes == 4) & (ndim == 2)\n",
|
||
" basis(xi) = [\n",
|
||
" (1-xi[1])*(1-xi[2])/4\n",
|
||
" (1+xi[1])*(1-xi[2])/4\n",
|
||
" (1+xi[1])*(1+xi[2])/4\n",
|
||
" (1-xi[1])*(1+xi[2])/4]\n",
|
||
" dbasis(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n",
|
||
" (1-xi[2])/4.0 -(1+xi[1])/4.0\n",
|
||
" (1+xi[2])/4.0 (1+xi[1])/4.0\n",
|
||
" -(1+xi[2])/4.0 (1-xi[1])/4.0]\n",
|
||
" return basis, dbasis\n",
|
||
" elseif (nnodes == 10) & (ndim == 3)\n",
|
||
" basis(xi) = [\n",
|
||
" (xi[1] + xi[2] + xi[3] - 1)*(2*xi[1] + 2*xi[2] + 2*xi[3] - 1)\n",
|
||
" xi[1]*(2*xi[1] - 1)\n",
|
||
" xi[2]*(2*xi[2] - 1)\n",
|
||
" xi[3]*(2*xi[3] - 1)\n",
|
||
" -4*xi[1]*(xi[1] + xi[2] + xi[3] - 1)\n",
|
||
" 4*xi[1]*xi[2]\n",
|
||
" -4*xi[2]*(xi[1] + xi[2] + xi[3] - 1)\n",
|
||
" -4*xi[3]*(xi[1] + xi[2] + xi[3] - 1)\n",
|
||
" 4*xi[1]*xi[3]\n",
|
||
" 4*xi[2]*xi[3]\n",
|
||
" ]\n",
|
||
" dbasis(xi) = [\n",
|
||
" 4*xi[1] + 4*xi[2] + 4*xi[3] - 3 4*xi[1] + 4*xi[2] + 4*xi[3] - 3 4*xi[1] + 4*xi[2] + 4*xi[3] - 3\n",
|
||
" 4*xi[1] - 1 0 0\n",
|
||
" 0 4*xi[2] - 1 0\n",
|
||
" 0 0 4*xi[3] - 1\n",
|
||
" -4*(2*xi[1] + xi[2] + xi[3] - 1) -4*xi[1] -4*xi[1]\n",
|
||
" 4*xi[2] 4*xi[1] 0\n",
|
||
" -4*xi[2] -4*(xi[1] + 2*xi[2] + xi[3] - 1) -4*xi[2]\n",
|
||
" -4*xi[3] -4*xi[3] -4*(xi[1] + xi[2] + 2*xi[3] - 1)\n",
|
||
" 4*xi[3] 0 4*xi[1]\n",
|
||
" 0 4*xi[3] 4*xi[2]\n",
|
||
" ]\n",
|
||
" return basis, dbasis\n",
|
||
" end\n",
|
||
" throw(\"Unknown function space, ndim=$ndim, nnodes=$nnodes\")\n",
|
||
"end\n",
|
||
"\n",
|
||
"\"\"\"\n",
|
||
"\"\"\"\n",
|
||
"function get_integration_scheme(el::Element, order=2)\n",
|
||
" ndim, nnodes = size(el.coordinates)\n",
|
||
" if (nnodes == 4) & (order == 2) & (ndim == 2)\n",
|
||
" ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
|
||
" iweights = [1, 1, 1, 1]\n",
|
||
" return ipoints, iweights\n",
|
||
" elseif (nnodes == 10) & (order == 2) & (ndim == 3) # c3d10\n",
|
||
" # from code aster documentation\n",
|
||
" a = 1/20*(5-sqrt(5))\n",
|
||
" b = 1/20*(5+3*sqrt(5))\n",
|
||
" ipoints = [a a a; a a b; a b a; b a a]\n",
|
||
" iweights = 1/24*[1 1 1 1]\n",
|
||
" return ipoints, iweights\n",
|
||
" end\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 9,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"assemble_element! (generic function with 2 methods)"
|
||
]
|
||
},
|
||
"execution_count": 9,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"function assemble_element!(ass::Assembly, el::Element, io=2)\n",
|
||
"\n",
|
||
" # Material properties\n",
|
||
" E = el.attributes[\"Young\"]\n",
|
||
" nu = el.attributes[\"Poisson\"]\n",
|
||
" mu = E/(2*(1+nu))\n",
|
||
" la = E*nu/((1+nu)*(1-2*nu))\n",
|
||
" la = 2*la*mu/(la + 2*mu)\n",
|
||
"\n",
|
||
" dofs = prod(size(el.coordinates))\n",
|
||
" X = el.coordinates\n",
|
||
" u = el.attributes[\"displacement\"]\n",
|
||
" R = el.attributes[\"displacement nodal force\"]\n",
|
||
" K = el.attributes[\"displacement tangent stiffness\"]\n",
|
||
"\n",
|
||
" gdofs = ass.gdofs[el.id]\n",
|
||
" #Logging.debug(\"Assemble element $(el.id) to gdofs $gdofs\")\n",
|
||
" basis, dbasis = get_shape_functions(el)\n",
|
||
" ipoints, iweights = get_integration_scheme(el, io)\n",
|
||
" calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
|
||
"\n",
|
||
" for i=1:dofs\n",
|
||
" for j=1:dofs\n",
|
||
" push!(ass.I, gdofs[i])\n",
|
||
" push!(ass.J, gdofs[j])\n",
|
||
" push!(ass.A, K[i,j])\n",
|
||
" end\n",
|
||
" push!(ass.i, gdofs[i])\n",
|
||
" push!(ass.b, R[i])\n",
|
||
" end\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Time to test again. From last test we know that correct solution is\n",
|
||
"\n",
|
||
" [0.0 -0.39914506095474317 -0.07228582695592449 0.0\n",
|
||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
|
||
"\n",
|
||
"This time we assemble global stiffness matrix in different order, 2 3 4 1"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 10,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"WARNING: the `=>` syntax is deprecated, use `-->` instead\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"one element assembly\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:30:DEBUG:root:Adding nodes to array\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Creating elements\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 1\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 3.090022136728999\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 2\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 0.32121316021535135\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 3\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 0.040431781939994194\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 4\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 0.0009291101052124042\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 5\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 1.5638899022213175e-7\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Starting iteration 6\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Solution norm = 1.0118539067290854e-14\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Converged in 6 iterations.\n",
|
||
"13-Aug 00:37:30:DEBUG:root:Displacement of element = \n",
|
||
"[-0.39914506095474334 -0.0722858269559246 0.0 0.0\n",
|
||
" -2.1779892317073504 -2.2222447544017645 0.0 0.0]\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"1 fact verified.\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"delayed_handler (generic function with 4 methods)"
|
||
]
|
||
},
|
||
"execution_count": 10,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"facts(\"one element assembly\") do\n",
|
||
" # Create model\n",
|
||
" #Logging.debug(\"Creating nodes\")\n",
|
||
" #n1 = Node(1)\n",
|
||
" #n2 = Node(2)\n",
|
||
" #n3 = Node(3)\n",
|
||
" #n4 = Node(4)\n",
|
||
" Logging.debug(\"Adding nodes to array\")\n",
|
||
" #nodes = [n1.id, n2.id, n3.id, n4.id]\n",
|
||
" node_ids = [1, 2, 3, 4]\n",
|
||
" coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n",
|
||
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||
" Logging.debug(\"Creating elements\")\n",
|
||
" el = Element(1, node_ids, coordinates, attributes)\n",
|
||
"\n",
|
||
" # Initialize elements ready for solution\n",
|
||
" el.attributes[\"displacement\"] = zeros(2, 4)\n",
|
||
" el.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
|
||
" el.attributes[\"displacement tangent stiffness\"] = zeros(8, 8)\n",
|
||
"\n",
|
||
" for i=1:10\n",
|
||
" Logging.debug(\"Starting iteration $i\")\n",
|
||
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())\n",
|
||
" ass.gdofs[el.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n",
|
||
" Logging.debug(\"Assembling\")\n",
|
||
" assemble_element!(ass, el)\n",
|
||
"\n",
|
||
" # Boundary conditions\n",
|
||
" F = [0 0; 0 -2; 0 0; 0 0]'\n",
|
||
" F = reshape(F, prod(size(F)))\n",
|
||
" free_dofs = [1, 2, 3, 4]\n",
|
||
"\n",
|
||
" # solution\n",
|
||
" K = sparse(ass.I, ass.J, ass.A)\n",
|
||
" R = full(sparsevec(ass.i, ass.b))\n",
|
||
" R = R - F\n",
|
||
" du = zeros(8) # must be determined from ass\n",
|
||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n",
|
||
"\n",
|
||
" Logging.debug(\"Solution norm = $(norm(du))\")\n",
|
||
"\n",
|
||
" # update solution back to elements\n",
|
||
" eldu = du[ass.gdofs[el.id]]\n",
|
||
" eldu = reshape(eldu, (2, round(Int, length(eldu)/2)))\n",
|
||
" el.attributes[\"displacement\"] += eldu\n",
|
||
" if norm(du) < 1.0e-9\n",
|
||
" Logging.debug(\"Converged in $i iterations.\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
" disp = el.attributes[\"displacement\"]\n",
|
||
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
|
||
" @fact norm(disp) => roughly(3.1292483947150043)\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Seems to be working. But we still need to handle boundary conditions more \"cleverly\" and generalize assembly to several elements (which is not problem).\n",
|
||
"\n",
|
||
"First of all, essential boundary conditions are nothing more than equality constraints saying that value for some degree of freedom is fixed. Elimination is just a special case when this value equals to zero. There is couple of different strategies to handle essential boundary conditions. One option is to force them using Lagrange multipliers which can also be used to create all kind of kinematic constraints also. (For example, contact can be considered as a kinematic constraint.) Another option is to manipulate matrix such a way that constraint is satisfied.\n",
|
||
"\n",
|
||
"Because we are now going \"bottom-up\", we develop something extremely simple that however deals with the problem:"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 11,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [],
|
||
"source": [
|
||
"type BC\n",
|
||
" dofs :: Array{Int64, 1}\n",
|
||
" values :: Array{Float64, 1}\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 12,
|
||
"metadata": {
|
||
"collapsed": false,
|
||
"scrolled": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"solve two element problem\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:31:DEBUG:root:Creating elements\n",
|
||
"13-Aug 00:37:31:INFO:root:create_ldof2gdofmap: dofs per node: 2\n",
|
||
"13-Aug 00:37:31:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],5=>[9,10],6=>[11,12],1=>[1,2])\n",
|
||
"13-Aug 00:37:32:INFO:root:solve!: dofs per node: 2\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Problem size = 12\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Starting iteration 1\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Added 5 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Solving system of equations. Total size = 16\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Solution norm du = 0.6015838690633517\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Starting iteration 2\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Added 5 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Solving system of equations. Total size = 16\n",
|
||
"13-Aug 00:37:32:DEBUG:root:Solution norm du = 0.013420417380980414\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Starting iteration 3\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solution norm du = 0.00032202957936854873\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Starting iteration 4\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solution norm du = 9.906677094801476e-8\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Starting iteration 5\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Solution norm du = 5.027624635820698e-15\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Converged in 5 iterations.\n",
|
||
"13-Aug 00:37:33:DEBUG:root:Displacement of element = \n",
|
||
"[-0.02442313597467864 -0.039356000063335075 0.021097993207232584 0.020877031423993653\n",
|
||
" -0.12673626841485705 -0.40433021969759375 -0.40656320177872923 -0.1275776940913048]\n"
|
||
]
|
||
},
|
||
{
|
||
"name": "stdout",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"0 facts verified.\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"delayed_handler (generic function with 4 methods)"
|
||
]
|
||
},
|
||
"execution_count": 12,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"\"\"\"\n",
|
||
"Create local dof to global dof mapping for given elements\n",
|
||
"\"\"\"\n",
|
||
"function create_ldof2gdofmap(elements; ndofs=2)\n",
|
||
" Logging.info(\"create_ldof2gdofmap: dofs per node: $ndofs\")\n",
|
||
"\n",
|
||
" all_node_ids = Int64[]\n",
|
||
" for el in elements\n",
|
||
" eldim, elnodes = size(el.coordinates)\n",
|
||
" for nid in el.node_ids\n",
|
||
" push!(all_node_ids, nid)\n",
|
||
" end\n",
|
||
" end\n",
|
||
" all_node_ids = unique(all_node_ids)\n",
|
||
" sort!(all_node_ids)\n",
|
||
"\n",
|
||
" # Assign global dof for each node\n",
|
||
" pdim = 1\n",
|
||
" ngdofs = Dict{Int64, Array{Int64,1}}()\n",
|
||
" for nid in all_node_ids\n",
|
||
" ngdofs[nid] = collect(pdim:pdim+ndofs-1)\n",
|
||
" pdim += ndofs\n",
|
||
" end\n",
|
||
"\n",
|
||
" return ngdofs\n",
|
||
"end\n",
|
||
"\n",
|
||
"function solve!(elements, dofmap, neumann_bcs, dirichlet_bcs; ndofs=2, max_iterations=10)\n",
|
||
"\n",
|
||
" dbc = \"lagrange\"\n",
|
||
" \n",
|
||
" Logging.info(\"solve!: dofs per node: $ndofs\")\n",
|
||
" pdim = length(dofmap)*ndofs\n",
|
||
" Logging.debug(\"Problem size = $pdim\")\n",
|
||
"\n",
|
||
" # Initialize elements ready for solution\n",
|
||
" for el in elements\n",
|
||
" eldim, elnodes = size(el.coordinates)\n",
|
||
" el.attributes[\"displacement\"] = zeros(ndofs, elnodes)\n",
|
||
" el.attributes[\"displacement nodal force\"] = zeros(ndofs, elnodes)\n",
|
||
" el.attributes[\"displacement tangent stiffness\"] = zeros(ndofs*elnodes, ndofs*elnodes)\n",
|
||
" end\n",
|
||
"\n",
|
||
" # Assign global dofs for elements\n",
|
||
" gdofs = Dict{Int64, Array{Int64,1}}()\n",
|
||
" for el in elements\n",
|
||
" gdofs[el.id] = Int64[]\n",
|
||
" for nid in el.node_ids\n",
|
||
" for ndof in dofmap[nid]\n",
|
||
" push!(gdofs[el.id], ndof)\n",
|
||
" end\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n",
|
||
"\n",
|
||
" for iter=1:max_iterations\n",
|
||
" Logging.debug(\"Starting iteration $iter\")\n",
|
||
" ass.I = []\n",
|
||
" ass.J = []\n",
|
||
" ass.A = []\n",
|
||
" ass.i = []\n",
|
||
" ass.b = []\n",
|
||
" \n",
|
||
" Logging.debug(\"Assembling\")\n",
|
||
" for el in elements\n",
|
||
" assemble_element!(ass, el)\n",
|
||
" #println(\"Element stiffness matrix\")\n",
|
||
" #dump(round(el.attributes[\"displacement tangent stiffness\"], 2))\n",
|
||
" end\n",
|
||
"\n",
|
||
" i = 1\n",
|
||
" if dbc == \"lagrange\"\n",
|
||
" Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n",
|
||
" # Dirichlet boundary conditions\n",
|
||
" for bc in dirichlet_bcs\n",
|
||
" for (dof, val) in zip(bc.dofs, bc.values)\n",
|
||
" #Logging.debug(\"dof $dof => $val\")\n",
|
||
" push!(ass.I, dof)\n",
|
||
" push!(ass.J, pdim+i)\n",
|
||
" push!(ass.A, 1)\n",
|
||
" push!(ass.I, pdim+i)\n",
|
||
" push!(ass.J, dof)\n",
|
||
" push!(ass.A, 1)\n",
|
||
" push!(ass.i, pdim+i)\n",
|
||
" push!(ass.b, 0)\n",
|
||
" i += 1\n",
|
||
" end\n",
|
||
" end\n",
|
||
" Logging.debug(\"Added $i Lagrange multipliers to model\")\n",
|
||
" end\n",
|
||
" i -= 1\n",
|
||
"\n",
|
||
" Logging.debug(\"Adding Neumann boundary conditions\")\n",
|
||
" F = zeros(pdim+i)\n",
|
||
" # Neumann boundary conditions\n",
|
||
" for bc in neumann_bcs\n",
|
||
" for (dof, val) in zip(bc.dofs, bc.values)\n",
|
||
" F[dof] += val\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" Logging.debug(\"Solving system of equations. Total size = $(pdim+i)\")\n",
|
||
" # solution\n",
|
||
" K = sparse(ass.I, ass.J, ass.A)\n",
|
||
" R = full(sparsevec(ass.i, ass.b))\n",
|
||
" R = R - F\n",
|
||
" #Logging.debug(dump(round(full(K), 2)))\n",
|
||
" #print_matrix(full(K))\n",
|
||
" #Logging.debug(dump(round(R', 1)))\n",
|
||
"\n",
|
||
" if dbc == \"eliminate\"\n",
|
||
" Logging.debug(\"Eliminating Dirichlet boundary conditions\")\n",
|
||
" throw(\"Implement this properly\")\n",
|
||
" free_dofs = [1, 2, 3, 4, 5, 6, 7, 8]\n",
|
||
" du = zeros(12)\n",
|
||
" Logging.debug(\"K norm = $(norm(full(K[free_dofs, free_dofs])))\")\n",
|
||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n",
|
||
" else\n",
|
||
" du = K \\ -R\n",
|
||
" end\n",
|
||
"\n",
|
||
" #du = reshape(du, 2, 6)\n",
|
||
" solnorm = norm(du[1:pdim])\n",
|
||
" #Logging.debug(\"du = $du\")\n",
|
||
" Logging.debug(\"Solution norm du = $solnorm\")\n",
|
||
"\n",
|
||
" # update solution back to elements\n",
|
||
" for el in elements\n",
|
||
" #Logging.debug(\"update element $(el.id)\")\n",
|
||
" eldisp = el.attributes[\"displacement\"]\n",
|
||
" #Logging.debug(\"displacement before update $(el.id) : \\n$eldisp\")\n",
|
||
" eldu = du[ass.gdofs[el.id]]\n",
|
||
" eldu = reshape(eldu, (ndofs, round(Int, length(eldu)/ndofs)))\n",
|
||
" #Logging.debug(\"eldu for element $(el.id) \\n$eldu\")\n",
|
||
" el.attributes[\"displacement\"] += eldu\n",
|
||
" eldisp = el.attributes[\"displacement\"]\n",
|
||
" #Logging.debug(\"displacement after update $(el.id) : \\n$eldisp\")\n",
|
||
"\n",
|
||
" end\n",
|
||
" if solnorm < 1.0e-9\n",
|
||
" Logging.debug(\"Converged in $iter iterations.\")\n",
|
||
" break\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
"end\n",
|
||
"\n",
|
||
"ENV[\"COLUMNS\"] = 160\n",
|
||
"\n",
|
||
"function test1():\n",
|
||
" facts(\"solve one element problem\") do\n",
|
||
" # Create model\n",
|
||
" Logging.debug(\"Creating nodes\")\n",
|
||
" node_ids = [1, 2, 3, 4]\n",
|
||
" coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n",
|
||
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||
" Logging.debug(\"Creating elements\")\n",
|
||
" el = Element(1, node_ids, coordinates, attributes)\n",
|
||
" elements = [el]\n",
|
||
" dofmap = create_ldof2gdofmap(elements)\n",
|
||
" Logging.debug(dofmap)\n",
|
||
" # Boundary conditions\n",
|
||
" # here we want to create nodal force for third dof, that is, node id 2, second dof\n",
|
||
" bc1 = BC([dofmap[2][2]], [-2.0])\n",
|
||
" # dirichlet bc, set dx=dy=0 on support\n",
|
||
" bc2 = BC([dofmap[3][1], dofmap[3][2], dofmap[4][1], dofmap[4][2]], [0.0, 0.0, 0.0, 0.0])\n",
|
||
" solve!(elements, dofmap, [bc1], [bc2]; max_iterations=7)\n",
|
||
" disp = elements[1].attributes[\"displacement\"]\n",
|
||
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
|
||
" @fact norm(disp) => roughly(3.1292483947150043)\n",
|
||
" end\n",
|
||
"end\n",
|
||
"\n",
|
||
"facts(\"solve two element problem\") do\n",
|
||
" # Create model\n",
|
||
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||
"\n",
|
||
" Logging.debug(\"Creating elements\")\n",
|
||
" nids1 = [5, 1, 3, 6]\n",
|
||
" coords1 = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'\n",
|
||
" el1 = Element(1, nids1, coords1, copy(attributes))\n",
|
||
" nids2 = [1, 2, 4, 3]\n",
|
||
" coords2 = [5.0 0.0; 10.0 0.0; 10.0 1.0; 5.0 1.0]'\n",
|
||
" el2 = Element(2, nids2, coords2, copy(attributes))\n",
|
||
" elements = [el1, el2]\n",
|
||
"\n",
|
||
" dofmap = create_ldof2gdofmap(elements)\n",
|
||
" Logging.debug(dofmap)\n",
|
||
" # Boundary conditions\n",
|
||
" # here we want to create nodal force for third dof, that is, node id 2, second dof\n",
|
||
" bc1 = BC([dofmap[4][2]], [-0.1])\n",
|
||
" # dirichlet bc, set dx=dy=0 on support\n",
|
||
" bc2 = BC([dofmap[5][1], dofmap[5][2], dofmap[6][1], dofmap[6][2]], [0.0, 0.0, 0.0, 0.0])\n",
|
||
" solve!(elements, dofmap, [bc1], [bc2]; max_iterations=7)\n",
|
||
" disp = elements[2].attributes[\"displacement\"]\n",
|
||
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
|
||
" #@fact norm(disp) => roughly(3.1292483947150043)\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Top-down design\n",
|
||
"\n",
|
||
"Next we see this problem from \"other direction\", by parsing ABAQUS .inp file and making 3d simulation."
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 13,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"\n",
|
||
"WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:11.\n",
|
||
"Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
|
||
"\n",
|
||
"WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:32.\n",
|
||
"Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
|
||
"13-Aug 00:37:37:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n",
|
||
"WARNING: beginswith is deprecated, use startswith instead.\n",
|
||
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
|
||
" in beginswith at deprecated.jl:30\n",
|
||
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n",
|
||
" in include_string at loading.jl:99\n",
|
||
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
|
||
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
|
||
" in anonymous at task.jl:365\n",
|
||
"while loading In[13], in expression starting on line 3\n",
|
||
"WARNING: beginswith is deprecated, use startswith instead.\n",
|
||
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
|
||
" in beginswith at deprecated.jl:30\n",
|
||
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n",
|
||
" in include_string at loading.jl:99\n",
|
||
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
|
||
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
|
||
" in anonymous at task.jl:365\n",
|
||
"while loading In[13], in expression starting on line 3\n",
|
||
"13-Aug 00:37:38:DEBUG:root:Found NODE section\n",
|
||
"13-Aug 00:37:38:DEBUG:root:Found ELEMENT section\n",
|
||
"WARNING: integer(s::AbstractString) is deprecated, use parse(Int,s) instead.\n",
|
||
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
|
||
" in integer at deprecated.jl:49\n",
|
||
" in map at abstractarray.jl:1251\n",
|
||
" in parse_element_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:59\n",
|
||
" in process_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:108\n",
|
||
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:117\n",
|
||
" in include_string at loading.jl:99\n",
|
||
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
|
||
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
|
||
" in anonymous at task.jl:365\n",
|
||
"while loading In[13], in expression starting on line 3\n",
|
||
"13-Aug 00:37:39:DEBUG:root:120 elements found\n",
|
||
"13-Aug 00:37:40:INFO:root:Creating ELSET Body1\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Creating node set SUPPORT\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Creating node set LOAD\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
|
||
"13-Aug 00:37:40:DEBUG:root:Creating node set TOP\n"
|
||
]
|
||
},
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"Dict{Any,Any} with 4 entries:\n",
|
||
" \"nodes\" => Dict{Any,Any}(288=>[97.5,7.5,10.0],11=>[92.5,2.5,5.0],134=>[45.0,10.0,0.0],158=>[2.5,2.5,0.0],160=>[7.5,7.5,0.0],215=>[60.0,0.0,5.0],29=>[2.5,7…\n",
|
||
" \"elements\" => Dict{Any,Any}(68=>[71,144,149,198,51,150,57,43,50,214],2=>[204,199,175,130,207,208,209,3,4,176],89=>[95,78,104,52,127,126,106,60,68,67],11=>[15…\n",
|
||
" \"elsets\" => Dict{Any,Any}(\"Body1\"=>[1,2,3,4,5,6,7,8,9,10 … 111,112,113,114,115,116,117,118,119,120])\n",
|
||
" \"nsets\" => Dict{Any,Any}(\"LOAD\"=>[82,84,87,179,197,246,249,256,257],\"SUPPORT\"=>[108,109,111,155,162,216,225,281,298],\"TOP\"=>[70,75,76,84,88,90,95,96,98,10…"
|
||
]
|
||
},
|
||
"execution_count": 13,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"using JuliaFEM.abaqus_reader\n",
|
||
"fid = open(\"../geometry/3d_beam/palkki.inp\")\n",
|
||
"model = JuliaFEM.abaqus_reader.parse_abaqus(fid)\n",
|
||
"close(fid)\n",
|
||
"model"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 14,
|
||
"metadata": {
|
||
"collapsed": false,
|
||
"scrolled": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:41:DEBUG:root:Creating elements\n",
|
||
"13-Aug 00:37:41:INFO:root:create_ldof2gdofmap: dofs per node: 3\n",
|
||
"13-Aug 00:37:41:INFO:root:solve!: dofs per node: 3\n",
|
||
"13-Aug 00:37:41:DEBUG:root:Problem size = 894\n",
|
||
"13-Aug 00:37:41:DEBUG:root:Starting iteration 1\n",
|
||
"13-Aug 00:37:41:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Added 28 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Solving system of equations. Total size = 921\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Solution norm du = 54.19642700242575\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Starting iteration 2\n",
|
||
"13-Aug 00:37:43:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Added 28 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Solving system of equations. Total size = 921\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Solution norm du = 1.68729144400063\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Starting iteration 3\n",
|
||
"13-Aug 00:37:45:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Added 28 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Solving system of equations. Total size = 921\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Solution norm du = 0.04175080278208098\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Starting iteration 4\n",
|
||
"13-Aug 00:37:47:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:48:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:48:DEBUG:root:Added 28 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:48:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:48:DEBUG:root:Solving system of equations. Total size = 921\n",
|
||
"13-Aug 00:37:49:DEBUG:root:Solution norm du = 2.4176157836844963e-5\n",
|
||
"13-Aug 00:37:49:DEBUG:root:Starting iteration 5\n",
|
||
"13-Aug 00:37:49:DEBUG:root:Assembling\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Added 28 Lagrange multipliers to model\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Adding Neumann boundary conditions\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Solving system of equations. Total size = 921\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Solution norm du = 1.4448231751500831e-11\n",
|
||
"13-Aug 00:37:50:DEBUG:root:Converged in 5 iterations.\n",
|
||
"13-Aug 00:37:50:INFO:root:Maximum absolute displacement in y direction: 5.245400568184194\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"function solve_3d_model()\n",
|
||
" Logging.debug(\"Creating elements\")\n",
|
||
" elements = Element[]\n",
|
||
" coordinates = zeros(3, 10)\n",
|
||
" for (elid, node_ids) in model[\"elements\"]\n",
|
||
" for (i, nid) in enumerate(node_ids)\n",
|
||
" coordinates[:,i] = model[\"nodes\"][nid]\n",
|
||
" end\n",
|
||
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||
" el = Element(elid, node_ids, coordinates, attributes)\n",
|
||
" push!(elements, el)\n",
|
||
" end\n",
|
||
"\n",
|
||
" # create \"dofmap\" so that we know how to assemble global stiffness matrix\n",
|
||
" dofmap = create_ldof2gdofmap(elements; ndofs=3)\n",
|
||
"\n",
|
||
" # Boundary conditions\n",
|
||
"\n",
|
||
" # dirichlet bc, set dx=dy=dz for all nodes in set SUPPORT\n",
|
||
" bc_support = BC(Int64[], Float64[])\n",
|
||
" for nid in model[\"nsets\"][\"SUPPORT\"]\n",
|
||
" for i=1:3\n",
|
||
" push!(bc_support.dofs, dofmap[nid][i])\n",
|
||
" push!(bc_support.values, 0.0)\n",
|
||
" end\n",
|
||
" end\n",
|
||
"\n",
|
||
" # force boundary condition, put -1 to 2nd dof for each node in set LOAD\n",
|
||
" bc_load = BC(Int64[], Float64[])\n",
|
||
" for nid in model[\"nsets\"][\"LOAD\"]\n",
|
||
" push!(bc_load.dofs, dofmap[nid][2])\n",
|
||
" push!(bc_load.values, -30.0)\n",
|
||
" end\n",
|
||
"\n",
|
||
" #solve!(elements, [bc1], [bc2]; max_iterations=7)\n",
|
||
" neumann_bcs = [bc_load]\n",
|
||
" dirichlet_bcs = [bc_support]\n",
|
||
" # ndofs = dimension of unknown field in nodes\n",
|
||
" solve!(elements, dofmap, neumann_bcs, dirichlet_bcs; ndofs=3, max_iterations=10)\n",
|
||
"\n",
|
||
" # Let's pick maximum absolute displacement in y direction\n",
|
||
" maxdisp = 0.0\n",
|
||
" for el in elements\n",
|
||
" eldisp = el.attributes[\"displacement\"]\n",
|
||
" eldispy = eldisp[2,:]\n",
|
||
" maxeldisp = maximum(abs(eldispy))\n",
|
||
" if maxeldisp > maxdisp\n",
|
||
" maxdisp = maxeldisp\n",
|
||
" end\n",
|
||
" end\n",
|
||
" Logging.info(\"Maximum absolute displacement in y direction: $maxdisp\")\n",
|
||
" return model, elements, dofmap\n",
|
||
"end\n",
|
||
"\n",
|
||
"model, elements, dofmap = solve_3d_model();"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"## Saving results to file"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 15,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"<Grid Name=\"Grid\">\n",
|
||
" <Time Value=\"0\"/>\n",
|
||
"</Grid>\n"
|
||
]
|
||
},
|
||
"execution_count": 15,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"xdoc, xmodel = JuliaFEM.xdmf.xdmf_new_model()\n",
|
||
"temporal_collection = JuliaFEM.xdmf.xdmf_new_temporal_collection(xmodel)\n",
|
||
"grid = JuliaFEM.xdmf.xdmf_new_grid(temporal_collection; time=0)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Save geometry to xdmf file"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 16,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"name": "stderr",
|
||
"output_type": "stream",
|
||
"text": [
|
||
"13-Aug 00:37:51:INFO:root:Number of nodes in model: 298\n",
|
||
"13-Aug 00:37:51:INFO:root:Number of elements in model: 120\n"
|
||
]
|
||
}
|
||
],
|
||
"source": [
|
||
"nnodes = length(model[\"nodes\"])\n",
|
||
"Logging.info(\"Number of nodes in model: $nnodes\")\n",
|
||
"node_ids = Int64[]\n",
|
||
"for nid in keys(model[\"nodes\"])\n",
|
||
" push!(node_ids, nid)\n",
|
||
"end\n",
|
||
"sort!(node_ids)\n",
|
||
"X = zeros(3, nnodes)\n",
|
||
"for (i, nid) in enumerate(node_ids)\n",
|
||
" X[:,i] = model[\"nodes\"][nid]\n",
|
||
"end\n",
|
||
"\n",
|
||
"nelements = length(model[\"elements\"])\n",
|
||
"Logging.info(\"Number of elements in model: $nelements\")\n",
|
||
"elmap = zeros(Int64, 11, nelements)\n",
|
||
"elmap[1,:] = 0x0026\n",
|
||
"for elid in 1:nelements\n",
|
||
" elmap[2:end,elid] = model[\"elements\"][elid]\n",
|
||
"end"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 17,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"true"
|
||
]
|
||
},
|
||
"execution_count": 17,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"using LightXML\n",
|
||
"\n",
|
||
"function xdmf_new_mesh(grid, X, elmap)\n",
|
||
" dim, nnodes = size(X)\n",
|
||
" geometry = new_child(grid, \"Geometry\")\n",
|
||
" set_attribute(geometry, \"Type\", \"XYZ\")\n",
|
||
" dataitem = new_child(geometry, \"DataItem\")\n",
|
||
" set_attribute(dataitem, \"DataType\", \"Float\")\n",
|
||
" set_attribute(dataitem, \"Dimensions\", \"$nnodes $dim\")\n",
|
||
" set_attribute(dataitem, \"Format\", \"XML\")\n",
|
||
" set_attribute(dataitem, \"Precision\", 8)\n",
|
||
" #add_text(dataitem, join(X, \" \"))\n",
|
||
" s = \"\\n\"\n",
|
||
" \n",
|
||
" for i=1:nnodes\n",
|
||
" s *= \"\\t\\t\" * join(X[:,i], \" \") * \"\\n\"\n",
|
||
" end\n",
|
||
" s *= \" \"\n",
|
||
" add_text(dataitem, s)\n",
|
||
"\n",
|
||
" elmap2 = copy(elmap)\n",
|
||
" elmap2[2:end,:] -= 1\n",
|
||
" dim, nelements = size(elmap2)\n",
|
||
"\n",
|
||
" topology = new_child(grid, \"Topology\")\n",
|
||
" #set_attribute(topology, \"Dimensions\", \"1\")\n",
|
||
" set_attribute(topology, \"TopologyType\", \"Mixed\")\n",
|
||
" set_attribute(topology, \"NumberOfElements\", nelements)\n",
|
||
" dataitem = new_child(topology, \"DataItem\")\n",
|
||
" set_attribute(dataitem, \"DataType\", \"Int\")\n",
|
||
" set_attribute(dataitem, \"Dimensions\", \"$nelements $dim\")\n",
|
||
" set_attribute(dataitem, \"Format\", \"XML\")\n",
|
||
" set_attribute(dataitem, \"Precision\", 8)\n",
|
||
" s = \"\\n\"\n",
|
||
" for i=1:nelements\n",
|
||
" s *= \"\\t\\t\" * join(elmap2[:,i], \" \") * \"\\n\"\n",
|
||
" end\n",
|
||
" add_text(dataitem, s)\n",
|
||
" #add_text(dataitem, join(elmap2, \" \"))\n",
|
||
" \n",
|
||
"end\n",
|
||
"\n",
|
||
"xdmf_new_mesh(grid, X, elmap)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 18,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"10435"
|
||
]
|
||
},
|
||
"execution_count": 18,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"JuliaFEM.xdmf.xdmf_save_model(xdoc, \"/tmp/3d_solid_model.xmf\")"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "markdown",
|
||
"metadata": {},
|
||
"source": [
|
||
"Save nodal data to model"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 19,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"Dict{Any,Any} with 298 entries:\n",
|
||
" 288 => [-0.11672879201601734,-4.945796660856531,0.34641446977739976]\n",
|
||
" 11 => [-0.10351195894967631,-4.756564623966246,0.3419475317395125]\n",
|
||
" 158 => [-0.004781967238459207,-0.15800625216143013,-0.03641750463465448]\n",
|
||
" 215 => [-0.009679264390442352,-3.1283115653077243,0.1060066992444004]\n",
|
||
" 134 => [-0.012153531002752209,-2.568061461984137,0.1358780416251414]\n",
|
||
" 160 => [0.045841800099574996,-0.3952996715282528,-0.04337647899598485]\n",
|
||
" 29 => [-0.019610173840903154,-0.13106872577808154,-0.04187494739380516]\n",
|
||
" 131 => [0.04911611452200464,-2.5861023637514955,0.15290704689462717]\n",
|
||
" 249 => [-0.18589774363886435,-4.932931637644884,0.27789363131184]\n",
|
||
" 207 => [0.05588675730587902,-2.563471323701039,0.0669746466494972]\n",
|
||
" 173 => [-0.017826777186855432,-3.4323824666471117,0.08923268461041232]\n",
|
||
" 289 => [-0.08168560138254058,-4.782381469728047,0.4566974667444831]\n",
|
||
" 74 => [-0.013759506382279562,-2.9054222587693794,0.1389413161160111]\n",
|
||
" 201 => [0.006183182231407451,-2.8511450493652664,0.14659955727801022]\n",
|
||
" 176 => [0.04229935999959756,-2.349387641785469,0.09405299762367403]\n",
|
||
" 57 => [-0.0022029054234204426,-3.062804411311773,0.14846103116563106]\n",
|
||
" 31 => [0.005313697788644969,-0.25237461591643895,-0.030307540603490793]\n",
|
||
" 285 => [-0.10859185872935265,-4.832314881155342,0.35832378735042353]\n",
|
||
" 70 => [0.042811232565748765,-2.7950866978216338,0.1421452814276251]\n",
|
||
" 33 => [-0.016450709743293015,-0.379272230318473,0.027549328283933208]\n",
|
||
" 252 => [-0.14470395638682446,-0.9284887677850233,0.043594035381458736]\n",
|
||
" 114 => [-0.14842327883584733,-0.3532980011491899,-0.018641651848308006]\n",
|
||
" 165 => [-0.04346921888336824,-3.2745129973662452,0.07316785266950793]\n",
|
||
" 96 => [0.1501458469524447,-3.8471140528652787,0.1489845291162686]\n",
|
||
" 133 => [0.09045793517058681,-2.3950191788705957,0.12799724500423504]\n",
|
||
" ⋮ => ⋮"
|
||
]
|
||
},
|
||
"execution_count": 19,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"nodaldisp = Dict()\n",
|
||
"for el in elements\n",
|
||
" for (i, nid) in enumerate(el.node_ids)\n",
|
||
" nodaldisp[nid] = el.attributes[\"displacement\"][:, i]\n",
|
||
" end\n",
|
||
"end\n",
|
||
"nodaldisp"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 20,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"3x298 Array{Float64,2}:\n",
|
||
" -0.174259 -0.175081 0.0359566 0.0279021 0.00925611 0.0379208 … -0.0248959 -0.116472 0.191156 0.193018 -0.0153311 -0.123282 0.0\n",
|
||
" -0.978534 -1.32327 -2.42646 -2.50175 -2.72095 -2.59357 -3.20696 -0.707364 -1.76025 -1.33439 -2.22191 -1.93398 0.0\n",
|
||
" 0.0535874 0.111758 0.125464 0.115461 0.145374 0.1251 0.0902774 0.00418633 0.368958 0.387435 0.185441 0.230546 0.0"
|
||
]
|
||
},
|
||
"execution_count": 20,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"u = zeros(3, nnodes)\n",
|
||
"for (i, nid) in enumerate(node_ids)\n",
|
||
" u[:,i] = nodaldisp[nid]\n",
|
||
"end\n",
|
||
"u"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 21,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"true"
|
||
]
|
||
},
|
||
"execution_count": 21,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"JuliaFEM.xdmf.xdmf_new_field(grid, \"Displacement\", \"nodes\", u)"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": 22,
|
||
"metadata": {
|
||
"collapsed": false
|
||
},
|
||
"outputs": [
|
||
{
|
||
"data": {
|
||
"text/plain": [
|
||
"28215"
|
||
]
|
||
},
|
||
"execution_count": 22,
|
||
"metadata": {},
|
||
"output_type": "execute_result"
|
||
}
|
||
],
|
||
"source": [
|
||
"JuliaFEM.xdmf.xdmf_save_model(xdoc, \"/tmp/3d_solid_model.xmf\")"
|
||
]
|
||
},
|
||
{
|
||
"cell_type": "code",
|
||
"execution_count": null,
|
||
"metadata": {
|
||
"collapsed": true
|
||
},
|
||
"outputs": [],
|
||
"source": []
|
||
}
|
||
],
|
||
"metadata": {
|
||
"kernelspec": {
|
||
"display_name": "Julia 0.4.0-dev",
|
||
"language": "julia",
|
||
"name": "julia-0.4"
|
||
},
|
||
"language_info": {
|
||
"name": "julia",
|
||
"version": "0.4.0"
|
||
}
|
||
},
|
||
"nbformat": 4,
|
||
"nbformat_minor": 0
|
||
}
|