mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-12 14:23:20 +00:00
lagrange elements can now easily be created using simple macro
This commit is contained in:
Jukka Aho
@@ -21,7 +21,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
|
||||
"Logger(root,DEBUG,PipeEndpoint(open, 0 bytes waiting),root)"
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
@@ -64,7 +64,18 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Creating element Point1: 1 node point element\n",
|
||||
"Creating element Seg2: 2 node linear line element\n",
|
||||
"Creating Lagrange basis for element Seg2. Number of basis functions: 2. Element dimension: 1\n",
|
||||
"Calculating inverse of A\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: Element"
|
||||
]
|
||||
@@ -82,9 +93,29 @@
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Element JuliaFEM.Seg2 created.\n",
|
||||
"Creating element Seg3: 3 node quadratic line element\n",
|
||||
"Creating Lagrange basis for element Seg3. Number of basis functions: 3. Element dimension: 1\n",
|
||||
"Calculating inverse of A\n",
|
||||
"Element JuliaFEM.Seg3 created.\n",
|
||||
"Creating element Quad4: 4 node bilinear quadrangle element\n",
|
||||
"Creating Lagrange basis for element Quad4. Number of basis functions: 4. Element dimension: 2\n",
|
||||
"Calculating inverse of A\n",
|
||||
"Element JuliaFEM.Quad4 created.\n",
|
||||
"Creating element Tet10: 10 node quadratic tetrahedron\n",
|
||||
"Creating Lagrange basis for element Tet10. Number of basis functions: 10. Element dimension: 3\n",
|
||||
"Calculating inverse of A\n",
|
||||
"Element JuliaFEM.Tet10 created.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"type Quad4 <: Element\n",
|
||||
"type MyQuad4 <: Element\n",
|
||||
" connectivity :: Array{Int, 1}\n",
|
||||
" fields :: Dict{Any, Any}\n",
|
||||
"end"
|
||||
@@ -107,7 +138,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Quad4"
|
||||
"MyQuad4"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
@@ -116,14 +147,14 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"Quad4(connectivity) = Quad4(connectivity, Dict{ASCIIString, Any}())"
|
||||
"MyQuad4(connectivity) = MyQuad4(connectivity, Dict{Any, Any}())"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Some basic charasteristics like number of nodes / connectivity points and dimension:"
|
||||
"Some basic charasteristics like number of basis funcitons and dimension:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -136,7 +167,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_element_dimension (generic function with 2 methods)"
|
||||
"get_element_dimension (generic function with 6 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 5,
|
||||
@@ -145,8 +176,8 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"JuliaFEM.get_number_of_nodes(el::Type{Quad4}) = 4\n",
|
||||
"JuliaFEM.get_element_dimension(el::Quad4) = 2"
|
||||
"JuliaFEM.get_number_of_basis_functions(el::Type{MyQuad4}) = 4\n",
|
||||
"JuliaFEM.get_element_dimension(el::MyQuad4) = 2"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -166,7 +197,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_dbasisdxi (generic function with 2 methods)"
|
||||
"get_dbasisdxi (generic function with 6 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 6,
|
||||
@@ -175,14 +206,14 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"function JuliaFEM.get_basis(el::Quad4, xi::Array{Float64,1})\n",
|
||||
"function JuliaFEM.get_basis(el::MyQuad4, xi)\n",
|
||||
" [(1-xi[1])*(1-xi[2])/4\n",
|
||||
" (1+xi[1])*(1-xi[2])/4\n",
|
||||
" (1+xi[1])*(1+xi[2])/4\n",
|
||||
" (1-xi[1])*(1+xi[2])/4]\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function JuliaFEM.get_dbasisdxi(el::Quad4, xi::Array{Float64,1})\n",
|
||||
"function JuliaFEM.get_dbasisdxi(el::MyQuad4, xi)\n",
|
||||
" [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n",
|
||||
" (1-xi[2])/4.0 -(1+xi[1])/4.0\n",
|
||||
" (1+xi[2])/4.0 (1+xi[1])/4.0\n",
|
||||
@@ -208,26 +239,43 @@
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"31-Aug 20:05:30:INFO:root:number of connectivity points (nodes) in this element: 4\n",
|
||||
"31-Aug 20:05:31:INFO:root:Constructing element..\n",
|
||||
"31-Aug 20:05:31:INFO:root:Element dimension: 2\n",
|
||||
"31-Aug 20:05:31:INFO:root:Setting scalar field [1 2 3 4] to element.\n",
|
||||
"31-Aug 20:05:31:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
|
||||
"31-Aug 20:05:31:INFO:root:Value: [2.5]\n",
|
||||
"31-Aug 20:05:31:INFO:root:Element Quad4 passed tests.\n"
|
||||
"10-Sep 20:16:49:INFO:root:Testing element MyQuad4\n",
|
||||
"10-Sep 20:16:49:INFO:root:number of basis functions in this element: 4\n",
|
||||
"10-Sep 20:16:49:INFO:root:Initializing element\n",
|
||||
"10-Sep 20:16:49:INFO:root:Element dimension: 2\n",
|
||||
"10-Sep 20:16:49:INFO:root:Setting scalar field [1 2 3 4] to element.\n",
|
||||
"10-Sep 20:16:50:INFO:root:Interpolating scalar field at [0.0,0.0]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"PipeEndpoint(open, 0 bytes waiting)"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"10-Sep 20:16:50:INFO:root:Value: [2.5]\n",
|
||||
"10-Sep 20:16:50:INFO:root:Element MyQuad4 passed tests.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: test_element\n",
|
||||
"test_element(Quad4)"
|
||||
"test_element(MyQuad4)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"If test_element passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:"
|
||||
"If `test_element` passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -239,9 +287,9 @@
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM: set_field, interpolate\n",
|
||||
"el1 = Quad4([1, 2, 3, 4])\n",
|
||||
"el1 = MyQuad4([1, 2, 3, 4])\n",
|
||||
"set_field(el1, :temperature, [1 2 3 4])\n",
|
||||
"set_field(el1, :coordinates, [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, :geometry, [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);"
|
||||
]
|
||||
},
|
||||
@@ -264,7 +312,8 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"interpolate(el1, :temperature, [0.0, 0.0]) # temperature at the middle poinf of the element, 1/4*(1+2+3+4)"
|
||||
"# temperature at the middle poinf of the element, 1/4*(1+2+3+4)\n",
|
||||
"interpolate(el1, :temperature, [0.0, 0.0])"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -289,12 +338,13 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"interpolate(el1, :coordinates, [0.0, 0.0]) # midpoint of element"
|
||||
"# geometry midpoint of element\n",
|
||||
"interpolate(el1, :geometry, [0.0, 0.0])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"execution_count": 12,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -305,12 +355,13 @@
|
||||
"1"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# interpolating scalar -> scalar.\n",
|
||||
"interpolate(el1, :\"heat coefficient\", [0.0, 0.0])"
|
||||
]
|
||||
},
|
||||
@@ -347,13 +398,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"execution_count": 15,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM: Equation, IntegrationPoint\n",
|
||||
"using JuliaFEM: Equation, IntegrationPoint, Quad4\n",
|
||||
"\n",
|
||||
"abstract Heat <: Equation"
|
||||
]
|
||||
@@ -367,7 +418,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"execution_count": 16,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -392,7 +443,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 14,
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -403,7 +454,7 @@
|
||||
"DC2D4"
|
||||
]
|
||||
},
|
||||
"execution_count": 14,
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -429,7 +480,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 15,
|
||||
"execution_count": 18,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -437,10 +488,10 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_lhs (generic function with 2 methods)"
|
||||
"has_lhs (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 15,
|
||||
"execution_count": 18,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -456,7 +507,8 @@
|
||||
" dNdX = get_dbasisdX(el, ip.xi)\n",
|
||||
" hc = interpolate(el, :\"temperature thermal conductivity\", ip.xi)\n",
|
||||
" return dNdX*hc*dNdX'\n",
|
||||
"end"
|
||||
"end\n",
|
||||
"JuliaFEM.has_lhs(eq::DC2D4) = true"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -468,7 +520,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"execution_count": 19,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -483,7 +535,7 @@
|
||||
" -1.0 -2.0 -1.0 4.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
@@ -491,12 +543,41 @@
|
||||
"source": [
|
||||
"using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
|
||||
"el = Quad4([1, 2, 3, 4])\n",
|
||||
"set_field(el, :coordinates, [0 0; 1 0; 1 1; 0 1]')\n",
|
||||
"set_field(el, :geometry, [0 0; 1 0; 1 1; 0 1]')\n",
|
||||
"set_field(el, :\"temperature thermal conductivity\", 6)\n",
|
||||
"eq = DC2D4(el)\n",
|
||||
"integrate_lhs(eq)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"If rhs or lhs is not defined, integration returns nothing."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 20,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"true"
|
||||
]
|
||||
},
|
||||
"execution_count": 20,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"isa(integrate_rhs(eq), Void)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
@@ -506,41 +587,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"31-Aug 20:05:35:INFO:root:number of connectivity points (nodes) in this element: 2\n",
|
||||
"31-Aug 20:05:35:INFO:root:Constructing element..\n",
|
||||
"31-Aug 20:05:35:INFO:root:Element dimension: 1\n",
|
||||
"31-Aug 20:05:35:INFO:root:Setting scalar field [1 2] to element.\n",
|
||||
"31-Aug 20:05:35:INFO:root:Interpolating scalar field at [0.0]\n",
|
||||
"31-Aug 20:05:35:INFO:root:Value: [1.5]\n",
|
||||
"31-Aug 20:05:35:INFO:root:Element Seg2 passed tests.\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"type Seg2 <: Element\n",
|
||||
" connectivity :: Array{Int, 1}\n",
|
||||
" fields :: Dict{Any, Any}\n",
|
||||
"end\n",
|
||||
"Seg2(connectivity) = Seg2(connectivity, Dict{ASCIIString, Any}())\n",
|
||||
"JuliaFEM.get_number_of_nodes(el::Type{Seg2}) = 2\n",
|
||||
"JuliaFEM.get_element_dimension(el::Seg2) = 1\n",
|
||||
"JuliaFEM.get_basis(el::Seg2, xi::Array{Float64,1}) = [0.5*(1-xi[1]), 0.5*(1+xi[1])]\n",
|
||||
"JuliaFEM.get_dbasisdxi(el::Seg2, xi::Array{Float64,1}) = [-0.5 0.5]'\n",
|
||||
"test_element(Seg2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 29,
|
||||
"execution_count": 21,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -548,19 +595,19 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"get_rhs (generic function with 2 methods)"
|
||||
"has_rhs (generic function with 2 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 29,
|
||||
"execution_count": 21,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using JuliaFEM: get_basis\n",
|
||||
"using JuliaFEM: get_basis, Seg2\n",
|
||||
"\n",
|
||||
"\"\"\"\n",
|
||||
"Diffusive heat transfer for 4-node bilinear element.\n",
|
||||
"Diffusive heat transfer for 2-node linear segment.\n",
|
||||
"\"\"\"\n",
|
||||
"type DC2D2 <: Heat\n",
|
||||
" element :: Seg2\n",
|
||||
@@ -583,12 +630,13 @@
|
||||
" N = get_basis(el, ip.xi)\n",
|
||||
" f = interpolate(el, :\"temperature flux\", ip.xi)\n",
|
||||
" return f*N\n",
|
||||
"end"
|
||||
"end\n",
|
||||
"JuliaFEM.has_rhs(eq::DC2D2) = true"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 30,
|
||||
"execution_count": 22,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -601,19 +649,30 @@
|
||||
" 50.0"
|
||||
]
|
||||
},
|
||||
"execution_count": 30,
|
||||
"execution_count": 22,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"el = Seg2([1, 2])\n",
|
||||
"set_field(el, :coordinates, [0.0 0.0; 0.0 1.0]')\n",
|
||||
"set_field(el, :geometry, [0.0 0.0; 0.0 1.0]')\n",
|
||||
"set_field(el, :\"temperature flux\", 100.0)\n",
|
||||
"eq = DC2D2(el)\n",
|
||||
"integrate_rhs(eq)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Defining own problem\n",
|
||||
"\n",
|
||||
"- takes a set of elements and maps corresponding equations for them\n",
|
||||
"- boundary conditions\n",
|
||||
"- solve!(problem) updates fields"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
@@ -621,18 +680,95 @@
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
"source": [
|
||||
"abstract Problem\n",
|
||||
"\n",
|
||||
"type HeatProblem <: Problem\n",
|
||||
" elements: Array{Any, 1}\n",
|
||||
" equations: Array{Any, 1}\n",
|
||||
" boundary_conditions: Array{Any, 1}\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Developing boundary conditions"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"abstract BC\n",
|
||||
"\n",
|
||||
"type DefaultBC <: BC\n",
|
||||
" \n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"bc = NodalBC()\n",
|
||||
"bc[:displacement, [1, 2, 3], 1:2] = 3.0"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"bc = WeakBC()\n",
|
||||
"bc[:displacement, [el1, el2, el3]] = (X) -> X[0] + X[1] - 2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"bc = MPCBC()\n",
|
||||
"bc[:displacement, "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"bc = MortarBC()\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 0.4.0-dev",
|
||||
"display_name": "Julia 0.5.0-dev",
|
||||
"language": "julia",
|
||||
"name": "julia-0.4"
|
||||
"name": "julia-0.5"
|
||||
},
|
||||
"language_info": {
|
||||
"name": "julia",
|
||||
"version": "0.4.0"
|
||||
"version": "0.5.0"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
|
||||
+137
-179
@@ -24,10 +24,10 @@ 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_number_of_basis_functions*
|
||||
- get_element_dimension *
|
||||
- get_basis *
|
||||
- get_dbasisdxi
|
||||
- get_dbasisdxi*
|
||||
- get_dbasisdX
|
||||
- get_field
|
||||
- set_field
|
||||
@@ -39,33 +39,41 @@ Which should work if element is defined following some rules. Functions marked w
|
||||
=#
|
||||
|
||||
# These must be implemented for your own element
|
||||
get_number_of_basis_functions(el::Type{Element}) = nothing
|
||||
get_element_dimension(el::Element) = nothing
|
||||
get_basis(el::Element, xi) = nothing
|
||||
""" Return partial derivatives of basis using ForwardDiff """
|
||||
function get_dbasisdxi(el::Element, xi)
|
||||
f(xi) = get_basis(el, xi)
|
||||
ForwardDiff.jacobian(f, xi)
|
||||
end
|
||||
function get_number_of_basis_functions(el::Type{Element})
|
||||
Logging.info("You really should define get_number_of_basis_functions")
|
||||
# this is hack, evaluate basis in some point and return length of vector.
|
||||
bf = get_basis([0.0, 0.0, 0.0])
|
||||
length(bf)
|
||||
end
|
||||
|
||||
abstract CG <: Element # Lagrange (continous Galerkin) element family
|
||||
get_dbasisdxi(el::Element, xi) = nothing
|
||||
|
||||
"""
|
||||
4 node bilinear quadrangle element
|
||||
Create new element with element_name to family element_family
|
||||
|
||||
X = [-1.0 1.0 1.0 -1.0
|
||||
-1.0 -1.0 1.0 1.0]
|
||||
|
||||
P(xi) = [1.0, xi[1], xi[2], xi[1]*xi[2]]
|
||||
|
||||
dP(xi) = [0.0 1.0 0.0 xi[2]
|
||||
0.0 0.0 1.0 xi[1]]
|
||||
Examples
|
||||
--------
|
||||
>>> @create_element(Seg2, CG, "2 node linear segment")
|
||||
"""
|
||||
macro create_element(element_name, element_family, element_description)
|
||||
print("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
|
||||
|
||||
#=
|
||||
|
||||
Start of example
|
||||
|
||||
Example how to create new element. This is commented because I use code
|
||||
generation for simple elements like Lagrage elements. Feel free to use
|
||||
code generation but elements can be of course created manually too!
|
||||
|
||||
type Quad4 <: CG
|
||||
connectivity :: Array{Int, 1}
|
||||
fields :: Dict{Any, Any}
|
||||
@@ -75,7 +83,7 @@ end
|
||||
Quad4(connectivity) = Quad4(connectivity, Dict{Any, Any}())
|
||||
|
||||
""" Return number of basis functions of this element. """
|
||||
get_number_of_basis_functions(el::Type{Quad4})
|
||||
get_number_of_basis_functions(el::Type{Quad4}) = 4
|
||||
|
||||
""" Return element dimension (length of xi vector). """
|
||||
get_element_dimension(el::Type{Quad4}) = 2
|
||||
@@ -95,10 +103,93 @@ function get_dbasisdxi(el::Quad4, xi)
|
||||
(1+xi[2])/4.0 (1+xi[1])/4.0
|
||||
-(1+xi[2])/4.0 (1-xi[1])/4.0]
|
||||
end
|
||||
# TODO: create_lagrange_element(:Quad4, X, P, dP)
|
||||
|
||||
End of example.
|
||||
|
||||
=#
|
||||
|
||||
abstract CG <: Element # Lagrange (continous Galerkin) element family
|
||||
|
||||
"""
|
||||
Given polynomial P and coordinates of reference element, calculate
|
||||
Lagrange basis function and partial derivatives.
|
||||
"""
|
||||
function calculate_lagrange_basis(P, X)
|
||||
dim, nbasis = size(X)
|
||||
A = zeros(nbasis, nbasis)
|
||||
for i=1:nbasis
|
||||
A[i,:] = P(X[:, i])
|
||||
end
|
||||
println("Calculating inverse of A")
|
||||
invA = inv(A)'
|
||||
basis(xi) = invA*P(xi)
|
||||
dbasisdxi = ForwardDiff.jacobian(basis)
|
||||
basis, dbasisdxi
|
||||
end
|
||||
|
||||
"""
|
||||
Assign Lagrange basis for element.
|
||||
"""
|
||||
macro create_lagrange_basis(element_name, X, P)
|
||||
|
||||
print("Creating Lagrange basis for element ", element_name, ". ")
|
||||
eltype = esc(element_name)
|
||||
|
||||
quote
|
||||
|
||||
global get_number_of_basis_functions, get_element_dimension
|
||||
global get_basis, get_dbasisdxi
|
||||
|
||||
dim = size($X, 1)
|
||||
nbasis = size($X, 2)
|
||||
print("Number of basis functions: ", nbasis, ". ")
|
||||
println("Element dimension: ", dim)
|
||||
|
||||
get_number_of_basis_functions(el::Type{$(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)
|
||||
println("Element ", $element_name, " created.")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# 0d Lagrange 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_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])
|
||||
|
||||
# 2d Lagrange elements
|
||||
|
||||
@create_element(Quad4, CG, "4 node bilinear quadrangle element")
|
||||
@create_lagrange_basis(Quad4,
|
||||
[-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,
|
||||
[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],
|
||||
(xi) -> [ 1.0, xi[1], xi[2], xi[3], xi[1]^2,
|
||||
xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
|
||||
|
||||
# Common element routines
|
||||
|
||||
|
||||
"""
|
||||
Test routine for element.
|
||||
|
||||
@@ -112,25 +203,30 @@ Raises
|
||||
This uses FactCheck and throws exception if element is not passing.
|
||||
"""
|
||||
function test_element(eltype)
|
||||
Logging.info("Testing element $eltype")
|
||||
local el
|
||||
n = get_number_of_nodes(eltype)
|
||||
Logging.info("number of connectivity points (nodes) in this element: $n")
|
||||
@fact n --> not(-1) """Unable to determine number of nodes for $eltype
|
||||
define a function 'get_number_of_nodes' which returns the number of nodes for this element."""
|
||||
n = get_number_of_basis_functions(eltype)
|
||||
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
|
||||
'get_number_of_basis_functions' which returns the number of nodes
|
||||
for this element."""
|
||||
|
||||
Logging.info("Constructing element..")
|
||||
Logging.info("Initializing element")
|
||||
try
|
||||
el = eltype(collect(1:n))
|
||||
catch
|
||||
Logging.error("""Unable to create element with default constructor
|
||||
define function $eltype(connectivity) which initializes this element.
|
||||
""")
|
||||
Logging.error("""
|
||||
Unable to create element with default constructor define function
|
||||
$eltype(connectivity) which initializes this element.""")
|
||||
return false
|
||||
end
|
||||
|
||||
dim = get_element_dimension(el)
|
||||
Logging.info("Element dimension: $dim")
|
||||
@fact dim --> not(-1) """Unable to get element dimension
|
||||
define function 'get_element_dimension' which return the dimension of this element (1, 2, 3)"""
|
||||
@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 = collect(1:n)'
|
||||
@@ -141,14 +237,16 @@ function test_element(eltype)
|
||||
try
|
||||
get_basis(el, zeros(dim))
|
||||
catch
|
||||
Logging.error("""Unable to evaluate basis, define function 'get_basis' for this element.
|
||||
""")
|
||||
Logging.error("""
|
||||
Unable to evaluate basis, define function 'get_basis' for
|
||||
this element.""")
|
||||
end
|
||||
try
|
||||
get_dbasisdxi(el, zeros(dim))
|
||||
catch
|
||||
Logging.error("""Unable to evaluate partial derivatives of basis, define function 'get_dbasisdxi' for this element.
|
||||
""")
|
||||
Logging.error("""
|
||||
Unable to evaluate partial derivatives of basis,
|
||||
define function 'get_dbasisdxi' for this element.""")
|
||||
end
|
||||
|
||||
xi = zeros(dim)
|
||||
@@ -207,143 +305,3 @@ function interpolate(el::Element, field, xi)
|
||||
return result
|
||||
end
|
||||
end
|
||||
|
||||
#=
|
||||
"""
|
||||
Create new Lagrange element
|
||||
|
||||
FIXME: this is not working
|
||||
|
||||
LoadError: error compiling anonymous: type definition not allowed inside a local scope
|
||||
|
||||
It's the for loop which is causing problems. See
|
||||
https://github.com/JuliaLang/julia/issues/10555
|
||||
|
||||
"""
|
||||
function create_lagrange_element(element_name, X, P, dP)
|
||||
|
||||
@eval begin
|
||||
|
||||
nnodes, dim = size(X)
|
||||
A = zeros(nnodes, nnodes)
|
||||
for i=1:nnodes
|
||||
A[i,:] = P(X[i,:])
|
||||
end
|
||||
invA = inv(A)'
|
||||
|
||||
type $element_name
|
||||
node_ids :: Array{Int, 1}
|
||||
fields :: Dict{ASCIIString, Any}
|
||||
end
|
||||
|
||||
function $element_name(node_ids)
|
||||
fields = Dict{ASCIIString, Any}()
|
||||
$element_name(node_ids, fields)
|
||||
end
|
||||
|
||||
function get_basis(el::$element_name, xi)
|
||||
invA*P(xi)
|
||||
end
|
||||
|
||||
function get_dbasisdxi(el::$element_name, xi)
|
||||
invA*dP(xi)
|
||||
end
|
||||
|
||||
$element_name
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
=#
|
||||
|
||||
# 0d Lagrange elements
|
||||
|
||||
#=
|
||||
"""
|
||||
1 node point element
|
||||
"""
|
||||
type Point1 <: CG
|
||||
node_ids :: Array{Int, 1}
|
||||
fields :: Dict{ASCIIString, Any}
|
||||
end
|
||||
function Point1(node_ids)
|
||||
fields = Dict{ASCIIString, Any}()
|
||||
Point1(node_ids, fields)
|
||||
end
|
||||
|
||||
# 1d Lagrange elements
|
||||
|
||||
"""
|
||||
2 node linear line element
|
||||
"""
|
||||
type Seg2 <: CG
|
||||
node_ids :: Array{Int, 1}
|
||||
fields :: Dict{ASCIIString, Any}
|
||||
end
|
||||
|
||||
# X = [-1.0 1.0]'
|
||||
# P = (xi) -> [1.0 xi[1]]'
|
||||
# dP = (xi) -> [0.0 1.0]'
|
||||
# create_lagrange_element(:Seg2, X, P, dP)
|
||||
|
||||
"""
|
||||
3 node quadratic line element
|
||||
"""
|
||||
type Seg2 <: CG
|
||||
node_ids :: Array{Int, 1}
|
||||
fields :: Dict{ASCIIString, Any}
|
||||
end
|
||||
#X = [-1.0 1.0 0.0]'
|
||||
#P = (xi) -> [1.0 xi[1] xi[1]^2]'
|
||||
#dP = (xi) -> [0.0 1.0 2*xi[1]]'
|
||||
#create_lagrange_element(:Seg3, X, P, dP)
|
||||
|
||||
# 2d Lagrange elements
|
||||
=#
|
||||
|
||||
# 3d Lagrange elements
|
||||
#=
|
||||
"""
|
||||
10 node quadratic tethahedron
|
||||
"""
|
||||
type Tet10 <: CG
|
||||
node_ids :: Array{Int, 1}
|
||||
fields :: Dict{ASCIIString, Any}
|
||||
end
|
||||
=#
|
||||
# X = [
|
||||
# 0.0 0.0 0.0
|
||||
# 1.0 0.0 0.0
|
||||
# 0.0 1.0 0.0
|
||||
# 0.0 0.0 1.0
|
||||
# 0.5 0.0 0.0
|
||||
# 0.5 0.5 0.0
|
||||
# 0.0 0.5 0.0
|
||||
# 0.0 0.0 0.5
|
||||
# 0.5 0.0 0.5
|
||||
# 0.0 0.5 0.5]
|
||||
# P(xi) = [
|
||||
# 1
|
||||
# xi[1]
|
||||
# xi[2]
|
||||
# xi[3]
|
||||
# xi[1]^2
|
||||
# xi[2]^2
|
||||
# xi[3]^2
|
||||
# xi[1]*xi[2]
|
||||
# xi[2]*xi[3]
|
||||
# xi[3]*xi[1]]
|
||||
# dP(xi) = [
|
||||
# 0 0 0
|
||||
# 1 0 0
|
||||
# 0 1 0
|
||||
# 0 0 1
|
||||
# 2*xi[1] 0 0
|
||||
# 0 2*xi[2] 0
|
||||
# 0 0 2*xi[3]
|
||||
# xi[2] xi[1] 0
|
||||
# 0 xi[3] xi[2]
|
||||
# xi[3] 0 xi[1]
|
||||
# ]
|
||||
#create_lagrange_element(:Tet10, X, P, dP)
|
||||
|
||||
Reference in New Issue
Block a user