diff --git a/README.rst b/README.rst index 27308c2..30a61ae 100644 --- a/README.rst +++ b/README.rst @@ -20,15 +20,15 @@ The JuliaFEM project develops open-source software for reliable, scalable, distr .. image:: https://coveralls.io/repos/JuliaFEM/JuliaFEM.jl/badge.svg?branch=master :target: https://coveralls.io/r/JuliaFEM/JuliaFEM.jl?branch=master -.. image:: http://juliafem.kapsi.fi/_static/notebooks-status.svg?ts=1 - :target: http://juliafem.kapsi.fi/tutorials/index.html +.. image:: http://juliafem.org/_static/notebooks-status.svg?ts=1 + :target: http://juliafem.org/tutorials/index.html -.. image:: http://juliafem.kapsi.fi/_static/doctests-status.svg?ts=1 +.. image:: http://juliafem.org/_static/doctests-status.svg?ts=1 -.. image:: http://juliafem.kapsi.fi/_static/unittests-status.svg?ts=1 +.. image:: http://juliafem.org/_static/unittests-status.svg?ts=1 -.. image:: http://juliafem.kapsi.fi/_static/lint-status.svg?ts=1 - :target: http://juliafem.kapsi.fi/quality/index.html#lint-report +.. image:: http://juliafem.org/_static/lint-status.svg?ts=1 + :target: http://juliafem.org/quality/index.html#lint-report .. image:: https://img.shields.io/github/issues/JuliaFEM/JuliaFEM.jl.svg :target: https://github.com/JuliaFEM/JuliaFEM.jl/issues diff --git a/deps/ReST.jl b/deps/ReST.jl new file mode 100644 index 0000000..dff65bf --- /dev/null +++ b/deps/ReST.jl @@ -0,0 +1,3 @@ +module ReST + +end diff --git a/docs/api/JuliaFEM.abaqus_reader.rst b/docs/api/JuliaFEM.abaqus_reader.rst new file mode 100644 index 0000000..dc71683 --- /dev/null +++ b/docs/api/JuliaFEM.abaqus_reader.rst @@ -0,0 +1,12 @@ +JuliaFEM.abaqus_reader +====================== + +Internal +-------- + + .. function:: add_handler(section, function_name) + +source: +------- +[JuliaFEM/src/abaqus_reader.jl:18](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/abaqus_reader.jl#L18) + diff --git a/docs/api/JuliaFEM.elasticity_solver.rst b/docs/api/JuliaFEM.elasticity_solver.rst new file mode 100644 index 0000000..4a8532b --- /dev/null +++ b/docs/api/JuliaFEM.elasticity_solver.rst @@ -0,0 +1,93 @@ +JuliaFEM.elasticity_solver +========================== + +Internal +-------- + + .. function:: assemble!(fe, eldofs_, I, V) + + Assemble global RHS to I,V ready for sparse format + + :param fe : local vector + :param eldofs_ : Array + :param I,V : Arrays for sparse matrix + :notes: eldofs can also be node ids for convenience. In that case dimension + is calculated and eldofs are "extended" to problem dimension. +source: +------- +[JuliaFEM/src/elasticity_solver.jl:171](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L171) + + .. function:: assemble!(ke, eldofs_, I, J, V) + + Assemble global stiffness matrix to I,J,V ready for sparse format + + :param ke : local matrix + :param eldofs_ : Array + :param I,J,V : Arrays for sparse matrix + :notes: eldofs can also be node ids for convenience. In that case dimension + is calculated and eldofs are "extended" to problem dimension. +source: +------- +[JuliaFEM/src/elasticity_solver.jl:130](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L130) + + .. function:: calc_local_matrices!(X, u, R, Kt, N, dNdchi, lambda_, mu_, ipoints, iweights) + +source: +------- +[JuliaFEM/src/elasticity_solver.jl:68](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L68) + + .. function:: eliminate_boundary_conditions(dirichletbc, I, J, V) + + Eliminate Dirichlet boundary conditions from matrix + + :param dirichletbc : array [dim x nnodes] + :param I, J, V : sparse matrix arrays + :returns: I, J, V : boundary conditions removed + :notes: pros: + - matrix assembly remains positive definite + cons: + - maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) + - if u != 0 in dirichlet boundary requires extra care + :raises: Exception, if displacement boundary conditions given, i.e. + DX=2 for some node, for example. +source: +------- +[JuliaFEM/src/elasticity_solver.jl:218](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L218) + + .. function:: eliminate_boundary_conditions(dirichletbc, I, V) + + Eliminate Dirichlet boundary conditions from vector + + :param dirichletbc : array [dim x nnodes] + :param I, V : sparse vector arrays + :returns: I, V : boundary conditions removed + :notes: pros: + - matrix assembly remains positive definite + cons: + - maybe inefficient because of extra sparse matrix operations. (It's hard to remove stuff from sparse matrix.) + - if u != 0 in dirichlet boundary requires extra care + :raises: Exception, if displacement boundary conditions given, i.e. + DX=2 for some node, for example. +source: +------- +[JuliaFEM/src/elasticity_solver.jl:257](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L257) + + .. function:: interpolate{T<:Real}(field::Array{T<:Real, 1}, basis::Function, ip) + + Interpolate field variable using basis functions f for point ip. + This function tries to be as general as possible and allows interpolating + lot of different fields. + + :param field : Array{Number, dim} + :param basis : Function + :param ip : Array{Number, 1} +source: +------- +[JuliaFEM/src/elasticity_solver.jl:30](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L30) + + .. function:: solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc, lambda, mu, N, dNdchi, ipoints, iweights) + +source: +------- +[JuliaFEM/src/elasticity_solver.jl:278](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/elasticity_solver.jl#L278) + diff --git a/docs/api/JuliaFEM.interfaces.rst b/docs/api/JuliaFEM.interfaces.rst new file mode 100644 index 0000000..f6c7d4b --- /dev/null +++ b/docs/api/JuliaFEM.interfaces.rst @@ -0,0 +1,16 @@ +JuliaFEM.interfaces +=================== + +Exported +-------- + + .. function:: solve_elasticity_interface!() + + This is generic interface that reads data from data model, solves elasticity + problem and updates model. + + :param model : to be defined +source: +------- +[JuliaFEM/src/interfaces.jl:22](https://github.com/JuliaFEM/JuliaFEM.jl/tree/08077a783a6459336cf07b02e323df0e9977c74f/src/interfaces.jl#L22) + diff --git a/docs/api/JuliaFEM.rst b/docs/api/JuliaFEM.rst new file mode 100644 index 0000000..1d6d74a --- /dev/null +++ b/docs/api/JuliaFEM.rst @@ -0,0 +1,3 @@ +JuliaFEM +======== + diff --git a/docs/api/JuliaFEM.xdmf.rst b/docs/api/JuliaFEM.xdmf.rst new file mode 100644 index 0000000..79f91c8 --- /dev/null +++ b/docs/api/JuliaFEM.xdmf.rst @@ -0,0 +1,3 @@ +JuliaFEM.xdmf +============= + diff --git a/docs/build_api.jl b/docs/build_api.jl index e60f10c..aac615e 100644 --- a/docs/build_api.jl +++ b/docs/build_api.jl @@ -41,6 +41,7 @@ const modules = append_list # run(`cp $main_folder/$file_ $this_folder`) cd(dirname(@__FILE__)) do +<<<<<<< HEAD # Run the doctests *before* we start to generate *any* documentation. # for m in modules # failures = failed(doctest(m)) @@ -50,11 +51,14 @@ cd(dirname(@__FILE__)) do # exit(1) # Bail when doctests fail. # end # end +======= +>>>>>>> origin/master # Generate and save the contents of docstrings as markdown files. index = Index() for mod in modules Lexicon.update!(index, save(joinpath(api_directory, "$(mod).rst"), mod)) end +<<<<<<< HEAD #save(joinpath(api_directory, "index.rst")) # Add a reminder not to edit the generated files. @@ -64,35 +68,11 @@ cd(dirname(@__FILE__)) do # all changes to the originating docstrings/files rather than these ones. # """) # end +======= + # save(joinpath(api_directory, "index.rst"), index; md_subheader = :category) +>>>>>>> origin/master # info("Adding all documentation changes in $(api_directory) to this commit.") # success(`git add $(api_directory)`) || exit(1) end - -# point of this? -#cd(dirname(dirname(@__FILE__))) do -# yaml = """ -# # This is automatically generated by docs/build.jl. Edit that file if you -# # wish to make any permenant changes. -# site_name: JuliaFEM.jl -# site_description: JuliaFEM.jl, open-source software for reliable, scalable, distributed Finite Element Method. -# repo_name: GitHub -# # docs_dir: 'docs' -# # site_dir: 'site' -# repo_url: https://github.com/JuliaFEM/JuliaFEM.jl -# pages: -# - Home: 'index.md' -# #- API Docs: -# - JuliaFEM: 'api/JuliaFEM.md' -# - JuliaFEM.elasticity_solver: 'api/JuliaFEM.elasticity_solver.md' -# theme: readthedocs -# """ -# -# # TODO: add the solutions if I figure out how to get them to render properly -# -# open("mkdocs.yml", "w") do f -# write(f, yaml) -# end -#end - diff --git a/docs/source/modules.rst b/docs/source/modules.rst new file mode 100644 index 0000000..7033b51 --- /dev/null +++ b/docs/source/modules.rst @@ -0,0 +1,6 @@ +src +=== + +.. toctree:: + :maxdepth: 4 + diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb index b1b8047..baccd33 100644 --- a/notebooks/2015-06-25-elasticity-solver-example.ipynb +++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb @@ -27,7 +27,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 1, "metadata": { "collapsed": false }, @@ -38,7 +38,7 @@ "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" ] }, - "execution_count": 8, + "execution_count": 1, "metadata": {}, "output_type": "execute_result" } @@ -63,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 2, "metadata": { "collapsed": false }, @@ -74,7 +74,7 @@ "calc_local_matrices! (generic function with 1 method)" ] }, - "execution_count": 9, + "execution_count": 2, "metadata": {}, "output_type": "execute_result" } @@ -159,12 +159,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "*Design principle 6*: we test our code. We use FactCheck in testing." + "*Design principle 6*: we test our code. We use FactCheck for testing." ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 3, "metadata": { "collapsed": true }, @@ -175,7 +175,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -191,16 +191,23 @@ "name": "stderr", "output_type": "stream", "text": [ - "28-Jul 22:15:35:DEBUG:root:Converged in 6 iterations.\n", - "28-Jul 22:15:35:DEBUG:root:solution vector: \n", + "29-Jul 23:12:24:DEBUG:root:Converged in 6 iterations.\n", + "29-Jul 23:12:24:DEBUG:root:solution vector: \n", " [0.0 -0.39914506095474317 -0.07228582695592449 0.0\n", " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n", - "28-Jul 22:15:35:DEBUG:root:norm of u: 3.1292483947150043\n", - "28-Jul 22:15:35:DEBUG:root:Converged in 6 iterations.\n", - "28-Jul 22:15:35:DEBUG:root:solution vector: \n", + "29-Jul 23:12:24:DEBUG:root:norm of u: 3.1292483947150043\n", + "29-Jul 23:12:25:DEBUG:root:Converged in 6 iterations.\n", + "29-Jul 23:12:25:DEBUG:root:solution vector: \n", " [0.0 0.7433248532717793 1.0485210147234858 0.0\n", " 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n", - "28-Jul 22:15:35:DEBUG:root:norm of u: 3.129248394715006\n" + "29-Jul 23:12:25:DEBUG:root:norm of u: 3.129248394715006\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2 facts verified.\n" ] }, { @@ -209,7 +216,7 @@ "delayed_handler (generic function with 4 methods)" ] }, - "execution_count": 11, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -289,56 +296,63 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "One element solutions are not particularly interesting so next step is to create function that assembles global matrix from local matrices." + "One element solutions are not particularly interesting so next step is to create function that assembles global matrix from local matrices. Some data types:" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 5, "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "2 facts verified.\n" - ] - } - ], + "outputs": [], "source": [ "type Node\n", " id :: Int\n", - " coordinates :: Array{Float64, 1}\n", - "end\n", - "\n", - "type Element\n", - " id :: Int\n", - " nodes :: Array{Node, 1}\n", - " attributes :: Dict{ASCIIString, Any}\n", - "end\n", - "\n", - "type FunctionSpace\n", - " family :: String\n", - " order :: Int\n", - "end\n", - "\n", - "type IntegrationScheme\n", - " points :: Array{Float64, 2}\n", - " weights :: Array{Float64, 1}\n", - "end\n", - "\n", - "type IJVMatrix\n", - " I :: Array{Int64, 1}\n", - " J :: Array{Int64, 1}\n", - " V :: Array{Float64, 1}\n", + " elements :: Array{Int64, 1}\n", "end" ] }, { "cell_type": "code", - "execution_count": 105, + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "type Element\n", + " id :: Int\n", + " nodes :: Array{Int64, 1}\n", + " coordinates :: Array{Float64, 2}\n", + " attributes :: Dict{ASCIIString, Any}\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "type Assembly\n", + " # LHS\n", + " I :: Array{Int64, 1}\n", + " J :: Array{Int64, 1}\n", + " A :: Array{Float64, 1}\n", + " # RHS\n", + " i :: Array{Int64, 1}\n", + " b :: Array{Float64, 1}\n", + " # global dofs for each element\n", + " gdofs :: Dict{Int64, Array{Int64, 1}}\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 19, "metadata": { "collapsed": false }, @@ -346,17 +360,29 @@ { "data": { "text/plain": [ - "Element2(1,[Node2(1,[1.0,2.0,3.0])],Dict(\"Young's modulus\"=>2.1e11,\"Poisson's ration\"=>0.3),\"Lagrange(1,2)\",\"FPG4\")" + "get_shape_functions (generic function with 1 method)" ] }, - "execution_count": 105, + "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function get_shape_functions(fs::FunctionSpace)\n", - " if (fs.family == \"Lagrange\") & (fs.order == 1)\n", + "\"\"\"\n", + "Return shape functions and their derivatives for a element.\n", + "\n", + "Parameters\n", + "----------\n", + "element::Element\n", + "\n", + "Returns\n", + "-------\n", + "tuple (basis, dbasis)\n", + "\"\"\"\n", + "function get_shape_functions(el::Element)\n", + " ndim, nnodes = size(el.coordinates)\n", + " if (nnodes == 4) & (ndim == 2)\n", " basis(xi) = [\n", " (1-xi[1])*(1-xi[2])/4\n", " (1+xi[1])*(1-xi[2])/4\n", @@ -369,36 +395,12 @@ " return basis, dbasis\n", " end\n", " throw(\"Unknown function space\")\n", - "end\n", - "\n", - "function get_integration_scheme(is::IntegrationScheme)\n", - " if is.\n", - "end\n", - "\n", - "function assemble!(IJV::IJVMatrix, element::Element, gdofs, fs::FunctionSpace, is::IntegrationScheme)\n", - " nnodes = length(element.nodes)\n", - " ndim = length(elements.nodes[1].coordinates)\n", - "\n", - " # Material properties\n", - " E = element.attributes[\"Young\"]\n", - " nu = element.attributes[\"Poisson\"]\n", - " mu = E/(2*(1+nu))\n", - " la = E*nu/((1+nu)*(1-2*nu))\n", - " la = 2*la*mu/(la + 2*mu)\n", - "\n", - " R = zeros(ndim, nnodes)\n", - " K = zeros(ndim*nnodes, ndim*nnodes)\n", - "\n", - " basis, dbasis = get_shape_functions(fs)\n", - "\n", - " calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n", - "\n", "end" ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -406,25 +408,30 @@ { "data": { "text/plain": [ - "IntegrationScheme(4x2 Array{Float64,2}:\n", - " -0.57735 -0.57735\n", - " 0.57735 -0.57735\n", - " 0.57735 0.57735\n", - " -0.57735 0.57735,[1.0,1.0,1.0,1.0])" + "get_integration_scheme (generic function with 2 methods)" ] }, - "execution_count": 13, + "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "is = IntegrationScheme(1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1], [1, 1, 1, 1])\n" + "\"\"\"\n", + "\"\"\"\n", + "function get_integration_scheme(el::Element, order=2)\n", + " ndim, nnodes = size(el.coordinates)\n", + " if (nnodes == 4) & (order == 2)\n", + " ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n", + " iweights = [1, 1, 1, 1]\n", + " return ipoints, iweights\n", + " end\n", + "end" ] }, { "cell_type": "code", - "execution_count": 111, + "execution_count": 48, "metadata": { "collapsed": false }, @@ -432,45 +439,183 @@ { "data": { "text/plain": [ - "Node2(2,[3.0,2.0,3.0])" + "assemble_element! (generic function with 2 methods)" ] }, - "execution_count": 111, + "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "e2 = Element2(2, [n1], Dict(\"Young's modulus\" => 90.0e9, \"Poisson's ration\" => 0.3), \"Lagrange(1,2)\", \"FPG4\")\n", - "e2.nodes[1]" + "function assemble_element!(ass::Assembly, el::Element, io=2)\n", + "\n", + " # Material properties\n", + " E = el.attributes[\"Young\"]\n", + " nu = el.attributes[\"Poisson\"]\n", + " mu = E/(2*(1+nu))\n", + " la = E*nu/((1+nu)*(1-2*nu))\n", + " la = 2*la*mu/(la + 2*mu)\n", + "\n", + " dofs = prod(size(el.coordinates))\n", + " X = el.coordinates\n", + " u = el.attributes[\"displacement\"]\n", + " R = el.attributes[\"displacement nodal force\"]\n", + " K = el.attributes[\"displacement tangent stiffness\"]\n", + "\n", + " gdofs = ass.gdofs[el.id]\n", + " basis, dbasis = get_shape_functions(el)\n", + " ipoints, iweights = get_integration_scheme(el, io)\n", + " calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n", + "\n", + " for i=1:dofs\n", + " for j=1:dofs\n", + " push!(ass.I, gdofs[i])\n", + " push!(ass.J, gdofs[j])\n", + " push!(ass.A, K[i,j])\n", + " end\n", + " push!(ass.i, gdofs[i])\n", + " push!(ass.b, R[i])\n", + " end\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Time to test again. From last test we know that correct solution is\n", + "\n", + " [0.0 -0.39914506095474317 -0.07228582695592449 0.0\n", + " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n", + "\n", + "This time we assemble global stiffness matrix in different order, 2 3 4 1" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 49, "metadata": { - "collapsed": true + "collapsed": false }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "one element assembly\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "30-Jul 00:08:18:DEBUG:root:Creating nodes\n", + "30-Jul 00:08:18:DEBUG:root:Creating elements\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 1\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 3.0814821107320176\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 2\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 0.32007464366194766\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 3\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 0.040279810888447135\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 4\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 0.000925649536063315\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 5\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 1.5582555024825817e-7\n", + "30-Jul 00:08:18:DEBUG:root:Starting iteration 6\n", + "30-Jul 00:08:18:DEBUG:root:Assembling\n", + "30-Jul 00:08:18:DEBUG:root:Solving\n", + "30-Jul 00:08:18:DEBUG:root:Solution norm = 1.0317166868587156e-14\n", + "30-Jul 00:08:18:DEBUG:root:Converged in 6 iterations.\n", + "30-Jul 00:08:18:DEBUG:root:Displacement of element = \n", + "[-0.39914506095474334 -0.07228582695592464 0.0 0.0\n", + " -2.1779892317073504 -2.222244754401764 0.0 0.0]\n" + ] + }, + { + "data": { + "text/plain": [ + "delayed_handler (generic function with 4 methods)" + ] + }, + "execution_count": 49, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 fact verified.\n" + ] + } + ], "source": [ - "n1 = new_node()\n", - "set_node_id(n1, 1)\n", - "set_node_coords(n1, [0.0, 0.0])\n", + "facts(\"one element assembly\") do\n", + " Logging.debug(\"Creating nodes\")\n", + " n1 = Node(1, Int64[])\n", + " n2 = Node(2, Int64[])\n", + " n3 = Node(3, Int64[])\n", + " n4 = Node(4, Int64[])\n", + " nodes = [n1.id, n2.id, n3.id, n4.id]\n", + " coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n", + " attributes = Dict(\n", + " \"Young\" => 90, \"Poisson\" => 0.25,\n", + " \"displacement\" => zeros(2, 4),\n", + " \"displacement nodal force\" => zeros(2, 4),\n", + " \"displacement tangent stiffness\" => zeros(8, 8))\n", + " Logging.debug(\"Creating elements\")\n", + " el = Element(1, nodes, coordinates, attributes)\n", + " for i=1:10\n", + " Logging.debug(\"Starting iteration $i\")\n", + " ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())\n", + " ass.gdofs[el.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n", + " Logging.debug(\"Assembling\")\n", + " assemble_element!(ass, el)\n", + " Logging.debug(\"Solving\")\n", + " du = zeros(2, 4) # must be determined from ass\n", + " F = [0 0; 0 -2; 0 0; 0 0]'\n", + " free_dofs = [1, 2, 3, 4]\n", "\n", - "n2 = new_node()\n", - "set_node_id(n2, 2)\n", - "set_node_coords(n2, [10.0, 0.0])\n", + " # solution\n", + " K = sparse(ass.I, ass.J, ass.A)\n", + " R = full(sparsevec(ass.i, ass.b))\n", + " R = reshape(R, (2, round(Int, length(R)/2)))\n", + " du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n", "\n", - "n3 = new_node()\n", - "set_node_id(n3, 3)\n", - "set_node_coords(n3, [10.0, 1.0])\n", + " Logging.debug(\"Solution norm = $(norm(du))\")\n", + " #Logging.debug(\"Solution increment = \\n$du\")\n", "\n", - "n4 = new_node()\n", - "set_node_id(n4, 4)\n", - "set_node_coords(n4, [0.0, 1.0])\n", - "\n", - "nodes = [n1, n2, n3, n4]\n", - "add_nodes(m, \"NALL\", nodes)" + " # update solution back to elements\n", + " eldu = du[ass.gdofs[el.id]]\n", + " eldu = reshape(eldu, (2, round(Int, length(eldu)/2)))\n", + " el.attributes[\"displacement\"] += eldu\n", + " if norm(du) < 1.0e-9\n", + " Logging.debug(\"Converged in $i iterations.\")\n", + " break\n", + " end\n", + " end\n", + " disp = el.attributes[\"displacement\"]\n", + " Logging.debug(\"Displacement of element = \\n$disp\")\n", + " @fact norm(disp) => roughly(3.1292483947150043)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Seems to be working. But we still need to handle boundary conditions more \"cleverly\" and generalize assembly to several elements (which is not problem)." ] }, {