mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
removed obsolete notebooks
This commit is contained in:
Jukka Aho
@@ -1,191 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Let's collect here examples how to parse mesh files (i.e. topologies; nodes and connectivity) from different kind of sources."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"f = open(\"../testdata/tripod.GiD.msh\");\n",
|
||||
"data = readlines(f)\n",
|
||||
"close(f)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Dict{Any,Any} with 0 entries"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"type Node\n",
|
||||
" node_id\n",
|
||||
" coords\n",
|
||||
"end\n",
|
||||
"type Element\n",
|
||||
" element_id\n",
|
||||
" node_ids\n",
|
||||
"end\n",
|
||||
"nodes = Dict()\n",
|
||||
"elements = Dict()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Entering section Coordinates\n",
|
||||
"Leaving section Coordinates\n",
|
||||
"Entering section Elements\n",
|
||||
"Leaving section Elements\n",
|
||||
"Mesh parsed\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"section = \"\"\n",
|
||||
"for line in data\n",
|
||||
" line = strip(line)\n",
|
||||
" if line == \"Coordinates\" || line == \"Elements\"\n",
|
||||
" section = line\n",
|
||||
" println(\"Entering section \",section)\n",
|
||||
" end\n",
|
||||
" if beginswith(line, \"end\")\n",
|
||||
" println(\"Leaving section \",section)\n",
|
||||
" section = \"\"\n",
|
||||
" end\n",
|
||||
" if section == \"Coordinates\"\n",
|
||||
" m = matchall(r\"[-0-9.]+\", line)\n",
|
||||
" if length(m) == 0\n",
|
||||
" continue\n",
|
||||
" end\n",
|
||||
" id = integer(m[1])\n",
|
||||
" coords = float(m[2:end])\n",
|
||||
" nodes[id] = Node(id, coords)\n",
|
||||
" end\n",
|
||||
" if section == \"Elements\"\n",
|
||||
" #println(line)\n",
|
||||
" m = matchall(r\"[0-9]+\", line)\n",
|
||||
" #println(m)\n",
|
||||
" if length(m) == 0\n",
|
||||
" continue\n",
|
||||
" end\n",
|
||||
" id = integer(m[1])\n",
|
||||
" #println(line)\n",
|
||||
" #println(m)\n",
|
||||
" #print(id)\n",
|
||||
" node_ids = m[2:end]\n",
|
||||
" #println(id,node_ids)\n",
|
||||
" elements[id] = Element(id, map(integer, node_ids))\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"println(\"Mesh parsed\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"10-element Array{Int64,1}:\n",
|
||||
" 1411\n",
|
||||
" 1527\n",
|
||||
" 1338\n",
|
||||
" 1581\n",
|
||||
" 1467\n",
|
||||
" 1431\n",
|
||||
" 1375\n",
|
||||
" 1494\n",
|
||||
" 1557\n",
|
||||
" 1455"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"elements[100].node_ids"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"10-element Array{Node,1}:\n",
|
||||
" Node(1411,[1.55541,4.93615,12.9066])\n",
|
||||
" Node(1527,[3.55077,4.78495,12.5057])\n",
|
||||
" Node(1338,[1.26995,2.70013,12.9378])\n",
|
||||
" Node(1581,[2.09046,5.60655,10.2598])\n",
|
||||
" Node(1467,[2.56484,4.86055,12.7445])\n",
|
||||
" Node(1431,[2.42483,3.74254,12.7719])\n",
|
||||
" Node(1375,[1.41281,3.81814,12.923]) \n",
|
||||
" Node(1494,[1.82293,5.27135,11.5832])\n",
|
||||
" Node(1557,[2.82062,5.19575,11.3827])\n",
|
||||
" Node(1455,[1.68021,4.15334,11.5988])"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"map(id -> nodes[id], elements[100].node_ids)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
File diff suppressed because one or more lines are too long
@@ -1,565 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#Assembly\n",
|
||||
"\n",
|
||||
"Let's make two element model and assemble it. We split the previous one element model to 2 quadrilaterals, make assembly and solve it. Small modifications to functions, I think it's better that they don't allocate memory but do in place operations. \n",
|
||||
"\n",
|
||||
"**TODO**\n",
|
||||
"- Tangent stiffness is calculated using forward finite difference. I think we should try ReverseDiffSparse for it's sparse matrix support, but I don't know how to use it. Or alternatively use FAD like before and assemble after linearization. It would be nice experiment to try linearization *after* assembly, would it work?\n",
|
||||
"- Verify calculations using some well known FEM software.\n",
|
||||
"\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": [
|
||||
"0x00000c95556709c2"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"tic()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Dict{Any,Any} with 0 entries"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"type Node\n",
|
||||
" coords\n",
|
||||
"end\n",
|
||||
"type Element\n",
|
||||
" node_ids\n",
|
||||
"end\n",
|
||||
"elements = Dict()\n",
|
||||
"nodes = Dict()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Element([5,6,3,2])"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"ndim = 2\n",
|
||||
"nnodes = 6\n",
|
||||
"\n",
|
||||
"nodes[1] = Node([0, 1, 0])\n",
|
||||
"nodes[2] = Node([5, 1, 0])\n",
|
||||
"nodes[3] = Node([10, 1, 0])\n",
|
||||
"nodes[4] = Node([0, 0, 0])\n",
|
||||
"nodes[5] = Node([5, 0, 0])\n",
|
||||
"nodes[6] = Node([10, 0, 0])\n",
|
||||
"elements[1] = Element([4, 5, 2, 1])\n",
|
||||
"elements[2] = Element([5, 6, 3, 2])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(36.0,24.0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# 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": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"calculate_internal_energy! (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function calculate_internal_energy!(X, u, Wint, dNdξ, λ, μ, dim=2)\n",
|
||||
" \"\"\"Calculate internal energy for a single element.\n",
|
||||
"\n",
|
||||
" Parameters\n",
|
||||
" ----------\n",
|
||||
" X : array [dim x nodes]\n",
|
||||
" u : array [dim x nodes]\n",
|
||||
" dNdξ : shape function derivatives\n",
|
||||
" λ : float\n",
|
||||
" μ : float\n",
|
||||
" dim : integer, optinal\n",
|
||||
"\n",
|
||||
" Returns\n",
|
||||
" -------\n",
|
||||
" Nothing, this is inplace function\n",
|
||||
" \n",
|
||||
" \"\"\"\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",
|
||||
" Wint[:,:] += w*J(ξ)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x4 Array{Float64,2}:\n",
|
||||
" -60.6667 -1.33333 133.333 -71.3333\n",
|
||||
" -24.0 -8.0 35.0 -3.0 "
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"Wint = zeros(2, 4)\n",
|
||||
"X = [0 0; 1 0; 1 1; 0 1]'\n",
|
||||
"u = [0 0; 0 0; 1 0; 0 0]'\n",
|
||||
"calculate_internal_energy!(X, u, Wint, dNdξ, λ, μ)\n",
|
||||
"Wint"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"assemble! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function assemble!(u, R)\n",
|
||||
" \"\"\" Assemble global residual vector R = T - F\n",
|
||||
" \"\"\"\n",
|
||||
" R[:] = 0.0\n",
|
||||
" u = reshape(u, ndim, nnodes)\n",
|
||||
" R = reshape(R, ndim, nnodes)\n",
|
||||
" \n",
|
||||
" Xe = zeros(ndim, 4)\n",
|
||||
" Winte = zeros(ndim, 4)\n",
|
||||
" \n",
|
||||
" # Internal forces, T\n",
|
||||
" for i=1:length(elements) # loop through elements\n",
|
||||
" Xe[:,:] = 0.0 # FIXME: how to efficiently empty array?\n",
|
||||
" el = elements[i]\n",
|
||||
" nids = el.node_ids\n",
|
||||
" for i=1:length(nids) # loop through nodes\n",
|
||||
" Xe[:,i] = nodes[nids[i]].coords[1:2]\n",
|
||||
" end\n",
|
||||
" Winte[:,:] = 0.0\n",
|
||||
" calculate_internal_energy!(Xe, u[:,nids], Winte, dNdξ, λ, μ)\n",
|
||||
" R[:,nids] += Winte\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # External forces, F\n",
|
||||
" # T - F = T + (-F)\n",
|
||||
" R[2, 3] += 2 # Force to the tip of härveli\n",
|
||||
"\n",
|
||||
" u = reshape(u, ndim*nnodes)\n",
|
||||
" R = reshape(R, ndim*nnodes)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1x12 Array{Float64,2}:\n",
|
||||
" 0.0 0.0 34.8853 16.2 140.715 78.6 … -64.1013 -36.0 -111.499 -56.8"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"R = zeros(ndim*nnodes)\n",
|
||||
"u = [0 0; 0 0; 1 0; 0 0; 0 0; 0 0]'\n",
|
||||
"#u = [0 0; 0 0; 1 0; 0 0]'\n",
|
||||
"assemble!(reshape(u, ndim*nnodes), R)\n",
|
||||
"R'"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We should use ReverseDiffSparse because of it's sparse matrix support"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# FIXME: I don't know how to get these working!\n",
|
||||
"#using ReverseDiffSparse\n",
|
||||
"#using ForwardDiff\n",
|
||||
"#Kt! = forwarddiff_jacobian!(assemble!, Float64, fadtype=:dual, n=12, m=12)\n",
|
||||
"#Kt!(reshape(u, 12), Kt)\n",
|
||||
"#Kt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Lin (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# So we go to plan B. FIXME: change this to analytical version\n",
|
||||
"\n",
|
||||
"function Lin(f!, h=1.0e-6)\n",
|
||||
"\n",
|
||||
" function D!(x, J)\n",
|
||||
" J[:,:] = 0\n",
|
||||
" N = length(x)\n",
|
||||
" Δx = zeros(N)\n",
|
||||
" y = zeros(N)\n",
|
||||
" Δy = zeros(N)\n",
|
||||
" f!(x, y) # Evaluate function f in x and store results to y\n",
|
||||
" for i=1:N\n",
|
||||
" Δx[:] = 0.0\n",
|
||||
" Δx[i] += h\n",
|
||||
" f!(x+Δx, Δy) # Evaluate function f in x+Δx and store results to Δy\n",
|
||||
" J[i, :] = (Δy-y) / h\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" return D!\n",
|
||||
"\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Handling homogeneous Dirichlet conditions, using elimination.\n",
|
||||
"\n",
|
||||
"**INFO**: We could try something like this: http://www.code-aster.org/V2/doc/default/en/man_r/r3/r3.03.01.pdf\n",
|
||||
"\n",
|
||||
"Here's an idea how to make a very simply elimination"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x6 Array{Int64,2}:\n",
|
||||
" 1 0 0 1 0 0\n",
|
||||
" 1 0 0 1 0 0"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fixed_dofs = integer(zeros(ndim, nnodes))\n",
|
||||
"fixed_dofs[:,1] = fixed_dofs[:,4] = 1\n",
|
||||
"fixed_dofs"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x6 Array{Int64,2}:\n",
|
||||
" 0 1 1 0 1 1\n",
|
||||
" 0 1 1 0 1 1"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"free_dofs = integer(ones(ndim, nnodes)) - fixed_dofs\n",
|
||||
"free_dofs"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1x8 Array{Int64,2}:\n",
|
||||
" 3 4 5 6 9 10 11 12"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"free_dofs = find(free_dofs)\n",
|
||||
"free_dofs'"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Starting Newton iterations\n",
|
||||
"Iteration 1, norm = 12.031651257897487\n",
|
||||
"Iteration 2, norm = 3.141275068747564\n",
|
||||
"Iteration 3, norm = 1.1530769747882408\n",
|
||||
"Iteration 4, norm = 0.2665784114781298\n",
|
||||
"Iteration 5, norm = 0.035146969229970175\n",
|
||||
"Iteration 6, norm = 0.003125584242156337\n",
|
||||
"Iteration 7, norm = 2.227371173806294e-6\n",
|
||||
"Iteration 8, norm = 2.1397676611257455e-10\n",
|
||||
"Converged.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x6 Array{Float64,2}:\n",
|
||||
" 0.0 -0.106192 -1.69115 0.0 -0.761252 -2.48604\n",
|
||||
" 0.0 -2.10509 -6.00728 0.0 -1.89221 -5.5914 "
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function solve!(u, free_dofs; max_iterations=10, eps=1.0e-7)\n",
|
||||
" ∇R = Lin(assemble!)\n",
|
||||
" Kt = zeros(ndim*nnodes, ndim*nnodes)\n",
|
||||
" R = zeros(size(u))\n",
|
||||
" println(\"Starting Newton iterations\")\n",
|
||||
" for i=1:max_iterations\n",
|
||||
" print(\"Iteration \",i, \", \")\n",
|
||||
" R[:] = 0.0\n",
|
||||
" assemble!(u, R) # Calculate internal energy in nodes and store results to R\n",
|
||||
" ∇R(u, Kt) # Linearize residual in u and save result to Kt\n",
|
||||
" # Solve !\n",
|
||||
" du = Kt[free_dofs, free_dofs] \\ -R[free_dofs]\n",
|
||||
" u[free_dofs] += du\n",
|
||||
" println(\"norm = \",norm(du))\n",
|
||||
" if norm(du) < eps\n",
|
||||
" println(\"Converged.\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" return u\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"u = zeros(ndim*nnodes)\n",
|
||||
"solve!(u, free_dofs)\n",
|
||||
"u = reshape(u, ndim, nnodes)\n",
|
||||
"u"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"elapsed time: 4.541729719 seconds\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4.541729719"
|
||||
]
|
||||
},
|
||||
"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
|
||||
}
|
||||
@@ -1,382 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Compare analytical and autodiffed stiffness matrix\n",
|
||||
"\n",
|
||||
"Here we compare how much autodiffed solution is slower than analytical.\n",
|
||||
"\n",
|
||||
"Author(s): Jukka Aho <jukka.aho@kapsi.fi>\n",
|
||||
"\n",
|
||||
"Last updated:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2015-06-15"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using Dates\n",
|
||||
"today()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"160"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using ForwardDiff\n",
|
||||
"ENV[\"COLUMNS\"] = 160"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(36.0,24.0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# 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",
|
||||
"a = 1/sqrt(3)\n",
|
||||
"ipoints = [[-a -a], [a -a], [a a], [-a a]]\n",
|
||||
"iweights = [1 1 1 1]\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": [
|
||||
"## Version using automatic differentiation"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"calc_local_matrices! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function calc_local_matrices!(X, u, R, Kt; dim=2)\n",
|
||||
" I = eye(dim)\n",
|
||||
" \n",
|
||||
" function calc_Wint!(u, Wint)\n",
|
||||
" for m = 1:length(iweights)\n",
|
||||
" w = iweights[m]\n",
|
||||
" ξ = ipoints[m, :]\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",
|
||||
" S = λ*trace(E)*I + 2*μ*E # PK2 stress tensor\n",
|
||||
" P = F*S # PK1 stress tensor\n",
|
||||
" Wint[:,:] += w*P*∇N*det(Jᵀ)\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # herlper for tangent stiffness matrix\n",
|
||||
" function R!(u, R)\n",
|
||||
" R[:] = 0\n",
|
||||
" calc_Wint!(reshape(u, 2, 4), reshape(R, 2, 4))\n",
|
||||
" #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))\n",
|
||||
" end\n",
|
||||
" Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=8, m=8)\n",
|
||||
"\n",
|
||||
" Kt[:,:] = Jacobian(reshape(u, 8))\n",
|
||||
" R!(reshape(u, 8), reshape(R, 8))\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Converged in 6 iterations\n",
|
||||
"[0.0 -0.3991450609547433 -0.07228582695592461 0.0\n",
|
||||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# validation\n",
|
||||
"X = [0 0; 10 0; 10 1; 0 1]'\n",
|
||||
"u = zeros(2,4)\n",
|
||||
"R = zeros(2,4)\n",
|
||||
"Kt = zeros(8,8)\n",
|
||||
"\n",
|
||||
"free_dofs = [3, 4, 5, 6]\n",
|
||||
"for i in 1:10\n",
|
||||
" calc_local_matrices!(X, u, R, Kt)\n",
|
||||
" R[2,3] += 2\n",
|
||||
" du = Kt[free_dofs, free_dofs] \\ -reshape(R, 8)[free_dofs]\n",
|
||||
" u[free_dofs] += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged in \", i, \" iterations\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"println(u)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function test_algo1(N=10000)\n",
|
||||
" for i=1:N\n",
|
||||
" calc_local_matrices!(X, u, R, Kt)\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"test_algo1()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"elapsed time: 11.336512664 seconds (2020893880 bytes allocated, 26.15% gc time)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"@time test_algo1()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Analytical tangent stiffness"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"calc_local_matrices2! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function calc_local_matrices2!(X, u, R, Kt; dim=2)\n",
|
||||
" I = eye(dim)\n",
|
||||
" R[:,:] = 0.0\n",
|
||||
" Kt[:,:] = 0.0\n",
|
||||
" N = 4 # number of shape functions\n",
|
||||
"\n",
|
||||
" dF = zeros(2, 2)\n",
|
||||
"\n",
|
||||
" for m = 1:length(iweights)\n",
|
||||
" w = iweights[m]\n",
|
||||
" ξ = ipoints[m, :]\n",
|
||||
" Jᵀ = X*dNdξ(ξ)\n",
|
||||
" detJ = det(Jᵀ)\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",
|
||||
" S = λ*trace(E)*I + 2*μ*E # PK2 stress tensor\n",
|
||||
" P = F*S # PK1 stress tensor\n",
|
||||
" R[:,:] += w*P*∇N*detJ\n",
|
||||
"\n",
|
||||
" for p = 1:N\n",
|
||||
" for i = 1:dim\n",
|
||||
" dF[:,:] = 0.0\n",
|
||||
" dF[i,:] = ∇N[:,p]\n",
|
||||
" dE = 1/2*(F'*dF + dF'*F)\n",
|
||||
" dS = λ*trace(dE)*I + 2*μ*dE\n",
|
||||
" dP = dF*S + F*dS\n",
|
||||
" for q = 1:N\n",
|
||||
" for j = 1:dim\n",
|
||||
" Kt[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*∇N[:,q])[1]*detJ\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Converged in 6 iterations\n",
|
||||
"[0.0 -0.39914506095474317 -0.07228582695592449 0.0\n",
|
||||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# validation\n",
|
||||
"X = [0 0; 10 0; 10 1; 0 1]'\n",
|
||||
"u = zeros(2,4)\n",
|
||||
"R = zeros(2,4)\n",
|
||||
"Kt = zeros(8,8)\n",
|
||||
"\n",
|
||||
"free_dofs = [3, 4, 5, 6]\n",
|
||||
"for i in 1:10\n",
|
||||
" calc_local_matrices2!(X, u, R, Kt)\n",
|
||||
" R[2,3] += 2\n",
|
||||
" du = Kt[free_dofs, free_dofs] \\ -reshape(R, 8)[free_dofs]\n",
|
||||
" u[free_dofs] += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged in \", i, \" iterations\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"println(u)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"function test_algo2(N=10000)\n",
|
||||
" for i=1:N\n",
|
||||
" calc_local_matrices2!(X, u, R, Kt)\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"test_algo2()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"elapsed time: 8.447644445 seconds (1554800080 bytes allocated, 29.45% gc time)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"@time test_algo2()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
File diff suppressed because one or more lines are too long
@@ -1,426 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Exporting data to Xdmf format\n",
|
||||
"\n",
|
||||
"Author(s): Jukka Aho <jukka.aho@kapsi.fi>\n",
|
||||
"\n",
|
||||
"Python example how to export data to Xdmf format using Python"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import sys\n",
|
||||
"import os\n",
|
||||
"sys.path.append(os.path.expanduser('~/opt/xdmf/lib/python'))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from Xdmf import *"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Export mesh"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"grid = XdmfUnstructuredGrid.New()\n",
|
||||
"grid.setTime(XdmfTime.New(123))\n",
|
||||
"geometry = XdmfGeometry.New()\n",
|
||||
"topology = XdmfTopology.New()\n",
|
||||
"geometry.setType(XdmfGeometryType.XYZ())\n",
|
||||
"topology.setType(XdmfTopologyType.Mixed())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"array([[ 0., 0., 0.],\n",
|
||||
" [ 1., 0., 0.],\n",
|
||||
" [ 0., 1., 0.],\n",
|
||||
" [ 0., 0., 1.]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"# one linear tetra\n",
|
||||
"nodes = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])\n",
|
||||
"nodes"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"It needs to be in 1d list:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"[0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"list(nodes.flatten())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"geometry.insertAsFloat32(0, list(nodes.flatten()))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Topology; connect nodes 1 - 4 to form tetra\n",
|
||||
"\n",
|
||||
"List is here:\n",
|
||||
"\n",
|
||||
"http://public.kiteware.com/pipermail/2013-July/028859.html"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"topology.insertAsInt32(0, [0x6, 1, 2, 3, 4]) # Not sure starts this from 0 or 1"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**update**, node numbering starts from 0."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"grid.setGeometry(geometry)\n",
|
||||
"grid.setTopology(topology)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Export fields"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Give scalar value to cells (i.e. elements)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"temperature = XdmfAttribute.New()\n",
|
||||
"temperature.setType(XdmfAttributeType.Scalar())\n",
|
||||
"temperature.setCenter(XdmfAttributeCenter.Cell())\n",
|
||||
"temperature.setName(\"Temperature field\")\n",
|
||||
"temperature.insertAsInt32(0, [56])\n",
|
||||
"grid.insert(temperature)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Give vector values to nodes"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"density = XdmfAttribute.New()\n",
|
||||
"density.setType(XdmfAttributeType.Vector())\n",
|
||||
"density.setCenter(XdmfAttributeCenter.Node())\n",
|
||||
"density.setName(\"Density\")\n",
|
||||
"density.insertAsFloat32(0, [1.3, 2.4, 4.5, 6.7])\n",
|
||||
"grid.insert(density)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Write to file"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"I think this is not mandatory but create temporal collection anyway:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"coll = XdmfGridCollection.New()\n",
|
||||
"coll.setType(XdmfGridCollectionType.Temporal())\n",
|
||||
"coll.insert(grid)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Create domain and write to disk"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"domain = XdmfDomain.New()\n",
|
||||
"domain.insert(coll)\n",
|
||||
"writer = XdmfWriter.New(\"testdata.xmf\")\n",
|
||||
"domain.accept(writer)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"testdata.h5 testdata.xmf\r\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"!ls testdata*"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"<?xml version=\"1.0\" encoding=\"utf-8\"?>\r\n",
|
||||
"<Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"2.1\">\r\n",
|
||||
" <Domain>\r\n",
|
||||
" <Grid CollectionType=\"Temporal\" GridType=\"Collection\" Name=\"Collection\">\r\n",
|
||||
" <Geometry Type=\"None\"/>\r\n",
|
||||
" <Topology Dimensions=\"0\" Type=\"NoTopology\"/>\r\n",
|
||||
" <Grid Name=\"Grid\">\r\n",
|
||||
" <Time Value=\"123\"/>\r\n",
|
||||
" <Geometry Type=\"XYZ\">\r\n",
|
||||
" <DataItem DataType=\"Float\" Dimensions=\"12\" Format=\"XML\" Precision=\"4\">0 0 0 1 0 0 0 1 0 0 0 1</DataItem>\r\n",
|
||||
" </Geometry>\r\n",
|
||||
" <Topology Dimensions=\"1\" Type=\"Mixed\">\r\n",
|
||||
" <DataItem DataType=\"Int\" Dimensions=\"5\" Format=\"XML\" Precision=\"4\">6 1 2 3 4</DataItem>\r\n",
|
||||
" </Topology>\r\n",
|
||||
" <Attribute Center=\"Cell\" Name=\"Temperature field\" Type=\"Scalar\">\r\n",
|
||||
" <DataItem DataType=\"Int\" Dimensions=\"1\" Format=\"XML\" Precision=\"4\">56</DataItem>\r\n",
|
||||
" </Attribute>\r\n",
|
||||
" <Attribute Center=\"Node\" Name=\"Density\" Type=\"Vector\">\r\n",
|
||||
" <DataItem DataType=\"Float\" Dimensions=\"4\" Format=\"XML\" Precision=\"4\">1.3 2.4000001 4.5 6.6999998</DataItem>\r\n",
|
||||
" </Attribute>\r\n",
|
||||
" </Grid>\r\n",
|
||||
" </Grid>\r\n",
|
||||
" </Domain>\r\n",
|
||||
"</Xdmf>\r\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"!cat testdata.xmf"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"This time Xdmf decided to write results directly to the xml file. If the model is bigger the correct format is"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
" <?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
|
||||
" <Xdmf xmlns:xi=\"http://www.w3.org/2001/XInclude\" Version=\"2.1\">\n",
|
||||
" <Domain>\n",
|
||||
" <Grid CollectionType=\"Temporal\" GridType=\"Collection\" Name=\"Collection\">\n",
|
||||
" <Geometry Type=\"None\"/>\n",
|
||||
" <Topology Dimensions=\"0\" Type=\"NoTopology\"/>\n",
|
||||
" <Grid Name=\"Grid\">\n",
|
||||
" <Time Value=\"0\"/>\n",
|
||||
" <Geometry Type=\"XYZ\">\n",
|
||||
" <DataItem DataType=\"Float\" Dimensions=\"1260999\" Format=\"HDF\" Precision=\"4\">model.h5:Data16</DataItem>\n",
|
||||
" </Geometry>\n",
|
||||
" <Topology Dimensions=\"292001\" Type=\"Mixed\">\n",
|
||||
" <DataItem DataType=\"Int\" Dimensions=\"3046831\" Format=\"HDF\" Precision=\"4\">model.h5:Data17</DataItem>\n",
|
||||
" </Topology>\n",
|
||||
" <Attribute Center=\"Cell\" Name=\"Temperature field\" Type=\"Scalar\">\n",
|
||||
" <DataItem DataType=\"Int\" Dimensions=\"292001\" Format=\"HDF\" Precision=\"4\">model.h5:Data18</DataItem>\n",
|
||||
" </Attribute>\n",
|
||||
" <Attribute Center=\"Node\" Name=\"Density\" Type=\"Vector\">\n",
|
||||
" <DataItem DataType=\"Float\" Dimensions=\"292001\" Format=\"HDF\" Precision=\"4\">model.h5:Data19</DataItem>\n",
|
||||
" </Attribute>\n",
|
||||
" </Grid>\n",
|
||||
" </Grid>\n",
|
||||
" </Domain>\n",
|
||||
" </Xdmf>\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Structure of h5 file\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" File(filename=model.h5, title='', mode='r+', root_uep='/', filters=Filters(complevel=0, shuffle=False, fletcher32=False, least_significant_digit=None)) \n",
|
||||
" / (RootGroup) '' \n",
|
||||
" /Data0 (EArray(1125135,)) ''\n",
|
||||
" atom := Float32Atom(shape=(), dflt=0.0)\n",
|
||||
" maindim := 0 \n",
|
||||
" flavor := 'numpy' \n",
|
||||
" byteorder := 'little' \n",
|
||||
" chunkshape := (1000,) \n",
|
||||
" /Data1 (EArray(2680009,)) ''\n",
|
||||
" atom := Int32Atom(shape=(), dflt=0)\n",
|
||||
" maindim := 0 \n",
|
||||
" flavor := 'numpy' \n",
|
||||
" byteorder := 'little' \n",
|
||||
" chunkshape := (1000,) \n",
|
||||
" ...\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" In [7]: d.root.Data0[:]\n",
|
||||
" Out[7]:\n",
|
||||
" array([ 302.10440063, 222.19999695, 20.38028908, ..., -230. ,\n",
|
||||
" 0. , 7.33333015], dtype=float32)\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.9"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
@@ -1,530 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Shape function and integration points\n",
|
||||
"\n",
|
||||
"Author(s): Jukka Aho\n",
|
||||
"\n",
|
||||
"**Abstract**: Shape functions and element descriptions used in JuliaFEM."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from sympy import *\n",
|
||||
"#init_printing()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"xi = DeferredVector(r\"xi\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## 1D shape function"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Linear 2-node segment (Lagrange family)\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ |\n",
|
||||
"| ----- | -------:|\n",
|
||||
"| $N_1$ | -1 |\n",
|
||||
"| $N_2$ | 1 |\n",
|
||||
"\n",
|
||||
"\\begin{equation}\n",
|
||||
" \\left(\\mathbf{P}\\boldsymbol{\\alpha}\\right)\\left(\\xi_1\\right) = \\alpha_{1}+\\alpha_{2}\\xi_{1}\n",
|
||||
"\\end{equation}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Matrix([\n",
|
||||
" [-xi[1]/2 + 1/2],\n",
|
||||
" [ xi[1]/2 + 1/2]]), Matrix([\n",
|
||||
" [-1/2],\n",
|
||||
" [ 1/2]]))"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"A = Matrix([[1, -1], [1, 1]])\n",
|
||||
"P = Matrix([1, xi[1]]).T\n",
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T]).T\n",
|
||||
"N, dN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Quadratic 3-node segment (Lagrange family)\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ |\n",
|
||||
"| ----- | -------:|\n",
|
||||
"| $N_1$ | -1 |\n",
|
||||
"| $N_2$ | 1 |\n",
|
||||
"| $N_3$ | 0 |\n",
|
||||
"\n",
|
||||
"\\begin{equation}\n",
|
||||
" \\left(\\mathbf{P}\\boldsymbol{\\alpha}\\right)\\left(\\xi_1\\right) = \\alpha_1 + \\alpha_2\\xi_1 + \\alpha_3\\xi_1^2\n",
|
||||
"\\end{equation}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Matrix([\n",
|
||||
" [xi[1]**2/2 - xi[1]/2],\n",
|
||||
" [xi[1]**2/2 + xi[1]/2],\n",
|
||||
" [ -xi[1]**2 + 1]]), Matrix([\n",
|
||||
" [xi[1] - 1/2],\n",
|
||||
" [xi[1] + 1/2],\n",
|
||||
" [ -2*xi[1]]]))"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"A = Matrix([[1, -1, (-1)**2],\n",
|
||||
" [1, 1, 1**2],\n",
|
||||
" [1, 0, 0**2]])\n",
|
||||
"P = Matrix([1, xi[1], xi[1]**2]).T\n",
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T]).T\n",
|
||||
"N, dN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### P-elements"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## 2D shape functions"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Linear triangle\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ | $\\xi_2$ |\n",
|
||||
"| ----- | -------:| -------:|\n",
|
||||
"| $N_1$ | 0 | 0 |\n",
|
||||
"| $N_2$ | 1 | 0 |\n",
|
||||
"| $N_3$ | 0 | 1 |"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Matrix([\n",
|
||||
" [-xi[1] - xi[2] + 1],\n",
|
||||
" [ xi[1]],\n",
|
||||
" [ xi[2]]]), Matrix([\n",
|
||||
" [-1, -1],\n",
|
||||
" [ 1, 0],\n",
|
||||
" [ 0, 1]]))"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"A = Matrix([[1, 0, 0], [1, 1, 0], [1, 0, 1]])\n",
|
||||
"P = Matrix([1, xi[1], xi[2]]).T\n",
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T]).T\n",
|
||||
"N, dN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Quadratic triangle\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ | $\\xi_2$ |\n",
|
||||
"| ----- | -------:| -------:|\n",
|
||||
"| $N_1$ | 0 | 0 |\n",
|
||||
"| $N_2$ | 1 | 0 |\n",
|
||||
"| $N_3$ | 0 | 1 |\n",
|
||||
"| $N_4$ | 1/2 | 0 |\n",
|
||||
"| $N_5$ | 1/2 | 1/2 |\n",
|
||||
"| $N_6$ | 0 | 1/2 |"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Matrix([\n",
|
||||
"[1, 0, 0, 0, 0, 0],\n",
|
||||
"[1, 1, 0, 1, 0, 0],\n",
|
||||
"[1, 0, 1, 0, 1, 0],\n",
|
||||
"[1, 1/2, 0, 1/4, 0, 0],\n",
|
||||
"[1, 1/2, 1/2, 1/4, 1/4, 1/4],\n",
|
||||
"[1, 0, 1/2, 0, 1/4, 0]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P = Matrix([1, xi[1], xi[2], xi[1]**2, xi[2]**2, xi[1]*xi[2]]).T\n",
|
||||
"A = Matrix([\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0}),\n",
|
||||
" P.subs({xi[1]: 1, xi[2]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 1}),\n",
|
||||
" P.subs({xi[1]: Rational(1,2), xi[2]: 0}),\n",
|
||||
" P.subs({xi[1]: Rational(1,2), xi[2]: Rational(1,2)}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: Rational(1,2)}),\n",
|
||||
" ])\n",
|
||||
"A"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Matrix([\n",
|
||||
" [2*xi[1]**2 + 4*xi[1]*xi[2] - 3*xi[1] + 2*xi[2]**2 - 3*xi[2] + 1],\n",
|
||||
" [ 2*xi[1]**2 - xi[1]],\n",
|
||||
" [ 2*xi[2]**2 - xi[2]],\n",
|
||||
" [ -4*xi[1]**2 - 4*xi[1]*xi[2] + 4*xi[1]],\n",
|
||||
" [ 4*xi[1]*xi[2]],\n",
|
||||
" [ -4*xi[1]*xi[2] - 4*xi[2]**2 + 4*xi[2]]]), Matrix([\n",
|
||||
" [ 4*xi[1] + 4*xi[2] - 3, 4*xi[1] + 4*xi[2] - 3],\n",
|
||||
" [ 4*xi[1] - 1, 0],\n",
|
||||
" [ 0, 4*xi[2] - 1],\n",
|
||||
" [-8*xi[1] - 4*xi[2] + 4, -4*xi[1]],\n",
|
||||
" [ 4*xi[2], 4*xi[1]],\n",
|
||||
" [ -4*xi[2], -4*xi[1] - 8*xi[2] + 4]]))"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T]).T\n",
|
||||
"N, dN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## 3D shape functions"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Linear tetrahedra, **tet4**\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ | $\\xi_2$ | $\\xi_2$ |\n",
|
||||
"| ----- | -------:| -------:| -------:|\n",
|
||||
"| $N_1$ | 0 | 0 | 0 |\n",
|
||||
"| $N_2$ | 1 | 0 | 0 |\n",
|
||||
"| $N_3$ | 0 | 1 | 0 |\n",
|
||||
"| $N_4$ | 0 | 0 | 1 |"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Matrix([\n",
|
||||
"[1, 0, 0, 0],\n",
|
||||
"[1, 1, 0, 0],\n",
|
||||
"[1, 0, 1, 0],\n",
|
||||
"[1, 0, 0, 1]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P = Matrix([1, xi[1], xi[2], xi[3]]).T\n",
|
||||
"A = Matrix([\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 1, xi[2]: 0, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 1, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 1}),\n",
|
||||
" ])\n",
|
||||
"A"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(Matrix([\n",
|
||||
" [-xi[1] - xi[2] - xi[3] + 1],\n",
|
||||
" [ xi[1]],\n",
|
||||
" [ xi[2]],\n",
|
||||
" [ xi[3]]]), Matrix([\n",
|
||||
" [-1, -1, -1],\n",
|
||||
" [ 1, 0, 0],\n",
|
||||
" [ 0, 1, 0],\n",
|
||||
" [ 0, 0, 1]]))"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T, N.diff(xi[3]).T]).T\n",
|
||||
"N, dN"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Quadratic Lagrange tetrahedral element, 10 nodes, **tet10**\n",
|
||||
"\n",
|
||||
"| | $\\xi_1$ | $\\xi_2$ | $\\xi_2$ |\n",
|
||||
"| ----- | -------:| -------:| -------:|\n",
|
||||
"| $N_1$ | 0 | 0 | 0 |\n",
|
||||
"| $N_2$ | 1 | 0 | 0 |\n",
|
||||
"| $N_3$ | 0 | 1 | 0 |\n",
|
||||
"| $N_4$ | 0 | 0 | 1 |\n",
|
||||
"| $N_5$ | 1/2 | 0 | 0 |\n",
|
||||
"| $N_6$ | 1/2 | 1/2 | 0 |\n",
|
||||
"| $N_7$ | 0 | 1/2 | 0 |\n",
|
||||
"| $N_8$ | 0 | 0 | 1/2 |\n",
|
||||
"| $N_9$ | 1/2 | 0 | 1/2 |\n",
|
||||
"| $N_{10}$ | 0 | 1/2 | 1/2 |"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Matrix([\n",
|
||||
"[1, 0, 0, 0, 0, 0, 0, 0, 0, 0],\n",
|
||||
"[1, 1, 0, 0, 1, 0, 0, 0, 0, 0],\n",
|
||||
"[1, 0, 1, 0, 0, 0, 1, 0, 0, 0],\n",
|
||||
"[1, 0, 0, 1, 0, 0, 0, 0, 1, 0],\n",
|
||||
"[1, 1/2, 0, 0, 1/4, 0, 0, 0, 0, 0],\n",
|
||||
"[1, 1/2, 1/2, 0, 1/4, 1/4, 1/4, 0, 0, 0],\n",
|
||||
"[1, 0, 1/2, 0, 0, 0, 1/4, 0, 0, 0],\n",
|
||||
"[1, 0, 0, 1/2, 0, 0, 0, 0, 1/4, 0],\n",
|
||||
"[1, 1/2, 0, 1/2, 1/4, 0, 0, 0, 1/4, 1/4],\n",
|
||||
"[1, 0, 1/2, 1/2, 0, 0, 1/4, 1/4, 1/4, 0]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"P = Matrix([1, xi[1], xi[2], xi[3], xi[1]**2, xi[1]*xi[2], xi[2]**2, xi[2]*xi[3], xi[3]**2, xi[1]*xi[3]]).T\n",
|
||||
"A = Matrix([\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 1, xi[2]: 0, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 1, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: 1}),\n",
|
||||
"\n",
|
||||
" P.subs({xi[1]: Rational(1,2), xi[2]: 0, xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: Rational(1,2), xi[2]: Rational(1,2), xi[3]: 0}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: Rational(1,2), xi[3]: 0}),\n",
|
||||
"\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: 0, xi[3]: Rational(1,2)}),\n",
|
||||
" P.subs({xi[1]: Rational(1,2), xi[2]: 0, xi[3]: Rational(1,2)}),\n",
|
||||
" P.subs({xi[1]: 0, xi[2]: Rational(1,2), xi[3]: Rational(1,2)}),\n",
|
||||
" ])\n",
|
||||
"A"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Matrix([\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]]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"N = (P*A.inv()).T\n",
|
||||
"dN = Matrix([N.diff(xi[1]).T, N.diff(xi[2]).T, N.diff(xi[3]).T]).T\n",
|
||||
"factor(N)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Matrix([\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]]])"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"factor(dN)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 2",
|
||||
"language": "python",
|
||||
"name": "python2"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 2
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython2",
|
||||
"version": "2.7.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
File diff suppressed because one or more lines are too long
@@ -1,548 +0,0 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Interpolation and integration algorithms\n",
|
||||
"\n",
|
||||
"Author(s): Jukka Aho\n",
|
||||
"\n",
|
||||
"**Abstract**: Some strategies to implement automatic differentiation. The number of different choises is caused by a fact that the linearization of function can be done before integration or vice versa, and functions can return values or do in-place modifications. There is probably performance differences between different strategies, but all of them should work."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM\n",
|
||||
"using ForwardDiff"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"\"Old good\" elasticity force equilibrium equation $R = T - F$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"calc_residual_vector_integrand (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n",
|
||||
" # Calculate dN/dX and interpolate material parameters\n",
|
||||
" dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n",
|
||||
" u = el.attributes[\"displacement\"]\n",
|
||||
" lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n",
|
||||
" mu = JuliaFEM.interpolate(el, \"mu\", xi)\n",
|
||||
"\n",
|
||||
" # Calculate residual force vector R = T - F\n",
|
||||
" gradu = u*dbasisdX\n",
|
||||
" F = I + gradu\n",
|
||||
" E = 1/2*(gradu' + gradu + gradu'*gradu)\n",
|
||||
" S = lambda*trace(E)*I + 2*mu*E\n",
|
||||
" P = F*S\n",
|
||||
" T = P*dbasisdX'\n",
|
||||
" return T\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Test case, already well known 2d elasticity in [0,10] x [0,1] grid."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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",
|
||||
"ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]'\n",
|
||||
"iweights = [1, 1, 1, 1]\n",
|
||||
"attributes = Dict()\n",
|
||||
"e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\n",
|
||||
"\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",
|
||||
"e.attributes[\"coordinates\"] = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
|
||||
"e.attributes[\"lambda\"] = la\n",
|
||||
"e.attributes[\"mu\"] = mu\n",
|
||||
"e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.5 0.0; 0.0 0.0]'\n",
|
||||
"e.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
|
||||
"e.attributes[\"displacement tangent stiffness\"] = zeros(8, 8);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Integration\n",
|
||||
"\n",
|
||||
"1. take element and function and return value\n",
|
||||
"2. take function and return function which can be integrated by passing element as a function\n",
|
||||
"3. do in-place integration, save values to target"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"integrate! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function integrate(f::Function, el::JuliaFEM.Element)\n",
|
||||
" target = []\n",
|
||||
" for m = 1:length(el.iweights)\n",
|
||||
" w = el.iweights[m]\n",
|
||||
" xi = el.ipoints[:, m]\n",
|
||||
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
|
||||
" push!(target, w*f(el, xi)*det(J))\n",
|
||||
" end\n",
|
||||
" return sum(target)\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function integrate(f::Function)\n",
|
||||
" function integrate(el::JuliaFEM.Element)\n",
|
||||
" target = []\n",
|
||||
" for m = 1:length(el.iweights)\n",
|
||||
" w = el.iweights[m]\n",
|
||||
" xi = el.ipoints[:, m]\n",
|
||||
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
|
||||
" push!(target, w*f(el, xi)*det(J))\n",
|
||||
" end\n",
|
||||
" return sum(target)\n",
|
||||
" end\n",
|
||||
" return integrate\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function integrate!(f::Function, el::JuliaFEM.Element, target)\n",
|
||||
" # set target to zero\n",
|
||||
" el.attributes[target][:] = 0.0\n",
|
||||
" for m = 1:length(el.iweights)\n",
|
||||
" w = el.iweights[m]\n",
|
||||
" xi = el.ipoints[:, m]\n",
|
||||
" J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n",
|
||||
" el.attributes[target][:,:] += w*f(el, xi)*det(J)\n",
|
||||
" end\n",
|
||||
"end\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x4 Array{Float64,2}:\n",
|
||||
" -38.2303 -72.8697 79.4912 31.6088\n",
|
||||
" -17.625 -28.475 37.7 8.4 "
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"integrate(calc_residual_vector_integrand, e)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x4 Array{Float64,2}:\n",
|
||||
" -38.2303 -72.8697 79.4912 31.6088\n",
|
||||
" -17.625 -28.475 37.7 8.4 "
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"calc_residual_vector = integrate(calc_residual_vector_integrand)\n",
|
||||
"calc_residual_vector(e)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x4 Array{Float64,2}:\n",
|
||||
" -38.2303 -72.8697 79.4912 31.6088\n",
|
||||
" -17.625 -28.475 37.7 8.4 "
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"integrate!(calc_residual_vector_integrand, e, \"displacement nodal force\")\n",
|
||||
"e.attributes[\"displacement nodal force\"]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Linearization\n",
|
||||
"\n",
|
||||
"1. take function, element and field, and return partial derivative\n",
|
||||
"2. take function and field, return function which takes element as argument\n",
|
||||
"3. do in-place linearization to target, requires function which takes element as argument\n",
|
||||
"\n",
|
||||
"In general linearization can be done before integration and vice versa, i.e.\n",
|
||||
"\n",
|
||||
" integrate(linearize(f, \"u\"))(e) <-> linearize(integrate(f), \"u\")(e)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"linearize! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function linearize(f::Function, el::JuliaFEM.Element, field::ASCIIString)\n",
|
||||
" dim, nnodes = size(el.attributes[field])\n",
|
||||
" function helper!(x, y)\n",
|
||||
" orig = copy(el.attributes[field])\n",
|
||||
" el.attributes[field] = reshape(x, dim, nnodes)\n",
|
||||
" y[:] = f(el)\n",
|
||||
" el.attributes[field] = copy(orig)\n",
|
||||
" end\n",
|
||||
" jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
|
||||
" return jac(el.attributes[field][:])\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function linearize(f::Function, field::ASCIIString)\n",
|
||||
" function jacobian(el::JuliaFEM.Element, args...)\n",
|
||||
" dim, nnodes = size(el.attributes[field])\n",
|
||||
" function helper!(x, y)\n",
|
||||
" orig = copy(el.attributes[field])\n",
|
||||
" el.attributes[field] = reshape(x, dim, nnodes)\n",
|
||||
" y[:] = f(el, args...)\n",
|
||||
" el.attributes[field] = copy(orig)\n",
|
||||
" end\n",
|
||||
" jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
|
||||
" return jac(el.attributes[field][:])\n",
|
||||
" end\n",
|
||||
" return jacobian\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function linearize!(f::Function, el::JuliaFEM.Element, field::ASCIIString, target::ASCIIString)\n",
|
||||
" el.attributes[target][:] = 0.0\n",
|
||||
" dim, nnodes = size(el.attributes[field])\n",
|
||||
" function helper!(x, y)\n",
|
||||
" orig = copy(el.attributes[field])\n",
|
||||
" el.attributes[field] = reshape(x, dim, nnodes)\n",
|
||||
" y[:] = f(el)\n",
|
||||
" el.attributes[field] = copy(orig)\n",
|
||||
" end\n",
|
||||
" jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
|
||||
" jac!(el.attributes[field][:], el.attributes[target])\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
|
||||
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
|
||||
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
|
||||
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
|
||||
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
|
||||
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
|
||||
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
|
||||
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"integrate(linearize(calc_residual_vector_integrand, \"displacement\"))(e)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
|
||||
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
|
||||
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
|
||||
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
|
||||
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
|
||||
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
|
||||
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
|
||||
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
|
||||
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
|
||||
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
|
||||
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
|
||||
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
|
||||
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
|
||||
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
|
||||
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"linearize(integrate(calc_residual_vector_integrand), e, \"displacement\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
|
||||
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
|
||||
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
|
||||
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
|
||||
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
|
||||
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
|
||||
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
|
||||
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"linearize!(integrate(calc_residual_vector_integrand), e, \"displacement\", \"displacement tangent stiffness\")\n",
|
||||
"e.attributes[\"displacement tangent stiffness\"]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n",
|
||||
" 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n",
|
||||
" 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n",
|
||||
" 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n",
|
||||
" -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n",
|
||||
" -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n",
|
||||
" -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n",
|
||||
" -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"integrate!(linearize(calc_residual_vector_integrand, \"displacement\"), e, \"displacement tangent stiffness\")\n",
|
||||
"e.attributes[\"displacement tangent stiffness\"]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Validations"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Converged in 6 iterations.\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": 14,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"free_dofs = [3, 4, 5, 6]\n",
|
||||
"u = zeros(2, 4)\n",
|
||||
"du = zeros(2, 4)\n",
|
||||
"F = [0 0; 0 0; 0 -2; 0 0]'\n",
|
||||
"for i=1:10\n",
|
||||
" e.attributes[\"displacement\"] = u\n",
|
||||
" K = linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)\n",
|
||||
" R = integrate(calc_residual_vector_integrand)(e)\n",
|
||||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
|
||||
" u += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged in $i iterations.\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"u # -2.222244754401764"
|
||||
]
|
||||
}
|
||||
],
|
||||
"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
|
||||
}
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
Reference in New Issue
Block a user