new commit

This commit is contained in:
ovainola
2015-07-30 20:43:26 +03:00
10 changed files with 402 additions and 141 deletions
+6 -6
View File
@@ -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
+3
View File
@@ -0,0 +1,3 @@
module ReST
end
+12
View File
@@ -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)
+93
View File
@@ -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)
+16
View File
@@ -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)
+3
View File
@@ -0,0 +1,3 @@
JuliaFEM
========
+3
View File
@@ -0,0 +1,3 @@
JuliaFEM.xdmf
=============
+7 -27
View File
@@ -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
+6
View File
@@ -0,0 +1,6 @@
src
===
.. toctree::
:maxdepth: 4
@@ -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)."
]
},
{