From cf8a580eedcbdd546fdf2642bfa20005d7d095e7 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 24 Sep 2015 21:05:04 +0300 Subject: [PATCH] added description how to formulate new problems --- .../2015-08-29-developing-juliafem.ipynb | 295 +++++++++--------- src/JuliaFEM.jl | 4 +- src/equations.jl | 5 +- src/problems.jl | 52 +++ 4 files changed, 215 insertions(+), 141 deletions(-) create mode 100644 src/problems.jl diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index 74686ea..b2e99d7 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -64,18 +64,7 @@ "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Creating element Point1: 1 node point element\n", - "Creating element Seg2: 2 node linear line element\n", - "Creating Lagrange basis for element Seg2. Number of basis functions: 2. Element dimension: 1\n", - "Calculating inverse of A\n" - ] - } - ], + "outputs": [], "source": [ "using JuliaFEM: Element" ] @@ -93,27 +82,7 @@ "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Element JuliaFEM.Seg2 created.\n", - "Creating element Seg3: 3 node quadratic line element\n", - "Creating Lagrange basis for element Seg3. Number of basis functions: 3. Element dimension: 1\n", - "Calculating inverse of A\n", - "Element JuliaFEM.Seg3 created.\n", - "Creating element Quad4: 4 node bilinear quadrangle element\n", - "Creating Lagrange basis for element Quad4. Number of basis functions: 4. Element dimension: 2\n", - "Calculating inverse of A\n", - "Element JuliaFEM.Quad4 created.\n", - "Creating element Tet10: 10 node quadratic tetrahedron\n", - "Creating Lagrange basis for element Tet10. Number of basis functions: 10. Element dimension: 3\n", - "Calculating inverse of A\n", - "Element JuliaFEM.Tet10 created.\n" - ] - } - ], + "outputs": [], "source": [ "type MyQuad4 <: Element\n", " connectivity :: Array{Int, 1}\n", @@ -167,7 +136,7 @@ { "data": { "text/plain": [ - "get_element_dimension (generic function with 6 methods)" + "get_element_dimension (generic function with 7 methods)" ] }, "execution_count": 5, @@ -197,7 +166,7 @@ { "data": { "text/plain": [ - "get_dbasisdxi (generic function with 6 methods)" + "get_dbasisdxi (generic function with 7 methods)" ] }, "execution_count": 6, @@ -239,12 +208,12 @@ "name": "stderr", "output_type": "stream", "text": [ - "10-Sep 20:16:49:INFO:root:Testing element MyQuad4\n", - "10-Sep 20:16:49:INFO:root:number of basis functions in this element: 4\n", - "10-Sep 20:16:49:INFO:root:Initializing element\n", - "10-Sep 20:16:49:INFO:root:Element dimension: 2\n", - "10-Sep 20:16:49:INFO:root:Setting scalar field [1 2 3 4] to element.\n", - "10-Sep 20:16:50:INFO:root:Interpolating scalar field at [0.0,0.0]\n" + "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" ] }, { @@ -256,14 +225,6 @@ "execution_count": 7, "metadata": {}, "output_type": "execute_result" - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "10-Sep 20:16:50:INFO:root:Value: [2.5]\n", - "10-Sep 20:16:50:INFO:root:Element MyQuad4 passed tests.\n" - ] } ], "source": [ @@ -284,12 +245,21 @@ "metadata": { "collapsed": false }, - "outputs": [], + "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, [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, :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);" ] }, @@ -339,12 +309,12 @@ ], "source": [ "# geometry midpoint of element\n", - "interpolate(el1, :geometry, [0.0, 0.0])" + "interpolate(el1, :Geometry, [0.0, 0.0])" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -355,7 +325,7 @@ "1" ] }, - "execution_count": 12, + "execution_count": 11, "metadata": {}, "output_type": "execute_result" } @@ -365,6 +335,15 @@ "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": {}, @@ -398,7 +377,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 12, "metadata": { "collapsed": false }, @@ -418,7 +397,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 13, "metadata": { "collapsed": false }, @@ -443,7 +422,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -454,7 +433,7 @@ "DC2D4" ] }, - "execution_count": 17, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -466,7 +445,7 @@ " 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))\n", + " set_field(el, :temperature, zeros(2, 4)) # assign new field \"temperature\" to element\n", " DC2D4(el, integration_points, [])\n", "end" ] @@ -480,7 +459,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -491,13 +470,13 @@ "has_lhs (generic function with 2 methods)" ] }, - "execution_count": 18, + "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: get_element, get_dbasisdX\n", + "using JuliaFEM: get_element, get_dbasisdX, has_lhs, has_rhs\n", "\n", "\"\"\"\n", "Left hand side defined in integration point\n", @@ -520,7 +499,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 16, "metadata": { "collapsed": false }, @@ -535,7 +514,7 @@ " -1.0 -2.0 -1.0 4.0" ] }, - "execution_count": 19, + "execution_count": 16, "metadata": {}, "output_type": "execute_result" } @@ -543,7 +522,7 @@ "source": [ "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n", "el = Quad4([1, 2, 3, 4])\n", - "set_field(el, :geometry, [0 0; 1 0; 1 1; 0 1]')\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)" @@ -558,7 +537,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -569,7 +548,7 @@ "true" ] }, - "execution_count": 20, + "execution_count": 17, "metadata": {}, "output_type": "execute_result" } @@ -587,7 +566,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -598,7 +577,7 @@ "has_rhs (generic function with 2 methods)" ] }, - "execution_count": 21, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } @@ -636,7 +615,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -649,14 +628,14 @@ " 50.0" ] }, - "execution_count": 22, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "el = Seg2([1, 2])\n", - "set_field(el, :geometry, [0.0 0.0; 0.0 1.0]')\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)" @@ -669,94 +648,132 @@ "## Defining own problem\n", "\n", "- takes a set of elements and maps corresponding equations for them\n", - "- boundary conditions\n", - "- solve!(problem) updates fields" + "- problem should have all required information in order to be solvable" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 20, "metadata": { - "collapsed": true + "collapsed": false }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "PlaneHeatProblem" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ - "abstract Problem\n", + "using JuliaFEM: Problem, get_equation, get_dimension\n", "\n", - "type HeatProblem <: Problem\n", - " elements: Array{Any, 1}\n", - " equations: Array{Any, 1}\n", - " boundary_conditions: Array{Any, 1}\n", - "end" + "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": [ - "Developing boundary conditions" + "Our solution procedure so far is therefore" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 23, "metadata": { - "collapsed": true + "collapsed": false }, - "outputs": [], + "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": [ - "abstract BC\n", + "using JuliaFEM: set_global_dofs!, get_global_dofs, add_element!, get_equations, get_matrix_dimension\n", "\n", - "type DefaultBC <: BC\n", - " \n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "bc = NodalBC()\n", - "bc[:displacement, [1, 2, 3], 1:2] = 3.0" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "bc = WeakBC()\n", - "bc[:displacement, [el1, el2, el3]] = (X) -> X[0] + X[1] - 2" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "bc = MPCBC()\n", - "bc[:displacement, " - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [ - "bc = MortarBC()\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" ] } ], diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 41f90a7..15b3de4 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -10,7 +10,9 @@ using Logging include("types.jl") # type definitions include("elements.jl") # elements include("equations.jl") # formulations -include("math.jl") # basic mathematical operations +include("problems.jl") # problems + +include("math.jl") # basic mathematical operations -- obsolete ..? include("elasticity_solver.jl") include("xdmf.jl") diff --git a/src/equations.jl b/src/equations.jl index 0ab54db..1676955 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -27,13 +27,16 @@ has_rhs(eq::Equation) = false get_rhs(eq::Equation, xi) = nothing get_element(eq::Equation) = eq.element get_integration_points(eq::Equation) = eq.integration_points + # couple convenient functions -- could make weak form definition easier +get_connectivity(eq::Equation) = get_connectivity(get_element(eq)) get_basis(eq::Equation, ip::IntegrationPoint) = get_basis(get_element(eq), ip.xi) get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq), ip.xi) interpolate(eq::Equation, field::Union(ASCIIString, Symbol), ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi) integrate_lhs(eq::Equation) = has_lhs(eq) ? integrate(eq, get_lhs) : nothing integrate_rhs(eq::Equation) = has_rhs(eq) ? integrate(eq, get_rhs) : nothing + """ Return determinant of Jacobian for numerical integration. """ @@ -72,7 +75,7 @@ function get_global_dofs(eq::Equation) eq.global_dofs end -function set_global_dofs(eq::Equation, dofs) +function set_global_dofs!(eq::Equation, dofs) eq.global_dofs = dofs end diff --git a/src/problems.jl b/src/problems.jl new file mode 100644 index 0000000..914df2a --- /dev/null +++ b/src/problems.jl @@ -0,0 +1,52 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +abstract Problem + +get_equations(pr::Problem) = pr.equations +get_dimension(pr::Type{Problem}) = nothing +get_equation(pr::Type{Problem}, el::Type{Element}) = nothing + +""" +Add new element to problem +""" +function add_element!(pr::Problem, el::Element) + eq = get_equation(typeof(pr), typeof(el)) + push!(pr.equations, eq(el)) +end + +""" +Return total number of basis functions in problem +""" +function get_number_of_basis_functions(pr::Problem) + conn = Int[] + for eq in get_equations(pr) + append!(conn, get_connectivity(eq)) + end + length(unique(conn)) +end + +""" +Problem matrix size dimension +""" +function get_matrix_dimension(pr::Problem) + get_dimension(typeof(pr))*get_number_of_basis_functions(pr) +end + +""" +Assign global dofs for element. This doesn't do any reordering. +""" +function set_global_dofs!(pr::Problem) + #ndim = get_dimension(pr)*get_number_of_basis_functions(pr) + #ndim = get_matrix_dimension(pr) + dim = get_dimension(typeof(pr)) + nconn = get_number_of_basis_functions(pr) + ndim = dim*nconn + Logging.debug("Problem (matrix) dimension: $ndim") + gdofs = reshape(collect(1:ndim), dim, nconn) + for eq in get_equations(pr) + lconn = get_connectivity(eq) + gconn = gdofs[:, lconn][:] + set_global_dofs!(eq, gconn) + end +end