Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-27 04:10:54 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
cf8a580eedcbdd546fdf2642bfa20005d7d095e7
JuliaFEM.jl/notebooks/2015-08-29-developing-juliafem.ipynb
T
Jukka Aho cf8a580eed added description how to formulate new problems
2015-09-24 21:05:04 +03:00

794 lines
20 KiB
Plaintext
Raw Blame History

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Developing JuliaFEM\n",
"\n",
"Author(s): Jukka Aho\n",
"\n",
"**Abstract**: General developer notes."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Logger(root,DEBUG,PipeEndpoint(open, 0 bytes waiting),root)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Logging\n",
"using FactCheck\n",
"Logging.configure(level=DEBUG)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Developing own element\n",
"\n",
"Finite element definition, from [FEniCS-book](https://bitbucket.org/fenics-project/fenics-book/src/7d3a80e7dda0fc279c7964dc6000d57942f11eb3/fenicsbook.cls?at=master) [Ciarlet, 2002]:\n",
"\n",
"- the domain $T$ is a bounded, closed subset of $\\mathbb{R}^d$ (for $d = 1, 2, 3, \\dots$) with nonempty interior and piecewise smooth boundary;\n",
"- the space $\\mathcal{V} = \\mathcal{V}(T)$ is a finite dimensional function space on $T$ of dimension $n$;\n",
"- the set of degrees of freedom (nodes) $\\mathcal{L} = \\{\\ell_1, \\ell_2,\\ldots, \\ell_{n}\\}$ is a basis for the dual space $\\mathcal{V}'$; that is, the space of bounded linear functionals on $\\mathcal{V}$.\n",
"\n",
"We extend this definition so that domain $T$ can also be empty.\n",
"\n",
"Minimum requirements for element:\n",
"- subclass from Element, if not wanting 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",
"- create proper constructor (see example).\n",
"\n",
"Test the element using ``test_element`` function. It it passes, then element implementation should be fine. As an example, we define 4 node quadrilateral element using linear Lagrange basis. We really don't care much how element is implemented as long it's interface is constructed with some rules. The interface is tested using `test_element` and it also gives information how to fix element if something is missing."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM: Element"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's one basic implementation. The actual element:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"type MyQuad4 <: Element\n",
" connectivity :: Array{Int, 1}\n",
" fields :: Dict{Any, Any}\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Default constructor, providing connectivity data needed in assembly"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"MyQuad4"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MyQuad4(connectivity) = MyQuad4(connectivity, Dict{Any, Any}())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Some basic charasteristics like number of basis funcitons and dimension:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_element_dimension (generic function with 7 methods)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.get_number_of_basis_functions(el::Type{MyQuad4}) = 4\n",
"JuliaFEM.get_element_dimension(el::MyQuad4) = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The most important, it's basis (we probably want to interpolate something with this element):"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_dbasisdxi (generic function with 7 methods)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function JuliaFEM.get_basis(el::MyQuad4, xi)\n",
" [(1-xi[1])*(1-xi[2])/4\n",
" (1+xi[1])*(1-xi[2])/4\n",
" (1+xi[1])*(1+xi[2])/4\n",
" (1-xi[1])*(1+xi[2])/4]\n",
"end\n",
"\n",
"function JuliaFEM.get_dbasisdxi(el::MyQuad4, xi)\n",
" [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n",
" (1-xi[2])/4.0 -(1+xi[1])/4.0\n",
" (1+xi[2])/4.0 (1+xi[1])/4.0\n",
" -(1+xi[2])/4.0 (1-xi[1])/4.0]\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we check that everything is well defined:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"24-Sep 21:04:04:INFO:root:Testing element MyQuad4\n",
"24-Sep 21:04:04:INFO:root:number of basis functions in this element: 4\n",
"24-Sep 21:04:04:INFO:root:Initializing element\n",
"24-Sep 21:04:04:INFO:root:Element dimension: 2\n",
"24-Sep 21:04:04:INFO:root:Setting scalar field [1 2 3 4] to element.\n",
"24-Sep 21:04:04:INFO:root:Interpolating scalar field at [0.0,0.0]\n"
]
},
{
"data": {
"text/plain": [
"PipeEndpoint(open, 0 bytes waiting)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: test_element\n",
"test_element(MyQuad4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If `test_element` passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"24-Sep 21:04:04:INFO:root:Value: [2.5]\n",
"24-Sep 21:04:04:INFO:root:Element MyQuad4 passed tests.\n"
]
}
],
"source": [
"using JuliaFEM: set_field, interpolate\n",
"el1 = MyQuad4([1, 2, 3, 4])\n",
"set_field(el1, :temperature, [1, 2, 3, 4])\n",
"set_field(el1, :Geometry, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]);\n",
"set_field(el1, :\"heat coefficient\", 1);"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.5"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# temperature at the middle poinf of the element, 1/4*(1+2+3+4)\n",
"interpolate(el1, :temperature, [0.0, 0.0])"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 5.0\n",
" 0.5\n",
" 0.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# geometry midpoint of element\n",
"interpolate(el1, :Geometry, [0.0, 0.0])"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolating scalar -> scalar.\n",
"interpolate(el1, :\"heat coefficient\", [0.0, 0.0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 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 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 formulation\n",
"\n",
"Let's consider Poisson equation\n",
"\\begin{align}\n",
"\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\n",
"u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
"\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}}\n",
"\\end{align}\n",
"\n",
"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_{\\Gamma_{\\mathrm{N}}}g v\\,\\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",
"- it needs to have lhs and rhs functions\n",
"- default constructor takes the element as input argument\n",
"\n",
"Now we have function `test_equation`, which we can use to test that everything is working as expected. \n",
"\n",
"Again thanks to multiple dispatch, you are free to code your weak form however you want as long as it returns lhs and rhs sides for element dofs. This kind of freedom gives good opportunities to wrap e.g. Fortran code from some other projects. And again we have some suggestions ad following these ideas you get a lot of stuff for free. First we look the left hand side of the equation, that is,\n",
"\\begin{equation}\n",
" \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM: Equation, IntegrationPoint, Quad4\n",
"\n",
"abstract Heat <: Equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our basic data type often looks something like this:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"\"\"\"\n",
"Diffusive heat transfer for 4-node bilinear element.\n",
"\"\"\"\n",
"type DC2D4 <: Heat\n",
" element :: Quad4\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
" global_dofs :: Array{Int64, 1}\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We must provide default constructor which takes element as input argument:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"DC2D4"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function DC2D4(el::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",
" set_field(el, :temperature, zeros(2, 4)) # assign new field \"temperature\" to element\n",
" DC2D4(el, integration_points, [])\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the actual implementation for $\\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x$:"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_lhs (generic function with 2 methods)"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_element, get_dbasisdX, has_lhs, has_rhs\n",
"\n",
"\"\"\"\n",
"Left hand side defined in integration point\n",
"\"\"\"\n",
"function JuliaFEM.get_lhs(eq::DC2D4, ip)\n",
" el = get_element(eq)\n",
" dNdX = get_dbasisdX(el, ip.xi)\n",
" hc = interpolate(el, :\"temperature thermal conductivity\", ip.xi)\n",
" return dNdX*hc*dNdX'\n",
"end\n",
"JuliaFEM.has_lhs(eq::DC2D4) = true"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And that's it. If we want to play with this formulation, we must create element and assign this equation for it:"
]
},
{
"cell_type": "code",
"execution_count": 16,
"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": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
"el = Quad4([1, 2, 3, 4])\n",
"set_field(el, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
"set_field(el, :\"temperature thermal conductivity\", 6)\n",
"eq = DC2D4(el)\n",
"integrate_lhs(eq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If rhs or lhs is not defined, integration returns nothing."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"isa(integrate_rhs(eq), Void)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next heat flux on boundary:"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_rhs (generic function with 2 methods)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_basis, Seg2\n",
"\n",
"\"\"\"\n",
"Diffusive heat transfer for 2-node linear segment.\n",
"\"\"\"\n",
"type DC2D2 <: Heat\n",
" element :: Seg2\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
" global_dofs :: Array{Int64, 1}\n",
"end\n",
"\n",
"function DC2D2(el::Seg2)\n",
" integration_points = [\n",
" IntegrationPoint([0], 2.0)]\n",
" set_field(el, :temperature, zeros(2, 1))\n",
" DC2D2(el, integration_points, [])\n",
"end\n",
"\n",
"\"\"\"\n",
"Right hand side defined in integration point\n",
"\"\"\"\n",
"function JuliaFEM.get_rhs(eq::DC2D2, ip)\n",
" el = get_element(eq)\n",
" N = get_basis(el, ip.xi)\n",
" f = interpolate(el, :\"temperature flux\", ip.xi)\n",
" return f*N\n",
"end\n",
"JuliaFEM.has_rhs(eq::DC2D2) = true"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 50.0\n",
" 50.0"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"el = Seg2([1, 2])\n",
"set_field(el, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
"set_field(el, :\"temperature flux\", 100.0)\n",
"eq = DC2D2(el)\n",
"integrate_rhs(eq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Defining own problem\n",
"\n",
"- takes a set of elements and maps corresponding equations for them\n",
"- problem should have all required information in order to be solvable"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"PlaneHeatProblem"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Problem, get_equation, get_dimension\n",
"\n",
"type PlaneHeatProblem <: Problem\n",
" equations :: Array{Any, 1}\n",
"end\n",
"PlaneHeatProblem() = PlaneHeatProblem([])"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_equation (generic function with 3 methods)"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.get_dimension(pr::Type{PlaneHeatProblem}) = 1\n",
"JuliaFEM.get_equation(pr::Type{PlaneHeatProblem}, el::Type{Quad4}) = DC2D4\n",
"JuliaFEM.get_equation(pr::Type{PlaneHeatProblem}, el::Type{Seg2}) = DC2D2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our solution procedure so far is therefore"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"24-Sep 21:04:17:DEBUG:root:Problem (matrix) dimension: 4\n"
]
},
{
"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",
"[50.0,50.0,0.0,0.0])"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: set_global_dofs!, get_global_dofs, add_element!, get_equations, get_matrix_dimension\n",
"\n",
"# create elements and add necessary properties like connectivity and geometry\n",
"el1 = Quad4([1, 2, 3, 4])\n",
"set_field(el1, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
"set_field(el1, :\"temperature thermal conductivity\", 6)\n",
"el2 = Seg2([1, 2])\n",
"set_field(el2, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
"set_field(el2, :\"temperature flux\", 100.0)\n",
"\n",
"problem = PlaneHeatProblem()\n",
"add_element!(problem, el1)\n",
"add_element!(problem, el2)\n",
"\n",
"# set global dofs for equations\n",
"set_global_dofs!(problem)\n",
"\n",
"n = get_matrix_dimension(problem)\n",
"\n",
"# integrate and assembly\n",
"A = zeros(n, n)\n",
"b = zeros(n)\n",
"for eq in get_equations(problem)\n",
" dofs = get_global_dofs(eq)\n",
" if has_lhs(eq)\n",
" A[dofs, dofs] += integrate_lhs(eq)\n",
" end\n",
" if has_rhs(eq)\n",
" b[dofs] += integrate_rhs(eq)\n",
" end\n",
"end\n",
"A, b"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.5.0-dev",
"language": "julia",
"name": "julia-0.5"
},
"language_info": {
"name": "julia",
"version": "0.5.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Reference in New Issue View Git Blame Copy Permalink