diff --git a/.gitignore b/.gitignore index eaebe9c..1795503 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ .ipynb_checkpoints docs/build/html *.swp +*.lnk diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index f5314b4..27f8ee4 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -25,31 +25,18 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n", - "WARNING: Base.String is deprecated, use AbstractString instead.\n" - ] - }, { "data": { "text/plain": [ - "Logger(root,DEBUG,PipeEndpoint(open, 0 bytes waiting),root)" + "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)" ] }, - "execution_count": 1, + "execution_count": 39, "metadata": {}, "output_type": "execute_result" } @@ -78,7 +65,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 40, "metadata": { "collapsed": false }, @@ -96,7 +83,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 41, "metadata": { "collapsed": false }, @@ -118,7 +105,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 42, "metadata": { "collapsed": false }, @@ -129,7 +116,7 @@ "MyQuad4" ] }, - "execution_count": 4, + "execution_count": 42, "metadata": {}, "output_type": "execute_result" } @@ -160,7 +147,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 43, "metadata": { "collapsed": false }, @@ -171,7 +158,7 @@ "get_element_dimension (generic function with 7 methods)" ] }, - "execution_count": 5, + "execution_count": 43, "metadata": {}, "output_type": "execute_result" } @@ -190,7 +177,7 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 44, "metadata": { "collapsed": false }, @@ -199,25 +186,8 @@ "name": "stderr", "output_type": "stream", "text": [ - "09-Oct 01:41:19:INFO:root:Testing element MyQuad4\n", - "09-Oct 01:41:19:INFO:root:number of basis functions in this element: 4\n", - "09-Oct 01:41:19:INFO:root:Initializing element\n", - "09-Oct 01:41:19:INFO:root:Element dimension: 2\n", - "09-Oct 01:41:19:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n", - "09-Oct 01:41:19:INFO:root:Pushing field to element.\n", - "09-Oct 01:41:19:INFO:root:Interpolating scalar field at [0.0,0.0]\n", - "09-Oct 01:41:19:INFO:root:Value: 2.5\n" + "20-loka 15:49:17:INFO:root:Testing element MyQuad4\n" ] - }, - { - "data": { - "text/plain": [ - "PipeEndpoint(open, 0 bytes waiting)" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" } ], "source": [ @@ -234,7 +204,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 45, "metadata": { "collapsed": false }, @@ -242,39 +212,41 @@ { "data": { "text/plain": [ - "2-element Array{JuliaFEM.Field{T},1}:\n", - " JuliaFEM.Field{Int64}(0.0,1,2)\n", - " JuliaFEM.Field{Int64}(1.0,1,3)" + "JuliaFEM.FieldSet(symbol(\"heat coefficient\"),JuliaFEM.Field[JuliaFEM.Field{Int64}(0.0,1,2),JuliaFEM.Field{Int64}(1.0,1,3)])" ] }, - "execution_count": 7, + "execution_count": 45, "metadata": {}, "output_type": "execute_result" - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "09-Oct 01:41:19:INFO:root:Element MyQuad4 passed tests.\n" - ] } ], "source": [ - "using JuliaFEM: new_fieldset!, add_field!, interpolate, dinterpolate\n", - "el1 = MyQuad4([1, 2, 3, 4])\n", - "new_fieldset!(el1, \"temperature\")\n", - "add_field!(el1, \"temperature\", Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n", - "add_field!(el1, \"temperature\", Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n", - "new_fieldset!(el1, \"geometry\")\n", - "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]))\n", - "new_fieldset!(el1, \"heat coefficient\")\n", - "add_field!(el1, \"heat coefficient\", Field(0.0, 2))\n", - "add_field!(el1, \"heat coefficient\", Field(1.0, 3))" + "using JuliaFEM: FieldSet, Field, interpolate, dinterpolate\n", + "element = MyQuad4([1, 2, 3, 4])\n", + "\n", + "geometry_field = Field(0.0, Vector[]) # Create empty field at time t=0.0\n", + "push!(geometry_field, [ 0.0, 0.0, 0.0]) # push some values for field\n", + "push!(geometry_field, [10.0, 0.0, 0.0])\n", + "push!(geometry_field, [10.0, 1.0, 0.0])\n", + "push!(geometry_field, [ 0.0, 1.0, 0.0])\n", + "geometry_fieldset = FieldSet(\"geometry\") # create fieldset \"geometry\"\n", + "push!(geometry_fieldset, geometry_field) # add field to fieldset\n", + "push!(element, geometry_fieldset) # add fieldset to element\n", + "\n", + "temperature_fieldset = FieldSet(\"temperature\")\n", + "push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n", + "push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n", + "push!(element, temperature_fieldset)\n", + "\n", + "heat_coefficient_fieldset = FieldSet(\"heat coefficient\")\n", + "push!(heat_coefficient_fieldset, Field(0.0, 2))\n", + "push!(heat_coefficient_fieldset, Field(1.0, 3))\n", + "push!(element, heat_coefficient_fieldset)" ] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 46, "metadata": { "collapsed": false }, @@ -285,19 +257,19 @@ "1.25" ] }, - "execution_count": 8, + "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# temperature at the middle poinf of the element, 1/4*(1+2+3+4) at t=0.5\n", - "interpolate(el1, \"temperature\", [0.0, 0.0], 0.5)" + "interpolate(element, \"temperature\", [0.0, 0.0], 0.5)" ] }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 47, "metadata": { "collapsed": false }, @@ -311,19 +283,19 @@ " 0.0" ] }, - "execution_count": 9, + "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# geometry midpoint of element\n", - "interpolate(el1, \"geometry\", [0.0, 0.0], 0.0)" + "interpolate(element, \"geometry\", [0.0, 0.0], 0.0)" ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 48, "metadata": { "collapsed": false }, @@ -337,19 +309,19 @@ " 0.0 0.0" ] }, - "execution_count": 10, + "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# interpolate derivatives works too\n", - "dinterpolate(el1, \"geometry\", [0.0, 0.0], 0.0)" + "dinterpolate(element, \"geometry\", [0.0, 0.0], 0.0)" ] }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 49, "metadata": { "collapsed": false }, @@ -360,19 +332,19 @@ "2.5" ] }, - "execution_count": 11, + "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# interpolating scalar -> scalar.\n", - "interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5)" + "interpolate(element, \"heat coefficient\", [0.0, 0.0], 0.5)" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 50, "metadata": { "collapsed": false }, @@ -383,13 +355,13 @@ "3" ] }, - "execution_count": 12, + "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)" + "interpolate(element, \"heat coefficient\", [0.0, 0.0], Inf)" ] }, { @@ -434,7 +406,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 51, "metadata": { "collapsed": false }, @@ -442,20 +414,21 @@ { "data": { "text/plain": [ - "get_unknown_field_name (generic function with 1 method)" + "get_unknown_field_name (generic function with 4 methods)" ] }, - "execution_count": 13, + "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: Equation, IntegrationPoint, Quad4\n", + "using JuliaFEM: Equation, IntegrationPoint, Quad4, get_unknown_field_name\n", "\n", "abstract Heat <: Equation\n", "\n", - "get_unknown_field_name(eq::Heat) = symbol(\"temperature\")" + "# it's important to define for which fieldset to save unknown values when solving\n", + "JuliaFEM.get_unknown_field_name(eq::Heat) = symbol(\"temperature\")" ] }, { @@ -467,7 +440,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 52, "metadata": { "collapsed": false }, @@ -492,7 +465,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 53, "metadata": { "collapsed": false }, @@ -503,20 +476,20 @@ "DC2D4" ] }, - "execution_count": 15, + "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "function DC2D4(el::Quad4)\n", + "function DC2D4(element::Quad4)\n", " integration_points = [\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n", " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n", - " new_fieldset!(el, \"temperature\") # assign new field \"temperature\" to element\n", - " DC2D4(el, integration_points, [])\n", + " push!(element, FieldSet(\"temperature\"))\n", + " DC2D4(element, integration_points, [])\n", "end" ] }, @@ -529,7 +502,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 54, "metadata": { "collapsed": false }, @@ -537,10 +510,10 @@ { "data": { "text/plain": [ - "has_lhs (generic function with 2 methods)" + "has_lhs (generic function with 4 methods)" ] }, - "execution_count": 16, + "execution_count": 54, "metadata": {}, "output_type": "execute_result" } @@ -569,7 +542,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 55, "metadata": { "collapsed": false }, @@ -584,25 +557,28 @@ " -1.0 -2.0 -1.0 4.0" ] }, - "execution_count": 17, + "execution_count": 55, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "using JuliaFEM: integrate, integrate_lhs, integrate_rhs, new_fieldset!, add_field!\n", - "el = Quad4([1, 2, 3, 4])\n", - "new_fieldset!(el, \"geometry\")\n", - "add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]))\n", - "new_fieldset!(el, \"temperature thermal conductivity\")\n", - "add_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n", - "eq = DC2D4(el)\n", - "integrate_lhs(eq, 1.0)" + "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n", + "element = Quad4([1, 2, 3, 4])\n", + "fieldset1 = FieldSet(\"geometry\")\n", + "field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n", + "push!(fieldset1, field1)\n", + "fieldset2 = FieldSet(\"temperature thermal conductivity\")\n", + "push!(fieldset2, Field(0.0, 6.0))\n", + "push!(element, fieldset1)\n", + "push!(element, fieldset2)\n", + "equation = DC2D4(element)\n", + "integrate_lhs(equation, 1.0)" ] }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 56, "metadata": { "collapsed": false }, @@ -613,13 +589,13 @@ "true" ] }, - "execution_count": 18, + "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "JuliaFEM.has_lhs(eq)" + "JuliaFEM.has_lhs(equation)" ] }, { @@ -631,7 +607,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 57, "metadata": { "collapsed": false }, @@ -642,13 +618,13 @@ "true" ] }, - "execution_count": 19, + "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "isa(integrate_rhs(eq, 1.0), Void)" + "integrate_rhs(equation, 1.0) == nothing" ] }, { @@ -660,7 +636,7 @@ }, { "cell_type": "code", - "execution_count": 20, + "execution_count": 58, "metadata": { "collapsed": false }, @@ -668,10 +644,10 @@ { "data": { "text/plain": [ - "has_rhs (generic function with 2 methods)" + "has_rhs (generic function with 4 methods)" ] }, - "execution_count": 20, + "execution_count": 58, "metadata": {}, "output_type": "execute_result" } @@ -688,11 +664,11 @@ " global_dofs :: Array{Int64, 1}\n", "end\n", "\n", - "function DC2D2(el::Seg2)\n", + "function DC2D2(element::Seg2)\n", " integration_points = [\n", - " IntegrationPoint([0], 2.0)]\n", - " new_fieldset!(el, \"temperature\")\n", - " DC2D2(el, integration_points, [])\n", + " IntegrationPoint([0.0], 2.0)]\n", + " push!(element, FieldSet(\"temperature\"))\n", + " DC2D2(element, integration_points, [])\n", "end\n", "\n", "\"\"\"\n", @@ -710,7 +686,7 @@ }, { "cell_type": "code", - "execution_count": 21, + "execution_count": 59, "metadata": { "collapsed": false }, @@ -723,19 +699,17 @@ " 50.0" ] }, - "execution_count": 21, + "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "el = Seg2([1, 2])\n", - "new_fieldset!(el, \"geometry\")\n", - "add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n", - "new_fieldset!(el, \"temperature flux\")\n", - "add_field!(el, \"temperature flux\", Field(0.0, 100.0))\n", - "eq = DC2D2(el)\n", - "integrate_rhs(eq, 1.0)" + "element = Seg2([1, 2])\n", + "push!(element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0,0.0], [0.0,1.0]])]))\n", + "push!(element, FieldSet(\"temperature flux\", [Field(0.0, 100.0)]))\n", + "equation = DC2D2(element)\n", + "integrate_rhs(equation, 1.0)" ] }, { @@ -749,7 +723,7 @@ }, { "cell_type": "code", - "execution_count": 22, + "execution_count": 60, "metadata": { "collapsed": false }, @@ -760,7 +734,7 @@ "PlaneHeatProblem" ] }, - "execution_count": 22, + "execution_count": 60, "metadata": {}, "output_type": "execute_result" } @@ -776,7 +750,7 @@ }, { "cell_type": "code", - "execution_count": 23, + "execution_count": 61, "metadata": { "collapsed": false }, @@ -784,10 +758,10 @@ { "data": { "text/plain": [ - "get_equation (generic function with 3 methods)" + "get_equation (generic function with 6 methods)" ] }, - "execution_count": 23, + "execution_count": 61, "metadata": {}, "output_type": "execute_result" } @@ -807,7 +781,7 @@ }, { "cell_type": "code", - "execution_count": 24, + "execution_count": 62, "metadata": { "collapsed": false }, @@ -816,14 +790,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "09-Oct 01:41:22:DEBUG:root:total dofs: 4\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "dofmap: Dict(4=>[3],2=>[1],3=>[2],5=>[4])\n" + "20-loka 15:49:17:DEBUG:root:total dofs: 4\n" ] }, { @@ -834,114 +801,40 @@ "\t[2, 1] = 300.0" ] }, - "execution_count": 24, + "execution_count": 62, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using JuliaFEM: get_connectivity, set_global_dofs!, get_global_dofs\n", - "using JuliaFEM: add_element!, get_equations, get_matrix_dimension\n", + "using JuliaFEM: add_element!, get_equations, get_matrix_dimension, calculate_global_dofs\n", + "using JuliaFEM: assign_global_dofs!, get_lhs, get_rhs\n", "\n", "# create elements and add necessary properties like connectivity and geometry\n", "el1 = Quad4([2, 3, 4, 5])\n", - "new_fieldset!(el1, \"geometry\")\n", - "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n", - "new_fieldset!(el1, \"temperature thermal conductivity\")\n", - "add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n", + "fs1geom = FieldSet(\"geometry\")\n", + "push!(fs1geom, Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n", + "fs1temp = FieldSet(\"temperature thermal conductivity\")\n", + "push!(fs1temp, Field(0.0, 6.0))\n", + "push!(el1, fs1geom)\n", + "push!(el1, fs1temp)\n", + "\n", "el2 = Seg2([2, 3])\n", - "new_fieldset!(el2, \"geometry\")\n", - "add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", - "new_fieldset!(el2, \"temperature flux\")\n", - "add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n", + "fs2geom = FieldSet(\"geometry\")\n", + "push!(fs2geom, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", + "fs2temp = FieldSet(\"temperature flux\")\n", + "push!(fs2temp, Field(1.0, 600.0))\n", + "push!(el2, fs2geom)\n", + "push!(el2, fs2temp)\n", "\n", "problem = PlaneHeatProblem()\n", - "add_element!(problem, el1)\n", - "add_element!(problem, el2)\n", - "\n", - "function JuliaFEM.get_connectivity(pr::Problem)\n", - " conn = Int[]\n", - " for eq in get_equations(pr)\n", - " el = get_element(eq)\n", - " append!(conn, get_connectivity(el))\n", - " end\n", - " conn = unique(conn)\n", - " return conn\n", - "end\n", - "\n", - "\"\"\"\n", - "Calculate global dofs for equations, maybe using some bandwidth\n", - "minimizing or fill reducing algorithm\n", - "\"\"\"\n", - "function calculate_global_dofs(pr::Problem)\n", - " conn = get_connectivity(pr)\n", - " dim = get_dimension(typeof(pr))\n", - " ndofs = dim*length(conn)\n", - " Logging.debug(\"total dofs: $ndofs\")\n", - "\n", - " mconn = maximum(conn)\n", - " gdofs = reshape(collect(1:mconn), dim, mconn)\n", - " dofmap = Dict{Int64, Array{Int64, 1}}()\n", - " for (i, c) in enumerate(conn)\n", - " dofmap[c] = gdofs[:, i]\n", - " end\n", - " return dofmap\n", - "end\n", - "\n", - "\"\"\"\n", - "Assign global dofs for equations.\n", - "\"\"\"\n", - "function assign_global_dofs!(pr::Problem, dofmap)\n", - " for eq in get_equations(pr)\n", - " el = get_element(eq)\n", - " c = get_connectivity(el)\n", - " #gdofs = [dofmap[ci] for ci in c]\n", - " gdofs = Int64[]\n", - " for ci in c\n", - " append!(gdofs, dofmap[ci])\n", - " end\n", - " set_global_dofs!(eq, gdofs)\n", - " end\n", - "end\n", + "push!(problem, el1)\n", + "push!(problem, el2)\n", "\n", "dofmap = calculate_global_dofs(problem)\n", - "println(\"dofmap: $dofmap\")\n", "assign_global_dofs!(problem, dofmap)\n", "\n", - "function get_lhs(pr::Problem, t::Float64)\n", - " I = Int64[]\n", - " J = Int64[]\n", - " V = Float64[]\n", - " dim = get_dimension(typeof(pr))\n", - " for eq in filter(has_lhs, get_equations(pr))\n", - " dofs = get_global_dofs(eq)\n", - " lhs = integrate_lhs(eq, t)\n", - " for (li, i) in enumerate(dofs)\n", - " for (lj, j) in enumerate(dofs)\n", - " push!(I, i)\n", - " push!(J, j)\n", - " push!(V, lhs[li, lj])\n", - " end\n", - " end\n", - " end\n", - " return I, J, V\n", - "end\n", - "\n", - "function get_rhs(pr::Problem, t::Float64)\n", - " I = Int64[]\n", - " V = Float64[]\n", - " dim = get_dimension(typeof(pr))\n", - " for eq in filter(has_rhs, get_equations(pr))\n", - " dofs = get_global_dofs(eq)\n", - " rhs = integrate_rhs(eq, t)\n", - " for (li, i) in enumerate(dofs)\n", - " push!(I, i)\n", - " push!(V, rhs[li])\n", - " end\n", - " end\n", - " return I, V\n", - "end\n", - "\n", "t = 1.0\n", "A = sparse(get_lhs(problem, t)...)\n", "b = sparsevec(get_rhs(problem, t)..., size(A, 1))" @@ -949,7 +842,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 63, "metadata": { "collapsed": false }, @@ -962,7 +855,7 @@ "\t[2, 1] = 100.0" ] }, - "execution_count": 25, + "execution_count": 63, "metadata": {}, "output_type": "execute_result" } @@ -984,7 +877,7 @@ }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 64, "metadata": { "collapsed": false }, @@ -992,118 +885,29 @@ { "data": { "text/plain": [ - "DirichletProblem" + "1-element Array{JuliaFEM.DirichletEquation,1}:\n", + " JuliaFEM.DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.FieldSet(symbol(\"reaction force\"),JuliaFEM.Field[]),:geometry=>JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])]))),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)" ] }, - "execution_count": 26, + "execution_count": 64, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "abstract BoundaryProblem <: Problem\n", + "using JuliaFEM: DirichletProblem\n", "\n", - "type DirichletProblem <: BoundaryProblem\n", - " equations :: Array{Equation, 1}\n", - "end\n", - "DirichletProblem() = DirichletProblem([])" - ] - }, - { - "cell_type": "code", - "execution_count": 27, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "get_equation (generic function with 4 methods)" - ] - }, - "execution_count": 27, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "abstract DirichletBC <: Equation\n", - "\n", - "get_unknown_field_name(eq::DirichletBC) = symbol(\"reaction force\")\n", - "\n", - "\"\"\"\n", - "Dirichlet boundary condition element for 2 node line segment\n", - "(plane problems).\n", - "\"\"\"\n", - "type DBC2D2 <: DirichletBC\n", - " element :: Seg2\n", - " integration_points :: Array{IntegrationPoint, 1}\n", - " global_dofs :: Array{Int64, 1}\n", - " fieldval :: Function\n", - "end\n", - "function DBC2D2(el::Seg2)\n", - " integration_points = [\n", - " IntegrationPoint([-sqrt(1/3)], 1.0),\n", - " IntegrationPoint([+sqrt(1/3)], 1.0)]\n", - " new_fieldset!(el, \"reaction force\")\n", - " fieldval(X, t) = 0.0\n", - " DBC2D2(el, integration_points, [], fieldval)\n", - "end\n", - "\n", - "function JuliaFEM.get_lhs(eq::DBC2D2, ip, t)\n", - " el = get_element(eq)\n", - " h = get_basis(el)(ip.xi)\n", - " return h*h'\n", - "end\n", - "\n", - "function JuliaFEM.get_rhs(eq::DBC2D2, ip, t)\n", - " el = get_element(eq)\n", - " h = get_basis(el, ip.xi)\n", - " f = eq.fieldval\n", - " X = interpolate(el, \"geometry\", ip.xi, t)\n", - " return h*f(X, t)\n", - "end\n", - "JuliaFEM.has_lhs(eq::DBC2D2) = true\n", - "JuliaFEM.has_rhs(eq::DBC2D2) = true\n", - "JuliaFEM.get_dimension(pr::Type{DirichletProblem}) = 1\n", - "JuliaFEM.get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2" - ] - }, - { - "cell_type": "code", - "execution_count": 28, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "1-element Array{JuliaFEM.Equation,1}:\n", - " DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[],:geometry=>JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)" - ] - }, - "execution_count": 28, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ "# create elements and add necessary properties like connectivity and geometry\n", "el3 = Seg2([4, 5])\n", - "new_fieldset!(el3, \"geometry\")\n", - "add_field!(el3, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", + "push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]])]))\n", "\n", "bc1 = DirichletProblem()\n", - "#get_equation(typeof(bc1), typeof(el3))\n", - "#DBC2D2(el3)\n", - "add_element!(bc1, el3)" + "push!(bc1, el3)" ] }, { "cell_type": "code", - "execution_count": 29, + "execution_count": 65, "metadata": { "collapsed": false }, @@ -1116,7 +920,7 @@ "\t[4, 1] = 0.0" ] }, - "execution_count": 29, + "execution_count": 65, "metadata": {}, "output_type": "execute_result" } @@ -1141,7 +945,7 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 66, "metadata": { "collapsed": false }, @@ -1156,7 +960,7 @@ " -1.0 -2.0 -1.0 4.0" ] }, - "execution_count": 30, + "execution_count": 66, "metadata": {}, "output_type": "execute_result" } @@ -1167,7 +971,7 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 67, "metadata": { "collapsed": false }, @@ -1182,7 +986,7 @@ " 0.0" ] }, - "execution_count": 31, + "execution_count": 67, "metadata": {}, "output_type": "execute_result" } @@ -1193,7 +997,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 68, "metadata": { "collapsed": false }, @@ -1208,7 +1012,7 @@ " 0.0 0.0 0.166667 0.333333" ] }, - "execution_count": 32, + "execution_count": 68, "metadata": {}, "output_type": "execute_result" } @@ -1219,7 +1023,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 69, "metadata": { "collapsed": false }, @@ -1234,7 +1038,7 @@ " 0.0" ] }, - "execution_count": 33, + "execution_count": 69, "metadata": {}, "output_type": "execute_result" } @@ -1245,7 +1049,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 70, "metadata": { "collapsed": false }, @@ -1278,7 +1082,7 @@ "\t[4, 8] = 0.333333" ] }, - "execution_count": 34, + "execution_count": 70, "metadata": {}, "output_type": "execute_result" } @@ -1289,7 +1093,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 71, "metadata": { "collapsed": false }, @@ -1304,7 +1108,7 @@ "\t[8, 1] = 0.0" ] }, - "execution_count": 35, + "execution_count": 71, "metadata": {}, "output_type": "execute_result" } @@ -1322,7 +1126,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 72, "metadata": { "collapsed": false }, @@ -1341,7 +1145,7 @@ " 0.0 0.0 0.166667 0.333333 0.0 0.0 0.0 0.0 " ] }, - "execution_count": 36, + "execution_count": 72, "metadata": {}, "output_type": "execute_result" } @@ -1352,7 +1156,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 73, "metadata": { "collapsed": false }, @@ -1361,7 +1165,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Non-zero rows: [1,2,3,4,7,8]\n" + "Non-zero rows: [1,2,3,4,7,8]" ] }, { @@ -1378,7 +1182,7 @@ " 600.0" ] }, - "execution_count": 37, + "execution_count": 73, "metadata": {}, "output_type": "execute_result" } @@ -1409,7 +1213,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 74, "metadata": { "collapsed": false }, @@ -1417,56 +1221,27 @@ { "data": { "text/plain": [ - "add_problem! (generic function with 1 method)" + "call (generic function with 1270 methods)" ] }, - "execution_count": 38, + "execution_count": 74, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "abstract Solver\n", + "using JuliaFEM: Solver, get_problems\n", "\n", - "\"\"\"\n", - "Simple solver for educational purposes.\n", - "\"\"\"\n", + "\"\"\" Simple solver for educational purposes. \"\"\"\n", "type SimpleSolver <: Solver\n", " problems\n", "end\n", "\n", - "\"\"\"\n", - "Default initializer\n", - "\"\"\"\n", + "\"\"\" Default initializer. \"\"\"\n", "function SimpleSolver()\n", " SimpleSolver(Problem[])\n", "end\n", "\n", - "function add_problem!(solver::Solver, problem::Problem)\n", - " push!(solver.problems, problem)\n", - "end" - ] - }, - { - "cell_type": "code", - "execution_count": 39, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "call (generic function with 1253 methods)" - ] - }, - "execution_count": 39, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "get_problems(s::Solver) = s.problems\n", "\"\"\"\n", "Call solver to solve a set of problems.\n", "\n", @@ -1516,7 +1291,7 @@ " element = get_element(equation)\n", " field_name = get_unknown_field_name(equation) # field we are solving\n", " field = Field(t, full(x1[gdofs])[:])\n", - " add_field!(element, field_name, field)\n", + " push!(element[field_name], field)\n", " end\n", "\n", " # update field for elements in problem 2\n", @@ -1525,7 +1300,7 @@ " element = get_element(equation)\n", " field_name = get_unknown_field_name(equation)\n", " field = Field(t, full(x2[gdofs]))\n", - " add_field!(element, field_name, field)\n", + " push!(element[field_name], field)\n", " end\n", "end" ] @@ -1541,41 +1316,7 @@ }, { "cell_type": "code", - "execution_count": 40, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/plain": [ - "5-element Array{Array{Float64,1},1}:\n", - " [-1.0,0.0]\n", - " [-0.5,0.5]\n", - " [0.0,1.0] \n", - " [0.5,1.5] \n", - " [1.0,2.0] " - ] - }, - "execution_count": 40, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "\"\"\"\n", - "Simple linspace extension to multidimensional values.\n", - "\"\"\"\n", - "function Base.linspace(X1, X2, n)\n", - " [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]\n", - "end\n", - "\n", - "linspace([-1.0, 0.0], [1.0, 2.0], 5)" - ] - }, - { - "cell_type": "code", - "execution_count": 41, + "execution_count": 75, "metadata": { "collapsed": false }, @@ -1584,59 +1325,45 @@ "name": "stderr", "output_type": "stream", "text": [ - "09-Oct 01:41:25:DEBUG:root:total dofs: 4\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Residual norm: 9.845568954283847e-14\n", - "Array{T,1}[[0.0,0.0],[0.25,0.0],[0.5,0.0],[0.75,0.0],[1.0,0.0]]\n", - "[100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003]\n" + "20-loka 15:49:17:DEBUG:root:total dofs: 4\n" ] } ], "source": [ - "using JuliaFEM: get_fieldset\n", - "\n", "# Define Problem 1:\n", "# - Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}\n", "# - Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁\n", "\n", "el1 = Quad4([1, 2, 3, 4])\n", - "new_fieldset!(el1, \"geometry\")\n", - "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n", - "new_fieldset!(el1, \"temperature thermal conductivity\")\n", - "add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n", + "push!(el1, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))\n", + "push!(el1, FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)]))\n", "\n", "el2 = Seg2([1, 2])\n", - "new_fieldset!(el2, \"geometry\")\n", - "add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]]))\n", + "push!(el2, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n", "\n", "# Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n", - "new_fieldset!(el2, \"temperature flux\")\n", - "add_field!(el2, \"temperature flux\", Field(0.0, 0.0))\n", - "add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n", + "load = FieldSet(\"temperature flux\")\n", + "push!(load, Field(0.0, 0.0))\n", + "push!(load, Field(1.0, 600.0))\n", + "push!(el2, load)\n", "\n", "problem1 = PlaneHeatProblem()\n", - "add_element!(problem1, el1)\n", - "add_element!(problem1, el2)\n", + "push!(problem1, el1)\n", + "push!(problem1, el2)\n", "\n", "# Define Problem 2:\n", "# - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n", "\n", "el3 = Seg2([3, 4])\n", - "new_fieldset!(el3, \"geometry\")\n", - "add_field!(el3, \"geometry\", Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]]))\n", + "push!(el3, FieldSet(\"geometry\", [Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]])]))\n", "\n", "problem2 = DirichletProblem()\n", - "add_element!(problem2, el3)\n", + "push!(problem2, el3)\n", "\n", "# Create a solver for a set of problems\n", "solver = SimpleSolver()\n", - "add_problem!(solver, problem1)\n", - "add_problem!(solver, problem2)\n", + "push!(solver, problem1)\n", + "push!(solver, problem2)\n", "\n", "# Solve problem at time t=1.0 and update fields\n", "call(solver, 1.0)\n", @@ -1650,6 +1377,41 @@ "println(T)" ] }, + { + "cell_type": "code", + "execution_count": 76, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Success :: (line:-1) :: fact was true\n", + " Expression: mean(T) --> roughly(100.0)\n", + " Expected: 100.0\n", + " Occurred: 100.00000000000003" + ] + }, + "execution_count": 76, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Residual norm: 9.845568954283847e-14\n", + "Array{T,1}[[0.0,0.0],[0.25,0.0],[0.5,0.0],[0.75,0.0],[1.0,0.0]]\n", + "[100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003]\n" + ] + } + ], + "source": [ + "@fact mean(T) --> roughly(100.0)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -1668,20 +1430,11 @@ "\n", "Any question about data structures, better ideas and improvements are very welcome." ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Julia 0.4.0-rc2", + "display_name": "Julia 0.4.0", "language": "julia", "name": "julia-0.4" }, diff --git a/notebooks/2015-09-15-2d-segmentation.ipynb b/notebooks/2015-09-15-2d-segmentation.ipynb index 8bfe810..8c5f147 100644 --- a/notebooks/2015-09-15-2d-segmentation.ipynb +++ b/notebooks/2015-09-15-2d-segmentation.ipynb @@ -8,7 +8,7 @@ "\n", "Author(s): Jukka Aho\n", "\n", - "**Abstract**: Contact segmentation in 2d" + "**Abstract**: Calculate mortar projection for 2d" ] }, { @@ -17,19 +17,10 @@ "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "WARNING: could not import Base.help into PyCall\n" - ] - } - ], + "outputs": [], "source": [ - "using JuliaFEM: PSeg, get_field, set_field, interpolate, new_field!, push_field!\n", - "using JuliaFEM: calculate_normals!, average_normals!, fit_derivative_field!\n", - "using JuliaFEM: set_degree, get_number_of_basis_functions, dinterpolate\n", + "using JuliaFEM: new_fieldset!, add_field!, dinterpolate\n", + "using JuliaFEM: Field, FieldSet, Seg2\n", "using PyPlot\n", "using ForwardDiff" ] @@ -38,12 +29,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Create some test boundaries" + "Create some test boundaries:" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 4, "metadata": { "collapsed": false }, @@ -54,7 +45,7 @@ "rlinspace (generic function with 1 method)" ] }, - "execution_count": 2, + "execution_count": 4, "metadata": {}, "output_type": "execute_result" } @@ -90,7 +81,31 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1-element Array{Array{T,1},1}:\n", + " [0,1]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "f = Field(0.0, Vector[])\n", + "push!(f.values, [0, 1])" + ] + }, + { + "cell_type": "code", + "execution_count": 9, "metadata": { "collapsed": false }, @@ -99,35 +114,29 @@ "name": "stdout", "output_type": "stream", "text": [ - "Converged. p = [4.2905787641613236,5.624156522861784,0.04706336239279979,6.561746550428803]\n" + "Converged. p = [4.2905787641613236,5.624156522861784,0.04706336239279979,6.561746550428803]" ] }, { - "data": { - "image/png": [ - "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVVX+//HX4X5HRU2JME0wM40k1LCLZt/MBLwEKDOa6ZTNTFdsUrMpv2XONFaj4zg2U7/SQY3xK2kKWpqW5TimJpaVElIhGqCCoSgXuZzfH2egkNsBgQ3nvJ+Px3moa++192fXFnmz1l7bZDabzYiIiIiIiNgAB6MLEBERERERaSkKOCIiIiIiYjMUcERERERExGYo4IiIiIiIiM1QwBEREREREZuhgCMiIiIiIjZDAUdERERERGyGAo6IiIiIiNgMBRwREREREbEZCjgiIiIiImIznJqy89dff83//u//kpqaSm5uLm5ubvTr14+HH36YX/7ylw32XblyJTNmzKhzW25uLt27d6/VnpeXx9atW7n66qtxd3dvSqkiIiIiImJDiouLyczMZPTo0XTt2rXe/ZoUcLKysjh//jz3338//v7+FBUVkZSUxNSpU8nMzOSZZ55p9BgLFiygd+/eNdp8fX3r3Hfr1q1MmTKlKSWKiIiIiIgNW716dYODKyaz2Wy+nBNUVlYSGhrKmTNnOHbsWL37VY3gfPbZZwwePNiqY+/evZtbbrmF1atX079//8sps0OKj49n8eLFRpchBtN9ILoHBHQfiIXuAwH7vQ+OHDnClClT+Pe//83w4cPr3a9JIzh1cXBwICAggMLCQqv2N5vNFBYW4uHhgaOjY4P7Vk1L69+/v9WhyJb4+vra5XVLTboPRPeAgO4DsdB9IKD7oLFHV5q1yEBRURF5eXl8++23LF68mK1btzJ79myr+o4cORJfX188PT0ZN24cGRkZzSlBRERExC4cOXKEoKAgXF1dSUlJwdXVlaCgII4cOWJ0aSLtUrNGcGbNmsXrr79uOYCTE0uXLmXmzJkN9vH09GT69OmMHDkSHx8fPvvsM/785z8THh5OamoqAQEBzSlFRERExGZ9/fXXhISEUF5eXt128eJFMjIyGDRoEIcOHbLLafwiDWlWwImPjyc2Npbs7GzWrFnDI488gru7O9OmTau3T0xMDDExMdV/joqKYvTo0dx2220sXLiQ1157rTmliIiIiNis8ePH1wg3P1deXk5UVBRHjx5t46pE2rdmBZx+/frRr18/AKZMmcLo0aN54okniI2NbdJyzsOHD2fo0KFs3769OWXYvLi4OKNLkHZA94HoHhDQfWCvsrKyLmu72CZ9PWjYZa+iBvD666/z61//mtTUVEJCQprUNzY2lg8//JC8vLxa21JTUwkNDeW2226rtZR0XFyc/ueKiIiITXNxcaGsrKze7c7Ozly8eLENKxJpG4mJiSQmJtZoO3v2LJ988gkHDhxocJGFy15FDSwv3QHLimpN9d1339GtW7cG91m8eLFdrxQhIiIi9slkMl3WdpGOqq7BjKrBj8Y0KZGcPn26VltZWRkJCQn4+fkxYMAAAHJyckhLS6sxZ7Suvlu2bCE1NZW77767KWWIiIiI2LySEnB0DGxwn8DAhreL2KMmjeDMnDmTwsJCbrvtNvz9/cnNzWXNmjWkp6ezYsWK6vfaPP300yQkJJCZmVn9Fy88PJzBgwcTGhqKr68vqampvPXWWwQGBjJv3ryWvzIRERGRDqq4GMaNg8rKTTg6DqKiou6FBh544IE2rkyk/WtSwJk8eTJvvvkmr732Gvn5+fj4+DB06FCWLVvGqFGjqvczmUy1hkwnT57M5s2b2bZtG0VFRfj7+/PQQw8xf/78RqeoiYiIiNiLoiKIioI9e+C99/rTo8choqKiyMrKwmw2YzKZuOqqq7jhhhuYN28egYGBei5Z5GdaZJGB1lI1z66xB4lEREREbMGFCxAZCfv2webNcPvt9e9bWVnJjBkzWLVqFYmJicTGxrZdoSIGsDYbtMgiAyIiIiJyec6fh7FjITUV3n8fbrml4f0dHBx48803qaio4Be/+AWOjo7ce++9bVOsSDvW9GXPRERERKRFFRbC3XfDwYOwdWvj4aaKo6MjK1euJDY2lsmTJ7Nhw4bWLVSkA1DAERERETHQ2bMwejR8+SVs2wbh4U3r7+joSEJCAhMnTiQ2NpZNmza1TqEiHYQCjoiIiIhBCgos4ebIEdi+HYYNa95xnJycWL16NePGjSM6OprNmze3bKEiHYgCjoiIiIgBfvwR/ud/ID3dEm7Cwi7veM7OziQmJjJ27FgmTpzI+++/3zKFinQwCjgiIiIibSw/H0aNgu+/hw8/BCtezm4VZ2dn1q5dy+jRoxk/fjzbtm1rmQOLdCAKOCIiIiJtKC/PEm6OH7eEm5CQlj2+i4sL69at484772TcuHHs2LGjZU8g0s4p4IiIiIi0kdOn4Y47IDsbPvoIBg1qnfO4urqSlJTEiBEjiIyMZN++fa1zIpF2SAFHREREpA2cPAkjR8KpU7BzJ1x/feuez83NjQ0bNhAfH8+AAQNa92Qi7Yhe9CkiIiLSynJzLSM3BQWWcHPttW1zXjc3NxYuXNg2JxNpJxRwRERERFpRdrYl3BQWWsJNcLDRFYnYNgUcERERkVbyww+WaWnFxfDxx9C3r9EVidg+BRwRERGRVnD8uCXclJVZwk2fPkZXJGIfFHBEREREWtixY5ZwYzZbpqX17m10RSL2Q6uoiYiIiLSg77+H228Hk0nhRsQICjgiIiIiLeS772DECHBysoSbXr2MrkjE/ijgiIiIiLSAjAzLyI2rq+WZm6uuMroiEfukgCMiIiJymdLTLeHG09MSbq680uiKROyX1QHn66+/JiYmhmuuuQZPT0/8/PwIDw9nzZo1VvUvKChg5syZdOvWDS8vL+644w4OHjzY7MJFRERE2oO0NEu48fW1TEvr2dPoikTsm9WrqGVlZXH+/Hnuv/9+/P39KSoqIikpialTp5KZmckzzzxTb9/KykrGjh3LoUOHmD17Nn5+fixfvpwRI0Zw4MAB+mpReBEREemADh+2vMSza1fYsQOuuMLoikTEZDabzc3tXFlZSWhoKGfOnOHYsWP17vd///d/TJ48maSkJCZOnAhAXl4ewcHBjBkzpt5RoNTUVEJDQzlw4ACDBw9ubpkiIiIiLe6rr2DUKEuo2bEDunUzuiIR22ZtNrisZ3AcHBwICAjA2dm5wf2SkpLo0aNHdbgB6Nq1K7GxsWzcuJGysrLLKUNERESkTR06ZHnPjb8/fPihwo1Ie9LkgFNUVEReXh7ffvstixcvZuvWrcyePbvBPgcPHqwzZYWFhVFUVER6enpTyxARERExxOefW8JNYKBl5KZrV6MrEpGfa3LAmTVrFt27dycoKIg5c+awdOlSZs6c2WCfnJwcetbxxF1VW3Z2dlPLEBEREWlzqamWZ2769IHt26FLF6MrEpFLWb3IQJX4+HhiY2PJzs5mzZo1PPLII7i7uzNt2rR6+5SUlODq6lqr3c3NDYDi4uKmliEiIiLSpvbvh7vuguBg2LoVOnUyuiIRqUuTA06/fv3o168fAFOmTGH06NE88cQTxMbG4u7uXmcfd3d3SktLa7WXlJRUbxcRERFpr/buhdGj4brr4L33LEtCi0j71OSAc6l7772XDz74gG+++YaQkJA69+nZs2ed09BycnIA8Pf3b/Ac8fHx+F7ylSQuLo64uLhmVi0iIiJinT17LOFm0CBLuPH2NroiEduXmJhIYmJijbazZ89a1feyA07V9DIHh/of5wkJCWHXrl2YzWZMJlN1+969e/H09CQ4OLjBcyxevFjLRIuIiEib270b7r4bBg+GzZvBy8voikTsQ12DGVXLRDfG6kUGTp8+XautrKyMhIQE/Pz8GDBgAGAZlUlLS6O8vLx6v+joaE6ePMn69eur2/Ly8li3bh2RkZGNLjMtIiIi0tY++cQychMWBlu2KNyIdBRWj+DMnDmTwsJCbrvtNvz9/cnNzWXNmjWkp6ezYsUKHB0dAXj66adJSEggMzOTwMBAwBJwhg0bxvTp0zl8+DB+fn4sX74cs9nM888/3zpXJiIiItJMH30EERFw882waRN4eBhdkYhYy+qAM3nyZN58801ee+018vPz8fHxYejQoSxbtoxRo0ZV72cymWpMQwPL9LUtW7bw1FNPsXTpUoqLixkyZAgJCQkEBQW13NWIiIiIXKYdOyAyEm65BTZuBK2FJNKxmMxms9noIupTNc/uwIEDegZHREREWt22bTBuHIwYARs2wH/faCEi7YC12aDJL/oUERERsUXvvQdRUTBqFLz7rsKNSEelgCMiIiJ2LyUFxo+3vMjznXegjveTi0gHoYAjIiIidm3jRpg4EcaOhaQkhRuRjk4BR0REROzWhg0QHW2ZmrZ2Lbi4GF2RiFwuBRwRERGxS0lJEBtrGb1JTAS9lk/ENijgiIiIiN1ZuxYmT4aYGFizRuFGxJYo4IiIiIhdeftt+MUvLJ9Vq8DJ6rcCikhHoIAjIiIidiMhAaZOtXxWrABHR6MrEpGWpoAjIiIidmHlSrj/fsvnrbcUbkRslQKOiIiI2Lw334QZM+DBB+GNN8BB3wGJ2Cz99RYRERGb9vrr8MAD8Otfw2uvKdyI2Dr9FRcRERGb9dpr8NBD8Oij8Le/KdyI2AP9NRcRERGbtGwZ/Pa38Pjj8Je/gMlkdEUi0hYUcERERMTm/OUvllGbJ5+ExYsVbkTsiQKOiIiI2JRXX4UnnoA5c+DllxVuROyNAo6IiIjYjEWL4He/g2eegT/+UeFGxB4p4IiIiIhN+MMfLKM2zz0HCxYo3IjYKwUcERER6fAWLLCM2jz/vOWjcCNivxRwREREpENISUnB2dkZk8lU/XF2dmby5BSeew5efNEyeiMi9s3qgLN//34eeeQRBgwYgJeXF7169WLSpEkcPXq00b4rV67EwcGhzs+pU6cu6wJERETE9m3cuJHIyEjKy8trtJeXl7N2bST33ZfCM88YVJyItCtO1u74pz/9iT179hATE8OgQYPIyclh2bJlDB48mE8//ZQBAwY0eowFCxbQu3fvGm2+vr5Nr1pERETsSnR0dIPb3357Av/8Z1kbVSMi7ZnVAefJJ58kLCwMJ6efukyaNImBAwfy0ksvsWrVqkaPMWbMGAYPHty8SkVERMRuXTpy09TtImI/rJ6idvPNN9cINwB9+/bluuuuIy0tzapjmM1mCgsLqaioaFqVIiIiIiIiVrisRQbMZjMnT56ka9euVu0/cuRIfH198fT0ZNy4cWRkZFzO6UVERERERGqweopaXdasWUN2djYvvvhig/t5enoyffp0Ro4ciY+PD5999hl//vOfCQ8PJzU1lYCAgMspQ0RERGzU2ZKz/PHff4TOwI/173fpLBMRsV/N/mqQlpbGww8/THh4ONOmTWtw35iYGGJiYqr/HBUVxejRo7nttttYuHAhr732WnPLEBERERtUXlnOGwfe4Lmdz1FUVkTcc3EkxifWu/+GDRvasDoRac+aNUUtNzeXsWPH0rlzZ5KSkjA1421aw4cPZ+jQoWzfvr05JYiIiIgNMpvNbDm6hUGvDeLhLQ8zNmgs6Y+k8/YTb5OcnFxrpKbqe5AdO3YYUa6ItENNHsE5e/YsY8aM4dy5c+zatYsePXo0++QBAQGkp6c3ul98fHyt5aTj4uKIi4tr9rlFRESkfTl08hC/2/Y7PvjuA0ZcPYLVE1czuOdPq69GRERQVlZzKejKykquvvpqlixZQkREBKNGjWrrskWkFSQmJpKYWHPU9uzZs1b1NZnNZrO1JyopKeGuu+7i4MGDbN++naFDhzat0kvcdNNNXLhwgSNHjtS5PTU1ldDQUA4cOKDlpUVERGxU7vlcnv3wWd48+CZ9u/TllbteITI40uoZIidOnKBPnz44OTmRnZ1Np06dWrliETGCtdnA6ilqFRUVTJo0ib1797Ju3bp6w01ubi5paWk11qM/ffp0rf22bNlCamoqd999t7UliIiIiA0pKivixU9epO/Svrxz5B2W3L2Er377FVH9opo0/T0gIIBVq1ZRXFzM8OHDW7FiEekImvSiz+TkZCIjI8nLy2P16tU1tk+ZMgWAuXPnkpCQQGZmJoGBgQCEh4czePBgQkND8fX1JTU1lbfeeovAwEDmzZvXgpcjIiIi7V2luZK3v3ybeTvmkXs+l0eHPMrvb/s9nd07N/uYkyZNYtOmTbz99ts89thjLF26tAUrFpGOxOqA88UXX2AymUhOTiY5ObnGNpPJVB1wTCZTrZ+6TJ48mc2bN7Nt2zaKiorw9/fnoYceYv78+XTr1q0FLkNEREQ6gl3HdvHktifZn72fif0n8qc7/0TfLn1b5NirVq1i9+7d/PWvf2Xs2LGMHj26RY4rIh1Lk57BaWt6BkdERMQ2fHvmW+Zsn8M7R94htGcoi0cv5tZet7b4ebKzs7n66qtxdHTkhx9+oEuXLi1+DhExRos/gyMiIiLSVD8W/8iTW5+k/9/6s/eHvSSMT2Dfg/taJdwA+Pv78/bbb1NSUqLncUTslAKOiIiItLiyijKW7l1K37/25R8H/sFztz/HN498w9QbpuJgat1vP6Kjo5k6dWr1S8lFxL4o4IiIiEiLMZvNbPpmE9e/dj3xW+OZeO1EMh7L4Pe3/R4PZ482q2PlypVcffXVLF++nPfee6/NzisixlPAERERkRbxee7njEoYxbh/jeMqn6tInZnKG1Fv0MOr+S8Fby4HBwf27NmDi4sLEydOJC8vr81rEBFjKOCIiIjIZckuzGbGxhkM/sdgcs7nkBKXwgdTP+CGHjcYWlePHj30PI6IHbJ6mWgRERGRn7tw8QKv/OcVFv1nER7OHiy7ZxkPDn4QZ0dno0urdu+993L//ffz8ccfU1JSgpubm9EliUgrU8ARERGRJqk0V5LwRQLPfPgMeUV5PD70cebdOo9Obp2MLq1Ob775JmCZtiYitk8BR0RERKy2M3Mns7bO4mDuQWIHxPLSqJfo3bm30WU1SMFGxL4o4IiIiEij0vPTmf3BbDZ+s5GhVw5l94zdhF8VbnRZIiK1KOCIiIhIvfKL8nnh4xdY/tly/L39Sbw3kUkDJmEymYwuTUSkTgo4IiIiUsvFiov8bd/feOGTF6iorGDByAU8PvRx3J3djS5NRKRBCjgiIiJSzWw2827au8zePpvvfvyOB258gBdGvsAVXlcYXZqIiFUUcERERASAA9kHmLVtFp8c+4TR14xmw6QNXN/9eqPLEhFpEgUcERERO3f87HGe+fAZVh1axYBuA3j/l+8zuu9oo8sSEWkWBRwRERE7df7ief707z/xyp5X8HH14R8R/2DGjTNwctC3ByLScekrmIiIiJ2pqKxg5ecr+f1Hv+fH4h958uYnmXPLHHxcfYwuTUTksingiIiI2JHt323nyW1PcujkIX4x8Bf84Y4/0KtTL6PLEhFpMQo4IiIiduDI6SM89cFTbD66mfCrwvn0V58yNGCo0WWJiLQ4B2t33L9/P4888ggDBgzAy8uLXr16MWnSJI4ePWpV/4KCAmbOnEm3bt3w8vLijjvu4ODBg80uXERERBp3+sJpHtnyCANfG8jh04dZF7OOf0//t8KNiNgsq0dw/vSnP7Fnzx5iYmIYNGgQOTk5LFu2jMGDB/Ppp58yYMCAevtWVlYyduxYDh06xOzZs/Hz82P58uWMGDGCAwcO0Ldv3xa5GBEREbEoLS9l6d6lLNy1EDNm/jjqjzw69FHcnNyMLk1EpFVZHXCefPJJwsLCcHL6qcukSZMYOHAgL730EqtWraq3b1JSEnv27CEpKYmJEycCEBsbS3BwMPPnz2fNmjWXcQkiIiJSxWw2k3Q4iTnb55B1Notf3/Rr5t8+n26e3YwuTUSkTVgdcG6++eZabX379uW6664jLS2twb5JSUn06NGjOtwAdO3aldjYWFavXk1ZWRnOzs5NKFtEREQutffEXmZtm8V/jv+HsUFj2fLLLVzb9VqjyxIRaVNWP4NTF7PZzMmTJ+natWuD+x08eJDBgwfXag8LC6OoqIj09PTLKUPE5hw5coSgoCBcXV1xcXHB1dWVoKAgjhw5YnRpItIOHSs4xi/e+QXD3hzG+Yvn+WDqB6T8IkXhRkTs0mUFnDVr1pCdnc2kSZMa3C8nJ4eePXvWaq9qy87OvpwyRGzK119/zaBBg8jIyODixYuUlZVx8eJFMjIyGDRokEKOiFQ7V3qOp7c/Tb9l/fgo8yPejHqT1Jmp3NnnTqNLExExTLOXiU5LS+Phhx8mPDycadOmNbhvSUkJrq6utdrd3CwPOhYXFze3DBGbM378eMrLy+vcVl5eTlRUlNWrF4qIbSqvLOetg2/x7EfPUlhayFPhTzF7+Gy8Xb2NLk1ExHDNCji5ubmMHTuWzp07k5SUhMlkanB/d3d3SktLa7WXlJRUbxcRi6ysrMvaLiK2bWvGVn73we/46tRXTB00lT+M+gMBPgFGlyUi0m40OeCcPXuWMWPGcO7cOXbt2kWPHj0a7dOzZ886p6Hl5OQA4O/v32D/+Ph4fH19a7TFxcURFxfXhMpFOgaz2XxZ20XENn116it+t+13bP12K7cG3sr+B/dzk/9NRpclItIqEhMTSUxMrNF29uxZq/o2KeCUlJQQGRlJRkYG27dv59prrXt4MSQkhF27dmE2m2uM9uzduxdPT0+Cg4Mb7L948eI6FykQsUWNjYiay8yk3pKKX4QffhF+eA7wbLSPiHRcpy6c4rmPnuON1Dfo07kP62PXM/7a8fp7LyI2ra7BjNTUVEJDQxvta/UiAxUVFUyaNIm9e/eybt06hg6t+w3Iubm5pKWl1XiGIDo6mpMnT7J+/frqtry8PNatW0dkZKSWiBb5mcDAwAa3B3QLwKWbC8cWHOOzgZ/xae9PSX8knfz386koqWijKkWkpaSkpODs7IzJZKr+ODs7s37jev6464/0XdqXtV+v5ZX/eYWvf/s1E/pPULgREWlAk170mZycTGRkJHl5eaxevbrG9ilTpgAwd+5cEhISyMzMrP5GLTo6mmHDhjF9+nQOHz6Mn58fy5cvx2w28/zzz7fg5Yh0fJs2bWLgwIFUVNQOK05OTmz5eAv9+/enoqSCgp0F5Kfkc2bzGbL/lo2DhwOd7+xsGd0Z64erf+3FPUSk/di4cSPjx4+v1V5eXs694++FX8LjUx7n2duexc/Dz4AKRUQ6HqsDzhdffIHJZCI5OZnk5OQa20wmU3XAqfrp0885ODiwZcsWnnrqKZYuXUpxcTFDhgwhISGBoKCgFrgMEdvRv39/Xn75ZWbNmlU9umkymQgMDGTTpk30798fAEc3R/zu9sPvbj/MfzVTdLiI/M355Cfnk/7rdKgEr1Cv6qls3oO9MTnop74i7Ul0dHSD2x3XOrJk9ZI2qkZExDaYzO34ieWqeXYHDhzQMzhiV6ZNm0Zqaipffvlls/qXnSnjzHtnLKM775+hvKAc5yuc8RtrCTud7+yMk3ezV4kXkRZizVSzdvzPtIhIm7I2G+g7HJF2pqKigi1btvDAAw80+xjOXZy54pdXcMUvr6CyrJJz/zlHfko++Sn55L6Vi8nFRKcRnapHd9x7a6l2kTZx5gx88gns3Ak7d+IMlBldk4iIjVHAEWln9u/fT15eHhERES1yPAdnBzrd3olOt3fimpevofjbYkvY2ZzPt09+S8ZjGXhc51Eddnxu9sHByer1R0SkIZcEGg4dArMZeveG22+n4osvjK5QRMTmKOCItDMpKSl06dKFYcOGtcrx3a9xJ+DxAAIeD6C8sJwft/1I/uZ8clfmcnzRcZw6O9FlTBf8IvzoMroLzl20yqGI1eoLNFdfDSNGQHw83H675c+Aw+rVVP5s1dFLOTnpn2kRkabSV06RdiYlJYV77rkHR0fHVj+Xk7cT3e7tRrd7u2GuNFP4WSH5yZapbKfePgWO4Bvui1+kZXTH41oPLU8r8nONBZonnoCRI6FXrzq7b9iwgcjIyHoPv2HDhtaoWkTEpingiLQjx48f54svvuDpp59u83ObHEz4DPHBZ4gPvRf0pvSHUsuqbJvzyZyfyXezv8Otj1v1VLZOt3XCwVVT2cTONDLljFmzLL/WE2guFRERQXJyMhMmTKjx/jiA9evXt9hUVRERe6KAI9KObN68GUdHR0aPHm10Kbhe6Yr/TH/8Z/pTUVxBwYcF5G/OJ29DHj8s/QFHL0c6/4/lnTtd7umCaw+9c0dsUBOnnDVHREQEZWU/LTXwz3/+k/vvv59Dhw4xYcKEy7wAERH7o2WiRdqRiIgILly4wEcffWR0KfUym81c+PJC9aps5z49B2bwDvO2jO5E+uEV4qWpbNIxNRZobr+9wSlnLaGyshIvLy86depEdnZ2q51HRKSj0TLRIh1MUVERO3bs4MUXXzS6lAaZTCa8BnnhNciLXvN6cfH0xep37hx/9TiZ8zNx8Xf56Z07ozrj6Nn6zxOJNEsLTzlrCQ4ODkRFRbF27VpSU1P1Az4RkSZSwBFpJz788ENKSko63Jx7l24u9LivBz3u60HlxUrO/vts9ehOzhs5mFxNdL7DMpXNb6wfbr3cjC5Z7FkbTDlrCYsWLWLt2rXMmTOHDz74wNBaREQ6GgUckXYiJSWFvn37EhwcbHQpzebg4kDnOzrT+Y7O9P1zX4rSi6rfuZPxeAZHHz6K50DPn965M9QHk6OmskkrusxVzowSGBhIUFAQO3fupLy8XMtFi4g0gb5iirQDZrOZlJQUYmJibOrZFY9gDzxmeXDVrKsoP1vOmW1nqkd2sv6YhZOfE373/Hcq212dce6kd+7IZWqHU86aa/bs2Tz44IO8/PLLhqysKCLSUWmRAZF24PPPP+fGG29k+/btjBo1yuhyWp25wsy5feeq37lz4csLmJxM+N7iWz2649HPw+gypSNoaIRm5EhLmGkHU86ao7KyEg8PD7p27cqJEyeMLkdExHBaZECkA0lJScHb25tbb73V6FLahMlW7mdcAAAgAElEQVTRhO/Nvvje7EufP/ShJKuk+p073//+e7793be4B7lXhx3fW3xxcNE7d4QOO+WsORwcHIiIiOCdd97hq6++4vrrrze6JBGRDkEBR6QduO6665g7dy4uLi5Gl2IIt0A3rvzNlVz5myupKKrgx+0/kr85n1NrT3Fi8QkcvR3pMrqL5Z07Y7rg0t0+/zvZJRuactYcL7/8Mu+88w5PPfUU7733ntHliIh0CAo4Iu3AxIkTjS6h3XD0cKRrVFe6RnXFbDZz/vPz1auypU1PA8BnqE/1O3c8B3ra1HNLdq++QNOrl2Vk5oknLCM1HXDKWXP07t2bPn36sH37di02ICJiJc35EJF2y2Qy4X2jN1c/ezWhe0MJzwmn35v9cLnShayXsvjshs/4tNenpP8mnfzN+VQUV9Ton5KSgrOzMyaTqfrj7OxMSkqKQVcktZw5A+++awkuISHQtStMmAAbN8KNN8Jbb8H330NmJqxYAfffbzfhpsqTTz5JeXk5ixcvNroUEZEOQYsMiEiHVFlaScEnBZaFCjbnU/JdCQ7uDnQeZXnnzqdOn3LvA/fW2z85ObnDvXPIJjQ25axqYQAbnXLWHJWVlbi7u3PFFVeQlZVldDkiIobRIgMiYtMcXB3o8j9d6PI/Xej7l74UfVNUHXbSH07n0YpHG+w/YcIEysrK2qha25OSksKECRMoLy+vbnNycmLDhg01g2Njq5y1kxdrtmcODg7cc889vPvuuxw5coT+/fsbXZKISLumgCMiHZ7JZMLzWk88r/Uk8KlAyn4sI6dLToN9ysvLGbR/P84mEy4ODjibTNW/d7n09z9rq/p9VR+XS/araqtvv5+fp77zOf53Ol17tXHjRsaPH1+rvby8nMjISJLnzSPiwgWbX+WsLS1atIhTp07h4KCZ5SIijWlSwLlw4QKLFi1i79697Nu3j4KCAlasWMG0adMa7bty5UpmzJhR57bc3Fy6d+/elFJEROrl3NkZM43Pvh3RqRNlZjMXKyurf71oNlNmNlNYUUFZeflP2y7Zr6qt7L99LlZWUtHoGZtwDZeEokuDUK1AVsf21gph0ZMmNVj7hD/8gTI7WeWsrQQFBbF7926jyxAR6RCaFHBOnz7NggUL6NWrFyEhIezcubPJP2VcsGABvXv3rtHm6+vbpGOIiLSEpUFBLXq8yv+Go7KfBaW6wlNdQamh8FR1zNI6jn3pfj8PZtaer8nBrKLhHuUA333XzP+KIiIil6dJAcff3796tOXAgQOEhYU1+YRjxozRggEi0uqcnJxqPB9S1/aW5mAy4Woy4drBphFVms2U/zcMlV4ayOoISrc38N9VRETEaE36F97FxaV6KllzF18zm80UFhbi4eGBo6Njs44hItKYDRs2EBkZ2eB2sXComp4GeOrrsoiIdHBt/mPGkSNH4uvri6enJ+PGjSMjI6OtSxAROxAREUFycnKtH6Q4OTlpiejL1Njol15GKSIiRmqzgOPp6cn06dNZvnw57777LrNnz2bHjh2Eh4dz4sSJtipDROxIREQEixYtAmDLli2YzWbKysoUbi5TY6NfGh0TEREjtVnAiYmJ4c0332TKlClERUXxwgsvsHXrVvLz81m4cGFblSEidmb//v0A3H777QZXYjuqRscuHanR6JiIiLQHhs4jGD58OEOHDmX79u1GliEiNuzIkSO4uLjg4eFhdCk2JSIigrKyMkaPHs0HH3xAZWWl0SWJiIgA7eBFnwEBAaSnpze4T3x8fK2lpOPi4oiLi2vN0kTEBpw4cQI/Pz+jy7BZffr0wWw2c+rUKb3PTEREWkxiYiKJiYk12s6ePWtVX8MDznfffUe3bt0a3Gfx4sVaWlpEmqWgoIAhQ4YYXYbNuvbaawE4cOAAY8aMMbgaERGxFXUNZqSmphIaGtpo31Z5Bic3N5e0tLQa76A4ffp0rf22bNlCamoqd999d2uUISJ27tSpU1RUVHD99dcbXYrNGjhwIABffvmlwZWIiIhYNHkEZ9myZRQUFJCdnQ3Apk2byMrKAuCxxx7Dx8eHuXPnkpCQQGZmJoGBgQCEh4czePBgQkND8fX1JTU1lbfeeovAwEDmzZvXgpckImKxY8cOAG6++WaDK7FdN910EwBpaWkGVyIiImLR5IDz6quvcuzYMQBMJhMbNmxg/fr1mEwm7rvvPnx8fDCZTJhMphr9Jk+ezObNm9m2bRtFRUX4+/vz0EMPMX/+/EanqImINMd//vMfAEaNGmVwJbar6mt+1b8LIiIiRmtywPn+++8b3WfFihWsWLGiRtuCBQtYsGBBU08nItJsX3zxBQ4ODtUjydI63Nzcqkf1RUREjNZm78EREWlr33//PT4+PkaXYfN8fX3Jy8szugwRERFAAUdEbFheXh5XXnml0WXYvG7dulFYWGh0GSIiIoACjojYqIsXL1JSUkK/fv2MLsXmBQQEUFpaanQZIiIigAKOiNioqgUGqlb5ktbTu3dvAD2HIyIi7YICjojYpI8//hiAESNGGFuIHejfvz9geQGbiIiI0RRwRMQm+fr6csUVVxAWFmZ0KTYvJCQE0Ms+RUSkfVDAERGb9MQTT5Cbm4uTU5NXw5cmqgo433zzjcGViIiIKOCIiMhl8vLyonfv3nrfkIiItAv60aaIiFy27777zugSREREAI3giIiInUpJScHZ2RmTyVT9cXZ2JiUlxejSRETkMijgiIiI3dm4cSORkZGUl5fXaC8vLycyMlIhR0SkA1PAERERuxMdHd3g9gkTJrRRJS1Lo1IiInoGR0RE7EBxMRw7BpmZlk95ubnB/cvLy5k4cSKenp54enri4eFR41drf+/o6Ngm1weWUanx48fXeS2RkZEkJycTERHRZvWIiBhFAUdERDq8SwPMpZ+TJ3/a15I5TFYcs5j8/HwuXLhQ/SkqKuLChQuUlJRYVZeLi0uDIak5oennx/n5MujWjEqVlZVZVbeISEemgCMiIu1eUwPMVVfB1VdD//4wZozl91WfK68EZ+fyS09Ry3vvvVfvtoqKCoqKiqoDz8/DT9Xv69r287YLFy60SICqCjyXPk90qca2i4jYCgUcERFpUEpKChMmTKjxDbKTkxMbNmxosSlPzQ0w114Ld99t+X3v3j8FmMbe7+rk5NTgN/yNvSDW0dERb29vvL29rbi6pquoqKC4uLhWIGooUL3wwguNHjctbTre3kPw8RmCp+dAHBxcWqV+EREjKeCIiEi9Wuq5jssdgenVyxJgevcGf//GA0xjNmzYQGRkZIPbjeTo6IiXlxdeXl5W97Em4Jw/f4iTJ1djNpdjMrni7X1jdeDx9h6Cu3tfTKbGp++JiLRnCjgiIlIva5/raG6Aue66n6aQ9epl+TUg4PIDTGMiIiJITk5u9ZGptmTNqNRNNx2goqKY8+c/p7BwH+fO7ePMmS388MPS/+7T6WeBJwxv7yG4uvZoq0sQEWkRTfon5MKFCyxatIi9e/eyb98+CgoKWLFiBdOmTbOqf0FBAbNnz2bDhg0UFxczZMgQXn31VW688cZmFS8iIq3Lmuc6evSoHWACAiyjLfU9A9PaAcYaERERNvXQvbWjUo6O7vj63oyv783V28rKzlBYuJ9z5/ZRWLiP7OzXKSt7EQBX16tqjPJ4e4fi5NQ6U/NERFpCk/6JOX36NAsWLKBXr16EhISwc+dOq4eyKysrGTt2LIcOHWL27Nn4+fmxfPlyRowYwYEDB+jbt2+zLkBERIz1m9/89AxMr17tJ8DYm8sZlXJ27kKXLqPp0mU0AGazmdLS49WB59y5fWRmvkBl5QXAhIdH//8GnqH4+ITpeR4RaVea9E+Qv78/ubm5dO/enQMHDhAWFmZ136SkJPbs2UNSUhITJ04EIDY2luDgYObPn8+aNWuaVrmIiLS4sooydh/fTUp6CsnpyeAIVDTcZ/78NilNrNBSo1Imkwk3t0Dc3ALp3t0yTdFsruDChcM1Rnr0PI+ItEdNCjguLi50794dsPx0pymSkpLo0aNHdbgB6Nq1K7GxsaxevZqysjKcnZ2bdEwREbl8Z4rP8H7G+ySnJ/N+xvsUlBTQw6sHEUERZDhnUFlRWW/fxlYbE9thMjni5TUQL6+B9Ow5A+C/z/McrA48NZ/n6Yy3d1h14PHxGYKLyxVGXoKI2Ik2+5fp4MGDDB48uFZ7WFgYr7/+Ounp6QwYMKCtyhERsVtms5m0vDRS0lNIOZrC7qzdVJgrGNxzMI8NeYyI4AhC/UNxMDkwbt24dr3amBjL8jxPOL6+4dVtLf08z5EjR4iKiiIrKwuz2YzJZCIwMJBNmzbRv3//Vr0+EemY2izg5OTkMGLEiFrtPXv2BCA7O1sBR0SklVysuMiuY7uqp559++O3uDu5M6rPKJaPXc7YoLFc6XNlrX62uNqYtK76n+fZW8/zPNfVGOWxPM9jmdHx9ddfExISUmuxi4yMDAYNGsShQ4cUckSkljYLOCUlJbi6utZqd3NzA6C4uLitShERsQt5RXm8d/Q9ktOT2frtVs6VnuNK7yuJCI5gSfAS7uh9Bx7OHo0ex9ZWG5O2VfN5nhjg58/z7OPcuf0UFu4lNzcBqMDBwQ0vrxC8vYcQEfGvelfyKy8vJyoqiqNHj7bh1YhIR9BmAcfd3Z3S0tJa7SUlJdXbRUSk+cxmM4dPHyY5PZmU9BT2nNhDpbmSMP8wnrz5SSKDIwnpEaIHv8VwNZ/n+RUAFRVF/32eZz+FhfvJz9/MiROnGjxOVlZWW5QrIh1MmwWcnj17kp2dXas9JycHsKzQVp/4+Hh8fX1rtMXFxREXF9eyRYqIdDCl5aV8fOxjkr9JJuVoCpkFmXg6e3Jnnzt5PeJ1xgaPpYeXXtQo7Z+jowe+vsPx9R1e3WYyuQD1jx42dcEjEek4EhMTSUxMrNF29uxZq/q2WcAJCQlh165d1Q8IVtm7dy+enp4EBwfX23fx4sV1LlAgImKPTl04xZajW0hOT2bbt9s4f/E8gb6BRARFEBEcwcjeI3FzcjO6TJHL1thoo0YjRWxXXYMZqamphIaGNtq3VQJObm4uBQUF9O3bt3oJ0ejoaJKSkli/fj333nsvAHl5eaxbt47IyEgtES0iUg+z2cyXp74k+ZtkktOT2ffDPgCGBgxl7vC5RPaLZGD3gfpmT2xOYGAgGRkZDW4XEblUkwPOsmXLKCgoqJ5utmnTpuo5sI899hg+Pj7MnTuXhIQEMjMzq7/4REdHM2zYMKZPn87hw4fx8/Nj+fLlmM1mnn/++Ra8JBGRjq+kvISPvv+o+nma4+eO4+XixehrRvObm37DmKAxdPfsbnSZIq1q06ZNDBo0qM6FBpycnNi0aZMBVYlIe9fkgPPqq69y7NgxwDI0vGHDBtavX4/JZOK+++7Dx8cHk8lU6yeJDg4ObNmyhaeeeoqlS5dSXFzMkCFDSEhIICgoqGWuRkSkA8spzKmeevbBdx9QVFZE7069GX/teCKDI7mt1224OtVejVLEVvXv359Dhw7pPTgi0iQmczt+Qq9qnt2BAwf0DI6IdCgpKSmNvjvGbDZzMPdg9btpPsv+DAeTAzcH3ExEcASRwZFc1+06TT0TERHB+mzQZosMiIjYi40bNzJ+/Pha7eXl5URGRvLsP57lZM+TpBxNIbswGx9XH+7uezePDXmMMUFj6OrR1YCqRUREbIMCjohIC4uOjm5w+4JfL6DvX/oSe10skf0iuTXwVpwdtdCKiIhIS1DAERFpYfW9eb2aGdIfSdfUMxERkVbgYHQBIiI2pbISa2KLwo2IiEjrUMAREWkJFRXwr3/BoEFospmIiIhxFHBERC5HeTmsXg3XXw9xcRAQgLOjY4Ndql6ALCIiIi1PAUdEpDnKymDFCujfH6ZOhaAg2LsX3n+ff737boNdN2zY0EZFioiI2B8FHBGRprh4Ed54A4KDYcYMGDgQDhyATZtgyBAAIiIiSE5OrjVS4+TkRHJycvV7cERERKTlaZ6EiIg1SkrgrbfgpZfgxAmIjraEmoED69w9IiKCsrKyNi5SRERENIIjItKQoiL4y1+gTx949FG49Vb46iv4v/+rN9yIiIiIcTSCIyJSl/Pn4e9/h5dfhvx8mDIF5s2zTE0TERGRdksBR0Tk5woL4W9/g1dfhYICmDbNEmz69DG6MhEREbGCAo6ICFjCzNKlsGQJXLhgWUBg7lzo1cvoykRERKQJFHBExL6dOWMJNX/5i2WFtAcfhNmzISDA6MpERESkGRRwRMQ+nT4Nf/4zLFsGFRXwm9/A734HPXsaXZmIiIhcBgUcEbEvJ0/CK6/A8uXg4AC//S08+SR07250ZSIiItICFHBExD788INlRbR//ANcXCA+3vLx8zO6MhEREWlBCjgiYtuOH4c//QneeAM8PGDOHHj8cejc2ejKREREpBUo4IiIbfr+e3jpJVixAnx84Lnn4JFHwNfX6MpERESkFTk0ZefS0lLmzJmDv78/Hh4eDBs2jO3btzfab+XKlTg4ONT5OXXqVLOLFxGpJSPDssRzUBBs2AAvvgiZmfDMMwo3IiIidqBJIzj3338/77zzDvHx8QQFBbFixQruuecePvroI4YPH95o/wULFtC7d+8abb76hkNEWkJaGixcCG+/bVkw4OWX4aGHLNPSRERExG5YHXD27dvH2rVreeWVV5g1axYAU6dO5frrr2f27Nns3r270WOMGTOGwYMHN79aEZFLffWVJdisXQv+/pZ32jzwALi7G12ZiIiIGMDqKWpJSUk4OTkxc+bM6jZXV1d+9atfsWfPHn744YdGj2E2myksLKSioqJ51YqIVPn8c4iOhoEDYc8ey7LP334Ljz6qcCMiImLHrA44Bw8eJDg4GC8vrxrtYWFhAHz++eeNHmPkyJH4+vri6enJuHHjyMjIaGK5ImIvjhw5QlBQEK6urri4uODq6kpQUBBH1q2DcePgxhvh4EHL6mjp6fDrX4Orq9Fli4iIiMGsnqKWk5NDzzre8F3Vlp2dXW9fT09Ppk+fzsiRI/Hx8eGzzz7jz3/+M+Hh4aSmphIQENCM0kXEVn399deEhIRQXl5eoz0jI4NBsbEc6tWL/itXwi9+Ac7OxhQpIiIi7ZLVAae4uBjXOn466ubmVr29PjExMcTExFT/OSoqitGjR3PbbbexcOFCXnvttabULCI2bvz48bXCTZVyIMrJiaPTprVtUSIiItIhWD1Fzd3dndLS0lrtJSUl1dubYvjw4QwdOtSqZaZFxL5kZWU1vP348TaqRERERDoaq0dwevbsWec0tJycHAD8/f2bfPKAgADS09Mb3S8+Pr7WctJxcXHExcU1+Zwi0v6ZzebL2i4iIiIdW2JiIomJiTXazp49a1VfqwPOjTfeyM6dOyksLMTb27u6fe/evQCEhIRYe6hq3333Hd26dWt0v8WLF2t5aRE7YjKZLmu7iIiIdGx1DWakpqYSGhraaF+rp6hFR0dTUVHB66+/Xt1WWlrKihUrGDZsGFdeeSUAubm5pKWl1Zg/f/r06VrH27JlC6mpqdx9993WliAidiIwMPCytouIiIj9snoEZ8iQIcTExPD0009z6tQprrnmGv75z3+SlZXFihUrqvebO3cuCQkJZGZmVn8TEh4ezuDBgwkNDcXX15fU1FTeeustAgMDmTdvXstflYh0aJs2bWLQoEF1LjTg5OTEpk2bDKhKREREOgKrAw5AQkICzz77LKtWreLHH3/khhtuICUlhVtuuaV6H5PJVGv6yOTJk9m8eTPbtm2jqKgIf39/HnroIebPn2/VFDURsS/9+/fn0KFDREVFkZWVhdlsxmQyERgYyKZNm+jfv7/RJYqIiEg7ZTK346d1q+bZHThwQM/giIiIiIjYMWuzgdXP4IiIiIiIiLR3CjgiIiIiImIzFHBERERERMRmKOCIiIiIiIjNUMARERERERGboYAjIiIiIiI2QwFHRERERERshgKOiIiIiIjYDAUcERERERGxGQo4IiIiIiJiMxRwRERERETEZijgiIiIiIiIzVDAERERERERm6GAIyIiIiIiNkMBR0REREREbIYCjoiIiIiI2AwFHBERERERsRkKOCIiIiIiYjOaFHBKS0uZM2cO/v7+eHh4MGzYMLZv325V34KCAmbOnEm3bt3w8vLijjvu4ODBg80qWkREREREpC5NCjj3338/ixcvZurUqSxduhRHR0fuuecedu/e3WC/yspKxo4dS2JiIo899hiLFi3i1KlTjBgxgoyMjMu6ABERERERkSpO1u64b98+1q5dyyuvvMKsWbMAmDp1Ktdffz2zZ89uMOQkJSWxZ88ekpKSmDhxIgCxsbEEBwczf/581qxZc5mXISIiIiIi0oQRnKSkJJycnJg5c2Z1m6urK7/61a/Ys2cPP/zwQ4N9e/ToUR1uALp27UpsbCwbN26krKysmeWLiIiIiIj8xOqAc/DgQYKDg/Hy8qrRHhYWBsDnn3/eYN/BgwfXag8LC6OoqIj09HRryxAREREREamX1QEnJyeHnj171mqvasvOzm6VviIiIiIiItayOuAUFxfj6upaq93Nza16e31KSkqa3VdERERERMRaVgccd3d3SktLa7WXlJRUb2+NviIiIiIiItayOuD07NmzzqlkOTk5APj7+7dKX4D4+HiioqJqfBITE60tvcOyh2uUxuk+EN0DAroPxEL3gYB93AeJiYm1vv+Pj4+3qq/Vy0TfeOON7Ny5k8LCQry9vavb9+7dC0BISEi9fUNCQti1axdmsxmTyVSjr6enJ8HBwQ2ee/HixXUuUmDrEhMTiYuLM7oMMZjuA9E9IKD7QCx0HwjYx30QFxdX6xpTU1MJDQ1ttK/VIzjR0dFUVFTw+uuvV7eVlpayYsUKhg0bxpVXXglAbm4uaWlplJeX1+h78uRJ1q9fX92Wl5fHunXriIyMxNnZ2doyRERERERE6mX1CM6QIUOIiYnh6aef5tSpU1xzzTX885//JCsrixUrVlTvN3fuXBISEsjMzCQwMBCwBJxhw4Yxffp0Dh8+jJ+fH8uXL8dsNvP888+3/FWJiIiIiIhdsjrgACQkJPDss8+yatUqfvzxR2644QZSUlK45ZZbqvcxmUw1pqEBODg4sGXLFp566imWLl1KcXExQ4YMISEhgaCgoJa5EhERERERsXtNCjiurq4sWrSIRYsW1bvPihUraozoVOnUqRNvvPEGb7zxhtXnq1o++siRI00p02acPXuW1NRUo8sQg+k+EN0DAroPxEL3gYD93gdVmaCxV8yYzGazuS0Kao41a9YwZcoUo8sQEREREZF2YvXq1fzyl7+sd3u7Djh5eXls3bqVq6++Wu/KERERERGxY8XFxWRmZjJ69Gi6du1a737tOuCIiIiIiIg0hdXLRIuIiIiIiLR3CjgiIiIiImIzFHBERERERMRmKOCIiIiIiIjNUMDpgB588EEcHByIjIw0uhRpQzt27GDGjBkEBwfj6enJNddcw4MPPkhubq7RpUkrKC0tZc6cOfj7++Ph4cGwYcPYvn270WVJG9q/fz+PPPIIAwYMwMvLi169ejFp0iSOHj1qdGlioIULF+Lg4MDAgQONLkXaWGpqKlFRUfj5+eHp6cnAgQP561//anRZ7ZJWUetgPvvsM8LDw3FycuLOO+9k06ZNRpckbeSmm26ioKCAmJgYgoKC+Pbbb1m2bBkeHh58/vnnXHHFFUaXKC0oLi6Od955h/j4eIKCglixYgX79+/no48+Yvjw4UaXJ20gOjqaPXv2EBMTw6BBg8jJyWHZsmWcP3+eTz/9lAEDBhhdorSxEydO0K9fPxwcHOjduzeHDh0yuiRpI9u2bSMyMpLQ0FAmTZqEl5cXGRkZmM1mXnrpJaPLa3cUcDoQs9nM8OHDGTBgANu3b2fgwIEKOHbk3//+N7fcckuNtl27dnH77bfzzDPPsGDBAoMqk5a2b98+hg0bxiuvvMKsWbMAy4jO9ddfT/fu3dm9e7fBFUpb2LNnD2FhYTg5OVW3ZWRkMHDgQKKjo1m1apWB1YkRJk+eTH5+PuXl5eTl5fHll18aXZK0gXPnzhEcHMwtt9xCUlKS0eV0CJqi1oGsWrWKw4cP8+KLL6Jcan8uDTcAt956K126dCEtLc2AiqS1JCUl4eTkxMyZM6vbXF1d+dWvfsWePXv44YcfDKxO2srNN99cI9wA9O3bl+uuu05/5+3QJ598wjvvvMOSJUswm82YTCajS5I28vbbb3Pq1CkWLlwIwIULF6isrDS4qvZNAaeDKCwsZM6cOcybN09TkaTa+fPnKSwsbPBtvtLxHDx4kODgYLy8vGq0h4WFAfD5558bUZa0A2azmZMnT+rvvJ2pqKjg0Ucf5cEHH9TURDu0fft2fHx8OH78OP369cPb2xtfX19++9vfUlpaanR57ZICTgfxwgsv4OnpSXx8vNGlSDuyZMkSysrKmDRpktGlSAvKycmhZ8+etdqr2rKzs9u6JGkn1qxZQ3Z2tv7O25m///3vZGVlaSqynTp69Cjl5eWMHz+eMWPGsH79embMmMHf//53pk+fbnR57ZJT47tISzKbzVanbTc3NwDS09NZunQp//rXv3B2dm7N8qSNNOc+uNQnn3zC888/z6RJkxgxYkQLVidGKy4uxtXVtVZ71b1QXFzc1iVJO5CWlsbDDz9MeHg406ZNM7ocaSP5+fk899xzPPfcc/j5+Rldjhjg/PnzFBUV8Zvf/IYlS5YAMH78eC5evMg//vEPXnjhBfr27Wtwle2LRnDa2Mcff4yHh4dVn/T0dAAef/xxhg8fzoQJEwyuXlpKc+6Dn0tLS2PChAkMGjSI//f//p8BVyCtyd3dvc4AXFJSUr1d7Etubi5jx46lc+fOJCUl6fkLO/L73/+erl278uijjxpdihik6mt+XFxcjfaqP3/66Rf7ko8AAAOmSURBVKdtXlN7pxGcNta/f39Wrlxp1b49evTgww8/ZOvWraxfv57MzMzqbeXl5RQVFXHs2DG6dOmCt7d36xQsraKp98HPHT9+nLvuuovOnTuzZcsWPD09W6FCMVLPnj3rnIaWk5MDgL+/f1uXJAY6e/YsY8aM4dy5c+zatavW1wSxXUePHuWNN95gyZIlnDhxorq9pKSEixcvcuzYMXx8fOjcubOBVUpr8/f35/Dhw7Wewe7evTsAP/74oxFltWsKOG3siiuu4L777rN6/6ysLAAmTpxYa1t2dja9e/dmyZIlPPbYYy1Wo7S+pt4HVfLz87nrrrsoKyvjo48+0oITNurGG29k586d/7+9u2dpJQijOH78AJbpdw2Cha/BWpKAYBMUQYVgFYKNIAl2USsljSDoB1CElIpoZZGgiAQCAcVOEMGX2C0WUVBIzK304jW3dCeO/18525xmWc7s7jyqVqufNi9KpZIkqb+/31Q0+Ozl5UWxWExXV1fK5/Pq6uoyHQk+qlQqent709zcXNPnvOu6SqVSWltbM5AOfhkcHFQ+n9f9/b06Ozs/1t83wgKBgKloLYs5OC3u7u5OZ2dnn9YajYZmZmbkOI4WFhbU3d2tjo4OQwnhl+fnZ0WjUV1eXuro6EgDAwOmI+GbvM/BWV1d1fz8vKS/c3ACgYCKxaLhhPBDvV7X+Pi4Dg8Ptb+/r5GREdOR4DPP83R6evrpk8RGo6HFxUU9PT1pfX1dwWCQk9Usd35+rlAopHg8rlwu97Eej8e1u7urm5sb3uz+g4LzQzmOo97eXgZ9/iJjY2M6ODhQIpH4cqhAe3u7RkdHzQTDt5iamtLe3p7S6bSCwaC2t7dVLpdVKBSazkSCfVKplDY2NhSLxTQxMfHl+vT0tIFUaAXhcFie5zHo8xdJJpPa3NzU5OSkhoaGdHx8rJ2dHWUyGa2srJiO13IoOD+U67rq6emh4Pwiruvq9va26ZBXx3F0fX1tIBW+y+vrq5aWlpTL5fT4+Ki+vj4tLy9reHjYdDT4JBKJ6OTkpOk939bWpnq9biAVWkEkEpHnebq4uDAdBT6p1WrKZrPa2trSw8ODHMfR7Owsvyj8BwUHAAAAgDU4JhoAAACANSg4AAAAAKxBwQEAAABgDQoOAAAAAGtQcAAAAABYg4IDAAAAwBoUHAAAAADWoOAAAAAAsAYFBwAAAIA1KDgAAAAArEHBAQAAAGANCg4AAAAAa1BwAAAAAFjjD8pR/tt+TMPDAAAAAElFTkSuQmCC" - ], - "text/plain": [ - "PyPlot.Figure(PyObject )" - ] - }, - "metadata": {}, - "output_type": "display_data" - }, - { - "data": { - "text/plain": [ - "(-0.1,3.6)" - ] - }, - "execution_count": 3, - "metadata": {}, - "output_type": "execute_result" + "ename": "LoadError", + "evalue": "LoadError: \"Empty set of fields.\"\nwhile loading In[9], in expression starting on line 60", + "output_type": "error", + "traceback": [ + "LoadError: \"Empty set of fields.\"\nwhile loading In[9], in expression starting on line 60", + "", + " in interpolate at C:\\Users\\jahx06\\.julia\\v0.4\\JuliaFEM\\src\\interpolate.jl:28", + " in dinterpolate at C:\\Users\\jahx06\\.julia\\v0.4\\JuliaFEM\\src\\elements.jl:208", + " in calculate_normals! at In[9]:36", + " in calculate_normals! at In[9]:33", + " [inlined code] from In[9]:61", + " in anonymous at no file:0" + ] }, { "name": "stdout", "output_type": "stream", "text": [ + "\n", "Converged. p = [2.110067133637596,1.2802355090457196,1.5675520985912952,-3.7904167887303712]\n" ] } @@ -161,14 +170,30 @@ "#x2 = linspace(0, 3, nm) + rand(nm)*0.1 + 0.4\n", "#y2 = 0.5*rand(nm) + 0.5\n", "\n", + "\"\"\"\n", + "Calculate element local normal field\n", + "\"\"\"\n", + "function calculate_normals!(el::Seg2, time, field_name=\"normal\")\n", + " new_fieldset!(el, field_name)\n", + " normal_field = Field(time, Vector[])\n", + " for xi in Vector[[-1.0], [1.0]]\n", + " tangent = dinterpolate(el, \"geometry\", xi, time)\n", + " normal = [0 -1; 1 0]*tangent\n", + " normal /= norm(normal)\n", + " push!(normal_field.values, normal)\n", + " end\n", + " #add_field!(el, field_name, normal_field)\n", + "end\n", + "\n", "function create_elements(X, sid=0)\n", " Γ = []\n", " nnodes = size(X, 2)\n", " nelements = nnodes-1\n", " for i=1:nelements\n", " con = sid+[i, i+1]\n", - " el = PSeg(con)\n", - " set_field(el, :Geometry, Vector[X[:, i], X[:, i+1]])\n", + " el = Seg2(con)\n", + " new_fieldset!(el, \"geometry\")\n", + " #add_field!(el, \"geometry\", Field(0.0, Vector[X[:, i], X[:, i+1]]))\n", " push!(Γ, el)\n", " end\n", " return Γ\n", @@ -176,10 +201,70 @@ "Γ₁ = create_elements([x1 y1]', 0)\n", "Γ₂ = create_elements([x2 y2]', nsl);\n", "\n", - "# calculate and average normals like se did in last notebook\n", "for el in [Γ₁; Γ₂]\n", - " calculate_normals!(el)\n", + " calculate_normals!(el, 0.0)\n", "end\n", + "\n", + "figure(figsize=(10, 3))\n", + "plot(x1, y1, \"-ko\")\n", + "plot(x2, y2, \"-ko\")\n", + "axis(\"equal\")\n", + "ylim(-0.1, 4.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Calculating surface normals\n", + "\n", + "- must have unique normal" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged. p = [4.2905787641613236,5.624156522861784,0.04706336239279979,6.561746550428803]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XtcVVX+//HX4X5HRU2JME0wM40k1LCLZt/MBLwEKDOa6ZTNTFdsUrMpv2XONFaj4zg2U7/SQY3xK2kKWpqW5TimJpaVElIhGqCCoSgXuZzfH2egkNsBgQ3nvJ+Px3moa++192fXFnmz1l7bZDabzYiIiIiIiNgAB6MLEBERERERaSkKOCIiIiIiYjMUcERERERExGYo4IiIiIiIiM1QwBEREREREZuhgCMiIiIiIjZDAUdERERERGyGAo6IiIiIiNgMBRwREREREbEZCjgiIiIiImIznJqy89dff83//u//kpqaSm5uLm5ubvTr14+HH36YX/7ylw32XblyJTNmzKhzW25uLt27d6/VnpeXx9atW7n66qtxd3dvSqkiIiIiImJDiouLyczMZPTo0XTt2rXe/ZoUcLKysjh//jz3338//v7+FBUVkZSUxNSpU8nMzOSZZ55p9BgLFiygd+/eNdp8fX3r3Hfr1q1MmTKlKSWKiIiIiIgNW716dYODKyaz2Wy+nBNUVlYSGhrKmTNnOHbsWL37VY3gfPbZZwwePNiqY+/evZtbbrmF1atX079//8sps0OKj49n8eLFRpchBtN9ILoHBHQfiIXuAwH7vQ+OHDnClClT+Pe//83w4cPr3a9JIzh1cXBwICAggMLCQqv2N5vNFBYW4uHhgaOjY4P7Vk1L69+/v9WhyJb4+vra5XVLTboPRPeAgO4DsdB9IKD7oLFHV5q1yEBRURF5eXl8++23LF68mK1btzJ79myr+o4cORJfX188PT0ZN24cGRkZzSlBRERExC4cOXKEoKAgXF1dSUlJwdXVlaCgII4cOWJ0aSLtUrNGcGbNmsXrr79uOYCTE0uXLmXmzJkN9vH09GT69OmMHDkSHx8fPvvsM/785z8THh5OamoqAQEBzSlFRERExGZ9/fXXhISEUF5eXt128eJFMjIyGDRoEIcOHbLLafwiDWlWwImPjyc2Npbs7GzWrFnDI488gru7O9OmTau3T0xMDDExMdV/joqKYvTo0dx2220sXLiQ1157rTmliIiIiNis8ePH1wg3P1deXk5UVBRHjx5t46pE2rdmBZx+/frRr18/AKZMmcLo0aN54okniI2NbdJyzsOHD2fo0KFs3769OWXYvLi4OKNLkHZA94HoHhDQfWCvsrKyLmu72CZ9PWjYZa+iBvD666/z61//mtTUVEJCQprUNzY2lg8//JC8vLxa21JTUwkNDeW2226rtZR0XFyc/ueKiIiITXNxcaGsrKze7c7Ozly8eLENKxJpG4mJiSQmJtZoO3v2LJ988gkHDhxocJGFy15FDSwv3QHLimpN9d1339GtW7cG91m8eLFdrxQhIiIi9slkMl3WdpGOqq7BjKrBj8Y0KZGcPn26VltZWRkJCQn4+fkxYMAAAHJyckhLS6sxZ7Suvlu2bCE1NZW77767KWWIiIiI2LySEnB0DGxwn8DAhreL2KMmjeDMnDmTwsJCbrvtNvz9/cnNzWXNmjWkp6ezYsWK6vfaPP300yQkJJCZmVn9Fy88PJzBgwcTGhqKr68vqampvPXWWwQGBjJv3ryWvzIRERGRDqq4GMaNg8rKTTg6DqKiou6FBh544IE2rkyk/WtSwJk8eTJvvvkmr732Gvn5+fj4+DB06FCWLVvGqFGjqvczmUy1hkwnT57M5s2b2bZtG0VFRfj7+/PQQw8xf/78RqeoiYiIiNiLoiKIioI9e+C99/rTo8choqKiyMrKwmw2YzKZuOqqq7jhhhuYN28egYGBei5Z5GdaZJGB1lI1z66xB4lEREREbMGFCxAZCfv2webNcPvt9e9bWVnJjBkzWLVqFYmJicTGxrZdoSIGsDYbtMgiAyIiIiJyec6fh7FjITUV3n8fbrml4f0dHBx48803qaio4Be/+AWOjo7ce++9bVOsSDvW9GXPRERERKRFFRbC3XfDwYOwdWvj4aaKo6MjK1euJDY2lsmTJ7Nhw4bWLVSkA1DAERERETHQ2bMwejR8+SVs2wbh4U3r7+joSEJCAhMnTiQ2NpZNmza1TqEiHYQCjoiIiIhBCgos4ebIEdi+HYYNa95xnJycWL16NePGjSM6OprNmze3bKEiHYgCjoiIiIgBfvwR/ud/ID3dEm7Cwi7veM7OziQmJjJ27FgmTpzI+++/3zKFinQwCjgiIiIibSw/H0aNgu+/hw8/BCtezm4VZ2dn1q5dy+jRoxk/fjzbtm1rmQOLdCAKOCIiIiJtKC/PEm6OH7eEm5CQlj2+i4sL69at484772TcuHHs2LGjZU8g0s4p4IiIiIi0kdOn4Y47IDsbPvoIBg1qnfO4urqSlJTEiBEjiIyMZN++fa1zIpF2SAFHREREpA2cPAkjR8KpU7BzJ1x/feuez83NjQ0bNhAfH8+AAQNa92Qi7Yhe9CkiIiLSynJzLSM3BQWWcHPttW1zXjc3NxYuXNg2JxNpJxRwRERERFpRdrYl3BQWWsJNcLDRFYnYNgUcERERkVbyww+WaWnFxfDxx9C3r9EVidg+BRwRERGRVnD8uCXclJVZwk2fPkZXJGIfFHBEREREWtixY5ZwYzZbpqX17m10RSL2Q6uoiYiIiLSg77+H228Hk0nhRsQICjgiIiIiLeS772DECHBysoSbXr2MrkjE/ijgiIiIiLSAjAzLyI2rq+WZm6uuMroiEfukgCMiIiJymdLTLeHG09MSbq680uiKROyX1QHn66+/JiYmhmuuuQZPT0/8/PwIDw9nzZo1VvUvKChg5syZdOvWDS8vL+644w4OHjzY7MJFRERE2oO0NEu48fW1TEvr2dPoikTsm9WrqGVlZXH+/Hnuv/9+/P39KSoqIikpialTp5KZmckzzzxTb9/KykrGjh3LoUOHmD17Nn5+fixfvpwRI0Zw4MAB+mpReBEREemADh+2vMSza1fYsQOuuMLoikTEZDabzc3tXFlZSWhoKGfOnOHYsWP17vd///d/TJ48maSkJCZOnAhAXl4ewcHBjBkzpt5RoNTUVEJDQzlw4ACDBw9ubpkiIiIiLe6rr2DUKEuo2bEDunUzuiIR22ZtNrisZ3AcHBwICAjA2dm5wf2SkpLo0aNHdbgB6Nq1K7GxsWzcuJGysrLLKUNERESkTR06ZHnPjb8/fPihwo1Ie9LkgFNUVEReXh7ffvstixcvZuvWrcyePbvBPgcPHqwzZYWFhVFUVER6enpTyxARERExxOefW8JNYKBl5KZrV6MrEpGfa3LAmTVrFt27dycoKIg5c+awdOlSZs6c2WCfnJwcetbxxF1VW3Z2dlPLEBEREWlzqamWZ2769IHt26FLF6MrEpFLWb3IQJX4+HhiY2PJzs5mzZo1PPLII7i7uzNt2rR6+5SUlODq6lqr3c3NDYDi4uKmliEiIiLSpvbvh7vuguBg2LoVOnUyuiIRqUuTA06/fv3o168fAFOmTGH06NE88cQTxMbG4u7uXmcfd3d3SktLa7WXlJRUbxcRERFpr/buhdGj4brr4L33LEtCi0j71OSAc6l7772XDz74gG+++YaQkJA69+nZs2ed09BycnIA8Pf3b/Ac8fHx+F7ylSQuLo64uLhmVi0iIiJinT17LOFm0CBLuPH2NroiEduXmJhIYmJijbazZ89a1feyA07V9DIHh/of5wkJCWHXrl2YzWZMJlN1+969e/H09CQ4OLjBcyxevFjLRIuIiEib270b7r4bBg+GzZvBy8voikTsQ12DGVXLRDfG6kUGTp8+XautrKyMhIQE/Pz8GDBgAGAZlUlLS6O8vLx6v+joaE6ePMn69eur2/Ly8li3bh2RkZGNLjMtIiIi0tY++cQychMWBlu2KNyIdBRWj+DMnDmTwsJCbrvtNvz9/cnNzWXNmjWkp6ezYsUKHB0dAXj66adJSEggMzOTwMBAwBJwhg0bxvTp0zl8+DB+fn4sX74cs9nM888/3zpXJiIiItJMH30EERFw882waRN4eBhdkYhYy+qAM3nyZN58801ee+018vPz8fHxYejQoSxbtoxRo0ZV72cymWpMQwPL9LUtW7bw1FNPsXTpUoqLixkyZAgJCQkEBQW13NWIiIiIXKYdOyAyEm65BTZuBK2FJNKxmMxms9noIupTNc/uwIEDegZHREREWt22bTBuHIwYARs2wH/faCEi7YC12aDJL/oUERERsUXvvQdRUTBqFLz7rsKNSEelgCMiIiJ2LyUFxo+3vMjznXegjveTi0gHoYAjIiIidm3jRpg4EcaOhaQkhRuRjk4BR0REROzWhg0QHW2ZmrZ2Lbi4GF2RiFwuBRwRERGxS0lJEBtrGb1JTAS9lk/ENijgiIiIiN1ZuxYmT4aYGFizRuFGxJYo4IiIiIhdeftt+MUvLJ9Vq8DJ6rcCikhHoIAjIiIidiMhAaZOtXxWrABHR6MrEpGWpoAjIiIidmHlSrj/fsvnrbcUbkRslQKOiIiI2Lw334QZM+DBB+GNN8BB3wGJ2Cz99RYRERGb9vrr8MAD8Otfw2uvKdyI2Dr9FRcRERGb9dpr8NBD8Oij8Le/KdyI2AP9NRcRERGbtGwZ/Pa38Pjj8Je/gMlkdEUi0hYUcERERMTm/OUvllGbJ5+ExYsVbkTsiQKOiIiI2JRXX4UnnoA5c+DllxVuROyNAo6IiIjYjEWL4He/g2eegT/+UeFGxB4p4IiIiIhN+MMfLKM2zz0HCxYo3IjYKwUcERER6fAWLLCM2jz/vOWjcCNivxRwREREpENISUnB2dkZk8lU/XF2dmby5BSeew5efNEyeiMi9s3qgLN//34eeeQRBgwYgJeXF7169WLSpEkcPXq00b4rV67EwcGhzs+pU6cu6wJERETE9m3cuJHIyEjKy8trtJeXl7N2bST33ZfCM88YVJyItCtO1u74pz/9iT179hATE8OgQYPIyclh2bJlDB48mE8//ZQBAwY0eowFCxbQu3fvGm2+vr5Nr1pERETsSnR0dIPb3357Av/8Z1kbVSMi7ZnVAefJJ58kLCwMJ6efukyaNImBAwfy0ksvsWrVqkaPMWbMGAYPHty8SkVERMRuXTpy09TtImI/rJ6idvPNN9cINwB9+/bluuuuIy0tzapjmM1mCgsLqaioaFqVIiIiIiIiVrisRQbMZjMnT56ka9euVu0/cuRIfH198fT0ZNy4cWRkZFzO6UVERERERGqweopaXdasWUN2djYvvvhig/t5enoyffp0Ro4ciY+PD5999hl//vOfCQ8PJzU1lYCAgMspQ0RERGzU2ZKz/PHff4TOwI/173fpLBMRsV/N/mqQlpbGww8/THh4ONOmTWtw35iYGGJiYqr/HBUVxejRo7nttttYuHAhr732WnPLEBERERtUXlnOGwfe4Lmdz1FUVkTcc3EkxifWu/+GDRvasDoRac+aNUUtNzeXsWPH0rlzZ5KSkjA1421aw4cPZ+jQoWzfvr05JYiIiIgNMpvNbDm6hUGvDeLhLQ8zNmgs6Y+k8/YTb5OcnFxrpKbqe5AdO3YYUa6ItENNHsE5e/YsY8aM4dy5c+zatYsePXo0++QBAQGkp6c3ul98fHyt5aTj4uKIi4tr9rlFRESkfTl08hC/2/Y7PvjuA0ZcPYLVE1czuOdPq69GRERQVlZzKejKykquvvpqlixZQkREBKNGjWrrskWkFSQmJpKYWHPU9uzZs1b1NZnNZrO1JyopKeGuu+7i4MGDbN++naFDhzat0kvcdNNNXLhwgSNHjtS5PTU1ldDQUA4cOKDlpUVERGxU7vlcnv3wWd48+CZ9u/TllbteITI40uoZIidOnKBPnz44OTmRnZ1Np06dWrliETGCtdnA6ilqFRUVTJo0ib1797Ju3bp6w01ubi5paWk11qM/ffp0rf22bNlCamoqd999t7UliIiIiA0pKivixU9epO/Svrxz5B2W3L2Er377FVH9opo0/T0gIIBVq1ZRXFzM8OHDW7FiEekImvSiz+TkZCIjI8nLy2P16tU1tk+ZMgWAuXPnkpCQQGZmJoGBgQCEh4czePBgQkND8fX1JTU1lbfeeovAwEDmzZvXgpcjIiIi7V2luZK3v3ybeTvmkXs+l0eHPMrvb/s9nd07N/uYkyZNYtOmTbz99ts89thjLF26tAUrFpGOxOqA88UXX2AymUhOTiY5ObnGNpPJVB1wTCZTrZ+6TJ48mc2bN7Nt2zaKiorw9/fnoYceYv78+XTr1q0FLkNEREQ6gl3HdvHktifZn72fif0n8qc7/0TfLn1b5NirVq1i9+7d/PWvf2Xs2LGMHj26RY4rIh1Lk57BaWt6BkdERMQ2fHvmW+Zsn8M7R94htGcoi0cv5tZet7b4ebKzs7n66qtxdHTkhx9+oEuXLi1+DhExRos/gyMiIiLSVD8W/8iTW5+k/9/6s/eHvSSMT2Dfg/taJdwA+Pv78/bbb1NSUqLncUTslAKOiIiItLiyijKW7l1K37/25R8H/sFztz/HN498w9QbpuJgat1vP6Kjo5k6dWr1S8lFxL4o4IiIiEiLMZvNbPpmE9e/dj3xW+OZeO1EMh7L4Pe3/R4PZ482q2PlypVcffXVLF++nPfee6/NzisixlPAERERkRbxee7njEoYxbh/jeMqn6tInZnKG1Fv0MOr+S8Fby4HBwf27NmDi4sLEydOJC8vr81rEBFjKOCIiIjIZckuzGbGxhkM/sdgcs7nkBKXwgdTP+CGHjcYWlePHj30PI6IHbJ6mWgRERGRn7tw8QKv/OcVFv1nER7OHiy7ZxkPDn4QZ0dno0urdu+993L//ffz8ccfU1JSgpubm9EliUgrU8ARERGRJqk0V5LwRQLPfPgMeUV5PD70cebdOo9Obp2MLq1Ob775JmCZtiYitk8BR0RERKy2M3Mns7bO4mDuQWIHxPLSqJfo3bm30WU1SMFGxL4o4IiIiEij0vPTmf3BbDZ+s5GhVw5l94zdhF8VbnRZIiK1KOCIiIhIvfKL8nnh4xdY/tly/L39Sbw3kUkDJmEymYwuTUSkTgo4IiIiUsvFiov8bd/feOGTF6iorGDByAU8PvRx3J3djS5NRKRBCjgiIiJSzWw2827au8zePpvvfvyOB258gBdGvsAVXlcYXZqIiFUUcERERASAA9kHmLVtFp8c+4TR14xmw6QNXN/9eqPLEhFpEgUcERERO3f87HGe+fAZVh1axYBuA3j/l+8zuu9oo8sSEWkWBRwRERE7df7ief707z/xyp5X8HH14R8R/2DGjTNwctC3ByLScekrmIiIiJ2pqKxg5ecr+f1Hv+fH4h958uYnmXPLHHxcfYwuTUTksingiIiI2JHt323nyW1PcujkIX4x8Bf84Y4/0KtTL6PLEhFpMQo4IiIiduDI6SM89cFTbD66mfCrwvn0V58yNGCo0WWJiLQ4B2t33L9/P4888ggDBgzAy8uLXr16MWnSJI4ePWpV/4KCAmbOnEm3bt3w8vLijjvu4ODBg80uXERERBp3+sJpHtnyCANfG8jh04dZF7OOf0//t8KNiNgsq0dw/vSnP7Fnzx5iYmIYNGgQOTk5LFu2jMGDB/Ppp58yYMCAevtWVlYyduxYDh06xOzZs/Hz82P58uWMGDGCAwcO0Ldv3xa5GBEREbEoLS9l6d6lLNy1EDNm/jjqjzw69FHcnNyMLk1EpFVZHXCefPJJwsLCcHL6qcukSZMYOHAgL730EqtWraq3b1JSEnv27CEpKYmJEycCEBsbS3BwMPPnz2fNmjWXcQkiIiJSxWw2k3Q4iTnb55B1Notf3/Rr5t8+n26e3YwuTUSkTVgdcG6++eZabX379uW6664jLS2twb5JSUn06NGjOtwAdO3aldjYWFavXk1ZWRnOzs5NKFtEREQutffEXmZtm8V/jv+HsUFj2fLLLVzb9VqjyxIRaVNWP4NTF7PZzMmTJ+natWuD+x08eJDBgwfXag8LC6OoqIj09PTLKUPE5hw5coSgoCBcXV1xcXHB1dWVoKAgjhw5YnRpItIOHSs4xi/e+QXD3hzG+Yvn+WDqB6T8IkXhRkTs0mUFnDVr1pCdnc2kSZMa3C8nJ4eePXvWaq9qy87OvpwyRGzK119/zaBBg8jIyODixYuUlZVx8eJFMjIyGDRokEKOiFQ7V3qOp7c/Tb9l/fgo8yPejHqT1Jmp3NnnTqNLExExTLOXiU5LS+Phhx8mPDycadOmNbhvSUkJrq6utdrd3CwPOhYXFze3DBGbM378eMrLy+vcVl5eTlRUlNWrF4qIbSqvLOetg2/x7EfPUlhayFPhTzF7+Gy8Xb2NLk1ExHDNCji5ubmMHTuWzp07k5SUhMlkanB/d3d3SktLa7WXlJRUbxcRi6ysrMvaLiK2bWvGVn73we/46tRXTB00lT+M+gMBPgFGlyUi0m40OeCcPXuWMWPGcO7cOXbt2kWPHj0a7dOzZ886p6Hl5OQA4O/v32D/+Ph4fH19a7TFxcURFxfXhMpFOgaz2XxZ20XENn116it+t+13bP12K7cG3sr+B/dzk/9NRpclItIqEhMTSUxMrNF29uxZq/o2KeCUlJQQGRlJRkYG27dv59prrXt4MSQkhF27dmE2m2uM9uzduxdPT0+Cg4Mb7L948eI6FykQsUWNjYiay8yk3pKKX4QffhF+eA7wbLSPiHRcpy6c4rmPnuON1Dfo07kP62PXM/7a8fp7LyI2ra7BjNTUVEJDQxvta/UiAxUVFUyaNIm9e/eybt06hg6t+w3Iubm5pKWl1XiGIDo6mpMnT7J+/frqtry8PNatW0dkZKSWiBb5mcDAwAa3B3QLwKWbC8cWHOOzgZ/xae9PSX8knfz386koqWijKkWkpaSkpODs7IzJZKr+ODs7s37jev6464/0XdqXtV+v5ZX/eYWvf/s1E/pPULgREWlAk170mZycTGRkJHl5eaxevbrG9ilTpgAwd+5cEhISyMzMrP5GLTo6mmHDhjF9+nQOHz6Mn58fy5cvx2w28/zzz7fg5Yh0fJs2bWLgwIFUVNQOK05OTmz5eAv9+/enoqSCgp0F5Kfkc2bzGbL/lo2DhwOd7+xsGd0Z64erf+3FPUSk/di4cSPjx4+v1V5eXs694++FX8LjUx7n2duexc/Dz4AKRUQ6HqsDzhdffIHJZCI5OZnk5OQa20wmU3XAqfrp0885ODiwZcsWnnrqKZYuXUpxcTFDhgwhISGBoKCgFrgMEdvRv39/Xn75ZWbNmlU9umkymQgMDGTTpk30798fAEc3R/zu9sPvbj/MfzVTdLiI/M355Cfnk/7rdKgEr1Cv6qls3oO9MTnop74i7Ul0dHSD2x3XOrJk9ZI2qkZExDaYzO34ieWqeXYHDhzQMzhiV6ZNm0Zqaipffvlls/qXnSnjzHtnLKM775+hvKAc5yuc8RtrCTud7+yMk3ezV4kXkRZizVSzdvzPtIhIm7I2G+g7HJF2pqKigi1btvDAAw80+xjOXZy54pdXcMUvr6CyrJJz/zlHfko++Sn55L6Vi8nFRKcRnapHd9x7a6l2kTZx5gx88gns3Ak7d+IMlBldk4iIjVHAEWln9u/fT15eHhERES1yPAdnBzrd3olOt3fimpevofjbYkvY2ZzPt09+S8ZjGXhc51Eddnxu9sHByer1R0SkIZcEGg4dArMZeveG22+n4osvjK5QRMTmKOCItDMpKSl06dKFYcOGtcrx3a9xJ+DxAAIeD6C8sJwft/1I/uZ8clfmcnzRcZw6O9FlTBf8IvzoMroLzl20yqGI1eoLNFdfDSNGQHw83H675c+Aw+rVVP5s1dFLOTnpn2kRkabSV06RdiYlJYV77rkHR0fHVj+Xk7cT3e7tRrd7u2GuNFP4WSH5yZapbKfePgWO4Bvui1+kZXTH41oPLU8r8nONBZonnoCRI6FXrzq7b9iwgcjIyHoPv2HDhtaoWkTEpingiLQjx48f54svvuDpp59u83ObHEz4DPHBZ4gPvRf0pvSHUsuqbJvzyZyfyXezv8Otj1v1VLZOt3XCwVVT2cTONDLljFmzLL/WE2guFRERQXJyMhMmTKjx/jiA9evXt9hUVRERe6KAI9KObN68GUdHR0aPHm10Kbhe6Yr/TH/8Z/pTUVxBwYcF5G/OJ29DHj8s/QFHL0c6/4/lnTtd7umCaw+9c0dsUBOnnDVHREQEZWU/LTXwz3/+k/vvv59Dhw4xYcKEy7wAERH7o2WiRdqRiIgILly4wEcffWR0KfUym81c+PJC9aps5z49B2bwDvO2jO5E+uEV4qWpbNIxNRZobr+9wSlnLaGyshIvLy86depEdnZ2q51HRKSj0TLRIh1MUVERO3bs4MUXXzS6lAaZTCa8BnnhNciLXvN6cfH0xep37hx/9TiZ8zNx8Xf56Z07ozrj6Nn6zxOJNEsLTzlrCQ4ODkRFRbF27VpSU1P1Az4RkSZSwBFpJz788ENKSko63Jx7l24u9LivBz3u60HlxUrO/vts9ehOzhs5mFxNdL7DMpXNb6wfbr3cjC5Z7FkbTDlrCYsWLWLt2rXMmTOHDz74wNBaREQ6GgUckXYiJSWFvn37EhwcbHQpzebg4kDnOzrT+Y7O9P1zX4rSi6rfuZPxeAZHHz6K50DPn965M9QHk6OmskkrusxVzowSGBhIUFAQO3fupLy8XMtFi4g0gb5iirQDZrOZlJQUYmJibOrZFY9gDzxmeXDVrKsoP1vOmW1nqkd2sv6YhZOfE373/Hcq212dce6kd+7IZWqHU86aa/bs2Tz44IO8/PLLhqysKCLSUWmRAZF24PPPP+fGG29k+/btjBo1yuhyWp25wsy5feeq37lz4csLmJxM+N7iWz2649HPw+gypSNoaIRm5EhLmGkHU86ao7KyEg8PD7p27cqJEyeMLkdExHBaZECkA0lJScHb25tbb73V6FLahMlW7mdcAAAgAElEQVTRhO/Nvvje7EufP/ShJKuk+p073//+e7793be4B7lXhx3fW3xxcNE7d4QOO+WsORwcHIiIiOCdd97hq6++4vrrrze6JBGRDkEBR6QduO6665g7dy4uLi5Gl2IIt0A3rvzNlVz5myupKKrgx+0/kr85n1NrT3Fi8QkcvR3pMrqL5Z07Y7rg0t0+/zvZJRuactYcL7/8Mu+88w5PPfUU7733ntHliIh0CAo4Iu3AxIkTjS6h3XD0cKRrVFe6RnXFbDZz/vPz1auypU1PA8BnqE/1O3c8B3ra1HNLdq++QNOrl2Vk5oknLCM1HXDKWXP07t2bPn36sH37di02ICJiJc35EJF2y2Qy4X2jN1c/ezWhe0MJzwmn35v9cLnShayXsvjshs/4tNenpP8mnfzN+VQUV9Ton5KSgrOzMyaTqfrj7OxMSkqKQVcktZw5A+++awkuISHQtStMmAAbN8KNN8Jbb8H330NmJqxYAfffbzfhpsqTTz5JeXk5ixcvNroUEZEOQYsMiEiHVFlaScEnBZaFCjbnU/JdCQ7uDnQeZXnnzqdOn3LvA/fW2z85ObnDvXPIJjQ25axqYQAbnXLWHJWVlbi7u3PFFVeQlZVldDkiIobRIgMiYtMcXB3o8j9d6PI/Xej7l74UfVNUHXbSH07n0YpHG+w/YcIEysrK2qha25OSksKECRMoLy+vbnNycmLDhg01g2Njq5y1kxdrtmcODg7cc889vPvuuxw5coT+/fsbXZKISLumgCMiHZ7JZMLzWk88r/Uk8KlAyn4sI6dLToN9ysvLGbR/P84mEy4ODjibTNW/d7n09z9rq/p9VR+XS/araqtvv5+fp77zOf53Ol17tXHjRsaPH1+rvby8nMjISJLnzSPiwgWbX+WsLS1atIhTp07h4KCZ5SIijWlSwLlw4QKLFi1i79697Nu3j4KCAlasWMG0adMa7bty5UpmzJhR57bc3Fy6d+/elFJEROrl3NkZM43Pvh3RqRNlZjMXKyurf71oNlNmNlNYUUFZeflP2y7Zr6qt7L99LlZWUtHoGZtwDZeEokuDUK1AVsf21gph0ZMmNVj7hD/8gTI7WeWsrQQFBbF7926jyxAR6RCaFHBOnz7NggUL6NWrFyEhIezcubPJP2VcsGABvXv3rtHm6+vbpGOIiLSEpUFBLXq8yv+Go7KfBaW6wlNdQamh8FR1zNI6jn3pfj8PZtaer8nBrKLhHuUA333XzP+KIiIil6dJAcff3796tOXAgQOEhYU1+YRjxozRggEi0uqcnJxqPB9S1/aW5mAy4Woy4drBphFVms2U/zcMlV4ayOoISrc38N9VRETEaE36F97FxaV6KllzF18zm80UFhbi4eGBo6Njs44hItKYDRs2EBkZ2eB2sXComp4GeOrrsoiIdHBt/mPGkSNH4uvri6enJ+PGjSMjI6OtSxAROxAREUFycnKtH6Q4OTlpiejL1Njol15GKSIiRmqzgOPp6cn06dNZvnw57777LrNnz2bHjh2Eh4dz4sSJtipDROxIREQEixYtAmDLli2YzWbKysoUbi5TY6NfGh0TEREjtVnAiYmJ4c0332TKlClERUXxwgsvsHXrVvLz81m4cGFblSEidmb//v0A3H777QZXYjuqRscuHanR6JiIiLQHhs4jGD58OEOHDmX79u1GliEiNuzIkSO4uLjg4eFhdCk2JSIigrKyMkaPHs0HH3xAZWWl0SWJiIgA7eBFnwEBAaSnpze4T3x8fK2lpOPi4oiLi2vN0kTEBpw4cQI/Pz+jy7BZffr0wWw2c+rUKb3PTEREWkxiYiKJiYk12s6ePWtVX8MDznfffUe3bt0a3Gfx4sVaWlpEmqWgoIAhQ4YYXYbNuvbaawE4cOAAY8aMMbgaERGxFXUNZqSmphIaGtpo31Z5Bic3N5e0tLQa76A4ffp0rf22bNlCamoqd999d2uUISJ27tSpU1RUVHD99dcbXYrNGjhwIABffvmlwZWIiIhYNHkEZ9myZRQUFJCdnQ3Apk2byMrKAuCxxx7Dx8eHuXPnkpCQQGZmJoGBgQCEh4czePBgQkND8fX1JTU1lbfeeovAwEDmzZvXgpckImKxY8cOAG6++WaDK7FdN910EwBpaWkGVyIiImLR5IDz6quvcuzYMQBMJhMbNmxg/fr1mEwm7rvvPnx8fDCZTJhMphr9Jk+ezObNm9m2bRtFRUX4+/vz0EMPMX/+/EanqImINMd//vMfAEaNGmVwJbar6mt+1b8LIiIiRmtywPn+++8b3WfFihWsWLGiRtuCBQtYsGBBU08nItJsX3zxBQ4ODtUjydI63Nzcqkf1RUREjNZm78EREWlr33//PT4+PkaXYfN8fX3Jy8szugwRERFAAUdEbFheXh5XXnml0WXYvG7dulFYWGh0GSIiIoACjojYqIsXL1JSUkK/fv2MLsXmBQQEUFpaanQZIiIigAKOiNioqgUGqlb5ktbTu3dvAD2HIyIi7YICjojYpI8//hiAESNGGFuIHejfvz9geQGbiIiI0RRwRMQm+fr6csUVVxAWFmZ0KTYvJCQE0Ms+RUSkfVDAERGb9MQTT5Cbm4uTU5NXw5cmqgo433zzjcGViIiIKOCIiMhl8vLyonfv3nrfkIiItAv60aaIiFy27777zugSREREAI3giIiInUpJScHZ2RmTyVT9cXZ2JiUlxejSRETkMijgiIiI3dm4cSORkZGUl5fXaC8vLycyMlIhR0SkA1PAERERuxMdHd3g9gkTJrRRJS1Lo1IiInoGR0RE7EBxMRw7BpmZlk95ubnB/cvLy5k4cSKenp54enri4eFR41drf+/o6Ngm1weWUanx48fXeS2RkZEkJycTERHRZvWIiBhFAUdERDq8SwPMpZ+TJ3/a15I5TFYcs5j8/HwuXLhQ/SkqKuLChQuUlJRYVZeLi0uDIak5oennx/n5MujWjEqVlZVZVbeISEemgCMiIu1eUwPMVVfB1VdD//4wZozl91WfK68EZ+fyS09Ry3vvvVfvtoqKCoqKiqoDz8/DT9Xv69r287YLFy60SICqCjyXPk90qca2i4jYCgUcERFpUEpKChMmTKjxDbKTkxMbNmxosSlPzQ0w114Ld99t+X3v3j8FmMbe7+rk5NTgN/yNvSDW0dERb29vvL29rbi6pquoqKC4uLhWIGooUL3wwguNHjctbTre3kPw8RmCp+dAHBxcWqV+EREjKeCIiEi9Wuq5jssdgenVyxJgevcGf//GA0xjNmzYQGRkZIPbjeTo6IiXlxdeXl5W97Em4Jw/f4iTJ1djNpdjMrni7X1jdeDx9h6Cu3tfTKbGp++JiLRnCjgiIlIva5/raG6Aue66n6aQ9epl+TUg4PIDTGMiIiJITk5u9ZGptmTNqNRNNx2goqKY8+c/p7BwH+fO7ePMmS388MPS/+7T6WeBJwxv7yG4uvZoq0sQEWkRTfon5MKFCyxatIi9e/eyb98+CgoKWLFiBdOmTbOqf0FBAbNnz2bDhg0UFxczZMgQXn31VW688cZmFS8iIq3Lmuc6evSoHWACAiyjLfU9A9PaAcYaERERNvXQvbWjUo6O7vj63oyv783V28rKzlBYuJ9z5/ZRWLiP7OzXKSt7EQBX16tqjPJ4e4fi5NQ6U/NERFpCk/6JOX36NAsWLKBXr16EhISwc+dOq4eyKysrGTt2LIcOHWL27Nn4+fmxfPlyRowYwYEDB+jbt2+zLkBERIz1m9/89AxMr17tJ8DYm8sZlXJ27kKXLqPp0mU0AGazmdLS49WB59y5fWRmvkBl5QXAhIdH//8GnqH4+ITpeR4RaVea9E+Qv78/ubm5dO/enQMHDhAWFmZ136SkJPbs2UNSUhITJ04EIDY2luDgYObPn8+aNWuaVrmIiLS4sooydh/fTUp6CsnpyeAIVDTcZ/78NilNrNBSo1Imkwk3t0Dc3ALp3t0yTdFsruDChcM1Rnr0PI+ItEdNCjguLi50794dsPx0pymSkpLo0aNHdbgB6Nq1K7GxsaxevZqysjKcnZ2bdEwREbl8Z4rP8H7G+ySnJ/N+xvsUlBTQw6sHEUERZDhnUFlRWW/fxlYbE9thMjni5TUQL6+B9Ow5A+C/z/McrA48NZ/n6Yy3d1h14PHxGYKLyxVGXoKI2Ik2+5fp4MGDDB48uFZ7WFgYr7/+Ounp6QwYMKCtyhERsVtms5m0vDRS0lNIOZrC7qzdVJgrGNxzMI8NeYyI4AhC/UNxMDkwbt24dr3amBjL8jxPOL6+4dVtLf08z5EjR4iKiiIrKwuz2YzJZCIwMJBNmzbRv3//Vr0+EemY2izg5OTkMGLEiFrtPXv2BCA7O1sBR0SklVysuMiuY7uqp559++O3uDu5M6rPKJaPXc7YoLFc6XNlrX62uNqYtK76n+fZW8/zPNfVGOWxPM9jmdHx9ddfExISUmuxi4yMDAYNGsShQ4cUckSkljYLOCUlJbi6utZqd3NzA6C4uLitShERsQt5RXm8d/Q9ktOT2frtVs6VnuNK7yuJCI5gSfAS7uh9Bx7OHo0ex9ZWG5O2VfN5nhjg58/z7OPcuf0UFu4lNzcBqMDBwQ0vrxC8vYcQEfGvelfyKy8vJyoqiqNHj7bh1YhIR9BmAcfd3Z3S0tJa7SUlJdXbRUSk+cxmM4dPHyY5PZmU9BT2nNhDpbmSMP8wnrz5SSKDIwnpEaIHv8VwNZ/n+RUAFRVF/32eZz+FhfvJz9/MiROnGjxOVlZWW5QrIh1MmwWcnj17kp2dXas9JycHsKzQVp/4+Hh8fX1rtMXFxREXF9eyRYqIdDCl5aV8fOxjkr9JJuVoCpkFmXg6e3Jnnzt5PeJ1xgaPpYeXXtQo7Z+jowe+vsPx9R1e3WYyuQD1jx42dcEjEek4EhMTSUxMrNF29uxZq/q2WcAJCQlh165d1Q8IVtm7dy+enp4EBwfX23fx4sV1LlAgImKPTl04xZajW0hOT2bbt9s4f/E8gb6BRARFEBEcwcjeI3FzcjO6TJHL1thoo0YjRWxXXYMZqamphIaGNtq3VQJObm4uBQUF9O3bt3oJ0ejoaJKSkli/fj333nsvAHl5eaxbt47IyEgtES0iUg+z2cyXp74k+ZtkktOT2ffDPgCGBgxl7vC5RPaLZGD3gfpmT2xOYGAgGRkZDW4XEblUkwPOsmXLKCgoqJ5utmnTpuo5sI899hg+Pj7MnTuXhIQEMjMzq7/4REdHM2zYMKZPn87hw4fx8/Nj+fLlmM1mnn/++Ra8JBGRjq+kvISPvv+o+nma4+eO4+XixehrRvObm37DmKAxdPfsbnSZIq1q06ZNDBo0qM6FBpycnNi0aZMBVYlIe9fkgPPqq69y7NgxwDI0vGHDBtavX4/JZOK+++7Dx8cHk8lU6yeJDg4ObNmyhaeeeoqlS5dSXFzMkCFDSEhIICgoqGWuRkSkA8spzKmeevbBdx9QVFZE7069GX/teCKDI7mt1224OtVejVLEVvXv359Dhw7pPTgi0iQmczt+Qq9qnt2BAwf0DI6IdCgpKSmNvjvGbDZzMPdg9btpPsv+DAeTAzcH3ExEcASRwZFc1+06TT0TERHB+mzQZosMiIjYi40bNzJ+/Pha7eXl5URGRvLsP57lZM+TpBxNIbswGx9XH+7uezePDXmMMUFj6OrR1YCqRUREbIMCjohIC4uOjm5w+4JfL6DvX/oSe10skf0iuTXwVpwdtdCKiIhIS1DAERFpYfW9eb2aGdIfSdfUMxERkVbgYHQBIiI2pbISa2KLwo2IiEjrUMAREWkJFRXwr3/BoEFospmIiIhxFHBERC5HeTmsXg3XXw9xcRAQgLOjY4Ndql6ALCIiIi1PAUdEpDnKymDFCujfH6ZOhaAg2LsX3n+ff737boNdN2zY0EZFioiI2B8FHBGRprh4Ed54A4KDYcYMGDgQDhyATZtgyBAAIiIiSE5OrjVS4+TkRHJycvV7cERERKTlaZ6EiIg1SkrgrbfgpZfgxAmIjraEmoED69w9IiKCsrKyNi5SRERENIIjItKQoiL4y1+gTx949FG49Vb46iv4v/+rN9yIiIiIcTSCIyJSl/Pn4e9/h5dfhvx8mDIF5s2zTE0TERGRdksBR0Tk5woL4W9/g1dfhYICmDbNEmz69DG6MhEREbGCAo6ICFjCzNKlsGQJXLhgWUBg7lzo1cvoykRERKQJFHBExL6dOWMJNX/5i2WFtAcfhNmzISDA6MpERESkGRRwRMQ+nT4Nf/4zLFsGFRXwm9/A734HPXsaXZmIiIhcBgUcEbEvJ0/CK6/A8uXg4AC//S08+SR07250ZSIiItICFHBExD788INlRbR//ANcXCA+3vLx8zO6MhEREWlBCjgiYtuOH4c//QneeAM8PGDOHHj8cejc2ejKREREpBUo4IiIbfr+e3jpJVixAnx84Lnn4JFHwNfX6MpERESkFTk0ZefS0lLmzJmDv78/Hh4eDBs2jO3btzfab+XKlTg4ONT5OXXqVLOLFxGpJSPDssRzUBBs2AAvvgiZmfDMMwo3IiIidqBJIzj3338/77zzDvHx8QQFBbFixQruuecePvroI4YPH95o/wULFtC7d+8abb76hkNEWkJaGixcCG+/bVkw4OWX4aGHLNPSRERExG5YHXD27dvH2rVreeWVV5g1axYAU6dO5frrr2f27Nns3r270WOMGTOGwYMHN79aEZFLffWVJdisXQv+/pZ32jzwALi7G12ZiIiIGMDqKWpJSUk4OTkxc+bM6jZXV1d+9atfsWfPHn744YdGj2E2myksLKSioqJ51YqIVPn8c4iOhoEDYc8ey7LP334Ljz6qcCMiImLHrA44Bw8eJDg4GC8vrxrtYWFhAHz++eeNHmPkyJH4+vri6enJuHHjyMjIaGK5ImIvjhw5QlBQEK6urri4uODq6kpQUBBH1q2DcePgxhvh4EHL6mjp6fDrX4Orq9Fli4iIiMGsnqKWk5NDzzre8F3Vlp2dXW9fT09Ppk+fzsiRI/Hx8eGzzz7jz3/+M+Hh4aSmphIQENCM0kXEVn399deEhIRQXl5eoz0jI4NBsbEc6tWL/itXwi9+Ac7OxhQpIiIi7ZLVAae4uBjXOn466ubmVr29PjExMcTExFT/OSoqitGjR3PbbbexcOFCXnvttabULCI2bvz48bXCTZVyIMrJiaPTprVtUSIiItIhWD1Fzd3dndLS0lrtJSUl1dubYvjw4QwdOtSqZaZFxL5kZWU1vP348TaqRERERDoaq0dwevbsWec0tJycHAD8/f2bfPKAgADS09Mb3S8+Pr7WctJxcXHExcU1+Zwi0v6ZzebL2i4iIiIdW2JiIomJiTXazp49a1VfqwPOjTfeyM6dOyksLMTb27u6fe/evQCEhIRYe6hq3333Hd26dWt0v8WLF2t5aRE7YjKZLmu7iIiIdGx1DWakpqYSGhraaF+rp6hFR0dTUVHB66+/Xt1WWlrKihUrGDZsGFdeeSUAubm5pKWl1Zg/f/r06VrH27JlC6mpqdx9993WliAidiIwMPCytouIiIj9snoEZ8iQIcTExPD0009z6tQprrnmGv75z3+SlZXFihUrqvebO3cuCQkJZGZmVn8TEh4ezuDBgwkNDcXX15fU1FTeeustAgMDmTdvXstflYh0aJs2bWLQoEF1LjTg5OTEpk2bDKhKREREOgKrAw5AQkICzz77LKtWreLHH3/khhtuICUlhVtuuaV6H5PJVGv6yOTJk9m8eTPbtm2jqKgIf39/HnroIebPn2/VFDURsS/9+/fn0KFDREVFkZWVhdlsxmQyERgYyKZNm+jfv7/RJYqIiEg7ZTK346d1q+bZHThwQM/giIiIiIjYMWuzgdXP4IiIiIiIiLR3CjgiIiIiImIzFHBERERERMRmKOCIiIiIiIjNUMARERERERGboYAjIiIiIiI2QwFHRERERERshgKOiIiIiIjYDAUcERERERGxGQo4IiIiIiJiMxRwRERERETEZijgiIiIiIiIzVDAERERERERm6GAIyIiIiIiNkMBR0REREREbIYCjoiIiIiI2AwFHBERERERsRkKOCIiIiIiYjOaFHBKS0uZM2cO/v7+eHh4MGzYMLZv325V34KCAmbOnEm3bt3w8vLijjvu4ODBg80qWkREREREpC5NCjj3338/ixcvZurUqSxduhRHR0fuuecedu/e3WC/yspKxo4dS2JiIo899hiLFi3i1KlTjBgxgoyMjMu6ABERERERkSpO1u64b98+1q5dyyuvvMKsWbMAmDp1Ktdffz2zZ89uMOQkJSWxZ88ekpKSmDhxIgCxsbEEBwczf/581qxZc5mXISIiIiIi0oQRnKSkJJycnJg5c2Z1m6urK7/61a/Ys2cPP/zwQ4N9e/ToUR1uALp27UpsbCwbN26krKysmeWLiIiIiIj8xOqAc/DgQYKDg/Hy8qrRHhYWBsDnn3/eYN/BgwfXag8LC6OoqIj09HRryxAREREREamX1QEnJyeHnj171mqvasvOzm6VviIiIiIiItayOuAUFxfj6upaq93Nza16e31KSkqa3VdERERERMRaVgccd3d3SktLa7WXlJRUb2+NviIiIiIiItayOuD07NmzzqlkOTk5APj7+7dKX4D4+HiioqJqfBITE60tvcOyh2uUxuk+EN0DAroPxEL3gYB93AeJiYm1vv+Pj4+3qq/Vy0TfeOON7Ny5k8LCQry9vavb9+7dC0BISEi9fUNCQti1axdmsxmTyVSjr6enJ8HBwQ2ee/HixXUuUmDrEhMTiYuLM7oMMZjuA9E9IKD7QCx0HwjYx30QFxdX6xpTU1MJDQ1ttK/VIzjR0dFUVFTw+uuvV7eVlpayYsUKhg0bxpVXXglAbm4uaWlplJeX1+h78uRJ1q9fX92Wl5fHunXriIyMxNnZ2doyRERERERE6mX1CM6QIUOIiYnh6aef5tSpU1xzzTX885//JCsrixUrVlTvN3fuXBISEsjMzCQwMBCwBJxhw4Yxffp0Dh8+jJ+fH8uXL8dsNvP888+3/FWJiIiIiIhdsjrgACQkJPDss8+yatUqfvzxR2644QZSUlK45ZZbqvcxmUw1pqEBODg4sGXLFp566imWLl1KcXExQ4YMISEhgaCgoJa5EhERERERsXtNCjiurq4sWrSIRYsW1bvPihUraozoVOnUqRNvvPEGb7zxhtXnq1o++siRI00p02acPXuW1NRUo8sQg+k+EN0DAroPxEL3gYD93gdVmaCxV8yYzGazuS0Kao41a9YwZcoUo8sQEREREZF2YvXq1fzyl7+sd3u7Djh5eXls3bqVq6++Wu/KERERERGxY8XFxWRmZjJ69Gi6du1a737tOuCIiIiIiIg0hdXLRIuIiIiIiLR3CjgiIiIiImIzFHBERERERMRmKOCIiIiIiIjNUMDpgB588EEcHByIjIw0uhRpQzt27GDGjBkEBwfj6enJNddcw4MPPkhubq7RpUkrKC0tZc6cOfj7++Ph4cGwYcPYvn270WVJG9q/fz+PPPIIAwYMwMvLi169ejFp0iSOHj1qdGlioIULF+Lg4MDAgQONLkXaWGpqKlFRUfj5+eHp6cnAgQP561//anRZ7ZJWUetgPvvsM8LDw3FycuLOO+9k06ZNRpckbeSmm26ioKCAmJgYgoKC+Pbbb1m2bBkeHh58/vnnXHHFFUaXKC0oLi6Od955h/j4eIKCglixYgX79+/no48+Yvjw4UaXJ20gOjqaPXv2EBMTw6BBg8jJyWHZsmWcP3+eTz/9lAEDBhhdorSxEydO0K9fPxwcHOjduzeHDh0yuiRpI9u2bSMyMpLQ0FAmTZqEl5cXGRkZmM1mXnrpJaPLa3cUcDoQs9nM8OHDGTBgANu3b2fgwIEKOHbk3//+N7fcckuNtl27dnH77bfzzDPPsGDBAoMqk5a2b98+hg0bxiuvvMKsWbMAy4jO9ddfT/fu3dm9e7fBFUpb2LNnD2FhYTg5OVW3ZWRkMHDgQKKjo1m1apWB1YkRJk+eTH5+PuXl5eTl5fHll18aXZK0gXPnzhEcHMwtt9xCUlKS0eV0CJqi1oGsWrWKw4cP8+KLL6Jcan8uDTcAt956K126dCEtLc2AiqS1JCUl4eTkxMyZM6vbXF1d+dWvfsWePXv44YcfDKxO2srNN99cI9wA9O3bl+uuu05/5+3QJ598wjvvvMOSJUswm82YTCajS5I28vbbb3Pq1CkWLlwIwIULF6isrDS4qvZNAaeDKCwsZM6cOcybN09TkaTa+fPnKSwsbPBtvtLxHDx4kODgYLy8vGq0h4WFAfD5558bUZa0A2azmZMnT+rvvJ2pqKjg0Ucf5cEHH9TURDu0fft2fHx8OH78OP369cPb2xtfX19++9vfUlpaanR57ZICTgfxwgsv4OnpSXx8vNGlSDuyZMkSysrKmDRpktGlSAvKycmhZ8+etdqr2rKzs9u6JGkn1qxZQ3Z2tv7O25m///3vZGVlaSqynTp69Cjl5eWMHz+eMWPGsH79embMmMHf//53pk+fbnR57ZJT47tISzKbzVanbTc3NwDS09NZunQp//rXv3B2dm7N8qSNNOc+uNQnn3zC888/z6RJkxgxYkQLVidGKy4uxtXVtVZ71b1QXFzc1iVJO5CWlsbDDz9MeHg406ZNM7ocaSP5+fk899xzPPfcc/j5+Rldjhjg/PnzFBUV8Zvf/IYlS5YAMH78eC5evMg//vEPXnjhBfr27Wtwle2LRnDa2Mcff4yHh4dVn/T0dAAef/xxhg8fzoQJEwyuXlpKc+6Dn0tLS2PChAkMGjSI//f//p8BVyCtyd3dvc4AXFJSUr1d7Etubi5jx46lc+fOJCUl6fkLO/L73/+erl278uijjxpdihik6mt+XFxcjfaqP3/66Rf7ko8AAAOmSURBVKdtXlN7pxGcNta/f39Wrlxp1b49evTgww8/ZOvWraxfv57MzMzqbeXl5RQVFXHs2DG6dOmCt7d36xQsraKp98HPHT9+nLvuuovOnTuzZcsWPD09W6FCMVLPnj3rnIaWk5MDgL+/f1uXJAY6e/YsY8aM4dy5c+zatavW1wSxXUePHuWNN95gyZIlnDhxorq9pKSEixcvcuzYMXx8fOjcubOBVUpr8/f35/Dhw7Wewe7evTsAP/74oxFltWsKOG3siiuu4L777rN6/6ysLAAmTpxYa1t2dja9e/dmyZIlPPbYYy1Wo7S+pt4HVfLz87nrrrsoKyvjo48+0oITNurGG29k586d/7+9u2dpJQijOH78AJbpdw2Cha/BWpKAYBMUQYVgFYKNIAl2USsljSDoB1CElIpoZZGgiAQCAcVOEMGX2C0WUVBIzK304jW3dCeO/18525xmWc7s7jyqVqufNi9KpZIkqb+/31Q0+Ozl5UWxWExXV1fK5/Pq6uoyHQk+qlQqent709zcXNPnvOu6SqVSWltbM5AOfhkcHFQ+n9f9/b06Ozs/1t83wgKBgKloLYs5OC3u7u5OZ2dnn9YajYZmZmbkOI4WFhbU3d2tjo4OQwnhl+fnZ0WjUV1eXuro6EgDAwOmI+GbvM/BWV1d1fz8vKS/c3ACgYCKxaLhhPBDvV7X+Pi4Dg8Ptb+/r5GREdOR4DPP83R6evrpk8RGo6HFxUU9PT1pfX1dwWCQk9Usd35+rlAopHg8rlwu97Eej8e1u7urm5sb3uz+g4LzQzmOo97eXgZ9/iJjY2M6ODhQIpH4cqhAe3u7RkdHzQTDt5iamtLe3p7S6bSCwaC2t7dVLpdVKBSazkSCfVKplDY2NhSLxTQxMfHl+vT0tIFUaAXhcFie5zHo8xdJJpPa3NzU5OSkhoaGdHx8rJ2dHWUyGa2srJiO13IoOD+U67rq6emh4Pwiruvq9va26ZBXx3F0fX1tIBW+y+vrq5aWlpTL5fT4+Ki+vj4tLy9reHjYdDT4JBKJ6OTkpOk939bWpnq9biAVWkEkEpHnebq4uDAdBT6p1WrKZrPa2trSw8ODHMfR7Owsvyj8BwUHAAAAgDU4JhoAAACANSg4AAAAAKxBwQEAAABgDQoOAAAAAGtQcAAAAABYg4IDAAAAwBoUHAAAAADWoOAAAAAAsAYFBwAAAIA1KDgAAAAArEHBAQAAAGANCg4AAAAAa1BwAAAAAFjjD8pR/tt+TMPDAAAAAElFTkSuQmCC", + "text/plain": [ + "PyPlot.Figure(PyObject )" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "(-0.1,3.6)" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Converged. p = [2.110067133637596,1.2802355090457196,1.5675520985912952,-3.7904167887303712]\n" + ] + } + ], + "source": [ + "# calculate and average normals like we did in last notebook\n", "average_normals!(Γ₁)\n", "average_normals!(Γ₂)\n", "\n", @@ -431,9 +516,7 @@ "outputs": [ { "data": { - "image/png": [ - "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4VPX99vF7kpkJ2diCyCabCQEFRVEQiaJAA/4gVFxRqiiVAlYfRagookVNJbiLVbAUFLHaYpFqUmRfg7K4pEgChIBhlS0BsyeTmXn+GBMJWckyZ2byfl1XLug5M+fcwUJy53vO55icTqdTAAAAAOAD/IwOAAAAAAD1hYIDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAADAg+3evVsREREKCAiQ1WpVQECAIiIitHv3bqOjAR7J5HQ6nUaHAAAAQHnJycnq3bu3iouLy+0zm83auXOnevToYUAywHNRcAAAADxURESE0tLSKt0fHh6uffv2uTER4PkoOAAAAB4qICBARUVFle63Wq0qLCx0YyLA83EPDgAAgIeq7ufQ/JwaKI+CAwAA4KFMJlOd9gONEQUHAADAAxUUSP7+Hat8TceOVe8HGiMKDgAAgIfJz5dGjpQcji/k72+u9HUPPfSQG1MB3oGCAwAA4EHy8qSYGGnLFunLL3vohx92Kjw8XFarVRaLRVarVZdeeqluu+02TZ8+XZ988onRkQGPUvmPBAAAAOBWubmucrN9u7R8uTRwoCT1qHAUtMPh0Lhx4/S73/1O/v7+uuuuu9yeF/BEFBwAAAAPkJMjDR8uffedtGKFFBVV9ev9/Py0YMEC2e123XvvvfL399ftt9/unrCAB+M5OAAAAAbLzpZuuUXaudNVbq6/vubvtdvtuu+++/Tpp59qyZIlGjVqVMMFBbwABQcAAMBAP//sKjfJydLKldJ11134MYqLizVmzBh99tlnWrp0qUaOHFn/QQEvQcEBAAAwyNmz0rBh0t690qpV0rXX1v5YNptN99xzj7744gstW7ZMw4cPr7+ggBeh4AAAABjgzBkpOlrav19avVrq06fux7TZbLrrrru0fPlyff755xo2bFjdDwp4GQoOAACAm2VkSL/5jXTokLRmjdS7d/0du6ioSHfccYdWrVqlL774QtHR0fV3cMALUHAAAADc6PRpacgQ6ehRae1a6Yor6v8chYWFuv3227V27VolJCRo8ODB9X8SwENRcAAAANzk1Clp8GDp+HFp3TqpZ8+GO1dBQYFGjRqljRs3asOGDerbt2/DnQzwIH5GBwAAAGgMTpyQbr5ZOnlS2rChYcuNJDVp0kTLli3T5MmTdfnllzfsyQAPwgoOAABAAzt+XBo0yDU1bd06qXt3oxMBvstsdAAAAABfduyYq9xkZ7tWbrp1MzoR4NsoOAAAAA3k6FHXZWn5+dLGjVJ4uNGJAN9HwQEAAGgAhw+7yo3N5io3XbsanQhoHCg4AAAA9ezgQVe5cTpdl6V16WJ0IqDxYIoaAABAPfrxR2ngQMlkotwARqDgAAAA1JMDB6SbbpLMZle56dTJ6ERA40PBAQAAqAdpaa6Vm4AA1z03l1xidCKgcaLgAAAA1FFqqqvcBAe7yk379kYnAhovCg4AAEAd7NnjKjfNmrkuS2vb1uhEQONGwQEAAKillBTXPTdhYdL69VKbNkYnAkDBAQAAqIVdu1yjoFu3dpWbiy82OhEAiYIDAABwwXbudJWbdu2kdeukiy4yOhGAEhQcAACAC5CU5Co3HTtKa9dKrVoZnQjAuSg4AAAANfTdd9KgQVLXrtKaNVLLlkYnAnA+Cg4AAEAN7NghDR4sRURIq1dLLVoYnQhARSg4AAAA1di2TfrNb6QePaRVq6TmzY1OBKAyFBwAAIAqfP21q9z07CmtXOl63g0Az0XBAQAAqMSWLVJ0tHTVVdKKFVJoqNGJAFSHggMAAFCBTZukoUOla6+Vli+XQkKMTgSgJig4AAAA51m/XrrlFum666SEBCk42OhEAGqKggMAAHCOtWul4cOlAQOk+HgpKMjoRAAuBAUHAADgF6tWSSNGSAMHSl98IQUGGp0IwIWi4AAAAEj68ktp5EjXs27+8x+pSROjEwGoDQoOAABo9BISpFtvdU1MW7pUCggwOhGA2qLgAACARu3zz6XbbnPdd/Pvf1NuAG9HwQEAAI3WsmXSHXe4Lk37178kq9XoRADqioIDAAAapX//W7rrLtfqzSefSBaL0YkA1AcKDgAAaHT+9S9p9Gjpzjulf/yDcgP4EgoOAABoVD7+WLr3XtfH4sWS2Wx0IgD1iYIDAAAajQ8/lO67z/Xx/vuSv7/RiQDUNwoOAABoFD74QHrgAdfHwoWUG8BXUXAAAIDPW7BAGjdOGj9emj9f8uM7IMBn8dcbAAD4tL/9TXroIWniRGnuXMoN4Ov4Kw4AAHzW3LnShAnSo49K77xDuQEaA/6aAwAAn/TXv0oPPyw99pj01luSyWR0IgDuQMEBAAA+5623XKs2U6ZIb7xBuQEaEwoOAADwKa+9Jj3+uDRtmvTKK5QboLGh4AAAAJ/x8svS1KnSM89Is2ZRboDGiIIDAAB8wksvuVZtnntOevFFyg3QWFFwAACA13vxRdeqzfPPuz4oN0DjRcEBAABeISEhQRaLRSaTqfTDYrFo9OgEPfecFBvrWr0B0LiZnE6n0+gQAAAAVfn888916623Vrr//vvjtWjRCDcmAuCpKDgAAMDjWSwWFRcXV7rfbDbLZrO5MREAT0XBAQAAHs9Ug5tq+JYGgMQ9OAAAAAB8CAUHAAAAgM+g4AAAAI/1c8HPemrNU1KLql9nNpvdEwiAx6PgAAAAj1PsKNbcHXMV/na43t7+tu557p4qX79s2TI3JQPg6Sg4AADAYzidTi3ft1xXzL1Cf1z+Rw2PGK7UR1L18eMfKz4+vtxKTcnwgbVr1xoRF4AHYooaAADwCDtP7NTUVVO1+sBq3dT5Jr0W/Zqubnt1le9xOBzq3LmzDh8+rDVr1mjw4MFuSgvAU1FwAACAoY7nHNez657Vgu8XKLxluF6NflUx3WJqNBpako4cOaKuXbvKbDbr2LFjat68eQMnBuDJuEQNAAAYIs+Wp9hNsQqfE66lu5fqzWFvatfDuzQycmSNy40kdejQQYsXL1Z+fr4GDBjQgIkBeANWcAAAgFs5nA59/MPHmr52uo7nHNejfR/VjBtnqEVgNaPSqjFmzBh9/PHHevTRRzVnzpx6SgvA21BwAACA22w+uFlTVk3RjmM7dFuP2zR7yGyFtwyvl2M7HA517dpVBw8e1IoVKzR06NB6OS4A70LBAQAADW5/5n5NWzNNS3cvVZ+2ffTG0Dd0Q6cb6v08x44dU+fOneXv76+jR4+qZcuW9X4OAJ6Ne3AAAECDOZN/RlNWTlGPd3po29Ft+vDWD7V9/PYGKTeS1K5dO3388ccqKCjgfhygkWIFBwAA1Dub3aa538zV8xufV2FxoZ6KekpP9H9CQZYgt5z//vvv1+LFi/Xwww/rnXfeccs5AXgGCg4AAKg3TqdT8anx+tPqPyktM03jeo/Ti4NeVJuQNm7N4XA4dOmllyo9PV3Lly/XLbfc4tbzAzAOBQcAANSLpONJemLlE1qfvl6DuwzWa9Gv6co2VxqW5/jx4+rUqZP8/Px0+PBhtWrVyrAsANyHe3AAAECdHMs+pnGfj9PV712tn3J+UsI9CVp932pDy40ktWnThvtxgEaIFRwAAFAruUW5evWrV/XyVy8ryBKk5296XuOvHi+Lv8XoaGU8+OCD2rhxo1JSUtSkSROj4wBoYBQcAABwQRxOhz7834d6Zt0zOp13Wo/1e0zTb5iu5k2aGx2tQg6HQ5Lk58eFK0BjQMEBAAA1tiF9g55Y+YS+P/697rr8LsUNjlOXFl2MjgUApcxGBwAAAJ4vNSNVT65+Up/v/Vz92vfTlnFbdP0l1xsdCwDKoeAAAIBKZeRl6IWNL+jdb95Vu9B2+uT2T3T35XfLZDIZHQ0AKsQlagAANEKLFi3SwIED1blz53L70tPTtXb9WmV1z9ILm16Q3WHX9Bum67F+jynQEuj+sABwASg4AAA0Qunp6Ro3bpwWLlxYpuT8+OOPGnnPSGUPy9Zhv8N66KqH9MLNL+jikIuNCwsAF4BxIgAANEKdO3fWwoULNW7cOKWnp0uS4rfF66r/u0q7rtul7uHd9b+J/9N7Me9RbgB4FVZwAABoxNLT03XH3XeoqGORfvjmB0WMi9Db97ytoeFDjY4GALXCkAEAABqpnKIcLTiwQEmOJNn/bdfvY3+veU/Pk9mPbw8AeC9WcAAAaGTsDrs+SPpAM9bPUOaxTLVZ00Y5R3JksVi0devWCgcPAIC34B4cAAAakTUH1ujqv12th+IfUt/gvuqzo482fr5RL730kk6cOKG77rqr9J4cAPBGFBwAABqB3ad2a8THI/Sbxb9RiDVEn0V/puwl2fr4w4/VuXNnjRs3TuHh4QoNDS0zeAAAvA0FBwAAH3Yq95QeWf6Ies3tpZRTKfr0zk+V+GCislKzyoyItlgsio2N1bp16zRp0iRt3LjR2OAAUEvcgwMAgA8qLC7UnG1z9JfNf5FTTs24YYYe7feompibVPoeh8Oha665RkFBQdq8ebNMJpMbEwNA/aDgAADgQ5xOp/6d8m9NWzNNh34+pInXTNSfB/5ZFwVfVKP3r1y5UsOGDVNCQoKGDx/ewGkBoP5RcAAA8BHbjmzTE6ue0FeHv9LwiOF6NfpVdW/V/YKO4XQ6NWjQIGVkZCgpKUl+flzNDsC78K8W4IF2796tiIgIBQQEyGq1KiAgQBEREdq9e7fR0QB4oINnD+repffqugXXKacoR6vvW62EexMuuNxIkslk0qxZs/TDDz/ok08+aYC0ANCwWMEBPExycrJ69+6t4uLicvvMZrN27typHj16GJAMgKfJKszSrM2z9MbWN9QisIX+MugvGnvlWPn7+df52Lfeeqt27typPXv2yGq11kNaAHAPCg7gYSIiIpSWllbp/vDwcO3bt8+NiQB4mmJHsRZ+v1DPrn9W2YXZmtJ/ip4c8KRCA0Lr7RzJycnq1auX3n77bf3xj3+st+MCQEOj4AAeJiAgQEVFRZXut1qtKiwsdGMiAJ5kZdpKTV09VbtO7tJ9V9ynlwa/pA5NOzTIuR544AGtWLFC+/fvV3BwcIOcAwDqG/fgAB6mup858DMJoHHadXKXhn00TMP+MUwtmrTQjvE79OGoDxus3EjSzJkzdebMGb311lsNdg4AqG8UHMDDVPfcCafNqe+ivtPBuIPK2ZVD4QF83Mnck5qYMFFXzrtS+8/s12d3faaND2zUNe2uafBzd+7cWZMmTdLs2bOVkZHR4OcDgPpAwQE8TMeOHavc3+GiDrJeZNXBFw/qm17faGuXrUp9JFUZKzJkL7C7KSWA+hAZGamJEyfKYrHIZDKVflgsFo3/w3hd1PEihc8J17+S/6VXf/Oqkh9O1qgeo9z6AM7p06fL4XBo9uzZbjsnANQF9+AAHmb37t3q1auX7PbyZeXcKWr2ArvObjirjIQMZf43UwXpBfIL8lOLIS0UNiJMYcPDFNAuwIDPAEBN/eEPf9D8+fMrf8E10mMvPqZnb3xWYUFh7gt2npkzZ2r27Nnat2+fOnRouEviAKA+UHAAD/TGG2/oiSeekMVikeS6bK1jx4764osvKhwR7XQ6lZeSp4z/ZigjPkM/f/Wz5JBC+oS4ys6IMIVeHSqTn/t+6gugehaLpcKR8CX8zf4qtlW+312ysrJ06aWXatSoUfrb3/5mdBwAqBIFB/BAY8eO1XfffacffvihVu+3ZdqU+WWma3VnRaaKzxbLcrFFYcNdZafFkBYyh5rrOTWAC1WTS8085cv0m2++qalTpyo5OVmRkZFGxwGASlFwAA9jt9vVpk0bPfTQQ5o1a1adj+ewOZT1VZYyEjKUkZChvD15MllNan5T89LVncAugfWQHEC1MjOlTZukDRukDRtk/d//ZKvmLZ7yZbqgoECRkZHq16+flixZYnQcAKgUBQfwMFu3blX//v2VmJioAQMG1Pvx8/fnu8rOfzN0dsNZOW1OBV0WVFp2mvZvKj8z80eAenFeodHOnZLTKXXpIg0cKP8PPpCjmkN40pfp999/X+PGjdOOHTt0zTUNP8UNAGqDggN4mBkzZmju3Lk6efKk/P39G/RcxdnFOrPqjOvenf9myHbSJnMLs1re0lJhI8LUcmhLWVpaGjQD4FMqKzSdO0s33eT6GDjQ9b9V/T04ZrNZNlt1azzuU1xcrCuuuEIdOnTQqlWrjI4DABWi4AAepnfv3urVq5cWL17s1vM6HU5lf5OtjHjXpWw5STmSv9Ts+mYKi3Gt7gR1D3LreFrA41VXaAYOlG6+WerUqcK3T5w4Ue+9916lh58wYYLmzZvXEMlrbdmyZbrtttu0Zs0aDR482Og4AFAOBQfwIIcPH1bHjh31z3/+U3fffbehWQqPFpau7JxZfUaOfIeadG1Seilb8xubyy+AS9nQyFRzyZluvtn1ayWF5nyRkZG6+eabtWDBgnIrOQ888IC++uor7d27t/4/jzpwOp267rrrJLkuqeWHHgA8DQUH8CDz5s3TI488otOnT6t58+ZGxyllz7fr7LqzrsKTkKHCw4XyD/FXi9+4nrnT8v9aKqANz9yBD7rAS87qQ9++fbVjxw6Fh4dr37599Xbc+rR+/XoNGjRIn332mUaNGmV0HAAog4IDeJARI0YoNzdX69evNzpKpZxOp3J/yC2dypa1NUtySqHXhrpWd2LCFNI7hJ/qwjvV8ZKz+lBQUKDAQNdkQ0/+Ej106FAdPnxYO3fulNnM2HkAnoOCA3iIvLw8hYWFKTY2VlOmTDE6To0VnSr69Zk7KzNlz7LL2s766zN3BreQf3DDDksAzjV+/HiNHTtWUVFR5fYlJiZq0aJFmj9/vmtDPV9yVl+aNGmiwsJCvfHGG3r88cfdeu6a+vbbb3XNNddo4cKFevDBB42OAwClKDiAh0hISFBMTIz27NnjtQ/RcxQ59HPiz6WrO/n78mUKMKnFINelbGHDw9SkUxOjY8LHJSYmKiYmRvHx8WVKTmJiomJGjFD8k08q6uRJt11yVhuLFy/W/fffL4vFoqKiIkOzVOXuu+/W119/rdTUVDVpwt9tAJ6BggN4iIkTJ2rt2rVKTU31mcu78lLzSp+58/Omn+Usdiq4V/Cvz9zp11Qmf9/4XOFZSkrO3FdeUdiRIwpMTlbMsmWKt9sVJbn1krPaKvl3IDc3V0FBQQanqVhqaqouu+wyvfLKK5o8ebLRcQBAEgUH8AhOp1OXXHKJ7rzzTr3xxhtGx2kQxT8XK3PVL5eyLc+U7bRN5jCzwv7vl0vZolvI0pxn7qCOzrnkLDEhQTfu369ASVY/P8VHRyvqnnsMueSsNq666iolJSWpR48eSklJMTpOpSZMmKClS5fqwIEDatq0qdFxAICCA3iCpKQkXXXVVY3muRJOu1NZ27NKn7mT+0OuTGaTmkU1K13dCYr0zJ9Yw8NUNRTg5pvVY8UK7fnpJ7388sv605/+ZHDYC+MtwwaOHj2q8PBwPfnkk3r++eeNjgMAFBzAE8TGxurll1/W6dOnZbVajY7jdgWHCkqfuXN27Vk5ChwKjAgsLTvNoprJz8ozd6ALmnKWmJiooUOHKi8vT0FBQVq5cmWFgwc8WUBAgIqKijRv3jxNmDDB6DiVevLJJ/Xuu+/qwIEDat26tdFxADRyFBzAA3z22Wfas2ePpk+fbnQUw9nz7Dqz5kzpM3eKjhXJP9RfLYe2dD1z55aWsrYuWwIXLVqkgQMHqnMFN4anp6dr48aNGjt2rJs+A9SrWk45K7kHZ8mSJYqOjtZ1112nPXv2lBs84On+/ve/a/z48bJarSosLDQ6TqUyMzPVtWtXjR07Vm+99ZbRcQA0chQcAB7L6XQqJymndCpb9o5sSVLTfk1Ln7kT3CtYBw8e1Lhx47Rw4cIyJSc9Pb3C7fBglRWaTp1+LTM33VTllLPzp6iFhYXJbDZr6dKlFU5X83TeMGxAkl566SXNnDlTqamp/H0DYCgKDgCvUXSiSBnLXZeynVl5RvYcuwIuCVDY8DDlXpurJz58Qu9/8L46d+6s0aNHa8eOHTp48KDsdnvpMcxms2JjY5WWlvbrs1AayAU9j6WxqsklZ9UUmvOd/+c+ZMgQrV27VoWFhdq+fbvX/bn37NlTycnJ6tWrl3bu3Gl0nErl5ubq0ksv1dChQ7Vo0SKj4wBoxCg4ALySo9Chs5vOugYV/DdDBQcKdCLghKb4TdHvon+nJd8v0Y+Hfqz0/XFxcZo2bVqDZqzyeSwevJLQoMXMgAdr/vWvf9Wjjz6qjz76SGPGjKm347pLXl6egoODJXn2sAFJevfdd/XII49o586d6tmzp9FxADRSFBwAXs/pdCpvb55Of3Fal02/THn2PIUpTBnKqPQ9ZrNZNputwbOdX2Y8vdxI1RezqKgoRUVFacaMGSouLi7dX+HqWDVTzjRwYIM/WDMrK0vNmjXTrbfeqmXLljXYeRpSybCB999/Xw888IDRcSpVVFSkHj16qFevXvrPf/5jdBwAjRQFB4DP2L59u/r166eZ02dq5kszq319r+3bZTGZZPXzk8VkKv299fzfn7Ot5Pcl77Ge97qSbee+bs/WrZo+erSihg7VV6tW6a+ffqr+UVGVns/fZDL8Ya9VFbNNmzbpmWeeqfS9cb/7naaFhVU75cydmjdvrqCgIB07dsyt560vc+bM0WOPPabWrVvrxIkTRsep0scff6wxY8boq6++Uv/+/Y2OA6ARouAA8Anp6ekaNWqUkpKStHfvXkVGRlb7nkdTU2VzOlXkcJT+WuR0lttmczpVVMk22y/vKXI4ZK/qZI895vpm//77pQcfrDbbueWpoiJUrpBVsP9CStj557OYTNq7bZueGT1a/YYM0fa1a0uL2eUtW6q4ioleZkm2BrzkrDYGDhyozZs3q6ioSGaz2dAsteXv7y+HwyG73S4/P88dm+5wOHTVVVepefPm2rBhg+FlHUDj453/ygPAOUqmpf32t79VSkqKunbtWqP3zYmIqNccjl/Kke28ovTVli26b/du2SSFfP65Xhk9Wpf171/j8lRyzMIKjn3+67LtdtmKiy+orFVazIKCpJ49tW7ZMunWW/WA1Spt3y7Zq6xyKpakAwfq9c+2rkaOHKlNmzbp888/1+233250nFoZPHiwVq9erfvvv18fffSR0XEq5efnp1mzZmn48OFauXKlhg0bZnQkAI0MKzgAvF7Jc3Di4uL09ddf63//+58sFkuZ+0PO5+57cHr27KnExERt3rzZ4+7BcTidKv6lDBWeU3y2bNmicSNHKj8vTwGBgXpz6VJd1r+/BrZoUe0xPe1Ly+nTp3XRRRfpjjvu0Keffmp0nFopGTbgrv/v1oXT6dTAgQOVnZ2tb7/91qNXnAD4Hv7FAeD1xo4dq86dOyslJUWXXXaZJCk2NrbK91S3vz6ce9+K2WyWn5+foqKiFB8fr5iYGCUmJjZ4hprw++XStGB/f7W0WHSx1aqD33yjSXfcoXlz50qS7rjtNj19773y27XL4LS106pVK4WGhmrr1q1GR6m1oKAgtWjRQsXFxR7z/53KmEwmzZo1S0lJSVqyZInRcQA0MhQcAD7B6XQqOTm5tOCkpaUpLi6u3PX/ZrNZcXFxSktLa/BMixYtKl2pyc7Olr+/vySVlhxPfVbIucVs9OjRklx/viXFrLqfxnvqPS49e/bU0aNH5XA4jI5Sa2+//bYk6d577zU4SfUGDBigESNGaMaMGR6/4gTAt1BwAPiEkydPKjMzU5dffrkkaf78+Zo2bZpiYmIkSS+88IKcTqdsNpumTZvmlgc9zp8/v/QytNzcXFksltJ9UVFRHvuwyXOLmdVqlZ+fn3766afSYtanT58q3++O1bHaGD58uJxOp1asWGF0lFobM2aM/Pz8dPjwYa8oan/5y1904MABLViwwOgoABoRCg4An5CSkiJJpSs4JZKSkiRJkyZNcnumc+Xm5spqtRqaoabOLWaSZLVaderUKUmuYnbllVcqLi6u3EqNO1fHauPBX6bX/eMf/zA4Sd3ceOONklwPZPV0V1xxhcaMGaPnn39eeXl5RscB0EgwZACAT3jnnXc0efJk5eXllfnGOzAwUAUFBYbf9N66dWv5+/vrp59+MjRHbYSFhSkwMFBHjhypcF9mZqbhf741FRISolatWik9Pd3oKLVW8uBSi8WioqIio+NU68CBA+revbteeOEFPfXUU0bHAdAIsIIDwCckJycrMjKy3KpCYWGhR0xwKiwsVFBQkNExaqVZs2bKzs6ucF+bNm0kuf78vUGPHj285vKuyjRt2lTNmzeXzWbTjh07jI5Tra5du2rChAmaPXu2zpw5Y3QcAI2A8V/1AaAenDtB7VxOp9MjikVRUZFH5KiNli1bKj8/v8J93bp1kyR99tln7oxUa7fccoscDoc2bNhgdJQ6ee211yRJd911l8FJaqZk0MDs2bONjgKgEaDgAPAJycnJpQMGzt0mSe3btzciUhnFxcUKCQkxOkatXHzxxZVOwSq5V+frr792Z6Ra+/3vfy9JWrx4scFJ6mbcuHEymUxKT0/3itWoiy++WJMnT9acOXN07Ngxo+MA8HEUHABe79SpUzp9+nS5FZx58+ZJkvr3729ErDLsdrtCQ0ONjlErJQXx7Nmz5fbdfffdkuSxgwXO16lTJwUGBmrTpk1GR6mz66+/XpL0xz/+0eAkNTN16lQFBgbqxRdfNDoKAB9HwQHg9SqboLZx40ZJ0kMPPeT2TOdzOp1q1qyZ0TFqpVOnTpIqvs+mQ4cOklQ6Zc0bREZG6uDBg0bHqLOEhARJ0sKFCw1OUjPNmjXT9OnTNX/+fO3bt8/oOAB8GAUHgNe0b8NdAAAXaklEQVRLTk6W2WxWREREme0l38QOGDDAiFilsrKyJEktWrQwNEdthYeHS5JSU1MrfU1ubq674tRZdHS07Ha7tmzZYnSUOmnevLmaNm2qoqIi7dy50+g4NfLwww+rbdu2eu6554yOAsCHUXAAeL2UlBR169atzIM0Jdc33SaTyaBUvyoZDe2tBScyMlKSa9xvRfz8/FRcXOzOSHUybtw4Sd5/H44kxcXFSZJGjRplcJKaCQwM1MyZM/XPf/5T33//vdFxAPgoCg4Ar1fZBDW73e4RD9csKTitWrUyOEntdO/eXZJ06NChCvdbrVaveQ6O5CpsAQEBWr9+vdFR6mzSpEkymUw6cOCAVwwbkKSxY8cqMjJS06dPNzoKAB9FwQHg9SqaoJaTkyPJ9SBKo504cUKS62Gf3qhJkyYymUyVTr/yxulwERER+vHHH42OUS/69u0rSZo8ebLBSWrGbDYrNjZWK1as8Ppx3QA8EwUHgFc7ffq0Tp48WW4FZ8GCBZKknj17GhGrjJIb8C+66CKDk9Se1WrVyZMnK9xX8rDP7777zp2R6mTIkCGy2WxelbkyJcMG3nvvPYOT1Nztt9+uPn366Omnn/aq1T8A3oGCA8Cr7d69W1L5CWr/+c9/JLm+kTLa6dOnJf1aBLxRUFBQpU+hL7mEreTP3Bs8+OCDkqQPPvjA2CD1oFWrVgoJCVFhYWHpREFPZzKZFBcXp61bt+qLL74wOg4AH0PBAeDVkpOT5e/vr27dupXZXrJqcv/99xsRq4zMzExJUrt27QxOUntNmzYtnQZ3vptuukmStG3bNjcmqpsrrrhCVqtVa9euNTpKvSh5toy3DBuQXKtogwcP1vTp02W3242OA8CHUHAAeLWUlBRFRESUGyZgt9sVGBioJk2aGJTsVyUFx1vvwZFc9zLl5+dXuK9klWz//v3ujFRnXbt29ZoHlFbn8ccfl8lkqnKUtyeaNWuWUlJS9NFHHxkdBYAPoeAA8GqVTVA7duyYxxSKs2fPSnKNU/ZWrVu3ls1mq3BfyaV33jQqWpJuvvlmFRUVVfgAU2909dVXS5KmTZtmcJKau/baa3X77bfrueeeU2FhodFxAPgI7/1qCwCqeIKaJGVnZ6tr164GJCovOzvbq8uNJLVv315Op7N0Ot35goODFRAQ4OZUdeNL9+FIvw4bmDNnjsFJLkxsbKyOHDniVUMSAHg27/6KC6BRy8zM1PHjx8ut4Ozbt09Op1NXXnmlQcnKys7OltlsNjpGnXTs2FHSr0Mdzte0aVNlZGS4M1KdjB8/XuvWrZMkvfrqqzKZTDKZTLJYLJo9e7bGjx9vcMIL16ZNGwUHB6ugoMCrLlXr3r27HnzwQcXGxio7O9voOAB8AAUHgNcq+Wb7/BWckm9cBwwY4PZMFcnJyZHFYjE6Rp1ceumlkqS9e/dWuL9169aVru54oi5duuipp54qt724uFhPPfWUwsPDDUhVd88++6wk6dZbbzU4yYX585//rKysLL3xxhtGRwHgAyg4ALxWcnKy/Pz8yk1Q27p1qyRp0KBBRsQqJz8/3+su3zpfZGSkJFV6U36HDh286h6KP//5z1XunzFjhpuS1J/IyMjSh5fu3r27zKrUxIkTS/8beqJLLrlEjzzyiF599dXSCYgAUFvefc0EgEYtJSVF4eHh5cpDcnKyzGazWrZsaVCysvLz8z1imltdlFwGePDgwQr3l9zvdOjQodLL2TxJfr508KCUnu76KC6u+uGSxcXFuu222xQcHKzg4GAFBQWV+bWmv/f393fL5ydJAwcOrPA+luLiYr333nuaMGGC27LUxtNPP6358+dr1qxZev31142OA8CLUXAAeK3KBgwcOnRILVq0MCBRxYqKijwqT20EBQXJZDLp2LFjFe7v0aOHJOn77783pOCcX2DO/zhx4tfXujqHqQbHzFdGRoZyc3NLP/Ly8pSbm6uCgoIa5bJarVWWpNqUpnOPc+69Xe+//36VWRYsWKB58+bVKLcRwsLCNHXqVMXGxurxxx/3yKIMwDtQcAB4rZSUlNJJWOfKzMxUr169DEhUsaKiIgUFBRkdo84sFotOnjxZZtv48eM1duxY9e7dW5L0ww8/6Le//a0kKTExUYsWLdL8+fPrfO4LLTCXXCJ17iz16CHdcovr9yUf7dtLFkv1I62//PLLSvfZ7Xbl5eWVFp5zy0/J7yvad+623NzceilQJYWnujHd3jDGe/LkyXr77bc1c+ZMLVy40Og4ALwUBQeAVzp79qyOHTtWboJaVlaWbDZb6YqCJyguLlZoaKjRMWpl/PjxCg8P14wZM1RcXKykpCSZTCaZzWbFxsbq+PHjiomJ0ZIlSySpdHpXYmKiYmJiFB8fX6Pz1LbAdO8uDRvm+n2XLr8WmOqG1pnN5iq/4a9u6p2/v79CQ0Mb7L+r3W5Xfn5+uUJUVaF64YUXqj3unj0PKjS0r5o27avg4F7y87NW+x53CgkJ0bPPPqvHH39cU6dOrfAZVwBQHQoOAK+UkpIiqfwEtQ0bNkiS+vXr5+5IlXI4HF5bcKqbNhYXF6dp06YpJiZGJpNJBw8eLFNuoqKiJNV9BaZTJ1eB6dJFateu+gJTndjY2Ao/r3P3G8nf318hISEKCQmp8XtqUnBycnbqxImP5HQWy2QKUGjoVaWFJzS0rwIDw2UyVX/5XkP6wx/+oNdff10zZszQZ599ZmgWAN7J5HQ6q77TEgA80N///ndNmDBBOTk5CgwMLN3+5JNP6pVXXtGuXbsqvD/H3RwOh/z9/TV69Gh98sknRse5YBaLpdqVDpvNprVrEzVkyA1q2rSVbLZijRoVL7s9qtoC06WLq7x07vzrrx061L3AVOf8lalzP5/Y2FilpaXVy6V17lTT/1Z2e75ycpKUnb1dWVnblZ29Xfn5ab+8pvk5hedahYb2VUBAG3d9CqUWL16s+++/X1u3bvWoH1YA8A4UHABe6YknnlB8fLz27dtXZnt0dLTWrFmj4uJi+fkZPwk/MzNTYWFhmjRpkt59912j41ywmvw0/+KLnb8UmOaSfpbJ9KY6dnys9JKx8z9qcgkZLtzEiRMrnKJWYsKECZUOGbDZMpWdvaO08GRlbZfN5rrfKiDgkjKrPKGhfWQ2N+yKpN1uV+/evdWqVSutW7fO8FUlAN6FLzEAvFJlE9T279+v4OBgjyg3kkqnjnn7FLWqTJokFRYm6s037br33slaunSmPvqoT+nlaXCP9evXa8KECVqwYEG5Vanf//73Wr9+faXvtVhaqmXLoWrZcqgkyel0qrDwcJnCk57+ghyOXEkmBQX1+KXw9FPTptfW+/08/v7+eumllzRy5EitXr1a0dHR9XZsAL6PggPAK6WkpOi+++4rt/3EiRNq166dAYkq9tNPP0mSLrroIoOT1IzNbtOWw1uUkJqg+NR4yV+Sver3DB7suudm1aovFRUVpQceuK3cPThoeHv37pWkehkFbTKZ1KRJRzVp0lGtW98hSXI67crNTSmz0tOQ9/OMGDFC119/vaZPn64hQ4Z4zA8tAHg+Cg4Ar5OVlaUjR46Um7DkcDiUm5ur8PBwg5KVVzJW2ZMLTmZ+plakrVB8arxWpK3Q2YKzahPSRiMiRijNkiaH3VHpe/38/MqVmaioKMXHx1NyfIzJ5K+QkF4KCemltm3HSdIv9/N8X1p4MjOX6+jROZIks7mFQkOvLS08TZv2ldV68QWcz6S4uDjdeOONWrp0qe68884G+bwA+B4KDgCvU9kEtaSkJEnSVVdd5fZMlTl16pQkqXXr1gYn+ZXT6dSe03uUkJqghH0J2nJoi+xOu65ue7X+X9//pxHdRqhPuz7yM/kpfGZ4ldPG+vTpo9dff71ciSkpOYsWLaLg+DB//0A1a3a9mjW7vnTb+ffzHDv2N9lsrql0F3o/z4EDB9S3b1+NGTNGY8aMkeQqPh07dtTcuXN19OhRjR07tmE/SQBeh4IDwOukpKTIZDIpMjKyzPaSewxuvPFGI2JVKCMjQ5LUtm1bQ3MU2Yu0+eDm0kvP9p/Zr0BzoAZ3Hax3h7+r4RHD1b5p+3LvS0tLU1xcXJXTxiorMFFRUZSbRqjy+3m2VXI/z2VlVnlc9/NYJLn+3mzfvr3cOdLS0hQdHa1Vq1a581MD4CUoOAC8TkpKirp06aKgoKAy23fs2CFJuuGGG4yIVaGSgtOmjftH7Z7OO60v932p+NR4rdy/UlmFWWof2l4juo3Qm93e1KAugxRkCaryGCWjkqdNm+aOyPBBZe/ncV1m9uv9PNuVlbVD2dnbdPz4h5Ls8vNropCQ3goN7asJE/5Z6XGdTqcmTZpUbpIiAFBwAHidyiao7d69W1artVzxMdKZM2ckSS1btmzwczmdTqWcSlF8arwSUhP09ZGv5XA6dG27azWl/xTFdItR7za9GbkLw5W9n+f3kiS7Pe+X+3l2KDt7hzIy/qsjR05WeZxDhw65Iy4AL0PBAeB1UlJSdM8995TbfuTIEYWFhRmQqHI///yzTCZTg02AKiwu1MaDGxW/N14J+xKUfjZdwZZgDek6RH8b8TcN7zZcbULcv3oEXCh//yA1azZAzZoNKN1mMlkl2Sp9D4/yA1ARCg4Ar5Kdna1Dhw5VuILz888/q2/fvgakqlxWVla9l5uTuSe1fN9yxafGa9X+VcopylHHZh01ImKERnQboZu73Kwm5ib1ek7ACNWtNrIaCaAiFBwAXmX37t2SVG5E9PHjx2W329WzZ08jYlUqOztbZnPd/ql1Op364eQPit8br/jUeG0/6rrpul+HfnpqwFOKiYxRr9a9+GYPPqdjx45KS0urcj8AnI+CA8CrlIyI7t69e5nta9eulST179/f7ZmqkpeXJ6v1wp/wXlBcoPU/ri+9n+Zw1mGFWEM09NKhmnTNJN0ScYtaB3vO6GmgIcydO1fR0dEVXopmMpk0d+5cA1IB8HQUHABeJTk5WV26dFFwcHCZ7V999ZUkafDgwUbEqtSFFJyfsn8qvfRs9YHVyrPlqUvzLrq1+62K6RajGzvdqABzQAMnBjzH0aNHtWrVKk2aNEmHDh2S0+ks9xwcADifyckdegC8yPDhw2UymZSQkFBm+4033qgtW7bIbrcblOxX48ePV3h4eJXPjpk/f76cTqe+P/596bNpvjn2jfxMfurfob9GdBuhmG4xuuyiy7j0DACAC0DBAeBVunTpojvvvFMvv/xyme2XXHKJcnJySscyG+mll17SM888U+n+sVPHKuDGACXsS9Cx7GNqGtBUw8KHaUTECN0ScYtaBbVyY1oAAHwLBQeA18jJyVFoaKg++OADjR07tsy+Jk2aKDw8XLt27TIo3a8sFkuZlZtyTFL4W+EaETFCMZExuqHjDbL4W9wXEAAAH8Y9OAC8xp49eySVn6BWVFSkwsJCRUZGGhGrnCrLjSQ5pdRHUrn0DACABtAwT54DgAZQMkGtR48eZbYnJiZKkq655hq3ZyrH4VBNagvlBgCAhkHBAeA1kpOT1alTJ4WEhJTZvmnTJknSTTfdZECqX9jt0j//KV1xhbjYDAAA41BwAHiNlJSUcpenSdK3334rSbr22mvdHUkqLpY++kjq2VO65x6pQwdZ/P2rfEtdH/wJAAAqR8EB4DUqKzj79u1TYGCge4uDzSa9/77Uo4d0331SRIS0bZu0YoWe/ctfqnxrbGysm0ICAND4UHAAeIW8vDz9+OOPuvzyy8vt++mnn9S6dWv3BCkqkubPl7p1k8aNk3r1kr79VvriC6lvX0lSWlqa4uLiyhUus9msuLg4paWluScrAACNENdJAPAKe/bskdPprHAFJzs7W3369GnYAAUF0sKFUlycdOSIdMcdrlLTq1e5l86fP1+SNG3atIbNBAAAymEFB4BXqGyC2t69e+V0OnXllVc2zInz8qS33pK6dpUefVS64QZp1y5pyZIKyw0AADAWKzgAvEJycrIuueQSNW3atMz2devWSZIGDBhQvyfMyZHmzZNeeUXKyJB+9ztp+nTXpWkAAMBjUXAAeIXKBgxs27ZNkjRo0KD6OVF2tvTOO9Jrr0lnz0pjx7qKTdeu9XN8AADQoCg4ALxCSkqKYmJiym1PTk6W2WxWy5Yt63aCs2elOXOkN9+UcnNdAwSeekrq1KluxwUAAG5FwQHg8fLz87V///4KJ6gdOnRILVq0qP3BMzNdpeatt1wT0saPl558UurQoQ6JAQCAUSg4ADxeySCBii5RO3PmjK644ooLP+ipU9Lrr0t//atkt0uTJklTp0pt29ZDYgAAYBQKDgCPVzJB7fyCk5WVJZvNVmHxqdSJE9Krr0rvviv5+UkPPyxNmSK56zk6AACgQVFwAHi85ORktW/fXs2aNSuzff369ZKkvr88YLNKR4+6JqK9955ktUqTJ7s+wsIaIjIAADAIBQeAx6tsglpiYqIkafDgwZW/+fBhafZsaf58KShImjZNeuwxqS737QAAAI9FwQHg8ZKTkzV8+PBy25OSkmQymRQZGVn+TT/+KMXFSe+/LzVtKj33nPTII9J5q0AAAMC3+BkdAACqUlBQoP3795eu4IwfP7505Wb//v0KCQmRn5/rn7LExESNv/tu14jniAhp2TIpNlZKT5eeeYZyAwBAI8AKDgCPlpqaKofDUVpwxo4dq5iYGMXHx+vEiRNq3769JCnxH/9QzIMPKt5mk9q0cd1vM2GC67I0AADQaFBwAHi05ORkSb9OUIuKilJ8fLxiYmKUl5eniDZtlDh4sGLWrVN8q1aKeu456aGHpMBAI2MDAACDcIkaAI+WkpKitm3blnmYZ1RUlF5+9FFJUvHmzYpZv17xU6Yo6sgR6dFHKTcAADRirOAA8DiLFi1S+/btNWnSJO3fv1+SFBAQoI4dO2ru44/r6N//rtZJSbJIWiXpzVdeUdSUKYZmBgAAnoEVHAAep23btoqOjlZaWpqcTqecTqeKioqUlpam6EceUftTpxT29NMKbt5cb775pmbGxpYOHgAAAI2byel0Oo0OAQDnioiIUFpaWqX727drp9y8PMXHxysqKkqJiYmlgweioqLcmBQAAHgaCg4AjxMQEKCioqIqX7N58+YyZYaSAwAAJAoOAA9ktVpls9kq3e/v76/i4uJy2xMTE7Vo0SLNnz+/IeMBAAAPRsEB4HGqW8GxWq0qLCx0YyIAAOAtGDIAwON07NixTvsBAEDjRcEB4HHmzp0rk8lU4T6TyaS5c+e6OREAAPAWFBwAHufo0aNatWqVwsPDZbVaZbFYZLVaFR4erlWrVuno0aNGRwQAAB6Ke3AAAAAA+AxWcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn/H/AZbnZlLoDlpFAAAAAElFTkSuQmCC" - ], + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzgAAAEeCAYAAABG5rCAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4VPX99vF7kpkJ2diCyCabCQEFRVEQiaJAA/4gVFxRqiiVAlYfRagookVNJbiLVbAUFLHaYpFqUmRfg7K4pEgChIBhlS0BsyeTmXn+GBMJWckyZ2byfl1XLug5M+fcwUJy53vO55icTqdTAAAAAOAD/IwOAAAAAAD1hYIDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAADAg+3evVsREREKCAiQ1WpVQECAIiIitHv3bqOjAR7J5HQ6nUaHAAAAQHnJycnq3bu3iouLy+0zm83auXOnevToYUAywHNRcAAAADxURESE0tLSKt0fHh6uffv2uTER4PkoOAAAAB4qICBARUVFle63Wq0qLCx0YyLA83EPDgAAgIeq7ufQ/JwaKI+CAwAA4KFMJlOd9gONEQUHAADAAxUUSP7+Hat8TceOVe8HGiMKDgAAgIfJz5dGjpQcji/k72+u9HUPPfSQG1MB3oGCAwAA4EHy8qSYGGnLFunLL3vohx92Kjw8XFarVRaLRVarVZdeeqluu+02TZ8+XZ988onRkQGPUvmPBAAAAOBWubmucrN9u7R8uTRwoCT1qHAUtMPh0Lhx4/S73/1O/v7+uuuuu9yeF/BEFBwAAAAPkJMjDR8uffedtGKFFBVV9ev9/Py0YMEC2e123XvvvfL399ftt9/unrCAB+M5OAAAAAbLzpZuuUXaudNVbq6/vubvtdvtuu+++/Tpp59qyZIlGjVqVMMFBbwABQcAAMBAP//sKjfJydLKldJ11134MYqLizVmzBh99tlnWrp0qUaOHFn/QQEvQcEBAAAwyNmz0rBh0t690qpV0rXX1v5YNptN99xzj7744gstW7ZMw4cPr7+ggBeh4AAAABjgzBkpOlrav19avVrq06fux7TZbLrrrru0fPlyff755xo2bFjdDwp4GQoOAACAm2VkSL/5jXTokLRmjdS7d/0du6ioSHfccYdWrVqlL774QtHR0fV3cMALUHAAAADc6PRpacgQ6ehRae1a6Yor6v8chYWFuv3227V27VolJCRo8ODB9X8SwENRcAAAANzk1Clp8GDp+HFp3TqpZ8+GO1dBQYFGjRqljRs3asOGDerbt2/DnQzwIH5GBwAAAGgMTpyQbr5ZOnlS2rChYcuNJDVp0kTLli3T5MmTdfnllzfsyQAPwgoOAABAAzt+XBo0yDU1bd06qXt3oxMBvstsdAAAAABfduyYq9xkZ7tWbrp1MzoR4NsoOAAAAA3k6FHXZWn5+dLGjVJ4uNGJAN9HwQEAAGgAhw+7yo3N5io3XbsanQhoHCg4AAAA9ezgQVe5cTpdl6V16WJ0IqDxYIoaAABAPfrxR2ngQMlkotwARqDgAAAA1JMDB6SbbpLMZle56dTJ6ERA40PBAQAAqAdpaa6Vm4AA1z03l1xidCKgcaLgAAAA1FFqqqvcBAe7yk379kYnAhovCg4AAEAd7NnjKjfNmrkuS2vb1uhEQONGwQEAAKillBTXPTdhYdL69VKbNkYnAkDBAQAAqIVdu1yjoFu3dpWbiy82OhEAiYIDAABwwXbudJWbdu2kdeukiy4yOhGAEhQcAACAC5CU5Co3HTtKa9dKrVoZnQjAuSg4AAAANfTdd9KgQVLXrtKaNVLLlkYnAnA+Cg4AAEAN7NghDR4sRURIq1dLLVoYnQhARSg4AAAA1di2TfrNb6QePaRVq6TmzY1OBKAyFBwAAIAqfP21q9z07CmtXOl63g0Az0XBAQAAqMSWLVJ0tHTVVdKKFVJoqNGJAFSHggMAAFCBTZukoUOla6+Vli+XQkKMTgSgJig4AAAA51m/XrrlFum666SEBCk42OhEAGqKggMAAHCOtWul4cOlAQOk+HgpKMjoRAAuBAUHAADgF6tWSSNGSAMHSl98IQUGGp0IwIWi4AAAAEj68ktp5EjXs27+8x+pSROjEwGoDQoOAABo9BISpFtvdU1MW7pUCggwOhGA2qLgAACARu3zz6XbbnPdd/Pvf1NuAG9HwQEAAI3WsmXSHXe4Lk37178kq9XoRADqioIDAAAapX//W7rrLtfqzSefSBaL0YkA1AcKDgAAaHT+9S9p9Gjpzjulf/yDcgP4EgoOAABoVD7+WLr3XtfH4sWS2Wx0IgD1iYIDAAAajQ8/lO67z/Xx/vuSv7/RiQDUNwoOAABoFD74QHrgAdfHwoWUG8BXUXAAAIDPW7BAGjdOGj9emj9f8uM7IMBn8dcbAAD4tL/9TXroIWniRGnuXMoN4Ov4Kw4AAHzW3LnShAnSo49K77xDuQEaA/6aAwAAn/TXv0oPPyw99pj01luSyWR0IgDuQMEBAAA+5623XKs2U6ZIb7xBuQEaEwoOAADwKa+9Jj3+uDRtmvTKK5QboLGh4AAAAJ/x8svS1KnSM89Is2ZRboDGiIIDAAB8wksvuVZtnntOevFFyg3QWFFwAACA13vxRdeqzfPPuz4oN0DjRcEBAABeISEhQRaLRSaTqfTDYrFo9OgEPfecFBvrWr0B0LiZnE6n0+gQAAAAVfn888916623Vrr//vvjtWjRCDcmAuCpKDgAAMDjWSwWFRcXV7rfbDbLZrO5MREAT0XBAQAAHs9Ug5tq+JYGgMQ9OAAAAAB8CAUHAAAAgM+g4AAAAI/1c8HPemrNU1KLql9nNpvdEwiAx6PgAAAAj1PsKNbcHXMV/na43t7+tu557p4qX79s2TI3JQPg6Sg4AADAYzidTi3ft1xXzL1Cf1z+Rw2PGK7UR1L18eMfKz4+vtxKTcnwgbVr1xoRF4AHYooaAADwCDtP7NTUVVO1+sBq3dT5Jr0W/Zqubnt1le9xOBzq3LmzDh8+rDVr1mjw4MFuSgvAU1FwAACAoY7nHNez657Vgu8XKLxluF6NflUx3WJqNBpako4cOaKuXbvKbDbr2LFjat68eQMnBuDJuEQNAAAYIs+Wp9hNsQqfE66lu5fqzWFvatfDuzQycmSNy40kdejQQYsXL1Z+fr4GDBjQgIkBeANWcAAAgFs5nA59/MPHmr52uo7nHNejfR/VjBtnqEVgNaPSqjFmzBh9/PHHevTRRzVnzpx6SgvA21BwAACA22w+uFlTVk3RjmM7dFuP2zR7yGyFtwyvl2M7HA517dpVBw8e1IoVKzR06NB6OS4A70LBAQAADW5/5n5NWzNNS3cvVZ+2ffTG0Dd0Q6cb6v08x44dU+fOneXv76+jR4+qZcuW9X4OAJ6Ne3AAAECDOZN/RlNWTlGPd3po29Ft+vDWD7V9/PYGKTeS1K5dO3388ccqKCjgfhygkWIFBwAA1Dub3aa538zV8xufV2FxoZ6KekpP9H9CQZYgt5z//vvv1+LFi/Xwww/rnXfeccs5AXgGCg4AAKg3TqdT8anx+tPqPyktM03jeo/Ti4NeVJuQNm7N4XA4dOmllyo9PV3Lly/XLbfc4tbzAzAOBQcAANSLpONJemLlE1qfvl6DuwzWa9Gv6co2VxqW5/jx4+rUqZP8/Px0+PBhtWrVyrAsANyHe3AAAECdHMs+pnGfj9PV712tn3J+UsI9CVp932pDy40ktWnThvtxgEaIFRwAAFAruUW5evWrV/XyVy8ryBKk5296XuOvHi+Lv8XoaGU8+OCD2rhxo1JSUtSkSROj4wBoYBQcAABwQRxOhz7834d6Zt0zOp13Wo/1e0zTb5iu5k2aGx2tQg6HQ5Lk58eFK0BjQMEBAAA1tiF9g55Y+YS+P/697rr8LsUNjlOXFl2MjgUApcxGBwAAAJ4vNSNVT65+Up/v/Vz92vfTlnFbdP0l1xsdCwDKoeAAAIBKZeRl6IWNL+jdb95Vu9B2+uT2T3T35XfLZDIZHQ0AKsQlagAANEKLFi3SwIED1blz53L70tPTtXb9WmV1z9ILm16Q3WHX9Bum67F+jynQEuj+sABwASg4AAA0Qunp6Ro3bpwWLlxYpuT8+OOPGnnPSGUPy9Zhv8N66KqH9MLNL+jikIuNCwsAF4BxIgAANEKdO3fWwoULNW7cOKWnp0uS4rfF66r/u0q7rtul7uHd9b+J/9N7Me9RbgB4FVZwAABoxNLT03XH3XeoqGORfvjmB0WMi9Db97ytoeFDjY4GALXCkAEAABqpnKIcLTiwQEmOJNn/bdfvY3+veU/Pk9mPbw8AeC9WcAAAaGTsDrs+SPpAM9bPUOaxTLVZ00Y5R3JksVi0devWCgcPAIC34B4cAAAakTUH1ujqv12th+IfUt/gvuqzo482fr5RL730kk6cOKG77rqr9J4cAPBGFBwAABqB3ad2a8THI/Sbxb9RiDVEn0V/puwl2fr4w4/VuXNnjRs3TuHh4QoNDS0zeAAAvA0FBwAAH3Yq95QeWf6Ies3tpZRTKfr0zk+V+GCislKzyoyItlgsio2N1bp16zRp0iRt3LjR2OAAUEvcgwMAgA8qLC7UnG1z9JfNf5FTTs24YYYe7feompibVPoeh8Oha665RkFBQdq8ebNMJpMbEwNA/aDgAADgQ5xOp/6d8m9NWzNNh34+pInXTNSfB/5ZFwVfVKP3r1y5UsOGDVNCQoKGDx/ewGkBoP5RcAAA8BHbjmzTE6ue0FeHv9LwiOF6NfpVdW/V/YKO4XQ6NWjQIGVkZCgpKUl+flzNDsC78K8W4IF2796tiIgIBQQEyGq1KiAgQBEREdq9e7fR0QB4oINnD+repffqugXXKacoR6vvW62EexMuuNxIkslk0qxZs/TDDz/ok08+aYC0ANCwWMEBPExycrJ69+6t4uLicvvMZrN27typHj16GJAMgKfJKszSrM2z9MbWN9QisIX+MugvGnvlWPn7+df52Lfeeqt27typPXv2yGq11kNaAHAPCg7gYSIiIpSWllbp/vDwcO3bt8+NiQB4mmJHsRZ+v1DPrn9W2YXZmtJ/ip4c8KRCA0Lr7RzJycnq1auX3n77bf3xj3+st+MCQEOj4AAeJiAgQEVFRZXut1qtKiwsdGMiAJ5kZdpKTV09VbtO7tJ9V9ynlwa/pA5NOzTIuR544AGtWLFC+/fvV3BwcIOcAwDqG/fgAB6mup858DMJoHHadXKXhn00TMP+MUwtmrTQjvE79OGoDxus3EjSzJkzdebMGb311lsNdg4AqG8UHMDDVPfcCafNqe+ivtPBuIPK2ZVD4QF83Mnck5qYMFFXzrtS+8/s12d3faaND2zUNe2uafBzd+7cWZMmTdLs2bOVkZHR4OcDgPpAwQE8TMeOHavc3+GiDrJeZNXBFw/qm17faGuXrUp9JFUZKzJkL7C7KSWA+hAZGamJEyfKYrHIZDKVflgsFo3/w3hd1PEihc8J17+S/6VXf/Oqkh9O1qgeo9z6AM7p06fL4XBo9uzZbjsnANQF9+AAHmb37t3q1auX7PbyZeXcKWr2ArvObjirjIQMZf43UwXpBfIL8lOLIS0UNiJMYcPDFNAuwIDPAEBN/eEPf9D8+fMrf8E10mMvPqZnb3xWYUFh7gt2npkzZ2r27Nnat2+fOnRouEviAKA+UHAAD/TGG2/oiSeekMVikeS6bK1jx4764osvKhwR7XQ6lZeSp4z/ZigjPkM/f/Wz5JBC+oS4ys6IMIVeHSqTn/t+6gugehaLpcKR8CX8zf4qtlW+312ysrJ06aWXatSoUfrb3/5mdBwAqBIFB/BAY8eO1XfffacffvihVu+3ZdqU+WWma3VnRaaKzxbLcrFFYcNdZafFkBYyh5rrOTWAC1WTS8085cv0m2++qalTpyo5OVmRkZFGxwGASlFwAA9jt9vVpk0bPfTQQ5o1a1adj+ewOZT1VZYyEjKUkZChvD15MllNan5T89LVncAugfWQHEC1MjOlTZukDRukDRtk/d//ZKvmLZ7yZbqgoECRkZHq16+flixZYnQcAKgUBQfwMFu3blX//v2VmJioAQMG1Pvx8/fnu8rOfzN0dsNZOW1OBV0WVFp2mvZvKj8z80eAenFeodHOnZLTKXXpIg0cKP8PPpCjmkN40pfp999/X+PGjdOOHTt0zTUNP8UNAGqDggN4mBkzZmju3Lk6efKk/P39G/RcxdnFOrPqjOvenf9myHbSJnMLs1re0lJhI8LUcmhLWVpaGjQD4FMqKzSdO0s33eT6GDjQ9b9V/T04ZrNZNlt1azzuU1xcrCuuuEIdOnTQqlWrjI4DABWi4AAepnfv3urVq5cWL17s1vM6HU5lf5OtjHjXpWw5STmSv9Ts+mYKi3Gt7gR1D3LreFrA41VXaAYOlG6+WerUqcK3T5w4Ue+9916lh58wYYLmzZvXEMlrbdmyZbrtttu0Zs0aDR482Og4AFAOBQfwIIcPH1bHjh31z3/+U3fffbehWQqPFpau7JxZfUaOfIeadG1Seilb8xubyy+AS9nQyFRzyZluvtn1ayWF5nyRkZG6+eabtWDBgnIrOQ888IC++uor7d27t/4/jzpwOp267rrrJLkuqeWHHgA8DQUH8CDz5s3TI488otOnT6t58+ZGxyllz7fr7LqzrsKTkKHCw4XyD/FXi9+4nrnT8v9aKqANz9yBD7rAS87qQ9++fbVjxw6Fh4dr37599Xbc+rR+/XoNGjRIn332mUaNGmV0HAAog4IDeJARI0YoNzdX69evNzpKpZxOp3J/yC2dypa1NUtySqHXhrpWd2LCFNI7hJ/qwjvV8ZKz+lBQUKDAQNdkQ0/+Ej106FAdPnxYO3fulNnM2HkAnoOCA3iIvLw8hYWFKTY2VlOmTDE6To0VnSr69Zk7KzNlz7LL2s766zN3BreQf3DDDksAzjV+/HiNHTtWUVFR5fYlJiZq0aJFmj9/vmtDPV9yVl+aNGmiwsJCvfHGG3r88cfdeu6a+vbbb3XNNddo4cKFevDBB42OAwClKDiAh0hISFBMTIz27NnjtQ/RcxQ59HPiz6WrO/n78mUKMKnFINelbGHDw9SkUxOjY8LHJSYmKiYmRvHx8WVKTmJiomJGjFD8k08q6uRJt11yVhuLFy/W/fffL4vFoqKiIkOzVOXuu+/W119/rdTUVDVpwt9tAJ6BggN4iIkTJ2rt2rVKTU31mcu78lLzSp+58/Omn+Usdiq4V/Cvz9zp11Qmf9/4XOFZSkrO3FdeUdiRIwpMTlbMsmWKt9sVJbn1krPaKvl3IDc3V0FBQQanqVhqaqouu+wyvfLKK5o8ebLRcQBAEgUH8AhOp1OXXHKJ7rzzTr3xxhtGx2kQxT8XK3PVL5eyLc+U7bRN5jCzwv7vl0vZolvI0pxn7qCOzrnkLDEhQTfu369ASVY/P8VHRyvqnnsMueSsNq666iolJSWpR48eSklJMTpOpSZMmKClS5fqwIEDatq0qdFxAICCA3iCpKQkXXXVVY3muRJOu1NZ27NKn7mT+0OuTGaTmkU1K13dCYr0zJ9Yw8NUNRTg5pvVY8UK7fnpJ7388sv605/+ZHDYC+MtwwaOHj2q8PBwPfnkk3r++eeNjgMAFBzAE8TGxurll1/W6dOnZbVajY7jdgWHCkqfuXN27Vk5ChwKjAgsLTvNoprJz8ozd6ALmnKWmJiooUOHKi8vT0FBQVq5cmWFgwc8WUBAgIqKijRv3jxNmDDB6DiVevLJJ/Xuu+/qwIEDat26tdFxADRyFBzAA3z22Wfas2ePpk+fbnQUw9nz7Dqz5kzpM3eKjhXJP9RfLYe2dD1z55aWsrYuWwIXLVqkgQMHqnMFN4anp6dr48aNGjt2rJs+A9SrWk45K7kHZ8mSJYqOjtZ1112nPXv2lBs84On+/ve/a/z48bJarSosLDQ6TqUyMzPVtWtXjR07Vm+99ZbRcQA0chQcAB7L6XQqJymndCpb9o5sSVLTfk1Ln7kT3CtYBw8e1Lhx47Rw4cIyJSc9Pb3C7fBglRWaTp1+LTM33VTllLPzp6iFhYXJbDZr6dKlFU5X83TeMGxAkl566SXNnDlTqamp/H0DYCgKDgCvUXSiSBnLXZeynVl5RvYcuwIuCVDY8DDlXpurJz58Qu9/8L46d+6s0aNHa8eOHTp48KDsdnvpMcxms2JjY5WWlvbrs1AayAU9j6WxqsklZ9UUmvOd/+c+ZMgQrV27VoWFhdq+fbvX/bn37NlTycnJ6tWrl3bu3Gl0nErl5ubq0ksv1dChQ7Vo0SKj4wBoxCg4ALySo9Chs5vOugYV/DdDBQcKdCLghKb4TdHvon+nJd8v0Y+Hfqz0/XFxcZo2bVqDZqzyeSwevJLQoMXMgAdr/vWvf9Wjjz6qjz76SGPGjKm347pLXl6egoODJXn2sAFJevfdd/XII49o586d6tmzp9FxADRSFBwAXs/pdCpvb55Of3Fal02/THn2PIUpTBnKqPQ9ZrNZNputwbOdX2Y8vdxI1RezqKgoRUVFacaMGSouLi7dX+HqWDVTzjRwYIM/WDMrK0vNmjXTrbfeqmXLljXYeRpSybCB999/Xw888IDRcSpVVFSkHj16qFevXvrPf/5jdBwAjRQFB4DP2L59u/r166eZ02dq5kszq319r+3bZTGZZPXzk8VkKv299fzfn7Ot5Pcl77Ge97qSbee+bs/WrZo+erSihg7VV6tW6a+ffqr+UVGVns/fZDL8Ya9VFbNNmzbpmWeeqfS9cb/7naaFhVU75cydmjdvrqCgIB07dsyt560vc+bM0WOPPabWrVvrxIkTRsep0scff6wxY8boq6++Uv/+/Y2OA6ARouAA8Anp6ekaNWqUkpKStHfvXkVGRlb7nkdTU2VzOlXkcJT+WuR0lttmczpVVMk22y/vKXI4ZK/qZI895vpm//77pQcfrDbbueWpoiJUrpBVsP9CStj557OYTNq7bZueGT1a/YYM0fa1a0uL2eUtW6q4ioleZkm2BrzkrDYGDhyozZs3q6ioSGaz2dAsteXv7y+HwyG73S4/P88dm+5wOHTVVVepefPm2rBhg+FlHUDj453/ygPAOUqmpf32t79VSkqKunbtWqP3zYmIqNccjl/Kke28ovTVli26b/du2SSFfP65Xhk9Wpf171/j8lRyzMIKjn3+67LtdtmKiy+orFVazIKCpJ49tW7ZMunWW/WA1Spt3y7Zq6xyKpakAwfq9c+2rkaOHKlNmzbp888/1+233250nFoZPHiwVq9erfvvv18fffSR0XEq5efnp1mzZmn48OFauXKlhg0bZnQkAI0MKzgAvF7Jc3Di4uL09ddf63//+58sFkuZ+0PO5+57cHr27KnExERt3rzZ4+7BcTidKv6lDBWeU3y2bNmicSNHKj8vTwGBgXpz6VJd1r+/BrZoUe0xPe1Ly+nTp3XRRRfpjjvu0Keffmp0nFopGTbgrv/v1oXT6dTAgQOVnZ2tb7/91qNXnAD4Hv7FAeD1xo4dq86dOyslJUWXXXaZJCk2NrbK91S3vz6ce9+K2WyWn5+foqKiFB8fr5iYGCUmJjZ4hprw++XStGB/f7W0WHSx1aqD33yjSXfcoXlz50qS7rjtNj19773y27XL4LS106pVK4WGhmrr1q1GR6m1oKAgtWjRQsXFxR7z/53KmEwmzZo1S0lJSVqyZInRcQA0MhQcAD7B6XQqOTm5tOCkpaUpLi6u3PX/ZrNZcXFxSktLa/BMixYtKl2pyc7Olr+/vySVlhxPfVbIucVs9OjRklx/viXFrLqfxnvqPS49e/bU0aNH5XA4jI5Sa2+//bYk6d577zU4SfUGDBigESNGaMaMGR6/4gTAt1BwAPiEkydPKjMzU5dffrkkaf78+Zo2bZpiYmIkSS+88IKcTqdsNpumTZvmlgc9zp8/v/QytNzcXFksltJ9UVFRHvuwyXOLmdVqlZ+fn3766afSYtanT58q3++O1bHaGD58uJxOp1asWGF0lFobM2aM/Pz8dPjwYa8oan/5y1904MABLViwwOgoABoRCg4An5CSkiJJpSs4JZKSkiRJkyZNcnumc+Xm5spqtRqaoabOLWaSZLVaderUKUmuYnbllVcqLi6u3EqNO1fHauPBX6bX/eMf/zA4Sd3ceOONklwPZPV0V1xxhcaMGaPnn39eeXl5RscB0EgwZACAT3jnnXc0efJk5eXllfnGOzAwUAUFBYbf9N66dWv5+/vrp59+MjRHbYSFhSkwMFBHjhypcF9mZqbhf741FRISolatWik9Pd3oKLVW8uBSi8WioqIio+NU68CBA+revbteeOEFPfXUU0bHAdAIsIIDwCckJycrMjKy3KpCYWGhR0xwKiwsVFBQkNExaqVZs2bKzs6ucF+bNm0kuf78vUGPHj285vKuyjRt2lTNmzeXzWbTjh07jI5Tra5du2rChAmaPXu2zpw5Y3QcAI2A8V/1AaAenDtB7VxOp9MjikVRUZFH5KiNli1bKj8/v8J93bp1kyR99tln7oxUa7fccoscDoc2bNhgdJQ6ee211yRJd911l8FJaqZk0MDs2bONjgKgEaDgAPAJycnJpQMGzt0mSe3btzciUhnFxcUKCQkxOkatXHzxxZVOwSq5V+frr792Z6Ra+/3vfy9JWrx4scFJ6mbcuHEymUxKT0/3itWoiy++WJMnT9acOXN07Ngxo+MA8HEUHABe79SpUzp9+nS5FZx58+ZJkvr3729ErDLsdrtCQ0ONjlErJQXx7Nmz5fbdfffdkuSxgwXO16lTJwUGBmrTpk1GR6mz66+/XpL0xz/+0eAkNTN16lQFBgbqxRdfNDoKAB9HwQHg9SqboLZx40ZJ0kMPPeT2TOdzOp1q1qyZ0TFqpVOnTpIqvs+mQ4cOklQ6Zc0bREZG6uDBg0bHqLOEhARJ0sKFCw1OUjPNmjXT9OnTNX/+fO3bt8/oOAB8GAUHgNe0b8NdAAAXaklEQVRLTk6W2WxWREREme0l38QOGDDAiFilsrKyJEktWrQwNEdthYeHS5JSU1MrfU1ubq674tRZdHS07Ha7tmzZYnSUOmnevLmaNm2qoqIi7dy50+g4NfLwww+rbdu2eu6554yOAsCHUXAAeL2UlBR169atzIM0Jdc33SaTyaBUvyoZDe2tBScyMlKSa9xvRfz8/FRcXOzOSHUybtw4Sd5/H44kxcXFSZJGjRplcJKaCQwM1MyZM/XPf/5T33//vdFxAPgoCg4Ar1fZBDW73e4RD9csKTitWrUyOEntdO/eXZJ06NChCvdbrVaveQ6O5CpsAQEBWr9+vdFR6mzSpEkymUw6cOCAVwwbkKSxY8cqMjJS06dPNzoKAB9FwQHg9SqaoJaTkyPJ9SBKo504cUKS62Gf3qhJkyYymUyVTr/yxulwERER+vHHH42OUS/69u0rSZo8ebLBSWrGbDYrNjZWK1as8Ppx3QA8EwUHgFc7ffq0Tp48WW4FZ8GCBZKknj17GhGrjJIb8C+66CKDk9Se1WrVyZMnK9xX8rDP7777zp2R6mTIkCGy2WxelbkyJcMG3nvvPYOT1Nztt9+uPn366Omnn/aq1T8A3oGCA8Cr7d69W1L5CWr/+c9/JLm+kTLa6dOnJf1aBLxRUFBQpU+hL7mEreTP3Bs8+OCDkqQPPvjA2CD1oFWrVgoJCVFhYWHpREFPZzKZFBcXp61bt+qLL74wOg4AH0PBAeDVkpOT5e/vr27dupXZXrJqcv/99xsRq4zMzExJUrt27QxOUntNmzYtnQZ3vptuukmStG3bNjcmqpsrrrhCVqtVa9euNTpKvSh5toy3DBuQXKtogwcP1vTp02W3242OA8CHUHAAeLWUlBRFRESUGyZgt9sVGBioJk2aGJTsVyUFx1vvwZFc9zLl5+dXuK9klWz//v3ujFRnXbt29ZoHlFbn8ccfl8lkqnKUtyeaNWuWUlJS9NFHHxkdBYAPoeAA8GqVTVA7duyYxxSKs2fPSnKNU/ZWrVu3ls1mq3BfyaV33jQqWpJuvvlmFRUVVfgAU2909dVXS5KmTZtmcJKau/baa3X77bfrueeeU2FhodFxAPgI7/1qCwCqeIKaJGVnZ6tr164GJCovOzvbq8uNJLVv315Op7N0Ot35goODFRAQ4OZUdeNL9+FIvw4bmDNnjsFJLkxsbKyOHDniVUMSAHg27/6KC6BRy8zM1PHjx8ut4Ozbt09Op1NXXnmlQcnKys7OltlsNjpGnXTs2FHSr0Mdzte0aVNlZGS4M1KdjB8/XuvWrZMkvfrqqzKZTDKZTLJYLJo9e7bGjx9vcMIL16ZNGwUHB6ugoMCrLlXr3r27HnzwQcXGxio7O9voOAB8AAUHgNcq+Wb7/BWckm9cBwwY4PZMFcnJyZHFYjE6Rp1ceumlkqS9e/dWuL9169aVru54oi5duuipp54qt724uFhPPfWUwsPDDUhVd88++6wk6dZbbzU4yYX585//rKysLL3xxhtGRwHgAyg4ALxWcnKy/Pz8yk1Q27p1qyRp0KBBRsQqJz8/3+su3zpfZGSkJFV6U36HDh286h6KP//5z1XunzFjhpuS1J/IyMjSh5fu3r27zKrUxIkTS/8beqJLLrlEjzzyiF599dXSCYgAUFvefc0EgEYtJSVF4eHh5cpDcnKyzGazWrZsaVCysvLz8z1imltdlFwGePDgwQr3l9zvdOjQodLL2TxJfr508KCUnu76KC6u+uGSxcXFuu222xQcHKzg4GAFBQWV+bWmv/f393fL5ydJAwcOrPA+luLiYr333nuaMGGC27LUxtNPP6358+dr1qxZev31142OA8CLUXAAeK3KBgwcOnRILVq0MCBRxYqKijwqT20EBQXJZDLp2LFjFe7v0aOHJOn77783pOCcX2DO/zhx4tfXujqHqQbHzFdGRoZyc3NLP/Ly8pSbm6uCgoIa5bJarVWWpNqUpnOPc+69Xe+//36VWRYsWKB58+bVKLcRwsLCNHXqVMXGxurxxx/3yKIMwDtQcAB4rZSUlNJJWOfKzMxUr169DEhUsaKiIgUFBRkdo84sFotOnjxZZtv48eM1duxY9e7dW5L0ww8/6Le//a0kKTExUYsWLdL8+fPrfO4LLTCXXCJ17iz16CHdcovr9yUf7dtLFkv1I62//PLLSvfZ7Xbl5eWVFp5zy0/J7yvad+623NzceilQJYWnujHd3jDGe/LkyXr77bc1c+ZMLVy40Og4ALwUBQeAVzp79qyOHTtWboJaVlaWbDZb6YqCJyguLlZoaKjRMWpl/PjxCg8P14wZM1RcXKykpCSZTCaZzWbFxsbq+PHjiomJ0ZIlSySpdHpXYmKiYmJiFB8fX6Pz1LbAdO8uDRvm+n2XLr8WmOqG1pnN5iq/4a9u6p2/v79CQ0Mb7L+r3W5Xfn5+uUJUVaF64YUXqj3unj0PKjS0r5o27avg4F7y87NW+x53CgkJ0bPPPqvHH39cU6dOrfAZVwBQHQoOAK+UkpIiqfwEtQ0bNkiS+vXr5+5IlXI4HF5bcKqbNhYXF6dp06YpJiZGJpNJBw8eLFNuoqKiJNV9BaZTJ1eB6dJFateu+gJTndjY2Ao/r3P3G8nf318hISEKCQmp8XtqUnBycnbqxImP5HQWy2QKUGjoVaWFJzS0rwIDw2UyVX/5XkP6wx/+oNdff10zZszQZ599ZmgWAN7J5HQ6q77TEgA80N///ndNmDBBOTk5CgwMLN3+5JNP6pVXXtGuXbsqvD/H3RwOh/z9/TV69Gh98sknRse5YBaLpdqVDpvNprVrEzVkyA1q2rSVbLZijRoVL7s9qtoC06WLq7x07vzrrx061L3AVOf8lalzP5/Y2FilpaXVy6V17lTT/1Z2e75ycpKUnb1dWVnblZ29Xfn5ab+8pvk5hedahYb2VUBAG3d9CqUWL16s+++/X1u3bvWoH1YA8A4UHABe6YknnlB8fLz27dtXZnt0dLTWrFmj4uJi+fkZPwk/MzNTYWFhmjRpkt59912j41ywmvw0/+KLnb8UmOaSfpbJ9KY6dnys9JKx8z9qcgkZLtzEiRMrnKJWYsKECZUOGbDZMpWdvaO08GRlbZfN5rrfKiDgkjKrPKGhfWQ2N+yKpN1uV+/evdWqVSutW7fO8FUlAN6FLzEAvFJlE9T279+v4OBgjyg3kkqnjnn7FLWqTJokFRYm6s037br33slaunSmPvqoT+nlaXCP9evXa8KECVqwYEG5Vanf//73Wr9+faXvtVhaqmXLoWrZcqgkyel0qrDwcJnCk57+ghyOXEkmBQX1+KXw9FPTptfW+/08/v7+eumllzRy5EitXr1a0dHR9XZsAL6PggPAK6WkpOi+++4rt/3EiRNq166dAYkq9tNPP0mSLrroIoOT1IzNbtOWw1uUkJqg+NR4yV+Sver3DB7suudm1aovFRUVpQceuK3cPThoeHv37pWkehkFbTKZ1KRJRzVp0lGtW98hSXI67crNTSmz0tOQ9/OMGDFC119/vaZPn64hQ4Z4zA8tAHg+Cg4Ar5OVlaUjR46Um7DkcDiUm5ur8PBwg5KVVzJW2ZMLTmZ+plakrVB8arxWpK3Q2YKzahPSRiMiRijNkiaH3VHpe/38/MqVmaioKMXHx1NyfIzJ5K+QkF4KCemltm3HSdIv9/N8X1p4MjOX6+jROZIks7mFQkOvLS08TZv2ldV68QWcz6S4uDjdeOONWrp0qe68884G+bwA+B4KDgCvU9kEtaSkJEnSVVdd5fZMlTl16pQkqXXr1gYn+ZXT6dSe03uUkJqghH0J2nJoi+xOu65ue7X+X9//pxHdRqhPuz7yM/kpfGZ4ldPG+vTpo9dff71ciSkpOYsWLaLg+DB//0A1a3a9mjW7vnTb+ffzHDv2N9lsrql0F3o/z4EDB9S3b1+NGTNGY8aMkeQqPh07dtTcuXN19OhRjR07tmE/SQBeh4IDwOukpKTIZDIpMjKyzPaSewxuvPFGI2JVKCMjQ5LUtm1bQ3MU2Yu0+eDm0kvP9p/Zr0BzoAZ3Hax3h7+r4RHD1b5p+3LvS0tLU1xcXJXTxiorMFFRUZSbRqjy+3m2VXI/z2VlVnlc9/NYJLn+3mzfvr3cOdLS0hQdHa1Vq1a581MD4CUoOAC8TkpKirp06aKgoKAy23fs2CFJuuGGG4yIVaGSgtOmjftH7Z7OO60v932p+NR4rdy/UlmFWWof2l4juo3Qm93e1KAugxRkCaryGCWjkqdNm+aOyPBBZe/ncV1m9uv9PNuVlbVD2dnbdPz4h5Ls8vNropCQ3goN7asJE/5Z6XGdTqcmTZpUbpIiAFBwAHidyiao7d69W1artVzxMdKZM2ckSS1btmzwczmdTqWcSlF8arwSUhP09ZGv5XA6dG27azWl/xTFdItR7za9GbkLw5W9n+f3kiS7Pe+X+3l2KDt7hzIy/qsjR05WeZxDhw65Iy4AL0PBAeB1UlJSdM8995TbfuTIEYWFhRmQqHI///yzTCZTg02AKiwu1MaDGxW/N14J+xKUfjZdwZZgDek6RH8b8TcN7zZcbULcv3oEXCh//yA1azZAzZoNKN1mMlkl2Sp9D4/yA1ARCg4Ar5Kdna1Dhw5VuILz888/q2/fvgakqlxWVla9l5uTuSe1fN9yxafGa9X+VcopylHHZh01ImKERnQboZu73Kwm5ib1ek7ACNWtNrIaCaAiFBwAXmX37t2SVG5E9PHjx2W329WzZ08jYlUqOztbZnPd/ql1Op364eQPit8br/jUeG0/6rrpul+HfnpqwFOKiYxRr9a9+GYPPqdjx45KS0urcj8AnI+CA8CrlIyI7t69e5nta9eulST179/f7ZmqkpeXJ6v1wp/wXlBcoPU/ri+9n+Zw1mGFWEM09NKhmnTNJN0ScYtaB3vO6GmgIcydO1fR0dEVXopmMpk0d+5cA1IB8HQUHABeJTk5WV26dFFwcHCZ7V999ZUkafDgwUbEqtSFFJyfsn8qvfRs9YHVyrPlqUvzLrq1+62K6RajGzvdqABzQAMnBjzH0aNHtWrVKk2aNEmHDh2S0+ks9xwcADifyckdegC8yPDhw2UymZSQkFBm+4033qgtW7bIbrcblOxX48ePV3h4eJXPjpk/f76cTqe+P/596bNpvjn2jfxMfurfob9GdBuhmG4xuuyiy7j0DACAC0DBAeBVunTpojvvvFMvv/xyme2XXHKJcnJySscyG+mll17SM888U+n+sVPHKuDGACXsS9Cx7GNqGtBUw8KHaUTECN0ScYtaBbVyY1oAAHwLBQeA18jJyVFoaKg++OADjR07tsy+Jk2aKDw8XLt27TIo3a8sFkuZlZtyTFL4W+EaETFCMZExuqHjDbL4W9wXEAAAH8Y9OAC8xp49eySVn6BWVFSkwsJCRUZGGhGrnCrLjSQ5pdRHUrn0DACABtAwT54DgAZQMkGtR48eZbYnJiZKkq655hq3ZyrH4VBNagvlBgCAhkHBAeA1kpOT1alTJ4WEhJTZvmnTJknSTTfdZECqX9jt0j//KV1xhbjYDAAA41BwAHiNlJSUcpenSdK3334rSbr22mvdHUkqLpY++kjq2VO65x6pQwdZ/P2rfEtdH/wJAAAqR8EB4DUqKzj79u1TYGCge4uDzSa9/77Uo4d0331SRIS0bZu0YoWe/ctfqnxrbGysm0ICAND4UHAAeIW8vDz9+OOPuvzyy8vt++mnn9S6dWv3BCkqkubPl7p1k8aNk3r1kr79VvriC6lvX0lSWlqa4uLiyhUus9msuLg4paWluScrAACNENdJAPAKe/bskdPprHAFJzs7W3369GnYAAUF0sKFUlycdOSIdMcdrlLTq1e5l86fP1+SNG3atIbNBAAAymEFB4BXqGyC2t69e+V0OnXllVc2zInz8qS33pK6dpUefVS64QZp1y5pyZIKyw0AADAWKzgAvEJycrIuueQSNW3atMz2devWSZIGDBhQvyfMyZHmzZNeeUXKyJB+9ztp+nTXpWkAAMBjUXAAeIXKBgxs27ZNkjRo0KD6OVF2tvTOO9Jrr0lnz0pjx7qKTdeu9XN8AADQoCg4ALxCSkqKYmJiym1PTk6W2WxWy5Yt63aCs2elOXOkN9+UcnNdAwSeekrq1KluxwUAAG5FwQHg8fLz87V///4KJ6gdOnRILVq0qP3BMzNdpeatt1wT0saPl558UurQoQ6JAQCAUSg4ADxeySCBii5RO3PmjK644ooLP+ipU9Lrr0t//atkt0uTJklTp0pt29ZDYgAAYBQKDgCPVzJB7fyCk5WVJZvNVmHxqdSJE9Krr0rvviv5+UkPPyxNmSK56zk6AACgQVFwAHi85ORktW/fXs2aNSuzff369ZKkvr88YLNKR4+6JqK9955ktUqTJ7s+wsIaIjIAADAIBQeAx6tsglpiYqIkafDgwZW/+fBhafZsaf58KShImjZNeuwxqS737QAAAI9FwQHg8ZKTkzV8+PBy25OSkmQymRQZGVn+TT/+KMXFSe+/LzVtKj33nPTII9J5q0AAAMC3+BkdAACqUlBQoP3795eu4IwfP7505Wb//v0KCQmRn5/rn7LExESNv/tu14jniAhp2TIpNlZKT5eeeYZyAwBAI8AKDgCPlpqaKofDUVpwxo4dq5iYGMXHx+vEiRNq3769JCnxH/9QzIMPKt5mk9q0cd1vM2GC67I0AADQaFBwAHi05ORkSb9OUIuKilJ8fLxiYmKUl5eniDZtlDh4sGLWrVN8q1aKeu456aGHpMBAI2MDAACDcIkaAI+WkpKitm3blnmYZ1RUlF5+9FFJUvHmzYpZv17xU6Yo6sgR6dFHKTcAADRirOAA8DiLFi1S+/btNWnSJO3fv1+SFBAQoI4dO2ru44/r6N//rtZJSbJIWiXpzVdeUdSUKYZmBgAAnoEVHAAep23btoqOjlZaWpqcTqecTqeKioqUlpam6EceUftTpxT29NMKbt5cb775pmbGxpYOHgAAAI2byel0Oo0OAQDnioiIUFpaWqX727drp9y8PMXHxysqKkqJiYmlgweioqLcmBQAAHgaCg4AjxMQEKCioqIqX7N58+YyZYaSAwAAJAoOAA9ktVpls9kq3e/v76/i4uJy2xMTE7Vo0SLNnz+/IeMBAAAPRsEB4HGqW8GxWq0qLCx0YyIAAOAtGDIAwON07NixTvsBAEDjRcEB4HHmzp0rk8lU4T6TyaS5c+e6OREAAPAWFBwAHufo0aNatWqVwsPDZbVaZbFYZLVaFR4erlWrVuno0aNGRwQAAB6Ke3AAAAAA+AxWcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn0HBAQAAAOAzKDgAAAAAfAYFBwAAAIDPoOAAAAAA8BkUHAAAAAA+g4IDAAAAwGdQcAAAAAD4DAoOAAAAAJ9BwQEAAADgMyg4AAAAAHwGBQcAAACAz6DgAAAAAPAZFBwAAAAAPoOCAwAAAMBnUHAAAAAA+AwKDgAAAACfQcEBAAAA4DMoOAAAAAB8BgUHAAAAgM+g4AAAAADwGRQcAAAAAD6DggMAAADAZ1BwAAAAAPgMCg4AAAAAn/H/AZbnZlLoDlpFAAAAAElFTkSuQmCC", "text/plain": [ "PyPlot.Figure(PyObject )" ] @@ -535,13 +618,15 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 0.5.0-dev", + "display_name": "Julia 0.4.0", "language": "julia", - "name": "julia-0.5" + "name": "julia-0.4" }, "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", "name": "julia", - "version": "0.5.0" + "version": "0.4.0" } }, "nbformat": 4, diff --git a/src/dirichlet.jl b/src/dirichlet.jl index 5451760..35c36ae 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -27,13 +27,13 @@ type DBC2D2 <: DirichletEquation global_dofs :: Array{Int64, 1} fieldval :: Function end -function DBC2D2(el::Seg2) +function DBC2D2(element::Seg2) integration_points = [ IntegrationPoint([-sqrt(1/3)], 1.0), IntegrationPoint([+sqrt(1/3)], 1.0)] - new_fieldset!(el, "reaction force") + push!(element, FieldSet("reaction force")) fieldval(X, t) = 0.0 - DBC2D2(el, integration_points, [], fieldval) + DBC2D2(element, integration_points, [], fieldval) end function get_lhs(eq::DBC2D2, ip, t) el = get_element(eq) diff --git a/src/elements.jl b/src/elements.jl index d8a75b3..568bb0a 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -15,40 +15,25 @@ using ForwardDiff abstract Element +""" Get FieldSet from element. """ +function Base.getindex(element::Element, field_name::Union{Symbol, ASCIIString}) + element.fields[symbol(field_name)] +end + +""" Add new FieldSet to element. """ +function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name::Union{Symbol, ASCIIString}) + fieldset.name = symbol(fieldset_name) + element.fields[fieldset.name] = fieldset +end +function Base.push!(element::Element, fieldset::FieldSet) + element[fieldset.name] = fieldset +end + + #= ELEMENT DEFINITIONS -TODO: rewrite instructions after notebook. - -Each element must have - -1. Connectivity information. How element is connected to other elements. - This is typically node ids in Lagrange elements. -2. Ability to store fields, in array of shape dim × nnodes, where dim is - dimension of field and nnodes is number of nodes of element. Note that - this is always 2d array. -3. Default constructor which takes connectivity as argument. -4. Basis functions and derivative of basis functions. - -These rules probably will change, but there's a test_element function which -tests element and that it obeys current rules. If test_element passes, -everything should be ok. I use Quad4 as an example element here. -Several functions are inherited from Element abstract type: - -- get_connectivity -- get_number_of_basis_functions* -- get_element_dimension * -- get_basis * -- get_dbasisdxi * -- get_dbasisdX -- get_field -- set_field -- interpolate -- ... - -Which should work if element is defined following some rules. Functions marked with asterisk * are the ones which must necessarily to implement by your own. - -Start of example ----------------- +Example +------- This is example how to create new element. This is commented because I use code generation for simple elements like Lagrage elements. Feel free to use @@ -94,12 +79,6 @@ End of example. get_number_of_basis_functions(el::Type{Element}) = nothing get_element_dimension(el::Type{Element}) = nothing -### LAGRANGE ELEMENTS ### -#include("lagrange.jl") - -### HIERARCHICAL P-ELEMENTS ### -#include("hierarchical.jl") - ### COMMON ELEMENT ROUTINES ### """ @@ -115,10 +94,10 @@ Raises ------ This uses FactCheck and throws exceptions if element is not passing all tests. """ -function test_element(eltype) - Logging.info("Testing element $eltype") - local el - n = get_number_of_basis_functions(eltype) +function test_element(element_type) + Logging.info("Testing element $element_type") + local element + n = get_number_of_basis_functions(element_type) Logging.info("number of basis functions in this element: $n") @fact n --> not(nothing) """ Unable to determine number of nodes for $eltype define a function @@ -127,7 +106,7 @@ function test_element(eltype) Logging.info("Initializing element") try - el = eltype(collect(1:n)) + element = element_type(collect(1:n)) catch Logging.error(""" Unable to create element with default constructor define function @@ -135,31 +114,31 @@ function test_element(eltype) return false end - dim = get_element_dimension(eltype) + dim = get_element_dimension(element_type) Logging.info("Element dimension: $dim") @fact dim --> not(nothing) """ Unable to get element dimension define function 'get_element_dimension' which return the dimension of this element (1, 2, 3)""" # try to interpolate some scalar field - fld = Field(0.0, collect(1:n)) - Logging.info("Creating new scalar field $fld") - Logging.info("Pushing field to element.") - new_fieldset!(el, "field1") - add_field!(el, "field1", fld) - fieldset = get_fieldset(el, "field1") - @fact fieldset[1] --> fld + field = Field(0.0, collect(1:n)) + Logging.info("Creating new scalar field $field") + fieldset = FieldSet("field1") + push!(fieldset, field) + push!(element, fieldset) + @fact element["field1"][1] --> fld + # evaluate basis functions at middle point of element mid = zeros(dim) try - get_basis(el)(mid) + get_basis(element)(mid) catch Logging.error(""" Unable to evaluate basis, define function 'get_basis' for this element.""") end try - get_dbasisdxi(el)(mid) + get_dbasisdxi(element)(mid) catch Logging.error(""" Unable to evaluate partial derivatives of basis, @@ -167,46 +146,36 @@ function test_element(eltype) end Logging.info("Interpolating scalar field at $mid") - #f(field, xi, t) = el(xi)*el[field](t) - #i = f(:field1, mid, 0.0) - i = interpolate(el, "field1", mid, 0.0) + i = interpolate(element, "field1", mid, 0.0) Logging.info("Value: $i") - Logging.info("Element $eltype passed tests.") + Logging.info("Element $element_type passed tests.") end + get_connectivity(el::Element) = el.connectivity -""" -Get basis functions of element. -""" +""" Get basis functions of element. """ get_basis(el::Element) = el.basis get_basis(el::Element, xi::Vector) = el.basis(xi) Base.call(el::Element, xi::Vector) = el.basis(xi) -""" -Get partial derivatives of basis functions of element. -""" +""" Get partial derivatives of basis functions of element. """ get_dbasisdxi(el::Element) = el.basis.dbasisdxi get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi) -""" -Interpolate field on element. -""" -function interpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number) - fieldset = get_fieldset(el, symbol(field_name)) - field = interpolate(fieldset, t) - basis = get_basis(el) +""" Interpolate field on element. """ +function interpolate(element::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, time::Number) + fieldset = element[field_name] + field = interpolate(fieldset, time) + basis = get_basis(element) interpolate(basis, field, xi) end -""" -Interpolate derivative of field on element. -""" -function dinterpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number) - #get_dbasisdxi(el, xi)*el[field](t) - fieldset = get_fieldset(el, symbol(field_name)) - field = interpolate(fieldset, t) - basis = get_basis(el) +""" Interpolate derivative of field on element. """ +function dinterpolate(element::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, time::Number) + fieldset = element[field_name] + field = interpolate(fieldset, time) + basis = get_basis(element) dinterpolate(basis, field, xi) end @@ -215,24 +184,24 @@ Get jacobian of element evaluated at point ξ on element in reference configurat Parameters ---------- -el :: Element +element :: Element xi :: Vector -geometry_field :: Any, optional -time :: Number + spatial coordinate +time :: Float64 + temporal coordinate +geometry_field :: optional Returns ------- Vector or Matrix - depending on element type + depending on element dimension """ -function get_jacobian(el::Element, xi, t, geometry_field=symbol("geometry")) - dinterpolate(el, geometry_field, xi, t) +function get_jacobian(element::Element, xi, time, geometry_field="geometry") + dinterpolate(element, geometry_field, xi, time) end -""" -Evaluate partial derivatives of basis, dbasis/dX -""" +""" Evaluate partial derivatives of basis, dbasis/dX, at some time t""" function get_dbasisdX(el::Element, xi, t) dbasisdxi = get_dbasisdxi(el, xi) J = get_jacobian(el, xi, t) @@ -240,70 +209,8 @@ function get_dbasisdX(el::Element, xi, t) end -""" Create new empty set of fields for element. """ -function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString}) - el.fields[symbol(field_name)] = FieldSet() -end -function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field) - new_fieldset!(el, symbol(field_name)) - add_field!(el, symbol(field_name), field) -end -""" Add new field to fieldset of element. """ -function add_field!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field) - push!(el.fields[symbol(field_name)], field) -end - - -""" Get fieldset. """ -function get_fieldset(el::Element, field_name::Union{Symbol, ASCIIString}) - el.fields[symbol(field_name)] -end -""" Get fieldset, convenient function. """ -function Base.getindex(el::Element, field_name::Union{Symbol, ASCIIString}) - get_fieldset(el, field_name) -end - - - - -""" -calculate "local" normals in elements, in a way that -n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1] -""" -function calculate_normals!(el::Element, t, field_name=symbol("normals")) - new_field!(el, field_name, Vector) - for xi in Vector[[-1.0], [1.0]] - t = dinterpolate(el, :Geometry, xi) - n = [0 -1; 1 0]*t - n /= norm(n) - push_field!(el, field_name, n) - end -end - -""" -Alter normal field such that normals of adjacent elements are averaged. -""" -function average_normals!(elements, normal_field=symbol("normals")) - d = Dict() - for el in elements - c = get_connectivity(el) - n = get_field(el, normal_field) - for (ci, ni) in zip(c, n) - d[ci] = haskey(d, ci) ? d[ci] + ni : ni - end - end - for (ci, ni) in d - d[ci] /= norm(d[ci]) - end - for el in elements - c = get_connectivity(el) - new_normals = [d[ci] for ci in c] - set_field(el, normal_field, new_normals) - end -end - # FIXME: These two needs integration -- maybe not in elements.jl ..? """ @@ -326,7 +233,7 @@ function fit_field!(el::Element, field, f, fixed_coeffs=Int[]) xi = Vector[ [0.0], [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], + [-1/3*sqrt(5 - 2*sqrt(10/7))], [ 1/3*sqrt(5 + 2*sqrt(10/7))], [-1/3*sqrt(5 + 2*sqrt(10/7))]] n = get_number_of_basis_functions(el) @@ -390,7 +297,7 @@ function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) xi = Vector[ [0.0], [ 1/3*sqrt(5 - 2*sqrt(10/7))], - [-1/3*sqrt(5 - 2*sqrt(10/7))], + [-1/3*sqrt(5 - 2*sqrt(10/7))], [ 1/3*sqrt(5 + 2*sqrt(10/7))], [-1/3*sqrt(5 + 2*sqrt(10/7))]] n = get_number_of_basis_functions(el) @@ -442,67 +349,3 @@ function fit_derivative_field!(el::Element, field, f, fixed_coeffs=Int[]) end -""" Find projection from slave nodes to master element. """ -function calc_projection_slave_nodes_to_master_element(sel, mel) - X1 = get_field(sel, :Geometry) - N1 = get_field(sel, :Normals) - X2(xi) = interpolate(mel, :Geometry, xi) - dX2(xi) = dinterpolate(mel, :Geometry, xi) - R(xi, k) = det([X2(xi) - X1[k] N1[k]]') - dR(xi, k) = det([dX2(xi) N1[k]]') - xi2 = Vector[[0.0], [0.0]] - for k=1:2 - xi = xi2[k] - for i=1:3 - dxi = -R(xi, k)/dR(xi, k) - xi += dxi - if abs(dxi) < 1.0e-9 - break - end - end - xi2[k] = xi - end - clamp!(xi2, -1, 1) - return xi2 -end - -""" Find projection from master nodes to slave element. """ -function calc_projection_master_nodes_to_slave_element(sel, mel) - X1(xi) = interpolate(sel, :Geometry, xi) - dX1(xi) = dinterpolate(sel, :Geometry, xi) - N1(xi) = interpolate(sel, :Normals, xi) - dN1(xi) = dinterpolate(sel, :Normals, xi) - X2 = get_field(mel, :Geometry) - R(xi, k) = det([X1(xi) - X2[k] N1(xi)]') - dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]') - xi1 = Vector[[0.0], [0.0]] - for k=1:2 - xi = xi1[k] - for i=1:3 - dxi = -R(xi, k)/dR(xi, k) - xi += dxi - if abs(dxi) < 1.0e-9 - break - end - end - xi1[k] = xi - end - clamp!(xi1, -1, 1) - return xi1 -end - -function has_projection(sel, mel) - xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) - l = abs(xi1[2]-xi1[1])[1] - return l > 1.0e-9 -end - -""" -Calculate projection between 1d boundary elements -""" -function calc_projection(sel, mel) - xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) - xi2 = calc_projection_slave_nodes_to_master_element(sel, mel) - return xi1, xi2 -end - diff --git a/src/equations.jl b/src/equations.jl index 7acd838..c454595 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -3,6 +3,11 @@ abstract Equation +function get_unknown_field_name(equation::Equation) + eqtype = typeof(equation) + error("define get_unknown_field_name for this equation type $eqtype") +end + """ Integration point diff --git a/src/heat.jl b/src/heat.jl index c45dd49..1007d89 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -36,19 +36,19 @@ type DC2D4 <: HeatEquation integration_points :: Array{IntegrationPoint, 1} global_dofs :: Array{Int64, 1} end -function DC2D4(el::Quad4) +function DC2D4(element::Quad4) integration_points = [ IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0), IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0), IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0), IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)] - new_fieldset!(el, "temperature") - DC2D4(el, integration_points, []) + push!(element, FieldSet("temperature")) + DC2D4(element, integration_points, []) end -function get_lhs(eq::DC2D4, ip, t) - el = get_element(eq) - dNdX = get_dbasisdX(el, ip.xi, t) - k = interpolate(el, "temperature thermal conductivity", ip.xi, t) +function get_lhs(equation::DC2D4, ip, time) + element = get_element(equation) + dNdX = get_dbasisdX(element, ip.xi, time) + k = interpolate(element, "temperature thermal conductivity", ip.xi, time) return dNdX*k*dNdX' end JuliaFEM.has_lhs(eq::DC2D4) = true @@ -59,17 +59,15 @@ type DC2D2 <: HeatEquation integration_points :: Array{IntegrationPoint, 1} global_dofs :: Array{Int64, 1} end -function DC2D2(el::Seg2) - integration_points = [IntegrationPoint([0.0], 1.0)] - new_fieldset!(el, "temperature") - DC2D2(el, integration_points, []) +function DC2D2(element::Seg2) + integration_points = [IntegrationPoint([0.0], 2.0)] + push!(element, FieldSet("temperature")) + DC2D2(element, integration_points, []) end -function get_rhs(eq::DC2D2, ip, t) - el = get_element(eq) - h = get_basis(el, ip.xi) - f = interpolate(el, "temperature flux", ip.xi, t) - println("h = $h") - println("f = $f") +function get_rhs(equation::DC2D2, ip, time) + element = get_element(equation) + h = get_basis(element, ip.xi) + f = interpolate(element, "temperature flux", ip.xi, time) return h*f end JuliaFEM.has_rhs(eq::DC2D2) = true diff --git a/src/interpolate.jl b/src/interpolate.jl index f4c66ee..3921dc5 100644 --- a/src/interpolate.jl +++ b/src/interpolate.jl @@ -31,7 +31,7 @@ function interpolate(fields::FieldSet, t::Number) return Field(t, fields[1].values) end if t >= fields[end].time - return fields[end] + return Field(t, fields[end].values) end i = length(fields) while fields[i].time >= t diff --git a/src/lagrange.jl b/src/lagrange.jl index 92901c4..754d9c9 100644 --- a/src/lagrange.jl +++ b/src/lagrange.jl @@ -36,14 +36,14 @@ macro create_lagrange_element(element_name, element_description, X, P) global get_number_of_basis_functions, get_element_dimension dim = size($X, 1) nbasis = size($X, 2) - h = calculate_lagrange_basis($P, $X) + basis = calculate_lagrange_basis($P, $X) type $eltype <: CG connectivity :: Array{Int, 1} basis :: Basis - fields :: Dict{Symbol, Array{Field, 1}} + fields :: Dict{Symbol, FieldSet} end function $eltype(connectivity, args...) - $eltype(connectivity, Basis(h), Dict()) + $eltype(connectivity, Basis(basis), Dict()) end get_element_description(el::Type{$eltype}) = $element_description get_number_of_basis_functions(el::Type{$eltype}) = nbasis diff --git a/src/mortar.jl b/src/mortar.jl new file mode 100644 index 0000000..e7e48f8 --- /dev/null +++ b/src/mortar.jl @@ -0,0 +1,101 @@ +""" +calculate "local" normals in elements, in a way that +n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1] +""" +function calculate_normals!(el::Element, t, field_name=symbol("normals")) + new_field!(el, field_name, Vector) + for xi in Vector[[-1.0], [1.0]] + t = dinterpolate(el, :Geometry, xi) + n = [0 -1; 1 0]*t + n /= norm(n) + push_field!(el, field_name, n) + end +end + +""" +Alter normal field such that normals of adjacent elements are averaged. +""" +function average_normals!(elements, normal_field=symbol("normals")) + d = Dict() + for el in elements + c = get_connectivity(el) + n = get_field(el, normal_field) + for (ci, ni) in zip(c, n) + d[ci] = haskey(d, ci) ? d[ci] + ni : ni + end + end + for (ci, ni) in d + d[ci] /= norm(d[ci]) + end + for el in elements + c = get_connectivity(el) + new_normals = [d[ci] for ci in c] + set_field(el, normal_field, new_normals) + end +end + + +""" Find projection from slave nodes to master element. """ +function calc_projection_slave_nodes_to_master_element(sel, mel) + X1 = get_field(sel, :Geometry) + N1 = get_field(sel, :Normals) + X2(xi) = interpolate(mel, :Geometry, xi) + dX2(xi) = dinterpolate(mel, :Geometry, xi) + R(xi, k) = det([X2(xi) - X1[k] N1[k]]') + dR(xi, k) = det([dX2(xi) N1[k]]') + xi2 = Vector[[0.0], [0.0]] + for k=1:2 + xi = xi2[k] + for i=1:3 + dxi = -R(xi, k)/dR(xi, k) + xi += dxi + if abs(dxi) < 1.0e-9 + break + end + end + xi2[k] = xi + end + clamp!(xi2, -1, 1) + return xi2 +end + +""" Find projection from master nodes to slave element. """ +function calc_projection_master_nodes_to_slave_element(sel, mel) + X1(xi) = interpolate(sel, :Geometry, xi) + dX1(xi) = dinterpolate(sel, :Geometry, xi) + N1(xi) = interpolate(sel, :Normals, xi) + dN1(xi) = dinterpolate(sel, :Normals, xi) + X2 = get_field(mel, :Geometry) + R(xi, k) = det([X1(xi) - X2[k] N1(xi)]') + dR(xi, k) = det([dX1(xi) N1(xi)]') + det([X1(xi) - X2[k] dN1(xi)]') + xi1 = Vector[[0.0], [0.0]] + for k=1:2 + xi = xi1[k] + for i=1:3 + dxi = -R(xi, k)/dR(xi, k) + xi += dxi + if abs(dxi) < 1.0e-9 + break + end + end + xi1[k] = xi + end + clamp!(xi1, -1, 1) + return xi1 +end + +function has_projection(sel, mel) + xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) + l = abs(xi1[2]-xi1[1])[1] + return l > 1.0e-9 +end + +""" +Calculate projection between 1d boundary elements +""" +function calc_projection(sel, mel) + xi1 = calc_projection_master_nodes_to_slave_element(sel, mel) + xi2 = calc_projection_slave_nodes_to_master_element(sel, mel) + return xi1, xi2 +end + diff --git a/src/problems.jl b/src/problems.jl index b9426ca..0087934 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -22,6 +22,10 @@ function add_element!(problem::Problem, element::Element) equation = get_equation(typeof(problem), typeof(element)) push!(problem.equations, equation(element)) end +function Base.push!(problem::Problem, element::Element) + equation = get_equation(typeof(problem), typeof(element)) + push!(problem.equations, equation(element)) +end """ Return total number of basis functions in problem @@ -59,14 +63,15 @@ function set_global_dofs!(pr::Problem) end end -function get_connectivity(pr::Problem) - conn = Int[] - for eq in get_equations(pr) - el = get_element(eq) - append!(conn, get_connectivity(el)) +""" Return unique list of connectivity (i.e. node ids). """ +function get_connectivity(problem::Problem) + connectivity = Int[] + for equation in get_equations(problem) + element = get_element(equation) + append!(connectivity, get_connectivity(element)) end - conn = unique(conn) - return conn + connectivity = unique(connectivity) + return connectivity end """ diff --git a/src/solvers.jl b/src/solvers.jl index daf11ef..58532b2 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -5,12 +5,13 @@ abstract Solver -""" -Add new problem to solver -""" +""" Add new problem to solver. """ function add_problem!(solver::Solver, problem::Problem) push!(solver.problems, problem) end +function Base.push!(solver::Solver, problem::Problem) + push!(solver.problems, problem) +end """ Get all problems assigned to solver @@ -53,12 +54,10 @@ function call(solver::SimpleSolver, t) # assemble problem 2 A2 = sparse(get_lhs(problem2, t)...) b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1)) - + # make one monolithic assembly A = [A1 A2; A2' zeros(A2)] b = [b1; b2] - dump(full(A)) - dump(full(b)) # solve problem nz = unique(rowvals(A)) @@ -80,7 +79,7 @@ function call(solver::SimpleSolver, t) element = get_element(equation) field_name = get_unknown_field_name(equation) # field we are solving field = Field(t, full(x1[gdofs])[:]) - add_field!(element, field_name, field) + push!(element[field_name], field) end # update field for elements in problem 2 @@ -89,7 +88,7 @@ function call(solver::SimpleSolver, t) element = get_element(equation) field_name = get_unknown_field_name(equation) field = Field(t, full(x2[gdofs])) - add_field!(element, field_name, field) + push!(element[field_name], field) end end diff --git a/src/types.jl b/src/types.jl index dcd1f0d..747b448 100644 --- a/src/types.jl +++ b/src/types.jl @@ -5,7 +5,7 @@ using ForwardDiff -""" Field is a fundamental type which holds some values in some time t """ +""" Field is a fundamental data type which holds some values in some time t """ type Field{T} time :: Float64 increment :: Int64 @@ -19,6 +19,10 @@ end function Base.length(f::Field) length(f.values) end +""" Push value to field. """ +function Base.push!(f::Field, value) + push!(f.values, value) +end """ Get field discrete value at point i. """ function Base.getindex(f::Field, i::Int64) f.values[i] @@ -43,17 +47,39 @@ function Base.(:+)(f1::Field, f2::Field) Field(f1.time, f1.values + f2.values) end -""" -FieldSet is array of fields, each field maybe having different time and/or increment. -""" -typealias FieldSet Array{Field, 1} -""" Multiply fieldset with some vector x. """ -Base.(:*)(x::Array{Float64, 1}, f::Array{Field}) = sum(x .* f) -""" Add new field to fieldset. """ -function add_field!(fs::FieldSet, field::Field) - push!(fs, field) + +""" FieldSet is set of fields, each field can have different time and/or increment. """ +type FieldSet + name :: Symbol + fields :: Array{Field, 1} end +""" Initializer for FieldSet. """ +function FieldSet(field_name) + FieldSet(Symbol(field_name), []) +end +function FieldSet() + FieldSet(Symbol("unknown field"), []) +end +""" Add new field to fieldset. """ +function Base.push!(fs::FieldSet, field::Field) + push!(fs.fields, field) +end +""" Multiply fieldset with some vector x. """ +Base.(:*)(x::Array{Float64, 1}, fs::FieldSet) = sum(x .* fs.fields) +""" Get length of a fieldset. """ +function Base.length(fieldset::FieldSet) + length(fieldset.fields) +end +""" Return ith field from fieldset. """ +function Base.getindex(fieldset::FieldSet, i::Int64) + fieldset.fields[i] +end +#""" Return last field from fieldset. """ +function Base.endof(fieldset::FieldSet) + length(fieldset) +end + """ Basis function. """ @@ -70,7 +96,7 @@ diff(h::Basis) = h.dbasisdxi derivative(h::Basis) = h.dbasisdxi -# convenient functions +# convenient functions -- maybe this is not correct place for them """ Evaluate basis function in point ξ. """ call(b::Basis, xi) = b.basis(xi) #""" Interpolate field (h*f)(ξ) """