From c3c8283f8a18aba039b054e9992d7b258dc4e882 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 30 Sep 2015 09:16:44 +0300 Subject: [PATCH] dirichlet --- .../2015-08-29-developing-juliafem.ipynb | 223 +++++++++++++++--- src/elements.jl | 3 + 2 files changed, 198 insertions(+), 28 deletions(-) diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index 6600a57..bb8ca11 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -188,14 +188,14 @@ "name": "stderr", "output_type": "stream", "text": [ - "29-Sep 00:06:50:INFO:root:Testing element MyQuad4\n", - "29-Sep 00:06:50:INFO:root:number of basis functions in this element: 4\n", - "29-Sep 00:06:50:INFO:root:Initializing element\n", - "29-Sep 00:06:50:INFO:root:Element dimension: 2\n", - "29-Sep 00:06:50:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n", - "29-Sep 00:06:50:INFO:root:Pushing field to element.\n", - "29-Sep 00:06:50:INFO:root:Interpolating scalar field at [0.0,0.0]\n", - "29-Sep 00:06:51:INFO:root:Value: 2.5\n" + "30-Sep 09:08:39:INFO:root:Testing element MyQuad4\n", + "30-Sep 09:08:39:INFO:root:number of basis functions in this element: 4\n", + "30-Sep 09:08:39:INFO:root:Initializing element\n", + "30-Sep 09:08:40:INFO:root:Element dimension: 2\n", + "30-Sep 09:08:40:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n", + "30-Sep 09:08:40:INFO:root:Pushing field to element.\n", + "30-Sep 09:08:40:INFO:root:Interpolating scalar field at [0.0,0.0]\n", + "30-Sep 09:08:40:INFO:root:Value: 2.5\n" ] }, { @@ -244,7 +244,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "29-Sep 00:06:51:INFO:root:Element MyQuad4 passed tests.\n" + "30-Sep 09:08:40:INFO:root:Element MyQuad4 passed tests.\n" ] } ], @@ -762,7 +762,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "29-Sep 00:06:53:DEBUG:root:Problem (matrix) dimension: 4\n" + "30-Sep 09:08:43:DEBUG:root:Problem (matrix) dimension: 4\n" ] }, { @@ -864,9 +864,7 @@ { "data": { "text/plain": [ - "2x4 Array{Float64,2}:\n", - " 0.0 0.0 0.333333 0.166667\n", - " 0.0 0.0 0.166667 0.333333" + "InterfaceProblem" ] }, "execution_count": 25, @@ -875,12 +873,15 @@ } ], "source": [ - "M = 1/6*[0 0 2 1; 0 0 1 2]" + "type InterfaceProblem <: Problem\n", + " equations :: Array{Equation, 1}\n", + "end\n", + "InterfaceProblem() = InterfaceProblem([])" ] }, { "cell_type": "code", - "execution_count": 26, + "execution_count": 48, "metadata": { "collapsed": false }, @@ -888,32 +889,198 @@ { "data": { "text/plain": [ - "6-element Array{Float64,1}:\n", - " 100.0 \n", - " 100.0 \n", - " 1.36187e-14\n", - " 1.77636e-15\n", - " 600.0 \n", - " 600.0 " + "get_equation (generic function with 4 methods)" ] }, - "execution_count": 26, + "execution_count": 48, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "[A M'; M zeros(2, 2)] \\ [b; 0; 0]" + "abstract DirichletBC <: Equation\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, args...)\n", + "function DBC2D2(el::Seg2)\n", + " integration_points = [\n", + " IntegrationPoint([-sqrt(1/3)], 1.0),\n", + " IntegrationPoint([+sqrt(1/3)], 1.0)]\n", + " new_field!(el, \"reaction force\")\n", + " fieldval(X, t) = 0.0\n", + "# if length(args) != 0\n", + "# for (k, v) in args\n", + "# if k == :fieldval\n", + "# fieldval = v\n", + "# end\n", + "# end\n", + "# end\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", + " println(\"basis h: $h\")\n", + " out = h*h'\n", + " return out\n", + "end\n", + "\n", + "function JuliaFEM.get_rhs(eq::DBC2D2, ip, t)\n", + " el = get_element(eq)\n", + " h = get_basis(el)\n", + " #f = eq.fieldval\n", + " #X = interpolate(el, :Geometry, ip.xi, t)\n", + " #return h(ip.xi, t)*f(X, t)\n", + " return [0.0, 0.0]\n", + "end\n", + "JuliaFEM.has_lhs(eq::DBC2D2) = true\n", + "JuliaFEM.has_rhs(eq::DBC2D2) = true\n", + "JuliaFEM.get_dimension(pr::Type{InterfaceProblem}) = 1\n", + "JuliaFEM.get_equation(pr::Type{InterfaceProblem}, el::Type{Seg2}) = DBC2D2" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 49, "metadata": { - "collapsed": true + "collapsed": false }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "30-Sep 09:15:17:DEBUG:root:Problem (matrix) dimension: 2\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lhs\n" + ] + }, + { + "ename": "LoadError", + "evalue": "LoadError: MethodError: `call` has no method matching call(::JuliaFEM.Basis, ::Array{Float64,1}, ::Float64)\nClosest candidates are:\n BoundsError(!Matched::Any...)\n TypeVar(!Matched::Any...)\n TypeConstructor(!Matched::Any...)\n ...\nwhile loading In[49], in expression starting on line 17", + "output_type": "error", + "traceback": [ + "LoadError: MethodError: `call` has no method matching call(::JuliaFEM.Basis, ::Array{Float64,1}, ::Float64)\nClosest candidates are:\n BoundsError(!Matched::Any...)\n TypeVar(!Matched::Any...)\n TypeConstructor(!Matched::Any...)\n ...\nwhile loading In[49], in expression starting on line 17", + "", + " in get_lhs at In[48]:32", + " in integrate at /home/jukka/.julia/v0.4/JuliaFEM/src/equations.jl:69", + " in integrate_lhs at /home/jukka/.julia/v0.4/JuliaFEM/src/equations.jl:36", + " [inlined code] from In[49]:21", + " in anonymous at no file:0" + ] + } + ], + "source": [ + "# create elements and add necessary properties like connectivity and geometry\n", + "el3 = Seg2([1, 2])\n", + "new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n", + "\n", + "bc1 = InterfaceProblem()\n", + "add_element!(bc1, el3)\n", + "\n", + "# set global dofs for equations\n", + "set_global_dofs!(bc1)\n", + "\n", + "n = get_matrix_dimension(bc1)\n", + "\n", + "# integrate and assembly\n", + "t = 1.0\n", + "A = zeros(n, n)\n", + "b = zeros(n)\n", + "for eq in get_equations(bc1)\n", + " dofs = get_global_dofs(eq)\n", + " println(\"lhs\")\n", + " if has_lhs(eq)\n", + " A[dofs, dofs] += integrate_lhs(eq, t)\n", + " end\n", + " println(\"rhs\")\n", + " if has_rhs(eq)\n", + " b[dofs] += integrate_rhs(eq, t)\n", + " end\n", + "end\n", + "A2, b2" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "DBC2D2" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "JuliaFEM.get_equation(typeof(bc1), typeof(el3))" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "DBC2D2(JuliaFEM.Seg2([1,2],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "DBC2D2(el3)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "JuliaFEM.Seg2([1,2],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[]))" + ] + }, + "execution_count": 30, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "el3" + ] } ], "metadata": { diff --git a/src/elements.jl b/src/elements.jl index cae93f4..9c54896 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -276,6 +276,9 @@ end function new_field!(el::Element, field_name::ASCIIString, field::Field) new_field!(el, Symbol(field_name), field) end +function new_field!(el::Element, field_name::ASCIIString) + new_field!(el, Symbol(field_name)) +end """ Push to existing set field of fields. """