mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-17 17:22:10 +00:00
notebook with different optoins...
This commit is contained in:
Jukka Aho
@@ -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
|
||||
}
|
||||
+53
-5
@@ -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
|
||||
|
||||
|
||||
Reference in New Issue
Block a user