From fff0173c2d9c637e58acc2ce71b9b48aa572a0f3 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 17 Dec 2015 15:33:51 +0200 Subject: [PATCH] mortar tests, bug finding --- REQUIRE | 8 +- notebooks/2015-12-13-3d-tie-contact.ipynb | 604 +++++++++++----------- src/JuliaFEM.jl | 6 + src/common.jl | 52 ++ src/core.jl | 25 +- src/directsolver.jl | 1 + src/dirichlet.jl | 7 +- src/elements.jl | 10 + src/integrate.jl | 13 + src/linear_elasticity.jl | 17 +- src/mortar.jl | 186 +++++-- src/preprocess_aster_reader.jl | 140 +++++ src/problems.jl | 3 + src/test.jl | 5 +- test/test_mortar.jl | 240 ++++++++- test/test_preprocess_aster_reader.jl | 47 +- 16 files changed, 993 insertions(+), 371 deletions(-) create mode 100644 src/common.jl diff --git a/REQUIRE b/REQUIRE index ee621de..138afda 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,9 +1,7 @@ julia 0.4 +BaseTestNext ForwardDiff -Logging LightXML -Docile -Lexicon -FactCheck HDF5 -Lint +JLD + diff --git a/notebooks/2015-12-13-3d-tie-contact.ipynb b/notebooks/2015-12-13-3d-tie-contact.ipynb index 06ca4b4..94be2fd 100644 --- a/notebooks/2015-12-13-3d-tie-contact.ipynb +++ b/notebooks/2015-12-13-3d-tie-contact.ipynb @@ -360,16 +360,15 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO: Found 4 element sets: LSYM23, LSYM12, LSYM13, LOWER_BLOCK_TO_UPPER_BLOCK\n", - "INFO: Found 4 element sets: USYM23, UPPER_BLOCK_TO_LOWER_BLOCK, USYM13, LOAD\n" + "INFO: Found 6 element sets: SYM23, SYM12, UPPER_TO_LOWER, LOAD, SYM13, LOWER_TO_UPPER\n" ] }, { "data": { "text/plain": [ "Dict{ASCIIString,Any} with 2 entries:\n", - " \"nodes\" => Dict(2=>[0.0,0.0,1.0],11=>[0.0,0.6666666666666666,1.0],39=>…\n", - " \"connectivity\" => Dict(68=>(:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[26,4,14,42]),2=…" + " \"nodes\" => Dict(68=>[1.0,0.0,0.5],2=>[0.0,0.0,0.16666666666666666],89=…\n", + " \"connectivity\" => Dict(11=>(:SE2,:OTHER,[4,13]),158=>(:QU4,:UPPER_TO_LOWER,[8…" ] }, "execution_count": 1, @@ -380,9 +379,15 @@ "source": [ "using JuliaFEM\n", "using JuliaFEM.Preprocess: parse_aster_med_file\n", - "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Quad4, Hex8, LinearSolver\n", - "mesh_lower = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/LOWER_BLOCK.med\")\n", - "mesh_upper = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/UPPER_BLOCK.med\")" + "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity,\n", + " Quad4, Hex8, LinearSolver, update\n", + "#mesh_lower = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/LOWER_BLOCK.med\")\n", + "#mesh_upper = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/UPPER_BLOCK.med\")\n", + "#mesh_lower = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCK_TWOELEM.med\", \"LOWER_BLOCK_SPARSE\")\n", + "#mesh_upper = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCK_TWOELEM.med\", \"UPPER_BLOCK_SPARSE\")\n", + "#mesh1 = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCK_SPARSE.med\", \"LOWER_BLOCK_SPARSE\")\n", + "#mesh2 = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCK_SPARSE.med\", \"UPPER_BLOCK_SPARSE\")\n", + "mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/BLOCKS.med\")" ] }, { @@ -393,19 +398,36 @@ }, "outputs": [ { - "data": { - "text/plain": [ - "1x75 Array{Int64,2}:\n", - " 1 2 3 4 5 6 7 8 9 10 11 12 … 67 68 69 70 71 72 73 74 75" - ] - }, - "execution_count": 2, - "metadata": {}, - "output_type": "execute_result" + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: created 39 elements.\n" + ] } ], "source": [ - "sort(collect(keys(mesh_lower[\"nodes\"])))'" + "# interior elements are of type HE8\n", + "field_problem = LinearElasticityProblem()\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " eltype == :HE8 || continue\n", + " element = Hex8(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " element[\"youngs modulus\"] = 900.0\n", + " element[\"poissons ratio\"] = 0.25\n", + " #info(\"lower: add element with connectivity $elcon\")\n", + " push!(field_problem, element)\n", + "end\n", + "\n", + "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == :LOAD) || continue\n", + " element = Quad4(elcon)\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " element[\"displacement traction force 3\"] = -100.0\n", + " push!(field_problem, element)\n", + "end\n", + "info(\"created $(length(field_problem.elements)) elements.\")" ] }, { @@ -414,232 +436,12 @@ "metadata": { "collapsed": false }, - "outputs": [ - { - "data": { - "text/plain": [ - "1x48 Array{Int64,2}:\n", - " 1 2 3 4 5 6 7 8 9 10 11 12 … 40 41 42 43 44 45 46 47 48" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "sort(collect(keys(mesh_upper[\"nodes\"])))'" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Need to renumber nodes." - ] - }, - { - "cell_type": "code", - "execution_count": 4, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "Dict{Int64,Int64} with 48 entries:\n", - " 2 => 77\n", - " 11 => 86\n", - " 39 => 114\n", - " 46 => 121\n", - " 25 => 100\n", - " 42 => 117\n", - " 29 => 104\n", - " 8 => 83\n", - " 20 => 95\n", - " 14 => 89\n", - " 31 => 106\n", - " 33 => 108\n", - " 18 => 93\n", - " 26 => 101\n", - " 35 => 110\n", - " 17 => 92\n", - " 44 => 119\n", - " 4 => 79\n", - " 37 => 112\n", - " 45 => 120\n", - " 13 => 88\n", - " 30 => 105\n", - " 1 => 76\n", - " 47 => 122\n", - " 32 => 107\n", - " ⋮ => ⋮" - ] - }, - "execution_count": 4, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "new_node_numbering = Dict{Int64, Int64}()\n", - "upper_node_ids = sort(collect(keys(mesh_upper[\"nodes\"])))\n", - "for (k, old_node_id) in enumerate(upper_node_ids)\n", - " new_node_numbering[old_node_id] = 75+k\n", - "end\n", - "new_node_numbering" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "1x48 Array{Int64,2}:\n", - " 76 77 78 79 80 81 82 83 84 85 … 117 118 119 120 121 122 123" - ] - }, - "execution_count": 5, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function renumber!(mesh, node_numbering)\n", - " old_nodes = mesh[\"nodes\"]\n", - " new_nodes = typeof(old_nodes)()\n", - " for (node_id, node_coords) in old_nodes\n", - " new_node_id = node_numbering[node_id]\n", - " new_nodes[new_node_id] = node_coords\n", - " end\n", - " mesh[\"nodes\"] = new_nodes\n", - " for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", - " new_elcon = [node_numbering[node_id] for node_id in elcon]\n", - " mesh[\"connectivity\"][elid] = (eltype, elset, new_elcon)\n", - " end\n", - " return mesh\n", - "end\n", - "\n", - "renumber!(mesh_upper, new_node_numbering)\n", - "sort(collect(keys(mesh_upper[\"nodes\"])))'" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "Dict{Int64,Tuple{Symbol,Symbol,Array{Int64,1}}} with 92 entries:\n", - " 68 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[101,79,89,117])\n", - " 2 => (:SE2,:OTHER,[77,84])\n", - " 89 => (:HE8,:OTHER,[114,102,103,115,121,118,119,123])\n", - " 11 => (:SE2,:OTHER,[90,80])\n", - " 39 => (:QU4,:USYM13,[93,106,90,80])\n", - " 46 => (:QU4,:LOAD,[96,109,108,95])\n", - " 85 => (:HE8,:OTHER,[94,93,114,116,107,106,121,120])\n", - " 25 => (:SE2,:OTHER,[99,83])\n", - " 55 => (:QU4,:OTHER,[95,113,112,96])\n", - " 42 => (:QU4,:USYM13,[106,92,81,90])\n", - " 66 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[116,88,76,94])\n", - " 58 => (:QU4,:OTHER,[113,101,100,112])\n", - " 29 => (:SE2,:OTHER,[79,101])\n", - " 59 => (:QU4,:OTHER,[87,79,101,113])\n", - " 8 => (:SE2,:OTHER,[88,76])\n", - " 74 => (:QU4,:OTHER,[99,83,103,119])\n", - " 90 => (:HE8,:OTHER,[89,117,101,79,105,122,113,87])\n", - " 57 => (:QU4,:OTHER,[112,100,83,99])\n", - " 20 => (:SE2,:OTHER,[96,95])\n", - " 78 => (:HE8,:OTHER,[104,120,122,105,85,110,108,86])\n", - " 14 => (:SE2,:OTHER,[92,91])\n", - " 31 => (:SE2,:OTHER,[103,102])\n", - " 70 => (:QU4,:OTHER,[98,119,118,97])\n", - " 33 => (:QU4,:USYM23,[88,104,84,76])\n", - " 52 => (:QU4,:LOAD,[111,92,91,110])\n", - " ⋮ => ⋮" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "mesh_upper[\"connectivity\"]" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": false - }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "INFO: created 59 elements.\n" - ] - } - ], - "source": [ - "# interior elements are of type HE8\n", - "field_problem = LinearElasticityProblem()\n", - "\n", - "for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n", - " eltype == :HE8 || continue\n", - " element = Hex8(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n", - " element[\"youngs modulus\"] = 900.0\n", - " element[\"poissons ratio\"] = 0.25\n", - " #info(\"lower: add element with connectivity $elcon\")\n", - " push!(field_problem, element)\n", - "end\n", - "\n", - "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", - " eltype == :HE8 || continue\n", - " element = Hex8(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", - " element[\"youngs modulus\"] = 900.0\n", - " element[\"poissons ratio\"] = 0.25\n", - " #info(\"upper: add element with connectivity $elcon\")\n", - " push!(field_problem, element)\n", - "end\n", - "\n", - "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n", - "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", - " (eltype == :QU4) && (elset == :LOAD) || continue\n", - " element = Quad4(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", - " element[\"displacement traction force\"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]\n", - " push!(field_problem, element)\n", - "end\n", - "info(\"created $(length(field_problem.elements)) elements.\")" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "INFO: created 44 boundary elements.\n" + "INFO: created 35 boundary elements.\n" ] } ], @@ -647,34 +449,19 @@ "# boundary conditions\n", "boundary_problem = DirichletProblem(\"displacement\", 3)\n", "\n", - "for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n", - " (eltype == :QU4) && (elset in [:LSYM23, :LSYM12, :LSYM13]) || continue\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n", " element = Quad4(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n", - " if elset == :LSYM23\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " if elset == :SYM23\n", " element[\"displacement 1\"] = 0.0\n", - " elseif elset == :LSYM12\n", + " elseif elset == :SYM12\n", " element[\"displacement 3\"] = 0.0\n", - " elseif elset == :LSYM13\n", + " elseif elset == :SYM13\n", " element[\"displacement 2\"] = 0.0\n", " end\n", " push!(boundary_problem, element)\n", "end\n", - "\n", - "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", - " (eltype == :QU4) && (elset in [:USYM23, :USYM12, :USYM13]) || continue\n", - " element = Quad4(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", - " if elset == :USYM23\n", - " element[\"displacement 1\"] = 0.0\n", - " elseif elset == :USYM12\n", - " element[\"displacement 3\"] = 0.0\n", - " elseif elset == :USYM13\n", - " element[\"displacement 2\"] = 0.0\n", - " end\n", - " push!(boundary_problem, element)\n", - "end\n", - "\n", "info(\"created $(length(boundary_problem.elements)) boundary elements.\")" ] }, @@ -689,39 +476,94 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 4, "metadata": { "collapsed": false }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "WARNING: using PyPlot.mesh in module Main conflicts with an existing identifier.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAG7CAYAAABEoAutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9UVXW+//HXhoOgCPRDM6hRFCVNC4U0/HGvP5prjqWN5Y+Y4GaatPo2y8ZmFHWmaY02/WC8F8euTVoT/oyxS5O1XDXOclplkmEDemcVeM2u4JQYVqZYQICf7x8qM8g5eA4/zjl8eD7WclGfvT+c92e/Yb/Y5+wDjjHGCAAAC4QEugAAANoLoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALCGK9AFtOSLL77Qzp07FR8fr+7duwe6HABAAFRXV6usrEy33nqrevXq1eK+QR1qO3fuVHp6eqDLAAAEgS1btuiee+5pcZ+gDrX4+HhJ5xYyZMiQgNSwaNEi5eTkBOSxg0VXPwZdff0Sx6Crr18K7DEoLS1Venp6Yya0JKhD7cJTjkOGDFFycnJAaoiJiQnYYweLrn4Muvr6JY5BV1+/FBzHwJuXobhRBABgDUINAGANQg0AYA1C7RLS0tICXULAdfVj0NXXL3EMuvr6pc5zDJxg/ntqxcXFSklJUVFRUcBfoAQABIYvWcCVGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBo+hdo333yjxx57TFOmTNEVV1yhkJAQbdy40ev5X3/9tTIzM9W7d2/17NlTkyZN0v79+30uGgAAd3wKtRMnTmjlypX63//9Xw0fPlyS5DiOV3PPnj2r2267TXl5eVq4cKGys7NVWVmpCRMm6PDhw75XDgDARVy+7BwXF6fjx4/rqquuUlFRkUaOHOn13Pz8fO3du1f5+fm68847JUmzZ89WYmKiHnvsMW3dutW3ygEAuIhPV2rdunXTVVddJUkyxvj0QPn5+br66qsbA02SevXqpdmzZ+u1115TXV2dT58PAICL+e1Gkf379ys5ObnZ+MiRI/Xtt9/q0KFD/irlkk6cOKH77lusoUNv03XXTdfQobfpvvsW68SJE4EuDW4UFBQoOnqYHGewHCdJjjNY0dHDVFBQEOjScBF61bmUlpZq0KAJCg+/Qd26pSg8/AYNGjRBpaWlgS7NI5+efmyLiooKTZgwodl4bGysJOnYsWMaOnSov8rxqLKyUmPG3K1PPnlCUrYkR9JZlZTs07vvztHevdvUu3fvAFeJC3bv3q3x4+dJ2iLpZl3oV1VVocaNy9CePRs1duzYwBYJSfSqs/noo480fPhdqq/foH/u1+HDhbrxxhn6299e1ZAhQwJbpBt+u1KrqalReHh4s/GIiAhJUnV1tb9KaVFW1m/OB1qqzjVROneYUvXJJ7/WkiXZgSsOzdx++/+TtFnN+zVa0ib94AcPBKo0XIRedS4//OFDqq/Plbt+1de/qOnTHwxccS3w25Va9+7dVVtb22y8pqamcXsw2LevROeu0Ny5We+887iKi/1ZEVpSVVWvc9907qSqqiqUfgWJqqoGtdyrenoVRMrLv1RL/Tp69Et/luM1v4VabGysjh071my8oqJC0rk7Kz1ZtGiRYmJimoylpaUpLS2tfYuUVF8fqn/8VHKxEB05EqqUlHZ/WLRauFrql9SPfgWNHmq5V+H0Kqh0U0v9MqZbhzxqXl6e8vLymoydOnXK6/l+C7Xhw4fr3XfflTGmyXvbCgsLFRkZqcTERI9zc3Jy3N5k0hFcrgZJRu6beVb9+zcoP98vpcALKSm1aqlfUrmKivxbE9xLSflWLfeqll4FkdTU71RX57lfjvNdhzyuuwuW4uJipXj5E0+HhNrx48f19ddfa+DAgXK5zj3EzJkzlZ+frz/+8Y+66667JElffPGF/vu//1vTpk1TWFhYR5Tis1GjrldJSaHcX3YXavz46+WnfIUXoqJcqqp6X+del7nY+4qKaqBfQSIqKvQSvXLRqyDSr9+VOnzYc7/69r3S3yV5x/jomWeeMStXrjQPPvigcRzH3HXXXWblypVm5cqV5tSpU8YYY+69917jOI4pLy9vnNfQ0GBGjx5toqKizIoVK8zatWvN0KFDTUxMjDl06JDbxyoqKjKSTFFRka9ltlplZaVJSJhopPeM1GAkc/7jeyYhYaKprKz0Wy24tD179hgpwUgFF/WrwEgJZs+ePYEuEefRq86lpKTEuFzXue2Xy3WdKSkp8VstvmSBz6EWHx9vHMcxjuOYkJAQExIS0vjfF0Js7ty5Tf7/gpMnT5r777/f9OrVy0RGRpqJEye2WGQgQs2Yc8E2d+7PTP/+U400zfTvP9XMnfszAi1I7dmzx0RFDTXSjUaaZqQbTVTUUE6SQegfvRpspGQjDaZXQaykpMQMHDjehIUNM1KyCQsbZgYOHO/XQDOmg0PNnwIVav94fGOkcx8R/LZsOdevLVsCXQkuhV51LoE+F/qSBfzpGQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1fAq12tpaZWVlKS4uTj169FBqaqp27drl1dyioiLdfvvtio2NVVRUlJKSkvTMM8/o7NmzrSocAICL+RRqc+fOVU5OjjIyMrRmzRqFhoZq6tSpKigoaHFeUVGRxowZo6NHj2rp0qX6z//8Tw0YMEAPP/ywHnnkkTYtAACAC1ze7rhv3z5t27ZNq1atagyijIwMDRs2TEuWLGkx2NatW6eQkBDt3r1bl112mSRpwYIFmjBhgjZs2KDVq1e3cRkAAPhwpZafny+Xy6XMzMzGsfDwcM2fP1979+7VZ5995nHu6dOnFR4erpiYmCbjV199tXr06NGKsgEAaM7rUNu/f78SExPVs2fPJuMjR46UJB04cMDj3IkTJ+r06dN64IEHdPDgQZWXl+u5557Tq6++qmXLlrWydAAAmvL66ceKigrFxsY2G78wduzYMY9zFyxYoI8++kjr1q3TCy+8IEkKDQ3V2rVrm1z5AQDQFl6HWnV1tcLDw5uNR0RENG73JCQkRAMGDNCUKVM0a9YsRURE6KWXXtKPf/xj9enTR3fccUcrSgcAoCmvQ6179+6qra1tNl5TU9O43ZOnnnpKa9as0eHDhxtfQ5s5c6YmTZqkhx56SLfffrtCQ0M9zl+0aFGz1+PS0tKUlpbmbfkAgE4gLy9PeXl5TcZOnTrl9XyvQy02NtbtU4wVFRWSpLi4OI9zn332Wd1yyy3NbgqZNm2afvrTn6q8vFwDBgzwOD8nJ0fJycnelgoA6KTcXbAUFxcrJSXFq/le3ygyYsQIHTp0SFVVVU3GCwsLJUnDhw/3OLeyslINDQ3Nxuvq6iRJ9fX13pYBAIBHXofazJkz1dDQoPXr1zeO1dbWKjc3V6mpqbrmmmskScePH9fBgwebBFViYqL+/Oc/66uvvmoca2ho0Msvv6zo6GglJCS0x1oAAF2c108/jho1SrNmzdKyZctUWVmphIQEbdy4UUePHlVubm7jfkuXLtWmTZtUVlamvn37No6lp6fr5ptvVmZmpiIiIpSXl6fi4mL9+te/bvH1NAAAvOV1qEnSpk2b9Oijj2rz5s06efKkkpKStGPHDo0bN65xH8dx5DhOk3k/+tGP1KtXLz355JP6zW9+o9OnT2vw4MFat26dFixY0D4rAQB0eT6FWnh4uLKzs5Wdne1xn9zc3CZXbhdMnjxZkydP9r1CAAC8xJ+eAQBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1Bzo7S0VIMGTVBq6g2SUpSaeoMGDZqg0tLSQJcGN7Zu3aqQkAFKT79O0o1KT79OISEDtHXr1kCXhos07VUSvQpyBQUFio4eppSUJEnTlZKSpOjoYSooKAh0aZ6ZIFZUVGQkmaKiIr895ocffmhcruuMtNdIZ41kjNRgpPeMy3WdKSkp8VstuLQNGzYYaYDbfkkDzJYtWwJdIs6jV53LO++8Y6QED/1KMHv27PFbLb5kgWOMMYEOVk+Ki4uVkpKioqIiJScn++UxBw2aoMOHn5Q02s3W9zRw4HJ9/PHbfqkFlxYSMkDGbJWnfjlOus6e/T9/lwU36FXnEh09TFVVz8tTv6KiMnX69Id+qcWXLHD5paJO5OjRLyWletiaqrKy0+KZkuBhTDe11C9jwuhXkKBXnUtVVaha6ldVVb0/y/EaoXaRc994joetIaqvv1bp6f6sCC2LVEv9kiLoV9C4VK/C6VVQ6adL9SsYEWoXcZzvJBm5b+ZZuVyfasMG/9YEz9LTv1FL/ZJqtGWLf2uCe5fuVS29CiLp6eW6VL+CEaF2kb59r9Thw+/L/fPI7ys+Plr33OPvquBJRsZ3MsZzvxynjn4FCXrVuTz4YIOqqjz3KyoqOOODW/ov8vrrv5PLdZ+k93TupxGd//ieXK55ev313wWuODSzefNKSely1y8p4/x2BAN61bm8+eY6SRly369/P789+ARn1AbQkCFD9Le/varp0x9Ueflp1dVdq7CwT9WvX7Ref/1VDRkyJNAl4p/cc/5H+4yM9POvh0ZK+kaO8502b17ZuB2BR686l7Fjx2rPno36wQ8yz98UEi6pVlFRLr355kaNHTs20CW61+FvMGiDQLxPrenjGyOd+4jgt2XLuX7xdqfgR686l0CfC33JAp5+BABYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWMOnUKutrVVWVpbi4uLUo0cPpaamateuXV7P37VrlyZNmqTLLrtM0dHRuummm/Tyyy/7XDQAAO74FGpz585VTk6OMjIytGbNGoWGhmrq1KkqKCi45Nzc3FzdeuutCg8P15NPPqlVq1bpX//1X/Xpp5+2ungAAP6Z17/QeN++fdq2bZtWrVqlRx55RJKUkZGhYcOGacmSJS0GW1lZmR566CEtXLhQOTk5ba8aAAA3vL5Sy8/Pl8vlUmZmZuNYeHi45s+fr7179+qzzz7zOPe5556TMUYrVqyQJJ05c0bGmDaUDQBAc16H2v79+5WYmKiePXs2GR85cqQk6cCBAx7n7tq1S4MHD9aOHTt07bXXKjo6Wr169dIvf/lLwg0A0G68fvqxoqJCsbGxzcYvjB07dszj3I8//lgul0vz5s1TVlaWkpKS9Morr+jxxx9XfX29nnjiiVaUDgBAU16HWnV1tcLDw5uNR0RENG735MLTjU8//bQWL14sSZoxY4a++uor/fa3v9Xy5cubXQECAOArr0Ote/fuqq2tbTZeU1PTuL2ludXV1UpLS2syfvfdd+tPf/qTDhw4oHHjxnmcv2jRIsXExDQZS0tLa/b5AACdW15envLy8pqMnTp1yuv5XodabGys26cYKyoqJElxcXEe58bFxemTTz5Rnz59moxfddVVkqSTJ0+2+Ng5OTlKTk72tlQAQCfl7oKluLhYKSkpXs33+kaRESNG6NChQ6qqqmoyXlhYKEkaPny4x7k33XSTjDHN3pN2ISR79+7tbRkAAHjkdajNnDlTDQ0NWr9+feNYbW2tcnNzlZqaqmuuuUaSdPz4cR08eFD19fWN+82ZM0eS9Pvf/75x7OzZs8rNzdWVV17pdQIDANASr59+HDVqlGbNmqVly5apsrJSCQkJ2rhxo44eParc3NzG/ZYuXapNmzaprKxMffv2lSTdcccduuWWW/Tkk0/qiy++0I033qjt27eroKBA69evV1hYWPuvDADQ5XgdapK0adMmPfroo9q8ebNOnjyppKQk7dixo8lNHo7jyHGcZnO3b9+uX/ziF9q2bZs2bNigwYMHa+vWrdzsAQBoNz6FWnh4uLKzs5Wdne1xn9zc3CZXbhdERkYqJyeHX5MFAOgw/OkZAIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1fAq12tpaZWVlKS4uTj169FBqaqp27drl84MuWLBAISEhmjZtms9zAQDwxKdQmzt3rnJycpSRkaE1a9YoNDRUU6dOVUFBgdef469//as2btyoiIgIOY7jc8EAAHjidajt27dP27Zt01NPPaWnn35a999/v9566y3169dPS5Ys8epzGGO0cOFC3XvvverTp0+riwYAwB2vQy0/P18ul0uZmZmNY+Hh4Zo/f7727t2rzz777JKfY/PmzSopKdHjjz8uY0zrKgYAwAOvQ23//v1KTExUz549m4yPHDlSknTgwIEW51dVVSkrK0vLly/nKg0A0CG8DrWKigrFxsY2G78wduzYsRbnr1ixQpGRkVq0aJGPJQIA4B2XtztWV1crPDy82XhERETjdk8OHTqkNWvW6A9/+IPCwsJaUSYAAJfm9ZVa9+7dVVtb22y8pqamcbsnDz/8sMaOHasZM2a0okQAALzj9ZVabGys26cYKyoqJElxcXFu57311lvauXOn/vjHP6qsrKxxvL6+Xt9++63Ky8t1xRVXKCoqyuNjL1q0SDExMU3G0tLSlJaW5m35AIBOIC8vT3l5eU3GTp065fV8r0NtxIgRevvtt1VVVdUkgAoLCyVJw4cPdzvv6NGjkqQ777yz2bZjx46pf//+Wr16tRYuXOjxsXNycpScnOxtqQCATsrdBUtxcbFSUlK8mu91qM2cOVOrVq3S+vXr9dOf/lTSud8wkpubq9TUVF1zzTWSpOPHj+vrr7/WwIED5XK5dMstt2j79u1NPpcxRpmZmYqPj9fPf/5zDRs2zNsyAADwyOtQGzVqlGbNmqVly5apsrJSCQkJ2rhxo44eParc3NzG/ZYuXapNmzaprKxMffv21fe+9z1973vfa/b5Hn74YfXp00fTp09vn5UAALo8r0NNkjZt2qRHH31Umzdv1smTJ5WUlKQdO3Zo3Lhxjfs4juPVr7/iV2QBANqbT6EWHh6u7OxsZWdne9wnNze3yZWbJ0eOHPHloQEAuCT+9AwAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahJobpaWlGjRoglJTkyVNV2pqsgYNmqDS0tJAlwY3tm7dqpCQAUpPT5I0XenpSQoJGaCtW7cGujRc5B+9GiwpRenpg+lVECsoKFB09DClpAyWlKSUlMGKjh6mgoKCQJfmmQliRUVFRpIpKiry22N++OGHxuW6zkh7jXTWSMZIDUZ6z7hc15mSkhK/1YJL27Bhg5EGuO2XNMBs2bIl0CXiPHrVubzzzjtGSvDQrwSzZ88ev9XiSxY4xhgT6GD15MKf8C4qKlJycrJfHnPQoAk6fPhJSaPdbH1PAwcu18cfv+2XWnBpISEDZMxWeeqX46Tr7Nn/83dZcINedS7R0cNUVfW8PPUrKipTp09/6JdafMkCn/6eWldw9OiXklI9bE1VWdlp8UxJ8DCmp1rqlzHd6FeQMKabWu5VT3oVRKqqGtRSv6qq6v1ZjtcItYuc+8bz9Fe5Q1Rff63S0/1ZEVoWr5b6JUXSr6ARqZZ7FU+vgkoPtdyvcD/W4j1C7SKO850kI/fNPCuX61Nt2ODfmuBZenqZWuqX9I22bPFrSfAgPf0btdyrMnoVRNLTv1XL/ar1b0FeItQu0rfvlTp8+H25fx75fcXHR+uee/xdFTzJyDgjYzz3y3G+o19BIiPju0v06gy9CiIPPhiqqirP/YqKCs744Jb+i7z++u/kct0n6T2d+2lE5z++J5drnl5//XeBKw7NbN68UlK63PVLyji/HcGAXnUub765TlKG3Pfr389vDz7BGbUBNGTIEP3tb69q+vQHVV7+perquiks7Dv163elXn/9VQ0ZMiTQJeKf3HP+R/uMjHQZEyYpQlKNHKdOmzevbNyOwGvaq54693pomRznDL0KQmPHjtWePRv1gx9knr9ppIekbxUVFao339yosWPHBrpE9zr8DQZtEIj3qTV9fGOkcx8R/LZsOdcv3u4U/OhV5xLoc6EvWcDTjwAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAaxBqAABrEGoAAGsQagAAa/gcarW1tcrKylJcXJx69Oih1NRU7dq165Lz/vKXv2jevHlKTExUZGSkEhIStGDBAh0/frxVhQMAcDGfQ23u3LnKyclRRkaG1qxZo9DQUE2dOlUFBQUtzsvKytLu3bt111136ZlnntHdd9+tl19+WSNGjNDnn3/e6gUAAHCBy5ed9+3bp23btmnVqlV65JFHJEkZGRkaNmyYlixZ0mKwrV69WuPGjWsyNmXKFI0fP17/9V//pZUrV7aifAAA/sGnK7X8/Hy5XC5lZmY2joWHh2v+/Pnau3evPvvsM49zLw40SfqXf/kXXXHFFTp48KAvZQAA4JZPobZ//34lJiaqZ8+eTcZHjhwpSTpw4IBPD37mzBlVVVWpV69ePs0DAMAdn0KtoqJCsbGxzcYvjB07dsynB1+9erXq6uo0Z84cn+YBAOCOT6FWXV2t8PDwZuMRERGN2721e/du/epXv9KcOXM0YcIEX8oAAMAtn24U6d69u2pra5uN19TUNG73xsGDBzVjxgzdeOONeuGFFy65/6JFixQTE9NkLC0tTWlpaV49HgCgc8jLy1NeXl6TsVOnTnk936dQi42NdfsUY0VFhSQpLi7ukp/j73//uyZPnqzLL79cb7zxhiIjIy85JycnR8nJyb6UCgDohNxdsBQXFyslJcWr+T49/ThixAgdOnRIVVVVTcYLCwslScOHD29x/pdffqnJkyerrq5OO3fuVJ8+fXx5eAAAWuRTqM2cOVMNDQ1av35941htba1yc3OVmpqqa665RpJ0/PhxHTx4UPX19Y37ffPNN5o6daoqKir0xhtvKCEhoZ2WAADAOT49/Thq1CjNmjVLy5YtU2VlpRISErRx40YdPXpUubm5jfstXbpUmzZtUllZmfr27StJuueee/TBBx9o3rx5+uijj/TRRx817h8VFaU77rijnZYEAOiqfAo1Sdq0aZMeffRRbd68WSdPnlRSUpJ27NjR5M3VjuPIcZwm8/7nf/5HjuPoxRdf1IsvvthkW3x8PKEGAGgzn0MtPDxc2dnZys7O9rhPbm5ukys3STpy5Ijv1QEA4AP+9AwAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGoQaAMAahBoAwBqEGgDAGj6FWm1trbKyshQXF6cePXooNTVVu3bt8mru119/rczMTPXu3Vs9e/bUpEmTtH///lYVDQCAOz6F2ty5c5WTk6OMjAytWbNGoaGhmjp1qgoKClqcd/bsWd12223Ky8vTwoULlZ2drcrKSk2YMEGHDx9u0wIAALjA5e2O+/bt07Zt27Rq1So98sgjkqSMjAwNGzZMS5YsaTHY8vPztXfvXuXn5+vOO++UJM2ePVuJiYl67LHHtHXr1jYuAwAAH67U8vPz5XK5lJmZ2TgWHh6u+fPna+/evfrss89anHv11Vc3Bpok9erVS7Nnz9Zrr72murq6VpbfMU6cOKH77lusmTNvkzRdM2fepvvuW6wTJ04EujS4UVpaqkGDJmju3GRJ0zV3brIGDZqg0tLSQJeGi9CrzqVTnguNl77//e+boUOHNhvftWuXcRzH7Nixw+PcgQMHmttuu63Z+AsvvGAcxzEffvih23lFRUVGkikqKvK2zDb7/PPPTULCRCPtNdJZIxkjNRhpr0lImGgqKyv9Vgsu7cMPPzQu13Vu+vWecbmuMyUlJYEuEefRq84lmM6FvmSB11csnTTCAAANj0lEQVRqFRUVio2NbTZ+YezYsWMdMtffsrJ+o08+eUJSqiTn/GiIpFR98smvtWRJduCKQzM//OFDqq/PVfN+jVZ9/YuaPv3BwBWHJuhV59JZz4Vev6ZWXV2t8PDwZuMRERGN2z2pqalp9Vx/27evRJKnZt2sd955XMXF/qwILSkv/1LnvuncSVV5+Wn6FSToVefyzjstnwv37Xvcn+V4zetQ6969u2pra5uN19TUNG7viLmStGjRIsXExDQZS0tLU1pa2iXr9lV9faj+8VPJxUJ05EioUlLa/WHRat3UUr/q6q6lX0GDXnUuLZ8Lz50r219eXp7y8vKajJ06dcrr+V6HWmxsrNunCSsqKiRJcXFxHTJXknJycpScnOxtqW3icjVIMnLfzLPq379B+fl+KQVeSE39TnV1nvsVFvap3n/f31XBHXrVucyc2aAjRzz369y5sv25u2ApLi5Wipc/8XgdaiNGjNDbb7+tqqoqRUVFNY4XFhZKkoYPH+5x7vDhw/Xuu+/KGCPH+ccBKiwsVGRkpBITE70to8ONGnW9SkoK5f5pkkKNH3+9/JSv8EK/flfq8OH3JY12s/V99esXTb+CBL3qXMaPv15Hjng+F44adb2/S/KOt3efFBYWGsdxzKpVqxrHampqzMCBA83o0aMbxyoqKkxpaampq6trHNu2bZtxHMfk5+c3jp04ccJcdtllJi0trV3ueGkvlZWV5+/4ee/8nT7/uEOLux+DT0lJyfk76gou6lcBd9QFGXrVuQTTudCXLPA61IwxZvbs2SYsLMwsWbLErFu3zowZM8Z069bNvPvuu4373HvvvcZxHFNeXt441tDQYEaPHm2ioqLMihUrzNq1a83QoUNNTEyMOXToULsspD1VVlaauXN/Zq6/fqpJTJxmrr9+qpk792cEWpAqKSkxAweON926DTNhYcmmW7dhZuDA8ZwkgxC96lyC5VzoSxZ4/fSjJG3atEmPPvqoNm/erJMnTyopKUk7duzQuHHjGvdxHKfJU4ySFBISojfeeEOLFy/WmjVrVF1drVGjRmnTpk0aNGhQ2y8321nv3r2Vm/ubQJcBLw0ZMkQff/x2oMuAF+hV59IZz4WOMcYEughPLrw4WFRU5LcbRQAAwcWXLOBPzwAArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQahdQl5eXqBLCLiufgy6+voljkFXX7/UeY4BoXYJnaWRHamrH4Ouvn6JY9DV1y91nmNAqAEArEGoAQCsQagBAKzhCnQBLamurpYklZaWBqyGU6dOqbi4OGCPHwy6+jHo6uuXOAZdff1SYI/BhQy4kAktcYwxpqMLaq2tW7cqPT090GUAAILAli1bdM8997S4T1CH2hdffKGdO3cqPj5e3bt3D3Q5AIAAqK6uVllZmW699Vb16tWrxX2DOtQAAPAFN4oAAKxBqAEArEGoAQCsQagBAKxBqAEArNElQ622tlZZWVmKi4tTjx49lJqaql27dnk19+uvv1ZmZqZ69+6tnj17atKkSdq/f38HV9z+WnsM/vKXv2jevHlKTExUZGSkEhIStGDBAh0/ftwPVbeftnwN/LMFCxYoJCRE06ZN64AqO1Zbj8GuXbs0adIkXXbZZYqOjtZNN92kl19+uQMrbl9tWX9RUZFuv/12xcbGKioqSklJSXrmmWd09uzZDq66/XzzzTd67LHHNGXKFF1xxRUKCQnRxo0bvZ4ftOdC0wXdfffdJiwszCxZssQ8//zzZsyYMSYsLMzs2bOnxXkNDQ1mzJgxpmfPnmbFihVm7dq1ZujQoSY6Otp8/PHHfqq+fbT2GKSkpJiEhASzdOlS8/vf/94sX77cREdHm6uvvtocP37cT9W3XWvX/88++OADExYWZrp3726mTZvWgdV2jLYcgxdffNGEhISYKVOmmGeffdasW7fOLFq0yPzHf/yHHypvH61d/1//+lfTrVs3c8MNN5jVq1eb9evXmx/+8IfGcRzz8MMP+6n6tjty5IhxHMfEx8ebiRMnGsdxzMaNG72aG8znwi4XaoWFhcZxnCbffDU1NWbgwIFmzJgxLc7dtm2bcRzHvPLKK41jJ06cMJdffrn50Y9+1GE1t7e2HIN333232dju3buN4zjmF7/4RbvX2hHasv4Lzp49a0aPHm3uv/9+Ex8f3+lCrS3H4MiRI6Z79+7mJz/5SUeX2WHasv4FCxaYiIgIc/LkySbj48ePNzExMR1Sb0eora01n3/+uTHmXFD7EmrBfC7scqG2ePFiExYWZqqqqpqMP/nkk8ZxHPPpp596nDtr1iwTGxvbbPyBBx4wkZGR5rvvvmv3ejtCW46BJ1deeaWZOXNme5XYodpj/Rs3bjQxMTHm+PHjpl+/fp0u1NpyDLKyskxERIQ5ffq0McaYqqoqc/bs2Q6tt721Zf1z5swxMTExzdY8Z84ct+eHzuCDDz7wKdSC+VzY5V5T279/vxITE9WzZ88m4yNHjpQkHThwoMW5ycnJzcZHjhypb7/9VocOHWrfYjtIW46BO2fOnFFVVdUlf31NsGjr+quqqpSVlaXly5erT58+HVZnR2rLMdi1a5cGDx6sHTt26Nprr1V0dLR69eqlX/7ylzKd5BcUtWX9EydO1OnTp/XAAw/o4MGDKi8v13PPPadXX31Vy5Yt69C6g0UwnwuD+rf0d4SKigrFxsY2G78wduzYsRbnTpgwocW5Q4cObZ9CO1BbjoE7q1evVl1dnebMmdMu9XW0tq5/xYoVioyM1KJFizqkPn9oyzH4+OOP5XK5NG/ePGVlZSkpKUmvvPKKHn/8cdXX1+uJJ57osLrbS1vWv2DBAn300Udat26dXnjhBUlSaGio1q5dq8zMzI4pOMgE87mwy4VadXW1wsPDm41HREQ0bvekpqam1XODSVuOwcV2796tX/3qV5ozZ47bL/Jg1Jb1Hzp0SGvWrNEf/vAHhYWFdViNHa0tx+DMmTMyxujpp5/W4sWLJUkzZszQV199pd/+9rdavnx5syugYNOW9YeEhGjAgAGaMmWKZs2apYiICL300kv68Y9/rD59+uiOO+7osLqDRTCfC7tcqHXv3l21tbXNxmtqahq3d8TcYNJe6zh48KBmzJihG2+8sfEn1s6gLet/+OGHNXbsWM2YMaPD6vOHtn4fVFdXKy0trcn43XffrT/96U86cOCAxo0b174Ft7O2rP+pp57SmjVrdPjwYfXo0UOSNHPmTE2aNEkPPfSQbr/9doWGhnZM4UEimM+FXe41tdjYWLdPLVRUVEiS4uLiOmRuMGmPdfz973/X5MmTdfnll+uNN95QZGRku9fZUVq7/rfeeks7d+7UwoULVVZW1vivvr5e3377rcrLy1VVVdWhtbeXtnwNXNh28euJV111lSTp5MmT7VVmh2nL+p999lndcsstjYF2wbRp03Ts2DGVl5e3b7FBKJjPhV0u1EaMGKFDhw41O/kUFhZKkoYPH+5x7vDhw1VcXNzsxfDCwkJFRkYqMTGx/QvuAG05BpL05ZdfavLkyaqrq9POnTs73c0SrV3/0aNHJUl33nmnBgwY0Pjv2LFjeuutt9S/f3/l5uZ2bPHtpC1fAzfddJOMMfr000+bjF84yfXu3budq21/bVl/ZWWlGhoamo3X1dVJkurr69ux0uAU1OfCgN13GSAX3p+yatWqxrEL708ZPXp041hFRYUpLS01dXV1jWMX3puRn5/fOHbixAlz2WWXmbS0NP8soB205RicOXPGjBo1ysTExJji4mK/1t1eWrv+o0ePmtdee63Jv+3bt5urrrrKjBo1yrz22mvmk08+8ft6WqMtXwPbt283juOYn//8541jDQ0NZty4caZXr16d4q0tbVn/DTfcYK688krz5ZdfNo7V19eblJQUExMTY+rr6/2ziHbU0i39ne1c2OVCzRhjZs+e3fibBNatW2fGjBljunXr1uSNxffee69xHMeUl5c3jjU0NJjRo0ebqKioJu+ij4mJMYcOHQrEUlqttcfgjjvuMI7jmPnz55vNmzc3+bd9+/ZALKVVWrt+dzrj+9SMadsx+P73v29CQkLMAw88YNauXWv+7d/+zTiOY55//nl/L6PVWrv+rVu3GsdxzMCBA012drZZs2aNGT16tHEcxzzxxBOBWEqrPfPMM2blypXmwQcfNI7jmLvuususXLnSrFy50pw6dcoY0/nOhV0y1GpqaszixYtNbGysiYiIMDfffLP585//3GSfuXPnmpCQkGbfzCdPnjT333+/6dWrl4mMjDQTJ040RUVF/iy/XbT2GMTHx5uQkBDjOE6zf/379/f3MlqtLV8DF+uMv1HEmLYdgzNnzpif/OQnJjY21oSHh5ukpCTz0ksv+bP8NmvL+nfu3GkmTJhgevfu3bj+9evX+7P8dhEfH9/4/RsSEtL4vf3Pa+5s50LHmE7ybkkAAC6hy90oAgCwF6EGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsMb/B+ruUpOSe0bgAAAAAElFTkSuQmCC", + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbUAAAG7CAYAAABEoAutAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X9UVXW+//HX4YcHRKASM6hRE2NyNEFJwh/3+qN7zTG1LMwY5Q5q0uo2y8bmKtqM0y2bfjDei0NTN50mAjPGLk3WctU4y2k1FlexAb2trz+u6RKYAgcrf2ABcfDz/cPxzCAcOgc4HPyc52MtVzPvvT+c92fzYb/Y+2zAYYwxAgDAAiGBbgAAgJ5CqAEArEGoAQCsQagBAKxBqAEArEGoAQCsQagBAKwRFugGOvPZZ59px44dGjZsmCIjIwPdDgAgABobG1VVVaXbbrtNcXFxne7bp0Ntx44dWrRoUaDbAAD0Aa+88ooWLlzY6T59OtSGDRsm6cJERo4cGZAeVqxYofz8/IC8dl8R7Mcg2OcvcQyCff5SYI/BoUOHtGjRIncmdKZPh9rFW44jR47UuHHjAtJDbGxswF67rwj2YxDs85c4BsE+f6lvHANv3obiQREAgDUINQCANQg1AIA1CLVvkJmZGegWAi7Yj0Gwz1/iGAT7/KXL5xg4+vLfU6usrFRqaqoqKioC/gYlACAwfMkCrtQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADW8CnUvvzySz366KOaOXOmrrrqKoWEhKioqMjr8adPn1ZOTo4GDRqkAQMGaPr06dq3b5/PTQMA0BGfQu3kyZNat26d/u///k8pKSmSJIfD4dXY8+fP6/bbb1dJSYmWL1+uvLw81dfXa+rUqTp69KjvnQMAcIkwX3ZOSEjQiRMndPXVV6uiokLjx4/3emxpaal2796t0tJS3XXXXZKke+65R0lJSXr00Ue1ZcsW3zoHAOASPl2p9evXT1dffbUkyRjj0wuVlpbqmmuucQeaJMXFxemee+7Rm2++qZaWFp8+HgAAl+q1B0X27duncePGtauPHz9eX331lY4cOdJbrXyjQ4cOaeoNN+gmp1Op/frpJqdTU2+4QYcOHQp0awgiZWVlGh0ToxsdDiU7HLrR4dDomBiVlZUFujUEiZMnT2rl4sW6fdQozf32t3X7qFFauXixTp48GejWPPLp9mN31NXVaerUqe3q8fHxkqTa2lqNGjWqt9rx6MCBA7o7JUUvu1y6RZJD0nlJ5UePat6YMXrjo480cuTIAHcJ2+3atUtLpkzRK1LbddjQoKzJk1X0wQeaNGlSYJuE1err63XvxIl68tgx5elva3DvwYNa8P772rp7twYNGhTgLtvrtSu1pqYmOZ3OdvWIiAhJUmNjY2+10qkH77xThS6X0nXhkyhdOEgTJL3kcumBuXMD1xyCxr/Onq3NUofrsFjS/d/9bqBaQ5D4eW6unjx2rN0aTJf0s2PHlLdqVeCa60SvXalFRkaqubm5Xb2pqcm9vS/4vKZG6R62pUv6vLpaqqzszZYQhEIbGjpdh66GBtYh/OrgH/+oPA/bbpH0xN69vdmO13ot1OLj41VbW9uuXldXJ+nCk5WerFixQrGxsW1qmZmZyszM7NkmJfUzRp5+SCFEUr+WFik1tcdfF/h7Q6VO16FTYh3Cr0LV+RoMdbn88rolJSUqKSlpUztz5ozX43st1FJSUvT+++/LGNPmZ9vKy8sVFRWlpKQkj2Pz8/M7fMjEH752OGTU8SfzvKSvw8OlPXt6pRcEr+rU1E7XYbMkVVT0ak8ILq0ZGTLHj3tcg61h/omPji5YKisrlerlN3F+6erEiRM6ffq0RowYobC/TjwjI0OlpaX67W9/q7vvvluS9Nlnn+m///u/NWfOHIWHh/ujFZ8NHDJEe44e1YQOtu2RNHDoUKmXAhbBqzU6WnsaGjyuw7DoaNYh/Oo7U6ao/PjxDm+Dl0v6Tlpab7fkFZ9D7Ze//KVOnz7tvpX41ltvqaamRpK0fPlyxcTEaPXq1SouLlZVVZWGDBki6UKopaena/HixTp48KAGDhyo559/XsYYPfbYYz04pe75r7fe0rwxY/TSXx8WCdGF70r2SFoSFqY33norsA0iKGx85x1lTZ6sYqndOvwXSUXvvBPI9hAEVuXlacH77+tnx47pFv1tDZZL+nFiorbmeXrHLcCMj4YNG2YcDodxOBwmJCTEhISEuP93dXW1McaY7OzsNv//olOnTpn77rvPxMXFmaioKDNt2jRTUVHh8bUqKiqMpE738YeDBw+aKSNGmNHh4WacZEaHh5spI0aYgwcP9mofCG4ffPCBGRUdbcZIZo5kxkhmVHS0+eCDDwLdGoJEfX29+bfsbDPr+uvNHMnMuv5682/Z2aa+vr5X+/AlCxzG+PirQXrRxfuoFRUVvfae2iUNXHgzvqKCWz0InC1bpEWLpFdekRYuDHQ3CEYBPhf6kgX86RkAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANXwKtebmZuXm5iohIUH9+/dXenq6du7c6dXYiooKzZ49W/Hx8YqOjlZycrKeffZZnT9/vkuNAwBwKZ9CLTs7W/n5+crKylJBQYFCQ0M1a9YslZWVdTquoqJCEydOVE1NjVavXq3//M//1PDhw/XQQw/p4Ycf7tYEAAC4KMzbHffu3autW7dq/fr17iDKysrS6NGjtWrVqk6DbePGjQoJCdGuXbt0xRVXSJKWLVumqVOn6uWXX9aGDRu6OQ0AAHy4UistLVVYWJhycnLcNafTqaVLl2r37t369NNPPY49e/asnE6nYmNj29SvueYa9e/fvwttAwDQntehtm/fPiUlJWnAgAFt6uPHj5ck7d+/3+PYadOm6ezZs7r//vt1+PBhVVdX64UXXtAbb7yhNWvWdLF1AADa8vr2Y11dneLj49vVL9Zqa2s9jl22bJkOHDigjRs36sUXX5QkhYaG6rnnnmtz5QcAQHd4HWqNjY1yOp3t6hEREe7tnoSEhGj48OGaOXOm5s+fr4iICL366qv6wQ9+oMGDB+uOO+7oQusAALTldahFRkaqubm5Xb2pqcm93ZOnn35aBQUFOnr0qPs9tIyMDE2fPl0PPvigZs+erdDQUI/jV6xY0e79uMzMTGVmZnrbPgDgMlBSUqKSkpI2tTNnzng93utQi4+P7/AWY11dnSQpISHB49jnn39et956a7uHQubMmaMf/ehHqq6u1vDhwz2Oz8/P17hx47xtFQBwmerogqWyslKpqalejff6QZGxY8fqyJEjamhoaFMvLy+XJKWkpHgcW19fr9bW1nb1lpYWSZLL5fK2DQAAPPI61DIyMtTa2qpNmza5a83NzSosLFR6erquvfZaSdKJEyd0+PDhNkGVlJSk3//+9/riiy/ctdbWVr322muKiYlRYmJiT8wFABDkvL79mJaWpvnz52vNmjWqr69XYmKiioqKVFNTo8LCQvd+q1evVnFxsaqqqjRkyBB3bdGiRbrllluUk5OjiIgIlZSUqLKyUj/72c86fT8NAABveR1qklRcXKy1a9dq8+bNOnXqlJKTk7V9+3ZNnjzZvY/D4ZDD4Wgz7nvf+57i4uL01FNP6ec//7nOnj2rG2+8URs3btSyZct6ZiYAgKDnU6g5nU7l5eUpLy/P4z6FhYVtrtwumjFjhmbMmOF7hwAAeIk/PQMAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwhk+h1tzcrNzcXCUkJKh///5KT0/Xzp07vR6/c+dOTZ8+XVdccYViYmJ0880367XXXvO5aQAAOuJTqGVnZys/P19ZWVkqKChQaGioZs2apbKysm8cW1hYqNtuu01Op1NPPfWU1q9fr3/8x3/UJ5980uXmAQD4e2He7rh3715t3bpV69ev18MPPyxJysrK0ujRo7Vq1apOg62qqkoPPvigli9frvz8/O53DQBAB7y+UistLVVYWJhycnLcNafTqaVLl2r37t369NNPPY594YUXZIzR448/Lkk6d+6cjDHdaBsAgPa8DrV9+/YpKSlJAwYMaFMfP368JGn//v0ex+7cuVM33nijtm/fruuuu04xMTGKi4vTT3/6U8INANBjvL79WFdXp/j4+Hb1i7Xa2lqPYz/++GOFhYVpyZIlys3NVXJysl5//XU98cQTcrlcevLJJ7vQOgAAbXkdao2NjXI6ne3qERER7u2eXLzd+Mwzz2jlypWSpHnz5umLL77QL37xCz3yyCPtrgABAPCV16EWGRmp5ubmdvWmpib39s7GNjY2KjMzs0393nvv1e9+9zvt379fkydP9jh+xYoVio2NbVPLzMxs9/EAAJe3kpISlZSUtKmdOXPG6/Feh1p8fHyHtxjr6uokSQkJCR7HJiQk6NixYxo8eHCb+tVXXy1JOnXqVKevnZ+fr3HjxnnbKgDgMtXRBUtlZaVSU1O9Gu/1gyJjx47VkSNH1NDQ0KZeXl4uSUpJSfE49uabb5Yxpt3PpF0MyUGDBnnbBgAAHnkdahkZGWptbdWmTZvctebmZhUWFio9PV3XXnutJOnEiRM6fPiwXC6Xe78FCxZIkn7961+7a+fPn1dhYaEGDhzodQIDANAZr28/pqWlaf78+VqzZo3q6+uVmJiooqIi1dTUqLCw0L3f6tWrVVxcrKqqKg0ZMkSSdMcdd+jWW2/VU089pc8++0xjxozRtm3bVFZWpk2bNik8PLznZwYACDpeh5okFRcXa+3atdq8ebNOnTql5ORkbd++vc1DHg6HQw6Ho93Ybdu26Sc/+Ym2bt2ql19+WTfeeKO2bNnCwx4AgB7jU6g5nU7l5eUpLy/P4z6FhYVtrtwuioqKUn5+Pr8mCwDgN/zpGQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQi1DpSVlWl0TIySU1M1V1JyaqpGx8SorKws0K0hiGzZskXDQ0KUvGjRhXW4aJGGh4Roy5YtgW4NQeLQoUOaesMNGpeerrmSxqWna+oNN+jQoUOBbs2jsEA30Nfs2rVLS6ZM0SuSbpHkkHReUnlDg7ImT1bRBx9o0qRJgW0S1isqKtLj2dl6VZesQ2O0aNEiSdLChQsD2CFsd+DAAd2dkqKXXa6/rcGWFpUfPap5Y8bojY8+0siRIwPdZjtcqV3iX2fP1mZJ6brwSZQuHKQJkool3f/d7waqNQSRxxYv1ivqeB1ulrQ2KytQrSFIPHjnnSp0uTpcgy+5XHpg7tzANdcJrtQu4WpoULqHbemSQhsaJG7/wM8GGNPpOuxnDOsQfnW2qqrTNfh5TU1vtuM1Qu0STv3tu5JLhUgaKkl/vf0D+Mswdb4OoyTWIfzqOnW+BvsZ04vdeI9Qu0SzJKOOP5nnJVVL0iuv9GZLCEJVixZ1ug6/lFiH8KtPsrNlXC6Pa/Brh6fICyxC7RJh0dHa09CgCR1s2yOpNTpa4g16+Nm5rCztMcbjOvza4WAdwq9i/v3ftefoUY9rcOCQIb3dkld4UOQSG995R1mS/kcXvhvRX//7P5L+5a/bAX9bt3mzFqnjdZj11+2AP/3XW29pcVhYh2twSViY/uuttwLXXCe4UrvEpEmTVPTBB8r57nflamiQUxduSYZFR6vonXd4nB+94uLj+ouyshRujCIkNUlqcTi0bvNmHueH340cOVJvfPSRHpg7V59XV6tfS4u+Dg/XwKFD9cZbb/XJx/klQq1DkyZN0v87e1aqrJRSU6WKCmncuEC3hSCzcOHCC+G1ZcuFh0JeeYVbjuhVI0eO1Hsff/y3c+GePX3+XMjtRwCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1fAq15uZm5ebmKiEhQf3791d6erp27tzp84suW7ZMISEhmjNnjs9jAQDwxKdQy87OVn5+vrKyslRQUKDQ0FDNmjVLZWVlXn+MP/3pTyoqKlJERIQcDofPDQMA4InXobZ3715t3bpVTz/9tJ555hndd999evfddzV06FCtWrXKq49hjNHy5cv1/e9/X4MHD+5y0wAAdMTrUCstLVVYWJhycnLcNafTqaVLl2r37t369NNPv/FjbN68WQcPHtQTTzwhY0zXOgYAwAOvQ23fvn1KSkrSgAED2tTHjx8vSdq/f3+n4xsaGpSbm6tHHnmEqzQAgF94HWp1dXWKj49vV79Yq62t7XT8448/rqioKK1YscLHFgEA8E6Ytzs2NjbK6XS2q0dERLi3e3LkyBEVFBToN7/5jcLDw7vQJgAA38zrK7XIyEg1Nze3qzc1Nbm3e/LQQw9p0qRJmjdvXhdaBADAO15fqcXHx3d4i7Gurk6SlJCQ0OG4d999Vzt27NBvf/tbVVVVuesul0tfffWVqqurddVVVyk6Otrja69YsUKxsbFtapmZmcrMzPS2fQDAZaCkpEQlJSVtamfOnPF6vNehNnbsWL333ntqaGhoE0Dl5eWSpJSUlA7H1dTUSJLuuuuudttqa2t1/fXXa8OGDVq+fLnH187Pz9e4ceO8bRUAcJnq6IKlsrJSqampXo33OtQyMjK0fv16bdq0ST/60Y8kXfgNI4WFhUpPT9e1114rSTpx4oROnz6tESNGKCwsTLfeequ2bdvW5mMZY5STk6Nhw4bpxz/+sUaPHu1tGwAAeOR1qKWlpWn+/Plas2aN6uvrlZiYqKKiItXU1KiwsNC93+rVq1VcXKyqqioNGTJE3/rWt/Stb32r3cd76KGHNHjwYM2dO7dnZgIACHpeh5okFRcXa+3atdq8ebNOnTql5ORkbd++XZMnT3bv43A4vPr1V/yKLABAT/Mp1JxOp/Ly8pSXl+dxn8LCwjZXbp4cP37cl5cGAOAb8adnAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANYg1AAA1iDUAADWINQAANbwOdSam5uVm5urhIQE9e/fX+np6dq5c+c3jvvDH/6gJUuWKCkpSVFRUUpMTNSyZct04sSJLjUOAMClfA617Oxs5efnKysrSwUFBQoNDdWsWbNUVlbW6bjc3Fzt2rVLd999t5599lnde++9eu211zR27Fj95S9/6fIEAAC4KMyXnffu3autW7dq/fr1evjhhyVJWVlZGj16tFatWtVpsG3YsEGTJ09uU5s5c6amTJmiX/7yl1q3bl0X2gcA4G98ulIrLS1VWFiYcnJy3DWn06mlS5dq9+7d+vTTTz2OvTTQJOkf/uEfdNVVV+nw4cO+tAEAQId8CrV9+/YpKSlJAwYMaFMfP368JGn//v0+vfi5c+fU0NCguLg4n8YBANARn0Ktrq5O8fHx7eoXa7W1tT69+IYNG9TS0qIFCxb4NA4AgI74FGqNjY1yOp3t6hEREe7t3tq1a5cee+wxLViwQFOnTvWlDQAAOuTTgyKRkZFqbm5uV29qanJv98bhw4c1b948jRkzRi+++OI37r9ixQrFxsa2qWVmZiozM9Or1wMAXB5KSkpUUlLSpnbmzBmvx/sUavHx8R3eYqyrq5MkJSQkfOPH+POf/6wZM2boyiuv1Ntvv62oqKhvHJOfn69x48b50ioA4DLU0QVLZWWlUlNTvRrv0+3HsWPH6siRI2poaGhTLy8vlySlpKR0Ov7zzz/XjBkz1NLSoh07dmjw4MG+vDwAAJ3yKdQyMjLU2tqqTZs2uWvNzc0qLCxUenq6rr32WknSiRMndPjwYblcLvd+X375pWbNmqW6ujq9/fbbSkxM7KEpAABwgU+3H9PS0jR//nytWbNG9fX1SkxMVFFRkWpqalRYWOjeb/Xq1SouLlZVVZWGDBkiSVq4cKE+/PBDLVmyRAcOHNCBAwfc+0dHR+uOO+7ooSkBAIKVT6EmScXFxVq7dq02b96sU6dOKTk5Wdu3b2/zw9UOh0MOh6PNuP/93/+Vw+HQSy+9pJdeeqnNtmHDhhFqAIBu8znUnE6n8vLylJeX53GfwsLCNlduknT8+HHfuwMAwAf86RkAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANQg1AIA1CDUAgDUINQCANXwKtebmZuXm5iohIUH9+/dXenq6du7c6dXY06dPKycnR4MGDdKAAQM0ffp07du3r0tNAwDQEZ9CLTs7W/n5+crKylJBQYFCQ0M1a9YslZWVdTru/Pnzuv3221VSUqLly5crLy9P9fX1mjp1qo4ePdqtCQAAcFGYtzvu3btXW7du1fr16/Xwww9LkrKysjR69GitWrWq02ArLS3V7t27VVpaqrvuukuSdM899ygpKUmPPvqotmzZ0s1pAADgw5VaaWmpwsLClJOT4645nU4tXbpUu3fv1qefftrp2GuuucYdaJIUFxene+65R2+++aZaWlq62L5/nDx5UisXL9btGRmaK+n2jAytXLxYJ0+eDHRrCCKHDh3S1Btu0LjsbM2VNC47W1NvuEGHDh0KdGsIEpfjudDrK7V9+/YpKSlJAwYMaFMfP368JGn//v269tprPY4dN25cu/r48eO1adMmHTlyRKNGjfKlb7+pr6/XvRMn6sljx5QnySHp/PHj2nv8uBa8/7627t6tQYMGBbpNWO7AgQO6OyVFL7tcukV/XYcul8qPHtW8MWP0xkcfaeTIkYFuExa7XM+FXl+p1dXVKT4+vl39Yq22ttYvY3vbz3Nz9eSxY0rXhU+idOEgpUv62bFjylu1KnDNIWg8eOedKnS52q3DCZJecrn0wNy5gWsOQeFyPRd6faXW2Ngop9PZrh4REeHe7klTU1OXx/a2g3v3Ks/DtlskPfHHP0qVlb3ZEoLQ2epqpXvYli7p8+pq1iH86uAf/9j5uXDv3t5sx2teh1pkZKSam5vb1Zuamtzb/TFWklasWKHY2Ng2tczMTGVmZn5j374Kdbnc35VcKkRS6PHjUmpqj78u8Peukzpdh/1aWliH8KtQdb4GQ10uv7xuSUmJSkpK2tTOnDnj9XivQy0+Pr7D24R1dXWSpISEBL+MlaT8/PwO35Pzh9awMBl1/Mk8L6n1+uul0tJe6QXB65P0dJmWFo/r8OvwcGnPnt5uC0GkNSND5vhxz+fCMK/jwycdXbBUVlYq1ctv4rzuauzYsXrvvffU0NCg6Ohod728vFySlJKS4nFsSkqK3n//fRlj5HD87RCVl5crKipKSUlJ3rbhd99JS1P5wYMd3vopl/SdKVOkXgpYBK+YoUO15+hRTehg2x5JA4cOZR3Cr74zZYrKjx/3fC5MS+vtlrxjvFReXm4cDodZv369u9bU1GRGjBhhJkyY4K7V1dWZQ4cOmZaWFndt69atxuFwmNLSUnft5MmT5oorrjCZmZkeX7OiosJIMhUVFd622W319fVmWmKi+R/JtErG/PW//yOZaYmJpr6+vtd6QfA6ePCg+XZYmCm7ZB2WSebbYWHm4MGDgW4RlutL50JfssDrK7W0tDTNnz9fa9asUX19vRITE1VUVKSamhoVFha691u9erWKi4tVVVWlIUOGSJIyMjKUnp6uxYsX6+DBgxo4cKCef/55GWP02GOP9XROd8ugQYO0dfdu5a1apSf27lWoy6XWsDB9Jy1NW/Py+uQjrLDPyJGq2G7gAAAK00lEQVQj9cZHH+mBuXP1eU2N+hmjrx0ODRwyRG+89RaP88PvLtdzoU83RYuLi7V27Vpt3rxZp06dUnJysrZv367Jkye793E4HG1uMUpSSEiI3n77ba1cuVIFBQVqbGxUWlqaiouLdcMNN/TMTHrQoEGD9PO/C2ogEEaOHKn3Pv440G0giF2O50KHMcYEuglPLr45WFFR0WsPigAA+hZfsoA/PQMAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqEGALAGoQYAsAahBgCwBqH2DUpKSgLdQsAF+zEI9vlLHINgn790+RwDQu0bXC6fSH8K9mMQ7POXOAbBPn/p8jkGhBoAwBqEGgDAGoQaAMAaYYFuoDONjY2SpEOHDgWshzNnzqiysjJgr98XBPsxCPb5SxyDYJ+/FNhjcDEDLmZCZxzGGOPvhrpqy5YtWrRoUaDbAAD0Aa+88ooWLlzY6T59OtQ+++wz7dixQ8OGDVNkZGSg2wEABEBjY6Oqqqp02223KS4urtN9+3SoAQDgCx4UAQBYg1ADAFiDUAMAWINQAwBYg1ADAFgjKEOtublZubm5SkhIUP/+/ZWenq6dO3d6Nfb06dPKycnRoEGDNGDAAE2fPl379u3zc8c9r6vH4A9/+IOWLFmipKQkRUVFKTExUcuWLdOJEyd6oeue05018PeWLVumkJAQzZkzxw9d+ld3j8HOnTs1ffp0XXHFFYqJidHNN9+s1157zY8d96zuzL+iokKzZ89WfHy8oqOjlZycrGeffVbnz5/3c9c958svv9Sjjz6qmTNn6qqrrlJISIiKioq8Ht9nz4UmCN17770mPDzcrFq1yvzqV78yEydONOHh4eaDDz7odFxra6uZOHGiGTBggHn88cfNc889Z0aNGmViYmLMxx9/3Evd94yuHoPU1FSTmJhoVq9ebX7961+bRx55xMTExJhrrrnGnDhxope6776uzv/vffjhhyY8PNxERkaaOXPm+LFb/+jOMXjppZdMSEiImTlzpnn++efNxo0bzYoVK8x//Md/9ELnPaOr8//Tn/5k+vXrZ2666SazYcMGs2nTJnPnnXcah8NhHnrooV7qvvuOHz9uHA6HGTZsmJk2bZpxOBymqKjIq7F9+VwYdKFWXl5uHA5Hmy++pqYmM2LECDNx4sROx27dutU4HA7z+uuvu2snT540V155pfne977nt557WneOwfvvv9+utmvXLuNwOMxPfvKTHu/VH7oz/4vOnz9vJkyYYO677z4zbNiwyy7UunMMjh8/biIjI80Pf/hDf7fpN92Z/7Jly0xERIQ5depUm/qUKVNMbGysX/r1h+bmZvOXv/zFGHMhqH0Jtb58Lgy6UFu5cqUJDw83DQ0NbepPPfWUcTgc5pNPPvE4dv78+SY+Pr5d/f777zdRUVHm66+/7vF+/aE7x8CTgQMHmoyMjJ5q0a96Yv5FRUUmNjbWnDhxwgwdOvSyC7XuHIPc3FwTERFhzp49a4wxpqGhwZw/f96v/fa07sx/wYIFJjY2tt2cFyxY0OH54XLw4Ycf+hRqfflcGHTvqe3bt09JSUkaMGBAm/r48eMlSfv37+907Lhx49rVx48fr6+++kpHjhzp2Wb9pDvHoCPnzp1TQ0PDN/76mr6iu/NvaGhQbm6uHnnkEQ0ePNhvffpTd47Bzp07deONN2r79u267rrrFBMTo7i4OP30pz+VuUx+QVF35j9t2jSdPXtW999/vw4fPqzq6mq98MILeuONN7RmzRq/9t1X9OVzYZ/+Lf3+UFdXp/j4+Hb1i7Xa2tpOx06dOrXTsaNGjeqZRv2oO8egIxs2bFBLS4sWLFjQI/35W3fn//jjjysqKkorVqzwS3+9oTvH4OOPP1ZYWJiWLFmi3NxcJScn6/XXX9cTTzwhl8ulJ5980m9995TuzH/ZsmU6cOCANm7cqBdffFGSFBoaqueee045OTn+abiP6cvnwqALtcbGRjmdznb1iIgI93ZPmpqaujy2L+nOMbjUrl279Nhjj2nBggUdLvK+qDvzP3LkiAoKCvSb3/xG4eHhfuvR37pzDM6dOydjjJ555hmtXLlSkjRv3jx98cUX+sUvfqFHHnmk3RVQX9Od+YeEhGj48OGaOXOm5s+fr4iICL366qv6wQ9+oMGDB+uOO+7wW999RV8+FwZdqEVGRqq5ubldvampyb3dH2P7kp6ax+HDhzVv3jyNGTPG/R3r5aA783/ooYc0adIkzZs3z2/99Ybufh00NjYqMzOzTf3ee+/V7373O+3fv1+TJ0/u2YZ7WHfm//TTT6ugoEBHjx5V//79JUkZGRmaPn26HnzwQc2ePVuhoaH+abyP6MvnwqB7Ty0+Pr7DWwt1dXWSpISEBL+M7Ut6Yh5//vOfNWPGDF155ZV6++23FRUV1eN9+ktX5//uu+9qx44dWr58uaqqqtz/XC6XvvrqK1VXV6uhocGvvfeU7qyBi9sufT/x6quvliSdOnWqp9r0m+7M//nnn9ett97qDrSL5syZo9raWlVXV/dss31QXz4XBl2ojR07VkeOHGl38ikvL5ckpaSkeBybkpKiysrKdm+Gl5eXKyoqSklJST3fsB905xhI0ueff64ZM2aopaVFO3bsuOwelujq/GtqaiRJd911l4YPH+7+V1tbq3fffVfXX3+9CgsL/dt8D+nOGrj55ptljNEnn3zSpn7xJDdo0KAe7rbndWf+9fX1am1tbVdvaWmRJLlcrh7stG/q0+fCgD13GSAXfz5l/fr17trFn0+ZMGGCu1ZXV2cOHTpkWlpa3LWLP5tRWlrqrp08edJcccUVJjMzs3cm0AO6cwzOnTtn0tLSTGxsrKmsrOzVvntKV+dfU1Nj3nzzzTb/tm3bZq6++mqTlpZm3nzzTXPs2LFen09XdGcNbNu2zTgcDvPjH//YXWttbTWTJ082cXFxl8WPtnRn/jfddJMZOHCg+fzzz901l8tlUlNTTWxsrHG5XL0ziR7U2SP9l9u5MOhCzRhj7rnnHvdvEti4caOZOHGi6devX5sfLP7+979vHA6Hqa6udtdaW1vNhAkTTHR0dJufoo+NjTVHjhwJxFS6rKvH4I477jAOh8MsXbrUbN68uc2/bdu2BWIqXdLV+Xfkcvw5NWO6dwz+6Z/+yYSEhJj777/fPPfcc+af//mfjcPhML/61a96expd1tX5b9myxTgcDjNixAiTl5dnCgoKzIQJE4zD4TBPPvlkIKbSZc8++6xZt26deeCBB4zD4TB33323WbdunVm3bp05c+aMMebyOxcGZag1NTWZlStXmvj4eBMREWFuueUW8/vf/77NPtnZ2SYkJKTdF/OpU6fMfffdZ+Li4kxUVJSZNm2aqaio6M32e0RXj8GwYcNMSEiIcTgc7f5df/31vT2NLuvOGrjU5fgbRYzp3jE4d+6c+eEPf2ji4+ON0+k0ycnJ5tVXX+3N9rutO/PfsWOHmTp1qhk0aJB7/ps2berN9nvEsGHD3F+/ISEh7q/tv5/z5XYudBhzmfy0JAAA3yDoHhQBANiLUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBYg1ADAFiDUAMAWINQAwBY4/8DyWCjnQNsR7QAAAAASUVORK5CYII=", + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.1,1.1)" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n", "\n", - "mortar_side = mesh_upper\n", - "nonmortar_side = mesh_lower\n", + "using PyPlot\n", "\n", + "mortar_surface = :LOWER_TO_UPPER\n", + "slave_surface = :UPPER_TO_LOWER\n", + "\n", + "fig1 = figure(figsize=(5, 5))\n", "master_elements = JuliaFEM.Core.Element[]\n", - "for (elid, (eltype, elset, elcon)) in mortar_side[\"connectivity\"]\n", - " (eltype == :QU4) && (elset == :UPPER_BLOCK_TO_LOWER_BLOCK) || continue\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == mortar_surface) || continue\n", " element = Quad4(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[mortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", + " X = element(\"geometry\", 0.0)\n", + " x = [X[mod(i, 4)+1][1] for i=1:5]\n", + " y = [X[mod(i, 4)+1][2] for i=1:5]\n", + " plot(x, y, \"-bo\")\n", " push!(master_elements, element)\n", "end\n", "\n", + "xlim(-0.1, 1.1)\n", + "ylim(-0.1, 1.1)\n", + "\n", + "fig1 = figure(figsize=(5, 5))\n", "contact_problem = MortarProblem(\"displacement\", 3)\n", - "for (elid, (eltype, elset, elcon)) in nonmortar_side[\"connectivity\"]\n", - " (eltype == :QU4) && (elset == :LOWER_BLOCK_TO_UPPER_BLOCK) || continue\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == slave_surface) || continue\n", " element = Quad4(elcon)\n", - " element[\"geometry\"] = Vector{Float64}[nonmortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n", + " update(element, \"geometry\", mesh[\"nodes\"])\n", " element[\"master elements\"] = master_elements\n", + " X = element(\"geometry\", 0.0)\n", + " x = [X[mod(i, 4)+1][1] for i=1:5]\n", + " y = [X[mod(i, 4)+1][2] for i=1:5]\n", + " plot(x, y, \"-ro\")\n", " calculate_normal_tangential_coordinates!(element, 0.0)\n", " push!(contact_problem, element)\n", - "end" + "end\n", + "xlim(-0.1, 1.1)\n", + "ylim(-0.1, 1.1)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -735,13 +577,13 @@ "INFO: Starting iteration 1\n", "INFO: Assembling field problems...\n", "INFO: Assembling body 1...\n", - "INFO: dim = 369\n", + "INFO: dim = 273\n", "INFO: Assembling boundary problems...\n", "INFO: Assembling boundary 1...\n", "INFO: Assembling boundary 2...\n", "INFO: dumping matrices to disk, file = matrices_tie_contact_3d_host_1_iteration_1.jld\n", "INFO: Solving system\n", - "INFO: UMFPACK: solved in 0.18637418746948242 seconds. norm = 0.7276121358181449\n", + "INFO: UMFPACK: solved in 0.21015095710754395 seconds. norm = 0.5778929077881374\n", "INFO: timing info for iteration:\n" ] }, @@ -751,7 +593,7 @@ "(1,true)" ] }, - "execution_count": 10, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, @@ -759,13 +601,13 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO: boundary assembly : 1.6767330169677734\n", - "INFO: field assembly : 1.1170330047607422\n", - "INFO: dump matrices to disk : 1.0959408283233643\n", - "INFO: solve problem : 0.3340458869934082\n", - "INFO: update element data : 0.022793054580688477\n", - "INFO: non-linear iteration : 4.246567964553833\n", - "INFO: solver finished in 4.391422986984253 seconds.\n" + "INFO: boundary assembly : 1.677372932434082\n", + "INFO: field assembly : 1.1189799308776855\n", + "INFO: dump matrices to disk : 1.127087116241455\n", + "INFO: solve problem : 0.3656301498413086\n", + "INFO: update element data : 0.024611949920654297\n", + "INFO: non-linear iteration : 4.313704013824463\n", + "INFO: solver finished in 4.458410978317261 seconds.\n" ] } ], @@ -785,7 +627,43 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 84 1.000000 0.500000 0.500000\n", + "2 89 0.500000 0.500000 0.500000\n", + "3 69 0.500000 0.000000 0.500000\n", + "4 68 1.000000 0.000000 0.500000\n" + ] + } + ], + "source": [ + "prob = MortarProblem(\"temperature\", 1)\n", + "ENV[\"COLUMNS\"] = 160\n", + "ass = JuliaFEM.Core.Assembly()\n", + "el = contact_problem.elements[1]\n", + "push!(prob, el)\n", + "JuliaFEM.Core.assemble!(ass, prob, el, 0.0)\n", + "for (i, X) in enumerate(el(\"geometry\", 0.0))\n", + " ci = JuliaFEM.Core.get_connectivity(el)[i]\n", + " @printf \"%d %d %f %f %f\\n\" [i; ci; X]...\n", + "end\n", + "C = sparse(ass.stiffness_matrix)*46656\n", + "nz1 = sort(unique(rowvals(C)))\n", + "nz2 = sort(unique(rowvals(C')))\n", + "C = C[nz1,nz2]\n", + "CC = full(C);" + ] + }, + { + "cell_type": "code", + "execution_count": 9, "metadata": { "collapsed": false }, @@ -794,30 +672,140 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO: nid near corner = 82\n", - "INFO: displacement X = [1.0,1.0,1.0], u = [0.02865803993668297,0.02865803993668297,-0.11229776250350555]\n" + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n", + "INFO: win\n" + ] + } + ], + "source": [ + "for me in el[\"master elements\"]\n", + " \n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "JuliaFEM.Core.Field{JuliaFEM.Core.Discrete,JuliaFEM.Core.Variable,JuliaFEM.Core.TimeInvariant}(Array{T,2}[\n", + "3x3 Array{Float64,2}:\n", + " 0.0 1.0 0.0\n", + " 0.0 0.0 1.0\n", + " 1.0 0.0 0.0,\n", + "\n", + "3x3 Array{Float64,2}:\n", + " 0.0 1.0 0.0\n", + " 0.0 0.0 1.0\n", + " 1.0 0.0 0.0,\n", + "\n", + "3x3 Array{Float64,2}:\n", + " 0.0 1.0 0.0\n", + " 0.0 0.0 1.0\n", + " 1.0 0.0 0.0,\n", + "\n", + "3x3 Array{Float64,2}:\n", + " 0.0 1.0 0.0\n", + " 0.0 0.0 1.0\n", + " 1.0 0.0 0.0])" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "el(\"normal-tangential coordinates\", 0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "4x13 Array{Float64,2}:\n", + " -784.0 -700.0 -28.0 -700.0 -28.0 -23.7 -0.5 -622.1 -23.7 1295.91 647.289 647.289 319.511\n", + " -224.0 -1064.0 -224.0 -200.0 -8.0 -180.6 -2.2 -919.6 -27.6 647.289 1285.51 319.511 597.689\n", + " -224.0 -200.0 -8.0 -1064.0 -224.0 -27.6 -2.2 -919.6 -180.6 647.289 319.511 1285.51 597.689\n", + " -64.0 -304.0 -64.0 -304.0 -64.0 -186.6 -8.6 -1224.2 -186.6 319.511 597.689 597.689 891.111" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "CC" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: nid near corner = 73\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.027468071147715937,0.02746807114771595,-0.11069400634503247]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Test Failed\n", + " Expression: " ] }, { "ename": "LoadError", - "evalue": "LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14", + "evalue": "LoadError: There was an error during testing\nwhile loading In[24], in expression starting on line 14", "output_type": "error", "traceback": [ - "LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14", + "LoadError: There was an error during testing\nwhile loading In[24], in expression starting on line 14", "", - " in error at ./error.jl:21", - " in default_handler at test.jl:30", - " in do_test at test.jl:53", - " [inlined code] from In[11]:20", + " in record at /home/jukka/.julia/v0.4/BaseTestNext/src/BaseTestNext.jl:290", + " in do_test at /home/jukka/.julia/v0.4/BaseTestNext/src/BaseTestNext.jl:192", + " [inlined code] from In[24]:20", " in anonymous at no file:0" ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "isapprox(u[i],known_value)\n" + ] } ], "source": [ "using JuliaFEM.Test\n", "\n", "nid = 0\n", - "for (nid, coords) in mesh_upper[\"nodes\"]\n", + "for (nid, coords) in mesh[\"nodes\"]\n", " if isapprox(coords, [1.0, 1.0, 1.0])\n", " info(\"nid near corner = $nid\")\n", " break\n", @@ -833,13 +821,39 @@ " X = element(\"geometry\", 0.0)\n", " u = element(\"displacement\", 0.0)\n", " info(\"displacement X = $(X[i]), u = $(u[i])\")\n", - " @test isapprox(-u[i], known_value)\n", + " @test isapprox(u[i], known_value)\n", "end" ] }, { "cell_type": "code", - "execution_count": 28, + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Array{Float64,1}:\n", + " 0.0274681\n", + " 0.0274681\n", + " -0.110694 " + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "[0.02952918548952522,0.0295291854895252,-0.11404135697110807]\n", + "[0.027468071147715937,0.02746807114771595,-0.11069400634503247]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, "metadata": { "collapsed": false }, @@ -848,7 +862,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "INFO: XDFM: ndim = 369\n" + "INFO: XDFM: ndim = 273\n" ] } ], @@ -879,7 +893,7 @@ "#elements\n", "JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements)\n", "JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, \"displacement\", disp)\n", - "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/foobar3.xmf\");" + "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/blocks.xmf\");" ] }, { diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 1ad400a..d0e2017 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -6,6 +6,8 @@ This is JuliaFEM -- Finite Element Package """ module JuliaFEM +# include("common.jl") + """ JuliaFEM Core module. """ module Core include("core.jl") @@ -19,6 +21,10 @@ include("api.jl") end module Preprocess +macro debug(msg) + haskey(ENV, "DEBUG") || return + return msg +end include("abaqus_reader.jl") include("preprocess_aster_reader.jl") end diff --git a/src/common.jl b/src/common.jl new file mode 100644 index 0000000..3f1c525 --- /dev/null +++ b/src/common.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 + +# common routines + +""" +A very simple debugging macro. It executes commands if environment variable DEBUG is set. + +Usage +----- +Instead of starting session `julia file.jl`, do `DEBUG=1 julia file.jl`. +Or set `export DEBUG=1` for your `.bashrc`. + +Running inside code +------------------- + +julia> @debug info("moimoi heihei") + +will get executed iff environment variable DEBUG is set. + +Examples +-------- +julia> @debug info("moimoi") +(empty) +julia> ENV["DEBUG"] = 1 +julia> @debug info("moimoi") +INFO: moimoi +julia> @debug begin +... info("matrix is") +... dump([1 2; 3 4]) +... end +INFO: matrix is +Array(Int64(2,2)) 2x2 Array{Int64,2}: + 1 2 + 3 4 + +""" +macro debug(msg) + haskey(ENV, "DEBUG") || return + return msg +end + +function set_debug_on!() + ENV["DEBUG"] = 1; +end + +function set_debug_off!() + pop!(ENV, "DEBUG"); +end + +export @debug, set_debug_on!, set_debug_off! + diff --git a/src/core.jl b/src/core.jl index b88fbbd..61dd385 100644 --- a/src/core.jl +++ b/src/core.jl @@ -1,31 +1,10 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +using JuliaFEM + import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof, append! -""" -A very simple debugging macro. It prints debug message if environment variable -JULIAFEM_DEBUG is found. - -Usage: instead of starting session `julia file.jl` do `DEBUG=1 julia file.jl`. -Or set `export DEBUG=1` for your `.bashrc`. -""" -macro debug(msg) - haskey(ENV, "DEBUG") || return -# return :( println("DEBUG: ", $msg) ) - return msg -end - -function set_debug_on!() - ENV["DEBUG"] = 1; -end - -function set_debug_off!() - pop!(ENV, "DEBUG"); -end - -export @debug, set_debug_on!, set_debug_off! - using ForwardDiff autodiffcache = ForwardDiffCache() # export derivative, jacobian, hessian diff --git a/src/directsolver.jl b/src/directsolver.jl index 5a961bd..7bc533b 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -143,6 +143,7 @@ end function call(solver::DirectSolver, time::Number=0.0) info("# of field problems: $(length(solver.field_problems))") info("# of boundary problems: $(length(solver.boundary_problems))") + (length(solver.field_problems) != 0) || error("no field problems defined.") timing = Dict{ASCIIString, Float64}() tic(timing, "solver") diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 82b560c..853c09f 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -3,8 +3,11 @@ abstract DirichletProblem <: AbstractProblem -function DirichletProblem(parent_field_name, parent_field_dim, dim=1, elements=[]) - return BoundaryProblem{DirichletProblem}(parent_field_name, parent_field_dim, dim, elements) +function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]) + return BoundaryProblem{DirichletProblem}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements) +end +function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]) + return BoundaryProblem{DirichletProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements) end function assemble!(assembly::Assembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Number) diff --git a/src/elements.jl b/src/elements.jl index 58639f7..13bfe20 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -264,3 +264,13 @@ function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time:: element["normal-tangential coordinates"] = ntcoords end +""" Pick values from nodes and set to element according to connectivity. """ +function update(element::Element, field_name::ASCIIString, data::Union{Vector, Dict}) + element[field_name] = [data[i] for i in get_connectivity(element)] +end +function update(elements::Vector{Element}, field_name::ASCIIString, data::Union{Vector, Dict}) +# info("update $field_name for $(length(elements)) elements.") + for element in elements + update(element, field_name, data) + end +end diff --git a/src/integrate.jl b/src/integrate.jl index 8d6ba29..0601eb0 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -78,6 +78,19 @@ function get_integration_points(::TriangularElement, ::Type{Val{2}}) ] end +function get_integration_points(::TriangularElement, ::Type{Val{4}}) + # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF + # FIXME: something wrong here with weights ..? + [ + IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801), + IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801), + IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801), + IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532), + IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532), + IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532) + ] +end + function get_integration_points(::TriangularElement, ::Type{Val{5}}) # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF # FIXME: something wrong here with weights ..? diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl index 629efdc..96d9d14 100644 --- a/src/linear_elasticity.jl +++ b/src/linear_elasticity.jl @@ -5,8 +5,8 @@ abstract LinearElasticityProblem <: ElasticityProblem -function LinearElasticityProblem(dim::Int=3, elements=[]) - return Problem{LinearElasticityProblem}(dim, elements) +function LinearElasticityProblem(name="linear elasticity", dim::Int=3, elements=[]) + return Problem{LinearElasticityProblem}(name, dim, elements) end """ Elasticity equations, general 3D case. """ @@ -57,13 +57,22 @@ function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, proble L = w*T*N*norm(cross(JT[:,1], JT[:,2])) add!(assembly.force_vector, gdofs, vec(L)) end + for dim in 1:problem.dim + if haskey(element, "displacement traction force $dim") + T = element("displacement traction force $dim", ip, time) + ldofs = gdofs[dim:problem.dim:end] + JT = transpose(J) + L = w*T*N*norm(cross(JT[:,1], JT[:,2])) + add!(assembly.force_vector, ldofs, vec(L)) + end + end end end abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem -function PlaneStressLinearElasticityProblem(dim::Int=2, elements=[]) - return Problem{PlaneStressLinearElasticityProblem}(dim, elements) +function PlaneStressLinearElasticityProblem(name="plane stress linear elasticity", dim::Int=2, elements=[]) + return Problem{PlaneStressLinearElasticityProblem}(name, dim, elements) end """ Elasticity equations, plane stress. """ diff --git a/src/mortar.jl b/src/mortar.jl index 5834cf6..b241323 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -3,6 +3,13 @@ # Mortar projection calculation for 2d + +macro debug(msg) + haskey(ENV, "DEBUG") || return + return msg +end + + """ Find projection from slave nodes to master element, i.e. find xi2 from master element corresponding to the xi1. """ @@ -617,7 +624,18 @@ function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix return theta end end - error("project_point_to_auxiliary_plane: did not converge in $max_iterations iterations!") + begin + info("projecting point from auxiliary plane back to surface didn't go very well.") + info("element type: $E") + info("element connectivity: $(get_connectivity(element))") + info("auxiliary plane: x0 = $x0, Q = $Q") + info("point coordinates on plane: $p") + info("element geometry: $x") + info("ph: $ph") + info("normal direction: $n") + info("parameter vector before giving up: $theta") + end + error("project_point_to_surface: did not converge in $max_iterations iterations!") end @@ -632,8 +650,12 @@ node_csys """ abstract MortarProblem <: AbstractProblem -function MortarProblem(parent_field_name, parent_field_dim, dim=1, elements=[]) - return BoundaryProblem{MortarProblem}(parent_field_name, parent_field_dim, dim, elements) +function MortarProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[]) + return BoundaryProblem{MortarProblem}("mortar problem", parent_field_name, parent_field_dim, dim, elements) +end + +function MortarProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[]) + return BoundaryProblem{MortarProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements) end # Mortar assembly @@ -688,6 +710,31 @@ end typealias MortarElements3D Union{Tri3, Quad4} +""" Find master elements from list of potential master elements. """ +function find_master_elements(slave_element::Element, time::Real) + x0, Q = create_auxiliary_plane(slave_element, time) + Sl = Vector{Float64}[] + for p in slave_element("geometry", time) + push!(Sl, project_point_to_auxiliary_plane(p, x0, Q)) + end + S = hcat(Sl...) + master_elements = Element[] + + for master_element in slave_element["master elements"] + M = Vector{Float64}[] + for p in master_element("geometry", time) + push!(M, project_point_to_auxiliary_plane(p, x0, Q)) + end + M = hcat(M...) + P, neighbours = clip_polygon(S, M) + isa(P, Void) && continue # no clipping + size(P, 2) < 3 && continue # shared edge, no contribution + push!(master_elements, master_element) + end + + return master_elements +end + function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element{E}, time::Real) field_dim = problem.parent_field_dim field_name = problem.parent_field_name @@ -702,11 +749,14 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro for p in slave_element("geometry", time) push!(Sl, project_point_to_auxiliary_plane(p, x0, Q)) end - @debug info("auxiliary plane coords and basis: origo = $x0") - @debug info("basis:") - @debug dump(round(Q, 3)) + @debug begin + info("auxiliary plane coords and basis: origo = $x0") + info("basis:") + dump(round(Q, 3)) + end #S = reshape([S...;], 2, size(slave_element)[2]) S = hcat(Sl...) + slave_geom = Field(Vector{Float64}[S[:,j] for j=1:size(S,2)]) for master_element in slave_element["master elements"] master_dofs = get_gdofs(master_element, field_dim) @@ -717,11 +767,16 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro end #M = reshape([M...;], 2, size(master_element)[2]) M = hcat(M...) + master_geom = Field(Vector{Float64}[M[:,j] for j=1:size(M,2)]) + P = nothing neighbours = nothing - @debug info("applying polygon clip algorithm, S & M = ") - @debug dump(round(S, 3)) - @debug dump(round(M, 3)) + @debug begin + info("applying polygon clip algorithm, S & M = ") + dump(round(S, 3)) + dump(round(M, 3)) + end + try P, neighbours = clip_polygon(S, M) catch @@ -735,53 +790,114 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro error("cannot continue") end isa(P, Void) && continue # no clipping - @debug info("polygon on auxilyary plane: ") - @debug dump(round(P, 3)) - C = calculate_polygon_centerpoint(P) - @debug info("center point = $C") + @debug begin + info("polygon coords on auxilyary plane: ") + dump(round(P, 3)) + end + if size(P, 2) < 3 + # shared edge but no shared volume. skipping + continue + info("this is not polygon at all.") + info("clipping S") + dump(S) + info("clipping M") + dump(M) + error("size(P, 2) < 3") + end + C = calculate_polygon_centerpoint(P) npts = size(P, 2) # number of vertices in polygon - @debug info("number of vectices in polygon: $npts") -# S = zeros(3, 3) -# M = zeros(3, 3) + + @debug begin + info("clip polygon info") + theta = project_point_from_plane_to_surface(C, x0, Q, slave_element, time) + CC = slave_element("geometry", theta[2:3], time) + info("center point on slave: $CC") + info("number of vectices in polygon: $npts") + on_slave = zeros(3, 0) + on_master = zeros(3, 0) + for i=1:size(P, 2) + theta = project_point_from_plane_to_surface(P[:,i], x0, Q, slave_element, time) + on_slave = [on_slave slave_element("geometry", theta[2:3], time)] + theta = project_point_from_plane_to_surface(P[:,i], x0, Q, master_element, time) + on_master = [on_master master_element("geometry", theta[2:3], time)] + end + info("polygon coords projected to slave element") + dump(round(on_slave, 3)) + info("polygon coords projected to master element") + dump(round(on_master, 3)) + end + for i=1:npts # loop vertices and create temporary integrate cells xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]] yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]] X = hcat(xvec, yvec)' - @debug info("cell $i, coords = ") - @debug dump(round(X, 3)) - geom = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)]) + + @debug begin + on_slave = zeros(3, 0) + on_master = zeros(3, 0) + for j=1:size(X, 2) + theta = project_point_from_plane_to_surface(X[:,j], x0, Q, slave_element, time) + on_slave = [on_slave slave_element("geometry", theta[2:3], time)] + theta = project_point_from_plane_to_surface(X[:,j], x0, Q, master_element, time) + on_master = [on_master master_element("geometry", theta[2:3], time)] + end + info("cell $i coords projected to slave element") + dump(round(on_slave, 3)) + info("cell $i coords projected to master element") + dump(round(on_master, 3)) + end + + # integration cell geometry, i.e., Tri3 + cell = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)]) + # info("geom = $geom") for ip in get_integration_points(Tri3, Val{5}) - # calculate determiant of jacobian - #dN = get_dbasis(E, ip.xi) - dN = get_dbasis(Tri3, ip.xi) - J = sum([kron(dN[:,j], geom[j]') for j=1:length(geom)]) - w = ip.weight*det(J) # gauss point in auxiliary plane #N = get_basis(E, ip.xi) N = get_basis(Tri3, ip.xi) - x = vec(N*geom) + xi = vec(N*cell) # xi defined in auxilary plane + #xi = ip.xi + # info("x = $x") # find projection of gauss point to master and slave elements - theta1 = project_point_from_plane_to_surface(x, x0, Q, slave_element, time) - theta2 = project_point_from_plane_to_surface(x, x0, Q, master_element, time) + theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time) + theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time) + xi_slave = theta1[2:3] + xi_master = theta2[2:3] + @debug begin + X_slave = slave_element("geometry", xi_slave, time) + X_master = master_element("geometry", xi_master, time) + info("integration point on slave: $xi_slave => $X_slave") + info("integration point on master: $xi_master => $X_master") + end # evaluate shape functions values in gauss point and add contribution to matrices - N1 = slave_element(theta1[2:3], time) - N2 = master_element(theta2[2:3], time) - Sm = w*N1'*N1 + N1 = slave_element(xi_slave, time) + N2 = master_element(xi_master, time) + + # calculate determiant of jacobian + dNC = get_dbasis(Tri3, ip.xi) + dNS = get_dbasis(Quad4, xi_slave) + dNM = get_dbasis(Quad4, xi_master) + JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)]) + JN = sum([kron(dNS[:,j], slave_geom[j]') for j=1:length(slave_geom)]) + JM = sum([kron(dNM[:,j], master_geom[j]') for j=1:length(master_geom)]) + wS = det(JN) + wM = det(JM) + wC = det(JC) + @debug info("weight S = $wS, weight M = $wM, weight C = $wC") + + Sm = ip.weight*N1'*N1*wC # FIXME: master side transpose -- why? - Mm = w*(N1'*N2)' + Mm = ip.weight*(N1'*N2)'*wC for k=1:field_dim sd = slave_dofs[k:field_dim:end] md = master_dofs[k:field_dim:end] add!(assembly.stiffness_matrix, sd, sd, Sm) add!(assembly.stiffness_matrix, sd, md, -Mm) -# info("sd = $sd") -# info("md = $md") end end +# info("breaking on first") +# break end -# info("S = \n$S") -# info("M = \n$M") end end diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 003a6b5..c297813 100644 --- a/src/preprocess_aster_reader.jl +++ b/src/preprocess_aster_reader.jl @@ -1,6 +1,8 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +using JuliaFEM + using HDF5 function aster_parse_nodes(section::ASCIIString; strip_characters=true) @@ -53,6 +55,144 @@ function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}}) return model end + +function aster_renumber_nodes_!(mesh, node_numbering) + old_nodes = mesh["nodes"] + new_nodes = typeof(old_nodes)() + for (node_id, node_coords) in old_nodes + new_node_id = node_numbering[node_id] + new_nodes[new_node_id] = node_coords + end + mesh["nodes"] = new_nodes + for (elid, (eltype, elset, elcon)) in mesh["connectivity"] + new_elcon = [node_numbering[node_id] for node_id in elcon] + mesh["connectivity"][elid] = (eltype, elset, new_elcon) + end +end + + +function aster_renumber_nodes!(mesh1, mesh2) + + reserved_node_ids = Set(collect(keys(mesh1["nodes"]))) + @debug info("already reserved node ids: $reserved_node_ids") + mesh2_node_numbering = Dict{Int64, Int64}() + + # find new node ids assigned for mesh 2 + k = 1 + for node_id in sort(collect(keys(mesh2["nodes"]))) + @debug info("mesh2: processing node $node_id") + # if node id is reserved in mesh 1, find new number + if node_id in reserved_node_ids + @debug info("node id conflict, $node_id already defined in mesh 1, renumbering") + while k in reserved_node_ids + k += 1 + end + @debug info("mesh2: node $node_id -> $k") + mesh2_node_numbering[node_id] = k + push!(reserved_node_ids, k) + else + mesh2_node_numbering[node_id] = node_id + end + end + @debug info("new node numering:") + @debug println(mesh2_node_numbering) + aster_renumber_nodes_!(mesh2, mesh2_node_numbering) + +#= + # create new nodes + mesh2_old_nodes = mesh2["nodes"] + mesh2_new_nodes = typeof(mesh2_old_nodes)() + for (node_id, node_coords) in mesh2_old_nodes + new_node_id = mesh2_node_numbering[node_id] + mesh2_new_nodes[new_node_id] = node_coords + end + mesh2["nodes"] = mesh2_new_nodes + + # update connectivity + for (elid, (eltype, elset, elcon)) in mesh2["connectivity"] + new_elcon = [mesh2_node_numbering[node_id] for node_id in elcon] + mesh2["connectivity"][elid] = (eltype, elset, new_elcon) + end +=# + +end + + +function aster_renumber_elements!(mesh1, mesh2) + + reserved_element_ids = Set(collect(keys(mesh1["connectivity"]))) + @debug info("already reserved element ids: $reserved_element_ids") + mesh2_element_numbering = Dict{Int64, Int64}() + + # find new element ids assigned for mesh 2 + k = 1 + for element_id in sort(collect(keys(mesh2["connectivity"]))) + @debug info("mesh2: processing element $element_id") + # if node id is reserved in mesh 1, find new number + if element_id in reserved_element_ids + @debug info("element id conflict, $element_id already defined in mesh 1, renumbering") + while k in reserved_element_ids + k += 1 + end + @debug info("mesh2: element $element_id -> $k") + mesh2_element_numbering[element_id] = k + push!(reserved_element_ids, k) + else + mesh2_element_numbering[element_id] = element_id + end + end + @debug info("element numbering for mesh 2:") + @debug info(mesh2_element_numbering) + + # create new elements + mesh2_old_elements = mesh2["connectivity"] + mesh2_new_elements = typeof(mesh2_old_elements)() + for (element_id, element_data) in mesh2_old_elements + new_element_id = mesh2_element_numbering[element_id] + mesh2_new_elements[new_element_id] = element_data + end + mesh2["connectivity"] = mesh2_new_elements + +end + + +function aster_combine_meshes(mesh1, mesh2) + + # check that meshes are ready to be combined + node_ids_mesh_1 = collect(keys(mesh1["nodes"])) + node_ids_mesh_2 = collect(keys(mesh2["nodes"])) + if length(intersect(node_ids_mesh_1, node_ids_mesh_2)) != 0 + error("nodes with same id number in both meshes, failed.") + end + element_ids_mesh_1 = collect(keys(mesh1["connectivity"])) + element_ids_mesh_2 = collect(keys(mesh2["connectivity"])) + if length(intersect(element_ids_mesh_1, element_ids_mesh_2)) != 0 + error("elements with same id number in both meshes, failed.") + end + @assert similar(mesh1) == similar(mesh2) + @assert similar(mesh1["nodes"]) == similar(mesh2["nodes"]) + @assert similar(mesh1["connectivity"]) == similar(mesh2["connectivity"]) + + new_mesh = similar(mesh1) + new_mesh["nodes"] = similar(mesh1["nodes"]) + new_mesh["connectivity"] = similar(mesh1["connectivity"]) + + for (node_id, node_coords) in mesh1["nodes"] + new_mesh["nodes"][node_id] = node_coords + end + for (node_id, node_coords) in mesh2["nodes"] + new_mesh["nodes"][node_id] = node_coords + end + for (element_id, element_data) in mesh1["connectivity"] + new_mesh["connectivity"][element_id] = element_data + end + for (element_id, element_data) in mesh2["connectivity"] + new_mesh["connectivity"][element_id] = element_data + end + return new_mesh +end + + """ Code Aster binary file (.med), which is exported from SALOME. """ diff --git a/src/problems.jl b/src/problems.jl index 1021c3f..7091b15 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -4,11 +4,13 @@ abstract AbstractProblem type Problem{T<:AbstractProblem} + name :: ASCIIString dim :: Int elements :: Vector{Element} end type BoundaryProblem{T<:AbstractProblem} + name :: ASCIIString parent_field_name :: ASCIIString parent_field_dim :: Int dim :: Int @@ -36,3 +38,4 @@ end function Base.push!(problem::AllProblems, element::Element) push!(problem.elements, element) end + diff --git a/src/test.jl b/src/test.jl index 015b1b9..b8e2441 100644 --- a/src/test.jl +++ b/src/test.jl @@ -1,7 +1,7 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md -using Base.Test +using BaseTestNext abstract TestResult @@ -159,5 +159,4 @@ function print_test_statistics() return passed, failed, errors, critical end -export @test, run_test, print_test_statistics - +export @test, @testset, @test_throws, run_test, print_test_statistics diff --git a/test/test_mortar.jl b/test/test_mortar.jl index e7b1434..980f0fe 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -5,7 +5,8 @@ module MortarTests using JuliaFEM.Test -using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, MortarProblem, Assembly, assemble! +using JuliaFEM.Core: Element, Seg2, Quad4, Tri3, Hex8, MortarProblem, Assembly, assemble!, + get_connectivity, update using JuliaFEM.Core: PlaneStressElasticityProblem, DirichletProblem, DirectSolver # 2d stuff @@ -19,6 +20,7 @@ using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane, calculate_normal_tangential_coordinates!, is_point_inside_convex_polygon +using JuliaFEM.Core: LinearElasticityProblem function get_test_2d_model() # this is hand calculated and given as an example in my thesis @@ -691,6 +693,240 @@ function test_assemble_3d_problem_quad4_3() @test isapprox(stiffness_matrix, B) end -test_assemble_3d_problem_quad4_3() +#test_assemble_3d_problem_quad4_3() + + +function test_3d_problem() + nodes = Vector{Float64}[ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.5], + [1.0, 0.0, 0.5], + [1.0, 1.0, 0.5], + [0.0, 1.0, 0.5], + + [0.0, 0.0, 0.5], + [1.0, 0.0, 0.5], + [1.0, 1.0, 0.5], + [0.0, 1.0, 0.5], + [0.0, 0.0, 1.0], + [1.0, 0.0, 1.0], + [1.0, 1.0, 1.0], + [0.0, 1.0, 1.0], + ] + el1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8]) + el2 = Hex8([9, 10, 11, 12, 13, 14, 15, 16]) + sym121 = Quad4([1, 2, 3, 4]) + sym131 = Quad4([1, 2, 6, 5]) + sym132 = Quad4([9, 10, 14, 13]) + sym231 = Quad4([4, 1, 5, 8]) + sym232 = Quad4([12, 9, 13, 16]) + force = Quad4([14, 15, 16, 13]) + l2u = Quad4([5, 6, 7, 8]) + u2l = Quad4([9, 10, 11, 12]) + elements = Element[el1, el2, sym121, sym131, sym132, sym231, sym232, force, l2u, u2l] + update(elements, "geometry", nodes) + el1["youngs modulus"] = el2["youngs modulus"] = 900.0 + el1["poissons ratio"] = el2["poissons ratio"] = 0.25 + sym121["displacement 3"] = 0.0 + sym131["displacement 2"] = sym132["displacement 2"] = 0.0 + sym231["displacement 1"] = sym232["displacement 1"] = 0.0 + force["displacement traction force 3"] = -100.0 + l2u["master elements"] = Element[u2l] + calculate_normal_tangential_coordinates!(l2u, 0.0) + fb = LinearElasticityProblem("two elastic blocks") + push!(fb, el1, el2, force) + bc = DirichletProblem("symmetry boundaries", "displacement", 3) + push!(bc, sym121, sym131, sym132, sym231, sym232) + tie = MortarProblem("tie contact between bodies", "displacement", 3) + push!(tie, l2u) + solver = DirectSolver("solution of elasticity problem") + push!(solver, fb) + push!(solver, bc) + push!(solver, tie) + solver.nonlinear_problem = false + solver.method = :UMFPACK + call(solver, 0.0) + X = el2("geometry", [1.0, 1.0, 1.0], 0.0) + u = el2("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at $X = $u") + @test isapprox(u, 1/36*[1, 1, -4]) +end +#test_3d_problem() + +#= +@testset "plane quad4 projector tests" begin + a = 1/2 + b = 1/3 + nodes = Dict{Int64, Vector{Float64}}( + 1 => [0.0, 0.0, 0.0], + 2 => [1/2, 0.0, 0.0], + 3 => [1.0, 0.0, 0.0], + 4 => [0.0, 1.0, 0.0], + 5 => [1/2, 1.0, 0.0], + 6 => [1.0, 1.0, 0.0], + + 7 => [0.0, 0.0, 0.0], + 8 => [1/3, 0.0, 0.0], + 9 => [2/3, 0.0, 0.0], + 10 => [1.0, 0.0, 0.0], + 11 => [0.0, 1/2, 0.0], + 12 => [1/3, 1/2, 0.0], + 13 => [2/3, 1/2, 0.0], + 14 => [1.0, 1/2, 0.0], + 15 => [0.0, 1.0, 0.0], + 16 => [1/3, 1.0, 0.0], + 17 => [2/3, 1.0, 0.0], + 18 => [1.0, 1.0, 0.0], + ) + sel1 = Quad4([1, 2, 5, 4]) + sel2 = Quad4([2, 3, 6, 5]) + mel1 = Quad4([7, 8, 12, 11]) + mel2 = Quad4([8, 9, 13, 12]) + mel3 = Quad4([9, 10, 14, 13]) + mel4 = Quad4([11, 12, 16, 15]) + mel5 = Quad4([12, 13, 17, 16]) + mel6 = Quad4([13, 14, 18, 17]) + update(Element[sel1, sel2, mel1, mel2, mel3, mel4, mel5, mel6], "geometry", nodes) + calculate_normal_tangential_coordinates!(sel1, 0.0) + calculate_normal_tangential_coordinates!(sel2, 0.0) + prob = MortarProblem("temperature", 1) + push!(prob, sel1) + push!(prob, sel2) + + sel1["master elements"] = [mel1, mel2, mel4, mel5] + sel2["master elements"] = [mel2, mel3, mel5, mel6] + stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)*2592*6 + info("interface matrix:") + dump(round(stiffness_matrix, 3)) + B = [ + 864 432 0 432 216 0 -420 -375 -15 0 -504 -450 -18 0 -84 -75 -3 0 + 432 1728 432 216 864 216 -120 -690 -690 -120 -144 -828 -828 -144 -24 -138 -138 -24 + 0 432 864 0 216 432 0 -15 -375 -420 0 -18 -450 -504 0 -3 -75 -84 + 432 216 0 864 432 0 -84 -75 -3 0 -504 -450 -18 0 -420 -375 -15 0 + 216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120 + 0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420 + ] + info("expected interface matrix:") + dump(round(B, 3)) + @test isapprox(stiffness_matrix, B) +end +=# + +@testset "plane quad4 projector master 3x3 slave 2x2" begin + a = 1/2 + b = 1/3 + nodes = Dict{Int64, Vector{Float64}}( + 1 => [0*a, 0*a, 0.0], + 2 => [1*a, 0*a, 0.0], + 3 => [2*a, 0*a, 0.0], + 4 => [0*a, 1*a, 0.0], + 5 => [1*a, 1*a, 0.0], + 6 => [2*a, 1*a, 0.0], + 7 => [0*a, 2*a, 0.0], + 8 => [1*a, 2*a, 0.0], + 9 => [2*a, 2*a, 0.0], + + 10 => [0*b, 0*b, 0.0], + 11 => [1*b, 0*b, 0.0], + 12 => [2*b, 0*b, 0.0], + 13 => [3*b, 0*b, 0.0], + 14 => [0*b, 1*b, 0.0], + 15 => [1*b, 1*b, 0.0], + 16 => [2*b, 1*b, 0.0], + 17 => [3*b, 1*b, 0.0], + 18 => [0*b, 2*b, 0.0], + 19 => [1*b, 2*b, 0.0], + 20 => [2*b, 2*b, 0.0], + 21 => [3*b, 2*b, 0.0], + 22 => [0*b, 3*b, 0.0], + 23 => [1*b, 3*b, 0.0], + 24 => [2*b, 3*b, 0.0], + 25 => [3*b, 3*b, 0.0], + ) + sel1 = Quad4([1, 2, 5, 4]) + sel2 = Quad4([2, 3, 6, 5]) + sel3 = Quad4([4, 5, 8, 7]) + sel4 = Quad4([5, 6, 9, 8]) + mel1 = Quad4([10, 11, 15, 14]) + mel2 = Quad4([11, 12, 16, 15]) + mel3 = Quad4([12, 13, 17, 16]) + mel4 = Quad4([14, 15, 19, 18]) + mel5 = Quad4([15, 16, 20, 19]) + mel6 = Quad4([16, 17, 21, 20]) + mel7 = Quad4([18, 19, 23, 22]) + mel8 = Quad4([19, 20, 24, 23]) + mel9 = Quad4([20, 21, 25, 24]) + update(Element[sel1, sel2, sel3, sel4, mel1, mel2, mel3, + mel4, mel5, mel6, mel7, mel8, mel9], "geometry", nodes) + calculate_normal_tangential_coordinates!(sel1, 0.0) + calculate_normal_tangential_coordinates!(sel2, 0.0) + calculate_normal_tangential_coordinates!(sel3, 0.0) + calculate_normal_tangential_coordinates!(sel4, 0.0) + prob = MortarProblem("temperature", 1) + push!(prob, sel1) + push!(prob, sel2) + push!(prob, sel3) + push!(prob, sel4) + + master_elements = [mel1, mel2, mel3, mel4, mel5, mel6, mel7, mel8, mel9] + sel1["master elements"] = master_elements + sel2["master elements"] = master_elements + sel3["master elements"] = master_elements + sel4["master elements"] = master_elements + B = sparse(assemble(prob, 0.0).stiffness_matrix, 25, 25)*46656 + B = full(B) + D = B[1:9,1:9] + M = B[1:9,10:end] + info("interface matrix D:") + dump(round(D, 3)) + info("interface matrix M:") + dump(round(M, 3)) + D_expected = [ + 1296 648 0 648 324 0 0 0 0 + 648 2592 648 324 1296 324 0 0 0 + 0 648 1296 0 324 648 0 0 0 + 648 324 0 2592 1296 0 648 324 0 + 324 1296 324 1296 5184 1296 324 1296 324 + 0 324 648 0 1296 2592 0 324 648 + 0 0 0 648 324 0 1296 648 0 + 0 0 0 324 1296 324 648 2592 648 + 0 0 0 0 324 648 0 648 1296] + + M_expected = [ + -784 -700 -28 0 -700 -625 -25 0 -28 -25 -1 0 0 0 0 0 + -224 -1288 -1288 -224 -200 -1150 -1150 -200 -8 -46 -46 -8 0 0 0 0 + 0 -28 -700 -784 0 -25 -625 -700 0 -1 -25 -28 0 0 0 0 + -224 -200 -8 0 -1288 -1150 -46 0 -1288 -1150 -46 0 -224 -200 -8 0 + -64 -368 -368 -64 -368 -2116 -2116 -368 -368 -2116 -2116 -368 -64 -368 -368 -64 + 0 -8 -200 -224 0 -46 -1150 -1288 0 -46 -1150 -1288 0 -8 -200 -224 + 0 0 0 0 -28 -25 -1 0 -700 -625 -25 0 -784 -700 -28 0 + 0 0 0 0 -8 -46 -46 -8 -200 -1150 -1150 -200 -224 -1288 -1288 -224 + 0 0 0 0 0 -1 -25 -28 0 -25 -625 -700 0 -28 -700 -784] + + + info("D - D_expected") + dump(D - D_expected) + info("M - M_expected") + dump(M - M_expected) + + @test isapprox(D, D_expected) + @test isapprox(M, M_expected) + #= + B = [ + 864 432 0 432 216 0 -420 -375 -15 0 -504 -450 -18 0 -84 -75 -3 0 + 432 1728 432 216 864 216 -120 -690 -690 -120 -144 -828 -828 -144 -24 -138 -138 -24 + 0 432 864 0 216 432 0 -15 -375 -420 0 -18 -450 -504 0 -3 -75 -84 + 432 216 0 864 432 0 -84 -75 -3 0 -504 -450 -18 0 -420 -375 -15 0 + 216 864 216 432 1728 432 -24 -138 -138 -24 -144 -828 -828 -144 -120 -690 -690 -120 + 0 216 432 0 432 864 0 -3 -75 -84 0 -18 -450 -504 0 -15 -375 -420 + ] + info("expected interface matrix:") + dump(round(B, 3)) + @test isapprox(stiffness_matrix, B) +=# +end end diff --git a/test/test_preprocess_aster_reader.jl b/test/test_preprocess_aster_reader.jl index a1a0dd8..59134e2 100644 --- a/test/test_preprocess_aster_reader.jl +++ b/test/test_preprocess_aster_reader.jl @@ -7,7 +7,8 @@ using JuliaFEM using JuliaFEM.Test #using JuliaFEM: parse -using JuliaFEM.Preprocess: aster_parse_nodes +using JuliaFEM.Preprocess: aster_parse_nodes, aster_renumber_nodes!, + aster_renumber_elements!, aster_combine_meshes function test_read_mesh() mesh = """ @@ -69,6 +70,48 @@ function test_parse_nodes() @test nodes[8] == Float64[0.0, 1.0, 1.0] @test length(nodes) == 8 end -test_parse_nodes() +#test_parse_nodes() + +@testset "test combining meshes" begin + + mesh1 = Dict{ASCIIString, Any}( + "nodes" => Dict{Int64, Vector{Float64}}( + 1 => [1.0, 2.0], + 2 => [2.0, 3.0] + ), + "connectivity" => Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}( + 1 => (:SE2, :GRP1, [1, 2]) + ), + ) + + mesh2 = Dict{ASCIIString, Any}( + "nodes" => Dict{Int64, Vector{Float64}}( + 1 => [3.0, 4.0], + 2 => [4.0, 5.0] + ), + "connectivity" => Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}( + 1 => (:SE2, :GRP1, [1, 2]) + ), + ) + + aster_renumber_nodes!(mesh1, mesh2) + @test length(mesh2["nodes"]) == 2 + @test 3 in keys(mesh2["nodes"]) + @test 4 in keys(mesh2["nodes"]) + @test mesh2["connectivity"][1] == (:SE2, :GRP1, [3, 4]) + + aster_renumber_elements!(mesh1, mesh2) + @test length(mesh2["connectivity"]) == 1 + @test mesh2["connectivity"][2] == (:SE2, :GRP1, [3, 4]) + + mesh = aster_combine_meshes(mesh1, mesh2) + @test length(mesh["nodes"]) == 4 + @test mesh["nodes"][1] == [1.0, 2.0] + @test mesh["nodes"][2] == [2.0, 3.0] + @test mesh["nodes"][3] == [3.0, 4.0] + @test mesh["nodes"][4] == [4.0, 5.0] + @test mesh["connectivity"][1] == (:SE2, :GRP1, [1, 2]) + @test mesh["connectivity"][2] == (:SE2, :GRP1, [3, 4]) +end end