From 66963fffe75d74f46c45c4ded4d5caf1d2fb3490 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 24 Aug 2015 01:14:03 +0300 Subject: [PATCH] Broken again. --- ...2015-06-25-elasticity-solver-example.ipynb | 797 ++++++++++-------- src/JuliaFEM.jl | 5 + src/elements.jl | 92 ++ src/math.jl | 19 +- 4 files changed, 546 insertions(+), 367 deletions(-) create mode 100644 src/elements.jl diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 70230a8..59e621c 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -41,16 +41,6 @@ "collapsed": false }, "outputs": [ - { - "data": { - "text/plain": [ - "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - }, { "name": "stderr", "output_type": "stream", @@ -69,12 +59,25 @@ " 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" + "while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n", + "24-Aug 00:56:48:INFO:root:loading types\n", + "24-Aug 00:56:48:INFO:root:loading elements\n" ] + }, + { + "data": { + "text/plain": [ + "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ "using JuliaFEM\n", + "using JuliaFEM: Assembly, Element, Quad4, get_integration_points, get_jacobian, get_basis, get_dbasisdxi, get_dbasisdX\n", "using Logging\n", "Logging.configure(level=DEBUG)" ] @@ -87,9 +90,14 @@ "\n", "*Design principle 4*: we don't use greek characters in code which is implemented to JuliaFEM. In notebooks they are ok.\n", "\n", - "*Design principle 5*: we use 4 space indentation like in Python.\n", - "\n", - "First we construct some type for our element which contains all relevant data. We don't care a much how every element is actually implemented as long as it follows some general rules how the interface is constructed. Our element implementation for continuum 3d element is" + "*Design principle 5*: we use 4 space indentation like in Python." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First we construct some type for our element which contains all relevant data. We don't care a much how every element is actually implemented as long as it follows some general rules how the interface is constructed. First we define our element family and it's basis functions, derivatives of them etc. These needs to be defined for each element type only once." ] }, { @@ -98,28 +106,6 @@ "metadata": { "collapsed": false }, - "outputs": [], - "source": [ - "abstract ContinuumElement <: Element\n", - "\n", - "\"\"\"\n", - "4-node bilinear plane stress element\n", - "\"\"\"\n", - "type CPS4 <: ContinuumElement\n", - " id :: Int\n", - " node_ids :: Array{Int, 1}\n", - " shape_functions :: FunctionSpace\n", - " integration_points :: Array{IntegrationPoint, 1}\n", - " attributes :: Dict{ASCIIString, Any}\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "metadata": { - "collapsed": false - }, "outputs": [ { "data": { @@ -127,15 +113,41 @@ "get_rhs (generic function with 1 method)" ] }, - "execution_count": 3, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ + "\"\"\"\n", + "Integrate f over element using Gaussian quadrature rules.\n", + "\n", + "Parameters\n", + "----------\n", + "el::Element\n", + " well defined element\n", + "f::Function\n", + " Function to integrate\n", + "\"\"\"\n", + "function integrate(el::Element, f::Function)\n", + " target = []\n", + " for ip in get_integration_points(el)\n", + " J = get_jacobian(el, ip.xi)\n", + " push!(target, ip.weight*f(el, ip)*det(J))\n", + " end\n", + " return sum(target)\n", + "end\n", + "\n", + "\"\"\"\n", + "Return left hand side of the equation Ax = b (i.e. A)\n", + "\"\"\"\n", "function get_lhs(el::Element)\n", " return None\n", "end\n", + "\n", + "\"\"\"\n", + "Return right hand side of the equation Ax = b (i.e. b)\n", + "\"\"\"\n", "function get_rhs(el::Element)\n", " return None\n", "end" @@ -145,7 +157,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next we need to think a bit of interface design. We give good defaults if element is constructed like this, so user doesn't have to provide everything (although it's totally possible). Here's the implementation how to calculate internal nodal forces for some integration point in continuum elements general:" + "Next we define our elements for mechanical problem. Because calculating internal and external energy for several types of elements follow same procedure, we construct whole family of mechanical elements which share common functions." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "#abstract Mechanical <: CG" ] }, { @@ -158,7 +181,7 @@ { "data": { "text/plain": [ - "get_lhs (generic function with 2 methods)" + "get_field (generic function with 1 method)" ] }, "execution_count": 4, @@ -166,14 +189,95 @@ "output_type": "execute_result" } ], + "source": [ + "\"\"\"\n", + "(4-node bilinear) plane stress element\n", + "\"\"\"\n", + "type CPS4 <: Quad4\n", + " id :: Int\n", + " node_ids :: Array{Int, 1}\n", + " coordinates :: Array{Float64, 2}\n", + " integration_points :: Array{IntegrationPoint, 1}\n", + " attributes :: Dict{ASCIIString, Any}\n", + "end\n", + "\n", + "function get_field(el::Quad4)\n", + " el.attributes[\"displacement\"]\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Each element must provide a standard way how it's initialized. For constructor we need only unique element id and it's connectivity data." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "CPS4" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function CPS4(element_id, node_ids)\n", + " integration_points = [\n", + " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0, Dict()),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0, Dict())]\n", + " attributes = Dict(\"displacement\" => zeros(2, 4))\n", + " coordinates = zeros(2, 4)\n", + " CPS4(element_id, node_ids, coordinates, integration_points, attributes)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's the actual basic implementation for mechanical elements\n", + "\n", + "Our task is: for given $\\mathbf{u}$ calculate $\\mathbf{R}(\\mathbf{u}) = \\mathbf{T}(\\mathbf{u}) - \\mathbf{F}(\\mathbf{u})$ and it's partial derivative with respect to $\\mathbf{u}$, i.e. $\\partial \\mathbf{R}(\\mathbf{u}) / \\partial \\mathbf{u}$. Here is our $\\mathbf{T}$:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "get_lhs (generic function with 2 methods)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "\"\"\"\n", "Calculate internal nodal forces for continuum element.\n", "\"\"\"\n", - "function Wint(el::ContinuumElement)\n", + "function Wint(el::Quad4)\n", "\n", - " dNdX(xi) = JuliaFEM.get_dbasisdX(el, xi)\n", - " #dNdX(xi) = el.shape_functions.dbasis(xi)\n", + " dNdX(xi) = get_dbasisdX(el, xi)\n", " # material\n", " lambda(xi) = interpolate(el, \"lambda\", xi)\n", " mu(xi) = interpolate(el, \"mu\", xi)\n", @@ -185,25 +289,11 @@ " S(xi, u) = lambda(xi)*trace(E(xi, u))*I + 2*mu(xi)*E(xi, u)\n", " P(xi, u) = F(xi, u)*S(xi, u)\n", " T(xi, u) = P(xi, u)*dNdX(xi)'\n", - "\n", - " function Wint_(el, ip)\n", - " xi = ip.xi\n", - " u = el.attributes[\"displacement\"]\n", - " return T(xi, u)\n", - " end\n", - " return integrate(Wint_, el)\n", - " #return integrate(T, el.integration_points, el.attributes[\"displacement\"])\n", + " integrate(el, (el, ip) -> T(ip.xi, get_field(el)))\n", "end\n", "\n", - "get_rhs(el::ContinuumElement) = -Wint(el)\n", - "get_lhs(el::ContinuumElement) = linearize(Wint, \"displacement\")(el)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Our task is: for given $\\mathbf{u}$ calculate $\\mathbf{R}(\\mathbf{u}) = \\mathbf{T}(\\mathbf{u}) - \\mathbf{F}(\\mathbf{u})$ and it's partial derivative with respect to $\\mathbf{u}$, i.e. $\\partial \\mathbf{R}(\\mathbf{u}) / \\partial \\mathbf{u}$. Here is our $\\mathbf{T}$:" + "get_rhs(el::Quad4) = -Wint(el) # rhs = -R = -(T-F)\n", + "get_lhs(el::Quad4) = linearize(Wint, \"displacement\")(el)" ] }, { @@ -217,7 +307,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 7, "metadata": { "collapsed": true }, @@ -228,7 +318,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -239,7 +329,7 @@ "get_test_element (generic function with 1 method)" ] }, - "execution_count": 6, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } @@ -247,47 +337,25 @@ "source": [ "function get_test_element()\n", " # set up one linear quadrangle element\n", - " basis(xi) = [\n", - " (1-xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1+xi[2])/4\n", - " (1-xi[1])*(1+xi[2])/4]\n", - "\n", - " dbasisdxi(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n", - " (1-xi[2])/4.0 -(1+xi[1])/4.0\n", - " (1+xi[2])/4.0 (1+xi[1])/4.0\n", - " -(1+xi[2])/4.0 (1-xi[1])/4.0]\n", - "\n", - " X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n", - " #J(xi) = X*dbasisdxi(xi)\n", - " #dbasisdX(xi) = dbasisdxi(xi)*inv(J(xi)')\n", - " shape_functions = FunctionSpace(basis, dbasisdxi)\n", - " integration_points = [\n", - " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0, Dict())]\n", - " attributes = Dict()\n", " element_id = 1\n", " node_ids = [1, 2, 3, 4]\n", - " el = CPS4(element_id, node_ids, shape_functions, integration_points, attributes)\n", - " #el.shape_functions.dbasis(xi) = JuliaFEM.get_dbasisdX(el, xi)\n", + " el = CPS4(element_id, node_ids)\n", + "\n", " E = 90.0\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", - " el.attributes[\"coordinates\"] = X\n", - " el.attributes[\"lambda\"] = la\n", - " el.attributes[\"mu\"] = mu\n", - " el.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n", + " X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n", + " set_coordinates(el, X)\n", + " set_material(el, la, mu)\n", " return el\n", "end" ] }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 9, "metadata": { "collapsed": false, "scrolled": false @@ -304,40 +372,40 @@ "name": "stderr", "output_type": "stream", "text": [ - "22-Aug 20:45:52:DEBUG:root:Iteration 1\n", - "22-Aug 20:46:00:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:02:DEBUG:root:Iteration 2\n", - "22-Aug 20:46:02:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:02:DEBUG:root:Iteration 3\n", - "22-Aug 20:46:02:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:03:DEBUG:root:Iteration 4\n", - "22-Aug 20:46:03:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:03:DEBUG:root:Iteration 5\n", - "22-Aug 20:46:03:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:03:DEBUG:root:Iteration 6\n", - "22-Aug 20:46:03:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:03:DEBUG:root:Converged in 6 iterations.\n", - "22-Aug 20:46:03:DEBUG:root:solution vector: \n", + "24-Aug 00:57:07:DEBUG:root:Iteration 1\n", + "24-Aug 00:57:10:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Iteration 2\n", + "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Iteration 3\n", + "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Iteration 4\n", + "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Iteration 5\n", + "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Iteration 6\n", + "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:13:DEBUG:root:Converged in 6 iterations.\n", + "24-Aug 00:57:14:DEBUG:root:solution vector: \n", " [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n", - "22-Aug 20:46:04:DEBUG:root:norm of u: 3.1292483947150047\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 1\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 2\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 3\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 4\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 5\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Iteration 6\n", - "22-Aug 20:46:04:DEBUG:root:Solving Ax = b\n", - "22-Aug 20:46:04:DEBUG:root:Converged in 6 iterations.\n", - "22-Aug 20:46:04:DEBUG:root:solution vector: \n", - " [0.0 0.7433248532717796 1.048521014723486 0.0\n", - " 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n", - "22-Aug 20:46:04:DEBUG:root:norm of u: 3.1292483947150056\n" + "24-Aug 00:57:14:DEBUG:root:norm of u: 3.1292483947150047\n", + "24-Aug 00:57:14:DEBUG:root:Iteration 1\n", + "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:14:DEBUG:root:Iteration 2\n", + "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:14:DEBUG:root:Iteration 3\n", + "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:14:DEBUG:root:Iteration 4\n", + "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:14:DEBUG:root:Iteration 5\n", + "24-Aug 00:57:15:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:15:DEBUG:root:Iteration 6\n", + "24-Aug 00:57:15:DEBUG:root:Solving Ax = b\n", + "24-Aug 00:57:15:DEBUG:root:Converged in 6 iterations.\n", + "24-Aug 00:57:15:DEBUG:root:solution vector: \n", + " [0.0 1.2578327758133292 1.5202505368695098 0.0\n", + " 0.0 -1.8223091343697626 -1.6224781337179326 0.0]\n", + "24-Aug 00:57:15:DEBUG:root:norm of u: 3.129248394715004\n" ] }, { @@ -353,7 +421,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 7, + "execution_count": 9, "metadata": {}, "output_type": "execute_result" } @@ -369,8 +437,9 @@ " free_dofs = [3, 4, 5, 6]\n", " for i=1:10\n", " Logging.debug(\"Iteration $i\")\n", - " A = get_lhs(e)\n", " b = get_rhs(e)\n", + " #Logging.debug(\"rhs = $b\")\n", + " A = get_lhs(e)\n", " Logging.debug(\"Solving Ax = b\")\n", " du[free_dofs] = A[free_dofs, free_dofs] \\ (b + F)[free_dofs]\n", "\n", @@ -389,11 +458,11 @@ " Logging.debug(\"norm of u: $(norm(u))\")\n", "\n", " # We rotate model a bit and make sure that norm remains same\n", - " phi = 30/180*pi\n", + " phi = 45/180*pi\n", " rmat = [\n", " cos(phi) -sin(phi)\n", " sin(phi) cos(phi)]\n", - " e.attributes[\"coordinates\"] = rmat*e.attributes[\"coordinates\"]\n", + " set_coordinates(e, rmat*get_coordinates(e))\n", " F = rmat*F\n", "\n", " e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n", @@ -427,7 +496,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 10, "metadata": { "collapsed": false }, @@ -435,16 +504,16 @@ { "data": { "text/plain": [ - "assemble_element! (generic function with 1 method)" + "assemble_rhs! (generic function with 1 method)" ] }, - "execution_count": 8, + "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function assemble_element!(ass::Assembly, el::Element)\n", + "function assemble_lhs!(ass::Assembly, el::Element)\n", "\n", " gdofs = ass.gdofs[el.id]\n", "\n", @@ -459,6 +528,11 @@ " end\n", " end\n", " end\n", + "end\n", + "\n", + "function assemble_rhs!(ass::Assembly, el::Element)\n", + "\n", + " gdofs = ass.gdofs[el.id]\n", " \n", " b = get_rhs(el)\n", " if !(b == None)\n", @@ -472,7 +546,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -480,19 +554,19 @@ { "data": { "text/plain": [ - "get_field (generic function with 1 method)" + "get_field (generic function with 2 methods)" ] }, - "execution_count": 9, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function update_field(el::ContinuumElement, du)\n", + "function update_field(el::Element, du)\n", " el.attributes[\"displacement\"][:] += du\n", "end\n", - "function get_field(el::ContinuumElement)\n", + "function get_field(el::Element)\n", " return el.attributes[\"displacement\"]\n", "end" ] @@ -504,46 +578,6 @@ "We also need to construct our element in somehow \"standard\" way. My proposal is: element id and node ids (connectivity) information. Of course other fields must also be provided. Here's example for CPS4 element:" ] }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "CPS4" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function CPS4(element_id, node_ids)\n", - " basis(xi) = [\n", - " (1-xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1-xi[2])/4\n", - " (1+xi[1])*(1+xi[2])/4\n", - " (1-xi[1])*(1+xi[2])/4]\n", - " dbasis(xi) = [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n", - " (1-xi[2])/4.0 -(1+xi[1])/4.0\n", - " (1+xi[2])/4.0 (1+xi[1])/4.0\n", - " -(1+xi[2])/4.0 (1-xi[1])/4.0]\n", - " shape_functions = FunctionSpace(basis, dbasis)\n", - " integration_points = [\n", - " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0, Dict()),\n", - " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0, Dict())]\n", - " attributes = Dict(\"displacement\" => zeros(2, 4))\n", - " CPS4(element_id, node_ids, shape_functions, integration_points, attributes)\n", - "end" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -553,7 +587,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -564,7 +598,7 @@ "set_attribute (generic function with 1 method)" ] }, - "execution_count": 11, + "execution_count": 12, "metadata": {}, "output_type": "execute_result" } @@ -584,16 +618,18 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 13, "metadata": { - "collapsed": true + "collapsed": false }, "outputs": [], "source": [ + "using JuliaFEM: CG\n", + "abstract Point0 <: CG\n", "\"\"\"\n", "1-node point force element for plane stress problems.\n", "\"\"\"\n", - "type CPS1 <: ContinuumElement\n", + "type CPS1 <: Point0\n", " id :: Int\n", " node_ids :: Array{Int, 1}\n", " attributes :: Dict{ASCIIString, Any}\n", @@ -602,7 +638,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -613,7 +649,7 @@ "get_lhs (generic function with 3 methods)" ] }, - "execution_count": 13, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -655,7 +691,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -671,26 +707,38 @@ "name": "stderr", "output_type": "stream", "text": [ - "22-Aug 20:46:07:DEBUG:root:Starting iteration 1\n", - "22-Aug 20:46:07:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 3.0900221367289986\n", - "22-Aug 20:46:08:DEBUG:root:Starting iteration 2\n", - "22-Aug 20:46:08:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 0.3212131602153472\n", - "22-Aug 20:46:08:DEBUG:root:Starting iteration 3\n", - "22-Aug 20:46:08:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 0.04043178193999703\n", - "22-Aug 20:46:08:DEBUG:root:Starting iteration 4\n", - "22-Aug 20:46:08:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 0.0009291101052105739\n", - "22-Aug 20:46:08:DEBUG:root:Starting iteration 5\n", - "22-Aug 20:46:08:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 1.5638899027804743e-7\n", - "22-Aug 20:46:08:DEBUG:root:Starting iteration 6\n", - "22-Aug 20:46:08:DEBUG:root:Assembling\n", - "22-Aug 20:46:08:DEBUG:root:Solution norm = 1.0464940956129567e-14\n", - "22-Aug 20:46:08:DEBUG:root:Converged in 6 iterations.\n", - "22-Aug 20:46:08:DEBUG:root:Displacement of element = \n", + "24-Aug 01:06:39:DEBUG:root:Starting iteration 1\n", + "24-Aug 01:06:39:DEBUG:root:Assembling\n", + "24-Aug 01:06:39:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:39:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 3.0900221367289986\n", + "24-Aug 01:06:40:DEBUG:root:Starting iteration 2\n", + "24-Aug 01:06:40:DEBUG:root:Assembling\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.3212131602153472\n", + "24-Aug 01:06:40:DEBUG:root:Starting iteration 3\n", + "24-Aug 01:06:40:DEBUG:root:Assembling\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.04043178193999703\n", + "24-Aug 01:06:40:DEBUG:root:Starting iteration 4\n", + "24-Aug 01:06:40:DEBUG:root:Assembling\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.0009291101052105739\n", + "24-Aug 01:06:40:DEBUG:root:Starting iteration 5\n", + "24-Aug 01:06:40:DEBUG:root:Assembling\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 1.5638899027804743e-7\n", + "24-Aug 01:06:40:DEBUG:root:Starting iteration 6\n", + "24-Aug 01:06:40:DEBUG:root:Assembling\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n", + "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n", + "24-Aug 01:06:40:DEBUG:root:Solution norm = 1.0464940956129567e-14\n", + "24-Aug 01:06:40:DEBUG:root:Converged in 6 iterations.\n", + "24-Aug 01:06:40:DEBUG:root:Displacement of element = \n", "[0.0 -0.39914506095474334 -0.0722858269559246 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n" ] @@ -708,7 +756,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 14, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -725,7 +773,7 @@ " mu = E/(2*(1+nu))\n", " la = E*nu/((1+nu)*(1-2*nu))\n", " la = 2*la*mu/(la + 2*mu)\n", - " set_attribute(el1, \"coordinates\", [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n", + " set_coordinates(el1, [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n", " set_attribute(el1, \"lambda\", la)\n", " set_attribute(el1, \"mu\", mu)\n", "\n", @@ -742,8 +790,10 @@ " ass.gdofs[el1.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n", " ass.gdofs[el2.id] = [5, 6]\n", "\n", - " for el in elements\n", - " assemble_element!(ass, el)\n", + " for (j, el) in enumerate(elements)\n", + " Logging.debug(\"Assembling element $j\")\n", + " assemble_lhs!(ass, el)\n", + " assemble_rhs!(ass, el)\n", " end\n", "\n", " # (Dirichlet) boundary conditions \"handled\"\n", @@ -782,13 +832,16 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "type MPC\n", + "abstract BoundaryCondition\n", + "abstract DirichletBC <: BoundaryCondition\n", + "\n", + "type MPC <: DirichletBC\n", " slave_dof :: Int64\n", " slave_value :: Float64\n", " master_dofs :: Array{Int64, 1}\n", @@ -799,7 +852,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -810,7 +863,7 @@ "MPC" ] }, - "execution_count": 16, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } @@ -826,7 +879,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -837,7 +890,7 @@ "create_ldof2gdofmap (generic function with 1 method)" ] }, - "execution_count": 17, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } @@ -870,9 +923,40 @@ "end" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we want to solve a _problem_. It's yet another container and if it's defined somewhat standard way default solver can solve it. Nothing stops user to write his/hers own solver." + ] + }, { "cell_type": "code", - "execution_count": 18, + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "type SolverParameters\n", + " max_iterations :: Int64\n", + " eps :: Float64\n", + "end\n", + "SolverParameters() = SolverParameters(10, 0.05)\n", + "\n", + "abstract Problem\n", + "abstract ElasticityProblem <: Problem\n", + "type PlaneStressProblem <: ElasticityProblem\n", + " elements :: Array{Element, 1}\n", + " dofmap :: Dict{Int64, Array{Int64,1}}() # a dict node_id : (dof1, dof2, ...)\n", + " dirichlet_bcs :: Array{BoundaryCondition, 1}\n", + " solver_parameters :: SolverParameters\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 22, "metadata": { "collapsed": false, "scrolled": false @@ -889,47 +973,53 @@ "name": "stderr", "output_type": "stream", "text": [ - "22-Aug 20:46:09:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n", - "22-Aug 20:46:09:INFO:root:solve!: dofs per node: 2\n", - "22-Aug 20:46:09:DEBUG:root:Problem size = 8\n", - "22-Aug 20:46:09:DEBUG:root:Starting iteration 1\n", - "22-Aug 20:46:09:DEBUG:root:Assembling\n", - "22-Aug 20:46:09:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 3.0900221367289444\n", - "22-Aug 20:46:10:DEBUG:root:Starting iteration 2\n", - "22-Aug 20:46:10:DEBUG:root:Assembling\n", - "22-Aug 20:46:10:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 0.32121316021534796\n", - "22-Aug 20:46:10:DEBUG:root:Starting iteration 3\n", - "22-Aug 20:46:10:DEBUG:root:Assembling\n", - "22-Aug 20:46:10:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 0.040431781940014504\n", - "22-Aug 20:46:10:DEBUG:root:Starting iteration 4\n", - "22-Aug 20:46:10:DEBUG:root:Assembling\n", - "22-Aug 20:46:10:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 0.0009291101052065917\n", - "22-Aug 20:46:10:DEBUG:root:Starting iteration 5\n", - "22-Aug 20:46:10:DEBUG:root:Assembling\n", - "22-Aug 20:46:10:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n", - "22-Aug 20:46:10:DEBUG:root:Starting iteration 6\n", - "22-Aug 20:46:10:DEBUG:root:Assembling\n", - "22-Aug 20:46:10:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:10:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "22-Aug 20:46:10:DEBUG:root:Solving system of equations. Total size = 12\n", - "22-Aug 20:46:10:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n", - "22-Aug 20:46:10:DEBUG:root:Converged in 6 iterations.\n", - "22-Aug 20:46:10:DEBUG:root:Displacement on upper right = \n", + "24-Aug 01:07:29:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n", + "24-Aug 01:07:29:INFO:root:solve!: dofs per node: 2\n", + "24-Aug 01:07:29:DEBUG:root:Problem size = 8\n", + "24-Aug 01:07:29:DEBUG:root:Starting iteration 1\n", + "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:29:DEBUG:root:Solution norm du = 3.0900221367289444\n", + "24-Aug 01:07:29:DEBUG:root:Starting iteration 2\n", + "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.32121316021534796\n", + "24-Aug 01:07:29:DEBUG:root:Starting iteration 3\n", + "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.040431781940014504\n", + "24-Aug 01:07:29:DEBUG:root:Starting iteration 4\n", + "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.0009291101052065917\n", + "24-Aug 01:07:30:DEBUG:root:Starting iteration 5\n", + "24-Aug 01:07:30:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:30:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:30:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:30:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:30:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:30:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n", + "24-Aug 01:07:30:DEBUG:root:Starting iteration 6\n", + "24-Aug 01:07:30:DEBUG:root:Assembling lhs\n", + "24-Aug 01:07:30:DEBUG:root:Assembling rhs\n", + "24-Aug 01:07:30:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "24-Aug 01:07:30:DEBUG:root:Added 4 Lagrange multipliers to model\n", + "24-Aug 01:07:30:DEBUG:root:Solving system of equations. Total size = 12\n", + "24-Aug 01:07:30:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n", + "24-Aug 01:07:30:DEBUG:root:Converged in 6 iterations.\n", + "24-Aug 01:07:30:DEBUG:root:Displacement on upper right = \n", "[-0.07228582695592467\n", " -2.222244754401765]\n" ] @@ -947,7 +1037,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 18, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -974,9 +1064,13 @@ " Logging.debug(\"Starting iteration $iter\")\n", " ass = JuliaFEM.Assembly(gdofs)\n", "\n", - " Logging.debug(\"Assembling\")\n", + " Logging.debug(\"Assembling lhs\")\n", " for el in elements\n", - " assemble_element!(ass, el)\n", + " assemble_lhs!(ass, el)\n", + " end\n", + " Logging.debug(\"Assembling rhs\")\n", + " for el in elements\n", + " assemble_rhs!(ass, el)\n", " end\n", "\n", " i = 0\n", @@ -1033,7 +1127,7 @@ " mu = E/(2*(1+nu))\n", " la = E*nu/((1+nu)*(1-2*nu))\n", " la = 2*la*mu/(la + 2*mu)\n", - " set_attribute(el1, \"coordinates\", [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n", + " set_coordinates(el1, [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n", " set_attribute(el1, \"lambda\", la)\n", " set_attribute(el1, \"mu\", mu)\n", "\n", @@ -1071,7 +1165,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -1080,17 +1174,17 @@ "name": "stderr", "output_type": "stream", "text": [ - "22-Aug 20:46:11:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", - "22-Aug 20:46:11:DEBUG:root:Found NODE section\n", - "22-Aug 20:46:11:DEBUG:root:Found ELEMENT section\n", - "22-Aug 20:46:11:DEBUG:root:120 elements found\n", - "22-Aug 20:46:12:INFO:root:Creating ELSET Body1\n", - "22-Aug 20:46:12:DEBUG:root:Found NSET section\n", - "22-Aug 20:46:12:DEBUG:root:Creating node set SUPPORT\n", - "22-Aug 20:46:12:DEBUG:root:Found NSET section\n", - "22-Aug 20:46:12:DEBUG:root:Creating node set LOAD\n", - "22-Aug 20:46:12:DEBUG:root:Found NSET section\n", - "22-Aug 20:46:12:DEBUG:root:Creating node set TOP\n" + "24-Aug 01:07:34:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n", + "24-Aug 01:07:34:DEBUG:root:Found NODE section\n", + "24-Aug 01:07:35:DEBUG:root:Found ELEMENT section\n", + "24-Aug 01:07:35:DEBUG:root:120 elements found\n", + "24-Aug 01:07:35:INFO:root:Creating ELSET Body1\n", + "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", + "24-Aug 01:07:35:DEBUG:root:Creating node set SUPPORT\n", + "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", + "24-Aug 01:07:35:DEBUG:root:Creating node set LOAD\n", + "24-Aug 01:07:35:DEBUG:root:Found NSET section\n", + "24-Aug 01:07:35:DEBUG:root:Creating node set TOP\n" ] }, { @@ -1103,7 +1197,7 @@ " \"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": 19, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } @@ -1123,13 +1217,15 @@ }, "outputs": [], "source": [ + "abstract Tet10 <: CG\n", + "\n", "\"\"\"\n", "Stress/displacement elements. 10-node quadratic tetrahedron.\n", "\"\"\"\n", - "type C3D10 <: ContinuumElement\n", + "type C3D10 <: Tet10\n", " id :: Int\n", " node_ids :: Array{Int, 1}\n", - " shape_functions :: FunctionSpace\n", + " coordinates :: Array{Float64, 2}\n", " integration_points :: Array{IntegrationPoint, 1}\n", " attributes :: Dict{ASCIIString, Any}\n", "end" @@ -1266,7 +1362,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 40, "metadata": { "collapsed": false, "scrolled": false @@ -1276,72 +1372,82 @@ "name": "stderr", "output_type": "stream", "text": [ - "22-Aug 20:46:13:DEBUG:root:Creating elements\n", - "22-Aug 20:46:13:DEBUG:root:Creating elements\n", - "22-Aug 20:46:13:INFO:root:solve!: dofs per node: 3\n", - "22-Aug 20:46:13:DEBUG:root:Problem size = 894\n", - "22-Aug 20:46:13:DEBUG:root:Starting iteration 1\n", - "22-Aug 20:46:13:DEBUG:root:Assembling\n", - "22-Aug 20:46:26:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:26:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:46:26:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:46:26:DEBUG:root:Solution norm du = 550.6462282437674\n", - "22-Aug 20:46:26:DEBUG:root:Starting iteration 2\n", - "22-Aug 20:46:26:DEBUG:root:Assembling\n", - "22-Aug 20:46:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:39:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:46:39:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:46:39:DEBUG:root:Solution norm du = 126.14054730775176\n", - "22-Aug 20:46:39:DEBUG:root:Starting iteration 3\n", - "22-Aug 20:46:39:DEBUG:root:Assembling\n", - "22-Aug 20:46:53:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:46:53:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:46:53:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:46:53:DEBUG:root:Solution norm du = 38.949840553368894\n", - "22-Aug 20:46:53:DEBUG:root:Starting iteration 4\n", - "22-Aug 20:46:54:DEBUG:root:Assembling\n", - "22-Aug 20:47:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:47:06:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:47:06:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:47:06:DEBUG:root:Solution norm du = 15.167069063650652\n", - "22-Aug 20:47:06:DEBUG:root:Starting iteration 5\n", - "22-Aug 20:47:06:DEBUG:root:Assembling\n", - "22-Aug 20:47:18:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:47:18:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:47:18:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:47:18:DEBUG:root:Solution norm du = 9.516311534304958\n", - "22-Aug 20:47:18:DEBUG:root:Starting iteration 6\n", - "22-Aug 20:47:18:DEBUG:root:Assembling\n", - "22-Aug 20:47:31:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:47:31:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:47:31:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:47:31:DEBUG:root:Solution norm du = 1.6222822043785954\n", - "22-Aug 20:47:31:DEBUG:root:Starting iteration 7\n", - "22-Aug 20:47:31:DEBUG:root:Assembling\n", - "22-Aug 20:47:43:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:47:43:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:47:43:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:47:43:DEBUG:root:Solution norm du = 0.09626397754176579\n", - "22-Aug 20:47:43:DEBUG:root:Starting iteration 8\n", - "22-Aug 20:47:43:DEBUG:root:Assembling\n", - "22-Aug 20:47:56:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:47:56:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:47:56:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:47:56:DEBUG:root:Solution norm du = 0.00026304537198068307\n", - "22-Aug 20:47:56:DEBUG:root:Starting iteration 9\n", - "22-Aug 20:47:56:DEBUG:root:Assembling\n", - "22-Aug 20:48:09:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:48:09:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:48:09:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:48:09:DEBUG:root:Solution norm du = 2.7245126807720425e-9\n", - "22-Aug 20:48:09:DEBUG:root:Starting iteration 10\n", - "22-Aug 20:48:09:DEBUG:root:Assembling\n", - "22-Aug 20:48:22:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "22-Aug 20:48:22:DEBUG:root:Added 27 Lagrange multipliers to model\n", - "22-Aug 20:48:22:DEBUG:root:Solving system of equations. Total size = 921\n", - "22-Aug 20:48:22:DEBUG:root:Solution norm du = 1.0919112177573261e-13\n", - "22-Aug 20:48:22:DEBUG:root:Converged in 10 iterations.\n", - "22-Aug 20:48:22:INFO:root:Maximum absolute displacement in y direction: 49.40459927455298\n" + "22-Aug 21:35:10:DEBUG:root:Creating elements\n", + "22-Aug 21:35:10:DEBUG:root:Creating elements\n", + "22-Aug 21:35:10:INFO:root:solve!: dofs per node: 3\n", + "22-Aug 21:35:10:DEBUG:root:Problem size = 894\n", + "22-Aug 21:35:10:DEBUG:root:Starting iteration 1\n", + "22-Aug 21:35:10:DEBUG:root:Assembling lhs\n", + "22-Aug 21:35:23:DEBUG:root:Assembling rhs\n", + "22-Aug 21:35:23:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:35:23:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:35:23:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:35:23:DEBUG:root:Solution norm du = 550.6462282437674\n", + "22-Aug 21:35:23:DEBUG:root:Starting iteration 2\n", + "22-Aug 21:35:23:DEBUG:root:Assembling lhs\n", + "22-Aug 21:35:36:DEBUG:root:Assembling rhs\n", + "22-Aug 21:35:36:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:35:36:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:35:36:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:35:36:DEBUG:root:Solution norm du = 126.14054730775176\n", + "22-Aug 21:35:36:DEBUG:root:Starting iteration 3\n", + "22-Aug 21:35:36:DEBUG:root:Assembling lhs\n", + "22-Aug 21:35:49:DEBUG:root:Assembling rhs\n", + "22-Aug 21:35:49:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:35:49:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:35:49:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:35:49:DEBUG:root:Solution norm du = 38.949840553368894\n", + "22-Aug 21:35:49:DEBUG:root:Starting iteration 4\n", + "22-Aug 21:35:49:DEBUG:root:Assembling lhs\n", + "22-Aug 21:36:03:DEBUG:root:Assembling rhs\n", + "22-Aug 21:36:04:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:36:04:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:36:04:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:36:04:DEBUG:root:Solution norm du = 15.167069063650652\n", + "22-Aug 21:36:04:DEBUG:root:Starting iteration 5\n", + "22-Aug 21:36:04:DEBUG:root:Assembling lhs\n", + "22-Aug 21:36:17:DEBUG:root:Assembling rhs\n", + "22-Aug 21:36:17:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:36:17:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:36:17:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:36:17:DEBUG:root:Solution norm du = 9.516311534304958\n", + "22-Aug 21:36:17:DEBUG:root:Starting iteration 6\n", + "22-Aug 21:36:17:DEBUG:root:Assembling lhs\n", + "22-Aug 21:36:32:DEBUG:root:Assembling rhs\n", + "22-Aug 21:36:32:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:36:32:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:36:32:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:36:32:DEBUG:root:Solution norm du = 1.6222822043785954\n", + "22-Aug 21:36:32:DEBUG:root:Starting iteration 7\n", + "22-Aug 21:36:32:DEBUG:root:Assembling lhs\n", + "22-Aug 21:36:46:DEBUG:root:Assembling rhs\n", + "22-Aug 21:36:48:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:36:48:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:36:48:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:36:48:DEBUG:root:Solution norm du = 0.09626397754176579\n", + "22-Aug 21:36:48:DEBUG:root:Starting iteration 8\n", + "22-Aug 21:36:48:DEBUG:root:Assembling lhs\n", + "22-Aug 21:37:07:DEBUG:root:Assembling rhs\n", + "22-Aug 21:37:07:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:37:07:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:37:07:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:37:07:DEBUG:root:Solution norm du = 0.00026304537198068307\n", + "22-Aug 21:37:07:DEBUG:root:Starting iteration 9\n", + "22-Aug 21:37:07:DEBUG:root:Assembling lhs\n", + "22-Aug 21:37:22:DEBUG:root:Assembling rhs\n", + "22-Aug 21:37:22:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:37:22:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:37:22:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:37:22:DEBUG:root:Solution norm du = 2.7245126807720425e-9\n", + "22-Aug 21:37:22:DEBUG:root:Starting iteration 10\n", + "22-Aug 21:37:22:DEBUG:root:Assembling lhs\n", + "22-Aug 21:37:35:DEBUG:root:Assembling rhs\n", + "22-Aug 21:37:35:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", + "22-Aug 21:37:35:DEBUG:root:Added 27 Lagrange multipliers to model\n", + "22-Aug 21:37:35:DEBUG:root:Solving system of equations. Total size = 921\n", + "22-Aug 21:37:35:DEBUG:root:Solution norm du = 1.0919112177573261e-13\n", + "22-Aug 21:37:35:DEBUG:root:Converged in 10 iterations.\n", + "22-Aug 21:37:35:INFO:root:Maximum absolute displacement in y direction: 49.40459927455298\n" ] } ], @@ -1914,15 +2020,6 @@ "source": [ "d.Image(\"/tmp/piston.png\")" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 016f0ce..4e12f3a 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -7,7 +7,10 @@ using Lexicon using Logging @Logging.configure(level=DEBUG) +Logging.info("loading types") include("types.jl") # type definitions +Logging.info("loading elements") +include("elements.jl") # elements include("math.jl") # basic mathematical operations include("elasticity_solver.jl") @@ -15,4 +18,6 @@ include("xdmf.jl") include("abaqus_reader.jl") include("interfaces.jl") +export set_coordinates, get_coordinates, set_material + end # module diff --git a/src/elements.jl b/src/elements.jl new file mode 100644 index 0000000..2eaef9c --- /dev/null +++ b/src/elements.jl @@ -0,0 +1,92 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +abstract CG <: Element # Lagrange element family +abstract Quad4 <: CG # 4 node quadrangle elements + + +""" +Evaluate basis functions in point xi. +""" +function get_basis(el::Quad4) + (xi) -> [(1-xi[1])*(1-xi[2])/4 + (1+xi[1])*(1-xi[2])/4 + (1+xi[1])*(1+xi[2])/4 + (1-xi[1])*(1+xi[2])/4] +end +function get_basis(el::Quad4, xi) + get_basis(el)(xi) +end + +""" +Evaluate partial derivatives of basis function w.r.t +dimensionless coordinate xi, i.e. dbasis/dxi +""" +function get_dbasisdxi(el::Quad4) + (xi) -> [-(1-xi[2])/4.0 -(1-xi[1])/4.0 + (1-xi[2])/4.0 -(1+xi[1])/4.0 + (1+xi[2])/4.0 (1+xi[1])/4.0 + -(1+xi[2])/4.0 (1-xi[1])/4.0] +end +function get_dbasisdxi(el::Quad4, xi) + get_dbasisdxi(el)(xi) +end + + + +""" +Get jacobian of element evaluated at point xi +""" +function get_jacobian(el::Element, xi) + dbasisdxi = get_dbasisdxi(el) + X = get_coordinates(el) + J = interpolate(X, dbasisdxi, xi)' + return J +end + +""" +Evaluate partial derivatives of basis function w.r.t +material description X, i.e. dbasis/dX +""" +function get_dbasisdX(el::CG) + function get_dbasisdX_(xi) + dbasisdxi = get_dbasisdxi(el, xi) + J = get_jacobian(el, xi) + dbasisdxi*inv(J) + end +end +function get_dbasisdX(el::CG, xi) + get_dbasisdX(el)(xi) +end + +""" +Return coordinates of element in array of size dim x nnodes +""" +function get_coordinates(el::Element) + # Make sure you define at least this field to your element if you want + # to build everything yourself + el.coordinates +end + +function set_coordinates(el::Element, coordinates) + el.coordinates = coordinates +end + +function set_material(el::Element, lambda, mu) + el.attributes["lambda"] = lambda + el.attributes["mu"] = mu +end + +""" +Get element id +""" +function get_element_id(el::Element) + el.id +end + +function get_integration_points(el::Element) + el.integration_points +end + + + diff --git a/src/math.jl b/src/math.jl index 0a962af..ba0ad99 100644 --- a/src/math.jl +++ b/src/math.jl @@ -60,26 +60,11 @@ function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip) return result end function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false) - return interpolate(e.attributes[field], derivative ? e.shape_functions.dbasis : e.shape_functions.basis, x) + basis = derivative ? get_dbasisdxi(e) : get_basis(e) + return interpolate(e.attributes[field], basis, x) end -""" - -""" -function get_basis(el::Element, xi) - return el.shape_functions.basis(xi) -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.shape_functions.dbasis(xi)*inv(J') - return dbasisdX -end - """ Linearize function f w.r.t some given field, i.e. calculate dR/du