mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-18 03:19:21 +00:00
time properly implemented to .. everything.
This commit is contained in:
Jukka Aho
@@ -13,11 +13,35 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 1,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"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)"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using Logging\n",
|
||||
"using FactCheck\n",
|
||||
@@ -49,7 +73,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -133,7 +157,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_element_dimension (generic function with 8 methods)"
|
||||
"get_element_dimension (generic function with 7 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
@@ -164,21 +188,25 @@
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"28-Sep 12:56:52:INFO:root:Testing element MyQuad4\n",
|
||||
"28-Sep 12:56:52:INFO:root:number of basis functions in this element: 4\n",
|
||||
"28-Sep 12:56:53:INFO:root:Initializing element\n",
|
||||
"28-Sep 12:56:53:INFO:root:Element dimension: 2\n",
|
||||
"28-Sep 12:56:53:INFO:root:Setting scalar field JuliaFEM.Field{Int64}(0.0,1,[1,2,3,4]) to element.\n"
|
||||
"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"
|
||||
]
|
||||
},
|
||||
{
|
||||
"ename": "LoadError",
|
||||
"evalue": "LoadError: MethodError: `convert` has no method matching convert(::Type{Array{JuliaFEM.Field{T},1}}, ::JuliaFEM.Field{Int64})\nThis may have arisen from a call to the constructor Array{JuliaFEM.Field{T},1}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T}(::Type{Array{T,1}}, !Matched::Range{T})\n convert{T,S,N}(::Type{Array{T,N}}, !Matched::SubArray{S,N,P<:AbstractArray{T,N},I<:Tuple{Vararg{Union{AbstractArray{T,1},Colon,Int64}}},LD})\n ...\nwhile loading In[6], in expression starting on line 2",
|
||||
"output_type": "error",
|
||||
"traceback": [
|
||||
"LoadError: MethodError: `convert` has no method matching convert(::Type{Array{JuliaFEM.Field{T},1}}, ::JuliaFEM.Field{Int64})\nThis may have arisen from a call to the constructor Array{JuliaFEM.Field{T},1}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T}(::Type{Array{T,1}}, !Matched::Range{T})\n convert{T,S,N}(::Type{Array{T,N}}, !Matched::SubArray{S,N,P<:AbstractArray{T,N},I<:Tuple{Vararg{Union{AbstractArray{T,1},Colon,Int64}}},LD})\n ...\nwhile loading In[6], in expression starting on line 2",
|
||||
""
|
||||
]
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"PipeEndpoint(open, 0 bytes waiting)"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
@@ -195,31 +223,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"24-Sep 21:04:04:INFO:root:Value: [2.5]\n",
|
||||
"24-Sep 21:04:04:INFO:root:Element MyQuad4 passed tests.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: set_field, interpolate\n",
|
||||
"el1 = MyQuad4([1, 2, 3, 4])\n",
|
||||
"set_field(el1, :temperature, [1, 2, 3, 4])\n",
|
||||
"set_field(el1, :Geometry, 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",
|
||||
"set_field(el1, :\"heat coefficient\", 1);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 9,
|
||||
"execution_count": 7,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -227,22 +231,59 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2.5"
|
||||
"2-element Array{JuliaFEM.Field{T},1}:\n",
|
||||
" JuliaFEM.Field{Int64}(0.0,1,2)\n",
|
||||
" JuliaFEM.Field{Int64}(1.0,1,3)"
|
||||
]
|
||||
},
|
||||
"execution_count": 9,
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"29-Sep 00:06:51:INFO:root:Element MyQuad4 passed tests.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: new_field!, push_field!, interpolate, dinterpolate\n",
|
||||
"el1 = MyQuad4([1, 2, 3, 4])\n",
|
||||
"new_field!(el1, :temperature, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
|
||||
"push_field!(el1, :temperature, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
|
||||
"new_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_field!(el1, \"heat coefficient\", Field(0.0, 2))\n",
|
||||
"push_field!(el1, \"heat coefficient\", Field(1.0, 3))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1.25"
|
||||
]
|
||||
},
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# temperature at the middle poinf of the element, 1/4*(1+2+3+4)\n",
|
||||
"interpolate(el1, :temperature, [0.0, 0.0])"
|
||||
"# 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)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"execution_count": 9,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -256,14 +297,39 @@
|
||||
" 0.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"execution_count": 9,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# geometry midpoint of element\n",
|
||||
"interpolate(el1, :Geometry, [0.0, 0.0])"
|
||||
"interpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"3x2 Array{Float64,2}:\n",
|
||||
" 5.0 0.0\n",
|
||||
" 0.0 0.5\n",
|
||||
" 0.0 0.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"dinterpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -276,7 +342,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"1"
|
||||
"([0.625,0.625,0.625,0.625],[0.75,0.75,0.75,0.75])"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
@@ -286,7 +352,7 @@
|
||||
],
|
||||
"source": [
|
||||
"# interpolating scalar -> scalar.\n",
|
||||
"interpolate(el1, :\"heat coefficient\", [0.0, 0.0])"
|
||||
"interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5), interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -304,10 +370,9 @@
|
||||
"source": [
|
||||
"## Developing own formulation\n",
|
||||
"\n",
|
||||
"Let's consider Poisson equation\n",
|
||||
"Let's consider a Laplace equation\n",
|
||||
"\\begin{align}\n",
|
||||
"\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\n",
|
||||
"u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
|
||||
"\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}}\n",
|
||||
"\\end{align}\n",
|
||||
"\n",
|
||||
@@ -399,7 +464,7 @@
|
||||
" 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",
|
||||
" set_field(el, :temperature, zeros(2, 4)) # assign new field \"temperature\" to element\n",
|
||||
" new_field!(el, :temperature) # assign new field \"temperature\" to element\n",
|
||||
" DC2D4(el, integration_points, [])\n",
|
||||
"end"
|
||||
]
|
||||
@@ -435,10 +500,11 @@
|
||||
"\"\"\"\n",
|
||||
"Left hand side defined in integration point\n",
|
||||
"\"\"\"\n",
|
||||
"function JuliaFEM.get_lhs(eq::DC2D4, ip)\n",
|
||||
"function JuliaFEM.get_lhs(eq::DC2D4, ip, t)\n",
|
||||
" el = get_element(eq)\n",
|
||||
" dNdX = get_dbasisdX(el, ip.xi)\n",
|
||||
" hc = interpolate(el, :\"temperature thermal conductivity\", ip.xi)\n",
|
||||
" dNdX = get_dbasisdX(el, ip.xi, t)\n",
|
||||
" fld = el[\"temperature thermal conductivity\"](t)\n",
|
||||
" hc = sum(el(ip.xi) * fld)\n",
|
||||
" return dNdX*hc*dNdX'\n",
|
||||
"end\n",
|
||||
"JuliaFEM.has_lhs(eq::DC2D4) = true"
|
||||
@@ -476,17 +542,10 @@
|
||||
"source": [
|
||||
"using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
|
||||
"el = Quad4([1, 2, 3, 4])\n",
|
||||
"set_field(el, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
|
||||
"set_field(el, :\"temperature thermal conductivity\", 6)\n",
|
||||
"new_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_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
|
||||
"eq = DC2D4(el)\n",
|
||||
"integrate_lhs(eq)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"If rhs or lhs is not defined, integration returns nothing."
|
||||
"integrate_lhs(eq, 1.0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -508,14 +567,14 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"isa(integrate_rhs(eq), Void)"
|
||||
"JuliaFEM.has_lhs(eq)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next heat flux on boundary:"
|
||||
"If rhs or lhs is not defined, integration returns nothing."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -528,7 +587,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"has_rhs (generic function with 2 methods)"
|
||||
"true"
|
||||
]
|
||||
},
|
||||
"execution_count": 18,
|
||||
@@ -536,6 +595,35 @@
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"isa(integrate_rhs(eq, 1.0), Void)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next heat flux on boundary:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"has_rhs (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: get_basis, Seg2\n",
|
||||
"\n",
|
||||
@@ -551,25 +639,25 @@
|
||||
"function DC2D2(el::Seg2)\n",
|
||||
" integration_points = [\n",
|
||||
" IntegrationPoint([0], 2.0)]\n",
|
||||
" set_field(el, :temperature, zeros(2, 1))\n",
|
||||
" new_field!(el, :temperature)\n",
|
||||
" DC2D2(el, integration_points, [])\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"\"\"\"\n",
|
||||
"Right hand side defined in integration point\n",
|
||||
"\"\"\"\n",
|
||||
"function JuliaFEM.get_rhs(eq::DC2D2, ip)\n",
|
||||
"function JuliaFEM.get_rhs(eq::DC2D2, ip, t)\n",
|
||||
" el = get_element(eq)\n",
|
||||
" N = get_basis(el, ip.xi)\n",
|
||||
" f = interpolate(el, :\"temperature flux\", ip.xi)\n",
|
||||
" return f*N\n",
|
||||
" ϕ = get_basis(el)\n",
|
||||
" f = el[\"temperature flux\"]\n",
|
||||
" return ϕ(ip.xi)*f(t)\n",
|
||||
"end\n",
|
||||
"JuliaFEM.has_rhs(eq::DC2D2) = true"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -582,17 +670,17 @@
|
||||
" 50.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"el = Seg2([1, 2])\n",
|
||||
"set_field(el, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
|
||||
"set_field(el, :\"temperature flux\", 100.0)\n",
|
||||
"new_field!(el, :Geometry, Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n",
|
||||
"new_field!(el, \"temperature flux\", Field(0.0, 100.0))\n",
|
||||
"eq = DC2D2(el)\n",
|
||||
"integrate_rhs(eq)"
|
||||
"integrate_rhs(eq, 1.0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -607,7 +695,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -618,7 +706,7 @@
|
||||
"PlaneHeatProblem"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -627,14 +715,14 @@
|
||||
"using JuliaFEM: Problem, get_equation, get_dimension\n",
|
||||
"\n",
|
||||
"type PlaneHeatProblem <: Problem\n",
|
||||
" equations :: Array{Any, 1}\n",
|
||||
" equations :: Array{Equation, 1}\n",
|
||||
"end\n",
|
||||
"PlaneHeatProblem() = PlaneHeatProblem([])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -645,7 +733,7 @@
|
||||
"get_equation (generic function with 3 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 21,
|
||||
"execution_count": 22,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -674,7 +762,7 @@
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"24-Sep 21:04:17:DEBUG:root:Problem (matrix) dimension: 4\n"
|
||||
"29-Sep 00:06:53:DEBUG:root:Problem (matrix) dimension: 4\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -687,7 +775,7 @@
|
||||
" -2.0 -1.0 4.0 -1.0\n",
|
||||
" -1.0 -2.0 -1.0 4.0,\n",
|
||||
"\n",
|
||||
"[50.0,50.0,0.0,0.0])"
|
||||
"[300.0,300.0,0.0,0.0])"
|
||||
]
|
||||
},
|
||||
"execution_count": 23,
|
||||
@@ -700,11 +788,11 @@
|
||||
"\n",
|
||||
"# create elements and add necessary properties like connectivity and geometry\n",
|
||||
"el1 = Quad4([1, 2, 3, 4])\n",
|
||||
"set_field(el1, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
|
||||
"set_field(el1, :\"temperature thermal conductivity\", 6)\n",
|
||||
"new_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_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
|
||||
"el2 = Seg2([1, 2])\n",
|
||||
"set_field(el2, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
|
||||
"set_field(el2, :\"temperature flux\", 100.0)\n",
|
||||
"new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
|
||||
"new_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
|
||||
"\n",
|
||||
"problem = PlaneHeatProblem()\n",
|
||||
"add_element!(problem, el1)\n",
|
||||
@@ -716,20 +804,108 @@
|
||||
"n = get_matrix_dimension(problem)\n",
|
||||
"\n",
|
||||
"# integrate and assembly\n",
|
||||
"t = 1.0\n",
|
||||
"A = zeros(n, n)\n",
|
||||
"b = zeros(n)\n",
|
||||
"for eq in get_equations(problem)\n",
|
||||
" dofs = get_global_dofs(eq)\n",
|
||||
" if has_lhs(eq)\n",
|
||||
" A[dofs, dofs] += integrate_lhs(eq)\n",
|
||||
" A[dofs, dofs] += integrate_lhs(eq, t)\n",
|
||||
" end\n",
|
||||
" if has_rhs(eq)\n",
|
||||
" b[dofs] += integrate_rhs(eq)\n",
|
||||
" b[dofs] += integrate_rhs(eq, t)\n",
|
||||
" end\n",
|
||||
"end\n",
|
||||
"A, b"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2-element Array{Float64,1}:\n",
|
||||
" 100.0\n",
|
||||
" 100.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 24,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"fdofs = [1, 2]\n",
|
||||
"A[fdofs, fdofs] \\ b[fdofs]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We still need to consider Dirichlet boundary conditions:\n",
|
||||
"\\begin{align}\n",
|
||||
"u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
|
||||
"\\end{align}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 25,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"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"
|
||||
]
|
||||
},
|
||||
"execution_count": 25,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"M = 1/6*[0 0 2 1; 0 0 1 2]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 26,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"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 "
|
||||
]
|
||||
},
|
||||
"execution_count": 26,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"[A M'; M zeros(2, 2)] \\ [b; 0; 0]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
@@ -742,13 +918,15 @@
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 0.5.0-dev",
|
||||
"display_name": "Julia 0.4.0-rc2",
|
||||
"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,
|
||||
|
||||
+7
-1
@@ -8,7 +8,13 @@ using Logging
|
||||
@Logging.configure(level=DEBUG)
|
||||
|
||||
include("types.jl") # type definitions
|
||||
include("elements.jl") # elements
|
||||
|
||||
### ELEMENTS ###
|
||||
include("elements.jl")
|
||||
include("lagrange.jl") # Lagrange elements
|
||||
#include("hierarchical.jl") # P-elements
|
||||
|
||||
|
||||
include("equations.jl") # formulations
|
||||
include("problems.jl") # problems
|
||||
|
||||
|
||||
+68
-27
@@ -90,16 +90,13 @@ End of example.
|
||||
|
||||
# These must be implemented for your own element
|
||||
get_number_of_basis_functions(el::Type{Element}) = nothing
|
||||
get_number_of_basis_functions(el::Element) = nothing
|
||||
get_element_dimension(el::Element) = nothing
|
||||
get_dbasisdxi(el::Element, xi) = nothing
|
||||
get_connectivity(el::Element) = el.connectivity
|
||||
get_element_dimension(el::Type{Element}) = nothing
|
||||
|
||||
### LAGRANGE ELEMENTS ###
|
||||
include("lagrange.jl")
|
||||
#include("lagrange.jl")
|
||||
|
||||
### HIERARCHICAL P-ELEMENTS ###
|
||||
include("hierarchical.jl")
|
||||
#include("hierarchical.jl")
|
||||
|
||||
### COMMON ELEMENT ROUTINES ###
|
||||
|
||||
@@ -144,14 +141,15 @@ function test_element(eltype)
|
||||
|
||||
# try to interpolate some scalar field
|
||||
fld = Field(0.0, collect(1:n))
|
||||
Logging.info("Pushing scalar field $fld to element.")
|
||||
Logging.info("Creating new scalar field $fld")
|
||||
Logging.info("Pushing field to element.")
|
||||
new_field!(el, :field1)
|
||||
push_field!(el, :field1, fld)
|
||||
@fact el[:field1][1] --> fld
|
||||
|
||||
mid = zeros(dim)
|
||||
try
|
||||
f = get_basis(el)(mid)
|
||||
get_basis(el)(mid)
|
||||
catch
|
||||
Logging.error("""
|
||||
Unable to evaluate basis, define function 'get_basis' for
|
||||
@@ -167,11 +165,12 @@ function test_element(eltype)
|
||||
|
||||
Logging.info("Interpolating scalar field at $mid")
|
||||
f(field, xi, t) = el(xi)*el[field](t)
|
||||
i = f(:field, mid, 0.0)
|
||||
i = f(:field1, mid, 0.0)
|
||||
Logging.info("Value: $i")
|
||||
Logging.info("Element $eltype passed tests.")
|
||||
end
|
||||
|
||||
get_connectivity(el::Element) = el.connectivity
|
||||
|
||||
"""
|
||||
Get basis functions of element.
|
||||
@@ -180,6 +179,32 @@ 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_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::Symbol, xi::Vector, t::Number)
|
||||
get_basis(el, xi)*el[field](t)
|
||||
end
|
||||
function interpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
|
||||
interpolate(el, Symbol(field), xi, t)
|
||||
end
|
||||
|
||||
"""
|
||||
Interpolate derivative of field on element.
|
||||
"""
|
||||
function dinterpolate(el::Element, field::Symbol, xi::Vector, t::Number)
|
||||
get_dbasisdxi(el, xi)*el[field](t)
|
||||
end
|
||||
function dinterpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
|
||||
dinterpolate(el, Symbol(field), xi, t)
|
||||
end
|
||||
|
||||
"""
|
||||
Get jacobian of element evaluated at point ξ on element in reference configuration.
|
||||
|
||||
@@ -188,6 +213,7 @@ Parameters
|
||||
el::Element
|
||||
xi::Vector
|
||||
geometry_field::Any, optional
|
||||
time::Number, optional, default=0.0
|
||||
|
||||
Returns
|
||||
-------
|
||||
@@ -198,8 +224,8 @@ Notes
|
||||
-----
|
||||
Big "J" comes from reference (undeformed) configuration.
|
||||
"""
|
||||
function get_Jacobian(el::Element, xi, geometry_field=:Geometry)
|
||||
dinterpolate(el, geometry_field, xi)
|
||||
function get_Jacobian(el::Element, xi, t, geometry_field=:Geometry)
|
||||
dinterpolate(el, geometry_field, xi, t)
|
||||
end
|
||||
|
||||
|
||||
@@ -210,11 +236,11 @@ Notes
|
||||
-----
|
||||
Small "j" comes from current (deformed) configuration.
|
||||
"""
|
||||
function get_jacobian(el::Element, xi, geometry_field=:Geometry, displacement_field=:displacement)
|
||||
function get_jacobian(el::Element, xi, t, geometry_field=:Geometry, displacement_field=:displacement)
|
||||
dbasisdxi = get_dbasisdxi(el, xi)
|
||||
X = get_field(el, geometry_field)
|
||||
u = get_field(el, displacement_field)
|
||||
j = (X+u)*dbasisdxi
|
||||
X = get_field(el, geometry_field)(t)
|
||||
u = get_field(el, displacement_field)(t)
|
||||
j = dbasisdxi*(X+u)
|
||||
return j
|
||||
end
|
||||
|
||||
@@ -222,9 +248,9 @@ end
|
||||
"""
|
||||
Evaluate partial derivatives of basis, dbasis/dX
|
||||
"""
|
||||
function get_dbasisdX(el::Element, xi)
|
||||
function get_dbasisdX(el::Element, xi, t)
|
||||
dbasisdxi = get_dbasisdxi(el, xi)
|
||||
J = get_Jacobian(el, xi)
|
||||
J = get_Jacobian(el, xi, t)
|
||||
dbasisdxi*inv(J)
|
||||
end
|
||||
|
||||
@@ -232,32 +258,48 @@ end
|
||||
"""
|
||||
Evaluate partial derivatives of basis, dbasis/dx
|
||||
"""
|
||||
function get_dbasisdx(el::Element, xi)
|
||||
function get_dbasisdx(el::Element, xi, t)
|
||||
dbasisdxi = get_dbasisdxi(el, xi)
|
||||
j = get_jacobian(el, xi)
|
||||
j = get_jacobian(el, xi, t)
|
||||
dbasisdxi*inv(j)
|
||||
end
|
||||
|
||||
|
||||
""" Create new empty field of some type. """
|
||||
function new_field!(el::Element, field_name)
|
||||
function new_field!(el::Element, field_name::Symbol)
|
||||
el.fields[field_name] = Field[]
|
||||
end
|
||||
function new_field!(el::Element, field_name::Symbol, field::Field)
|
||||
new_field!(el, field_name)
|
||||
push_field!(el, field_name, field)
|
||||
end
|
||||
function new_field!(el::Element, field_name::ASCIIString, field::Field)
|
||||
new_field!(el, Symbol(field_name), field)
|
||||
end
|
||||
|
||||
|
||||
""" Push to existing set field of fields. """
|
||||
function push_field!(el::Element, field_name, field::Field)
|
||||
function push_field!(el::Element, field_name::Symbol, field::Field)
|
||||
push!(el.fields[field_name], field)
|
||||
end
|
||||
function push_field!(el::Element, field_name::ASCIIString, field::Field)
|
||||
push_field!(el, Symbol(field_name), field)
|
||||
end
|
||||
|
||||
|
||||
""" Get field variable. """
|
||||
function get_field(el::Element, field_name)
|
||||
function get_field(el::Element, field_name::Symbol)
|
||||
el.fields[field_name]
|
||||
end
|
||||
function Base.getindex(el::Element, field_name)
|
||||
el.fields[field_name]
|
||||
function get_field(el::Element, field_name::ASCIIString)
|
||||
el.fields[Symbol(field_name)]
|
||||
end
|
||||
function Base.getindex(el::Element, field_name::Union{ASCIIString, Symbol})
|
||||
get_field(el, field_name)
|
||||
end
|
||||
|
||||
|
||||
#=
|
||||
"""
|
||||
Evaluate some field in point ξ on element using basis functions.
|
||||
|
||||
@@ -296,8 +338,6 @@ function interpolate(el::Element, field, xis::Array{Vector, 1})
|
||||
map(interpolate_, xis)
|
||||
end
|
||||
|
||||
"""
|
||||
"""
|
||||
function dinterpolate(el::Element, field, xi::Number)
|
||||
dinterpolate(el, field, [xi])
|
||||
end
|
||||
@@ -309,12 +349,13 @@ function dinterpolate(el::Element, field, xi::Vector)
|
||||
end
|
||||
return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
|
||||
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, field_name=:Normals)
|
||||
function calculate_normals!(el::Element, t, field_name=:Normals)
|
||||
new_field!(el, field_name, Vector)
|
||||
for xi in Vector[[-1.0], [1.0]]
|
||||
t = dinterpolate(el, :Geometry, xi)
|
||||
|
||||
+12
-12
@@ -32,23 +32,23 @@ get_integration_points(eq::Equation) = eq.integration_points
|
||||
get_connectivity(eq::Equation) = get_connectivity(get_element(eq))
|
||||
get_basis(eq::Equation, ip::IntegrationPoint) = get_basis(get_element(eq), ip.xi)
|
||||
get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq), ip.xi)
|
||||
interpolate(eq::Equation, field::Union(ASCIIString, Symbol), ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
|
||||
integrate_lhs(eq::Equation) = has_lhs(eq) ? integrate(eq, get_lhs) : nothing
|
||||
integrate_rhs(eq::Equation) = has_rhs(eq) ? integrate(eq, get_rhs) : nothing
|
||||
interpolate(eq::Equation, field::Union{ASCIIString, Symbol}, ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
|
||||
integrate_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
|
||||
integrate_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
|
||||
|
||||
|
||||
"""
|
||||
Return determinant of Jacobian for numerical integration.
|
||||
"""
|
||||
function get_detJ(eq::Equation, ip::IntegrationPoint)
|
||||
function get_detJ(eq::Equation, ip::IntegrationPoint, t::Float64)
|
||||
el = get_element(eq)
|
||||
get_detJ(el, ip)
|
||||
get_detJ(el, ip, t)
|
||||
end
|
||||
function get_detJ(el::Element, ip::IntegrationPoint)
|
||||
J = get_detJ(el, ip.xi)
|
||||
function get_detJ(el::Element, ip::IntegrationPoint, t::Float64)
|
||||
get_detJ(el, ip.xi, t)
|
||||
end
|
||||
function get_detJ(el::Element, xi::Vector)
|
||||
J = get_Jacobian(el, xi)
|
||||
function get_detJ(el::Element, xi::Vector, t::Float64)
|
||||
J = get_Jacobian(el, xi, t)
|
||||
s = size(J)
|
||||
return s[1] == s[2] ? det(J) : norm(J)
|
||||
end
|
||||
@@ -63,10 +63,10 @@ eq::Equation
|
||||
f::Function
|
||||
Function to integrate
|
||||
"""
|
||||
function integrate(eq::Equation, f::Function)
|
||||
function integrate(eq::Equation, f::Function, t::Float64)
|
||||
target = []
|
||||
for ip in get_integration_points(eq)
|
||||
push!(target, ip.weight*f(eq, ip)*get_detJ(eq, ip))
|
||||
push!(target, ip.weight*f(eq, ip, t)*get_detJ(eq, ip, t))
|
||||
end
|
||||
return sum(target)
|
||||
end
|
||||
@@ -80,4 +80,4 @@ function set_global_dofs!(eq::Equation, dofs)
|
||||
end
|
||||
|
||||
# Equations for heat problems
|
||||
include("heat_equations.jl")
|
||||
#include("heat_equations.jl")
|
||||
|
||||
+30
-56
@@ -5,31 +5,9 @@
|
||||
|
||||
abstract CG <: Element
|
||||
|
||||
"""
|
||||
Create new element with element_name to family element_family
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> @create_element(Seg2, CG, "2 node linear segment")
|
||||
"""
|
||||
macro create_element(element_name, element_family, element_description)
|
||||
# Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
|
||||
eltype = esc(element_name)
|
||||
elfam = esc(element_family)
|
||||
quote
|
||||
global get_element_description
|
||||
type $eltype <: $elfam
|
||||
connectivity :: Array{Int, 1}
|
||||
fields :: Dict{Any, Any}
|
||||
end
|
||||
$eltype(connectivity) = $eltype(connectivity, Dict{Any, Any}())
|
||||
get_element_description(el::Type{$eltype}) = $element_description
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
Given polynomial P and coordinates of reference element, calculate
|
||||
Lagrange basis function and partial derivatives.
|
||||
Lagrange basis functions
|
||||
"""
|
||||
function calculate_lagrange_basis(P, X)
|
||||
dim, nbasis = size(X)
|
||||
@@ -40,70 +18,66 @@ function calculate_lagrange_basis(P, X)
|
||||
# Logging.debug("Calculating inverse of A")
|
||||
invA = inv(A)'
|
||||
basis(xi) = invA*P(xi)
|
||||
dbasisdxi = ForwardDiff.jacobian(basis)
|
||||
basis, dbasisdxi
|
||||
basis
|
||||
end
|
||||
|
||||
"""
|
||||
Assign Lagrange basis for element.
|
||||
Create new Lagrange element
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> @create_lagrange_element(Seg2, "2 node linear segment", X, P)
|
||||
"""
|
||||
macro create_lagrange_basis(element_name, X, P)
|
||||
|
||||
# Logging.debug("Creating Lagrange basis for element ", element_name, ". ")
|
||||
macro create_lagrange_element(element_name, element_description, X, P)
|
||||
# Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
|
||||
eltype = esc(element_name)
|
||||
|
||||
quote
|
||||
|
||||
global get_element_description
|
||||
global get_number_of_basis_functions, get_element_dimension
|
||||
global get_basis, get_dbasisdxi
|
||||
|
||||
dim = size($X, 1)
|
||||
nbasis = size($X, 2)
|
||||
# Logging.debug("Number of basis functions: ", nbasis, ". ")
|
||||
# Logging.debug("Element dimension: ", dim)
|
||||
|
||||
get_number_of_basis_functions(el::Type{$(esc(element_name))}) = nbasis
|
||||
get_number_of_basis_functions(el::$(esc(element_name))) = nbasis
|
||||
get_element_dimension(el::$(esc(element_name))) = dim
|
||||
|
||||
basis, dbasisdxi = calculate_lagrange_basis($P, $X)
|
||||
get_basis(el::$eltype, xi) = basis(xi)
|
||||
get_dbasisdxi(el::$eltype, xi) = dbasisdxi(xi)
|
||||
# Logging.debug("Element ", $element_name, " created.")
|
||||
h = calculate_lagrange_basis($P, $X)
|
||||
type $eltype <: CG
|
||||
connectivity :: Array{Int, 1}
|
||||
basis :: Basis
|
||||
fields :: Dict{Symbol, Array{Field, 1}}
|
||||
end
|
||||
function $eltype(connectivity, args...)
|
||||
$eltype(connectivity, Basis(h), Dict())
|
||||
end
|
||||
get_element_description(el::Type{$eltype}) = $element_description
|
||||
get_number_of_basis_functions(el::Type{$eltype}) = nbasis
|
||||
get_element_dimension(el::Type{$eltype}) = dim
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# 0d Lagrange element
|
||||
|
||||
@create_element(Point1, CG, "1 node point element")
|
||||
#@create_element(Point1, CG, "1 node point element")
|
||||
|
||||
# 1d Lagrange elements
|
||||
|
||||
@create_element(Seg2, CG, "2 node linear line element")
|
||||
@create_lagrange_basis(Seg2, [-1.0 1.0], (xi) -> [1.0, xi[1]])
|
||||
@create_lagrange_element(Seg2, "2 node linear line element",
|
||||
[-1.0 1.0], (xi) -> [1.0, xi[1]])
|
||||
|
||||
@create_element(Seg3, CG, "3 node quadratic line element")
|
||||
@create_lagrange_basis(Seg3, [-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
|
||||
@create_lagrange_element(Seg3, "3 node quadratic line element",
|
||||
[-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
|
||||
|
||||
# 2d Lagrange elements
|
||||
|
||||
@create_element(Tri3, CG, "3 node bilinear triangle element")
|
||||
@create_lagrange_basis(Tri3,
|
||||
@create_lagrange_element(Tri3, "3 node bilinear triangle element",
|
||||
[0.0 1.0 0.0
|
||||
0.0 0.0 1.0],
|
||||
(xi) -> [1.0, xi[1], xi[2]])
|
||||
|
||||
@create_element(Quad4, CG, "4 node bilinear quadrangle element")
|
||||
@create_lagrange_basis(Quad4,
|
||||
@create_lagrange_element(Quad4, "4 node bilinear quadrangle element",
|
||||
[-1.0 1.0 1.0 -1.0
|
||||
-1.0 -1.0 1.0 1.0],
|
||||
(xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2]])
|
||||
|
||||
# 3d Lagrange elements
|
||||
|
||||
@create_element(Tet10, CG, "10 node quadratic tetrahedron")
|
||||
@create_lagrange_basis(Tet10,
|
||||
@create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
|
||||
[0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
|
||||
0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5
|
||||
0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.5 0.5],
|
||||
|
||||
+2
-84
@@ -1,69 +1,8 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
## This module contains math stuff, including interpolation, integration, linearization, ...
|
||||
|
||||
using ForwardDiff
|
||||
|
||||
#export interpolate, integrate, linearize
|
||||
|
||||
"""
|
||||
Interpolate field variable using basis functions f for point ip.
|
||||
This function tries to be as general as possible and allows interpolating
|
||||
lot of different fields.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
field :: Array{Number, dim}
|
||||
Field variable
|
||||
basis :: Function
|
||||
Basis functions
|
||||
ip :: Array{Number, 1}
|
||||
Point to interpolate
|
||||
"""
|
||||
function interpolate(field::Float64, basis::Function, ip::Array{Float64,1})
|
||||
# dummy function, unable to interpolate scalar value!
|
||||
return field
|
||||
end
|
||||
function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip)
|
||||
result = dot(field, basis(ip))
|
||||
return result
|
||||
end
|
||||
function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
|
||||
m, n = size(field)
|
||||
bip = basis(ip)
|
||||
tmp = size(bip)
|
||||
if length(tmp) == 1
|
||||
ndim = 1
|
||||
nnodes = tmp[1]
|
||||
else
|
||||
ndim, nnodes = size(bip)
|
||||
end
|
||||
if ndim == 1
|
||||
if n == nnodes
|
||||
result = field * bip
|
||||
elseif m == nnodes
|
||||
result = field' * bip
|
||||
end
|
||||
else
|
||||
if n == nnodes
|
||||
result = bip' * field
|
||||
elseif m == nnodes
|
||||
result = bip' * field'
|
||||
end
|
||||
end
|
||||
if length(result) == 1
|
||||
result = result[1]
|
||||
end
|
||||
return result
|
||||
end
|
||||
#function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
|
||||
# basis = derivative ? get_dbasisdxi(e) : get_basis(e)
|
||||
# return interpolate(e.attributes[field], basis, x)
|
||||
#end
|
||||
|
||||
|
||||
|
||||
"""
|
||||
Linearize function f w.r.t some given field, i.e. calculate dR/du
|
||||
|
||||
@@ -78,7 +17,6 @@ Returns
|
||||
-------
|
||||
Array{Float64, 2}
|
||||
jacobian / "tangent stiffness matrix"
|
||||
|
||||
"""
|
||||
function linearize(f::Function, el::Element, field::ASCIIString)
|
||||
dim, nnodes = size(el.attributes[field])
|
||||
@@ -92,6 +30,7 @@ function linearize(f::Function, el::Element, field::ASCIIString)
|
||||
return jac(el.attributes[field][:])
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
This version returns another function which can be then evaluated against field
|
||||
"""
|
||||
@@ -111,6 +50,7 @@ function linearize(f::Function, field::ASCIIString)
|
||||
return jacobian
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
In-place version, no additional garbage collection.
|
||||
"""
|
||||
@@ -128,8 +68,6 @@ function linearize!(f::Function, el::Element, field::ASCIIString, target::ASCIIS
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
"""
|
||||
This version returns a function which must be operated with element e
|
||||
"""
|
||||
@@ -157,25 +95,6 @@ function integrate!(f::Function, el::Element, target)
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
Evaluate field in point xi using basis functions.
|
||||
"""
|
||||
function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})
|
||||
f = get_field(el, field)
|
||||
if !isa(f, Array)
|
||||
# This is scalar, nothing to interpolate
|
||||
return f
|
||||
end
|
||||
basis = get_basis(el, xi)
|
||||
dim, nnodes = size(f)
|
||||
result = zeros(dim)
|
||||
for i=1:nnodes
|
||||
result += basis[i]*f[:,i]
|
||||
end
|
||||
return result
|
||||
end
|
||||
|
||||
|
||||
function linearize(eq::Equation, f::Function, field::ASCIIString)
|
||||
function jacobian(eq::Equation, args...)
|
||||
el = get_element(eq)
|
||||
@@ -194,4 +113,3 @@ function linearize(eq::Equation, f::Function, field::ASCIIString)
|
||||
return jacobian
|
||||
end
|
||||
|
||||
|
||||
|
||||
+11
-3
@@ -10,7 +10,7 @@ using ForwardDiff
|
||||
type Field{T}
|
||||
time :: Float64
|
||||
increment :: Int64
|
||||
values :: Array{T, 1}
|
||||
values :: T
|
||||
end
|
||||
|
||||
""" Initialize field. """
|
||||
@@ -25,8 +25,16 @@ Base.length(f::Field) = length(f.values)
|
||||
Base.getindex(f::Field, i::Int64) = f.values[i]
|
||||
|
||||
""" Interpolate field h(ξ)*f = x*f """
|
||||
Base.(:*)(x::Array{Float64, 1}, f::Field) = sum(x .* f.values)
|
||||
Base.(:*)(x::Array{Float64, 2}, f::Field) = sum([f[i]*x[i,:] for i in 1:length(f)])
|
||||
function interpolate{T}(x::Vector, f::Field{Vector{T}})
|
||||
sum([f[i]*x[i] for i in 1:length(f)])
|
||||
end
|
||||
function interpolate{T}(x::Matrix, f::Field{Vector{T}})
|
||||
sum([f[i]*x[i,:] for i in 1:length(f)])
|
||||
end
|
||||
function interpolate(x::Vector, f::Field)
|
||||
f.values*x
|
||||
end
|
||||
Base.(:*)(x::Union{Vector, Matrix}, f::Field) = interpolate(x, f)
|
||||
|
||||
""" Interpolate field (h*f)(ξ) """
|
||||
Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
|
||||
|
||||
@@ -2,13 +2,7 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using FactCheck
|
||||
using JuliaFEM: test_element
|
||||
|
||||
facts("test set and interpolate field variable") do
|
||||
el = Quad4(1, [1, 2, 3, 4])
|
||||
set_coordinates(el, [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')
|
||||
set_field(el, "displacement", [0.0 0.0; 0.0 0.0; 0.5 0.0; 0.0 0.0]'')
|
||||
fval = interpolate(el, "displacement", [0.0, 1.0])
|
||||
Logging.debug(fval)
|
||||
@fact fval --> roughly([0.25 0.0]')
|
||||
end
|
||||
|
||||
using JuliaFEM: Quad4
|
||||
test_element(Quad4)
|
||||
|
||||
@@ -57,4 +57,8 @@ facts("test fields and interpolation") do
|
||||
@fact diff(h)([0.0, 0.0])*X --> [0.5 0.0; 0.0 0.5]
|
||||
@fact (diff(h)*X)([0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
|
||||
|
||||
b = Basis((xi) -> [1/2*(1-xi[1]), 1/2*(1+xi[1])])
|
||||
f = Field(0.0, 100.0)
|
||||
@fact b(0.0) * f --> [50.0, 50.0]
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user