Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-28 20:46:58 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
e0d974b3a9f3f1ea06fe38d4eff571c728eb6821
JuliaFEM.jl/notebooks/2015-08-29-developing-juliafem.ipynb
T
Jukka Aho e0d974b3a9 data types defined
2015-10-09 01:44:08 +03:00

1698 lines
44 KiB
Plaintext
Raw Blame History

This file contains ambiguous Unicode characters
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 general guidelines how to implement things to JuliaFEM.\n",
"\n",
"In short, we have\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 which needs to be solve. One element can have several equations but not vice versa. Typical equations: Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc.\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.\n",
"\n",
"Moreover, we have (will have) **`test_element`**, **`test_equation`**, **`test_problem`** and **`test_solver`** which can be used to test that implementation has all necessary things defined."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n",
"WARNING: Base.String is deprecated, use AbstractString instead.\n"
]
},
{
"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",
"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": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuliaFEM: Element, Field, FieldSet, Basis"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here's the basic implementation for element:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"type MyQuad4 <: Element\n",
" connectivity :: Array{Int, 1}\n",
" basis :: Basis\n",
" fields :: Dict{Symbol, 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": 4,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"MyQuad4"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function MyQuad4(connectivity)\n",
" h(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",
" dh(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",
" 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."
]
},
{
"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::Type{MyQuad4}) = 2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we check that everything is well defined:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"09-Oct 01:41:19:INFO:root:Testing element MyQuad4\n",
"09-Oct 01:41:19:INFO:root:number of basis functions in this element: 4\n",
"09-Oct 01:41:19:INFO:root:Initializing element\n",
"09-Oct 01:41:19:INFO:root:Element dimension: 2\n",
"09-Oct 01:41:19:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n",
"09-Oct 01:41:19:INFO:root:Pushing field to element.\n",
"09-Oct 01:41:19:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
"09-Oct 01:41:19:INFO:root:Value: 2.5\n"
]
},
{
"data": {
"text/plain": [
"PipeEndpoint(open, 0 bytes waiting)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"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": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{JuliaFEM.Field{T},1}:\n",
" JuliaFEM.Field{Int64}(0.0,1,2)\n",
" JuliaFEM.Field{Int64}(1.0,1,3)"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"09-Oct 01:41:19:INFO:root:Element MyQuad4 passed tests.\n"
]
}
],
"source": [
"using JuliaFEM: new_fieldset!, add_field!, interpolate, dinterpolate\n",
"el1 = MyQuad4([1, 2, 3, 4])\n",
"new_fieldset!(el1, \"temperature\")\n",
"add_field!(el1, \"temperature\", Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
"add_field!(el1, \"temperature\", Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
"new_fieldset!(el1, \"geometry\")\n",
"add_field!(el1, \"geometry\", Field(0.0, 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",
"new_fieldset!(el1, \"heat coefficient\")\n",
"add_field!(el1, \"heat coefficient\", Field(0.0, 2))\n",
"add_field!(el1, \"heat coefficient\", Field(1.0, 3))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1.25"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# temperature at the middle poinf of the element, 1/4*(1+2+3+4) at t=0.5\n",
"interpolate(el1, \"temperature\", [0.0, 0.0], 0.5)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Float64,1}:\n",
" 5.0\n",
" 0.5\n",
" 0.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# geometry midpoint of element\n",
"interpolate(el1, \"geometry\", [0.0, 0.0], 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3x2 Array{Float64,2}:\n",
" 5.0 0.0\n",
" 0.0 0.5\n",
" 0.0 0.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolate derivatives works too\n",
"dinterpolate(el1, \"geometry\", [0.0, 0.0], 0.0)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2.5"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# interpolating scalar -> scalar.\n",
"interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"3"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)"
]
},
{
"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 equation\n",
"\n",
"Let's consider a Laplace equation\n",
"\\begin{align}\n",
"\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\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",
"- provide the name of the unknown field variable trying to solve\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": 13,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_unknown_field_name (generic function with 1 method)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Equation, IntegrationPoint, Quad4\n",
"\n",
"abstract Heat <: Equation\n",
"\n",
"get_unknown_field_name(eq::Heat) = symbol(\"temperature\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our basic data type often looks something like this:"
]
},
{
"cell_type": "code",
"execution_count": 14,
"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": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"DC2D4"
]
},
"execution_count": 15,
"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",
" new_fieldset!(el, \"temperature\") # 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": 16,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_lhs (generic function with 2 methods)"
]
},
"execution_count": 16,
"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, t)\n",
" el = get_element(eq)\n",
" dNdX = get_dbasisdX(el, ip.xi, t)\n",
" hc = interpolate(el, \"temperature thermal conductivity\", ip.xi, t)\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": 17,
"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": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: integrate, integrate_lhs, integrate_rhs, new_fieldset!, add_field!\n",
"el = Quad4([1, 2, 3, 4])\n",
"new_fieldset!(el, \"geometry\")\n",
"add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]))\n",
"new_fieldset!(el, \"temperature thermal conductivity\")\n",
"add_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
"eq = DC2D4(el)\n",
"integrate_lhs(eq, 1.0)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"JuliaFEM.has_lhs(eq)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If rhs or lhs is not defined, integration returns nothing."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"isa(integrate_rhs(eq, 1.0), Void)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next heat flux on boundary:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"has_rhs (generic function with 2 methods)"
]
},
"execution_count": 20,
"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",
" new_fieldset!(el, \"temperature\")\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, t)\n",
" el = get_element(eq)\n",
" h = get_basis(el, ip.xi)\n",
" #f = el[\"temperature flux\"]\n",
" f = interpolate(el, \"temperature flux\", ip.xi, t)\n",
" return h*f\n",
"end\n",
"JuliaFEM.has_rhs(eq::DC2D2) = true"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 50.0\n",
" 50.0"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"el = Seg2([1, 2])\n",
"new_fieldset!(el, \"geometry\")\n",
"add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n",
"new_fieldset!(el, \"temperature flux\")\n",
"add_field!(el, \"temperature flux\", Field(0.0, 100.0))\n",
"eq = DC2D2(el)\n",
"integrate_rhs(eq, 1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Defining own problem\n",
"\n",
"- main object: takes a set of elements and maps corresponding field equations to them"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"PlaneHeatProblem"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: Problem, get_equation, get_dimension\n",
"\n",
"type PlaneHeatProblem <: Problem\n",
" equations :: Array{Equation, 1}\n",
"end\n",
"PlaneHeatProblem() = PlaneHeatProblem([])"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_equation (generic function with 3 methods)"
]
},
"execution_count": 23,
"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": 24,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"09-Oct 01:41:22:DEBUG:root:total dofs: 4\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"dofmap: Dict(4=>[3],2=>[1],3=>[2],5=>[4])\n"
]
},
{
"data": {
"text/plain": [
"4x1 sparse matrix with 2 Float64 entries:\n",
"\t[1, 1] = 300.0\n",
"\t[2, 1] = 300.0"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using JuliaFEM: get_connectivity, set_global_dofs!, get_global_dofs\n",
"using JuliaFEM: add_element!, get_equations, get_matrix_dimension\n",
"\n",
"# create elements and add necessary properties like connectivity and geometry\n",
"el1 = Quad4([2, 3, 4, 5])\n",
"new_fieldset!(el1, \"geometry\")\n",
"add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
"new_fieldset!(el1, \"temperature thermal conductivity\")\n",
"add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
"el2 = Seg2([2, 3])\n",
"new_fieldset!(el2, \"geometry\")\n",
"add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
"new_fieldset!(el2, \"temperature flux\")\n",
"add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
"\n",
"problem = PlaneHeatProblem()\n",
"add_element!(problem, el1)\n",
"add_element!(problem, el2)\n",
"\n",
"function JuliaFEM.get_connectivity(pr::Problem)\n",
" conn = Int[]\n",
" for eq in get_equations(pr)\n",
" el = get_element(eq)\n",
" append!(conn, get_connectivity(el))\n",
" end\n",
" conn = unique(conn)\n",
" return conn\n",
"end\n",
"\n",
"\"\"\"\n",
"Calculate global dofs for equations, maybe using some bandwidth\n",
"minimizing or fill reducing algorithm\n",
"\"\"\"\n",
"function calculate_global_dofs(pr::Problem)\n",
" conn = get_connectivity(pr)\n",
" dim = get_dimension(typeof(pr))\n",
" ndofs = dim*length(conn)\n",
" Logging.debug(\"total dofs: $ndofs\")\n",
"\n",
" mconn = maximum(conn)\n",
" gdofs = reshape(collect(1:mconn), dim, mconn)\n",
" dofmap = Dict{Int64, Array{Int64, 1}}()\n",
" for (i, c) in enumerate(conn)\n",
" dofmap[c] = gdofs[:, i]\n",
" end\n",
" return dofmap\n",
"end\n",
"\n",
"\"\"\"\n",
"Assign global dofs for equations.\n",
"\"\"\"\n",
"function assign_global_dofs!(pr::Problem, dofmap)\n",
" for eq in get_equations(pr)\n",
" el = get_element(eq)\n",
" c = get_connectivity(el)\n",
" #gdofs = [dofmap[ci] for ci in c]\n",
" gdofs = Int64[]\n",
" for ci in c\n",
" append!(gdofs, dofmap[ci])\n",
" end\n",
" set_global_dofs!(eq, gdofs)\n",
" end\n",
"end\n",
"\n",
"dofmap = calculate_global_dofs(problem)\n",
"println(\"dofmap: $dofmap\")\n",
"assign_global_dofs!(problem, dofmap)\n",
"\n",
"function get_lhs(pr::Problem, t::Float64)\n",
" I = Int64[]\n",
" J = Int64[]\n",
" V = Float64[]\n",
" dim = get_dimension(typeof(pr))\n",
" for eq in filter(has_lhs, get_equations(pr))\n",
" dofs = get_global_dofs(eq)\n",
" lhs = integrate_lhs(eq, t)\n",
" for (li, i) in enumerate(dofs)\n",
" for (lj, j) in enumerate(dofs)\n",
" push!(I, i)\n",
" push!(J, j)\n",
" push!(V, lhs[li, lj])\n",
" end\n",
" end\n",
" end\n",
" return I, J, V\n",
"end\n",
"\n",
"function get_rhs(pr::Problem, t::Float64)\n",
" I = Int64[]\n",
" V = Float64[]\n",
" dim = get_dimension(typeof(pr))\n",
" for eq in filter(has_rhs, get_equations(pr))\n",
" dofs = get_global_dofs(eq)\n",
" rhs = integrate_rhs(eq, t)\n",
" for (li, i) in enumerate(dofs)\n",
" push!(I, i)\n",
" push!(V, rhs[li])\n",
" end\n",
" end\n",
" return I, V\n",
"end\n",
"\n",
"t = 1.0\n",
"A = sparse(get_lhs(problem, t)...)\n",
"b = sparsevec(get_rhs(problem, t)..., size(A, 1))"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x1 sparse matrix with 2 Float64 entries:\n",
"\t[1, 1] = 100.0\n",
"\t[2, 1] = 100.0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"fdofs = [1, 2]\n",
"A[fdofs, fdofs] \\ b[fdofs]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We still need to consider Dirichlet boundary conditions:\n",
"\\begin{align}\n",
"u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"DirichletProblem"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abstract BoundaryProblem <: Problem\n",
"\n",
"type DirichletProblem <: BoundaryProblem\n",
" equations :: Array{Equation, 1}\n",
"end\n",
"DirichletProblem() = DirichletProblem([])"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"get_equation (generic function with 4 methods)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abstract DirichletBC <: Equation\n",
"\n",
"get_unknown_field_name(eq::DirichletBC) = symbol(\"reaction force\")\n",
"\n",
"\"\"\"\n",
"Dirichlet boundary condition element for 2 node line segment\n",
"(plane problems).\n",
"\"\"\"\n",
"type DBC2D2 <: DirichletBC\n",
" element :: Seg2\n",
" integration_points :: Array{IntegrationPoint, 1}\n",
" global_dofs :: Array{Int64, 1}\n",
" fieldval :: Function\n",
"end\n",
"function DBC2D2(el::Seg2)\n",
" integration_points = [\n",
" IntegrationPoint([-sqrt(1/3)], 1.0),\n",
" IntegrationPoint([+sqrt(1/3)], 1.0)]\n",
" new_fieldset!(el, \"reaction force\")\n",
" fieldval(X, t) = 0.0\n",
" DBC2D2(el, integration_points, [], fieldval)\n",
"end\n",
"\n",
"function JuliaFEM.get_lhs(eq::DBC2D2, ip, t)\n",
" el = get_element(eq)\n",
" h = get_basis(el)(ip.xi)\n",
" return h*h'\n",
"end\n",
"\n",
"function JuliaFEM.get_rhs(eq::DBC2D2, ip, t)\n",
" el = get_element(eq)\n",
" h = get_basis(el, ip.xi)\n",
" f = eq.fieldval\n",
" X = interpolate(el, \"geometry\", ip.xi, t)\n",
" return h*f(X, t)\n",
"end\n",
"JuliaFEM.has_lhs(eq::DBC2D2) = true\n",
"JuliaFEM.has_rhs(eq::DBC2D2) = true\n",
"JuliaFEM.get_dimension(pr::Type{DirichletProblem}) = 1\n",
"JuliaFEM.get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"1-element Array{JuliaFEM.Equation,1}:\n",
" DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[],:geometry=>JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# create elements and add necessary properties like connectivity and geometry\n",
"el3 = Seg2([4, 5])\n",
"new_fieldset!(el3, \"geometry\")\n",
"add_field!(el3, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
"\n",
"bc1 = DirichletProblem()\n",
"#get_equation(typeof(bc1), typeof(el3))\n",
"#DBC2D2(el3)\n",
"add_element!(bc1, el3)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x1 sparse matrix with 2 Float64 entries:\n",
"\t[3, 1] = 0.0\n",
"\t[4, 1] = 0.0"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"assign_global_dofs!(bc1, dofmap)\n",
"\n",
"# integrate and assembly\n",
"t = 1.0\n",
"A2 = sparse(get_lhs(bc1, t)...)\n",
"b2 = sparsevec(get_rhs(bc1, t)...)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Now we have two problems defined, "
]
},
{
"cell_type": "code",
"execution_count": 30,
"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": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"full(A)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x1 Array{Float64,2}:\n",
" 300.0\n",
" 300.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"full(b)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x4 Array{Float64,2}:\n",
" 0.0 0.0 0.0 0.0 \n",
" 0.0 0.0 0.0 0.0 \n",
" 0.0 0.0 0.333333 0.166667\n",
" 0.0 0.0 0.166667 0.333333"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"full(A2)"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"4x1 Array{Float64,2}:\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"full(b2)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 sparse matrix with 24 Float64 entries:\n",
"\t[1, 1] = 4.0\n",
"\t[2, 1] = -1.0\n",
"\t[3, 1] = -2.0\n",
"\t[4, 1] = -1.0\n",
"\t[1, 2] = -1.0\n",
"\t[2, 2] = 4.0\n",
"\t[3, 2] = -1.0\n",
"\t[4, 2] = -2.0\n",
"\t[1, 3] = -2.0\n",
"\t[2, 3] = -1.0\n",
"\t⋮\n",
"\t[8, 3] = 0.166667\n",
"\t[1, 4] = -1.0\n",
"\t[2, 4] = -2.0\n",
"\t[3, 4] = -1.0\n",
"\t[4, 4] = 4.0\n",
"\t[7, 4] = 0.166667\n",
"\t[8, 4] = 0.333333\n",
"\t[3, 7] = 0.333333\n",
"\t[4, 7] = 0.166667\n",
"\t[3, 8] = 0.166667\n",
"\t[4, 8] = 0.333333"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Atot = [A A2; A2 zeros(A2)]"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x1 sparse matrix with 4 Float64 entries:\n",
"\t[1, 1] = 300.0\n",
"\t[2, 1] = 300.0\n",
"\t[7, 1] = 0.0\n",
"\t[8, 1] = 0.0"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"btot = [b; b2]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Problem is that now are total matrix Atot has zero rows which needs to be removed."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"8x8 Array{Float64,2}:\n",
" 4.0 -1.0 -2.0 -1.0 0.0 0.0 0.0 0.0 \n",
" -1.0 4.0 -1.0 -2.0 0.0 0.0 0.0 0.0 \n",
" -2.0 -1.0 4.0 -1.0 0.0 0.0 0.333333 0.166667\n",
" -1.0 -2.0 -1.0 4.0 0.0 0.0 0.166667 0.333333\n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n",
" 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n",
" 0.0 0.0 0.333333 0.166667 0.0 0.0 0.0 0.0 \n",
" 0.0 0.0 0.166667 0.333333 0.0 0.0 0.0 0.0 "
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"full(Atot)"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Non-zero rows: [1,2,3,4,7,8]\n"
]
},
{
"data": {
"text/plain": [
"8x1 Array{Float64,2}:\n",
" 100.0\n",
" 100.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 600.0\n",
" 600.0"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"r = unique(rowvals(Atot))\n",
"println(\"Non-zero rows: $r\")\n",
"xtot = zeros(btot)\n",
"F = lufact(Atot[r,r])\n",
"s = full(btot[r])\n",
"xtot[r] = F \\ s\n",
"full(xtot)"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Developing own solver\n",
"\n",
"Last part. Defining own solver.\n",
"\n",
"- takes a set of problems (typically main problem + boundary problems)\n",
"- solves them, updates fields"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"add_problem! (generic function with 1 method)"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"abstract Solver\n",
"\n",
"\"\"\"\n",
"Simple solver for educational purposes.\n",
"\"\"\"\n",
"type SimpleSolver <: Solver\n",
" problems\n",
"end\n",
"\n",
"\"\"\"\n",
"Default initializer\n",
"\"\"\"\n",
"function SimpleSolver()\n",
" SimpleSolver(Problem[])\n",
"end\n",
"\n",
"function add_problem!(solver::Solver, problem::Problem)\n",
" push!(solver.problems, problem)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"call (generic function with 1253 methods)"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"get_problems(s::Solver) = s.problems\n",
"\"\"\"\n",
"Call solver to solve a set of problems.\n",
"\n",
"This is simple serial solver for demonstration purposes. It handles the most\n",
"common situation, i.e., some main field problem and it's Dirichlet boundary.\n",
"\"\"\"\n",
"function call(solver::SimpleSolver, t)\n",
" problems = get_problems(solver)\n",
" problem1 = problems[1]\n",
" problem2 = problems[2]\n",
"\n",
" # calculate order of degrees of freedom in global matrix\n",
" # and set the ordering to problems\n",
" dofmap = calculate_global_dofs(problem1)\n",
" assign_global_dofs!(problem1, dofmap)\n",
" assign_global_dofs!(problem2, dofmap)\n",
"\n",
" # assemble problem 1\n",
" A1 = sparse(get_lhs(problem1, t)...)\n",
" b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))\n",
"\n",
" # assemble problem 2\n",
" A2 = sparse(get_lhs(problem2, t)...)\n",
" b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))\n",
" \n",
" # make one monolithic assembly\n",
" A = [A1 A2; A2' zeros(A2)]\n",
" b = [b1; b2]\n",
"\n",
" # solve problem\n",
" nz = unique(rowvals(A))\n",
" x = zeros(b)\n",
" x[nz] = lufact(Atot[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",
" # check residual\n",
" R1 = A1*x1 - b1\n",
" R2 = A2*x2 - b2\n",
" println(\"Residual norm: $(norm(R1+R2))\")\n",
"\n",
" # update field for elements in problem 1\n",
" for equation in get_equations(problem1)\n",
" gdofs = get_global_dofs(equation)\n",
" element = get_element(equation)\n",
" field_name = get_unknown_field_name(equation) # field we are solving\n",
" field = Field(t, full(x1[gdofs])[:])\n",
" add_field!(element, field_name, field)\n",
" end\n",
"\n",
" # update field for elements in problem 2\n",
" for equation in get_equations(problem2)\n",
" gdofs = get_global_dofs(equation)\n",
" element = get_element(equation)\n",
" field_name = get_unknown_field_name(equation)\n",
" field = Field(t, full(x2[gdofs]))\n",
" add_field!(element, field_name, field)\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary\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": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Array{Float64,1},1}:\n",
" [-1.0,0.0]\n",
" [-0.5,0.5]\n",
" [0.0,1.0] \n",
" [0.5,1.5] \n",
" [1.0,2.0] "
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"\"\"\"\n",
"Simple linspace extension to multidimensional values.\n",
"\"\"\"\n",
"function Base.linspace(X1, X2, n)\n",
" [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]\n",
"end\n",
"\n",
"linspace([-1.0, 0.0], [1.0, 2.0], 5)"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"09-Oct 01:41:25:DEBUG:root:total dofs: 4\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Residual norm: 9.845568954283847e-14\n",
"Array{T,1}[[0.0,0.0],[0.25,0.0],[0.5,0.0],[0.75,0.0],[1.0,0.0]]\n",
"[100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003]\n"
]
}
],
"source": [
"using JuliaFEM: get_fieldset\n",
"\n",
"# Define Problem 1:\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",
"el1 = Quad4([1, 2, 3, 4])\n",
"new_fieldset!(el1, \"geometry\")\n",
"add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
"new_fieldset!(el1, \"temperature thermal conductivity\")\n",
"add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
"\n",
"el2 = Seg2([1, 2])\n",
"new_fieldset!(el2, \"geometry\")\n",
"add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]]))\n",
"\n",
"# Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n",
"new_fieldset!(el2, \"temperature flux\")\n",
"add_field!(el2, \"temperature flux\", Field(0.0, 0.0))\n",
"add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
"\n",
"problem1 = PlaneHeatProblem()\n",
"add_element!(problem1, el1)\n",
"add_element!(problem1, el2)\n",
"\n",
"# Define Problem 2:\n",
"# - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n",
"\n",
"el3 = Seg2([3, 4])\n",
"new_fieldset!(el3, \"geometry\")\n",
"add_field!(el3, \"geometry\", Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]]))\n",
"\n",
"problem2 = DirichletProblem()\n",
"add_element!(problem2, el3)\n",
"\n",
"# Create a solver for a set of problems\n",
"solver = SimpleSolver()\n",
"add_problem!(solver, problem1)\n",
"add_problem!(solver, problem2)\n",
"\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 = linspace([-1.0], [1.0], 5)\n",
"X = interpolate(el2, \"geometry\", xi, 1.0)\n",
"T = interpolate(el2, \"temperature\", xi, 1.0)\n",
"println(X)\n",
"println(T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## TODO\n",
"\n",
"- dynamics\n",
"- different basis for trial and test (Petrov-Galerkin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Frequently asked questions\n",
"\n",
"Any question about data structures, better ideas and improvements are very welcome."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.0-rc2",
"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