From 1fa7f259fdb4dcb609db34a094f8937b79b6906c Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 19 Aug 2015 18:50:43 +0300 Subject: [PATCH] updated. --- ...2015-06-25-elasticity-solver-example.ipynb | 667 +++++++++--------- src/abaqus_reader.jl | 6 +- src/math.jl | 30 +- src/types.jl | 21 +- 4 files changed, 354 insertions(+), 370 deletions(-) diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 2c4110f..92a818b 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -36,7 +36,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 1, "metadata": { "collapsed": false }, @@ -47,9 +47,30 @@ "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" ] }, - "execution_count": 2, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", + " in int32 at deprecated.jl:49\n", + " in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n", + " in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n", + " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n", + " in anonymous at task.jl:365\n", + "while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n", + "WARNING: int32(x) is deprecated, use Int32(x) instead.\n", + " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", + " in int32 at deprecated.jl:49\n", + " in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n", + " in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n", + " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n", + " in anonymous at task.jl:365\n", + "while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n" + ] } ], "source": [ @@ -73,7 +94,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -84,39 +105,35 @@ "jacobian (generic function with 1 method)" ] }, - "execution_count": 3, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n", - " # Calculate dN/dX\n", - " dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n", - "\n", - " # Calculate residual force vector R(u) = T(u) - F(u)\n", - " function T(u)\n", - " # kinematics\n", - " gradu = u*dbasisdX\n", - " F = I + gradu\n", - " E = 1/2*(gradu' + gradu + gradu'*gradu)\n", - " # constitutive equation\n", - " lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n", - " mu = JuliaFEM.interpolate(el, \"mu\", xi)\n", - " S = lambda*trace(E)*I + 2*mu*E\n", - " P = F*S\n", - " T = P*dbasisdX'\n", - " return T\n", - " end\n", - "\n", + "function Wint_integrand(el, ip)\n", + " xi = ip.xi\n", + " dbasisdX = JuliaFEM.get_dbasisdX(el, ip)\n", " u = el.attributes[\"displacement\"]\n", - " R(u) = T(u)\n", - " return R(u) # we skip calculating external force vector for now.\n", + "\n", + " # kinematics\n", + " gradu = u*dbasisdX\n", + " F = I + gradu\n", + " E = 1/2*(gradu' + gradu + gradu'*gradu)\n", + "\n", + " # constitutive equation\n", + " lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n", + " mu = JuliaFEM.interpolate(el, \"mu\", xi)\n", + " S = lambda*trace(E)*I + 2*mu*E\n", + " P = F*S\n", + " T = P*dbasisdX'\n", + " return T\n", "end\n", "\n", - "calc_residual_vector = JuliaFEM.integrate(calc_residual_vector_integrand)\n", + "Wint = JuliaFEM.integrate(Wint_integrand)\n", + "\n", "# calculate partial derivatives of R with respect to field \"displacement\"\n", - "calc_tangent_stiffness = JuliaFEM.linearize(calc_residual_vector, \"displacement\")" + "jacobian = JuliaFEM.linearize(Wint, \"displacement\")" ] }, { @@ -130,7 +147,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 3, "metadata": { "collapsed": true }, @@ -141,7 +158,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 4, "metadata": { "collapsed": false, "scrolled": false @@ -158,16 +175,16 @@ "name": "stderr", "output_type": "stream", "text": [ - "17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n", - "17-Aug 18:21:25:DEBUG:root:solution vector: \n", + "19-Aug 18:39:44:DEBUG:root:Converged in 6 iterations.\n", + "19-Aug 18:39:44:DEBUG:root:solution vector: \n", " [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n", - "17-Aug 18:21:25:DEBUG:root:norm of u: 3.1292483947150047\n", - "17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n", - "17-Aug 18:21:25:DEBUG:root:solution vector: \n", + "19-Aug 18:39:44:DEBUG:root:norm of u: 3.1292483947150047\n", + "19-Aug 18:39:45:DEBUG:root:Converged in 6 iterations.\n", + "19-Aug 18:39:45:DEBUG:root:solution vector: \n", " [0.0 0.7433248532717796 1.048521014723486 0.0\n", " 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n", - "17-Aug 18:21:26:DEBUG:root:norm of u: 3.1292483947150056\n" + "19-Aug 18:39:45:DEBUG:root:norm of u: 3.1292483947150056\n" ] }, { @@ -183,7 +200,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 6, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -200,10 +217,15 @@ " (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.0, 1.0, 1.0, 1.0]\n", + " integration_points = [\n", + " JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n", + " JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n", + " JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0, Dict()),\n", + " JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0, Dict())]\n", " attributes = Dict()\n", - " e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\n", + " element_id = 1\n", + " node_ids = [1, 2, 3, 4]\n", + " e = JuliaFEM.Element(element_id, node_ids, basis, dbasis, integration_points, attributes)\n", " E = 90.0\n", " nu = 0.25\n", " mu = E/(2*(1+nu))\n", @@ -225,8 +247,8 @@ "\n", " free_dofs = [3, 4, 5, 6]\n", " for i=1:10\n", - " R = calc_residual_vector(e)\n", - " K = calc_tangent_stiffness(e)\n", + " R = Wint(e)\n", + " K = jacobian(e)\n", " du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n", " e.attributes[\"displacement\"] += du\n", " if norm(du) < 1.0e-9\n", @@ -253,8 +275,8 @@ " e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n", " du = zeros(2, 4)\n", " for i=1:10\n", - " R = calc_residual_vector(e)\n", - " K = calc_tangent_stiffness(e)\n", + " R = Wint(e)\n", + " K = jacobian(e)\n", " du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n", " e.attributes[\"displacement\"] += du\n", " if norm(du) < 1.0e-9\n", @@ -269,16 +291,9 @@ "end" ] }, - { - "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. Some data types:" - ] - }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -286,35 +301,21 @@ { "data": { "text/plain": [ - "assemble_element! (generic function with 2 methods)" + "assemble_element! (generic function with 1 method)" ] }, - "execution_count": 9, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function assemble_element!(ass::Assembly, el::Element)\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", + "function assemble_element!(ass::JuliaFEM.Assembly, el::JuliaFEM.Element)\n", "\n", " gdofs = ass.gdofs[el.id]\n", - " #Logging.debug(\"Assemble element $(el.id) to gdofs $gdofs\")\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", + " R = Wint(el)\n", + " K = jacobian(el)\n", + " dofs = length(R)\n", "\n", " for i=1:dofs\n", " for j=1:dofs\n", @@ -342,18 +343,11 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING: the `=>` syntax is deprecated, use `-->` instead\n" - ] - }, { "name": "stdout", "output_type": "stream", @@ -365,30 +359,29 @@ "name": "stderr", "output_type": "stream", "text": [ - "15-Aug 17:03:55:DEBUG:root:Adding nodes to array\n", - "15-Aug 17:03:55:DEBUG:root:Creating elements\n", - "15-Aug 17:03:55:DEBUG:root:Starting iteration 1\n", - "15-Aug 17:03:55:DEBUG:root:Assembling\n", - "15-Aug 17:03:57:DEBUG:root:Solution norm = 3.090022136728999\n", - "15-Aug 17:03:57:DEBUG:root:Starting iteration 2\n", - "15-Aug 17:03:57:DEBUG:root:Assembling\n", - "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.32121316021535135\n", - "15-Aug 17:03:57:DEBUG:root:Starting iteration 3\n", - "15-Aug 17:03:57:DEBUG:root:Assembling\n", - "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.040431781939994194\n", - "15-Aug 17:03:57:DEBUG:root:Starting iteration 4\n", - "15-Aug 17:03:57:DEBUG:root:Assembling\n", - "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.0009291101052124042\n", - "15-Aug 17:03:57:DEBUG:root:Starting iteration 5\n", - "15-Aug 17:03:57:DEBUG:root:Assembling\n", - "15-Aug 17:03:57:DEBUG:root:Solution norm = 1.5638899022213175e-7\n", - "15-Aug 17:03:58:DEBUG:root:Starting iteration 6\n", - "15-Aug 17:03:58:DEBUG:root:Assembling\n", - "15-Aug 17:03:58:DEBUG:root:Solution norm = 1.0118539067290854e-14\n", - "15-Aug 17:03:58:DEBUG:root:Converged in 6 iterations.\n", - "15-Aug 17:03:58:DEBUG:root:Displacement of element = \n", - "[-0.39914506095474334 -0.0722858269559246 0.0 0.0\n", - " -2.1779892317073504 -2.2222447544017645 0.0 0.0]\n" + "19-Aug 18:39:46:DEBUG:root:Adding nodes to array\n", + "19-Aug 18:39:46:DEBUG:root:Starting iteration 1\n", + "19-Aug 18:39:46:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 3.0900221367289986\n", + "19-Aug 18:39:47:DEBUG:root:Starting iteration 2\n", + "19-Aug 18:39:47:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 0.3212131602153472\n", + "19-Aug 18:39:47:DEBUG:root:Starting iteration 3\n", + "19-Aug 18:39:47:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 0.04043178193999703\n", + "19-Aug 18:39:47:DEBUG:root:Starting iteration 4\n", + "19-Aug 18:39:47:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 0.0009291101052105739\n", + "19-Aug 18:39:47:DEBUG:root:Starting iteration 5\n", + "19-Aug 18:39:47:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 1.5638899027804743e-7\n", + "19-Aug 18:39:47:DEBUG:root:Starting iteration 6\n", + "19-Aug 18:39:47:DEBUG:root:Assembling\n", + "19-Aug 18:39:47:DEBUG:root:Solution norm = 1.0464940956129567e-14\n", + "19-Aug 18:39:47:DEBUG:root:Converged in 6 iterations.\n", + "19-Aug 18:39:47:DEBUG:root:Displacement of element = \n", + "[0.0 -0.39914506095474334 -0.0722858269559246 0.0\n", + " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n" ] }, { @@ -404,7 +397,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 10, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -412,41 +405,26 @@ "source": [ "facts(\"one element assembly\") do\n", " # Create model\n", - " #Logging.debug(\"Creating nodes\")\n", - " #n1 = Node(1)\n", - " #n2 = Node(2)\n", - " #n3 = Node(3)\n", - " #n4 = Node(4)\n", " Logging.debug(\"Adding nodes to array\")\n", - " #nodes = [n1.id, n2.id, n3.id, n4.id]\n", - " node_ids = [1, 2, 3, 4]\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, node_ids, 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", + " el = get_test_element()\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 = JuliaFEM.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", "\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", + " F = [0.0 0.0; 0.0 0.0; 0.0 -2.0; 0.0 0.0]'\n", + " F = F[:]\n", + " free_dofs = [3, 4, 5, 6]\n", "\n", " # solution\n", " K = sparse(ass.I, ass.J, ass.A)\n", " R = full(sparsevec(ass.i, ass.b))\n", " R = R - F\n", - " du = zeros(8) # must be determined from ass\n", + " du = zeros(8)\n", " du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n", "\n", " Logging.debug(\"Solution norm = $(norm(du))\")\n", @@ -462,7 +440,7 @@ " end\n", " disp = el.attributes[\"displacement\"]\n", " Logging.debug(\"Displacement of element = \\n$disp\")\n", - " @fact norm(disp) => roughly(3.1292483947150043)\n", + " @fact norm(disp) --> roughly(3.1292483947150043)\n", "end" ] }, @@ -479,7 +457,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -493,83 +471,18 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 8, "metadata": { - "collapsed": false, - "scrolled": false + "collapsed": false }, "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "solve two element problem\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "15-Aug 17:03:59:DEBUG:root:Creating elements\n", - "15-Aug 17:03:59:INFO:root:create_ldof2gdofmap: dofs per node: 2\n", - "15-Aug 17:04:00:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],5=>[9,10],6=>[11,12],1=>[1,2])\n", - "15-Aug 17:04:04:INFO:root:solve!: dofs per node: 2\n", - "15-Aug 17:04:04:DEBUG:root:Problem size = 12\n", - "15-Aug 17:04:04:DEBUG:root:Starting iteration 1\n", - "15-Aug 17:04:04:DEBUG:root:Assembling\n", - "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n", - "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n", - "15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n", - "15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.6015838690633517\n", - "15-Aug 17:04:05:DEBUG:root:Starting iteration 2\n", - "15-Aug 17:04:05:DEBUG:root:Assembling\n", - "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n", - "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n", - "15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n", - "15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.013420417380980414\n", - "15-Aug 17:04:05:DEBUG:root:Starting iteration 3\n", - "15-Aug 17:04:05:DEBUG:root:Assembling\n", - "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n", - "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n", - "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n", - "15-Aug 17:04:06:DEBUG:root:Solution norm du = 0.00032202957936854873\n", - "15-Aug 17:04:06:DEBUG:root:Starting iteration 4\n", - "15-Aug 17:04:06:DEBUG:root:Assembling\n", - "15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n", - "15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n", - "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n", - "15-Aug 17:04:06:DEBUG:root:Solution norm du = 9.906677094801476e-8\n", - "15-Aug 17:04:06:DEBUG:root:Starting iteration 5\n", - "15-Aug 17:04:06:DEBUG:root:Assembling\n", - "15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n", - "15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n", - "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n", - "15-Aug 17:04:06:DEBUG:root:Solution norm du = 5.027624635820698e-15\n", - "15-Aug 17:04:06:DEBUG:root:Converged in 5 iterations.\n", - "15-Aug 17:04:06:DEBUG:root:Displacement of element = \n", - "[-0.02442313597467864 -0.039356000063335075 0.021097993207232584 0.020877031423993653\n", - " -0.12673626841485705 -0.40433021969759375 -0.40656320177872923 -0.1275776940913048]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "0 facts verified.\n" - ] - }, { "data": { "text/plain": [ - "delayed_handler (generic function with 4 methods)" + "create_ldof2gdofmap (generic function with 1 method)" ] }, - "execution_count": 12, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } @@ -578,18 +491,17 @@ "\"\"\"\n", "Create local dof to global dof mapping for given elements\n", "\"\"\"\n", - "function create_ldof2gdofmap(elements; ndofs=2)\n", - " Logging.info(\"create_ldof2gdofmap: dofs per node: $ndofs\")\n", + "function create_ldof2gdofmap(elements, field)\n", "\n", + " ndofs = size(elements[1].attributes[field], 1)\n", + " \n", " all_node_ids = Int64[]\n", " for el in elements\n", - " eldim, elnodes = size(el.coordinates)\n", " for nid in el.node_ids\n", " push!(all_node_ids, nid)\n", " end\n", " end\n", " all_node_ids = unique(all_node_ids)\n", - " sort!(all_node_ids)\n", "\n", " # Assign global dof for each node\n", " pdim = 1\n", @@ -600,24 +512,104 @@ " end\n", "\n", " return ngdofs\n", - "end\n", - "\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false, + "scrolled": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "solve one element problem\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "19-Aug 18:49:12:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n", + "19-Aug 18:49:13:INFO:root:solve!: dofs per node: 2\n", + "19-Aug 18:49:13:DEBUG:root:Problem size = 8\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 1\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 3.0900221367289444\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 2\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.32121316021534796\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 3\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.040431781940014504\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 4\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.0009291101052065917\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 5\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n", + "19-Aug 18:49:13:DEBUG:root:Starting iteration 6\n", + "19-Aug 18:49:13:DEBUG:root:Assembling\n", + "19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n", + "19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n", + "19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n", + "19-Aug 18:49:13:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n", + "19-Aug 18:49:13:DEBUG:root:Converged in 6 iterations.\n", + "19-Aug 18:49:13:DEBUG:root:Displacement of element = \n", + "[0.0 -0.39914506095474345 -0.07228582695592467 0.0\n", + " 0.0 -2.177989231707351 -2.222244754401765 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": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ "function solve!(elements, dofmap, neumann_bcs, dirichlet_bcs; ndofs=2, max_iterations=10)\n", "\n", - " dbc = \"lagrange\"\n", - " \n", " Logging.info(\"solve!: dofs per node: $ndofs\")\n", " pdim = length(dofmap)*ndofs\n", " Logging.debug(\"Problem size = $pdim\")\n", "\n", - " # Initialize elements ready for solution\n", - " for el in elements\n", - " eldim, elnodes = size(el.coordinates)\n", - " el.attributes[\"displacement\"] = zeros(ndofs, elnodes)\n", - " el.attributes[\"displacement nodal force\"] = zeros(ndofs, elnodes)\n", - " el.attributes[\"displacement tangent stiffness\"] = zeros(ndofs*elnodes, ndofs*elnodes)\n", - " end\n", - "\n", " # Assign global dofs for elements\n", " gdofs = Dict{Int64, Array{Int64,1}}()\n", " for el in elements\n", @@ -629,7 +621,7 @@ " end\n", " end\n", "\n", - " ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n", + " ass = JuliaFEM.Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n", "\n", " for iter=1:max_iterations\n", " Logging.debug(\"Starting iteration $iter\")\n", @@ -642,30 +634,26 @@ " Logging.debug(\"Assembling\")\n", " for el in elements\n", " assemble_element!(ass, el)\n", - " #println(\"Element stiffness matrix\")\n", - " #dump(round(el.attributes[\"displacement tangent stiffness\"], 2))\n", " end\n", "\n", " i = 1\n", - " if dbc == \"lagrange\"\n", - " Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n", - " # Dirichlet boundary conditions\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", + " Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n", + " # Dirichlet boundary conditions\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", - " Logging.debug(\"Added $i Lagrange multipliers to model\")\n", " end\n", + " Logging.debug(\"Added $i Lagrange multipliers to model\")\n", " i -= 1\n", "\n", " Logging.debug(\"Adding Neumann boundary conditions\")\n", @@ -682,20 +670,8 @@ " K = sparse(ass.I, ass.J, ass.A)\n", " R = full(sparsevec(ass.i, ass.b))\n", " R = R - F\n", - " #Logging.debug(dump(round(full(K), 2)))\n", - " #print_matrix(full(K))\n", - " #Logging.debug(dump(round(R', 1)))\n", "\n", - " if dbc == \"eliminate\"\n", - " Logging.debug(\"Eliminating Dirichlet boundary conditions\")\n", - " throw(\"Implement this properly\")\n", - " free_dofs = [1, 2, 3, 4, 5, 6, 7, 8]\n", - " du = zeros(12)\n", - " Logging.debug(\"K norm = $(norm(full(K[free_dofs, free_dofs])))\")\n", - " du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n", - " else\n", - " du = K \\ -R\n", - " end\n", + " du = K \\ -R\n", "\n", " #du = reshape(du, 2, 6)\n", " solnorm = norm(du[1:pdim])\n", @@ -725,55 +701,27 @@ "\n", "ENV[\"COLUMNS\"] = 160\n", "\n", - "function test1():\n", - " facts(\"solve one element problem\") do\n", - " # Create model\n", - " Logging.debug(\"Creating nodes\")\n", - " node_ids = [1, 2, 3, 4]\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, node_ids, coordinates, attributes)\n", - " elements = [el]\n", - " dofmap = create_ldof2gdofmap(elements)\n", - " Logging.debug(dofmap)\n", - " # Boundary conditions\n", - " # here we want to create nodal force for third dof, that is, node id 2, second dof\n", - " bc1 = BC([dofmap[2][2]], [-2.0])\n", - " # dirichlet bc, set dx=dy=0 on support\n", - " bc2 = BC([dofmap[3][1], dofmap[3][2], dofmap[4][1], dofmap[4][2]], [0.0, 0.0, 0.0, 0.0])\n", - " solve!(elements, dofmap, [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\n", - "end\n", - "\n", - "facts(\"solve two element problem\") do\n", + "facts(\"solve one element problem\") do\n", " # Create model\n", - " attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n", - "\n", - " Logging.debug(\"Creating elements\")\n", - " nids1 = [5, 1, 3, 6]\n", - " coords1 = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'\n", - " el1 = Element(1, nids1, coords1, copy(attributes))\n", - " nids2 = [1, 2, 4, 3]\n", - " coords2 = [5.0 0.0; 10.0 0.0; 10.0 1.0; 5.0 1.0]'\n", - " el2 = Element(2, nids2, coords2, copy(attributes))\n", - " elements = [el1, el2]\n", - "\n", - " dofmap = create_ldof2gdofmap(elements)\n", + " el = get_test_element()\n", + " elements = [el]\n", + " # Initialize elements ready for solution\n", + " for el in elements\n", + " eldim, elnodes = size(el.attributes[\"coordinates\"])\n", + " el.attributes[\"displacement\"] = zeros(2, elnodes)\n", + " end\n", + " dofmap = create_ldof2gdofmap(elements, \"displacement\")\n", " Logging.debug(dofmap)\n", " # Boundary conditions\n", " # here we want to create nodal force for third dof, that is, node id 2, second dof\n", - " bc1 = BC([dofmap[4][2]], [-0.1])\n", + " bc1 = BC([dofmap[3][2]], [-2.0])\n", " # dirichlet bc, set dx=dy=0 on support\n", - " bc2 = BC([dofmap[5][1], dofmap[5][2], dofmap[6][1], dofmap[6][2]], [0.0, 0.0, 0.0, 0.0])\n", + " bc2 = BC([dofmap[1][1], dofmap[1][2], dofmap[4][1], dofmap[4][2]], [0.0, 0.0, 0.0, 0.0])\n", " solve!(elements, dofmap, [bc1], [bc2]; max_iterations=7)\n", - " disp = elements[2].attributes[\"displacement\"]\n", + " disp = elements[1].attributes[\"displacement\"]\n", " Logging.debug(\"Displacement of element = \\n$disp\")\n", - " #@fact norm(disp) => roughly(3.1292483947150043)\n", - "end" + " @fact norm(disp) --> roughly(3.1292483947150043)\n", + "end\n" ] }, { @@ -787,7 +735,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -796,68 +744,44 @@ "name": "stderr", "output_type": "stream", "text": [ - "\n", - "WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:11.\n", - "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", - "15-Aug 17:04:20:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", - "WARNING: beginswith is deprecated, use startswith instead.\n", - " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", - " in beginswith at deprecated.jl:30\n", - " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n", - " in include_string at loading.jl:99\n", - " in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n", - " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n", - " in anonymous at task.jl:365\n", - "while loading In[13], in expression starting on line 3\n", - "WARNING: beginswith is deprecated, use startswith instead.\n", - " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", - " in beginswith at deprecated.jl:30\n", - " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n", - " in include_string at loading.jl:99\n", - " in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n", - " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n", - " in anonymous at task.jl:365\n", - "while loading In[13], in expression starting on line 3\n", - "15-Aug 17:04:21:DEBUG:root:Found NODE section\n", - "15-Aug 17:04:22:DEBUG:root:Found ELEMENT section\n", - "WARNING: integer(s::AbstractString) is deprecated, use parse(Int,s) instead.\n", - " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", - " in integer at deprecated.jl:49\n", - " in map at abstractarray.jl:1251\n", - " in parse_element_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:59\n", - " in process_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:108\n", - " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:117\n", - " in include_string at loading.jl:99\n", - " in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n", - " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n", - " in anonymous at task.jl:365\n", - "while loading In[13], in expression starting on line 3\n", - "15-Aug 17:04:27:DEBUG:root:120 elements found\n", - "15-Aug 17:04:27:INFO:root:Creating ELSET Body1\n", - "15-Aug 17:04:27:DEBUG:root:Found NSET section\n", - "15-Aug 17:04:28:DEBUG:root:Creating node set SUPPORT\n", - "15-Aug 17:04:28:DEBUG:root:Found NSET section\n", - "15-Aug 17:04:28:DEBUG:root:Creating node set LOAD\n", - "15-Aug 17:04:28:DEBUG:root:Found NSET section\n", - "15-Aug 17:04:28:DEBUG:root:Creating node set TOP\n" + "19-Aug 18:49:21:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", + "WARNING: beginswith is deprecated, use startswith instead.\n" ] }, { - "data": { - "text/plain": [ - "Dict{Any,Any} with 4 entries:\n", - " \"nodes\" => Dict{Any,Any}(288=>[97.5,7.5,10.0],11=>[92.5,2.5,5.0],134=>[45.0,10.0,0.0],158=>[2.5,2.5,0.0],160=>[7.5,7.5,0.0],215=>[60.0,0.0,5.0],29=>[2.5,7…\n", - " \"elements\" => Dict{Any,Any}(68=>[71,144,149,198,51,150,57,43,50,214],2=>[204,199,175,130,207,208,209,3,4,176],89=>[95,78,104,52,127,126,106,60,68,67],11=>[15…\n", - " \"elsets\" => Dict{Any,Any}(\"Body1\"=>[1,2,3,4,5,6,7,8,9,10 … 111,112,113,114,115,116,117,118,119,120])\n", - " \"nsets\" => Dict{Any,Any}(\"LOAD\"=>[82,84,87,179,197,246,249,256,257],\"SUPPORT\"=>[108,109,111,155,162,216,225,281,298],\"TOP\"=>[70,75,76,84,88,90,95,96,98,10…" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" + "ename": "LoadError", + "evalue": "LoadError: MethodError: `convert` has no method matching convert(::Type{SubString{ASCIIString}}, ::Dict{Any,Any})\nThis may have arisen from a call to the constructor SubString{ASCIIString}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T<:AbstractString}(::Type{T<:AbstractString}, !Matched::AbstractArray{Char,1})\n convert{T<:AbstractString}(::Type{SubString{T<:AbstractString}}, !Matched::T<:AbstractString)\n ...\nwhile loading In[17], in expression starting on line 3", + "output_type": "error", + "traceback": [ + "LoadError: MethodError: `convert` has no method matching convert(::Type{SubString{ASCIIString}}, ::Dict{Any,Any})\nThis may have arisen from a call to the constructor SubString{ASCIIString}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T<:AbstractString}(::Type{T<:AbstractString}, !Matched::AbstractArray{Char,1})\n convert{T<:AbstractString}(::Type{SubString{T<:AbstractString}}, !Matched::T<:AbstractString)\n ...\nwhile loading In[17], in expression starting on line 3", + "", + " in setindex! at dict.jl:615", + " in parse_header at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:40", + " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:120" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", + " in beginswith at deprecated.jl:30\n", + " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:115\n", + " in include_string at loading.jl:99\n", + " in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n", + " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n", + " in anonymous at task.jl:365\n", + "while loading In[17], in expression starting on line 3\n", + "WARNING: beginswith is deprecated, use startswith instead.\n", + " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n", + " in beginswith at deprecated.jl:30\n", + " in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:115\n", + " in include_string at loading.jl:99\n", + " in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n", + " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n", + " in anonymous at task.jl:365\n", + "while loading In[17], in expression starting on line 3\n" + ] } ], "source": [ @@ -1280,6 +1204,53 @@ "JuliaFEM.xdmf.xdmf_save_model(xdoc, \"/tmp/3d_solid_model.xmf\")" ] }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0,0.0])" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u = zeros(10)\n", + "u1 = u[1:5]\n", + "u, u1" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[1.0,1.0,1.0,1.0,1.0])" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "u1[:] = 1\n", + "u, u1" + ] + }, { "cell_type": "code", "execution_count": null, diff --git a/src/abaqus_reader.jl b/src/abaqus_reader.jl index 4791762..3e52099 100644 --- a/src/abaqus_reader.jl +++ b/src/abaqus_reader.jl @@ -8,7 +8,9 @@ using Logging VERSION < v"0.4-" && using Docile -eldims = Dict({"C3D10" => 10}) +eldims = Dict( + "C3D10" => 10, + "C3D4" => 4) global handlers = Dict() @@ -29,7 +31,7 @@ end function parse_header(header_line) args = map(s -> strip(s), split(header_line, ",")) args[1] = strip(args[1], '*') - d = Dict({"section" => args[1]}) + d = Dict("section" => args[1]) options = Dict() for k in args[2:end] args2 = split(k, "=") diff --git a/src/math.jl b/src/math.jl index b839cb6..d3dc9bf 100644 --- a/src/math.jl +++ b/src/math.jl @@ -73,9 +73,9 @@ end """ Return partial derivatives of shape functions w.r.t X using chain rule. """ -function get_dbasisdX(el::Element, xi) - J = interpolate(el, "coordinates", xi; derivative=true) - dbasisdX = el.dbasis(xi)*inv(J') +function get_dbasisdX(el::Element, ip) + J = interpolate(el, "coordinates", ip.xi; derivative=true) + dbasisdX = el.dbasis(ip.xi)*inv(J') return dbasisdX end @@ -156,11 +156,9 @@ f::Function """ function integrate(f::Function, el::JuliaFEM.Element) target = [] - for m = 1:length(el.iweights) - w = el.iweights[m] - xi = el.ipoints[:, m] - J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true) - push!(target, w*f(el, xi)*det(J)) + for ip in el.integration_points + J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true) + push!(target, ip.weight*f(el, ip)*det(J)) end return sum(target) end @@ -171,11 +169,9 @@ This version returns a function which must be operated with element e function integrate(f::Function) function integrate(el::JuliaFEM.Element) target = [] - for m = 1:length(el.iweights) - w = el.iweights[m] - xi = el.ipoints[:, m] - J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true) - push!(target, w*f(el, xi)*det(J)) + for ip in el.integration_points + J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true) + push!(target, ip.weight*f(el, ip)*det(J)) end return sum(target) end @@ -188,10 +184,8 @@ This version saves results inplace to target, garbage collection free function integrate!(f::Function, el::JuliaFEM.Element, target) # set target to zero el.attributes[target][:] = 0.0 - for m = 1:length(el.iweights) - w = el.iweights[m] - xi = el.ipoints[:, m] - J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true) - el.attributes[target][:,:] += w*f(el, xi)*det(J) + for ip in el.integration_points + J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true) + el.attributes[target][:,:] += ip.weight*f(el, ip)*det(J) end end diff --git a/src/types.jl b/src/types.jl index 05d38e9..68e663b 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,6 +1,22 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +""" +Integration point + +xi :: Array{Float64, 1} + (dimensionless) coordinates of integration point +weight :: Float64 + Integration weight +attributes :: Dict{ASCIIString, Any} + This is used to save internal variables of IP needed e.g. for incremental + material models. +""" +type IntegrationPoint + xi :: Array{Float64, 1} + weight :: Float64 + attributes :: Dict{ASCIIString, Any} +end type Element id :: Int @@ -8,9 +24,10 @@ type Element node_ids :: Array{Int, 1} basis :: Function dbasis :: Function + integration_points :: Array{IntegrationPoint, 1} attributes :: Dict{ASCIIString, Any} - ipoints :: Array{Float64, 2} - iweights :: Array{Float64, 1} +# ipoints :: Array{Float64, 2} +# iweights :: Array{Float64, 1} end