diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index b27ccad..02f5773 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -4,45 +4,39 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Solving elasticity problems using JuliaFEM\n", + "# Solving elasticity equations using JuliaFEM\n", "\n", "Author(s): Jukka Aho\n", "\n", - "**Abstract**: Elasticity equations design notes." + "**Abstract**: Elasticity equations design notes.\n", + "\n", + "##Weak form\n", + "\n", + "Given function spaces\n", + "\\begin{align}\n", + "\\boldsymbol{\\mathcal{U}} & =\\left\\{ \\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\boldsymbol{u}\\left(\\boldsymbol{X},t\\right)=\\hat{\\boldsymbol{u}}\\left(\\boldsymbol{X},t\\right)\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\\\\\n", + "\\boldsymbol{\\mathcal{V}} & =\\left\\{ \\delta\\boldsymbol{u}\\in H^{1}\\left(\\Omega\\right)|\\delta\\boldsymbol{u}\\left(\\boldsymbol{X}\\right)=0\\text{ on }\\Gamma_{\\mathrm{u}}\\right\\} ,\n", + "\\end{align}\n", + "find $\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{U}}$ such that\n", + "\\begin{equation}\n", + "\\delta\\mathcal{W}:=\\int_{\\Omega_{0}}\\rho_{0}\\ddot{\\boldsymbol{u}}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}+\\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\hat{\\boldsymbol{b}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}-\\int_{\\Gamma_{\\sigma}}\\hat{\\boldsymbol{t}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}A_{0} =0 \\qquad\\forall\\delta\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{V}}\n", + "\\end{equation}" ] }, { "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" } @@ -58,7 +52,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 3, "metadata": { "collapsed": false }, @@ -69,7 +63,7 @@ "CPS4" ] }, - "execution_count": 2, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -107,7 +101,7 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -115,79 +109,58 @@ { "data": { "text/plain": [ - "8-element Array{Float64,1}:\n", - " -189.267 \n", - " 26.1073 \n", - " 5.7554e-13 \n", - " 2.6894e-12 \n", - " -8.73968e-13\n", - " 20.0 \n", - " 189.267 \n", - " -46.1073 " + "has_rhs (generic function with 4 methods)" ] }, - "execution_count": 40, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "function get_lhs_and_rhs(equation::CPS4, ip, time)\n", + " # boilerplate code start\n", " element = get_element(equation)\n", - " N = FEM.get_basis(element)\n", - " dN = FEM.diff(N)(ip.xi)\n", + " geometry = element[\"geometry\"](time)\n", + " basis = FEM.get_basis(element)\n", + " dbasis = FEM.diff(basis)(ip.xi)\n", + " grad(u) = dbasis*u*inv(dbasis*geometry)\n", + " # boilerplate code end -- replace with a macro?\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", + " # interpolate fields in temporal dimension\n", + " young = element[\"youngs modulus\"](time)\n", + " poisson = element[\"poissons ratio\"](time)\n", + " displacement = element[\"displacement\"](time)\n", "\n", - " # material\n", - " young = interpolate(N, young, ip)\n", - " poisson = interpolate(N, poisson, ip)\n", + " # interpolate material in spatial dimension\n", + " young = interpolate(basis, young, ip)\n", + " poisson = interpolate(basis, 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", + " lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d\n", "\n", - " gradw = dN*inv(dN*X)\n", - " grad(u) = dN*u*inv(dN*(X+u))\n", - " Grad(u) = dN*u*inv(dN*X)\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", - " F = I + Grad(f) # deformation gradient\n", - " C = F'*F\n", - " E = 1/2*(C - I) # strain\n", + " function W(data::Vector)\n", + " Wint = 0.0\n", + " # create new field u, similar to field displacement, and fill it with data\n", + " u = similar(displacement, data)\n", + " F = I + grad(u) # deformation gradient\n", + " E = 1/2*(F'*F - I) # strain\n", " S = 2*mu*E + lambda*trace(E)*I # stress\n", - " P = F*S\n", - " return (P*gradw')[:]\n", + " Wint += 1/2*trace(S*E')\n", + "\n", + " Wext = 0.0\n", + " # any volume load?\n", + " if haskey(element, \"displacement volume load\")\n", + " b = interpolate(element, \"displacement volume load\", ip, time)\n", + " δu = interpolate(basis, u, ip)\n", + " Wext += dot(b, δu)\n", + " end\n", + " return Wint - Wext\n", " end\n", "\n", - " function calc_Wint(data::Vector)\n", - " # create new field similar to displacement and fill it with data\n", - " δu = similar(u, data)\n", - " # kinematics\n", - " F = I + Grad(δu) # deformation gradient\n", - " C = F'*F\n", - " E = 1/2*(C - I) # strain\n", - " #E = 1/2*(Grad(δu)' + Grad(δu) + Grad(δu)'*Grad(δu))\n", - " S = 2*mu*E + lambda*trace(E)*I # stress\n", - " δF = Grad(δu)\n", - " #δɛ = 1/2*(δF' + δF)\n", - " #δE = F'*δɛ*F\n", - " #δE = 1/2*(F'*δF + δF'*F)\n", - " #δE = 1/2*(Grad(δu)' + Grad(δu) + Grad(δu)'*Grad(δu))\n", - " δE = 1/2*E # ...???? energy 2x\n", - "\n", - " return trace(S*δE')\n", - " end\n", - "\n", - " R = ForwardDiff.gradient(calc_Wint, u[:])\n", - " Kt = ForwardDiff.hessian(calc_Wint, u[:])\n", - " return Kt, -R\n", + " R = ForwardDiff.gradient(W, displacement[:])\n", + " K = ForwardDiff.hessian(W, displacement[:])\n", + " return K, -R\n", "end\n", "\n", "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n", @@ -198,14 +171,12 @@ "end\n", "\n", "JuliaFEM.has_lhs(equation::CPS4) = true\n", - "JuliaFEM.has_rhs(equation::CPS4) = true\n", - "\n", - "JuliaFEM.integrate_rhs(equation, Inf)" + "JuliaFEM.has_rhs(equation::CPS4) = true" ] }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -214,19 +185,90 @@ "name": "stdout", "output_type": "stream", "text": [ + "testing primary field with point load versus code aster solution\n", + "increment 1, norm = 5.77653, Wint = 0.000, Wext = 81.525, |Wint-Wext| = 81.52516\n", + "increment 2, norm = 0.99988, Wint = 173.895, Wext = 79.037, |Wint-Wext| = 94.85823\n", + "increment 3, norm = 0.28354, Wint = 87.084, Wext = 82.175, |Wint-Wext| = 4.90953\n", + "increment 4, norm = 0.07071, Wint = 82.548, Wext = 83.100, |Wint-Wext| = 0.55141\n", + "increment 5, norm = 0.00082, Wint = 83.116, Wext = 83.109, |Wint-Wext| = 0.00644\n", + "increment 6, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", + "increment 7, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", + "increment 8, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", + "increment 9, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", + "increment 10, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", + "elapsed time: 4.257759666 seconds\n", + "1 fact verified.\n" + ] + }, + { + "data": { + "text/plain": [ + "delayed_handler (generic function with 4 methods)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "facts(\"testing primary field with point load versus code aster solution\") do\n", + " element = Quad4([1, 2, 3, 4])\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(\"youngs modulus\", [Field(0.0, 500.0)]))\n", + " push!(element, FieldSet(\"poissons ratio\", [Field(0.0, 0.3)]))\n", + " equation = CPS4(element)\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", + " fd = [3, 4, 5, 6]\n", + " f = zeros(8)\n", + " f[6] = -20.0\n", + " tic()\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", + " Wint = (-b'*u)[1]\n", + " Wext = (f'*u)[1]\n", + " W = abs(Wint-Wext)\n", + " @printf(\"increment %2d, norm = %8.5f, Wint = %8.3f, Wext = %8.3f, |Wint-Wext| = %8.5f\\n\", i, norm(du), Wint, Wext, W)\n", + " end\n", + " toc()\n", + " # verified using Code Aster.\n", + " @fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-4.15546385452579E+00)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "testing primary field with volume load versus code aster solution\n", "increment " ] }, { "data": { "text/plain": [ - "Success :: (line:-1) :: fact was true\n", - " Expression: (interpolate(element,\"displacement\",[1.0,1.0],Inf))[2] --> roughly(-4.15546385452579)\n", - " Expected: -4.15546385452579\n", - " Occurred: -4.155463854525783" + "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 41, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, @@ -234,47 +276,57 @@ "name": "stdout", "output_type": "stream", "text": [ - "1, norm = 5.77653, Wint = 0.000, Wext = 81.525, |Wint-Wext| = 81.52516\n", - "increment 2, norm = 0.99988, Wint = 173.895, Wext = 79.037, |Wint-Wext| = 94.85823\n", - "increment 3, norm = 0.28354, Wint = 87.084, Wext = 82.175, |Wint-Wext| = 4.90953\n", - "increment 4, norm = 0.07071, Wint = 82.548, Wext = 83.100, |Wint-Wext| = 0.55141\n", - "increment 5, norm = 0.00082, Wint = 83.116, Wext = 83.109, |Wint-Wext| = 0.00644\n", - "increment 6, norm = 0.00000, Wint = 83.109, Wext = 83.109, |Wint-Wext| = 0.00000\n", - "83.10927709051566\n" + "1, norm = 14.44128, Wint = -509.167, Wext = 0.000, |Wint-Wext| = 509.16667\n", + "increment 2, norm = 4.01742, Wint = 4335.728, Wext = 0.000, |Wint-Wext| = 4335.72801\n", + "increment 3, norm = 1.54645, Wint = 705.913, Wext = 0.000, |Wint-Wext| = 705.91292\n", + "increment 4, norm = 1.12361, Wint = 60.492, Wext = 0.000, |Wint-Wext| = 60.49208\n", + "increment 5, norm = 0.79119, Wint = -1.486, Wext = 0.000, |Wint-Wext| = 1.48555\n", + "increment 6, norm = 0.12733, Wint = 5.743, Wext = 0.000, |Wint-Wext| = 5.74331\n", + "increment 7, norm = 0.00725, Wint = 0.080, Wext = 0.000, |Wint-Wext| = 0.07997\n", + "increment 8, norm = 0.00001, Wint = 0.000, Wext = 0.000, |Wint-Wext| = 0.00045\n", + "increment 9, norm = 0.00000, Wint = 0.000, Wext = 0.000, |Wint-Wext| = 0.00000\n", + "increment 10, norm = 0.00000, Wint = 0.000, Wext = 0.000, |Wint-Wext| = 0.00000\n", + "elapsed time: 0.068575305 seconds\n", + "1 fact verified.\n" ] } ], "source": [ - "element = Quad4([1, 2, 3, 4])\n", - "push!(element, FieldSet(\"young\", [Field(0.0, 500.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.3)]))\n", - "equation = CPS4(element)\n", + "facts(\"testing primary field with volume load versus code aster solution\") do\n", + " element = Quad4([1, 2, 3, 4])\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(\"youngs modulus\", [Field(0.0, 500.0)]))\n", + " push!(element, FieldSet(\"poissons ratio\", [Field(0.0, 0.3)]))\n", + " push!(element, FieldSet(\"displacement volume load\",\n", + " [Field(0.0, Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]])]))\n", + " equation = CPS4(element)\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", - "fd = [3, 4, 5, 6]\n", - "f = zeros(8)\n", - "f[6] = -20.0\n", - "for i=1:6\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", - " Wint = (-b'*u)[1]\n", - " Wext = (f'*u)[1]\n", - " W = abs(Wint-Wext)\n", - " @printf(\"increment %2d, norm = %8.5f, Wint = %8.3f, Wext = %8.3f, |Wint-Wext| = %8.5f\\n\", i, norm(du), Wint, Wext, W)\n", - "end\n", - "println((f'*u)[1])\n", - "@fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-4.15546385452579E+00)\n", - "#interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2], (f'*u)[1]" + " 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", + " fd = [3, 4, 5, 6]\n", + " f = zeros(8)\n", + " f[6] = -20.0*0\n", + " tic()\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", + " Wint = (-b'*u)[1]\n", + " Wext = (f'*u)[1]\n", + " W = abs(Wint-Wext)\n", + " @printf(\"increment %2d, norm = %8.5f, Wint = %8.3f, Wext = %8.3f, |Wint-Wext| = %8.5f\\n\", i, norm(du), Wint, Wext, W)\n", + " end\n", + " toc()\n", + " # verified using Code Aster.\n", + " @fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-8.77303119819776E+00)\n", + "end" ] }, { diff --git a/src/elements.jl b/src/elements.jl index 63dcb8d..f882ba3 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -186,6 +186,11 @@ function dinterpolate(element::Element, field_name, ip::IntegrationPoint, time:: dinterpolate(element, field_name, ip.xi, time) end +""" Check does fieldset exist. """ +function Base.haskey(element::Element, what) + haskey(element.fields, symbol(what)) +end + """ Get jacobian of element evaluated at point ξ on element in reference configuration.