Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-10 13:17:42 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
b985ebf50cf9dd62b800cfaff3dfda9cd31dca1f
JuliaFEM.jl/docs/tutorials/2015-08-29-developing-juliafem.ipynb
T
Jukka Aho b985ebf50c data structures, new testing concept
2015-11-01 18:44:50 +02:00

1832 lines
66 KiB
Plaintext
Raw Blame History

This file contains invisible Unicode characters
This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Developing JuliaFEM\n",
"\n",
"Author(s): Jukka Aho\n",
"\n",
"**Abstract**: Developer notes. In this notebook we give the general guidelines how to implement basic fundamentals to JuliaFEM. The notebook is rapidly updated and should always represent the newest \"style\" how to develop things to JuliaFEM. At this phase of project interfaces are changing rapidly. For this reason we aim to develop couple of simple functions which aim to test that elements, equations, solvers etc. work as expected, in the sense that their interface hasn't any design problems. (Yes, we do believe in testing and follow also other practices proven to be good, like continuous integration). The aim of this document is not to be the guide to \"other contributors\" but is more like an guide line for all contributors, including ourselves, how things should be done.\n",
"\n",
"For demonstrational purposes we use the Poisson equation in our examples to demonstrate the main concepts of the package. Poisson equation has some analogies to solid mechanics. Consider for example $(EAu')' + q = 0$, which is nothing more than a Poisson equation defined in 1d.\n",
"\n",
"Because this is the second important document in this source repository, any questions or suggestions araising from this notebook are very welcome and desirable and should be adressed to the our issue log: https://github.com/JuliaFEM/JuliaFEM.jl/issues. Easiest way to make world better is to register GitHub (takes 1 min) and drop an issue. And the most important document? It's the \"how to contribute\" guide, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/CONTRIBUTING.rst.\n",
"\n",
"In this notebook briefly we have defined the following \"objects\":\n",
"\n",
"**Element**: object which holds basis functions and fields. One is able to interpolate fields using element. Introduction of new basis functions needs new elements. Typical elements: Lagrange elements, hierarchical elements, etc.\n",
"\n",
"**Equation**: object which defines some field equation, like Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc., which needs to be solve. This is the physics we are trying to solve by approximating the equations in element area by some interpolation function provided by element.\n",
"\n",
"**Problem**: object which maps equations to elements. For example, HeatProblem which maps Poisson equation to Lagrange elements or ElasticityProblem which maps elasticity equation to Lagrange elements.\n",
"\n",
"**Solver**: object which takes one or more problems, solves them using some strategy (iterative methods, multigrid, direct methods, ...) and updates corresponding fields to elements. The is the section where any numerical crunching happens and is targeted to high-performance computing. \n",
"\n",
"These concepts are quite well separated so that development process can focus just one of thing of interest. In this notebook we will give instructions to all of the sections describing the whole development process from element level to global assembly and solution.\n",
"\n",
"We have\n",
"- [x] **`test_element`**\n",
"- [ ] **`test_equation`**\n",
"- [ ] **`test_problem`** \n",
"- [ ] **`test_solver`**\n",
"\n",
"which can be used to test that implementation has all necessary things defined.\n",
"\n",
"\n",
"Poisson equation is used in this notebook to demonstrate all phases of development process. It's weak form is: find $u\\in\\mathcal{U}$ such that\n",
"\\begin{equation}\n",
" \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x=\\int_{\\Omega}fv\\,\\mathrm{d}x+\\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s\\quad\\forall v\\in\\mathcal{V}.\n",
"\\end{equation}\n",
"\n",
"## Definitions / nomenclature\n",
"\n",
"- weighted residual form: differential equation is multiplied by a weight function and integrating it.\n",
"- weak form, \"principle of virtual work\": typically established by partial integration of a weighted residual form. $\\delta\\mathcal{W}_{\\mathrm{int}}=\\delta\\mathcal{W}_{\\mathrm{ext}}$.\n",
"- variational form, \"principle of minimum potential energy\": there exists some functional or \"potential function\" $\\Pi$ we are minimizing"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Developing own element\n",
"\n",
"Element contains basis functions and one or several fieldsets.\n",
"\n",
"Minimum requirements for element:\n",
"- subclass it from Element, unless want to implement everything by youself\n",
"- define basis and partial derivatives of it, because we need to interpolate over it\n",
"- give connectivity information, how this element is connected to other elements\n",
"\n",
"Test the element using ``test_element`` function. It it passes, the element implementation should be fine."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM: Element, Basis, FieldSet"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the basic implementation for element:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"type MyQuad4 <: Element\n",
" connectivity :: Array{Int, 1}\n",
" basis :: Basis\n",
" fields :: FieldSet\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Default constructor takes only connectivity information as input argument. Most important thing here is to define `Basis` which is used to interpolate fields."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"MyQuad4"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function MyQuad4(connectivity)\n",
" h(xi) = 1/4*[(1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) (1+xi[1])*(1+xi[2]) (1-xi[1])*(1+xi[2])]\n",
" dh(xi) = 1/4*[\n",
" -(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2])\n",
" -(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])]\n",
" basis = Basis(h, dh)\n",
" MyQuad4(connectivity, basis, Dict())\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also some basic information like number of basis functions and element dimension is needed. We can set them directly to the `Base.size`:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"size (generic function with 81 methods)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Base.size(element::Type{MyQuad4}) = (2, 4) # 2 => (ξ₁, ξ₂), 4 => 4 basis functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here comes one important thing. We always define our \"things\" so that the first index is dimension, like $(x, y, z)$ or $(\\xi_1, \\xi_2, \\xi_3)$ and second index is basis function number / node id or something similar to that. To motivate this, consider the following example:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Int64,2}:\n",
" 1 3 5 7\n",
" 2 4 6 8"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = [1 2; 3 4; 5 6; 7 8]'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, if we consider $x$ as of some field e.g. geometry, our coordinates of first node is $(1, 2)$, second is $(3, 4)$ and so on. Typically on vector field problems the global assembly is something like $(u_1, v_1, u_2, v_2, \\ldots, )$. If fields are defined this way we can easily flatten matrix to vector and back:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8-element Array{Int64,1}:\n",
" 1\n",
" 2\n",
" 3\n",
" 4\n",
" 5\n",
" 6\n",
" 7\n",
" 8"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dim = size(x)\n",
"x2 = vec(x)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x4 Array{Int64,2}:\n",
" 1 3 5 7\n",
" 2 4 6 8"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"reshape(x2, dim)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Also:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 2 3 4 5 6 7 8 "
]
}
],
"source": [
"for i=1:8\n",
" print(\"$(x[i]) \")\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So when keeping this order, things automatically go right in Julia.\n",
"\n",
"Next we check that everything is well defined:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: Testing element MyQuad4\n",
"INFO: element dimension: 2 x 4\n",
"INFO: Initializing element\n",
"INFO: basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n",
"INFO: field val at [0.0,0.0]: 0.0\n"
]
},
{
"ename": "LoadError",
"evalue": "LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
"output_type": "error",
"traceback": [
"LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
"",
" in call at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:147",
" in test_element at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:57"
]
}
],
"source": [
"using JuliaFEM: test_element\n",
"test_element(MyQuad4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If `test_element` passes, elements *interface* should be well defined. After building element, one can interpolate things in it. Here we create three new fieldsets `temperature`, `geometry` and `heat coefficient`, add some values for them in time $t=0.0$ and $t=1.0$ and interpolate:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Dict{ASCIIString,JuliaFEM.FieldSet} with 4 entries:\n",
" \"geometry\" => JuliaFEM.FieldSet(\"geometry\",JuliaFEM.Field[JuliaFEM.Fi…\n",
" \"heat coefficient\" => JuliaFEM.FieldSet(\"heat coefficient\",JuliaFEM.Field[Jul…\n",
" \"displacement\" => JuliaFEM.FieldSet(\"displacement\",JuliaFEM.Field[JuliaFE…\n",
" \"temperature\" => JuliaFEM.FieldSet(\"temperature\",JuliaFEM.Field[JuliaFEM…"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: FieldSet, Field\n",
"element = MyQuad4([1, 2, 3, 4])\n",
"\n",
"geometry_field = Field(0.0, Vector[]) # Create empty field at time t=0.0\n",
"push!(geometry_field, [ 0.0, 0.0]) # push some values for field\n",
"push!(geometry_field, [ 1.0, 0.0])\n",
"push!(geometry_field, [ 1.0, 1.0])\n",
"push!(geometry_field, [ 0.0, 1.0])\n",
"geometry_fieldset = FieldSet(\"geometry\") # create fieldset \"geometry\"\n",
"push!(geometry_fieldset, geometry_field) # add field to fieldset\n",
"element[\"geometry\"] = geometry_fieldset # add fieldset to element\n",
"\n",
"temperature_fieldset = FieldSet(\"temperature\")\n",
"push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
"push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
"element[\"temperature\"] = temperature_fieldset\n",
"\n",
"displacement_fieldset = FieldSet(\"displacement\")\n",
"push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))\n",
"push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))\n",
"element[\"displacement\"] = displacement_fieldset\n",
"\n",
"heat_coefficient_fieldset = FieldSet(\"heat coefficient\")\n",
"push!(heat_coefficient_fieldset, Field(0.0, 2))\n",
"push!(heat_coefficient_fieldset, Field(1.0, 3))\n",
"element[\"heat coefficient\"] = heat_coefficient_fieldset\n",
"element.fields"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.5\n",
" 0.5"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# geometry midpoint of element\n",
"using JuliaFEM: get_basis, grad\n",
"basis = get_basis(element)\n",
"basis(\"geometry\", [0.0, 0.0], 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x2 Array{Float64,2}:\n",
" 0.0 1.0"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolate derivatives works too\n",
"dbasis = grad(basis)\n",
"dbasis(\"temperature\", [0.0, 0.0], 0.5) # temperature gradient at mid point of element at time t=0.5"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.5"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolating scalar -> scalar.\n",
"basis(\"heat coefficient\", [0.0, 0.0], 0.5)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# set time to Inf to get very last field value and -Inf to get first one.\n",
"basis(\"heat coefficient\", [0.0, 0.0], Inf)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x2 Array{Float64,2}:\n",
" 0.132813 0.0703125\n",
" 0.0703125 0.0078125"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# some continuum mechanics\n",
"∇u = dbasis(\"displacement\", [0.0, 0.0], 1.0)\n",
"ɛ = 1/2*(∇u + ∇u')\n",
"Ω = 1/2*(∇u - ∇u')\n",
"X = basis(\"geometry\", [0.0, 0.0], 1.0)\n",
"F = I + ∇u\n",
"C = F'*F\n",
"E = 1/2*(F'*F - I) # Green-Lagrange strain tensor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Basically one could already write local stiffness matrices:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"(\n",
"4x4 Array{Float64,2}:\n",
" 4.0 -1.0 -2.0 -1.0\n",
" -1.0 4.0 -1.0 -2.0\n",
" -2.0 -1.0 4.0 -1.0\n",
" -1.0 -2.0 -1.0 4.0,\n",
"\n",
"4x4 Array{Float64,2}:\n",
" 4.0 2.0 1.0 2.0\n",
" 2.0 4.0 2.0 1.0\n",
" 1.0 2.0 4.0 2.0\n",
" 2.0 1.0 2.0 4.0,\n",
"\n",
"4x1 Array{Float64,2}:\n",
" 1.0\n",
" 1.0\n",
" 1.0\n",
" 1.0)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"weights = [1.0, 1.0, 1.0, 1.0]\n",
"integration_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]\n",
"\n",
"M = zeros(4, 4)\n",
"K = zeros(4, 4)\n",
"f = zeros(4)\n",
"k = 6\n",
"rho = 36\n",
"q = 4\n",
"\n",
"for (w, xi) in zip(weights, integration_points)\n",
" N = get_basis(element)\n",
" ∇N = grad(N) # gradient is with respect to \"geometry\" field\n",
" detJ = det(N)(xi)\n",
" M += w*rho*N(xi)'*N(xi)*detJ # mass matrix\n",
" K += w*k*∇N(xi)'*∇N(xi)*detJ # stiffness matrix\n",
" f += w*N(xi)'*q*detJ # force vector\n",
"end\n",
"K, M, f"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"but we have something better than that. Keep on reading!\n",
"\n",
"### Summary of developing own elements\n",
"\n",
"Element itself if not calculating anything but only stores fields and basis functions so that the fields can be interpolated. We will provide command `test_element` which will ensure that everything necessary is defined. While a lot of things needs to be defined, by subclassing from `Element` most of these are already defined, thanks to multiple dispatch."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Developing own equation\n",
"\n",
"Let's consider a Poisson equation\n",
"\\begin{align}\n",
"-\\nabla \\cdot \\left( k \\nabla u\\right) &= f && \\text{on } \\Omega \\\\\n",
"\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}} \\\\\n",
"u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
"\\end{align}\n",
"\n",
"Weak form is, find $u\\in\\mathcal{U}$ such that\n",
"\\begin{equation}\n",
" \\int_{\\Omega}k \\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Omega}fv\\,\\mathrm{d}x + \\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s \\quad\\forall v\\in\\mathcal{V}.\n",
"\\end{equation}\n",
"\n",
"Minimum requirements for equation: \n",
"- subclass from Equation, if not want to implement from scratch\n",
"- provide the name of the unknown field variable trying to solve\n",
"- default constructor takes the element as input argument\n",
"\n",
"This time we will have function `test_equation`, which we can use to test that everything is working as expected. "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM: Equation, IntegrationPoint, Quad4, Seg2\n",
"\n",
"abstract Heat <: Equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Basic data type needs to contain element, integration points and global degrees of freedom for global assembly:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"size (generic function with 76 methods)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" Diffusive heat transfer for 4-node bilinear element. \"\"\"\n",
"type DC2D4 <: Heat\n",
" element :: Quad4\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"function DC2D4(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",
" if !haskey(element, \"temperature\")\n",
" element[\"temperature\"] = FieldSet()\n",
" end\n",
" DC2D4(element, integration_points)\n",
"end\n",
"Base.size(equation::DC2D4) = (1, 4)\n",
"\n",
"\"\"\" Diffusive heat transfer for 2-node linear segment. \"\"\"\n",
"type DC2D2 <: Heat\n",
" element :: Seg2\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"function DC2D2(element::Seg2)\n",
" integration_points = [\n",
" IntegrationPoint([0.0], 2.0)]\n",
" if !haskey(element, \"temperature\")\n",
" element[\"temperature\"] = FieldSet()\n",
" end\n",
" DC2D2(element, integration_points)\n",
"end\n",
"Base.size(equation::DC2D2) = (1, 2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the actual implementation for $\\int_{\\Omega}k \\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Omega}fv\\,\\mathrm{d}x + \\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s$.\n",
"\n",
"The low-level way is to write equations for local matrix representation directly:"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_force_vector (generic function with 3 methods)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_basis, get_element\n",
"\n",
"\"\"\" Mass matrix for dynamical problems. Not used in this example. \"\"\"\n",
"function JuliaFEM.get_mass_matrix(equation::DC2D4, ip, time)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" ρ = basis(\"density\")\n",
" return ρ(ip, time) * basis(ip, time)'*basis(ip, time)\n",
"end\n",
"\"\"\" Left hand side defined in integration point. \"\"\"\n",
"function JuliaFEM.get_stiffness_matrix(equation::DC2D4, ip, time)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" dbasis = grad(basis)\n",
" k = basis(\"temperature thermal conductivity\")\n",
" return k(ip, time) * dbasis(ip, time)'*dbasis(ip, time)\n",
"end\n",
"\"\"\" Right hand side defined in integration point. \"\"\"\n",
"function JuliaFEM.get_force_vector(equation::DC2D4, ip, time)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" f = basis(\"temperature load\")\n",
" return f(ip, time)*basis(ip, time)'\n",
"end\n",
"JuliaFEM.has_mass_matrix(equation::DC2D4) = true\n",
"JuliaFEM.has_stiffness_matrix(equation::DC2D4) = true\n",
"JuliaFEM.has_force_vector(equation::DC2D4) = true\n",
"\n",
"\"\"\" Right hand side defined in integration point. \"\"\"\n",
"function JuliaFEM.get_force_vector(equation::DC2D2, ip, time)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" g = basis(\"temperature flux\")\n",
" return g(ip, time)*basis(ip, time)'\n",
"end\n",
"JuliaFEM.has_force_vector(equation::DC2D2) = true"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It might at start to feel a bit ackward to pass integration point and time all around as function arguments, so a little reasoning why this is essential might be needed.\n",
"\n",
"First of all, equations are written in integration or gauss quadrature points for a certain reason. The reason is that for some incremental constitutive models internal variables are needed, which are saved to FieldSets inside integration points. We can always access these variables through `ip.fields` in the same way we can access all other field variables defined inside element via `element.fields`.\n",
"\n",
"Secondly, starting point for our equation is that everything can be time-dependent. FieldSet is, like name suggests, a set of fields with different times and/or increments. By transferring time parameter inside equation makes it possible to access not only current field values but earlier values too. In this example field variables thermal conductivity $k$ and load $f$ are constant in time but we are planning to set flux $g$ such a way that it is linear ramp from $g(t) = 6t$.\n",
"\n",
"If we want to try out this formulation, we must create element and assign this equation for it:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# if no time is given when defining fields it defaults to 0.0.\n",
"# create volume element\n",
"element = Quad4([1, 2, 3, 4])\n",
"fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
"fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n",
"fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n",
"fieldset4 = FieldSet(\"density\", [Field(36.0)])\n",
"element[\"geometry\"] = fieldset1\n",
"element[\"temperature thermal conductivity\"] = fieldset2\n",
"element[\"temperature load\"] = fieldset3\n",
"element[\"density\"] = fieldset4\n",
"equation = DC2D4(element)\n",
"\n",
"# create boundary element with \n",
"boundary_element = Seg2([1, 2])\n",
"boundary_element[\"geometry\"] = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])])\n",
"# linear ramp from 1 to 6 in time 0 to 1\n",
"boundary_element[\"temperature flux\"] = FieldSet(\"temperature flux\", [Field(0.0, 0.0), Field(1.0, 6.0)])\n",
"boundary_equation = DC2D2(boundary_element);"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x4 Array{Float64,2}:\n",
" 4.0 2.0 1.0 2.0\n",
" 2.0 4.0 2.0 1.0\n",
" 1.0 2.0 4.0 2.0\n",
" 2.0 1.0 2.0 4.0"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n",
"\n",
"local_assembly = initialize_local_assembly()\n",
"calculate_local_assembly!(local_assembly, equation, \"temperature\", 0.0)\n",
"local_assembly.mass_matrix"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x4 Array{Float64,2}:\n",
" 4.0 -1.0 -2.0 -1.0\n",
" -1.0 4.0 -1.0 -2.0\n",
" -2.0 -1.0 4.0 -1.0\n",
" -1.0 -2.0 -1.0 4.0"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly.stiffness_matrix"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x1 Array{Float64,2}:\n",
" 3.0\n",
" 3.0\n",
" 3.0\n",
" 3.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly.force_vector"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's set homogeneous Dirichlet boundary condition ($u=0$) on boundary on $\\Gamma_\\mathrm{D} = \\{(x,y) | 0 \\leq x \\leq 1, y=1\\}$. With source term $f=12$ and $k=6$ we should have constant $T=1$ on free boundary. According to my math we should have $u\\left(x,y\\right)=-\\frac{1}{6}\\left(\\frac{1}{2}fx^{2}-fx\\right)$ which equals to one if $x=1$."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = local_assembly.stiffness_matrix\n",
"b = local_assembly.force_vector\n",
"u = zeros(4)\n",
"fdofs = [1,2]\n",
"u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Notice that we created a separate boundary element which integrates heat flux. Defining boundary elements goes identically with other elements. With flux $g=6$ we should have constant $T=1$ on free boundary. This should equal to $u(x,y) = x$, i.e. 1d bar stretched from other end while the other one is fixed."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly = initialize_local_assembly()\n",
"calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\")\n",
"b = zeros(4)\n",
"b[fdofs] = local_assembly.force_vector\n",
"u = zeros(4)\n",
"u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some quick notes. Passing time all around is not strictly needed and defaults to $t=+\\inf$, when the most recent field is picked from `FieldSet`. Also when defining `Field` without time, it defaults to $t=0$. It's however good practice to interpolate fields inside functions using time parameter which makes the equations more general. If we want to solve the last equation in time $t=0.5$ we should get $T=0.5$:"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.5\n",
" 0.5"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly = initialize_local_assembly()\n",
"calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\", 0.5) # <-- time added\n",
"b = zeros(4)\n",
"b[fdofs] = local_assembly.force_vector\n",
"u = zeros(4)\n",
"u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Some alternative ways to define equations - principle of minimum potential energy\n",
"\n",
"Fundamentally we are always looking for a solution of systems of equation $Ax = b$, despite the fact that mathematicians doesn't like this \"engineering approach\" at all. For some equations it may be too hard to write left hand side, i.e., stiffness matrix $A$. (At least for me). For this reason JuliaFEM supports automatic differentiation using package called `ForwardDiff`, which makes the analytical linearization unnecessary.\n",
"\n",
"Let's consider the following functional\n",
"\\begin{equation}\n",
"\\min\\, J\\left(u\\right)=\\int_{\\Omega}(k+6u)\\left(\\nabla u \\cdot \\nabla u\\right) \\, \\mathrm{d}x-Pu,\n",
"\\end{equation}\n",
"where $P$ contains point loads at the free corners of the domain. This is nothing more but a redefinition of the earlier example problem defined such that the source term $k$, which in previous example was constant, is no more constant, but is replaced a term which depend from field $u$: $q(u) = k + 6u$.\n",
"\n",
"This functional doesn't actually have any or very little physical meaning, but it's non-linear due to the source term $q(u) = k+6u$ and therefore suits to our needs in demonstration purposes. In mechanical world this could be some sort of elastic spring with a spring \"constant\" depending on displacement (non-linear elasticity problem), which is compressed using point force $P$. The system is still *conservative*, i.e. it has a potential energy and can be solved using standard variational methods, finding a minimum of potential energy in suitable finite set of functions.\n",
"\n",
"In JuliaFEM we simply define it's potential energy functional $\\Pi(u)$ and let the `ForwardDiff`-package to handle the rest of the hard things. Because the problem is nonlinear, linearization and iterative solution method (=Newton's method) is needed. Accurate solution at the free end is\n",
"\\begin{equation}\n",
"u = \\frac{2}{3}.\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"size (generic function with 77 methods)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" Diffusive heat transfer for 4-node bilinear element, with a nonlinear source term. \"\"\"\n",
"type DC2D4NL <: Heat\n",
" element :: Quad4\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"function DC2D4NL(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",
" if !haskey(element, \"temperature\")\n",
" element[\"temperature\"] = FieldSet()\n",
" # Initial configuration needs to be defined if using autodiff\n",
" push!(element[\"temperature\"], Field([0.0, 0.0, 0.0, 0.0]))\n",
" end\n",
" DC2D4NL(element, integration_points)\n",
"end\n",
"Base.size(equation::DC2D4NL) = (1, 4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This time we don't write stiffness matrix, etc. explicitly but just write potential energy and let `ForwardDiff` take care of rest of the stuff."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_potential_energy (generic function with 2 methods)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" Calculate a potential Π = Wint - Wext of system. \"\"\"\n",
"function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" k = basis(\"temperature thermal conductivity\", ip, time)\n",
" f = basis(\"temperature load\", ip, time)\n",
" T = basis(\"temperature\", ip, time, variation)\n",
" ∇T = grad(basis)(\"temperature\", ip, time, variation)\n",
" Wint = (k + 6*T) * 1/2*vecdot(∇T, ∇T)\n",
" Wext = f*T\n",
" return Wint - Wext\n",
"end\n",
"JuliaFEM.has_potential_energy(eq::DC2D4NL) = true"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solving non-linear variational problem is rather complicated thing due to the nature of non-linearity. Here's a heavily commented version of simple Newton iteration."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"increment 1, |du| = 1.41421, |r| = 2.121\n",
"increment 2, |du| = 0.42426, |r| = 2.121\n",
"increment 3, |du| = 0.04657, |r| = 2.121\n",
"increment 4, |du| = 0.00057, |r| = 2.121\n",
"increment 5, |du| = 0.00000, |r| = 2.121\n",
"elapsed time: 0.617104273 seconds\n",
"error: 1.5543122344752192e-15\n",
"temperature at free end: [0.6666666666666682,0.6666666666666682], should be 0.6666666666666666\n"
]
}
],
"source": [
"function run_simulation()\n",
" # create model -- start\n",
" element = Quad4([1, 2, 3, 4])\n",
" fieldset1 = FieldSet(\"geometry\")\n",
" field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
" push!(fieldset1, field1)\n",
" fieldset2 = FieldSet(\"temperature thermal conductivity\")\n",
" push!(fieldset2, Field(0.0, 6.0))\n",
" fieldset3 = FieldSet(\"temperature load\")\n",
" push!(fieldset3, Field(0.0, 0*[12.0, 12.0, 12.0, 12.0]))\n",
" fieldset4 = FieldSet(\"temperature nodal load\")\n",
" push!(fieldset4, Field(0.0, [3.0, 3.0, 0.0, 0.0])) # <-- P is defined here\n",
" element[\"geometry\"] = fieldset1\n",
" element[\"temperature thermal conductivity\"] = fieldset2\n",
" element[\"temperature load\"] = fieldset3\n",
" element[\"temperature nodal load\"] = fieldset4\n",
" equation = DC2D4NL(element)\n",
" # create model -- end\n",
"\n",
" T0 = element[\"temperature\"][1] # initial temperature field, we need something to \"variate\"\n",
" la = initialize_local_assembly() # create workspace for local matrices\n",
" T = zeros(4) # create workspace for solution vector\n",
" ΔT = zeros(4) # \n",
" fd = [1, 2]\n",
" tic()\n",
" # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n",
" for i=1:5\n",
" calculate_local_assembly!(la, equation, \"temperature\") # calculate local matrices\n",
" ΔT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd] # <-- note sign convention, more on this below\n",
" T = T + ΔT # add increment to previous value\n",
" # create a new field \"similar\" to field T0 (i.e., same dimension of field variable with new data)\n",
" new_field = similar(T0, T)\n",
" new_field.time = 1.0\n",
" new_field.increment = i\n",
" push!(element[\"temperature\"], new_field) # add new field to \"temperature\" fieldset of element\n",
" # print some convergence information\n",
" @printf(\"increment %2d, |du| = %8.5f, |r| = %8.3f\\n\", i, norm(ΔT), norm(b[fd]))\n",
" end\n",
" toc()\n",
" err = element[\"temperature\"][end].values[1] - 2/3\n",
" println(\"error: $err\")\n",
" blaa = element[\"temperature\"][end].values[1:2]\n",
" acc = 2/3\n",
" println(\"temperature at free end: $blaa, should be $acc\")\n",
"end\n",
"run_simulation()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### A short note about sign convention\n",
"\n",
"When solving a non-linear equations like $\\mathbf{r}\\left(\\mathbf{u}\\right) = 0$ it is necessary to iterate, i.e. solve several linearized solutions using e.g. Newton iterations until residual is small enough. Using truncated Taylor series the solution algorithm is in general\n",
"\\begin{equation}\n",
"\\mathbf{r}\\left(\\mathbf{u}+\\Delta\\mathbf{u}\\right)=\\mathbf{r}\\left(\\mathbf{u}\\right)+\\frac{\\partial\\mathbf{r}\\left(\\mathbf{u}\\right)}{\\partial\\mathbf{u}}\\Delta\\mathbf{u}+\\cdots\\approx\\mathbf{r}\\left(\\mathbf{u}\\right)+\\mathbf{K}\\Delta\\mathbf{u}=0\n",
"\\end{equation}\n",
"\\begin{equation}\n",
"\\mathbf{K}\\Delta\\mathbf{u}=-\\mathbf{r}\\left(\\mathbf{u}\\right)\n",
"\\end{equation}\n",
"with a *negative* sign in front of residual vector.\n",
"\n",
"We do however want to solve something like\n",
"\\begin{equation}\n",
"\\mathbf{A}\\Delta\\mathbf{x}^{\\left(i+1\\right)}=\\mathbf{b},\\qquad\\mathbf{x}^{\\left(i+1\\right)}=\\mathbf{x}^{\\left(i\\right)}+\\Delta\\mathbf{x}^{\\left(i+1\\right)},\n",
"\\end{equation}\n",
"where new increment is *added* to the old result. If this doesn't say anything, compare the equation to the 1d Newton iteration (https://en.wikipedia.org/wiki/Newton's_method)\n",
"\\begin{equation}\n",
"x_{n+1}=x_{n}-\\frac{f\\left(x_{n}\\right)}{f'\\left(x_{n}\\right)}\n",
"\\end{equation}\n",
"\n",
"Our convention is that the negative sign of residual vector is positive. It makes sense because this way the \"normal\" linear system force vector is like it's used to be in school books. Typically the residual vector is defined $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$. In our case the residual vector must be multiplied by $-1$ after the linearization so that the increment can be *added* to old result. In some text books this has already be taken into account by defining right hand side as $\\mathbf{r}=\\mathbf{f}_{\\mathrm{ext}}-\\mathbf{f}_{\\mathrm{int}}$ while linearization is still done for vector $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}$. Quite confusing!\n",
"\n",
"Just keep in mind that there is always minus somewhere and it should be. In our variational implementation the sign of residual vector is automatically changed so that next increment is $\\mathbf{A} \\Delta\\mathbf{x}=\\mathbf{b}$. Like in the simple example above.\n",
"\n",
"### Principle of virtual work\n",
"\n",
"Defining variational form or \"energy form\", or in other words, using principle of minimum potential energy, requires system to be *conservative*, i.e. there is some energy functional describing the system. Indeed it's a very elegant approach to model e.g. hyperelasticity and other systems without a loss of energy. When system has energy dissipation, principle of minimum potential energy cannot be used. Possible situations when this happens include e.g. material plasticity or friction. In these situations the principle of virtual work is useful.\n",
"\n",
"In this situation we find equation $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$ which needs to be solved. Again we let the `ForwardDiff` do the linearization of the right hand side. We demonstrate this by constructing geometrically non-linear elasticity equation for 4-node bilinear element. Virtual work now looks\n",
"\n",
"\\begin{equation}\n",
"\\delta\\mathcal{W} := \\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\boldsymbol{b}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0} = 0\n",
"\\end{equation}\n",
"where\n",
"\\begin{equation}\n",
"\\delta\\mathcal{W}_{\\mathrm{int}}=\\int\\mathbf{S}:\\delta\\mathbf{E}\\,\\mathrm{d}V_{0}=\\int_{\\Omega_{0}}\\mathbf{P}\\cdot\\mbox{Grad}\\delta\\mathbf{u}\\,\\mathrm{d}V_{0}\n",
"\\end{equation}\n",
"and $\\boldsymbol{b}_{0}$ is a volume load defined in reference configuration."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_residual_vector (generic function with 3 methods)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_field\n",
"\n",
"abstract Elasticity <: Equation\n",
"abstract PlaneElasticity <: Elasticity\n",
"abstract PlaneStressElasticity <: PlaneElasticity\n",
"\n",
"\"\"\" Plane stress formulation for 4-node bilinear element. \"\"\"\n",
"type CPS4 <: PlaneStressElasticity\n",
" element :: Quad4\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"\n",
"function CPS4(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",
" if !haskey(element, \"displacement\")\n",
" element[\"displacement\"] = FieldSet()\n",
" # initial field must be defined if using autodiff\n",
" u0 = Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n",
" push!(element[\"displacement\"], u0)\n",
" end\n",
" CPS4(element, integration_points)\n",
"end\n",
"\n",
"JuliaFEM.size(eq::CPS4) = (2, 4)\n",
"\n",
"function JuliaFEM.get_residual_vector(equation::CPS4, ip, time; variation=nothing)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" dbasis = grad(basis)\n",
"\n",
" # material parameters\n",
" E = basis(\"youngs modulus\", ip, time)\n",
" ν = basis(\"poissons ratio\", ip, time)\n",
" μ = E/(2*(1+ν))\n",
" λ = E*ν/((1+ν)*(1-2*ν))\n",
" λ = 2*λ*μ/(λ + 2*μ) # <- correction for 2d\n",
"\n",
" # elasticity formulation\n",
" u = basis(\"displacement\", ip, time, variation)\n",
" ∇u = dbasis(\"displacement\", ip, time, variation)\n",
" F = I + ∇u\n",
" # F = F_e*F_p\n",
" b = basis(\"displacement volume load\", ip, time)\n",
" E = 1/2*(F'*F - I)\n",
" S = λ*trace(E)*I + 2*μ*E\n",
" J = det(F)\n",
" σ = J^-1 * F*S*F'\n",
" #push!(ip[\"cauchy stress\"], Field(time, σ))\n",
" P = F*S\n",
"\n",
" # residual vector\n",
" r_int = P*dbasis(ip,time)\n",
" r_ext = b*basis(ip,time)\n",
" r = r_int - r_ext\n",
" return vec(r)\n",
"end\n",
"JuliaFEM.has_residual_vector(equation::CPS4) = true"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There's a tiny Newton solver which can be used to test single element models to verify formulations. It makes the life a bit easier."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"displacement at tip: -8.773031198197748\n"
]
},
{
"data": {
"text/plain": [
"Success :: (line:-1) :: fact was true\n",
" Expression: disp --> roughly(-8.77303119819776)\n",
" Expected: -8.77303119819776\n",
" Occurred: -8.773031198197748"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_integration_points, solve!\n",
"\n",
"function run_simulation_2()\n",
" # create model -- start\n",
" # fields can be created less verbose way. just keep in mind that if time is not defined\n",
" # it defaults to 0.0. fields can be set \"directly\" to the element but it creates a fieldset\n",
" # anyway\n",
" element = Quad4([1, 2, 3, 4])\n",
" element[\"geometry\"] = Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])\n",
" element[\"youngs modulus\"] = Field(500.0)\n",
" element[\"poissons ratio\"] = Field(0.3)\n",
" element[\"displacement volume load\"] = Field(Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]])\n",
" equation = CPS4(element)\n",
" # create model -- end\n",
"\n",
" free_dofs = [3, 4, 5, 6]\n",
" solve!(equation, \"displacement\", free_dofs) # launch a newton solver for single element\n",
" disp = get_basis(element)(\"displacement\", [1.0, 1.0])[2]\n",
" println(\"displacement at tip: $disp\")\n",
" # verified using Code Aster.\n",
" @fact disp --> roughly(-8.77303119819776E+00)\n",
"end\n",
"run_simulation_2()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Manually assembling local matrices\n",
"\n",
"There might still be situations where all the above methods writing field equations are just not enough. This situation can happen for example when calculating tangent stiffness matrix analytically for a nonlinear problem. In this case all matrices are usually calculated at the same time. Or maybe for educational purposes it's important to show how matrices are actually calculated. Or for debugging. \n",
"\n",
"Anyway, it's possible to override `calculate_local_assembly!` function for your own equation and after that access all the low level stuff. Here's an example how to do that. "
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Success :: (line:-1) :: fact was true\n",
" Expression: la.stiffness_matrix[free_dofs,free_dofs] \\ la.force_vector[free_dofs] --> roughly([1.0,1.0])\n",
" Expected: [1.0,1.0]\n",
" Occurred: [1.0000000000000002,1.0000000000000002]"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: LocalAssembly, get_integration_points, initialize_local_assembly!\n",
"\n",
"function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::DC2D4,\n",
" unknown_field_name::ASCIIString, time::Number=Inf,\n",
" problem=nothing) # <= notice problem\n",
"\n",
" # zero all workspace. always start with this command, otherwise\n",
" # local matrices are added to the old ones.\n",
" initialize_local_assembly!(assembly, equation)\n",
"\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" dbasis = grad(basis)\n",
" detJ = det(basis)\n",
" for ip in get_integration_points(equation)\n",
" w = ip.weight * detJ(ip)\n",
" # evaluate fields in integration point\n",
" ρ = basis(\"density\", ip, time)\n",
" k = basis(\"temperature thermal conductivity\", ip, time)\n",
" f = basis(\"temperature load\", ip, time)\n",
" # evaluate basis functions and gradient in integration point\n",
" N = basis(ip, time)\n",
" dN = dbasis(ip, time)\n",
" # do assembly\n",
" assembly.mass_matrix += w * ρ*N'*N\n",
" assembly.stiffness_matrix += w * k*dN'*dN\n",
" assembly.force_vector += w * N'*f\n",
" end\n",
"end\n",
"\n",
"function test_local_assembly()\n",
" # create volume element\n",
" element = Quad4([1, 2, 3, 4])\n",
" # fields can be created even more directly:\n",
" element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
" element[\"temperature thermal conductivity\"] = 6.0\n",
" element[\"temperature load\"] = [12.0, 12.0, 12.0, 12.0]\n",
" element[\"density\"] = 36.0\n",
" equation = DC2D4(element)\n",
" la = initialize_local_assembly()\n",
" calculate_local_assembly!(la, equation, \"temperature\")\n",
" free_dofs = [1, 2]\n",
" @fact la.stiffness_matrix[free_dofs, free_dofs] \\ la.force_vector[free_dofs] --> roughly([1.0, 1.0])\n",
"end\n",
"\n",
"test_local_assembly()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that function takes also problem as input, which allows to access to rest of the equations and elements. From that in the next section."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Defining own problem\n",
"\n",
"The main purpose of `Problem` is to create a mapping between element types and field equations, so that solver knows how to construct equations. Loosely speaking `Problem` is a set of equations and instructions how they can be created from elements. Every time new element is pushed to problem, element mapping creates corresponding equations ready for integration.\n",
"\n",
"Notice that there is no any limits of how many different equations can be set to one element, so it's totally possible for example solve displacement and temperature equations with same element and even so that problems are loosely coupled, i.e. using temperature field to calculate stresses in mechanical problem."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"PlaneHeatProblem"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Problem, get_equations, get_connectivity\n",
"using JuliaFEM: get_element, get_unknown_field_dimension\n",
"\n",
"type PlaneHeatProblem <: Problem\n",
" unknown_field_name :: ASCIIString\n",
" unknown_field_dimension :: Int\n",
" equations :: Array{Equation, 1}\n",
" element_mapping :: Dict{DataType, DataType}\n",
"end\n",
"function PlaneHeatProblem(equations=[])\n",
" element_mapping = Dict(\n",
" Quad4 => DC2D4,\n",
" Seg2 => DC2D2)\n",
" PlaneHeatProblem(\"temperature\", 1, equations, element_mapping)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Developing own solver\n",
"\n",
"Solver is the part where all numerical crunching happens. Typically solver takes a set of problems. Each problem has all necessary information to get assembled. Solver then constructs the final problem to be solved. Typically this problem is some field problem + dirichlet boundaries but it can also be e.g. mixed problem. One essential task of solver is after solution to update corresponding fields to elements. Here we give a simple (yet parallel solution inside one machine) demo solver which solves a typical problem of one field and one dirichlet boundary. It isn't optimized any way but is appropriate as simple test solver."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"call (generic function with 1319 methods)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Solver, get_problems, get_unknown_field_name, get_gdofs\n",
"using JuliaFEM: initialize_global_assembly, calculate_global_assembly!\n",
"\n",
"\"\"\" Simple solver for educational purposes. \"\"\"\n",
"type SimpleSolver <: Solver\n",
" problems :: Array{Problem, 1}\n",
"end\n",
"\n",
"\"\"\" Default initializer. \"\"\"\n",
"function SimpleSolver()\n",
" SimpleSolver([])\n",
"end\n",
"\n",
"\"\"\"\n",
"Call solver to solve a set of problems.\n",
"\n",
"This is a simple direct solver for demonstration purposes. It handles the\n",
"common situation, i.e., some main field problem and it's Dirichlet boundary.\n",
"\n",
" Au + C'λ = f\n",
" Cu = g\n",
"\n",
"\"\"\"\n",
"function call(solver::SimpleSolver, time::Number=Inf)\n",
" p1, p2 = get_problems(solver)\n",
"\n",
" ga1 = initialize_global_assembly(p1)\n",
" calculate_global_assembly!(ga1, p1)\n",
" ga2 = initialize_global_assembly(p2)\n",
" calculate_global_assembly!(ga2, p2)\n",
"\n",
" A1 = ga1.stiffness_matrix\n",
" b1 = ga1.force_vector\n",
" A2 = ga2.stiffness_matrix\n",
" b2 = ga2.force_vector\n",
" \n",
" # create a saddle point problem\n",
" A = [A1 A2; A2' zeros(A2)]\n",
" b = [b1; b2]\n",
"\n",
" # solve problem\n",
" nz = unique(rowvals(A)) # here we remove any zero rows\n",
" x = zeros(b)\n",
" x[nz] = lufact(A[nz,nz]) \\ full(b[nz])\n",
"\n",
" # get \"problem-wise\" solution vectors\n",
" x1 = x[1:length(b1)]\n",
" x2 = x[length(b1)+1:end]\n",
"\n",
" # update field for elements in problem 1\n",
" for equation in get_equations(p1)\n",
" element = get_element(equation)\n",
" field_name = get_unknown_field_name(p1)\n",
" gdofs = get_gdofs(p1, equation)\n",
" element_solution = full(x1[gdofs])\n",
" field = Field(time, element_solution)\n",
" push!(element[field_name], field)\n",
" end\n",
"\n",
" # update field for elements in problem 2 (Dirichlet boundary)\n",
" for equation in get_equations(p2)\n",
" element = get_element(equation)\n",
" field_name = get_unknown_field_name(p2)\n",
" gdofs = get_gdofs(p2, equation)\n",
" element_solution = full(x2[gdofs])\n",
" field = Field(time, element_solution)\n",
" push!(element[field_name], field)\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be seen that we could already use parallel processing to assemble each problem concurrently. Moreover, LU factorization is already parallelized inside one machine, thanks to BLAS.\n",
"\n",
"Next we collect all things together and solve Poisson equation in domain $\\Omega = [0,1]\\times[0,1]$ with some Neumann boundary on $\\Gamma_1$ and Dirichlet boundary on $\\Gamma_2$."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"28-Oct 04:27:34:INFO:root:assembling problem for temperature\n",
"28-Oct 04:27:34:INFO:root:dimension of unknown field: 1, problem dofs: 4\n",
"28-Oct 04:27:35:INFO:root:assembling problem for reaction force\n",
"28-Oct 04:27:35:INFO:root:dimension of unknown field: 1, problem dofs: 4\n",
"28-Oct 04:27:35:INFO:root:Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n"
]
},
{
"data": {
"text/plain": [
"Success :: (line:-1) :: fact was true\n",
" Expression: T --> roughly(100.0)\n",
" Expected: 100.0\n",
" Occurred: 100.00000000000003"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" Define Problem 1:\n",
"\n",
"- Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}\n",
"- Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁\n",
"\"\"\"\n",
"function get_heatproblem()\n",
" el1 = Quad4([1, 2, 3, 4])\n",
" # these might look like normal values but believe me, they\n",
" # are fields with temporal and spatial dimension\n",
" el1[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
" el1[\"temperature thermal conductivity\"] = 6.0\n",
" el1[\"temperature load\"] = zeros(4)\n",
" el1[\"density\"] = 36.0\n",
"\n",
" el2 = Seg2([1, 2])\n",
" el2[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0]]\n",
" # Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n",
" # yet another simplification, if field is given as a tuple,\n",
" # multiple fields are created. there is 1 second time step between\n",
" # each field. So the following is basically same as\n",
" # fieldset = FieldSet(\"temperature flux\")\n",
" # field1 = Field(0.0, 0.0)\n",
" # field2 = Field(1.0, 600.0)\n",
" # push!(fieldset, field1)\n",
" # push!(fieldset, field2)\n",
" # element[\"temperature flux\"] = fieldset\n",
" el2[\"temperature flux\"] = (0.0, 600.0)\n",
"\n",
" problem1 = PlaneHeatProblem()\n",
" push!(problem1, el1)\n",
" push!(problem1, el2)\n",
" return problem1\n",
"end\n",
"\n",
"using JuliaFEM: DirichletProblem\n",
"\n",
"\"\"\" Define Problem 2:\n",
" - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n",
"\"\"\"\n",
"function get_boundaryproblem()\n",
" el3 = Seg2([3, 4])\n",
" el3[\"geometry\"] = Vector[[1.0, 1.0], [0.0, 1.0]]\n",
" problem2 = DirichletProblem(1)\n",
" push!(problem2, el3)\n",
" return problem2\n",
"end\n",
"\n",
"function main()\n",
" problem1 = get_heatproblem()\n",
" problem2 = get_boundaryproblem()\n",
" # Create a solver for a set of problems\n",
" solver = SimpleSolver()\n",
" push!(solver, problem1)\n",
" push!(solver, problem2)\n",
" # Solve problem at time t=1.0 and update fields\n",
" call(solver, 1.0)\n",
"\n",
" # Postprocess.\n",
" # Interpolate temperature field along boundary of Γ₁ at time t=1.0\n",
" xi = [0.0, -1.0]\n",
" el2 = get_element(problem1.equations[2])\n",
" basis = get_basis(el2)\n",
" X = basis(\"geometry\", xi, 1.0)\n",
" T = basis(\"temperature\", xi, 1.0)\n",
" Logging.info(\"Temperature at point X = $X is T = $T\")\n",
" @fact T --> roughly(100.0)\n",
"end\n",
"\n",
"main()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this notebook the basic instructions how to develop JuliaFEM has been given. The most imporant concepts has been considered; how to develop own element with own basis, several ways how to define own equation, and how to finally assemble and calculate the problem using solver. Any comments and/or discussion about technical details, theory, programming, or from life in general is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Advanced stuff\n",
"\n",
"In last section of this tutorial we consider some of the more advanced things which may araise when developing own models.\n",
"\n",
"### Boundary element access to parent element + overriding equations in problems\n",
"\n",
"This kind of situation might happen when one is almost happy for some problem setting, but would like to change just one or two equations from it. For example boundary equation is not satisfying all the requirements and one would like to test something new. \n",
"\n",
"### Accessing integration points\n",
"\n",
"### Fields as a function of something.\n",
"- statistical variables\n",
"- field dependent from another field\n",
"- field dependent from time or spatial domain etc.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"using JuliaFEM: get_default_integration_points\n",
"\"\"\" 2-node radiation boundary element. \"\"\"\n",
"type DC2D2RAD <: Heat\n",
" element :: Seg2\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"function DC2D2RAD(element::Seg2)\n",
" integration_points = [\n",
" IntegrationPoint([0.0], 2.0)]\n",
" if !haskey(element, \"temperature\")\n",
" element[\"temperature\"] = FieldSet()\n",
" push!(element[\"temperature\"], Field([0.0, 0.0]))\n",
" end\n",
" DC2D2RAD(element, integration_points)\n",
"end\n",
"Base.size(equation::DC2D2RAD) = (1, 2)\n",
"\n",
"\"\"\" Calculate potential energy caused by radiation.\n",
"https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_constant\n",
"\"\"\"\n",
"function JuliaFEM.get_potential_energy(equation::DC2D2RAD, ip, time; variation=nothing)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" eps = basis(\"emissivity\", ip, time)\n",
" #sig = basis(\"stefan-boltzmann constant\", ip, time)\n",
" sig = 5.670367e-8 # i guess stefan-boltzmann constant is constant ;)\n",
" T = basis(\"temperature\", ip, time, variation)\n",
" T_ext = basis(\"temperature external\", ip, time)\n",
" q = eps*sig*((T_ext+273.15)^4 - (T+273.15)^4)\n",
" println(ForwardDiff.value(q*T))\n",
" return q\n",
"end\n",
"JuliaFEM.has_potential_energy(equation::DC2D2RAD) = true\n",
"\n",
"Defining new equation mapping to old problem is one line command. Here we replace `DC2D2` $\\rightarrow$ `DC2DCRAD`\n",
"\n",
"function run_radiation_model()\n",
" # this is the same as before\n",
" element = Quad4([1, 2, 3, 4])\n",
" fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
" fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n",
" fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n",
" fieldset4 = FieldSet(\"density\", [Field(36.0)])\n",
" push!(element, fieldset1)\n",
" push!(element, fieldset2)\n",
" push!(element, fieldset3)\n",
" push!(element, fieldset4)\n",
"\n",
" # create boundary element\n",
" boundary_element = Seg2([1, 2])\n",
" push!(boundary_element, FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])]))\n",
" push!(boundary_element, FieldSet(\"emissivity\", [Field(0.5)]))\n",
" push!(boundary_element, FieldSet(\"temperature external\", [Field(20.0)]))\n",
"\n",
" # set initial conditions\n",
" push!(element, FieldSet(\"temperature\", [Field([0.0, 0.0, 0.0, 0.0])]))\n",
" push!(boundary_element, FieldSet(\"temperature\", [Field([0.0, 1.0])]))\n",
" \n",
" # create problem, change element mapping\n",
" problem = PlaneHeatProblem()\n",
" problem[Seg2] = DC2D2RAD # Seg2 was previous DC2D2\n",
" push!(problem, element)\n",
" push!(problem, boundary_element)\n",
"\n",
" # run our \"unit test solver\"\n",
" free_dofs = [3, 4]\n",
" solve!(problem, free_dofs; max_iterations=4, dump_matrices=true)\n",
" basis = get_basis(boundary_element)\n",
" T = basis(\"temperature\", [0.0])\n",
" println(\"Temperature at the midpoint of element: $T\")\n",
"end\n",
"\n",
"run_radiation_model()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.0",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Reference in New Issue View Git Blame Copy Permalink