Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-19 19:53:51 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
4e780b3ff54da7cd4efd8aa892ae1402bdd8a65a
JuliaFEM.jl/docs/tutorials/2015-08-29-developing-juliafem.ipynb
T
Jukka Aho 9b8842f5f3 updated developers guide + tests
2015-11-11 00:54:19 +02:00

1974 lines
81 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 aiming 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 the Poisson equation is used in examples to demonstrate the main concepts of the package. Poisson equation has some analogies to solid mechanics. Consider for example $(EAu')' + q = 0$, what 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 the 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 the following types are briefly described:\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 holding some field equation, like Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc., \"solvable\". This is the physics we are solving by approximating the equations in the element area by basis functionsprovided by element.\n",
"\n",
"**Problem**: object which maps equations to elements. For example, HeatProblem maps Poisson equation to Lagrange elements and ElasticityProblem maps elasticity equation to Lagrange elements.\n",
"\n",
"**Solver**: object 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 to point 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",
"what 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]: 2.5\n",
"INFO: derivative of basis at [0.0,0.0]:\n",
"[-0.5 0.5 0.5 -0.5\n",
" -0.5 -0.5 0.5 0.5]\n",
"INFO: field val at [0.0,0.0]: [0.0 2.0]\n",
"INFO: Element MyQuad4 passed tests.\n"
]
}
],
"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 fields `temperature`, `geometry` and `heat coefficient`, add some values for them in time $t=0.0$ and $t=1.0$ and interpolate. The concept of fields is described in [another notebook](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-06-14-data-structures.ipynb), and we are not particularly focusing on them here. Here's some examples how to store data to element."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Dict{ASCIIString,JuliaFEM.Field} with 4 entries:\n",
" \"geometry\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
" \"heat coefficient\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
" \"displacement\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
" \"temperature\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Field, TimeStep, Increment\n",
"element = MyQuad4([1, 2, 3, 4])\n",
"\n",
"# create one timestep for t=0\n",
"# is equivalent to \n",
"# element[\"geometry\"] = (0.0 => Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])\n",
"element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
"\n",
"# create two timesteps, one for t=0.0 and another for t=1.0\n",
"element[\"temperature\"] = (0.0, [0.0, 0.0, 0.0, 0.0]), (1.0, [1.0, 2.0, 3.0, 4.0])\n",
"\n",
"# create two timesteps, one is zeros at time t=0.0 and another with t=1.0\n",
"element[\"displacement\"] = (\n",
" 0.0 => Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]],\n",
" 1.0 => Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]])\n",
"\n",
"# create two timesteps, one for t=0.0 and another for t=1.0\n",
"element[\"heat coefficient\"] = (0.0, 2), (1.0, 3)\n",
"element.fields"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.5\n",
" 0.5"
]
},
"execution_count": 11,
"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": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1x2 Array{Float64,2}:\n",
" 0.0 1.0"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolate derivatives works too\n",
"dbasis = grad(basis)\n",
"# temperature gradient at mid point of element at time t=0.5\n",
"dbasis(\"temperature\", [0.0, 0.0], 0.5)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.5"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolating scalar -> scalar.\n",
"basis(\"heat coefficient\", [0.0, 0.0], 0.5)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 14,
"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": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x2 Array{Float64,2}:\n",
" 0.132813 0.0703125\n",
" 0.0703125 0.0078125"
]
},
"execution_count": 15,
"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": 16,
"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": 16,
"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",
" basis = get_basis(element)\n",
" dbasis = grad(basis) # gradient is with respect to \"geometry\" field\n",
" detJ = det(basis)\n",
" M += w*rho*basis(xi)'*basis(xi)*detJ(xi) # mass matrix\n",
" K += w*k*dbasis(xi)'*dbasis(xi)*detJ(xi) # stiffness matrix\n",
" f += w*basis(xi)'*q*detJ(xi) # 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": 17,
"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": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"size (generic function with 83 methods)"
]
},
"execution_count": 18,
"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\"] = []\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\"] = []\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": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_force_vector (generic function with 3 methods)"
]
},
"execution_count": 19,
"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",
" rho = basis(\"density\", ip, time)\n",
" return rho * 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\", ip, time)\n",
" return k * 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\", ip, time)\n",
" return f*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\", ip, time)\n",
" return g*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": 20,
"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",
"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",
"\n",
"# create boundary element with \n",
"boundary_element = Seg2([1, 2])\n",
"boundary_element[\"geometry\"] = 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\"] = (0.0 => 0.0, 1.0 => 6.0)\n",
"boundary_equation = DC2D2(boundary_element);"
]
},
{
"cell_type": "code",
"execution_count": 21,
"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": 21,
"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": 22,
"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": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly.stiffness_matrix"
]
},
{
"cell_type": "code",
"execution_count": 23,
"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": 23,
"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": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 24,
"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": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 1.0\n",
" 1.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"local_assembly = initialize_local_assembly()\n",
"calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\", 1.0)\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=0$. 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": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.5\n",
" 0.5"
]
},
"execution_count": 26,
"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": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"size (generic function with 84 methods)"
]
},
"execution_count": 27,
"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",
" # Initial configuration needs to be defined if using autodiff\n",
" element[\"temperature\"] = zeros(4)\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": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_potential_energy (generic function with 2 methods)"
]
},
"execution_count": 28,
"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",
" c = basis(\"temperature nonlinearity coefficient\", ip, time)\n",
" gradT = grad(basis)(\"temperature\", ip, time, variation)\n",
" Wint = 1/2 * k*(1 + c*T) * vecdot(gradT, gradT)\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": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"increment 1, |du| = 1.4142135624\n",
"increment 2, |du| = 0.4242640687\n",
"increment 3, |du| = 0.0465655685\n",
"increment 4, |du| = 0.0005747960\n",
"increment 5, |du| = 0.0000000876\n",
"increment 6, |du| = 0.0000000000\n",
"increment 7, |du| = 0.0000000000\n",
"elapsed time: 1.142529838 seconds\n",
"error: 0.0\n"
]
}
],
"source": [
"function run_simulation()\n",
" # create model -- start\n",
" element = Quad4([1, 2, 3, 4])\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\"] = [0.0, 0.0, 0.0, 0.0]\n",
" element[\"temperature nodal load\"] = [3.0, 3.0, 0.0, 0.0]\n",
" element[\"temperature nonlinearity coefficient\"] = [1.0, 1.0, 1.0, 1.0]\n",
" equation = DC2D4NL(element)\n",
" # create model -- end\n",
"\n",
" la = initialize_local_assembly() # create workspace for local matrices\n",
" T = zeros(4) # create workspace for solution vector\n",
" dT = zeros(4) # \n",
" fd = [1, 2] # free dofs\n",
" tic()\n",
" # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n",
" for i=1:10\n",
" calculate_local_assembly!(la, equation, \"temperature\") # calculate local matrices\n",
" dT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd] # <-- note sign convention, more on this below\n",
" T += dT\n",
" push!(element[\"temperature\"], T) # add new increment to model\n",
" # print some convergence information\n",
" @printf(\"increment %2d, |du| = %12.10f\\n\", i, norm(dT))\n",
" err = last(element[\"temperature\"])[1] - 2/3\n",
" isapprox(err, 0.0) && break\n",
" end\n",
" toc()\n",
" err = last(element[\"temperature\"])[1] - 2/3\n",
" println(\"error: $err\")\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",
"\\Rightarrow \\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_{\\Omega_{0}}\\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": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_residual_vector (generic function with 3 methods)"
]
},
"execution_count": 30,
"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",
" # initial field must be defined if using autodiff\n",
" element[\"displacement\"] = zeros(2, 4)\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",
" 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": [
"We have 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": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"displacement at tip: -8.773031198197748\n"
]
}
],
"source": [
"using JuliaFEM.Test\n",
"using JuliaFEM: get_integration_points, solve!\n",
"\n",
"function run_simulation_2()\n",
" # create model -- start\n",
" element = Quad4([1, 2, 3, 4])\n",
" element[\"geometry\"] = Vector[[0.0,0.0], [10.0,0.0], [10.0,1.0], [0.0,1.0]]\n",
" element[\"youngs modulus\"] = 500.0\n",
" element[\"poissons ratio\"] = 0.3\n",
" element[\"displacement volume load\"] = 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",
" @test isapprox(disp, -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": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: solution = [1.0000000000000002,1.0000000000000002]\n"
]
}
],
"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=0.0,\n",
" problem=nothing)\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",
" solution = la.stiffness_matrix[free_dofs, free_dofs] \\ la.force_vector[free_dofs]\n",
" info(\"solution = $solution\")\n",
" @test isapprox(solution, [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. More info about 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 \"equations\" can be set to one element, so it's totally possible for example solve displacement and temperature equations using 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": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"PlaneHeatProblem"
]
},
"execution_count": 33,
"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": 34,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"call (generic function with 1271 methods)"
]
},
"execution_count": 34,
"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=0.0)\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)) # 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 = full(x[1:length(b1)])\n",
" x2 = full(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",
" local_solution = reshape(x1[gdofs], size(equation))\n",
" push!(element[field_name], local_solution)\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",
" local_solution = reshape(x2[gdofs], size(equation))\n",
" push!(element[field_name], local_solution)\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": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: assembling problem for temperature\n",
"INFO: dimension of unknown field: 1, problem dofs: 4\n",
"INFO: assembling problem for reaction force\n",
"INFO: dimension of unknown field: 1, problem dofs: 4\n",
"INFO: Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n"
]
}
],
"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",
" 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",
" el2[\"temperature flux\"] = (0.0 => 0.0, 1.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",
" # 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)[1]\n",
" info(\"Temperature at point X = $X is T = $T\")\n",
" @test isapprox(T, 100.0)\n",
"end\n",
"\n",
"main()"
]
},
{
"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"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_potential_energy (generic function with 3 methods)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\" 2-node radiation boundary element. \"\"\"\n",
"type DC2D2NL <: Heat\n",
" element :: Seg2\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
"end\n",
"\n",
"function DC2D2NL(element::Seg2)\n",
" # https://pomax.github.io/bezierinfo/legendre-gauss.html\n",
" integration_points = [\n",
" IntegrationPoint([0.0], 128/225),\n",
" IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),\n",
" IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),\n",
" IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),\n",
" IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900)\n",
" ]\n",
" if !haskey(element, \"temperature\")\n",
" element[\"temperature\"] = zeros(2)\n",
" end\n",
" DC2D2NL(element, integration_points)\n",
"end\n",
"\n",
"function Base.size(equation::DC2D2NL)\n",
" return (1, 2)\n",
"end\n",
"\n",
"function JuliaFEM.get_potential_energy(equation::DC2D2NL, ip, time; variation=nothing)\n",
" element = get_element(equation)\n",
" basis = get_basis(element)\n",
" T = basis(\"temperature\", ip, time, variation)[1]\n",
" T_ext = basis(\"temperature external\", ip, time)[1]\n",
" coeff = basis(\"temperature coefficient\", ip, time)[1]\n",
" q0 = coeff*(T_ext^4 - T^4)\n",
" Wint = 0.0\n",
" Wext = q0*T\n",
" return Wint - Wext\n",
"end\n",
"\n",
"function JuliaFEM.has_potential_energy(equation::DC2D2NL)\n",
" return true\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Defining new equation mapping to old problem is one line command. Here's the equation `DC2D2` is replaced to $\\rightarrow$ `DC2DCNL`"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"INFO: residual norm: 0.7071067811865476\n",
"INFO: residual norm: 2.2097084800638154e-9\n",
"INFO: residual norm: 2.09345658541016e-16\n",
"INFO: Temperature at the midpoint of boundary element: 0.4999999984375\n"
]
}
],
"source": [
"function run_radiation_model()\n",
"\n",
" # create elements\n",
" N = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
" element1 = Quad4([1, 2, 3, 4])\n",
" element1[\"geometry\"] = Vector[N[1], N[2], N[3], N[4]]\n",
" element1[\"temperature thermal conductivity\"] = 6.0\n",
" element1[\"temperature load\"] = [0.0, 0.0, 0.0, 0.0]\n",
" element1[\"density\"] = 36.0\n",
" element1[\"temperature nonlinearity coefficient\"] = [0.0, 0.0, 0.0, 0.0]\n",
" element2 = Seg2([1, 2])\n",
" element2[\"geometry\"] = Vector[N[1], N[2]]\n",
" #element2[\"emissivity\"] = 0.5\n",
" element2[\"temperature coefficient\"] = 3.0e-8 # ~ 5.7e-8 * 0.5\n",
" element2[\"temperature external\"] = 100.0\n",
"\n",
" # set initial conditions\n",
" element1[\"temperature\"] = zeros(4)\n",
" element2[\"temperature\"] = zeros(2)\n",
" # create equations\n",
"\n",
" equation1 = DC2D4NL(element1)\n",
" equation2 = DC2D2NL(element2)\n",
" #equations = [equation1, equation2]\n",
"\n",
" # create workspace for local matrices\n",
" la1 = initialize_local_assembly(equation1)\n",
" la2 = initialize_local_assembly(equation2)\n",
"\n",
" #solve!(problem, free_dofs; max_iterations=5, dump_matrices=true)\n",
"\n",
" T = zeros(4)\n",
" dT = zeros(4)\n",
" free_dofs = [1, 2]\n",
" K = zeros(4,4)\n",
" f = zeros(4)\n",
" for i=1:20\n",
" fill!(K, 0.0)\n",
" fill!(f, 0.0)\n",
" calculate_local_assembly!(la1, equation1, \"temperature\")\n",
" calculate_local_assembly!(la2, equation2, \"temperature\")\n",
"\n",
" K += la1.stiffness_matrix\n",
" f += la1.force_vector\n",
" #info(\"f(linear) = \\n$(la1.force_vector)\")\n",
" #info(\"K(linear) = \\n$(la1.stiffness_matrix[free_dofs, free_dofs])\")\n",
"\n",
" K[free_dofs, free_dofs] += la2.stiffness_matrix\n",
" f[free_dofs] += la2.force_vector\n",
" #info(\"K(nonlinear) = \\n$(la2.stiffness_matrix)\")\n",
" #info(\"f(nonlinear) = \\n$(la2.force_vector)\")\n",
"\n",
" dT[free_dofs] = K[free_dofs, free_dofs] \\ f[free_dofs]\n",
" T += dT\n",
"\n",
" push!(element1[\"temperature\"], T)\n",
" push!(element2[\"temperature\"], T[free_dofs])\n",
" T_mid = get_basis(element2)(\"temperature\", [0.0])[1]\n",
" info(\"residual norm: $(norm(dT))\")\n",
" norm(dT) < 1.0e-12 && break\n",
" end\n",
" #solve!(problem, free_dofs; max_iterations=50)\n",
"\n",
" T_mid = get_basis(element2)(\"temperature\", [0.0])[1]\n",
" info(\"Temperature at the midpoint of boundary element: $T_mid\")\n",
"end\n",
"\n",
"run_radiation_model()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Accessing integration points\n",
"\n",
"It's possible to save internal variables to integration points. This is described [constitutive modelling tutorial](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb).\n",
"\n",
"### Statistical fields\n",
"\n",
"It has been already discussed in [data structures](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-06-14-data-structures.ipynb) how to make fields dependent from other fields. Another idea is to make field dependent from external resources (data from logging system / internet / random.org). For example, sampling data from normal distribution:"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"first (generic function with 13 methods)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"type RandomField <: JuliaFEM.DiscreteField\n",
" mu :: Float64\n",
" std :: Float64\n",
"end\n",
"Base.first(field::RandomField) = Increment(randn(4).*field.std^2 + field.mu)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4-element JuliaFEM.Increment{Float64}:\n",
" 9.18727\n",
" 10.2647 \n",
" 7.93236\n",
" 12.0819 "
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"RandomField(10.0, 2.0)(0.0)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"function calc_stochastic(N=500)\n",
" solutions = Vector[]\n",
" element = Quad4([1, 2, 3, 4])\n",
" element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
" element[\"temperature load\"] = [12.0, 12.0, 12.0, 12.0]\n",
" element[\"density\"] = 36.0\n",
" element[\"temperature thermal conductivity\"] = RandomField(6.0, 1.0)\n",
" equation = DC2D4(element)\n",
" la = initialize_local_assembly()\n",
" fd = [1,2]\n",
" for i=1:N\n",
" calculate_local_assembly!(la, equation, \"temperature\", 0.0)\n",
" push!(solutions, la.stiffness_matrix[fd, fd] \\ la.force_vector[fd])\n",
" end\n",
" return solutions\n",
"end\n",
"solutions = calc_stochastic();"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.997586\n",
" 0.998368"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mean(solutions)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"using PyPlot"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAADUCAYAAAAcN2ODAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X1UVHX+B/D3DMIwIpIMIIMco9UfiS7Es2lmGavhimXrilH+2FTQNRIXc43M8mxsWCoLu9LxKX/ppsuKT7s+lG2amlYrKGPH3wbJnvUJHBVEEXF4/v7+8DeT08A4I3fmDvl+ncPJvvd7733PJ/XTvXMfFEIIASIiIuqUUu4AREREroyNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKNkoiIiIr2CiJiIisYKMkIiKygo2SiIjICpsbZWNjI5YsWYLExET4+vpCqVRi48aNVtdpbW3F0KFDoVQqkZeX1+mc9evXIywsDGq1GqGhoSgsLLTvExARETmQzY2ypqYGOTk5+O677xAZGQkAUCgUVtdZuXIlLly40OXcNWvWID09HeHh4SgsLMSIESOQmZmJZcuW2fMZiIiIHMbmRhkUFIRLly7hzJkzWL58+V3nX7lyBTk5OcjOzu50ucFgwBtvvIGkpCQUFxdj5syZ2LhxI1588UXk5OTg+vXrtn8KIiIiB7G5UXp4eCAgIAAAYMsLR7KzszFkyBC8+OKLnS4/ePAg6urq8PLLL5uNZ2RkoLGxEXv37rU1GhERkcM45GKekpIS/PnPf0ZBQUGXc3Q6HQAgNjbWbDw6OhpKpRInT550RDQiIiK7SN4ohRCYO3cunn/+eQwfPrzLeXq9Hm5ubvDz8zMb9/DwgEajwcWLF6WORkREZLdeUm9ww4YN+N///V/s2LHD6jyDwQAPD49Ol6lUKhgMBovx2tpafPrppwgJCYFarZYkLxER9TwGgwFnz57F008/bXHAJTVJG+WNGzfw+uuvY+HChRgwYIDVuWq1Gi0tLZ0ua2pq6rQRfvrpp5g2bZokWYmIqOfbtGlTl9fCSEXSRrlixQq0trYiOTkZZ8+eBQBUVVUBAOrq6nD27FkMGDAA7u7u0Gq1aG9vR21trdn/DbS0tKCurg5BQUEW2w8JCQFwuzBhYWFSRneIrKws5Ofnyx3DJj0pK9Cz8jKrYzCrY/SUrOXl5Zg2bZqpLziSpI3ywoULuHbtGoYNG2axLDc3F7m5uTh58iQiIiJM92KWlpZi/PjxpnnHjx9HR0eHafmdjEeZYWFhiI6OljK6Q/j4+PSInEDPygr0rLzM6hjM6hg9KSsAp3wNJ2mjzMzMxHPPPWc2dvnyZcyePRvTp0/Hs88+a+r+CQkJ8PX1xapVq8wa5apVq+Dl5YUJEyZIGY2IiOie2NUoCwsLcf36ddMVqbt27cL58+cB3G6SUVFRiIqKMlvHeAp22LBheOaZZ0zjnp6eyMnJQUZGBpKTkzFu3DgcOXIEmzdvRm5uLh544IHufC4iIiJJ2NUo8/LycO7cOQC3H0m3c+dO7NixAwqFAqmpqejbt69dO58zZw7c3d2Rl5eHXbt2YeDAgSgoKEBmZqZd2yGS0q1bt1BRUSHZ9urr61FWVtatbQwZMgS9e/eWKBER2cOuRnnmzBm7dxASEoKOjo4ul6elpSEtLc3u7fYEKSkpckewWU/KCjg2b0VFBWJiYiTdZne3d+LECad8b9STfh8wq2P0pKzOohC2PI/ORZSVlSEmJsZpf2nQ/cn4+wzYBEDuq6vLAUzj73miH3BmP5D8gQNEPx5hANiciO53fHEzERGRFWyUREREVrBREhERWcFGSUREZAUbJRERkRVslERERFawURIREVlhc6NsbGzEkiVLkJiYCF9fXyiVSmzcuNFsjhACGzZswDPPPIOBAweiT58+CA8PxzvvvIPm5uZOt7t+/XqEhYVBrVYjNDQUhYWF3ftEREREErK5UdbU1CAnJwffffed6RVYCoXCbE5jYyNmzJiBq1evYs6cOfjjH/+I+Ph4LFmyxOwNIUZr1qxBeno6wsPDUVhYiBEjRiAzMxPLli3r5sciIiKShs1P5gkKCsKlS5cQEBCAEydOIC4uzmKOSqXCV199hUcffdQ0NnPmTISEhGDJkiU4cOAAEhISAAAGgwFvvPEGkpKSUFxcbJrb0dGBnJwczJo1i28QISIi2dl8ROnh4YGAgAAAt0+xdsbd3d2sSRpNmjQJAMzeyHDw4EHU1dXh5ZdfNpubkZGBxsZG7N2719ZoREREDuOUi3kuXboEAPDz8zON6XQ6AEBsbKzZ3OjoaCiVSpw8edIZ0YiIiKxySqNctmwZfHx8zL6n1Ov1cHNzM2uewO0jV41GY3o5NBERkZwc/vaQ3NxcHDhwAKtWrTJ7sbPBYICHh0en66hUKhgMBkdHIyIiuiuHNsotW7bgzTffRFpaGmbPnm22TK1Wo6WlpdP1mpqaoFaru9xuVlYWfHx8zMZSUlL4wlEioh+hoqIiFBUVmY3V19c7bf8Oa5SfffYZUlNTkZSUhNWrV1ss12q1aG9vR21trdnp15aWFtTV1SEoKKjLbefn5/MltkRE94nODoS+f8G64znkO8pjx47hueeeQ3x8PIqLi6FUWu4mKioKAFBaWmo2fvz4cXR0dJju1SQiIpKT5I2yvLwcEyZMwE9+8hPs2bMHKpWq03lPPfUUfH19sWrVKrPxVatWwcvLCxMmTJA6GhERkd3sOvVaWFiI69evm65I3bVrF86fPw8AyMzMhEKhwNNPP43r169j4cKF2L17t9n6gwcPNt1n6enpiZycHGRkZCA5ORnjxo3DkSNHsHnzZuTm5vJhA0RE5BLsapR5eXk4d+4cgNuPr9u5cyd27NgBhUKB1NRUdHR0oKqqCgqFAtnZ2Rbrv/TSS2YPJJgzZw7c3d2Rl5eHXbt2YeDAgSgoKEBmZmY3PxYREZE07GqUZ86cueucjo4OuwKkpaUhLS3NrnXox+fWrVtmT26SU3l5udwRiMiFOPw+SiJbVFRUOO0KNiIie7BRkovZBCBM5gwfA3hT5gxE5CrYKMnFhAGQ+x5Znnolou855VmvREREPRUbJRERkRVslERERFawURIREVnBRklERGQFGyUREZEVNjfKxsZGLFmyBImJifD19YVSqcTGjRs7nVteXo7ExER4e3tDo9EgNTUVtbW1nc5dv349wsLCoFarERoaisLCwnv7JERERA5gc6OsqalBTk4OvvvuO9MrsBQKhcW8qqoqjB49Gv/5z3+wdOlSLFiwAHv37sXYsWPR2tpqNnfNmjVIT09HeHg4CgsLMWLECGRmZmLZsmXd/FhERETSsPmBA0FBQbh06RICAgJw4sQJxMXFdTovNzcXBoMBOp0OwcHBAID4+HiMHTsWGzZsQHp6OgDAYDDgjTfeQFJSEoqLiwEAM2fOREdHB3JycjBr1iy+QYSIiGRn8xGlh4cHAgICAABCiC7nbd++HUlJSaYmCQAJCQkIDQ01NUQAOHjwIOrq6vDyyy+brZ+RkYHGxkbs3bvX5g9BRETkKJJezFNdXY2amhrExsZaLIuLi4NOpzP9u/HXP5wbHR0NpVKJkydPShmNiIjonkjaKPV6PQBAq9VaLNNqtairqzN9T6nX6+Hm5gY/Pz+zeR4eHtBoNKaXQxMREclJ0kZpMBgAACqVymKZp6en2RyDwQAPD49Ot6NSqUzziIiI5CTp20PUajUAoLm52WJZU1OT2Ry1Wo2WlpZOt9PU1GSa15msrCz4+PiYjaWkpCAlJeWechMRkesqKipCUVGR2Vh9fb3T9i9pozSecjWegr2TXq+HRqOBu7u7aW57eztqa2vNTr+2tLSgrq4OQUFBXe4nPz8f0dFyv4qJiIicobMDobKyMqe97F3SU68DBgyAv78/SktLLZaVlJSY7r8EgKioKACwmHv8+HF0dHSYzSUiIpKL5I+wmzx5Mvbs2YOqqirT2IEDB1BZWYkpU6aYxp566in4+vpi1apVZuuvWrUKXl5emDBhgtTRiIiI7GbXqdfCwkJcv37ddEXqrl27cP78eQBAZmYm+vbti0WLFmHr1q0YM2YM5s2bh4aGBixfvhwRERGYPn26aVuenp7IyclBRkYGkpOTMW7cOBw5cgSbN29Gbm4uHzZAREQuwa5GmZeXh3PnzgG4/fi6nTt3YseOHVAoFEhNTUXfvn0RHByMw4cPY/78+cjOzoZKpUJSUhLy8vJM308azZkzB+7u7sjLy8OuXbswcOBAFBQUIDMzU7pPSERE1A12NcozZ87YNG/o0KHYt2+fTXPT0tKQlpZmTwwiIiKn4Wu2iIiIrGCjJCIisoKNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKNkoiIiIr2CiJiIiskPQ1W0bHjx9HTk4OSktLUV9fj4EDB+KFF17AggULzN4zWV5ejqysLHz55Zfw8PDAhAkT8Ic//MHstVtE97fbLzAvLy+XOYe5IUOGoHfv3nLHIHIKyRvlqVOnMGrUKAQFBeE3v/kNfH198dVXX2HJkiU4ceIE/va3vwEAqqqqMHr0aPTr1w9Lly5FQ0MDVqxYgVOnTqGkpMTiubBE96ezAIBp06bJG+MHTpw4wXfC0n1D8kb517/+FS0tLdi7dy/CwsIA3H6ea0dHB/785z+jvr4ePj4+yM3NhcFggE6nQ3BwMAAgPj4eY8eOxYYNG5Ceni51NKIebBOAMLlDACgH4FpNm8jRJG+UxlOrAQEBZuOBgYFwc3ODh4cHAGD79u1ISkoyNUkASEhIQGhoKIqLi9koicyEAeARHJEcJL+YZ8aMGejfvz9mzpyJb775BhcuXMCWLVuwevVqZGZmQq1Wo7q6GjU1NYiNjbVYPy4uDjqdTupYRERE90TyI8qgoCB8+eWX+PnPf46oqCjT+OLFi/H2228DAPR6PQBAq9VarK/ValFXV4fW1lZ+T0lERLKTvFFevnwZ48ePBwCsW7cOGo0Ge/bswTvvvIP+/fsjIyMDBsPtK/lUKpXF+p6engAAg8HARklERLKTvFHm5OSguroap0+fRlBQEABg0qRJ6OjowGuvvYaUlBTT95jNzc0W6zc1NQGA2W0kP5SVlQUfHx+zsZSUFKSkpEj1MYiIyEUUFRWhqKjIbKy+vt5p+5e8UR49ehRRUVGmJmk0ceJEbNiwASdPnsTDDz8M4PtTsHfS6/XQaDRWjybz8/N5aToR0X2iswOhsrIyxMTEOGX/kl/M09raivb29k7HAaCtrQ0DBgyAv78/SktLLeaVlJQgMjJS6lhERET3RPJGGR0djbKyMlRWVpqNFxUVwc3NDREREQCAyZMnY8+ePaiqqjLNOXDgACorKzFlyhSpYxEREd0TyU+9/va3v8X27dvx+OOP45VXXoGvry/27NmDffv2IT09HYGBgQCARYsWYevWrRgzZgzmzZuHhoYGLF++HBEREZg+fbrUsYiIiO6J5I0yIiIChw4dwpIlS7B8+XI0NTXhJz/5CXJzc7Fw4ULTvODgYBw+fBjz589HdnY2VCoVkpKSkJeXx6tdneDWrVuoqKiQO4aJqz3LlIjIyCEPRY+Pj8cnn3xy13lDhw7Fvn37HBGB7qKiosJpX4QTEfVkDmmU1JO4yjNEPwbwptwhiIgssFHe91zlGaI89UpErokvbiYiIrKCjZKIiMgKNkoiIiIr2CiJiIisYKMkIiKygo2SiIjICoc1yrKyMjzzzDPQaDTw8vJCeHg4Vq5caTanvLwciYmJ8Pb2hkajQWpqKmprax0ViYiIyG4OuY/yH//4ByZOnIiYmBi89dZb6NOnD/7973+jurraNKeqqgqjR49Gv379sHTpUjQ0NGDFihU4deoUSkpK+Bg7IiJyCZI3yhs3biA1NRUTJ07Etm3bupyXm5sLg8EAnU6H4OBgALcffTd27Fhs2LAB6enpUkcjIiKym+SnXv/yl7/gypUreOeddwAAjY2N6OjosJi3fft2JCUlmZokACQkJCA0NBTFxcVSxyIiIronkjfK/fv3o2/fvrhw4QIefvhheHt7w8fHBy+//DKam5sBANXV1aipqUFsbKzF+nFxcdDpdFLHIiIiuieSN8rKykq0tbVh0qRJGD9+PHbs2IEZM2Zg9erVpvdM6vV6AIBWq7VYX6vVoq6uDq2trVJHIyIispvk31HevHkTt27dwpw5c1BQUAAAmDRpElpaWrBmzRq8/fbbMBgMAACVSmWxvqenJwDAYDDwgh4iIpKd5EeUarUaAJCSkmI2bvz3f/7zn6Y5xlOxd2pqajLbDhERkZwkP6IMCgrCt99+i/79+5uNBwQEAACuXbuGoKAgAN+fgr2TXq+HRqOxejSZlZUFHx8fs7GUlBSL5kxERD1fUVERioqKzMbq6+udtn/JG2VsbCz279+Pqqoq/Nd//Zdp/OLFiwAAf39/BAUFwd/fH6WlpRbrl5SUIDIy0uo+8vPzER3tCu9QJCIiR+vsQKisrAwxMTFO2b/kp16Tk5MBAOvXrzcb/+CDD+Du7o4nn3wSADB58mTs2bMHVVVVpjkHDhxAZWUlpkyZInUsIiKieyL5EWVkZCRmzJiB//mf/0FbWxtGjx6NQ4cOYdu2bVi0aBECAwMBAIsWLcLWrVsxZswYzJs3Dw0NDVi+fDkiIiJMV8cSERHJzSGPsFu9ejUGDhyIDz/8EDt37kRISAgKCgqQmZlpmhMcHIzDhw9j/vz5yM7OhkqlQlJSEvLy8ni1KxERuQyHNMpevXrhrbfewltvvWV13tChQ7Fv3z5HRCAiIpIEX7NFRERkhUOOKInox+r2w0LKy8tlzvG9IUOGoHfv3nLHoB8xNkoissNZAMC0adPkjXGHEydO8HYxcig2SiK6B5sAhMmcoRyA6zRs+vFioySiexAGgEdxdH/gxTxERERWsFESERFZwUZJRERkBRslERGRFU5plO+88w6USiXCw8MtlpWXlyMxMRHe3t7QaDRITU1FbW2tM2IRERHdlcOveq2qqkJubi68vLygUCgslo0ePRr9+vXD0qVL0dDQgBUrVuDUqVMoKSnhM1+JiEh2Dm+UCxYswMiRI9HW1mZxpJibmwuDwQCdTofg4GAAQHx8PMaOHYsNGzYgPT3d0fGIiIiscuip1y+++ALbt29HQUEBhBAWR5Tbt29HUlKSqUkCQEJCAkJDQ1FcXOzIaERERDZxWKNsb2/H3LlzkZ6ejmHDhlksr66uRk1NDWJjYy2WxcXFQafTOSoaERGRzRx26nX16tU4f/48Pv/8806X6/V6AIBWq7VYptVqUVdXh9bWVn5PSUREsnLIEeXVq1dN76PUaDSdzjEYbr+FQKVSWSzz9PQ0m0NERCQXhzTKxYsXw8/PD3Pnzu1yjlqtBgA0NzdbLGtqajKbQ0REJBfJT71WVlZi3bp1KCgoQFVVlWm8qakJLS0tOHfuHPr27Ws65Wo8BXsnvV4PjUbT5WnXrKws+Pj4mI2lpKQgJSVFwk9CRESuoKioCEVFRWZj9fX1Ttu/5I2yuroaHR0dyMzMRGZmpsXyhx56CL/5zW/whz/8Af7+/igtLbWYU1JSgsjIyC73kZ+fz/fPERHdJzo7ECorK0NMTIxT9i95owwPD8fOnTvNbgURQmDx4sW4efMm/vjHP2LQoEEAgMmTJ2Pjxo2oqqoy3SJy4MABVFZW4tVXX5U6GhERkd0kb5QajQbPPvusxXh+fj4A4JlnnjGNLVq0CFu3bsWYMWMwb948NDQ0YPny5YiIiMD06dOljkZERGQ3pz0UXaFQWDxwIDg4GIcPH8agQYOQnZ2NFStWICkpCZ999hlvCyEiIpfg8EfYGR08eLDT8aFDh2Lfvn3OikFERGQXvmaLiIjICjZKIiIiK9goiYiIrGCjJCIisoKNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKyRtlaWkpXnnlFQwbNgx9+vTBgw8+iKlTp6KystJibnl5ORITE+Ht7Q2NRoPU1FTU1tZKHYmIiOieSf4Iu/feew9ff/01pkyZgoiICOj1ehQWFiI6Ohr//Oc/MWzYMABAVVUVRo8ejX79+mHp0qVoaGjAihUrcOrUKZSUlPBZr0RE5BIkb5Svvvoq4uLi0KvX95ueOnUqwsPD8e677+Kjjz4CAOTm5sJgMECn05lesRUfH4+xY8diw4YNSE9PlzoaERGR3SQ/9TpixAizJgkAgwcPxtChQ1FRUWEa2759O5KSkkxNEgASEhIQGhqK4uJiqWMRERHdE6dczCOEwOXLl+Hn5wcAqK6uRk1NDWJjYy3mxsXFQafTOSMWERHRXTmlUW7evBkXL17E1KlTAQB6vR4AoNVqLeZqtVrU1dWhtbXVGdGIiIiscnijrKioQEZGBkaOHIlf/epXAACDwQAAUKlUFvM9PT3N5hAREcnJoY3y0qVLmDBhAvr164dt27ZBoVAAANRqNQCgubnZYp2mpiazOURERHKS/KpXo/r6eowfPx43btzAkSNHEBgYaFpmPOVqPAV7J71eD41GY/X2kKysLPj4+JiNpaSkICUlRaL0RETkKoqKilBUVGQ2Vl9f77T9O6RRNjU1YeLEifj3v/+N/fv3Y8iQIWbLBwwYAH9/f5SWllqsW1JSgsjISKvbz8/PR3R0tKSZiYjINXV2IFRWVoaYmBin7F/yU6/t7e2YOnUqjh07hq1bt2L48OGdzps8eTL27NmDqqoq09iBAwdQWVmJKVOmSB2LiIjonjjkgQO7d+/GxIkTUVtbi02bNpktnzZtGgBg0aJF2Lp1K8aMGYN58+ahoaEBy5cvR0REBKZPny51LCIionsieaP85ptvoFAosHv3buzevdtsmUKhMDXK4OBgHD58GPPnz0d2djZUKhWSkpKQl5fHx9cRkQ1uXxlfXl4ucw5zQ4YMQe/eveWOQRKSvFEePHjQ5rlDhw7Fvn37pI5ARPeFswC+P0vlKk6cOMFrKH5kHHbVKxGRc2wCECZ3CADlAFyraZM02CiJqIcLA8AjOHIcvriZiIjICjZKIiIiK3jqlYhIEq53FS6vwJUGG6UT3bp1y+ydnHJypT/MRD8OZwG41lW4vAJXGmyUTlRRUeG0Ry4RkVxc4SpcXoErJTZKWbjCH6SPAbwpcwaiHyNehftjw0YpC1f4g8RTr0REtpD1qtfm5ma89tprCAoKQu/evfHoo49i//79ckYiIiIyI2ujfOmll5Cfn4///u//xp/+9Ce4ubnh5z//Ob788ks5Y0nmh+9Pc209KSvQs/Iyq2N8JXcAO/Scuvasv7ecQ7ZGWVJSgi1btuDdd9/Fe++9h7S0NHz++ed48MEHsXDhQrliSapn/YbrSVmBnpWXWR3ja7kD2KHn1LVn/b3lHLI1ym3btqFXr16YNWuWaUylUmHmzJn4+uuvUV1dLVc0IiIiE9ku5tHpdAgNDUWfPn3MxuPi4gAAJ0+exIABA7q1j7///e/46iv5Ts98++23eO2110z/funSJdmyENH95N4fflBfX4+ysjKpAwHouQ9AkK1R6vV6aLVai3Hj2MWLF7u9j/379+P999eiV68Hu72te9HWdgn5+TuN/4bW1jOy5CCi+81ZAPf+8ANH3e/dUx+AIFujNBgMUKlUFuOenp6m5Z2tA9j+f0lXrlyBUhmA1tYF3UjaHe+jtTXj/399HYDx6PJjyH97hvGCKWOWKgCb5YtjkeduHJnX3ix3052sUme5m7tldXYea+r+/5+ukAWwXhtn//kyZpkJwPKAxLotAKZKGwd6AOslfSKYcVud9QrJCZkMGzZM/OxnP7MY/9e//iUUCoVYu3atxbJNmzYJAPzhD3/4wx/+CABi06ZNDu9Xsh1RarXaTk+v6vV6AEBQUJDFsqeffhqbNm1CSEgI1Gq1wzMSEZFrMhgMOHv2LJ5++mmH70u2RhkVFYVDhw6hoaEB3t7epvFjx44BACIjIy3W8fPzw4svvui0jERE5Loee+wxp+xHtttDfvnLX6K9vR1r1641jTU3N+PDDz/Eo48+2u0rXomIiKQg2xFlfHw8pkyZgtdffx1XrlzBoEGDsHHjRpw/fx4ffvihXLGIiIjMKIQQQq6dNzc3480338SmTZtw7do1PPLII8jJycHYsWPlikRERGRG1kZJRETk6mR9KDoREZGrc3qj7M6rtfbv34+EhAQEBATA29sbjzzyCFauXImOjg6LuV999RVGjRoFLy8vaLVazJs3D42NjS6X9cknn4RSqbT4GT9+vNOyfvbZZ6Za+fr6YsqUKTh37lync+Wuq61ZpaprY2MjlixZgsTERPj6+kKpVGLjxo02r3/9+nXMmjUL/v7+6NOnD5566inodLpO53a3ts7KKkVtu5P10qVLyM7OxpgxY+Dt7Q2lUonDhw93OV/OutqTVe66HjhwADNmzEBoaCi8vLwwaNAgpKend/noTTnrak9WSf4ucPidmj/w/PPPC3d3d7Fw4UKxbt06MXLkSOHu7i6OHj1qdb1PPvlEKBQKER4eLgoKCsTatWvFpEmThEKhEPPmzTObq9PphKenp4iJiRFr1qwRixcvFp6enmL8+PEul/WJJ54QAwcOFJs3bzb7OXjwoFOy7t69WyiVShEfHy9Wrlwpfv/73wt/f38RHBwsampqzObKXVd7skpV1zNnzgiFQiFCQkLEmDFjhEKhEBs3brRp3fb2djFy5EjRp08f8fbbb4v3339fDBs2TPTt21dUVlaazZWits7KKkVtu5P14MGDQqFQiIcffliMHDlSKBQKcfjw4U7nyl1Xe7LKXdeYmBgxaNAgkZ2dLdavXy8WLVok+vbtKwIDA8WlS5fM5spdV3uySlFXpzbKY8eOCYVCIfLy8kxjTU1NYvDgwWLkyJFW133hhReEp6enuHbtmtn4E088IXx8fMzGxo8fLwYMGCAaGhpMYx988IFQKBTiH//4h0tlfeKJJ0R4eLhNmRyRdejQoSI0NFS0traaxr755hvh5uYmXn31VbO5ctfVnqxS1FUIIZqbm8Xly5eFEEIcP37crj/MW7ZsEQqFQmzfvt00VlNTI/r16yfnUErKAAAIPklEQVReeOEFs7lS1NZZWaWobXeyNjQ0mP5sbd261Wrzkbuu9mSVu65HjhyxGPviiy+EQqEQixcvNhuXu672ZJWirk499dqdV2up1WqoVCr4+PiYjQcGBpo9jf7GjRvYv38/pk2bZvZmktTUVPTp0wfFxcUuk9VICIH29nbcvHnTpmxSZa2rq0N5eTmee+459Or1/Z1CERERGDJkCP7617+axuSuqz1ZjbpbVwDw8PBAQECAaXv22LZtGwIDA/GLX/zCNObn54fk5GT8/e9/R2trKwDpauuMrEbdrW13svbp0wcPPPDAXee5Ql1tzWokZ11HjRplMfb444/D19cXFRUVpjFXqKutWY26W1enNkpbXq3Vlblz56KjowOzZ89GRUUFzp07h9WrV2Pnzp14/fXXTfNOnTqFtrY2xMbGmq3v7u6OyMjILr8fkiOr0enTp+Hl5YW+fftCq9XirbfeQltbm005u5O1ubkZADp9HGDv3r2h1+tx5coVAPLX1Zasly9fNhvvbl27S6fTdfqmhLi4ONy6dQunT58GIF1tnZHVSO7a2sIV6movV6vrzZs30dDQAD8/P9OYq9a1s6xG3a2rUx840J1Xaz3yyCP4/PPPMXHiRHzwwQcAADc3N7z//vtmRyfGZ8V2tp/AwEAcPXrUZbICwODBg5GQkIDw8HA0NjZi69at+P3vf4/Tp093epQkZdb+/fvjgQcesKjJ1atX8e233wIAqqurERAQIHtdbcl68eJF9O/fH4A0de0uvV6PJ5980mL8zs86bNgwyWrbHbZmBVyjtrZwhbrawxXrWlBQgNbWVkyd+v3bRFy1rp1lBaSpq1Mb5b28WsuooqICEyZMwIMPPojly5fD09MTf/nLX/DKK6+gf//+ePbZZ8220dV+bH0lizOyAjA1UqMXX3wRs2fPxrp165CVlYXhw4c7LKtSqcTs2bPx3nvvYdGiRZg+fTpu3LiBhQsXorW1FUII07py19WerIA0de2upqYmmz6rVLXtDluzAq5RW1u4Ql3t4Wp1/eKLL/C73/0OU6dONfufKFesa1dZAWnq6tRTr2q12nQK7U5NTU2m5V1ZsGABevXqhUOHDmHatGn45S9/iR07dmDUqFHIyMgw3XZh3EZX+7H17dqOzNre3m5136+++iqA25dAOzrr22+/jZkzZ2LZsmV4+OGHERcXBw8PD8ycORMATKdIXaGutmbtir117S5bP6tUte2O7vx3AZxfW1u4Ql27S666VlRU4LnnnkNERIRFo3G1ulrL2hV76+rURnkvr9YyOnr0KJ566imL/wgTJ07ExYsXcfbsWdM+7tzmD/djbR/OytrVPYpGwcHBAG5fwOLorO7u7li3bh0uXryII0eO4PTp0/jkk09w/fp1uLm5YfDgwaZ93LnNH+7HGXW1NWtX7K1rd9n6WaWqbXd0578L4Pza2sIV6tpdctT1woULGDduHPr164ePP/4YXl5eZstdqa53y9oVe+vq1EYZFRWF06dPo6GhwWzc2qu1jNra2jo9EjNejWf8YvanP/0pevXqhdLSUrN5LS0tOHnypNV9ODtrV/7zn/8AAPz9/R2e1SggIACPPfYYBg8ejPb2dhw6dAjDhw83NXtXqKutWbtib127KzIyEmVlZRZX9B07dgxeXl4IDQ0FIF1tnZG1K86urS1coa7d5ey6Xr16FePGjUNrays+/fRT03f+d3KVutqStSt217VbN5fYyXgP3YoVK0xjxnvoRowYYRrT6/WivLzc7H65UaNGCY1GI65evWoaa2trEzExMcLHx0e0tbWZxsePHy+CgoI6vcfn008/dZmsN27cEE1NTWb77ejoEFOnThVKpVLodDqHZ+3Mu+++KxQKhdixY4fZuNx1tTWrVHX9odLS0i7v9eosr/HexG3btpnGampqxAMPPCBSUlLM1peits7I6oja2pv1TrbcRylnXW3N6gp1vXnzpoiPjxc+Pj6irKzM6rblrqutWaWqq9OfzJOcnGx6KsuaNWvEyJEjhYeHh9kNpL/61a+EQqEQ586dM4198sknQqlUisGDB4tly5aJP/3pT2LEiBFCoVCI3Nxcs32UlZUJT09PER0dLVatWiXeeOMNoVarRWJioktlPXjwoAgMDBTz588X77//vlixYoV47LHHhEKhEL/+9a+dkvWjjz4SkyZNEvn5+WLt2rUiOTlZKBQKMWvWLIt9yF1XW7NKWVchhFi5cqXIyckRc+bMEQqFQkyePFnk5OSInJwcUV9f32Xe9vZ2MWLECOHt7W32tBsfHx9x+vRps31IVVtHZ5WytveaVQhhmvf8888LhUIhZs6caRpzpbramtUV6vrss8+a8n300UdmP3/7299cqq62ZpWqrk5vlE1NTeK3v/2t0Gq1wtPTUwwfPtziSQ4vvfSSUCqVFr/h9u3bJx5//HHh5eUlVCqVeOSRR8TatWs73c/Ro0fFY489JtRqtejfv7+YO3euuHnzpktlPXPmjEhOThYPPfSQUKvVwsvLS8TFxXX5mRyRtaSkRDzxxBPC19dXqNVqERUVZXX/ctbV1qxS1lUIIUJCQoRCoRAKhUIolUqhVCpNvzbm6+r3wbVr10RaWprw8/MTXl5eYsyYMeLEiROd7keK2jo6q5S17U5W47w7/2n89Q/JXVdbsrpCXUNCQszy3fnz0EMPWexHzrramlWquvI1W0RERFbwNVtERERWsFESERFZwUZJRERkBRslERGRFWyUREREVrBREhERWcFGSUREZAUbJRERkRVslERERFawURIREVnBRklERGQFGyUREZEV/wep8Wkcr6uWHQAAAABJRU5ErkJggg==",
"text/plain": [
"PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f28c0d89590>)"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"figure(figsize=(5, 2))\n",
"plt[:hist](Float64[d[1] for d in solutions]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What is left is to code own MCMC sampler and start doing Bayesian FEM."
]
},
{
"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 etc. is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues"
]
}
],
"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