diff --git a/notebooks/2015-08-15-interpolation-integration-and-linearization-strategies.ipynb b/notebooks/2015-08-15-interpolation-integration-and-linearization-strategies.ipynb new file mode 100644 index 0000000..bf1c07f --- /dev/null +++ b/notebooks/2015-08-15-interpolation-integration-and-linearization-strategies.ipynb @@ -0,0 +1,548 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Interpolation and integration algorithms\n", + "\n", + "Author(s): Jukka Aho\n", + "\n", + "**Abstract**: Some strategies to implement automatic differentiation. The number of different choises is caused by a fact that the linearization of function can be done before integration or vice versa, and functions can return values or do in-place modifications. There is probably performance differences between different strategies, but all of them should work." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM\n", + "using ForwardDiff" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\"Old good\" elasticity force equilibrium equation $R = T - F$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "calc_residual_vector_integrand (generic function with 1 method)" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n", + " # Calculate dN/dX and interpolate material parameters\n", + " dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n", + " u = el.attributes[\"displacement\"]\n", + " lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n", + " mu = JuliaFEM.interpolate(el, \"mu\", xi)\n", + "\n", + " # Calculate residual force vector R = T - F\n", + " gradu = u*dbasisdX\n", + " F = I + gradu\n", + " E = 1/2*(gradu' + gradu + gradu'*gradu)\n", + " S = lambda*trace(E)*I + 2*mu*E\n", + " P = F*S\n", + " T = P*dbasisdX'\n", + " return T\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Test case, already well known 2d elasticity in [0,10] x [0,1] grid." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "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", + "ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]'\n", + "iweights = [1, 1, 1, 1]\n", + "attributes = Dict()\n", + "e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\n", + "\n", + "E = 90\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", + "\n", + "e.attributes[\"coordinates\"] = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n", + "e.attributes[\"lambda\"] = la\n", + "e.attributes[\"mu\"] = mu\n", + "e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.5 0.0; 0.0 0.0]'\n", + "e.attributes[\"displacement nodal force\"] = zeros(2, 4)\n", + "e.attributes[\"displacement tangent stiffness\"] = zeros(8, 8);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integration\n", + "\n", + "1. take element and function and return value\n", + "2. take function and return function which can be integrated by passing element as a function\n", + "3. do in-place integration, save values to target" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "integrate! (generic function with 1 method)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function integrate(f::Function, el::JuliaFEM.Element)\n", + " target = []\n", + " for m = 1:length(el.iweights)\n", + " w = el.iweights[m]\n", + " xi = el.ipoints[:, m]\n", + " J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n", + " push!(target, w*f(el, xi)*det(J))\n", + " end\n", + " return sum(target)\n", + "end\n", + "\n", + "function integrate(f::Function)\n", + " function integrate(el::JuliaFEM.Element)\n", + " target = []\n", + " for m = 1:length(el.iweights)\n", + " w = el.iweights[m]\n", + " xi = el.ipoints[:, m]\n", + " J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n", + " push!(target, w*f(el, xi)*det(J))\n", + " end\n", + " return sum(target)\n", + " end\n", + " return integrate\n", + "end\n", + "\n", + "function integrate!(f::Function, el::JuliaFEM.Element, target)\n", + " # set target to zero\n", + " el.attributes[target][:] = 0.0\n", + " for m = 1:length(el.iweights)\n", + " w = el.iweights[m]\n", + " xi = el.ipoints[:, m]\n", + " J = JuliaFEM.interpolate(el, \"coordinates\", xi; derivative=true)\n", + " el.attributes[target][:,:] += w*f(el, xi)*det(J)\n", + " end\n", + "end\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Float64,2}:\n", + " -38.2303 -72.8697 79.4912 31.6088\n", + " -17.625 -28.475 37.7 8.4 " + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "integrate(calc_residual_vector_integrand, e)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Float64,2}:\n", + " -38.2303 -72.8697 79.4912 31.6088\n", + " -17.625 -28.475 37.7 8.4 " + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "calc_residual_vector = integrate(calc_residual_vector_integrand)\n", + "calc_residual_vector(e)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "2x4 Array{Float64,2}:\n", + " -38.2303 -72.8697 79.4912 31.6088\n", + " -17.625 -28.475 37.7 8.4 " + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "integrate!(calc_residual_vector_integrand, e, \"displacement nodal force\")\n", + "e.attributes[\"displacement nodal force\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Linearization\n", + "\n", + "1. take function, element and field, and return partial derivative\n", + "2. take function and field, return function which takes element as argument\n", + "3. do in-place linearization to target, requires function which takes element as argument\n", + "\n", + "In general linearization can be done before integration and vice versa, i.e.\n", + "\n", + " integrate(linearize(f, \"u\"))(e) <-> linearize(integrate(f), \"u\")(e)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "linearize! (generic function with 1 method)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function linearize(f::Function, el::JuliaFEM.Element, field::ASCIIString)\n", + " dim, nnodes = size(el.attributes[field])\n", + " function helper!(x, y)\n", + " orig = copy(el.attributes[field])\n", + " el.attributes[field] = reshape(x, dim, nnodes)\n", + " y[:] = f(el)\n", + " el.attributes[field] = copy(orig)\n", + " end\n", + " jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n", + " return jac(el.attributes[field][:])\n", + "end\n", + "\n", + "function linearize(f::Function, field::ASCIIString)\n", + " function jacobian(el::JuliaFEM.Element, args...)\n", + " dim, nnodes = size(el.attributes[field])\n", + " function helper!(x, y)\n", + " orig = copy(el.attributes[field])\n", + " el.attributes[field] = reshape(x, dim, nnodes)\n", + " y[:] = f(el, args...)\n", + " el.attributes[field] = copy(orig)\n", + " end\n", + " jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n", + " return jac(el.attributes[field][:])\n", + " end\n", + " return jacobian\n", + "end\n", + "\n", + "function linearize!(f::Function, el::JuliaFEM.Element, field::ASCIIString, target::ASCIIString)\n", + " el.attributes[target][:] = 0.0\n", + " dim, nnodes = size(el.attributes[field])\n", + " function helper!(x, y)\n", + " orig = copy(el.attributes[field])\n", + " el.attributes[field] = reshape(x, dim, nnodes)\n", + " y[:] = f(el)\n", + " el.attributes[field] = copy(orig)\n", + " end\n", + " jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n", + " jac!(el.attributes[field][:], el.attributes[target])\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n", + " 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n", + " 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n", + " 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n", + " -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n", + " -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n", + " -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n", + " -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "integrate(linearize(calc_residual_vector_integrand, \"displacement\"))(e)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n", + " 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n", + " 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n", + " 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n", + " -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n", + " -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n", + " -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n", + " -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n", + " 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n", + " 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n", + " 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n", + " -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n", + " -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n", + " -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n", + " -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linearize(integrate(calc_residual_vector_integrand), e, \"displacement\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n", + " 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n", + " 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n", + " 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n", + " -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n", + " -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n", + " -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n", + " -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "linearize!(integrate(calc_residual_vector_integrand), e, \"displacement\", \"displacement tangent stiffness\")\n", + "e.attributes[\"displacement tangent stiffness\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "8x8 Array{Float64,2}:\n", + " 149.721 55.55 84.679 36.65 … -55.55 -136.278 -36.65 \n", + " 55.55 329.69 42.75 167.935 -172.941 -42.8 -324.684\n", + " 84.679 42.75 185.321 105.05 -123.05 -73.522 -24.75 \n", + " 36.65 167.935 105.05 340.54 -344.759 -24.8 -163.716\n", + " -98.122 -55.5 -196.478 -116.9 135.8 76.233 36.6 \n", + " -55.55 -172.941 -123.05 -344.759 … 352.922 42.8 164.778\n", + " -136.278 -42.8 -73.522 -24.8 42.8 133.567 24.8 \n", + " -36.65 -324.684 -24.75 -163.716 164.778 24.8 323.622" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "integrate!(linearize(calc_residual_vector_integrand, \"displacement\"), e, \"displacement tangent stiffness\")\n", + "e.attributes[\"displacement tangent stiffness\"]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Validations" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged in 6 iterations.\n" + ] + }, + { + "data": { + "text/plain": [ + "2x4 Array{Float64,2}:\n", + " 0.0 -0.399145 -0.0722858 0.0\n", + " 0.0 -2.17799 -2.22224 0.0" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "free_dofs = [3, 4, 5, 6]\n", + "u = zeros(2, 4)\n", + "du = zeros(2, 4)\n", + "F = [0 0; 0 0; 0 -2; 0 0]'\n", + "for i=1:10\n", + " e.attributes[\"displacement\"] = u\n", + " K = linearize(integrate(calc_residual_vector_integrand), \"displacement\")(e)\n", + " R = integrate(calc_residual_vector_integrand)(e)\n", + " du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n", + " u += du\n", + " if norm(du) < 1.0e-9\n", + " println(\"Converged in $i iterations.\")\n", + " break\n", + " end\n", + "end\n", + "u # -2.222244754401764" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.0-dev", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "name": "julia", + "version": "0.4.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/math.jl b/src/math.jl index df1e010..5f0864b 100644 --- a/src/math.jl +++ b/src/math.jl @@ -1,6 +1,10 @@ # 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, ... +""" + using JuliaFEM using ForwardDiff @@ -59,16 +63,23 @@ function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; deriva end +""" + +""" function get_basis(el::Element, xi) return el.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.dbasis(xi)*inv(J) return dbasisdX end + """ Linearize function f w.r.t some given field, i.e. calculate dR/du @@ -78,14 +89,35 @@ f::Function (possibly) nonlinear function to linearize field::ASCIIString field variable + +Returns +------- +Array{Float64, 2} + jacobian / "tangent stiffness matrix" + +""" +function linearize(f::Function, el::JuliaFEM.Element, field::ASCIIString) + dim, nnodes = size(el.attributes[field]) + function helper!(x, y) + orig = copy(el.attributes[field]) + el.attributes[field] = reshape(x, dim, nnodes) + y[:] = f(el) + el.attributes[field] = copy(orig) + end + jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) + return jac(el.attributes[field][:]) +end + +""" +This version returns another function which can be then evaluated against field """ function linearize(f::Function, field::ASCIIString) - function jacobian(el::Element, xi) + function jacobian(el::JuliaFEM.Element, args...) dim, nnodes = size(el.attributes[field]) function helper!(x, y) orig = copy(el.attributes[field]) el.attributes[field] = reshape(x, dim, nnodes) - y[:] = f(el, xi) + y[:] = f(el, args...) el.attributes[field] = copy(orig) end jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) @@ -94,6 +126,23 @@ function linearize(f::Function, field::ASCIIString) return jacobian end +""" +In-place version, no additional garbage collection. +""" +function linearize!(f::Function, el::JuliaFEM.Element, field::ASCIIString, target::ASCIIString) + el.attributes[target][:] = 0.0 + dim, nnodes = size(el.attributes[field]) + function helper!(x, y) + orig = copy(el.attributes[field]) + el.attributes[field] = reshape(x, dim, nnodes) + y[:] = f(el) + el.attributes[field] = copy(orig) + end + jac! = ForwardDiff.forwarddiff_jacobian!(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes) + jac!(el.attributes[field][:], el.attributes[target]) +end + + """ Integrate f over element using Gaussian quadrature rules. @@ -138,12 +187,11 @@ This version saves results inplace to target, garbage collection free """ function integrate!(f::Function, el::JuliaFEM.Element, target) # set target to zero - target[:] = 0.0 + el.attributes[target][:] = 0.0 for m = 1:length(el.iweights) w = el.iweights[m] xi = el.ipoints[:, m] J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true) - target[:,:] += w*f(el, xi)*det(J) + el.attributes[target][:,:] += w*f(el, xi)*det(J) end end -