Skip to content

Commit

Permalink
add example computing stiffness
Browse files Browse the repository at this point in the history
  • Loading branch information
@KristofferC
KristofferC committed Mar 16, 2016
1 parent f3aadb6 commit 4b27076
Showing 1 changed file with 275 additions and 0 deletions.
275 changes: 275 additions & 0 deletions examples/stiffness_example.ipynb
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,275 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example creating the stiffness for a linear elasticity element using tensors or matrices"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuAFEM\n",
"using ForwardDiff"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"E = 200e9\n",
"ν = 0.3\n",
"λ = E*ν / ((1+ν) * (1 - 2ν))\n",
"μ = E / (2(1+ν))\n",
"δ(i,j) = i == j ? 1.0 : 0.0\n",
"g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))\n",
"\n",
"# Create a random symmetric material stiffness\n",
"C = rand(SymmetricTensor{4, 2})\n",
"\n",
"Ee = [C[1,1,1,1] C[2,2,1,1] C[1,2,1,1];\n",
" C[1,1,2,2] C[2,2,2,2] C[1,2,2,2];\n",
" C[1,1,1,2] C[2,2,1,2] C[1,2,1,2]];"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false,
"scrolled": true
},
"outputs": [],
"source": [
"function_space = Lagrange{2, JuAFEM.Square, 1}()\n",
"quad_rule = get_gaussrule(Dim{2}, JuAFEM.Square(), 2)\n",
"fe_values = FEValues(Float64, quad_rule, function_space);\n",
"\n",
"x = [0. 1 1 0;\n",
" 0 0 1 1]\n",
"x_vec = reinterpret(Vec{2, Float64}, x, (4,));"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Stiffness"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"ke_element_mat! (generic function with 1 method)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ke_element_mat!{T,dim}(Ke, X::Vector{Vec{2, T}}, fe_values::FEValues{dim}, Ee, B, DB, BDB)\n",
" n_basefuncs = n_basefunctions(get_functionspace(fe_values))\n",
" @assert length(X) == n_basefuncs\n",
" \n",
" reinit!(fe_values, X)\n",
"\n",
" for q_point in 1:length(JuAFEM.points(get_quadrule(fe_values)))\n",
" \n",
" for i in 1:n_basefuncs\n",
" dNdx = shape_gradient(fe_values, q_point, i)[1]\n",
" dNdy = shape_gradient(fe_values, q_point, i)[2]\n",
" B[1, 2*i - 1] = dNdx\n",
" B[2, 2*i - 0] = dNdy\n",
" B[3, 2*i - 0] = dNdx\n",
" B[3, 2*i - 1] = dNdy\n",
" end\n",
" \n",
" A_mul_B!(DB, Ee, B)\n",
" At_mul_B!(BDB, B, DB)\n",
" scale!(BDB, detJdV(fe_values, q_point))\n",
" for p in 1:size(Ke,1)\n",
" for q in 1:size(Ke,2)\n",
" Ke[p, q] += BDB[p, q]\n",
" end\n",
" end\n",
" end\n",
" \n",
" return Ke\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 0.404717 0.477524 -0.0143671 … -0.370061 -0.113279 -0.0871745\n",
" 0.469551 0.430669 -0.0626017 -0.300082 -0.0425479 -0.0237196\n",
" -0.244121 -0.293388 0.105868 0.185925 0.0217783 0.155135 \n",
" 0.203806 0.115946 -0.0489529 -0.246533 0.154102 0.0389071\n",
" -0.277071 -0.370061 -0.113279 0.477524 -0.0143671 -0.0202885\n",
" -0.364402 -0.300082 -0.0425479 … 0.430669 -0.0626017 -0.106868 \n",
" 0.116475 0.185925 0.0217783 -0.293388 0.105868 -0.0476721\n",
" -0.308955 -0.246533 0.154102 0.115946 -0.0489529 0.0916802"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"B = zeros(3, 8)\n",
"DB = zeros(3,8)\n",
"BDB = zeros(8,8)\n",
"Ke = zeros(8,8)\n",
"fill!(Ke, 0.0)\n",
"ke_element_mat!(Ke, x_vec, fe_values, Ee, B, DB, BDB)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"ke_elementt! (generic function with 1 method)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function ke_elementt!{T,dim}(Ke, X::Vector{Vec{2, T}}, fe_values::FEValues{dim}, D)\n",
" n_basefuncs = n_basefunctions(get_functionspace(fe_values))\n",
" @assert length(X) == n_basefuncs\n",
" reinit!(fe_values, X)\n",
" for q_point in 1:length(JuAFEM.points(get_quadrule(fe_values)))\n",
" for i in 1:n_basefuncs\n",
" for j in 1:n_basefuncs\n",
" ∇ϕi = shape_gradient(fe_values, q_point, i)\n",
" ∇ϕj = shape_gradient(fe_values, q_point, j)\n",
" Ke_e = (∇ϕj ⊗ ∇ϕi) ⊡ D * detJdV(fe_values, q_point)\n",
" Ke[dim*(i-1) + 1:dim*i, dim*(j-1) + 1:dim*j] += Ke_e\n",
" end\n",
" end\n",
" end\n",
" return Ke\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 0.519594 0.336348 -0.129244 … 0.00159777 -0.175752 \n",
" 0.336348 0.319262 0.070602 -0.175752 -0.135126 \n",
" -0.129244 0.070602 -0.00900904 0.136655 0.0208988 \n",
" 0.070602 0.00453899 0.0842508 0.0208988 -0.0724996 \n",
" -0.391947 -0.231198 0.00159777 -0.129244 0.070602 \n",
" -0.231198 -0.188675 -0.175752 … 0.070602 0.00453899\n",
" 0.00159777 -0.175752 0.136655 -0.00900904 0.0842508 \n",
" -0.175752 -0.135126 0.0208988 0.0842508 0.203087 "
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Note the minor transpose!\n",
"C2 = Tensor{4, 2}((i,l,k,j) -> C[i,j,k,l])\n",
"x = [0. 1 1 0;\n",
" 0 0 1 1]\n",
"x_vec = reinterpret(Vec{2, Float64}, x, (4,))\n",
"Ke2 = zeros(8,8)\n",
"ke_elementt!(Ke2, x_vec, fe_values, C2)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"0.5443978289896796"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(Ke - Ke2) / norm(Ke)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.3",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}

0 comments on commit 4b27076

Please sign in to comment.