From f5f25a82c3133b4a2f9c48b60a4c7739a29a0f34 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 22 Oct 2015 07:00:13 +0300 Subject: [PATCH] elasticity equations --- ...2015-06-25-elasticity-solver-example.ipynb | 1174 +++++------------ src/interpolate.jl | 9 + src/types.jl | 36 + 3 files changed, 339 insertions(+), 880 deletions(-) diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index 1a77429..8a1c2d2 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -13,36 +13,18 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - " likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n" - ] - }, { "data": { "text/plain": [ "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)" ] }, - "execution_count": 1, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -50,14 +32,15 @@ "source": [ "using ForwardDiff\n", "using JuliaFEM: Quad4, Field, FieldSet, IntegrationPoint\n", - "using JuliaFEM: interpolate, get_element, get_dbasisdX\n", + "using JuliaFEM: interpolate, get_element, get_dbasisdX, dinterpolate\n", "using Logging\n", + "using FactCheck\n", "Logging.configure(level=DEBUG)" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -68,7 +51,7 @@ "CPS4" ] }, - "execution_count": 2, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -95,9 +78,18 @@ "end" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Geometrically non-linear Total Lagrange formulation\n", + "\n", + "### Method 1, tensor notation, automatic differentiation" + ] + }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -105,37 +97,68 @@ { "data": { "text/plain": [ - "has_lhs (generic function with 4 methods)" + "1x8 Array{Float64,2}:\n", + " -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0" ] }, - "execution_count": 11, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n", - " # just a quick test\n", + "function get_lhs_and_rhs(equation::CPS4, ip, time)\n", " element = get_element(equation)\n", - " dNdX = get_dbasisdX(element, ip, time)\n", - " E = interpolate(element, \"young\", ip, time)\n", - " nu = interpolate(element, \"poisson\", ip, time)\n", - " C = E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 1/2*(1-nu)]\n", - " B = zeros(3, 8)\n", - " for i=1:4\n", - " B[1, 2*(i-1)+1] = dNdX[i, 1]\n", - " B[2, 2*(i-1)+2] = dNdX[i, 2]\n", - " B[3, 2*(i-1)+1] = dNdX[i, 2]\n", - " B[3, 2*(i-1)+2] = dNdX[i, 1]\n", + " N = FEM.get_basis(element)\n", + " dN = FEM.diff(N)(ip.xi)\n", + "\n", + " # fields at time \n", + " X = element[\"geometry\"](time)\n", + " u = element[\"displacement\"](time)\n", + " young = element[\"young\"](time)\n", + " poisson = element[\"poisson\"](time)\n", + "\n", + " # material\n", + " young = interpolate(N, young, ip)\n", + " poisson = interpolate(N, poisson, ip)\n", + " mu = young/(2*(1+poisson))\n", + " lambda = young*poisson/((1+poisson)*(1-2*poisson))\n", + " lambda = 2*lambda*mu/(lambda + 2*mu)\n", + "\n", + " invJ = inv(dN*X)\n", + " gradw = dN*invJ\n", + "\n", + " function calc_R(data::Vector)\n", + " # create new field similar to displacement and fill it with data\n", + " f = similar(u, data)\n", + " # kinematics\n", + " gradu = dN*f*invJ\n", + " F = I + gradu # deformation gradient\n", + " E = 1/2*(F'*F - I) # strain\n", + " S = 2*mu*E + lambda*trace(E)*I # stress\n", + " return (F*S*gradw')[:]\n", " end\n", - " return B'*C*B\n", + "\n", + " R = calc_R(u[:])\n", + " Kt = ForwardDiff.jacobian(calc_R, u[:])\n", + " return Kt, -R\n", "end\n", - "JuliaFEM.has_lhs(equation::CPS4) = true" + "\n", + "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n", + " get_lhs_and_rhs(equation, ip, time)[1]\n", + "end\n", + "function JuliaFEM.get_rhs(equation::CPS4, ip, time)\n", + " get_lhs_and_rhs(equation, ip, time)[2]\n", + "end\n", + "\n", + "JuliaFEM.has_lhs(equation::CPS4) = true\n", + "JuliaFEM.has_rhs(equation::CPS4) = true\n", + "JuliaFEM.integrate_rhs(equation, 0.0)'" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -154,119 +177,95 @@ " 3.0 -319.4 15.0 -160.6 -3.0 158.8 -15.0 321.2" ] }, - "execution_count": 12, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], + "source": [ + "JuliaFEM.integrate_lhs(equation, 0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "increment 1, norm = 3.0900221367289094\n" + ] + }, + { + "data": { + "text/plain": [ + "2-element Array{Float64,1}:\n", + " -0.399145\n", + " -2.17799 " + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "increment 2, norm = 0.32121316021534246\n", + "increment 3, norm = 0.040431781940025155\n", + "increment 4, norm = 0.0009291101052125905\n", + "increment 5, norm = 1.5638899099729665e-7\n", + "increment 6, norm = 1.0201033584846351e-14\n", + "increment 7, norm = 1.3679680786658739e-15\n", + "increment 8, norm = 2.463346396682122e-15\n", + "increment 9, norm = 9.485831749489368e-16\n", + "increment 10, norm = 1.3372608112297287e-15\n" + ] + } + ], "source": [ "element = Quad4([1, 2, 3, 4])\n", "push!(element, FieldSet(\"young\", [Field(0.0, 90.0)]))\n", "push!(element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])]))\n", "push!(element, FieldSet(\"poisson\", [Field(0.0, 0.25)]))\n", "equation = CPS4(element)\n", - "JuliaFEM.integrate_lhs(equation, 1.0)" - ] - }, - { - "cell_type": "code", - "execution_count": 32, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "1x8 Array{Float64,2}:\n", - " 0.0 0.0 -0.158243 -2.17729 0.166432 -2.18057 0.0 0.0" - ] - }, - "execution_count": 32, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "F = [0.0, 0.0, 0.0, 0.0, 0.0, -2.0, 0.0, 0.0]\n", - "free_dofs = [3, 4, 5, 6]\n", - "u = zeros(8)\n", - "K = JuliaFEM.integrate_lhs(equation, 1.0)\n", - "u[free_dofs] = K[free_dofs, free_dofs] \\ F[free_dofs]\n", - "u'" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Potential energy:" - ] - }, - { - "cell_type": "code", - "execution_count": 35, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "-2.180569912309414" - ] - }, - "execution_count": 35, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "Wint(u::Vector) = 1/2*u'*K*u\n", - "Wext(u::Vector) = F'*u\n", - "Π(u::Vector) = (Wint(u) - Wext(u))[1]\n", - "Π(u)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "On the other hand:" - ] - }, - { - "cell_type": "code", - "execution_count": 49, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "1x8 Array{Float64,2}:\n", - " 0.0 0.0 -0.158243 -2.17729 0.166432 -2.18057 0.0 0.0" - ] - }, - "execution_count": 49, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "r = ForwardDiff.gradient(Π)\n", - "J = ForwardDiff.jacobian(r)\n", + "\n", + "u0 = Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n", + "push!(element[\"displacement\"], u0)\n", "u = zeros(8)\n", "du = zeros(8)\n", - "du[free_dofs] = J(u)[free_dofs, free_dofs] \\ -r(u)[free_dofs]\n", - "u += du\n", - "u'" + "fd = [3, 4, 5, 6]\n", + "f = zeros(8)\n", + "f[6] = -2.0\n", + "for i=1:10\n", + " A = JuliaFEM.integrate_lhs(equation, 1.0)\n", + " b = JuliaFEM.integrate_rhs(equation, 1.0)\n", + " du[fd] = A[fd,fd] \\ (b[fd]+f[fd])\n", + " u += du\n", + " new_field = similar(u0, u)\n", + " new_field.time = 1.0\n", + " new_field.increment = i\n", + " push!(element[\"displacement\"], new_field)\n", + " println(\"increment $i, norm = $(norm(du))\")\n", + "end\n", + "interpolate(element, \"displacement\", [1.0, -1.0], 1.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Method 2, Voigt notation, analytical linearization" ] }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 7, "metadata": { "collapsed": false }, @@ -274,269 +273,91 @@ { "data": { "text/plain": [ - "get_lhs (generic function with 2 methods)" + "1x8 Array{Float64,2}:\n", + " -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0 -0.0" ] }, - "execution_count": 25, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "\"\"\"\n", - "Calculate internal nodal forces for continuum element.\n", - "\"\"\"\n", - "function Wint(eq::CPS4)\n", - " el = eq.element\n", - " dNdX(xi) = get_dbasisdX(el, xi)\n", + "function get_lhs_and_rhs(equation::CPS4, ip, time)\n", + " element = get_element(equation)\n", + "\n", + " # fields\n", + " X = element[\"geometry\"](time)\n", + " u = element[\"displacement\"](time)\n", + " young = element[\"young\"](time)\n", + " poisson = element[\"poisson\"](time)\n", + "\n", " # material\n", - " lambda(xi) = interpolate(el, \"lambda\", xi)\n", - " mu(xi) = interpolate(el, \"mu\", xi)\n", + " N = FEM.get_basis(element)\n", + " young = interpolate(N, young, ip)\n", + " poisson = interpolate(N, poisson, ip)\n", + "\n", + " dN = FEM.diff(N)(ip.xi)\n", + " invJ = inv(dN*X)\n", + " dNdX = dN*invJ\n", + "\n", " # kinematics\n", - " Grad(xi, u) = u*dNdX(xi)\n", - " F(xi, u) = I + Grad(xi, u)\n", - " E(xi, u) = 1/2*(Grad(xi, u)' + Grad(xi, u) + Grad(xi, u)'*Grad(xi, u))\n", - " # constitutive equation\n", - " 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(get_element(eq), \"displacement\")))\n", + " gradu = dN*u*invJ\n", + " F = I + gradu # deformation gradient\n", + " E = 1/2*(F'*F - I) # GL strain tensor\n", + " E = [E[1,1], E[2,2], E[1,2]*2] # go to Voigt\n", + "\n", + " # constitutive equations\n", + " D = young/(1-poisson^2) * [1 poisson 0; poisson 1 0; 0 0 1/2*(1-poisson)]\n", + " S = D*E\n", + " T = zeros(4, 4)\n", + " T[1,1] = S[1]\n", + " T[2,2] = S[2]\n", + " T[1,2] = T[2,1] = S[3]\n", + " T[3:4,3:4] = T[1:2,1:2]\n", + "\n", + " # linear part\n", + " B_L = zeros(3, 8)\n", + " for i=1:4\n", + " B_L[1, 2*(i-1)+1] = F[1,1]*dNdX[i,1]\n", + " B_L[1, 2*(i-1)+2] = F[2,1]*dNdX[i,1]\n", + " B_L[2, 2*(i-1)+1] = F[1,2]*dNdX[i,2]\n", + " B_L[2, 2*(i-1)+2] = F[2,2]*dNdX[i,2]\n", + " B_L[3, 2*(i-1)+1] = F[1,1]*dNdX[i,2] + F[1,2]*dNdX[i,1]\n", + " B_L[3, 2*(i-1)+2] = F[2,1]*dNdX[i,2] + F[2,2]*dNdX[i,1]\n", + " end\n", + " K_L = B_L'*D*B_L\n", + "\n", + " # nonlinear part\n", + " B_NL = zeros(4, 8)\n", + " for i=1:4\n", + " B_NL[1, 2*(i-1)+1] = dNdX[i,1]\n", + " B_NL[2, 2*(i-1)+1] = dNdX[i,2]\n", + " B_NL[3, 2*(i-1)+2] = dNdX[i,1]\n", + " B_NL[4, 2*(i-1)+2] = dNdX[i,2]\n", + " end\n", + " K_NL = B_NL'*T*B_NL\n", + "\n", + " fint = B_L'*S\n", + "\n", + " R = fint\n", + " Kt = K_L + K_NL\n", + " return Kt, -R\n", "end\n", "\n", - "get_rhs(eq::CPS4) = -Wint(eq) # rhs = -R = -(T-F)\n", - "get_lhs(eq::CPS4) = linearize(eq, Wint, \"displacement\")(eq)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "That was our geometrically nonlinear elasticity solver. Note how we used automatic differentiation to linearize residual vector.\n", - "\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Design principle 6*: we test our code. We use FactCheck for testing." - ] - }, - { - "cell_type": "code", - "execution_count": 26, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "using FactCheck" - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "test rhs\n", - "1 fact verified.\n" - ] - }, - { - "data": { - "text/plain": [ - "delayed_handler (generic function with 4 methods)" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function get_test_equation()\n", - " # set up one linear quadrangle element\n", - " node_ids = [1, 2, 3, 4]\n", - " el = Quad4(node_ids)\n", - " eq = CPS4(el)\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", - " 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", - " set_field(el, \"coordinates\", X)\n", - " set_field(el, \"lambda\", la)\n", - " set_field(el, \"mu\", mu)\n", - " return eq\n", + "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n", + " get_lhs_and_rhs(equation, ip, time)[1]\n", "end\n", - "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" + "function JuliaFEM.get_rhs(equation::CPS4, ip, time)\n", + " get_lhs_and_rhs(equation, ip, time)[2]\n", + "end\n", + "\n", + "JuliaFEM.integrate_rhs(equation, 0.0)'" ] }, { "cell_type": "code", - "execution_count": 28, - "metadata": { - "collapsed": false, - "scrolled": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "test solve one element model\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "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", - "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", - "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": [ - "facts(\"test solve one element model\") do\n", - "\n", - " eq = get_test_equation()\n", - " el = get_element(eq)\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", - "\n", - " free_dofs = [3, 4, 5, 6]\n", - " for i=1:10\n", - " Logging.debug(\"Iteration $i\")\n", - " b = get_rhs(eq)\n", - " A = get_lhs(eq)\n", - " Logging.debug(\"Solving Ax = b\")\n", - " du[free_dofs] = A[free_dofs, free_dofs] \\ (b + F)[free_dofs]\n", - "\n", - " field = get_field(el, \"displacement\") + du\n", - " set_field(el, \"displacement\", field)\n", - " if norm(du) < 1.0e-9\n", - " Logging.debug(\"Converged in $i iterations.\")\n", - " break\n", - " end\n", - " end\n", - "\n", - " # Tested against Elmer solution\n", - " u = get_field(el, \"displacement\")\n", - " Logging.debug(\"solution vector: \\n $u\")\n", - " @fact u[2, 3] --> roughly(-2.222244754401764)\n", - " norm1 = norm(u)\n", - " Logging.debug(\"norm of u: $(norm(u))\")\n", - "\n", - " # We rotate model a bit and make sure that norm remains same\n", - " phi = 45/180*pi\n", - " rmat = [\n", - " cos(phi) -sin(phi)\n", - " sin(phi) cos(phi)]\n", - " set_field(el, \"coordinates\", rmat*get_field(el, \"coordinates\"))\n", - " F = rmat*F\n", - "\n", - " set_field(el, \"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", - " Logging.debug(\"Iteration $i\")\n", - " A = get_lhs(eq)\n", - " b = get_rhs(eq)\n", - " Logging.debug(\"Solving Ax = b\")\n", - " du[free_dofs] = A[free_dofs, free_dofs] \\ (b + F)[free_dofs]\n", - "\n", - " field = get_field(el, \"displacement\") + du\n", - " set_field(el, \"displacement\", field)\n", - " if norm(du) < 1.0e-9\n", - " Logging.debug(\"Converged in $i iterations.\")\n", - " break\n", - " end\n", - " end\n", - " u = get_field(el, \"displacement\")\n", - " Logging.debug(\"solution vector: \\n $u\")\n", - " Logging.debug(\"norm of u: $(norm(u))\")\n", - " @fact norm(u) --> roughly(norm1) \n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Assembly procedure is now very general because we always just have to call `get_lhs` and `get_rhs` to get corresponding vectors and matrices from element." - ] - }, - { - "cell_type": "code", - "execution_count": 29, + "execution_count": 8, "metadata": { "collapsed": false }, @@ -544,49 +365,24 @@ { "data": { "text/plain": [ - "assemble_rhs! (generic function with 1 method)" + "8x8 Array{Float64,2}:\n", + " 123.2 15.0 56.8 -3.0 -61.6 -15.0 -118.4 3.0\n", + " 15.0 321.2 3.0 158.8 -15.0 -160.6 -3.0 -319.4\n", + " 56.8 3.0 123.2 -15.0 -118.4 -3.0 -61.6 15.0\n", + " -3.0 158.8 -15.0 321.2 3.0 -319.4 15.0 -160.6\n", + " -61.6 -15.0 -118.4 3.0 123.2 15.0 56.8 -3.0\n", + " -15.0 -160.6 -3.0 -319.4 15.0 321.2 3.0 158.8\n", + " -118.4 -3.0 -61.6 15.0 56.8 3.0 123.2 -15.0\n", + " 3.0 -319.4 15.0 -160.6 -3.0 158.8 -15.0 321.2" ] }, - "execution_count": 29, + "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function assemble_lhs!(ass::Assembly, eq::Equation)\n", - "\n", - " el::Element = get_element(eq)\n", - " elid = get_field(el, \"id\")\n", - " gdofs = ass.gdofs[elid]\n", - "\n", - " A = get_lhs(eq)\n", - " if !(A == None)\n", - " ii, jj = size(A)\n", - " for i=1:ii\n", - " for j=1:jj\n", - " push!(ass.I, gdofs[i])\n", - " push!(ass.J, gdofs[j])\n", - " push!(ass.A, A[i,j])\n", - " end\n", - " end\n", - " end\n", - "end\n", - "\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", - " \n", - " b = get_rhs(eq)\n", - " if !(b == None)\n", - " for i=1:length(b)\n", - " push!(ass.i, gdofs[i])\n", - " push!(ass.b, b[i])\n", - " end\n", - " end\n", - "end" + "JuliaFEM.integrate_lhs(equation, 0.0)" ] }, { @@ -635,465 +431,6 @@ " 0.0 -2.1779892317073504 -2.222244754401764 0.0]" ] }, - { - "cell_type": "code", - "execution_count": 31, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "two element assembly\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "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" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1 fact verified.\n" - ] - }, - { - "data": { - "text/plain": [ - "delayed_handler (generic function with 4 methods)" - ] - }, - "execution_count": 31, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "facts(\"two element assembly\") do\n", - " # set up element 1\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_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 = 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[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, 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", - " free_dofs = [3, 4, 5, 6]\n", - "\n", - " # solution\n", - " A = sparse(ass.I, ass.J, ass.A)\n", - " b = full(sparsevec(ass.i, ass.b))\n", - " du = zeros(8)\n", - " du[free_dofs] = A[free_dofs, free_dofs] \\ b[free_dofs]\n", - "\n", - " Logging.debug(\"Solution norm = $(norm(du))\")\n", - "\n", - " # update solution back to elements\n", - " for el in elements\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, \"displacement\")\n", - " Logging.debug(\"Displacement of element = \\n$disp\")\n", - " @fact norm(disp) --> roughly(3.1292483947150043)\n", - "end" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Seems to be working. But we still need to handle Dirichlet boundary conditions somewhat more generally. In very general form Dirichlet bc can be expressed as $\\mathbf{B}\\mathbf{u} = \\mathbf{d}$ for variational problem and $\\mathbf{B}\\mathbf{u} \\leq \\mathbf{d}$ for variational inequality problems. In contact mechanics typically nodes are divided to several sets, one with slave nodes (can be eliminated), master nodes, and all other nodes. $\\mathbf{B}$ is usually something $\\mathbf{B}=\\begin{bmatrix}\\mathbf{0} & \\mathbf{D} & -\\mathbf{M}\\end{bmatrix}^\\mathrm{T}$. In normal Dirichlet boundary condition this simplifies to something $\\mathbf{D}\\mathbf{d}_\\mathcal{S} = \\mathbf{0}$, where nodes in set $\\mathcal{S}$ are known to be slave nodes. Finally, contact is nothing more than multi point constraint. Dirichlet boundary conditions can be easily eliminated if $\\mathbf{D}$ is diagonal: $\\mathbf{D}\\mathbf{d}_\\mathcal{S} = \\mathbf{M}\\mathbf{d}_\\mathcal{M} \\Rightarrow \\mathbf{d}_\\mathcal{S} = \\mathbf{D}^{-1}\\mathbf{M}\\mathbf{d}_\\mathcal{M} = \\mathbf{P}\\mathbf{d}_\\mathcal{M}$." - ] - }, - { - "cell_type": "code", - "execution_count": 12, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "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", - " master_values :: Array{Float64, 1}\n", - " constant :: Float64\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 13, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "MPC" - ] - }, - "execution_count": 13, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "\"\"\"\n", - "Convenient function, to set some dof=0\n", - "\"\"\"\n", - "function MPC(dof)\n", - " MPC(dof, 1.0, Int64[], Float64[], 0.0)\n", - "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, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "create_ldof2gdofmap (generic function with 1 method)" - ] - }, - "execution_count": 17, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "\"\"\"\n", - "Create local dof to global dof mapping for given elements\n", - "\"\"\"\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", - " 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", - "\n", - " # Assign global dof for each node\n", - " pdim = 1\n", - " ngdofs = Dict{Int64, Array{Int64,1}}()\n", - " for nid in all_node_ids\n", - " ngdofs[nid] = collect(pdim:pdim+ndofs-1)\n", - " pdim += ndofs\n", - " end\n", - "\n", - " return ngdofs\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 18, - "metadata": { - "collapsed": false, - "scrolled": false - }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "solve one element problem\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "24-Aug 18:34:38:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n", - "24-Aug 18:34:38:INFO:root:solve!: dofs per node: 2\n", - "24-Aug 18:34:39:DEBUG:root:Problem size = 8\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 1\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 3.0900221367289444\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 2\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.32121316021534796\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 3\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.040431781940014504\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 4\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.0009291101052065917\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 5\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n", - "24-Aug 18:34:39:DEBUG:root:Starting iteration 6\n", - "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n", - "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n", - "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n", - "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n", - "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n", - "24-Aug 18:34:39:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n", - "24-Aug 18:34:39:DEBUG:root:Converged in 6 iterations.\n", - "24-Aug 18:34:39:DEBUG:root:Displacement on upper right = \n", - "[-0.07228582695592467\n", - " -2.222244754401765]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "1 fact verified.\n" - ] - }, - { - "data": { - "text/plain": [ - "delayed_handler (generic function with 4 methods)" - ] - }, - "execution_count": 18, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function solve!(elements, dofmap, dirichlet_bcs; ndofs=2, max_iterations=10)\n", - "\n", - " Logging.info(\"solve!: dofs per node: $ndofs\")\n", - " pdim = length(dofmap)*ndofs\n", - " Logging.debug(\"Problem size = $pdim\")\n", - "\n", - " # Assign global dofs for element ids\n", - " gdofs = Dict{Int64, Array{Int64,1}}()\n", - " for el in elements\n", - " gdofs[el.id] = Int64[]\n", - " for nid in el.node_ids\n", - " for ndof in dofmap[nid]\n", - " push!(gdofs[el.id], ndof)\n", - " end\n", - " end\n", - " end\n", - "\n", - " for iter=1:max_iterations\n", - " Logging.debug(\"Starting iteration $iter\")\n", - " ass = JuliaFEM.Assembly(gdofs)\n", - "\n", - " Logging.debug(\"Assembling lhs\")\n", - " for el in elements\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", - " Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n", - " # Dirichlet boundary conditions (this leads to a saddle point problem)\n", - " for bc in dirichlet_bcs\n", - " i += 1\n", - " #Logging.debug(\"lock dof $(bc.slave_dof), matrix row $(pdim+i)\")\n", - " push!(ass.I, bc.slave_dof)\n", - " push!(ass.J, pdim+i)\n", - " push!(ass.A, bc.slave_value)\n", - " push!(ass.I, pdim+i)\n", - " push!(ass.J, bc.slave_dof)\n", - " push!(ass.A, bc.slave_value)\n", - " push!(ass.i, pdim+i)\n", - " push!(ass.b, bc.constant)\n", - " end\n", - " Logging.debug(\"Added $i Lagrange multipliers to model\")\n", - " Logging.debug(\"Solving system of equations. Total size = $(pdim+i)\")\n", - "\n", - " # solution\n", - " A = sparse(ass.I, ass.J, ass.A)\n", - " b = full(sparsevec(ass.i, ass.b))\n", - " du = A \\ b\n", - "\n", - " solnorm = norm(du[1:pdim])\n", - " Logging.debug(\"Solution norm du = $solnorm\")\n", - "\n", - " # update solution back to elements\n", - " for el in elements\n", - " eldu = du[ass.gdofs[el.id]]\n", - " update_field(el, eldu)\n", - "\n", - " end\n", - " if solnorm < 1.0e-9\n", - " Logging.debug(\"Converged in $iter iterations.\")\n", - " break\n", - " end\n", - " end\n", - "\n", - "end\n", - "\n", - "ENV[\"COLUMNS\"] = 160\n", - "\n", - "facts(\"solve one element problem\") do\n", - "\n", - " # set up element 1\n", - " element_id = 1\n", - " node_ids = [1, 2, 3, 4]\n", - " el1 = CPS4(element_id, node_ids)\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", - "\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", - "\n", - " elements = [el1, el2]\n", - "\n", - " dofmap = create_ldof2gdofmap(elements, \"displacement\")\n", - " Logging.debug(dofmap)\n", - " # Boundary conditions\n", - " # dirichlet bc, set dx=dy=0 on support\n", - " mpc1 = MPC(dofmap[1][1]) # node 1, dx=0\n", - " mpc2 = MPC(dofmap[1][2]) # node 1, dy\n", - " mpc3 = MPC(dofmap[4][1]) # node 4, dx\n", - " mpc4 = MPC(dofmap[4][2]) # node 4, dy\n", - " dbcs = [mpc1, mpc2, mpc3, mpc4]\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, dbcs; max_iterations=10)\n", - " Logging.debug(\"Displacement on upper right = \\n$(get_field(el2))\")\n", - " disp = get_field(el1)\n", - " @fact norm(disp) --> roughly(3.1292483947150043)\n", - "end\n" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -1465,6 +802,83 @@ "cell_type": "markdown", "metadata": {}, "source": [ + " 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", + "\n", "22.8.2015 13-14 seconds/assembly." ] }, diff --git a/src/interpolate.jl b/src/interpolate.jl index 721d7b2..b156823 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -51,7 +51,16 @@ function interpolate(fields::FieldSet, t::Number) return f end +function call(fieldset::FieldSet, time::Number) + interpolate(fieldset, time) +end + +function interpolate(basis::Basis, field::Field, ip::IntegrationPoint) + interpolate(basis, field, ip.xi) +end + function dinterpolate(N::Basis, u::Field, xi::Array{Float64, 1}) dN = diff(N) dN(xi)*u end + diff --git a/src/types.jl b/src/types.jl index 82fa169..1a2a93b 100644 --- a/src/types.jl +++ b/src/types.jl @@ -47,6 +47,42 @@ function Base.(:+)(f1::Field, f2::Field) Field(f1.time, f1.values + f2.values) end +""" Return data from field as a long array. + +Examples +-------- +>>> f = Field(0.0, Vector[[1.0, 2.0], [3.0, 4.0]]) +>>> f[:] +[1.0, 2.0, 3.0, 4.0] + +""" +function Base.getindex(field::Field, c::Colon) + [field.values...;] +end + +""" Return field similar to input but with new data in it. + +Examples +-------- +>>> f = Field(0.5, Vector[[1.0, 2.0], [3.0, 4.0]]) +>>> similar(f, ones(4)) +JuliaFEM.Field{Array{Array{T,1},1}}(0.5,1,Array{T,1}[[1.0,1.0],[1.0,1.0]]) + +""" +function Base.similar(field::Field, data::Vector) + new_field = Field(field.time, similar(field.values)) + data = reshape(data, round(Int, length(data)/length(field)), length(field)) + for i=1:length(new_field) + new_field.values[i] = data[:,i] + end + new_field +end + + + + + + """ FieldSet is set of fields, each field can have different time and/or increment. """