mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
nice progress, still unfinished
This commit is contained in:
Jukka Aho
@@ -27,7 +27,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -38,7 +38,7 @@
|
||||
"Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -63,7 +63,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -74,7 +74,7 @@
|
||||
"calc_local_matrices! (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -159,12 +159,12 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"*Design principle 6*: we test our code. We use FactCheck in testing."
|
||||
"*Design principle 6*: we test our code. We use FactCheck for testing."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
@@ -175,7 +175,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -191,16 +191,23 @@
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"28-Jul 22:15:35:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"28-Jul 22:15:35:DEBUG:root:solution vector: \n",
|
||||
"29-Jul 23:12:24:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"29-Jul 23:12:24:DEBUG:root:solution vector: \n",
|
||||
" [0.0 -0.39914506095474317 -0.07228582695592449 0.0\n",
|
||||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
|
||||
"28-Jul 22:15:35:DEBUG:root:norm of u: 3.1292483947150043\n",
|
||||
"28-Jul 22:15:35:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"28-Jul 22:15:35:DEBUG:root:solution vector: \n",
|
||||
"29-Jul 23:12:24:DEBUG:root:norm of u: 3.1292483947150043\n",
|
||||
"29-Jul 23:12:25:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"29-Jul 23:12:25:DEBUG:root:solution vector: \n",
|
||||
" [0.0 0.7433248532717793 1.0485210147234858 0.0\n",
|
||||
" 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n",
|
||||
"28-Jul 22:15:35:DEBUG:root:norm of u: 3.129248394715006\n"
|
||||
"29-Jul 23:12:25:DEBUG:root:norm of u: 3.129248394715006\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"2 facts verified.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -209,7 +216,7 @@
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -289,56 +296,63 @@
|
||||
"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."
|
||||
"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": 12,
|
||||
"execution_count": 5,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"2 facts verified.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"type Node\n",
|
||||
" id :: Int\n",
|
||||
" coordinates :: Array{Float64, 1}\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"type Element\n",
|
||||
" id :: Int\n",
|
||||
" nodes :: Array{Node, 1}\n",
|
||||
" attributes :: Dict{ASCIIString, Any}\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"type FunctionSpace\n",
|
||||
" family :: String\n",
|
||||
" order :: Int\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"type IntegrationScheme\n",
|
||||
" points :: Array{Float64, 2}\n",
|
||||
" weights :: Array{Float64, 1}\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"type IJVMatrix\n",
|
||||
" I :: Array{Int64, 1}\n",
|
||||
" J :: Array{Int64, 1}\n",
|
||||
" V :: Array{Float64, 1}\n",
|
||||
" elements :: Array{Int64, 1}\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 105,
|
||||
"execution_count": 6,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"type Element\n",
|
||||
" id :: Int\n",
|
||||
" nodes :: 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": [],
|
||||
"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": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -346,17 +360,29 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Element2(1,[Node2(1,[1.0,2.0,3.0])],Dict(\"Young's modulus\"=>2.1e11,\"Poisson's ration\"=>0.3),\"Lagrange(1,2)\",\"FPG4\")"
|
||||
"get_shape_functions (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 105,
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function get_shape_functions(fs::FunctionSpace)\n",
|
||||
" if (fs.family == \"Lagrange\") & (fs.order == 1)\n",
|
||||
"\"\"\"\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",
|
||||
" 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",
|
||||
@@ -369,36 +395,12 @@
|
||||
" return basis, dbasis\n",
|
||||
" end\n",
|
||||
" throw(\"Unknown function space\")\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function get_integration_scheme(is::IntegrationScheme)\n",
|
||||
" if is.\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function assemble!(IJV::IJVMatrix, element::Element, gdofs, fs::FunctionSpace, is::IntegrationScheme)\n",
|
||||
" nnodes = length(element.nodes)\n",
|
||||
" ndim = length(elements.nodes[1].coordinates)\n",
|
||||
"\n",
|
||||
" # Material properties\n",
|
||||
" E = element.attributes[\"Young\"]\n",
|
||||
" nu = element.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",
|
||||
" R = zeros(ndim, nnodes)\n",
|
||||
" K = zeros(ndim*nnodes, ndim*nnodes)\n",
|
||||
"\n",
|
||||
" basis, dbasis = get_shape_functions(fs)\n",
|
||||
"\n",
|
||||
" calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
|
||||
"\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -406,25 +408,30 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"IntegrationScheme(4x2 Array{Float64,2}:\n",
|
||||
" -0.57735 -0.57735\n",
|
||||
" 0.57735 -0.57735\n",
|
||||
" 0.57735 0.57735\n",
|
||||
" -0.57735 0.57735,[1.0,1.0,1.0,1.0])"
|
||||
"get_integration_scheme (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"is = IntegrationScheme(1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1], [1, 1, 1, 1])\n"
|
||||
"\"\"\"\n",
|
||||
"\"\"\"\n",
|
||||
"function get_integration_scheme(el::Element, order=2)\n",
|
||||
" ndim, nnodes = size(el.coordinates)\n",
|
||||
" if (nnodes == 4) & (order == 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",
|
||||
" end\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 111,
|
||||
"execution_count": 48,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -432,45 +439,183 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Node2(2,[3.0,2.0,3.0])"
|
||||
"assemble_element! (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 111,
|
||||
"execution_count": 48,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"e2 = Element2(2, [n1], Dict(\"Young's modulus\" => 90.0e9, \"Poisson's ration\" => 0.3), \"Lagrange(1,2)\", \"FPG4\")\n",
|
||||
"e2.nodes[1]"
|
||||
"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",
|
||||
" 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": null,
|
||||
"execution_count": 49,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"one element assembly\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"30-Jul 00:08:18:DEBUG:root:Creating nodes\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Creating elements\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 1\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 3.0814821107320176\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 2\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 0.32007464366194766\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 3\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 0.040279810888447135\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 4\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 0.000925649536063315\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 5\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 1.5582555024825817e-7\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Starting iteration 6\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Assembling\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solving\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Solution norm = 1.0317166868587156e-14\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"30-Jul 00:08:18:DEBUG:root:Displacement of element = \n",
|
||||
"[-0.39914506095474334 -0.07228582695592464 0.0 0.0\n",
|
||||
" -2.1779892317073504 -2.222244754401764 0.0 0.0]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 49,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"1 fact verified.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"n1 = new_node()\n",
|
||||
"set_node_id(n1, 1)\n",
|
||||
"set_node_coords(n1, [0.0, 0.0])\n",
|
||||
"facts(\"one element assembly\") do\n",
|
||||
" Logging.debug(\"Creating nodes\")\n",
|
||||
" n1 = Node(1, Int64[])\n",
|
||||
" n2 = Node(2, Int64[])\n",
|
||||
" n3 = Node(3, Int64[])\n",
|
||||
" n4 = Node(4, Int64[])\n",
|
||||
" nodes = [n1.id, n2.id, n3.id, n4.id]\n",
|
||||
" coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n",
|
||||
" attributes = Dict(\n",
|
||||
" \"Young\" => 90, \"Poisson\" => 0.25,\n",
|
||||
" \"displacement\" => zeros(2, 4),\n",
|
||||
" \"displacement nodal force\" => zeros(2, 4),\n",
|
||||
" \"displacement tangent stiffness\" => zeros(8, 8))\n",
|
||||
" Logging.debug(\"Creating elements\")\n",
|
||||
" el = Element(1, nodes, coordinates, attributes)\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",
|
||||
" Logging.debug(\"Solving\")\n",
|
||||
" du = zeros(2, 4) # must be determined from ass\n",
|
||||
" F = [0 0; 0 -2; 0 0; 0 0]'\n",
|
||||
" free_dofs = [1, 2, 3, 4]\n",
|
||||
"\n",
|
||||
"n2 = new_node()\n",
|
||||
"set_node_id(n2, 2)\n",
|
||||
"set_node_coords(n2, [10.0, 0.0])\n",
|
||||
" # solution\n",
|
||||
" K = sparse(ass.I, ass.J, ass.A)\n",
|
||||
" R = full(sparsevec(ass.i, ass.b))\n",
|
||||
" R = reshape(R, (2, round(Int, length(R)/2)))\n",
|
||||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
|
||||
"\n",
|
||||
"n3 = new_node()\n",
|
||||
"set_node_id(n3, 3)\n",
|
||||
"set_node_coords(n3, [10.0, 1.0])\n",
|
||||
" Logging.debug(\"Solution norm = $(norm(du))\")\n",
|
||||
" #Logging.debug(\"Solution increment = \\n$du\")\n",
|
||||
"\n",
|
||||
"n4 = new_node()\n",
|
||||
"set_node_id(n4, 4)\n",
|
||||
"set_node_coords(n4, [0.0, 1.0])\n",
|
||||
"\n",
|
||||
"nodes = [n1, n2, n3, n4]\n",
|
||||
"add_nodes(m, \"NALL\", nodes)"
|
||||
" # 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)."
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
Reference in New Issue
Block a user