From 2fb1587436c0e7fb42391bcd1b33e99ce3ec43a7 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 27 Aug 2015 08:18:20 +0300 Subject: [PATCH] Gets better all the time. --- ...2015-06-25-elasticity-solver-example.ipynb | 740 ++++++++---------- src/elements.jl | 11 +- src/math.jl | 64 +- src/types.jl | 3 +- 4 files changed, 413 insertions(+), 405 deletions(-) diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index a89aad2..e9208e6 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": 8, + "execution_count": 23, "metadata": { "collapsed": false }, @@ -47,15 +47,16 @@ "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" ] }, - "execution_count": 8, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ + "using ForwardDiff\n", "using JuliaFEM\n", - "using JuliaFEM: Assembly, Element, Quad4, get_jacobian, get_basis\n", - "using JuliaFEM: get_dbasisdxi, get_dbasisdX, get_field, set_field\n", + "using JuliaFEM: Assembly, Element, Equation, Point1, Quad4, IntegrationPoint\n", + "using JuliaFEM: get_field, set_field, get_element, linearize, integrate, interpolate, get_dbasisdX\n", "using Logging\n", "Logging.configure(level=DEBUG)" ] @@ -87,7 +88,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 24, "metadata": { "collapsed": false }, @@ -98,13 +99,12 @@ "CPS4" ] }, - "execution_count": 9, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "abstract Equation\n", "abstract Elasticity <: Equation\n", "\n", "\"\"\"\n", @@ -120,69 +120,11 @@ " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n", + " set_field(el, \"displacement\", zeros(2, 4))\n", " CPS4(el, integration_points)\n", "end" ] }, - { - "cell_type": "code", - "execution_count": 10, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "get_rhs (generic function with 2 methods)" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "using ForwardDiff\n", - "\n", - "get_integration_points(eq::Equation) = eq.integration_points\n", - "\n", - "\"\"\"\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(eq::Equation, f::Function)\n", - " target = []\n", - " for ip in get_integration_points(eq)\n", - " J = get_jacobian(eq.element, ip.xi)\n", - " push!(target, ip.weight*f(eq, ip)*det(J))\n", - " end\n", - " return sum(target)\n", - "end\n", - "\n", - "\n", - "\n", - "\"\"\"\n", - "Return left hand side of the equation Ax = b (i.e. A)\n", - "\"\"\"\n", - "function get_lhs(eq::Equation)\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(eq::Equation)\n", - " return None\n", - "end" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -194,7 +136,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 25, "metadata": { "collapsed": false }, @@ -205,55 +147,12 @@ "get_lhs (generic function with 2 methods)" ] }, - "execution_count": 15, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "#using JuliaFEM: linearize\n", - "\n", - "\"\"\"\n", - "Evaluate field in point xi using basis functions.\n", - "\"\"\"\n", - "function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})\n", - " f = get_field(el, field)\n", - " if !isa(f, Array)\n", - " # This is scalar, nothing to interpolate\n", - " return f\n", - " end\n", - " basis = get_basis(el, xi)\n", - " dim, nnodes = size(f)\n", - " result = zeros(dim)\n", - " for i=1:nnodes\n", - " result += basis[i]*f[:,i]\n", - " end\n", - " return result\n", - "end\n", - "\n", - "\n", - "StressEquation = Union()\n", - "\n", - "get_element(eq::Equation) = eq.element\n", - "\n", - "function linearize(f::Function, field::ASCIIString)\n", - " function jacobian(eq::Equation, args...)\n", - " el = get_element(eq)\n", - " fld = get_field(el, field)\n", - " dim, nnodes = size(fld)\n", - " function helper!(x, y)\n", - " orig = copy(fld)\n", - " set_field(el, field, reshape(x, dim, nnodes))\n", - " y[:] = f(eq, args...)\n", - " set_field(el, field, orig)\n", - " end\n", - " jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n", - " return jac(fld[:])\n", - " end\n", - " return jacobian\n", - "end\n", - "\n", - "\n", "\"\"\"\n", "Calculate internal nodal forces for continuum element.\n", "\"\"\"\n", @@ -271,11 +170,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", - " integrate(eq, (eq, ip) -> T(ip.xi, get_field(eq.element, \"displacement\")))\n", + " integrate(eq, (eq, ip) -> T(ip.xi, get_field(get_element(eq), \"displacement\")))\n", "end\n", "\n", "get_rhs(eq::CPS4) = -Wint(eq) # rhs = -R = -(T-F)\n", - "get_lhs(eq::CPS4) = linearize(Wint, \"displacement\")(eq)" + "get_lhs(eq::CPS4) = linearize(eq, Wint, \"displacement\")(eq)" ] }, { @@ -295,7 +194,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 26, "metadata": { "collapsed": true }, @@ -306,20 +205,26 @@ }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "test rhs\n", + "1 fact verified.\n" + ] + }, { "data": { "text/plain": [ - "CPS4(JuliaFEM.Quad4([1,2,3,4],Dict{ASCIIString,Any}(\"coordinates\"=>2x4 Array{Float64,2}:\n", - " 0.0 10.0 10.0 0.0\n", - " 0.0 0.0 1.0 1.0,\"lambda\"=>24.0,\"mu\"=>36.0)),[JuliaFEM.IntegrationPoint([-0.5773502691896258,-0.5773502691896258],1.0,Dict{ASCIIString,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896258,-0.5773502691896258],1.0,Dict{ASCIIString,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896258,0.5773502691896258],1.0,Dict{ASCIIString,Any}()),JuliaFEM.IntegrationPoint([-0.5773502691896258,0.5773502691896258],1.0,Dict{ASCIIString,Any}())])" + "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 28, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } @@ -344,59 +249,20 @@ " set_field(el, \"mu\", mu)\n", " return eq\n", "end\n", - "eq = get_test_equation()" + "facts(\"test rhs\") do\n", + " eq = get_test_equation()\n", + " utest = zeros(2, 4)\n", + " utest[1, 3] = 0.5\n", + " set_field(get_element(eq), \"displacement\", utest)\n", + " rhs = get_rhs(eq)\n", + " #println(rhs)\n", + " @fact rhs[2, 4] --> roughly(-8.4)\n", + "end" ] }, { "cell_type": "code", - "execution_count": 34, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "getindex (generic function with 132 methods)" - ] - }, - "execution_count": 34, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "# get_element(eq::Equation) = eq.element\n", - "#get_field(get_element(eq), \"coordinates\") += 1\n", - "#get_field(get_element(eq), \"coordinates\")\n", - "#getindex(eq::Equation, field) = eq.element.fields[field]\n", - "#setindex!(eq::Equation, Any, field) = " - ] - }, - { - "cell_type": "code", - "execution_count": 36, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "ename": "LoadError", - "evalue": "LoadError: MethodError: `setindex!` has no method matching setindex!(::CPS4, ::Array{Float64,2}, ::ASCIIString)\nwhile loading In[36], in expression starting on line 1", - "output_type": "error", - "traceback": [ - "LoadError: MethodError: `setindex!` has no method matching setindex!(::CPS4, ::Array{Float64,2}, ::ASCIIString)\nwhile loading In[36], in expression starting on line 1", - "" - ] - } - ], - "source": [ - "eq[\"coordinates\"] += 1.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, + "execution_count": 28, "metadata": { "collapsed": false, "scrolled": false @@ -413,40 +279,58 @@ "name": "stderr", "output_type": "stream", "text": [ - "26-Aug 01:36:30:DEBUG:root:Iteration 1\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 2\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 3\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 4\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 5\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 6\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Converged in 6 iterations.\n", - "26-Aug 01:36:30:DEBUG:root:solution vector: \n", + "27-Aug 08:16:52:DEBUG:root:Iteration 1\n", + "27-Aug 08:16:53:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Iteration 2\n", + "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Iteration 3\n", + "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Iteration 4\n", + "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Iteration 5\n", + "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Iteration 6\n", + "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:55:DEBUG:root:Converged in 6 iterations.\n", + "27-Aug 08:16:56:DEBUG:root:solution vector: \n", " [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n", - "26-Aug 01:36:30:DEBUG:root:norm of u: 3.1292483947150047\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 1\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 2\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 3\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 4\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 5\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Iteration 6\n", - "26-Aug 01:36:30:DEBUG:root:Solving Ax = b\n", - "26-Aug 01:36:30:DEBUG:root:Converged in 6 iterations.\n", - "26-Aug 01:36:30:DEBUG:root:solution vector: \n", + "27-Aug 08:16:56:DEBUG:root:norm of u: 3.1292483947150047\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 1\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 2\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 3\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 4\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 5\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Iteration 6\n", + "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n", + "27-Aug 08:16:56:DEBUG:root:Converged in 6 iterations.\n", + "27-Aug 08:16:56:DEBUG:root:solution vector: \n", " [0.0 1.2578327758133292 1.5202505368695098 0.0\n", - " 0.0 -1.8223091343697626 -1.6224781337179326 0.0]\n" + " 0.0 -1.8223091343697626 -1.6224781337179326 0.0]\n", + "27-Aug 08:16:56:DEBUG:root:norm of u: 3.129248394715004\n" ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2 facts verified.\n" + ] + }, + { + "data": { + "text/plain": [ + "delayed_handler (generic function with 4 methods)" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -454,7 +338,7 @@ "\n", " eq = get_test_equation()\n", " el = get_element(eq)\n", - " set_field(el, \"displacement\", zeros(2, 4))\n", + " #set_field(el, \"displacement\", zeros(2, 4))\n", " F = [0.0 0.0; 0.0 0.0; 0.0 -2.0; 0.0 0.0]'\n", "\n", " du = zeros(2, 4)\n", @@ -522,7 +406,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 29, "metadata": { "collapsed": false }, @@ -533,17 +417,19 @@ "assemble_rhs! (generic function with 1 method)" ] }, - "execution_count": 9, + "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function assemble_lhs!(ass::Assembly, el::Element)\n", + "function assemble_lhs!(ass::Assembly, eq::Equation)\n", "\n", - " gdofs = ass.gdofs[el.id]\n", + " el::Element = get_element(eq)\n", + " elid = get_field(el, \"id\")\n", + " gdofs = ass.gdofs[elid]\n", "\n", - " A = get_lhs(el)\n", + " A = get_lhs(eq)\n", " if !(A == None)\n", " ii, jj = size(A)\n", " for i=1:ii\n", @@ -556,11 +442,14 @@ " end\n", "end\n", "\n", - "function assemble_rhs!(ass::Assembly, el::Element)\n", + "function assemble_rhs!(ass::Assembly, eq::Equation)\n", + "\n", + " el::Element = get_element(eq)\n", + " elid = get_field(el, \"id\")\n", + " gdofs = ass.gdofs[elid]\n", "\n", - " gdofs = ass.gdofs[el.id]\n", " \n", - " b = get_rhs(el)\n", + " b = get_rhs(eq)\n", " if !(b == None)\n", " for i=1:length(b)\n", " push!(ass.i, gdofs[i])\n", @@ -572,7 +461,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 30, "metadata": { "collapsed": false }, @@ -580,122 +469,30 @@ { "data": { "text/plain": [ - "get_field (generic function with 2 methods)" + "get_lhs (generic function with 2 methods)" ] }, - "execution_count": 10, + "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function update_field(el::Element, du)\n", - " el.attributes[\"displacement\"][:] += du\n", - "end\n", - "function get_field(el::Element)\n", - " return el.attributes[\"displacement\"]\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Almost done. Also we must define interface how to set coordinates, material properties etc. for element." - ] - }, - { - "cell_type": "code", - "execution_count": 11, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "set_attribute (generic function with 1 method)" - ] - }, - "execution_count": 11, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function set_attribute(el::Element, field, value)\n", - " el.attributes[field] = value\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This might look a bit cumbersome at this point, but notice that everything is very general so far. We really don't have to define all this stuff when creating new elements if we follow some general construct. To demonstrate that, we define 0-dimensional \"point force\" element. Keep in mind, we need to have some functions defined how to operate using element, but they are already defined to superclass of element. So we need to only take care of construcor, rhs and (maybe) lhs. " - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": { - "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 <: Point0\n", - " id :: Int\n", - " node_ids :: Array{Int, 1}\n", - " attributes :: Dict{ASCIIString, Any}\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "get_lhs (generic function with 3 methods)" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "\"\"\"\n", - "Constructor of CPS1\n", - "\"\"\"\n", - "function CPS1(element_id, node_ids)\n", - " attributes = Dict(\n", - " \"displacement\" => zeros(2, 1),\n", - " \"displacement nodal load\" => zeros(2, 1))\n", - " CPS1(element_id, node_ids, attributes)\n", + "type CPS1 <: Elasticity\n", + " element :: Point1\n", + " integration_points :: Array{IntegrationPoint, 1}\n", "end\n", - "function get_rhs(el::CPS1)\n", - " return el.attributes[\"displacement nodal load\"]\n", + "function CPS1(el::Point1)\n", + " integration_points = []\n", + " set_field(el, \"displacement\", zeros(2, 1))\n", + " set_field(el, \"displacement nodal load\", zeros(2, 1))\n", + " CPS1(el, integration_points)\n", "end\n", - "function get_lhs(el::CPS1)\n", - " return None\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "And that's basically everything needs to be defined. During assembling get_lhs returns nothing and does not assemble anything to stiffness matrix (of course in case of follower point force the direction depents on the normal ..) and get_rhs adds simple point force pointing to some direction." + "get_rhs(eq::CPS1) = get_field(get_element(eq), \"displacement nodal load\")\n", + "get_lhs(eq::CPS1) = None" ] }, { @@ -710,7 +507,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 31, "metadata": { "collapsed": false }, @@ -719,45 +516,45 @@ "name": "stdout", "output_type": "stream", "text": [ - "one element assembly\n" + "two element assembly\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ - "24-Aug 18:33:33:DEBUG:root:Starting iteration 1\n", - "24-Aug 18:33:33:DEBUG:root:Assembling\n", - "24-Aug 18:33:33:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:33:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 3.0900221367289986\n", - "24-Aug 18:33:34:DEBUG:root:Starting iteration 2\n", - "24-Aug 18:33:34:DEBUG:root:Assembling\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 0.3212131602153472\n", - "24-Aug 18:33:34:DEBUG:root:Starting iteration 3\n", - "24-Aug 18:33:34:DEBUG:root:Assembling\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 0.04043178193999703\n", - "24-Aug 18:33:34:DEBUG:root:Starting iteration 4\n", - "24-Aug 18:33:34:DEBUG:root:Assembling\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 0.0009291101052105739\n", - "24-Aug 18:33:34:DEBUG:root:Starting iteration 5\n", - "24-Aug 18:33:34:DEBUG:root:Assembling\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 1.5638899027804743e-7\n", - "24-Aug 18:33:34:DEBUG:root:Starting iteration 6\n", - "24-Aug 18:33:34:DEBUG:root:Assembling\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 1\n", - "24-Aug 18:33:34:DEBUG:root:Assembling element 2\n", - "24-Aug 18:33:34:DEBUG:root:Solution norm = 1.0464940956129567e-14\n", - "24-Aug 18:33:34:DEBUG:root:Converged in 6 iterations.\n", - "24-Aug 18:33:34:DEBUG:root:Displacement of element = \n", + "27-Aug 08:17:00:DEBUG:root:Starting iteration 1\n", + "27-Aug 08:17:00:DEBUG:root:Assembling\n", + "27-Aug 08:17:00:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:00:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 3.0900221367289986\n", + "27-Aug 08:17:01:DEBUG:root:Starting iteration 2\n", + "27-Aug 08:17:01:DEBUG:root:Assembling\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.3212131602153472\n", + "27-Aug 08:17:01:DEBUG:root:Starting iteration 3\n", + "27-Aug 08:17:01:DEBUG:root:Assembling\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.04043178193999703\n", + "27-Aug 08:17:01:DEBUG:root:Starting iteration 4\n", + "27-Aug 08:17:01:DEBUG:root:Assembling\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.0009291101052105739\n", + "27-Aug 08:17:01:DEBUG:root:Starting iteration 5\n", + "27-Aug 08:17:01:DEBUG:root:Assembling\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 1.5638899027804743e-7\n", + "27-Aug 08:17:01:DEBUG:root:Starting iteration 6\n", + "27-Aug 08:17:01:DEBUG:root:Assembling\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n", + "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n", + "27-Aug 08:17:01:DEBUG:root:Solution norm = 1.0464940956129567e-14\n", + "27-Aug 08:17:01:DEBUG:root:Converged in 6 iterations.\n", + "27-Aug 08:17:01:DEBUG:root:Displacement of element = \n", "[0.0 -0.39914506095474334 -0.0722858269559246 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n" ] @@ -775,44 +572,47 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 14, + "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "facts(\"one element assembly\") do\n", + "facts(\"two element assembly\") do\n", " # set up element 1\n", - " element_id = 1\n", - " node_ids = [1, 2, 3, 4]\n", - " el1 = CPS4(element_id, node_ids)\n", + " el1 = Quad4([1, 2, 3, 4])\n", + " eq1 = CPS4(el1)\n", " # assign properties to element, e.g. coordinates, material properties, ...\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", - " 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", + " set_field(el1, \"coordinates\", [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n", + " set_field(el1, \"lambda\", la)\n", + " set_field(el1, \"mu\", mu)\n", + " set_field(el1, \"id\", 1)\n", "\n", " # set up element 2\n", - " el2 = CPS1(2, [3]) # Create nodal point force element with id 2 for node 3\n", - " set_attribute(el2, \"displacement nodal load\", [0.0, -2.0])\n", + " el2 = Point1([3])\n", + " eq2 = CPS1(el2) # Create nodal point force element with id 2 for node 3\n", + " set_field(el2, \"displacement nodal load\", [0.0, -2.0])\n", + " set_field(el2, \"id\", 2)\n", "\n", " elements = [el1, el2]\n", + " equations = [eq1, eq2]\n", "\n", " for i=1:10\n", " Logging.debug(\"Starting iteration $i\")\n", " Logging.debug(\"Assembling\")\n", " ass = Assembly()\n", - " ass.gdofs[el1.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n", - " ass.gdofs[el2.id] = [5, 6]\n", + " ass.gdofs[get_field(el1, \"id\")] = [1, 2, 3, 4, 5, 6, 7, 8]\n", + " ass.gdofs[get_field(el2, \"id\")] = [5, 6]\n", "\n", - " for (j, el) in enumerate(elements)\n", - " Logging.debug(\"Assembling element $j\")\n", - " assemble_lhs!(ass, el)\n", - " assemble_rhs!(ass, el)\n", + " for (j, eq) in enumerate(equations)\n", + " Logging.debug(\"Assembling lhs and rhs from equation $j\")\n", + " assemble_lhs!(ass, eq)\n", + " assemble_rhs!(ass, eq)\n", " end\n", "\n", " # (Dirichlet) boundary conditions \"handled\"\n", @@ -828,15 +628,18 @@ "\n", " # update solution back to elements\n", " for el in elements\n", - " eldu = du[ass.gdofs[el.id]]\n", - " update_field(el, eldu)\n", + " eldu = du[ass.gdofs[get_field(el, \"id\")]]\n", + " #tmp = get_field(el, \"displacement\")\n", + " #set_field(el, \"displacement\", tmp+eldu)\n", + " #update_field(el, eldu)\n", + " el.fields[\"displacement\"][:] += eldu\n", " end\n", " if norm(du) < 1.0e-9\n", " Logging.debug(\"Converged in $i iterations.\")\n", " break\n", " end\n", " end \n", - " disp = get_field(el1)\n", + " disp = get_field(el1, \"displacement\")\n", " Logging.debug(\"Displacement of element = \\n$disp\")\n", " @fact norm(disp) --> roughly(3.1292483947150043)\n", "end" @@ -851,7 +654,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -871,7 +674,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 13, "metadata": { "collapsed": false }, @@ -882,7 +685,7 @@ "MPC" ] }, - "execution_count": 16, + "execution_count": 13, "metadata": {}, "output_type": "execute_result" } @@ -896,6 +699,25 @@ "end" ] }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "abstract Problem\n", + "\n", + "type PlaneStressProblem <: Problem\n", + " equations :: Array{Equation, 1}\n", + " boundary_conditions :: Array{BoundaryCondition, 1}\n", + "\n", + "# dofmap :: Dict{Int64, Array{Int64,1}}() # a dict node_id : (dof1, dof2, ...)\n", + "# solver_parameters :: SolverParameters\n", + "end" + ] + }, { "cell_type": "code", "execution_count": 17, @@ -942,37 +764,6 @@ "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": 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": 18, @@ -2039,6 +1830,157 @@ "source": [ "d.Image(\"/tmp/piston.png\")" ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "using ForwardDiff" + ] + }, + { + "cell_type": "code", + "execution_count": 51, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);" + ] + }, + { + "cell_type": "code", + "execution_count": 52, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "5-element Array{Float64,1}:\n", + " 0.903617\n", + " 0.97582 \n", + " 0.36872 \n", + " 0.233261\n", + " 0.213665" + ] + }, + "execution_count": 52, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "x = rand(5)" + ] + }, + { + "cell_type": "code", + "execution_count": 53, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2.5477827220775753" + ] + }, + "execution_count": 53, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 54, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "5-element Array{Float64,1}:\n", + " 0.906762\n", + " 0.860405\n", + " 1.35153 \n", + " 1.59159 \n", + " 1.64713 " + ] + }, + "execution_count": 54, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g = ForwardDiff.gradient(f);\n", + "g(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 55, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Function" + ] + }, + "execution_count": 55, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "typeof(g)" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "LoadError", + "evalue": "LoadError: MethodError: `grad` has no method matching grad(::Function)\nwhile loading In[56], in expression starting on line 1", + "output_type": "error", + "traceback": [ + "LoadError: MethodError: `grad` has no method matching grad(::Function)\nwhile loading In[56], in expression starting on line 1", + "" + ] + } + ], + "source": [ + "j = ForwardDiff.jacobian(g)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/src/elements.jl b/src/elements.jl index 3371115..72ba279 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -3,13 +3,16 @@ abstract Element +get_element(eq::Equation) = eq.element + """ Get jacobian of element evaluated at point xi """ function get_jacobian(el::Element, xi) - dbasisdxi(xi) = get_dbasisdxi(el, xi) + dbasisdxi = get_dbasisdxi(el, xi) X = get_field(el, "coordinates") - J = interpolate(X, dbasisdxi, xi)' + #J = interpolate(X, dbasisdxi, xi)' + J = X*dbasisdxi return J end @@ -103,6 +106,10 @@ type Point1 <: CG node_ids :: Array{Int, 1} fields :: Dict{ASCIIString, Any} end +function Point1(node_ids) + fields = Dict{ASCIIString, Any}() + Point1(node_ids, fields) +end # 1d Lagrange elements diff --git a/src/math.jl b/src/math.jl index 642b380..6a4f0e5 100644 --- a/src/math.jl +++ b/src/math.jl @@ -1,9 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -""" -This module contains math stuff, including interpolation, integration, linearization, ... -""" +## This module contains math stuff, including interpolation, integration, linearization, ... using ForwardDiff @@ -185,3 +183,63 @@ function integrate!(f::Function, el::Element, target) end end +get_integration_points(eq::Equation) = eq.integration_points + +""" +Integrate f over element using Gaussian quadrature rules. + +Parameters +---------- +el::Element + well defined element +f::Function + Function to integrate +""" +function integrate(eq::Equation, f::Function) + target = [] + for ip in get_integration_points(eq) + J = get_jacobian(eq.element, ip.xi) + push!(target, ip.weight*f(eq, ip)*det(J)) + end + return sum(target) +end + + +""" +Evaluate field in point xi using basis functions. +""" +function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1}) + f = get_field(el, field) + if !isa(f, Array) + # This is scalar, nothing to interpolate + return f + end + basis = get_basis(el, xi) + dim, nnodes = size(f) + result = zeros(dim) + for i=1:nnodes + result += basis[i]*f[:,i] + end + return result +end + + +function linearize(eq::Equation, f::Function, field::ASCIIString) + function jacobian(eq::Equation, args...) + el = get_element(eq) + fld = get_field(el, field) + dim, nnodes = size(fld) + function helper(x::Vector) + orig = copy(fld) + set_field(el, field, reshape(x, dim, nnodes)) + y = f(eq, args...) + set_field(el, field, orig) + return y[:] + end + jac = ForwardDiff.jacobian(helper) + return jac(fld[:]) + end + return jacobian +end + + diff --git a/src/types.jl b/src/types.jl index 75c3255..c75a291 100644 --- a/src/types.jl +++ b/src/types.jl @@ -1,7 +1,8 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -export IntegrationPoint, Assembly +abstract Equation + """ Integration point