mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 09:54:55 +00:00
updated notebook
This commit is contained in:
Jukka Aho
@@ -191,16 +191,16 @@
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"29-Jul 23:12:24:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"29-Jul 23:12:24:DEBUG:root:solution vector: \n",
|
||||
"01-Aug 14:05:35:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"01-Aug 14:05:36:DEBUG:root:solution vector: \n",
|
||||
" [0.0 -0.39914506095474317 -0.07228582695592449 0.0\n",
|
||||
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\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",
|
||||
"01-Aug 14:05:36:DEBUG:root:norm of u: 3.1292483947150043\n",
|
||||
"01-Aug 14:05:37:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"01-Aug 14:05:37:DEBUG:root:solution vector: \n",
|
||||
" [0.0 0.7433248532717793 1.0485210147234858 0.0\n",
|
||||
" 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n",
|
||||
"29-Jul 23:12:25:DEBUG:root:norm of u: 3.129248394715006\n"
|
||||
"01-Aug 14:05:37:DEBUG:root:norm of u: 3.129248394715006\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -352,7 +352,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -363,7 +363,7 @@
|
||||
"get_shape_functions (generic function with 1 method)"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -400,7 +400,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -411,7 +411,7 @@
|
||||
"get_integration_scheme (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -431,7 +431,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 48,
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -442,7 +442,7 @@
|
||||
"assemble_element! (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 48,
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -494,7 +494,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 49,
|
||||
"execution_count": 11,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -510,58 +510,53 @@
|
||||
"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",
|
||||
"01-Aug 14:07:09:DEBUG:root:Creating nodes\n",
|
||||
"01-Aug 14:07:09:DEBUG:root:Creating elements\n",
|
||||
"01-Aug 14:07:09:DEBUG:root:Starting iteration 1\n",
|
||||
"01-Aug 14:07:09:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 3.090022136728999\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Starting iteration 2\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 0.3212131602153504\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Starting iteration 3\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 0.040431781940005365\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Starting iteration 4\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 0.0009291101052064257\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Starting iteration 5\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 1.5638899025378415e-7\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Starting iteration 6\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Assembling\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Solution norm = 1.0355045356935844e-14\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Converged in 6 iterations.\n",
|
||||
"01-Aug 14:07:10:DEBUG:root:Displacement of element = \n",
|
||||
"[-0.39914506095474334 -0.07228582695592464 0.0 0.0\n",
|
||||
" -2.1779892317073504 -2.222244754401764 0.0 0.0]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"1 fact verified.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 49,
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"1 fact verified.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"facts(\"one element assembly\") do\n",
|
||||
" # Create model\n",
|
||||
" Logging.debug(\"Creating nodes\")\n",
|
||||
" n1 = Node(1, Int64[])\n",
|
||||
" n2 = Node(2, Int64[])\n",
|
||||
@@ -569,32 +564,35 @@
|
||||
" 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",
|
||||
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||||
" Logging.debug(\"Creating elements\")\n",
|
||||
" el = Element(1, nodes, coordinates, attributes)\n",
|
||||
"\n",
|
||||
" # Initialize elements ready for solution\n",
|
||||
" el.attributes[\"displacement\"] = zeros(2, 4)\n",
|
||||
" el.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
|
||||
" el.attributes[\"displacement tangent stiffness\"] = zeros(8, 8)\n",
|
||||
"\n",
|
||||
" for i=1:10\n",
|
||||
" Logging.debug(\"Starting iteration $i\")\n",
|
||||
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())\n",
|
||||
" ass.gdofs[el.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n",
|
||||
" Logging.debug(\"Assembling\")\n",
|
||||
" assemble_element!(ass, el)\n",
|
||||
" Logging.debug(\"Solving\")\n",
|
||||
" du = zeros(2, 4) # must be determined from ass\n",
|
||||
"\n",
|
||||
" # Boundary conditions\n",
|
||||
" F = [0 0; 0 -2; 0 0; 0 0]'\n",
|
||||
" F = reshape(F, prod(size(F)))\n",
|
||||
" free_dofs = [1, 2, 3, 4]\n",
|
||||
"\n",
|
||||
" # solution\n",
|
||||
" K = sparse(ass.I, ass.J, ass.A)\n",
|
||||
" R = full(sparsevec(ass.i, ass.b))\n",
|
||||
" R = reshape(R, (2, round(Int, length(R)/2)))\n",
|
||||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
|
||||
" R = R - F\n",
|
||||
" du = zeros(8) # must be determined from ass\n",
|
||||
" du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n",
|
||||
"\n",
|
||||
" Logging.debug(\"Solution norm = $(norm(du))\")\n",
|
||||
" #Logging.debug(\"Solution increment = \\n$du\")\n",
|
||||
"\n",
|
||||
" # update solution back to elements\n",
|
||||
" eldu = du[ass.gdofs[el.id]]\n",
|
||||
@@ -615,269 +613,253 @@
|
||||
"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)."
|
||||
"Seems to be working. But we still need to handle boundary conditions more \"cleverly\" and generalize assembly to several elements (which is not problem).\n",
|
||||
"\n",
|
||||
"First of all, essential boundary conditions are nothing more than equality constraints saying that value for some degree of freedom is fixed. Elimination is just a special case when this value equals to zero. There is couple of different strategies to handle essential boundary conditions. One option is to force them using Lagrange multipliers which can also be used to create all kind of kinematic constraints also. (For example, contact can be considered as a kinematic constraint.) Another option is to manipulate matrix such a way that constraint is satisfied.\n",
|
||||
"\n",
|
||||
"Because we are now going \"bottom-up\", we develop something extremely simple that however deals with the problem:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"type BC\n",
|
||||
" dofs :: Array{Int64, 1}\n",
|
||||
" values :: Array{Float64, 1}\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 49,
|
||||
"metadata": {
|
||||
"collapsed": false,
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Converged\n"
|
||||
"solve one element problem\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"01-Aug 20:51:27:DEBUG:root:Creating nodes\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Creating elements\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Problem size = 8\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 1\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 3.0900221367289444\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 2\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 0.32121316021534363\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 3\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 0.04043178194002483\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 4\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 0.0009291101052060104\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 5\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 1.563889898905983e-7\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Starting iteration 6\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Assembling\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Dirichlet boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 5 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 6 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 7 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:dof 8 => 0.0\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Added 4 Lagrange multipliers to model\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Adding Neumann boundary conditions\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solving system of equations. Total size = 12\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Solution norm = 1.2782462771683917e-14\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Converged in 4 iterations.\n",
|
||||
"01-Aug 20:51:27:DEBUG:root:Displacement of element = \n",
|
||||
"[-0.3991450609547439 -0.07228582695592495 0.0 0.0\n",
|
||||
" -2.1779892317073513 -2.2222447544017654 0.0 0.0]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"1 fact verified.\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"
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"execution_count": 49,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"#(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
"# la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()\n",
|
||||
"function solve_one_element()\n",
|
||||
"function solve!(elements, neumann_bcs, dirichlet_bcs; dofs=2, max_iterations=10)\n",
|
||||
"\n",
|
||||
" # create nodes separately ...\n",
|
||||
" # Initialize elements ready for solution\n",
|
||||
" for el in elements\n",
|
||||
" el.attributes[\"displacement\"] = zeros(2, 4)\n",
|
||||
" el.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
|
||||
" el.attributes[\"displacement tangent stiffness\"] = zeros(8, 8)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # .. or create somewhat simpler syntax\n",
|
||||
" # nodes = Dict(1 => [0.0, 0.0], 2 => [10.0, 0.0], 3 => [10.0, 1.0], 4 => [0.0, 1.0])\n",
|
||||
" # add_nodes(m, \"NALL\", nodes)\n",
|
||||
" # Assign global dofs for elements\n",
|
||||
" gdofs = Dict{Int64,Array{Int64,1}}()\n",
|
||||
" pdim = 1\n",
|
||||
" for el in elements\n",
|
||||
" edofs = length(el.nodes)*dofs-1\n",
|
||||
" gdofs[el.id] = pdim:pdim+edofs\n",
|
||||
" pdim += edofs\n",
|
||||
" end\n",
|
||||
" Logging.debug(\"Problem size = $pdim\")\n",
|
||||
"\n",
|
||||
" # create element (hard way)\n",
|
||||
" e = new_element()\n",
|
||||
" set_element_id(1)\n",
|
||||
" set_node_ids(e, [1, 2, 3, 4])\n",
|
||||
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n",
|
||||
"\n",
|
||||
" # we can set function spaces and integration schema for element-wise ...\n",
|
||||
" set_function_space(e, \"Lagrange(1)\") # use linear Lagrange function space to approximate unknown field\n",
|
||||
" set_integration_schema(e, \"FPG4\") # use four integration points\n",
|
||||
" # or for model as a \"default value\"\n",
|
||||
" # set_function_space(m, \"Lagrange\", 1)\n",
|
||||
" # set_integration_schema(m, \"FPG4\")\n",
|
||||
"\n",
|
||||
" # set necessary material parameters\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",
|
||||
" # we can add attribute for each element node ...\n",
|
||||
" #add_attribute(e, \"lambda\", [la, la, la, la])\n",
|
||||
" #add_attribute(e, \"mu\", [mu, mu, mu, mu])\n",
|
||||
" # or just one for each element ...\n",
|
||||
" add_attribute(e, \"lambda\", lambda)\n",
|
||||
" add_attribute(e, \"mu\", mu)\n",
|
||||
" # .. or just for model as a default value\n",
|
||||
" # add_attribute(m, \"lambda\", lambda)\n",
|
||||
" # add_attribute(m, \"mu\", mu)\n",
|
||||
" # note that we don't assign attributes to nodes here so we can describe discontinous fields\n",
|
||||
" # in attributes\n",
|
||||
"\n",
|
||||
" # or we can use convenient syntax\n",
|
||||
" # e = new_element(element_id=1, node_ids=[1, 2, 3, 4], function_space=\"Lagrange(1)\",\n",
|
||||
" # integration_schema=\"FPG4\", attributes=Dict(\"lambda\" => lambda, \"mu\" => mu))\n",
|
||||
"\n",
|
||||
" elements = [e]\n",
|
||||
" add_elements(m, \"EALL\", elements)\n",
|
||||
"\n",
|
||||
" # boundary conditions are always assignet to sets\n",
|
||||
"\n",
|
||||
" # element boundary conditions\n",
|
||||
" \n",
|
||||
" # add load for element surface S1 in normal-tangential coordinate system\n",
|
||||
" # add_attribute(e, \"displacement S1 normal load\", 1)\n",
|
||||
" # add load for element surface S1 from -1 to -2\n",
|
||||
" # add_attribute(e, \"displacement S1 load\", [-1, -2])\n",
|
||||
"\n",
|
||||
" # nodal boundary conditions\n",
|
||||
" # add neumann boundary condition to node 3, in y direction\n",
|
||||
" add_attribute(n3, \"displacement 2 load\", -2)\n",
|
||||
" # add dirichlet boundary condition to nodes 1 and 2 (encastre)\n",
|
||||
" add_attribute(n1, \"displacement\", [0.0, 0.0])\n",
|
||||
" # or\n",
|
||||
" # add_attribute(n1, \"displacement 1\", 0.0)\n",
|
||||
" # add_attribute(n1, \"displacement 2\", 0.0)\n",
|
||||
" add_attribute(n2, \"displacement\", [0.0, 0.0])\n",
|
||||
" # or, for example, add dirichlet boundary conditions in directions 1 and 3 for node\n",
|
||||
" #add_attribute(n2, \"displacement 1,3\", [0.0, 1.0])\n",
|
||||
" # or fix all dofs of a nodes in nodeset \"SUPPORT\"\n",
|
||||
" # support_bc = add_nodeset(m, \"SUPPORT\", [n1, n2])\n",
|
||||
" # add_attribute(support_bc, \"displacement\", 0)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" # Everything is very general so far. Datamodel is well defined and in this point\n",
|
||||
" # we can save or load it to disk\n",
|
||||
"\n",
|
||||
" # save_model(m, \"mymodel\") # save model to xml/h5 (Xdmf)\n",
|
||||
" # m = load_model(\"mymodel\") # load model from file\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # END OF MODEL DEFINITON\n",
|
||||
" \n",
|
||||
" # Now we kick in elasticity iterations and solve displacement\n",
|
||||
" # field but of course it could be something else too\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" # m1, m2 = make_domain_decomposition(parts=2, keep_in_one_domain=[\"CONTACT_BOUNDARY\"])\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # In newton iteration, there might be situations where we just want to update RHS\n",
|
||||
" # like in radiation problems, it makes no sense to update stiffness matrix in that\n",
|
||||
" # case. For this reason local element matrix and force vector can be updated both\n",
|
||||
" # or just one of them. This time we don't have anything nonlinear in RHS so we can\n",
|
||||
" # calculate it outside of iteration loop\n",
|
||||
" p = new_problem(\"displacement\")\n",
|
||||
" \n",
|
||||
" RHS = zeros(length(nodes)*dofs)\n",
|
||||
" \n",
|
||||
" for i=1:10\n",
|
||||
" for i=1:max_iterations\n",
|
||||
" Logging.debug(\"Starting iteration $i\")\n",
|
||||
" ass.I = []\n",
|
||||
" ass.J = []\n",
|
||||
" ass.A = []\n",
|
||||
" ass.i = []\n",
|
||||
" ass.b = []\n",
|
||||
" \n",
|
||||
" solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
" la, mu, N, dNdξ, ipoints, iweights)\n",
|
||||
" u += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged\")\n",
|
||||
" Logging.debug(\"Assembling\")\n",
|
||||
" for el in elements\n",
|
||||
" assemble_element!(ass, el)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" Logging.debug(\"Adding Dirichlet boundary conditions\")\n",
|
||||
" # Dirichlet boundary conditions\n",
|
||||
" i = 1\n",
|
||||
" for bc in dirichlet_bcs\n",
|
||||
" for (dof, val) in zip(bc.dofs, bc.values)\n",
|
||||
" Logging.debug(\"dof $dof => $val\")\n",
|
||||
" push!(ass.I, dof)\n",
|
||||
" push!(ass.J, pdim+i)\n",
|
||||
" push!(ass.A, 1)\n",
|
||||
" push!(ass.I, pdim+i)\n",
|
||||
" push!(ass.J, dof)\n",
|
||||
" push!(ass.A, 1)\n",
|
||||
" push!(ass.i, pdim+i)\n",
|
||||
" push!(ass.b, 0)\n",
|
||||
" i += 1\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" i -= 1\n",
|
||||
" Logging.debug(\"Added $i Lagrange multipliers to model\")\n",
|
||||
"\n",
|
||||
" Logging.debug(\"Adding Neumann boundary conditions\")\n",
|
||||
" F = zeros(pdim+i)\n",
|
||||
" # Neumann boundary conditions\n",
|
||||
" for bc in neumann_bcs\n",
|
||||
" for (dof, val) in zip(bc.dofs, bc.values)\n",
|
||||
" F[dof] += val\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" Logging.debug(\"Solving system of equations. Total size = $(pdim+i)\")\n",
|
||||
" # solution\n",
|
||||
" K = sparse(ass.I, ass.J, ass.A)\n",
|
||||
" R = full(sparsevec(ass.i, ass.b))\n",
|
||||
" R = R - F\n",
|
||||
" #println(full(K))\n",
|
||||
" #println(R)\n",
|
||||
"\n",
|
||||
" du = K \\ -R\n",
|
||||
"\n",
|
||||
" solnorm = norm(du[1:pdim])\n",
|
||||
" Logging.debug(\"Solution norm = $solnorm\")\n",
|
||||
"\n",
|
||||
" # update solution back to elements\n",
|
||||
" for el in elements\n",
|
||||
" eldu = du[ass.gdofs[el.id]]\n",
|
||||
" eldu = reshape(eldu, (2, round(Int, length(eldu)/2)))\n",
|
||||
" el.attributes[\"displacement\"] += eldu\n",
|
||||
" end\n",
|
||||
" if solnorm < 1.0e-9\n",
|
||||
" Logging.debug(\"Converged in $i iterations.\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" return u\n",
|
||||
"\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"u"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"\n",
|
||||
"#### Parameters\n",
|
||||
"X : Array{Float64, 2}\n",
|
||||
" \n",
|
||||
"\"\"\"\n",
|
||||
"function solve_elasticity_increment!(X, u, du, elmap, nodalloads,\n",
|
||||
" dirichletbc, lambda, mu, N, dNdchi, ipoints,\n",
|
||||
" iweights)\n",
|
||||
" if length(size(elmap)) == 1\n",
|
||||
" # quick hack for just one element\n",
|
||||
" elmap = elmap''\n",
|
||||
" end\n",
|
||||
" nelnodes, nelements = size(elmap)\n",
|
||||
" dim, nnodes = size(u)\n",
|
||||
" dofs = dim*nelnodes\n",
|
||||
"\n",
|
||||
" Imat = Int64[]\n",
|
||||
" Jmat = Int64[]\n",
|
||||
" Vmat = Float64[]\n",
|
||||
" Ivec = Int64[]\n",
|
||||
" Vvec = Float64[]\n",
|
||||
"\n",
|
||||
" # FIXME: different number of nodes/element\n",
|
||||
" R = zeros(dim, nelnodes)\n",
|
||||
" Kt = zeros(dofs, dofs)\n",
|
||||
"\n",
|
||||
" # this can be parallelized\n",
|
||||
" for i in 1:nelements\n",
|
||||
" eldofs = elmap[:,i]\n",
|
||||
" calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdchi,\n",
|
||||
" lambda[eldofs], mu[eldofs], ipoints, iweights)\n",
|
||||
" assemble!(Kt, eldofs, Imat, Jmat, Vmat)\n",
|
||||
" assemble!(R, eldofs, Ivec, Vvec)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # add additional neumann boundary conditions to force vector\n",
|
||||
" for (i, nodal_load) in enumerate(nodalloads)\n",
|
||||
" if nodal_load == 0\n",
|
||||
" continue\n",
|
||||
" end\n",
|
||||
" push!(Ivec, i)\n",
|
||||
" push!(Vvec, -nodal_load)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # Create sparse matrix and vector\n",
|
||||
" A = sparse(Imat, Jmat, Vmat)\n",
|
||||
" b = sparsevec(Ivec, Vvec)\n",
|
||||
"\n",
|
||||
" # Remove dirichlet boundary conditions\n",
|
||||
" free_dofs = find(isnan(dirichletbc))\n",
|
||||
" #Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat)\n",
|
||||
" b = b[free_dofs]\n",
|
||||
" A = A[free_dofs, free_dofs]\n",
|
||||
"\n",
|
||||
" # solution\n",
|
||||
" du[free_dofs] = lufact(A) \\ -full(b)\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function one_elem_fixture()\n",
|
||||
" X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
|
||||
" elmap = [1; 2; 3; 4]\n",
|
||||
" nodalloads = [0 0; 0 0; 0 -2; 0 0]'\n",
|
||||
" dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0]'\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",
|
||||
" la = la*ones(1, 4)\n",
|
||||
" mu = mu*ones(1, 4)\n",
|
||||
" u = zeros(2, 4)\n",
|
||||
" du = zeros(2, 4)\n",
|
||||
"\n",
|
||||
" N(xi) = [\n",
|
||||
" (1-xi[1])*(1-xi[2])/4\n",
|
||||
" (1+xi[1])*(1-xi[2])/4\n",
|
||||
" (1+xi[1])*(1+xi[2])/4\n",
|
||||
" (1-xi[1])*(1+xi[2])/4\n",
|
||||
" ]\n",
|
||||
"\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",
|
||||
" ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
|
||||
" iweights = [1 1 1 1]\n",
|
||||
"\n",
|
||||
" return (X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
" la, mu, N, dNdξ, ipoints, iweights)\n",
|
||||
"facts(\"solve one element problem\") do\n",
|
||||
" # Create model\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(\"Young\" => 90, \"Poisson\" => 0.25)\n",
|
||||
" Logging.debug(\"Creating elements\")\n",
|
||||
" el = Element(1, nodes, coordinates, attributes)\n",
|
||||
" elements = [el]\n",
|
||||
" # Boundary conditions\n",
|
||||
" bc1 = BC([4], [-2.0]) # force boundary condition, third dof -2\n",
|
||||
" bc2 = BC([5, 6, 7, 8], [0.0, 0.0, 0.0, 0.0]) # dirichlet bc, set dx=dy=0 on support\n",
|
||||
" solve!(elements, [bc1], [bc2]; max_iterations=7)\n",
|
||||
" disp = elements[1].attributes[\"displacement\"]\n",
|
||||
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
|
||||
" @fact norm(disp) => roughly(3.1292483947150043)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"u3d = [u; 0 0 0 0] # extend to 3d vector field\n",
|
||||
"X3d = [X; 0 0 0 0]\n",
|
||||
"elmap2 = [0x5; elmap]'';"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
|
||||
Reference in New Issue
Block a user