mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-23 11:02:36 +00:00
- problem can be now represented using potential energy or residual
force vector, autodiff takes care of linearization - elasticity equations are now solved using e.g. principle of minimum potential energy. syntax is quite good, see notebook. - updated how to interpolate fields, by introducing function spaces. syntax is now good. still have to figure out how to do time derivatives - etc. etc. tutorial is broken at the moment, i took of get_lhs and get_rhs because they didn't really work.
This commit is contained in:
Jukka Aho
File diff suppressed because it is too large
Load Diff
@@ -10,7 +10,7 @@
|
||||
"\n",
|
||||
"**Abstract**: Elasticity equations design notes.\n",
|
||||
"\n",
|
||||
"##Weak form\n",
|
||||
"###Weak form\n",
|
||||
"\n",
|
||||
"Given function spaces\n",
|
||||
"\\begin{align}\n",
|
||||
@@ -20,6 +20,38 @@
|
||||
"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}\n",
|
||||
"\n",
|
||||
"### Some formulas\n",
|
||||
"\\begin{align}\n",
|
||||
"J & =\\det\\left(F\\right)\\\\\n",
|
||||
"I_{c} & =\\mbox{tr}\\left(C\\right)\\\\\n",
|
||||
"\\mathbf{C} & =\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}\\\\\n",
|
||||
"\\mathbf{F} & =\\mathbf{I}+\\nabla\\mathbf{u}\\\\\n",
|
||||
"\\mathbf{E} & =\\frac{1}{2}\\left(\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}-\\mathbf{I}\\right)\n",
|
||||
"\\end{align}\n",
|
||||
"\n",
|
||||
"### Potential energy\n",
|
||||
"\n",
|
||||
"\\begin{equation}\n",
|
||||
"\\underset{u\\in\\boldsymbol{\\mathcal{U}}}{\\min}\\Pi\\left(\\mathbf{u}\\right)\n",
|
||||
"\\end{equation}\n",
|
||||
"\\begin{equation}\n",
|
||||
"\\Pi\\left(\\mathbf{u}\\right)=\\int_{\\Omega}\\psi\\left(\\mathbf{u}\\right)-\\int_{\\Omega}\\hat{\\mathbf{b}}_{0}\\cdot\\mathbf{u}-\\int_{\\Gamma_{\\sigma}}\\hat{\\mathbf{t}}_{0}\\cdot\\mathbf{u}\\,\\mathrm{d}A_{0}\n",
|
||||
"\\end{equation}\n",
|
||||
"\n",
|
||||
"### Material models\n",
|
||||
"\n",
|
||||
"https://en.wikipedia.org/wiki/Strain_energy_density_function\n",
|
||||
"\n",
|
||||
"Saint-Venant-Kirchhoff model https://en.wikipedia.org/wiki/Hyperelastic_material\n",
|
||||
"\\begin{equation}\n",
|
||||
"\\psi\\left(\\mathbf{E}\\right)=\\frac{\\lambda}{2}\\left[\\mbox{tr}\\left(\\mathbf{E}\\right)\\right]^{2}+\\mu\\mbox{tr}\\left(\\mathbf{E}^2\\right)\n",
|
||||
"\\end{equation}\n",
|
||||
"\n",
|
||||
"neo-Hookean material https://en.wikipedia.org/wiki/Neo-Hookean_solid\n",
|
||||
"\\begin{equation}\n",
|
||||
"\\psi=\\frac{\\mu}{2}\\left(I_{c}-3\\right)-\\mu\\ln\\left(J\\right)+\\frac{\\lambda}{2}\\ln\\left(J\\right)^{2}\n",
|
||||
"\\end{equation}"
|
||||
]
|
||||
},
|
||||
@@ -43,8 +75,9 @@
|
||||
],
|
||||
"source": [
|
||||
"using ForwardDiff\n",
|
||||
"using JuliaFEM: Quad4, Field, FieldSet, IntegrationPoint\n",
|
||||
"using JuliaFEM: interpolate, get_element, get_dbasisdX, dinterpolate\n",
|
||||
"using JuliaFEM: Quad4, Field, FieldSet, IntegrationPoint, Equation, LocalAssembly\n",
|
||||
"using JuliaFEM: get_element, get_basis, grad, get_integration_points\n",
|
||||
"using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n",
|
||||
"using Logging\n",
|
||||
"using FactCheck\n",
|
||||
"Logging.configure(level=DEBUG)"
|
||||
@@ -60,7 +93,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"CPS4"
|
||||
"size (generic function with 63 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
@@ -69,10 +102,12 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"abstract Elasticity <: JuliaFEM.Equation\n",
|
||||
"abstract Elasticity <: Equation\n",
|
||||
"abstract PlaneElasticity <: Elasticity\n",
|
||||
"abstract PlaneStressElasticity <: PlaneElasticity\n",
|
||||
"\n",
|
||||
"JuliaFEM.get_unknown_field_name(equation::Elasticity) = \"displacement\"\n",
|
||||
"\n",
|
||||
"\"\"\" Plane stress formulation for 4-node bilinear element. \"\"\"\n",
|
||||
"type CPS4 <: PlaneStressElasticity\n",
|
||||
" element :: Quad4\n",
|
||||
@@ -87,7 +122,9 @@
|
||||
" IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n",
|
||||
" push!(element, FieldSet(\"displacement\"))\n",
|
||||
" CPS4(element, integration_points, [])\n",
|
||||
"end"
|
||||
"end\n",
|
||||
"\n",
|
||||
"JuliaFEM.size(eq::CPS4) = 8\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -101,7 +138,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"execution_count": 25,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -109,74 +146,83 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"has_rhs (generic function with 4 methods)"
|
||||
"has_potential_energy (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"execution_count": 25,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function get_lhs_and_rhs(equation::CPS4, ip, time)\n",
|
||||
" # boilerplate code start\n",
|
||||
"\"\"\"\n",
|
||||
"Calculate internal energy of system. This can be\n",
|
||||
"used to define own material models.\n",
|
||||
"\n",
|
||||
"Parameters\n",
|
||||
"----------\n",
|
||||
"equation\n",
|
||||
" field equation we are solving\n",
|
||||
"ip\n",
|
||||
" integration point, which can be used to access fields\n",
|
||||
"time\n",
|
||||
" current time\n",
|
||||
"F\n",
|
||||
" deformation gradient\n",
|
||||
"\n",
|
||||
"Returns\n",
|
||||
"-------\n",
|
||||
"Internal energy of system.\n",
|
||||
"\"\"\"\n",
|
||||
"function calculate_internal_energy(equation::Equation, ip::IntegrationPoint, time::Number, F::Matrix)\n",
|
||||
" element = get_element(equation)\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",
|
||||
" basis = get_basis(element)\n",
|
||||
"\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",
|
||||
" # interpolate material in spatial dimension\n",
|
||||
" young = interpolate(basis, young, ip)\n",
|
||||
" poisson = interpolate(basis, poisson, ip)\n",
|
||||
" # material parameters\n",
|
||||
" young = basis(\"youngs modulus\", ip, time)\n",
|
||||
" poisson = basis(\"poissons ratio\", ip, time)\n",
|
||||
" mu = young/(2*(1+poisson))\n",
|
||||
" lambda = young*poisson/((1+poisson)*(1-2*poisson))\n",
|
||||
" lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d\n",
|
||||
"\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",
|
||||
" 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",
|
||||
" if isa(equation, PlaneStressElasticity)\n",
|
||||
" lambda = 2*lambda*mu/(lambda + 2*mu) # <- correction for 2d\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" R = ForwardDiff.gradient(W, displacement[:])\n",
|
||||
" K = ForwardDiff.hessian(W, displacement[:])\n",
|
||||
" return K, -R\n",
|
||||
" E = 1/2*(F'*F - I) # strain\n",
|
||||
" Wint = 1/2*lambda*trace(E)^2 + mu*trace(E*E')\n",
|
||||
" # alternative way to calculate this:\n",
|
||||
" # S = lambda*trace(E)*I + 2*mu*E\n",
|
||||
" # Wint = 1/2*trace(S*E')\n",
|
||||
" return Wint\n",
|
||||
"end\n",
|
||||
"\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",
|
||||
"function JuliaFEM.get_potential_energy(equation::Elasticity, ip, time; variation=nothing)\n",
|
||||
" element = get_element(equation)\n",
|
||||
" basis = get_basis(element)\n",
|
||||
" dbasis = grad(basis)\n",
|
||||
"\n",
|
||||
" u = basis(\"displacement\", ip, time, variation)\n",
|
||||
" ∇u = dbasis(\"displacement\", ip, time, variation)\n",
|
||||
" F = I + ∇u # deformation gradient\n",
|
||||
"\n",
|
||||
" # internal energy\n",
|
||||
" Wint = calculate_internal_energy(equation, ip, time, F)\n",
|
||||
"\n",
|
||||
" # external energy -- any volume load?\n",
|
||||
" Wext = 0.0\n",
|
||||
" if haskey(element, \"displacement volume load\")\n",
|
||||
" b = basis(\"displacement volume load\", ip, time)\n",
|
||||
" Wext += dot(b, u)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" return Wint - Wext\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"JuliaFEM.has_lhs(equation::CPS4) = true\n",
|
||||
"JuliaFEM.has_rhs(equation::CPS4) = true"
|
||||
"JuliaFEM.has_potential_energy(equation::CPS4) = true"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"execution_count": 26,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -185,19 +231,8 @@
|
||||
"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"
|
||||
"testing primary field with nodal load versus code aster solution\n",
|
||||
"increment "
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -206,50 +241,65 @@
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
"execution_count": 26,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
" 1, |du| = 5.77653\n",
|
||||
"increment 2, |du| = 0.99988\n",
|
||||
"increment 3, |du| = 0.28354\n",
|
||||
"increment 4, |du| = 0.07071\n",
|
||||
"increment 5, |du| = 0.00082\n",
|
||||
"increment 6, |du| = 0.00000\n",
|
||||
"increment 7, |du| = 0.00000\n",
|
||||
"increment 8, |du| = 0.00000\n",
|
||||
"increment 9, |du| = 0.00000\n",
|
||||
"increment 10, |du| = 0.00000\n",
|
||||
"elapsed time: 0.03522289 seconds\n",
|
||||
"1 fact verified.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"facts(\"testing primary field with point load versus code aster solution\") do\n",
|
||||
"facts(\"testing primary field with nodal 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 nodal load\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, -20.0], [0.0, 0.0]])]))\n",
|
||||
" equation = CPS4(element)\n",
|
||||
"\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",
|
||||
" la = initialize_local_assembly(equation)\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",
|
||||
" la = initialize_local_assembly(equation, la)\n",
|
||||
" calculate_local_assembly!(la, equation)\n",
|
||||
" du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[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",
|
||||
" @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\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",
|
||||
" @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-4.15546385452579E+00)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"execution_count": 27,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -259,7 +309,7 @@
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"testing primary field with volume load versus code aster solution\n",
|
||||
"increment "
|
||||
"increment "
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -268,7 +318,7 @@
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"execution_count": 27,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -276,17 +326,17 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"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, |du| = 14.44128\n",
|
||||
"increment 2, |du| = 4.01742\n",
|
||||
"increment 3, |du| = 1.54645\n",
|
||||
"increment 4, |du| = 1.12361\n",
|
||||
"increment 5, |du| = 0.79119\n",
|
||||
"increment 6, |du| = 0.12733\n",
|
||||
"increment 7, |du| = 0.00725\n",
|
||||
"increment 8, |du| = 0.00001\n",
|
||||
"increment 9, |du| = 0.00000\n",
|
||||
"increment 10, |du| = 0.00000\n",
|
||||
"elapsed time: 0.004752839 seconds\n",
|
||||
"1 fact verified.\n"
|
||||
]
|
||||
}
|
||||
@@ -306,26 +356,22 @@
|
||||
" u = zeros(8)\n",
|
||||
" du = zeros(8)\n",
|
||||
" fd = [3, 4, 5, 6]\n",
|
||||
" f = zeros(8)\n",
|
||||
" f[6] = -20.0*0\n",
|
||||
" la = initialize_local_assembly(equation)\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",
|
||||
" la = initialize_local_assembly(equation, la)\n",
|
||||
" calculate_local_assembly!(la, equation)\n",
|
||||
" du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[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",
|
||||
" @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\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",
|
||||
" @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-8.77303119819776E+00)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
@@ -333,104 +379,14 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Method 2, Voigt notation, analytical linearization"
|
||||
"### Method 2, Voigt notation, analytical linearization\n",
|
||||
"\n",
|
||||
"I think I don't need to mention which way is more elegant. Here's the linearization of system is done manually anyway, maybe it has better performance."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 170,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_rhs (generic function with 6 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 170,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"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",
|
||||
" 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",
|
||||
" 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",
|
||||
"\n",
|
||||
" print('.')\n",
|
||||
" \n",
|
||||
" return Kt, -R\n",
|
||||
"end\n",
|
||||
"\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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 172,
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -439,24 +395,17 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
".."
|
||||
"testing primary field with nodal load versus code aster solution\n",
|
||||
"increment "
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"8x8 Array{Float64,2}:\n",
|
||||
" 119.461 13.486 61.8587 … -26.0246 -113.693 2.80625\n",
|
||||
" 13.486 309.576 15.6443 -148.518 -3.10973 -316.14 \n",
|
||||
" 61.8587 15.6443 130.663 -53.2776 -60.0524 1.57406\n",
|
||||
" 9.73233 155.082 36.0593 -303.877 1.51267 -163.427 \n",
|
||||
" -67.6267 -26.0206 -132.47 63.0312 56.539 10.2937 \n",
|
||||
" -26.0246 -148.518 -53.2776 … 299.439 16.271 152.956 \n",
|
||||
" -113.693 -3.10973 -60.0524 16.271 117.207 -14.674 \n",
|
||||
" 2.80625 -316.14 1.57406 152.956 -14.674 326.611 "
|
||||
"delayed_handler (generic function with 4 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 172,
|
||||
"execution_count": 24,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
@@ -464,70 +413,143 @@
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
".."
|
||||
" 1, |du| = 5.77653\n",
|
||||
"increment 2, |du| = 1.05083\n",
|
||||
"increment 3, |du| = 0.39176\n",
|
||||
"increment 4, |du| = 0.21527\n",
|
||||
"increment 5, |du| = 0.16285\n",
|
||||
"increment 6, |du| = 0.13115\n",
|
||||
"increment 7, |du| = 0.10774\n",
|
||||
"increment 8, |du| = 0.08994\n",
|
||||
"increment 9, |du| = 0.07617\n",
|
||||
"increment 10, |du| = 0.06533\n",
|
||||
"elapsed time: 0.004335419 seconds\n",
|
||||
" Failure :: (line:-1) :: fact was false\n",
|
||||
" Expression: ((get_basis(element))(\"displacement\",[1.0,1.0]))[2] --> roughly(-4.15546385452579)\n",
|
||||
" Expected: -5.091745430627231 ≅ -4.15546385452579\n",
|
||||
"Out of 1 total fact:\n",
|
||||
" Failed: 1\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"JuliaFEM.integrate_lhs(equation, Inf)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 169,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2.927103942720631"
|
||||
]
|
||||
},
|
||||
"execution_count": 169,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"(1/2*u'*JuliaFEM.integrate_lhs(equation, 0.0)*u)[1]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 30,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_lhs (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 30,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"\"\"\"\n",
|
||||
"1-node point force element for plane stress problems.\n",
|
||||
"\"\"\"\n",
|
||||
"type CPS1 <: Elasticity\n",
|
||||
" element :: Point1\n",
|
||||
"\"\"\" Plane stress formulation for 4-node bilinear element, manual formulation. \"\"\"\n",
|
||||
"type CPS4M <: PlaneStressElasticity\n",
|
||||
" element :: Quad4\n",
|
||||
" integration_points :: Array{IntegrationPoint, 1}\n",
|
||||
" global_dofs :: Array{Int64, 1}\n",
|
||||
"end\n",
|
||||
"function CPS1(el::Point1)\n",
|
||||
" integration_points = []\n",
|
||||
" set_field(el, \"displacement\", zeros(2, 1))\n",
|
||||
" set_field(el, \"displacement nodal load\", zeros(2, 1))\n",
|
||||
" CPS1(el, integration_points)\n",
|
||||
"function CPS4M(element::Quad4)\n",
|
||||
" integration_points = [\n",
|
||||
" IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
|
||||
" IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
|
||||
" IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n",
|
||||
" IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n",
|
||||
" push!(element, FieldSet(\"displacement\"))\n",
|
||||
" CPS4M(element, integration_points, [])\n",
|
||||
"end\n",
|
||||
"get_rhs(eq::CPS1) = get_field(get_element(eq), \"displacement nodal load\")\n",
|
||||
"get_lhs(eq::CPS1) = None"
|
||||
"\n",
|
||||
"JuliaFEM.size(eq::CPS4M) = 8\n",
|
||||
"\n",
|
||||
"function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::CPS4M, time::Number=Inf)\n",
|
||||
" initialize_local_assembly(assembly, equation)\n",
|
||||
" element = get_element(equation)\n",
|
||||
" basis = get_basis(element)\n",
|
||||
" dbasis = grad(basis)\n",
|
||||
" detJ = det(basis)\n",
|
||||
" ndofs = size(equation)\n",
|
||||
" nnodes = round(Int, ndofs/2)\n",
|
||||
"\n",
|
||||
" B_L = zeros(3, ndofs)\n",
|
||||
" B_NL = zeros(4, ndofs)\n",
|
||||
"\n",
|
||||
" for ip in get_integration_points(equation)\n",
|
||||
"\n",
|
||||
" fill!(B_L, 0.0)\n",
|
||||
" fill!(B_NL, 0.0)\n",
|
||||
"\n",
|
||||
" u = basis(\"displacement\", ip, time)\n",
|
||||
" ∇u = dbasis(\"displacement\", ip, time)\n",
|
||||
" young = basis(\"youngs modulus\", ip, time)\n",
|
||||
" poisson = basis(\"poissons ratio\", ip, time)\n",
|
||||
"\n",
|
||||
" # kinematics\n",
|
||||
" F = I + ∇u # 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 -- calculate stress\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",
|
||||
" dNdX = dbasis(ip, time)\n",
|
||||
"\n",
|
||||
" # linear part\n",
|
||||
" for i=1:nnodes\n",
|
||||
" B_L[1, 2*(i-1)+1] = F[1,1]*dNdX[1,i]\n",
|
||||
" B_L[1, 2*(i-1)+2] = F[2,1]*dNdX[1,i]\n",
|
||||
" B_L[2, 2*(i-1)+1] = F[1,2]*dNdX[2,i]\n",
|
||||
" B_L[2, 2*(i-1)+2] = F[2,2]*dNdX[2,i]\n",
|
||||
" B_L[3, 2*(i-1)+1] = F[1,1]*dNdX[2,i] + F[1,2]*dNdX[1,i]\n",
|
||||
" B_L[3, 2*(i-1)+2] = F[2,1]*dNdX[2,i] + F[2,2]*dNdX[1,i]\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # nonlinear part\n",
|
||||
" for i=1:nnodes\n",
|
||||
" B_NL[1, 2*(i-1)+1] = dNdX[1,i]\n",
|
||||
" B_NL[2, 2*(i-1)+1] = dNdX[2,i]\n",
|
||||
" B_NL[3, 2*(i-1)+2] = dNdX[1,i]\n",
|
||||
" B_NL[4, 2*(i-1)+2] = dNdX[2,i]\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" s = ip.weight*detJ(ip)\n",
|
||||
" assembly.stiffness_matrix += s*B_L'*D*B_L + s*B_NL'*T*B_NL\n",
|
||||
" assembly.force_vector += -s*B_L'*S\n",
|
||||
"\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" if haskey(element, \"displacement nodal load\")\n",
|
||||
" assembly.force_vector += element[\"displacement nodal load\"](time)[:]\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"facts(\"testing primary field with nodal 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 nodal load\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, -20.0], [0.0, 0.0]])]))\n",
|
||||
" equation = CPS4M(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",
|
||||
" \n",
|
||||
" fd = [3, 4, 5, 6]\n",
|
||||
" la = initialize_local_assembly(equation)\n",
|
||||
" tic()\n",
|
||||
" for i=1:10\n",
|
||||
" la = initialize_local_assembly(equation, la)\n",
|
||||
" calculate_local_assembly!(la, equation)\n",
|
||||
" du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[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",
|
||||
" @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n",
|
||||
" end\n",
|
||||
" toc()\n",
|
||||
" # verified using Code Aster. Noh, toimi se eilen.\n",
|
||||
" @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-4.15546385452579E+00)\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
|
||||
@@ -10,6 +10,9 @@ using Lexicon
|
||||
using Logging
|
||||
@Logging.configure(level=DEBUG)
|
||||
|
||||
using ForwardDiff
|
||||
autodiffcache = ForwardDiffCache()
|
||||
|
||||
""" Simple linspace extension to arrays.
|
||||
|
||||
Examples
|
||||
|
||||
+116
-71
@@ -8,20 +8,18 @@ Related notebooks
|
||||
2015-08-29-developing-juliafem.ipynb
|
||||
=#
|
||||
|
||||
using JuliaFEM: interpolate
|
||||
using FactCheck
|
||||
using ForwardDiff
|
||||
|
||||
|
||||
abstract Element
|
||||
|
||||
""" Get FieldSet from element. """
|
||||
function Base.getindex(element::Element, field_name::Union{Symbol, ASCIIString})
|
||||
function Base.getindex(element::Element, field_name)
|
||||
element.fields[symbol(field_name)]
|
||||
end
|
||||
|
||||
""" Add new FieldSet to element. """
|
||||
function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name::Union{Symbol, ASCIIString})
|
||||
function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name)
|
||||
fieldset.name = symbol(fieldset_name)
|
||||
element.fields[fieldset.name] = fieldset
|
||||
end
|
||||
@@ -126,103 +124,150 @@ function test_element(element_type)
|
||||
fieldset = FieldSet("field1")
|
||||
push!(fieldset, field)
|
||||
push!(element, fieldset)
|
||||
@fact element["field1"][1] --> fld
|
||||
push!(element, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))
|
||||
|
||||
# evaluate basis functions at middle point of element
|
||||
mid = zeros(dim)
|
||||
try
|
||||
get_basis(element)(mid)
|
||||
basis = get_basis(element)
|
||||
val1 = basis(mid, 0.0)
|
||||
Logging.info("basis at $mid: $val1")
|
||||
val2 = basis("field1", mid, 0.0)
|
||||
Logging.info("field val at $mid: $val2")
|
||||
catch
|
||||
Logging.error("""
|
||||
Unable to evaluate basis, define function 'get_basis' for
|
||||
this element.""")
|
||||
end
|
||||
try
|
||||
get_dbasisdxi(element)(mid)
|
||||
basis = get_basis(element)
|
||||
dbasis = grad(basis)
|
||||
val3 = dbasis(mid, 0.0)
|
||||
Logging.info("derivative of basis at $mid: $val3")
|
||||
val4 = dbasis("field1", mid, 0.0)
|
||||
Logging.info("field val at $mid: $val4")
|
||||
catch
|
||||
Logging.error("""
|
||||
Unable to evaluate partial derivatives of basis,
|
||||
define function 'get_dbasisdxi' for this element.""")
|
||||
end
|
||||
|
||||
Logging.info("Interpolating scalar field at $mid")
|
||||
i = interpolate(element, "field1", mid, 0.0)
|
||||
Logging.info("Value: $i")
|
||||
Logging.info("Element $element_type passed tests.")
|
||||
end
|
||||
|
||||
|
||||
get_connectivity(el::Element) = el.connectivity
|
||||
|
||||
""" Get basis functions of element. """
|
||||
get_basis(el::Element) = el.basis
|
||||
get_basis(el::Element, xi::Vector) = el.basis(xi)
|
||||
Base.call(el::Element, xi::Vector) = el.basis(xi)
|
||||
|
||||
""" Get partial derivatives of basis functions of element. """
|
||||
get_dbasisdxi(el::Element) = el.basis.dbasisdxi
|
||||
get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi)
|
||||
get_dbasisdxi(el::Element, ip::IntegrationPoint) = el.basis.dbasisdxi(ip.xi)
|
||||
|
||||
""" Interpolate field on element. """
|
||||
function interpolate(element::Element, field_name, xi::Vector, time::Number)
|
||||
fieldset = element[field_name]
|
||||
field = interpolate(fieldset, time)
|
||||
basis = get_basis(element)
|
||||
interpolate(basis, field, xi)
|
||||
end
|
||||
function interpolate(element::Element, field_name, ip::IntegrationPoint, time::Number)
|
||||
interpolate(element, field_name, ip.xi, time)
|
||||
function get_connectivity(el::Element)
|
||||
el.connectivity
|
||||
end
|
||||
|
||||
""" Interpolate derivative of field on element. """
|
||||
function dinterpolate(element::Element, field_name, xi::Vector, time::Number)
|
||||
fieldset = element[field_name]
|
||||
field = interpolate(fieldset, time)
|
||||
basis = get_basis(element)
|
||||
dinterpolate(basis, field, xi)
|
||||
type MixedFunctionSpace
|
||||
element1 :: Element
|
||||
element2 :: Element
|
||||
end
|
||||
function dinterpolate(element::Element, field_name, ip::IntegrationPoint, time::Number)
|
||||
dinterpolate(element, field_name, ip.xi, time)
|
||||
|
||||
type FunctionSpace
|
||||
element :: Element
|
||||
end
|
||||
|
||||
type GradientFunctionSpace
|
||||
element :: Element
|
||||
end
|
||||
|
||||
function grad(u::FunctionSpace)
|
||||
GradientFunctionSpace(u.element)
|
||||
end
|
||||
|
||||
""" Evaluate field on element function space. """
|
||||
function call(u::FunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
|
||||
f = !isa(variation, Void) ? variation : u.element[field_name](t)
|
||||
if length(f) == 1
|
||||
return f.values
|
||||
end
|
||||
h = u.element.basis.basis(xi)
|
||||
return h*f
|
||||
end
|
||||
|
||||
""" If basis is called without a field, return basis functions evaluated at that point. """
|
||||
function call(u::FunctionSpace, xi::Vector, t::Number=Inf)
|
||||
return u.element.basis.basis(xi)'
|
||||
end
|
||||
|
||||
""" Evaluate gradient of field on element function space. """
|
||||
function call(gradu::GradientFunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
|
||||
f = !isa(variation, Void) ? variation : gradu.element[field_name](t)
|
||||
X = gradu.element["geometry"](t)
|
||||
b = gradu.element.basis.dbasisdxi(xi)
|
||||
return b*f*inv(b*X)
|
||||
end
|
||||
|
||||
""" If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """
|
||||
function call(gradu::GradientFunctionSpace, xi::Vector, t::Number=Inf)
|
||||
X = gradu.element["geometry"](t)
|
||||
b = gradu.element.basis.dbasisdxi(xi)
|
||||
return (b*inv(b*X))'
|
||||
end
|
||||
|
||||
# on-line functions to get api more easy to use, ip -> xi.ip
|
||||
call(u::FunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t)
|
||||
call(u::FunctionSpace, ip::IntegrationPoint) = call(u, ip.xi)
|
||||
call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t)
|
||||
call(u::GradientFunctionSpace, ip::IntegrationPoint) = call(u, ip.xi)
|
||||
|
||||
""" Return field from function space. """
|
||||
function get_field(u::FunctionSpace, field_name, time=Inf)
|
||||
return u.element[field_name](time)
|
||||
end
|
||||
|
||||
""" Return field from function space. """
|
||||
function get_field(u::FunctionSpace, field_name, time=Inf, variation=nothing)
|
||||
return !isa(variation, Void) ? variation : u.element[field_name](time)
|
||||
end
|
||||
|
||||
""" Return fieldset from function space. """
|
||||
function get_fieldset(u::FunctionSpace, field_name)
|
||||
return u.element[field_name]
|
||||
end
|
||||
|
||||
# i think these will be the most called functions.
|
||||
call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
|
||||
call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
|
||||
|
||||
function jacobian(u::FunctionSpace, xi, t)
|
||||
u.element.basis.dbasisdxi(xi)*u.element["geometry"](t)
|
||||
end
|
||||
|
||||
function jacobian(u::FunctionSpace, ip::IntegrationPoint, t::Number)
|
||||
jacobian(u, ip.xi, t)
|
||||
end
|
||||
|
||||
function jacobian(u::FunctionSpace, xi)
|
||||
jacobian(u, xi, Inf)
|
||||
end
|
||||
|
||||
function LinAlg.det(u::FunctionSpace)
|
||||
function detJ(args...)
|
||||
J = jacobian(u, args...)
|
||||
m, n = size(J)
|
||||
return m == n ? det(J) : norm(J)
|
||||
end
|
||||
return detJ
|
||||
end
|
||||
|
||||
function get_basis(element::Element)
|
||||
return FunctionSpace(element)
|
||||
end
|
||||
|
||||
Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
|
||||
Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
|
||||
Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
|
||||
Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
|
||||
|
||||
|
||||
""" 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.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
element :: Element
|
||||
xi :: Vector
|
||||
spatial coordinate
|
||||
time :: Float64
|
||||
temporal coordinate
|
||||
geometry_field :: optional
|
||||
|
||||
Returns
|
||||
-------
|
||||
Vector or Matrix
|
||||
depending on element dimension
|
||||
|
||||
"""
|
||||
function get_jacobian(element::Element, xi, time, geometry_field="geometry")
|
||||
dinterpolate(element, geometry_field, xi, time)
|
||||
end
|
||||
|
||||
""" Evaluate partial derivatives of basis, dbasis/dX, at some time t"""
|
||||
function get_dbasisdX(el::Element, xi, t)
|
||||
dbasisdxi = get_dbasisdxi(el, xi)
|
||||
J = get_jacobian(el, xi, t)
|
||||
dbasisdxi*inv(J)
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
# FIXME: These two needs integration -- maybe not in elements.jl ..?
|
||||
"""
|
||||
|
||||
+151
-45
@@ -3,68 +3,174 @@
|
||||
|
||||
abstract Equation
|
||||
|
||||
abstract Assembly
|
||||
|
||||
""" Local element assembly. """
|
||||
type LocalAssembly <: Assembly
|
||||
ndofs :: Int
|
||||
mass_matrix :: Matrix
|
||||
stiffness_matrix :: Matrix
|
||||
force_vector :: Matrix
|
||||
potential_energy# :: Union{Array, Float64}
|
||||
residual_vector :: Vector
|
||||
end
|
||||
|
||||
function LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector::Matrix)
|
||||
LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector[:])
|
||||
end
|
||||
|
||||
""" Initialize workspace for local assembly. """
|
||||
function LocalAssembly(equation::Equation)
|
||||
ndofs = size(equation)
|
||||
mass_matrix = zeros(ndofs, ndofs)
|
||||
stiffness_matrix = zeros(ndofs, ndofs)
|
||||
force_vector = zeros(ndofs, 1)
|
||||
potential_energy = 0.0
|
||||
residual_vector = zeros(ndofs)
|
||||
return LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector,
|
||||
potential_energy, residual_vector)
|
||||
end
|
||||
|
||||
function initialize_local_assembly(equation::Equation)
|
||||
LocalAssembly(equation)
|
||||
end
|
||||
|
||||
function initialize_local_assembly(equation::Equation, assembly::LocalAssembly)
|
||||
if size(equation) != assembly.ndofs
|
||||
# if problem size changes, automatically initialize new work space
|
||||
return initialize_local_assembly(equation)
|
||||
end
|
||||
# otherwise, empty workspace ready for next iteration
|
||||
fill!(assembly.mass_matrix, 0.0)
|
||||
fill!(assembly.stiffness_matrix, 0.0)
|
||||
fill!(assembly.force_vector, 0.0)
|
||||
assembly.potential_energy = 0.0
|
||||
fill!(assembly.residual_vector, 0.0)
|
||||
return assembly
|
||||
end
|
||||
function initialize_local_assembly(assembly::LocalAssembly, equation::Equation)
|
||||
initialize_local_assembly(equation, assembly)
|
||||
end
|
||||
|
||||
function get_unknown_field_name(equation::Equation)
|
||||
eqtype = typeof(equation)
|
||||
error("define get_unknown_field_name for this equation type $eqtype")
|
||||
end
|
||||
|
||||
|
||||
has_lhs(eq::Equation) = false
|
||||
get_lhs(eq::Equation, xi) = nothing
|
||||
has_rhs(eq::Equation) = false
|
||||
get_rhs(eq::Equation, xi) = nothing
|
||||
get_element(eq::Equation) = eq.element
|
||||
get_integration_points(eq::Equation) = eq.integration_points
|
||||
|
||||
# couple convenient functions -- could make weak form definition easier
|
||||
get_connectivity(eq::Equation) = get_connectivity(get_element(eq))
|
||||
get_basis(eq::Equation, ip::IntegrationPoint) = get_basis(get_element(eq), ip.xi)
|
||||
get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq), ip.xi)
|
||||
interpolate(eq::Equation, field::Union{ASCIIString, Symbol}, ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
|
||||
integrate_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
|
||||
integrate_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
|
||||
get_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
|
||||
get_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
|
||||
has_mass_matrix(equation::Equation) = false
|
||||
get_mass_matrix(equation::Equation, ip, time) = nothing
|
||||
has_stiffness_matrix(equation::Equation) = false
|
||||
get_stiffness_matrix(equation::Equation, ip, time) = nothing
|
||||
has_force_vector(equation::Equation) = false
|
||||
get_force_vector(equation::Equation, ip, time) = nothing
|
||||
has_residual_vector(equation::Equation) = false
|
||||
get_residual_vector(equation::Equation, ip, time) = nothing
|
||||
has_potential_energy(equation::Equation) = false
|
||||
get_potential_energy(equation::Equation, ip, time) = nothing
|
||||
get_element(equation::Equation) = equation.element
|
||||
get_number_of_dofs(equation::Equation) = nothing
|
||||
get_integration_points(equation::Equation) = equation.integration_points
|
||||
|
||||
|
||||
"""
|
||||
Return determinant of Jacobian for numerical integration.
|
||||
"""
|
||||
function get_detJ(eq::Equation, ip::IntegrationPoint, t::Float64)
|
||||
el = get_element(eq)
|
||||
get_detJ(el, ip, t)
|
||||
end
|
||||
function get_detJ(el::Element, ip::IntegrationPoint, t::Float64)
|
||||
get_detJ(el, ip.xi, t)
|
||||
end
|
||||
function get_detJ(el::Element, xi::Vector, t::Float64)
|
||||
J = get_jacobian(el, xi, t)
|
||||
s = size(J)
|
||||
return s[1] == s[2] ? det(J) : norm(J)
|
||||
end
|
||||
""" Return a local assembly for element. """
|
||||
function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, time::Number=Inf)
|
||||
|
||||
"""
|
||||
Integrate f over element
|
||||
initialize_local_assembly(assembly, equation) # zero all
|
||||
|
||||
Parameters
|
||||
----------
|
||||
eq::Equation
|
||||
element = get_element(equation)
|
||||
basis = get_basis(element)
|
||||
detJ = det(basis)
|
||||
field_name = get_unknown_field_name(equation)
|
||||
|
||||
f::Function
|
||||
Function to integrate
|
||||
"""
|
||||
function integrate(eq::Equation, f::Function, t::Float64)
|
||||
target = []
|
||||
for ip in get_integration_points(eq)
|
||||
push!(target, ip.weight*f(eq, ip, t)*get_detJ(eq, ip, t))
|
||||
# 1. if equations are defined we just integrate them
|
||||
if has_mass_matrix(equation) || has_stiffness_matrix(equation) || has_force_vector(equation)
|
||||
for ip in get_integration_points(equation)
|
||||
s = ip.weight*detJ(ip)
|
||||
if has_mass_matrix(equation)
|
||||
assembly.mass_matrix += s*get_mass_matrix(equation, ip, time)
|
||||
end
|
||||
if has_stiffness_matrix(equation)
|
||||
assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time)
|
||||
end
|
||||
if has_force_vector(equation)
|
||||
assembly.force_vector += s*get_force_vector(equation, ip, time)[:]
|
||||
end
|
||||
# external loads -- if any nodal loads is defined add to force vector
|
||||
if haskey(element, "$field_name nodal load")
|
||||
assembly.force_vector += element["$field_name nodal load"](time)[:]
|
||||
end
|
||||
end
|
||||
end
|
||||
return sum(target)
|
||||
|
||||
# 2. variational / energy form - user has defined some potential energy / variational form
|
||||
if has_potential_energy(equation)
|
||||
field_name = get_unknown_field_name(equation)
|
||||
element = get_element(equation)
|
||||
field = element[field_name](time)
|
||||
function potential_energy(data::Vector)
|
||||
# calculate potential energy for some setting. this is needed by forwarddiff
|
||||
assembly.potential_energy = 0.0
|
||||
df = similar(field, data)
|
||||
# integrate potential energy
|
||||
for ip in get_integration_points(equation)
|
||||
dw = get_potential_energy(equation, ip, time; variation=df)
|
||||
assembly.potential_energy += ip.weight * dw * detJ(ip)
|
||||
end
|
||||
# external energy -- if any nodal loads is defined, decrease from potential energy
|
||||
if haskey(element, "$field_name nodal load")
|
||||
P = element["$field_name nodal load"](time)
|
||||
assembly.potential_energy -= dot(P[:], df[:])
|
||||
end
|
||||
if isa(assembly.potential_energy, Array)
|
||||
return assembly.potential_energy[1]
|
||||
end
|
||||
return assembly.potential_energy
|
||||
end
|
||||
hessian, allresults = ForwardDiff.hessian(potential_energy, field[:],
|
||||
AllResults, cache=autodiffcache)
|
||||
assembly.stiffness_matrix += hessian
|
||||
assembly.force_vector -= ForwardDiff.gradient(allresults) # <--- minus explained in tutorial
|
||||
assembly.potential_energy = ForwardDiff.value(allresults)
|
||||
end
|
||||
|
||||
# 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v
|
||||
if has_residual_vector(equation)
|
||||
field_name = get_unknown_field_name(equation)
|
||||
element = get_element(equation)
|
||||
field = element[field_name](time)
|
||||
function residual_vector(data::Vector)
|
||||
fill!(assembly.residual_vector, 0.0)
|
||||
df = similar(field, data)
|
||||
# integrate W
|
||||
for ip in get_integration_points(equation)
|
||||
dr = get_residual_vector(equation, ip, time; variation=df)
|
||||
assembly.residual_vector += ip.weight*dr*detJ(ip)
|
||||
end
|
||||
# external loads -- if any nodal loads is defined, remove from residual
|
||||
if haskey(element, "$field_name nodal load")
|
||||
assembly.residual_vector -= element["$field_name nodal load"](time)[:]
|
||||
end
|
||||
return assembly.residual_vector
|
||||
end
|
||||
jacobian, allresults = ForwardDiff.jacobian(residual_vector, field[:],
|
||||
AllResults, cache=autodiffcache)
|
||||
assembly.stiffness_matrix += jacobian
|
||||
assembly.force_vector -= ForwardDiff.value(allresults) # <-- minus explained in tutorial
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function calculate_local_assembly!(equation::Equation, assembly::LocalAssembly, time::Number=Inf)
|
||||
calculate_local_assembly!(assembly, equation)
|
||||
end
|
||||
|
||||
|
||||
""" Get global degrees of freedom for this element. """
|
||||
function get_global_dofs(eq::Equation)
|
||||
eq.global_dofs
|
||||
end
|
||||
|
||||
""" Set global degrees of freedom for this element. """
|
||||
function set_global_dofs!(eq::Equation, dofs)
|
||||
eq.global_dofs = dofs
|
||||
end
|
||||
|
||||
+2
-3
@@ -59,8 +59,7 @@ 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
|
||||
function dinterpolate(basis::Basis, u::Field, xi::Array{Float64, 1})
|
||||
basis.dbasisdxi(xi)*u
|
||||
end
|
||||
|
||||
|
||||
+14
-7
@@ -70,12 +70,17 @@ 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)
|
||||
fdim = round(Int, length(data)/length(field)) # dimension of field variable
|
||||
if fdim == 1
|
||||
new_field = Field(field.time, data)
|
||||
return new_field
|
||||
end
|
||||
new_field = Field(field.time, similar(field.values))
|
||||
data = reshape(data, round(Int, length(data)/length(field)), length(field))
|
||||
data = reshape(data, fdim, length(field))
|
||||
for i=1:length(new_field)
|
||||
new_field.values[i] = data[:,i]
|
||||
end
|
||||
new_field
|
||||
return new_field
|
||||
end
|
||||
|
||||
|
||||
@@ -117,7 +122,6 @@ function Base.endof(fieldset::FieldSet)
|
||||
end
|
||||
|
||||
|
||||
|
||||
""" Basis function. """
|
||||
type Basis
|
||||
basis :: Function
|
||||
@@ -128,9 +132,9 @@ function Basis(basis)
|
||||
Basis(basis, ForwardDiff.jacobian(basis))
|
||||
end
|
||||
""" Get partial derivative of basis function. """
|
||||
diff(h::Basis) = h.dbasisdxi
|
||||
derivative(h::Basis) = h.dbasisdxi
|
||||
|
||||
function grad(basis::Basis)
|
||||
(ip) -> basis.dbasisdxi(ip.xi)
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
@@ -158,7 +162,10 @@ end
|
||||
|
||||
# convenient functions -- maybe this is not correct place for them
|
||||
""" Evaluate basis function in point ξ. """
|
||||
call(b::Basis, xi) = b.basis(xi)
|
||||
call(b::Basis, xi::Vector) = b.basis(xi)
|
||||
call(b::Basis, ip::IntegrationPoint) = b.basis(ip.xi)
|
||||
Base.(:*)(basis::Basis, fs::FieldSet) = (xi, t) -> basis(xi)*fs(t)
|
||||
|
||||
#""" Interpolate field (h*f)(ξ) """
|
||||
#Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
|
||||
#""" Interpolate from set of fields with basis b, i.e. f(t) = b(t)*[f1, f2] """
|
||||
|
||||
+33
-1
@@ -2,7 +2,7 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using FactCheck
|
||||
using JuliaFEM: Element, Basis, FieldSet
|
||||
using JuliaFEM: Element, Basis, Field, FieldSet, FunctionSpace
|
||||
|
||||
""" Prototype element
|
||||
|
||||
@@ -51,3 +51,35 @@ facts("test adding fieldsets and fields to element") do
|
||||
@fact fields[2] --> field2
|
||||
end
|
||||
|
||||
facts("interpolation of fields in some function space") do
|
||||
|
||||
element = MockElement([1, 2, 3, 4])
|
||||
fieldset1 = FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])
|
||||
fieldset2 = FieldSet("constant scalar field", [Field(0.0, 1.0)])
|
||||
fieldset3 = FieldSet("scalar field", [Field(0.0, [1.0, 2.0, 3.0, 4.0])])
|
||||
fieldset4 = FieldSet("vector field 1", [Field(0.0, Vector[[1.0], [2.0], [3.0], [4.0]])])
|
||||
fieldset5 = FieldSet("vector field 2", [Field(0.0, Vector[[1.0, 5.0], [2.0, 6.0], [3.0, 7.0], [4.0, 8.0]])])
|
||||
fieldset6 = FieldSet("vector field 3", [Field(0.0, Vector[[1.0, 5.0, 9.0], [2.0, 6.0, 10.0], [3.0, 7.0, 11.0], [4.0, 8.0, 12.0]])])
|
||||
fieldset7 = FieldSet("tensor field 1", [Field(0.0, Matrix[[1.0 5.0; 9.0 13.0], [2.0 6.0; 10.0 14.0], [3.0 7.0; 11.0 15.0], [4.0 8.0; 12.0 16.0]])])
|
||||
|
||||
push!(element, fieldset1)
|
||||
push!(element, fieldset2)
|
||||
push!(element, fieldset3)
|
||||
push!(element, fieldset4)
|
||||
push!(element, fieldset5)
|
||||
push!(element, fieldset6)
|
||||
push!(element, fieldset7)
|
||||
|
||||
xi = [0.0, 0.0]
|
||||
t = 0.0
|
||||
u = FunctionSpace(element)
|
||||
v = FunctionSpace(element)
|
||||
|
||||
@fact v("constant scalar field", xi, t) --> 1.0
|
||||
@fact v("scalar field", xi, t) --> 1/4*(1+2+3+4)
|
||||
@fact v("vector field 1", xi, t) --> [1/4*(1+2+3+4)]
|
||||
@fact v("vector field 2", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8]
|
||||
@fact v("vector field 3", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8, 9+10+11+12]
|
||||
@fact v("tensor field 1", xi, t) --> 1/4*[1+2+3+4 5+6+7+8; 9+10+11+12 13+14+15+16]
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,50 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
|
||||
using FactCheck
|
||||
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
|
||||
geometry_field = Field(0.0, Vector[]) # Create empty field at time t=0.0
|
||||
push!(geometry_field, [ 0.0, 0.0]) # push some values for field
|
||||
push!(geometry_field, [ 1.0, 0.0])
|
||||
push!(geometry_field, [ 1.0, 1.0])
|
||||
push!(geometry_field, [ 0.0, 1.0])
|
||||
geometry_fieldset = FieldSet("geometry") # create fieldset "geometry"
|
||||
push!(geometry_fieldset, geometry_field) # add field to fieldset
|
||||
push!(element, geometry_fieldset) # add fieldset to element
|
||||
|
||||
temperature_fieldset = FieldSet("temperature")
|
||||
push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))
|
||||
push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))
|
||||
push!(element, temperature_fieldset)
|
||||
|
||||
displacement_fieldset = FieldSet("displacement")
|
||||
push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))
|
||||
push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))
|
||||
push!(element, displacement_fieldset)
|
||||
|
||||
facts("basic continuum interpolations") do
|
||||
# from my old home works
|
||||
basis = get_basis(element)
|
||||
dbasis = grad(basis)
|
||||
@fact basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0) --> [9/16, 1/2]
|
||||
gradu = dbasis("displacement", [0.0, 0.0], 1.0)
|
||||
epsilon = 1/2*(gradu + gradu')
|
||||
rotation = 1/2*(gradu - gradu')
|
||||
X = basis("geometry", [0.0, 0.0], 1.0)
|
||||
k = 0.25
|
||||
epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
|
||||
rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
|
||||
@fact epsilon --> roughly(epsilon_wanted)
|
||||
@fact rotation --> roughly(rotation_wanted)
|
||||
F = I + gradu
|
||||
@fact F --> [X[2]*k+1 X[1]*k; 0 1]
|
||||
C = F'*F
|
||||
@fact C --> [(X[2]*k+1)^2 (X[2]*k+1)*X[1]*k; (X[2]*k+1)*X[1]*k X[1]^2*k^2+1]
|
||||
E = 1/2*(F'*F - I)
|
||||
@fact E --> [1/2*(X[2]*k + 1)^2-1/2 1/2*(X[2]*k+1)*X[1]*k; 1/2*(X[2]*k + 1)*X[1]*k 1/2*X[1]^2*k^2]
|
||||
U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I)
|
||||
#@fact U --> roughly([1.24235 0.13804; 0.13804 1.02149])
|
||||
end
|
||||
Reference in New Issue
Block a user