diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index dd8e6f5..e509b7b 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -16,6 +16,13 @@ "collapsed": false }, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Docile: upgrading cache to 0.0.2.\n" + ] + }, { "name": "stderr", "output_type": "stream", @@ -25,7 +32,10 @@ "Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n", "\n", "WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:32.\n", - "Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n" + "Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n", + "Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n", + "Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n", + "Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n" ] }, { @@ -50,7 +60,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -66,6 +76,28 @@ "## 2d beam with linear elements" ] }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Array{Float64,2}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "typeof(Float64[1 2; 2 3])" + ] + }, { "cell_type": "code", "execution_count": 2, @@ -94,6 +126,64 @@ } ], "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", @@ -130,17 +220,137 @@ " la, mu, N, dNdξ, ipoints, iweights)\n", "end\n", "\n", - "(X, u, du, elmap, nodalloads, dirichletbc,\n", - " la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()\n", + "#(X, u, du, elmap, nodalloads, dirichletbc,\n", + "# la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()\n", + "function solve_one_element()\n", + "\n", + " m = new_model()\n", + "\n", + " # create nodes separately ...\n", + " n1 = new_node()\n", + " set_node_id(n1, 1)\n", + " set_node_coords(n1, [0.0, 0.0])\n", + "\n", + " n2 = new_node()\n", + " set_node_id(n2, 2)\n", + " set_node_coords(n2, [10.0, 0.0])\n", " \n", - "for i=1:10\n", - " JuliaFEM.elasticity_solver.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", - " break\n", + " n3 = new_node()\n", + " set_node_id(n3, 3)\n", + " set_node_coords(n3, [10.0, 1.0])\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)\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", + "\n", + " # create element (hard way)\n", + " e = new_element()\n", + " set_element_id(1)\n", + " set_node_ids(e, [1, 2, 3, 4])\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", + " \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", + " break\n", + " end\n", " end\n", + " return u\n", "end\n", "\n", "u"