Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-09 04:36:34 +00:00
Code Issues Packages Projects Releases Wiki Activity

updated.

Browse Source
This commit is contained in:
Jukka Aho
2015-08-19 18:50:43 +03:00
parent 584269d33a
commit 1fa7f259fd
4 changed files with 354 additions and 370 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/2015-06-25-elasticity-solver-example.ipynb
+319 -348
View File
@@ -36,7 +36,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 1,
"metadata": {
"collapsed": false
},
@@ -47,9 +47,30 @@
"Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
]
},
"execution_count": 2,
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in int32 at deprecated.jl:49\n",
" in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
" in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
" in anonymous at task.jl:365\n",
"while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n",
"WARNING: int32(x) is deprecated, use Int32(x) instead.\n",
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in int32 at deprecated.jl:49\n",
" in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
" in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
" in anonymous at task.jl:365\n",
"while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n"
]
}
],
"source": [
@@ -73,7 +94,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"metadata": {
"collapsed": false
},
@@ -84,39 +105,35 @@
"jacobian (generic function with 1 method)"
]
},
"execution_count": 3,
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n",
" # Calculate dN/dX\n",
" dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n",
"\n",
" # Calculate residual force vector R(u) = T(u) - F(u)\n",
" function T(u)\n",
" # kinematics\n",
" gradu = u*dbasisdX\n",
" F = I + gradu\n",
" E = 1/2*(gradu' + gradu + gradu'*gradu)\n",
" # constitutive equation\n",
" lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n",
" mu = JuliaFEM.interpolate(el, \"mu\", xi)\n",
" S = lambda*trace(E)*I + 2*mu*E\n",
" P = F*S\n",
" T = P*dbasisdX'\n",
" return T\n",
" end\n",
"\n",
"function Wint_integrand(el, ip)\n",
" xi = ip.xi\n",
" dbasisdX = JuliaFEM.get_dbasisdX(el, ip)\n",
" u = el.attributes[\"displacement\"]\n",
" R(u) = T(u)\n",
" return R(u) # we skip calculating external force vector for now.\n",
"\n",
" # kinematics\n",
" gradu = u*dbasisdX\n",
" F = I + gradu\n",
" E = 1/2*(gradu' + gradu + gradu'*gradu)\n",
"\n",
" # constitutive equation\n",
" lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n",
" mu = JuliaFEM.interpolate(el, \"mu\", xi)\n",
" S = lambda*trace(E)*I + 2*mu*E\n",
" P = F*S\n",
" T = P*dbasisdX'\n",
" return T\n",
"end\n",
"\n",
"calc_residual_vector = JuliaFEM.integrate(calc_residual_vector_integrand)\n",
"Wint = JuliaFEM.integrate(Wint_integrand)\n",
"\n",
"# calculate partial derivatives of R with respect to field \"displacement\"\n",
"calc_tangent_stiffness = JuliaFEM.linearize(calc_residual_vector, \"displacement\")"
"jacobian = JuliaFEM.linearize(Wint, \"displacement\")"
]
},
{
@@ -130,7 +147,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 3,
"metadata": {
"collapsed": true
},
@@ -141,7 +158,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 4,
"metadata": {
"collapsed": false,
"scrolled": false
@@ -158,16 +175,16 @@
"name": "stderr",
"output_type": "stream",
"text": [
"17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n",
"17-Aug 18:21:25:DEBUG:root:solution vector: \n",
"19-Aug 18:39:44:DEBUG:root:Converged in 6 iterations.\n",
"19-Aug 18:39:44:DEBUG:root:solution vector: \n",
" [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n",
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
"17-Aug 18:21:25:DEBUG:root:norm of u: 3.1292483947150047\n",
"17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n",
"17-Aug 18:21:25:DEBUG:root:solution vector: \n",
"19-Aug 18:39:44:DEBUG:root:norm of u: 3.1292483947150047\n",
"19-Aug 18:39:45:DEBUG:root:Converged in 6 iterations.\n",
"19-Aug 18:39:45:DEBUG:root:solution vector: \n",
" [0.0 0.7433248532717796 1.048521014723486 0.0\n",
" 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n",
"17-Aug 18:21:26:DEBUG:root:norm of u: 3.1292483947150056\n"
"19-Aug 18:39:45:DEBUG:root:norm of u: 3.1292483947150056\n"
]
},
{
@@ -183,7 +200,7 @@
"delayed_handler (generic function with 4 methods)"
]
},
"execution_count": 6,
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
@@ -200,10 +217,15 @@
" (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",
" ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]'\n",
" iweights = [1.0, 1.0, 1.0, 1.0]\n",
" integration_points = [\n",
" JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n",
" JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n",
" JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0, Dict()),\n",
" JuliaFEM.IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0, Dict())]\n",
" attributes = Dict()\n",
" e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\n",
" element_id = 1\n",
" node_ids = [1, 2, 3, 4]\n",
" e = JuliaFEM.Element(element_id, node_ids, basis, dbasis, integration_points, attributes)\n",
" E = 90.0\n",
" nu = 0.25\n",
" mu = E/(2*(1+nu))\n",
@@ -225,8 +247,8 @@
"\n",
" free_dofs = [3, 4, 5, 6]\n",
" for i=1:10\n",
" R = calc_residual_vector(e)\n",
" K = calc_tangent_stiffness(e)\n",
" R = Wint(e)\n",
" K = jacobian(e)\n",
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
" e.attributes[\"displacement\"] += du\n",
" if norm(du) < 1.0e-9\n",
@@ -253,8 +275,8 @@
" e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n",
" du = zeros(2, 4)\n",
" for i=1:10\n",
" R = calc_residual_vector(e)\n",
" K = calc_tangent_stiffness(e)\n",
" R = Wint(e)\n",
" K = jacobian(e)\n",
" du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
" e.attributes[\"displacement\"] += du\n",
" if norm(du) < 1.0e-9\n",
@@ -269,16 +291,9 @@
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One element solutions are not particularly interesting so next step is to create function that assembles global matrix from local matrices. Some data types:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 5,
"metadata": {
"collapsed": false
},
@@ -286,35 +301,21 @@
{
"data": {
"text/plain": [
"assemble_element! (generic function with 2 methods)"
"assemble_element! (generic function with 1 method)"
]
},
"execution_count": 9,
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function assemble_element!(ass::Assembly, el::Element)\n",
"\n",
" # Material properties\n",
" E = el.attributes[\"Young\"]\n",
" nu = el.attributes[\"Poisson\"]\n",
" mu = E/(2*(1+nu))\n",
" la = E*nu/((1+nu)*(1-2*nu))\n",
" la = 2*la*mu/(la + 2*mu)\n",
"\n",
" dofs = prod(size(el.coordinates))\n",
" X = el.coordinates\n",
" u = el.attributes[\"displacement\"]\n",
" R = el.attributes[\"displacement nodal force\"]\n",
" K = el.attributes[\"displacement tangent stiffness\"]\n",
"function assemble_element!(ass::JuliaFEM.Assembly, el::JuliaFEM.Element)\n",
"\n",
" gdofs = ass.gdofs[el.id]\n",
" #Logging.debug(\"Assemble element $(el.id) to gdofs $gdofs\")\n",
" basis, dbasis = get_shape_functions(el)\n",
" ipoints, iweights = get_integration_scheme(el, io)\n",
" calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
" R = Wint(el)\n",
" K = jacobian(el)\n",
" dofs = length(R)\n",
"\n",
" for i=1:dofs\n",
" for j=1:dofs\n",
@@ -342,18 +343,11 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: the `=>` syntax is deprecated, use `-->` instead\n"
]
},
{
"name": "stdout",
"output_type": "stream",
@@ -365,30 +359,29 @@
"name": "stderr",
"output_type": "stream",
"text": [
"15-Aug 17:03:55:DEBUG:root:Adding nodes to array\n",
"15-Aug 17:03:55:DEBUG:root:Creating elements\n",
"15-Aug 17:03:55:DEBUG:root:Starting iteration 1\n",
"15-Aug 17:03:55:DEBUG:root:Assembling\n",
"15-Aug 17:03:57:DEBUG:root:Solution norm = 3.090022136728999\n",
"15-Aug 17:03:57:DEBUG:root:Starting iteration 2\n",
"15-Aug 17:03:57:DEBUG:root:Assembling\n",
"15-Aug 17:03:57:DEBUG:root:Solution norm = 0.32121316021535135\n",
"15-Aug 17:03:57:DEBUG:root:Starting iteration 3\n",
"15-Aug 17:03:57:DEBUG:root:Assembling\n",
"15-Aug 17:03:57:DEBUG:root:Solution norm = 0.040431781939994194\n",
"15-Aug 17:03:57:DEBUG:root:Starting iteration 4\n",
"15-Aug 17:03:57:DEBUG:root:Assembling\n",
"15-Aug 17:03:57:DEBUG:root:Solution norm = 0.0009291101052124042\n",
"15-Aug 17:03:57:DEBUG:root:Starting iteration 5\n",
"15-Aug 17:03:57:DEBUG:root:Assembling\n",
"15-Aug 17:03:57:DEBUG:root:Solution norm = 1.5638899022213175e-7\n",
"15-Aug 17:03:58:DEBUG:root:Starting iteration 6\n",
"15-Aug 17:03:58:DEBUG:root:Assembling\n",
"15-Aug 17:03:58:DEBUG:root:Solution norm = 1.0118539067290854e-14\n",
"15-Aug 17:03:58:DEBUG:root:Converged in 6 iterations.\n",
"15-Aug 17:03:58:DEBUG:root:Displacement of element = \n",
"[-0.39914506095474334 -0.0722858269559246 0.0 0.0\n",
" -2.1779892317073504 -2.2222447544017645 0.0 0.0]\n"
"19-Aug 18:39:46:DEBUG:root:Adding nodes to array\n",
"19-Aug 18:39:46:DEBUG:root:Starting iteration 1\n",
"19-Aug 18:39:46:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 3.0900221367289986\n",
"19-Aug 18:39:47:DEBUG:root:Starting iteration 2\n",
"19-Aug 18:39:47:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 0.3212131602153472\n",
"19-Aug 18:39:47:DEBUG:root:Starting iteration 3\n",
"19-Aug 18:39:47:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 0.04043178193999703\n",
"19-Aug 18:39:47:DEBUG:root:Starting iteration 4\n",
"19-Aug 18:39:47:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 0.0009291101052105739\n",
"19-Aug 18:39:47:DEBUG:root:Starting iteration 5\n",
"19-Aug 18:39:47:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 1.5638899027804743e-7\n",
"19-Aug 18:39:47:DEBUG:root:Starting iteration 6\n",
"19-Aug 18:39:47:DEBUG:root:Assembling\n",
"19-Aug 18:39:47:DEBUG:root:Solution norm = 1.0464940956129567e-14\n",
"19-Aug 18:39:47:DEBUG:root:Converged in 6 iterations.\n",
"19-Aug 18:39:47:DEBUG:root:Displacement of element = \n",
"[0.0 -0.39914506095474334 -0.0722858269559246 0.0\n",
" 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n"
]
},
{
@@ -404,7 +397,7 @@
"delayed_handler (generic function with 4 methods)"
]
},
"execution_count": 10,
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
@@ -412,41 +405,26 @@
"source": [
"facts(\"one element assembly\") do\n",
" # Create model\n",
" #Logging.debug(\"Creating nodes\")\n",
" #n1 = Node(1)\n",
" #n2 = Node(2)\n",
" #n3 = Node(3)\n",
" #n4 = Node(4)\n",
" Logging.debug(\"Adding nodes to array\")\n",
" #nodes = [n1.id, n2.id, n3.id, n4.id]\n",
" node_ids = [1, 2, 3, 4]\n",
" coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n",
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
" Logging.debug(\"Creating elements\")\n",
" el = Element(1, node_ids, coordinates, attributes)\n",
"\n",
" # Initialize elements ready for solution\n",
" el.attributes[\"displacement\"] = zeros(2, 4)\n",
" el.attributes[\"displacement nodal force\"] = zeros(2, 4)\n",
" el.attributes[\"displacement tangent stiffness\"] = zeros(8, 8)\n",
"\n",
" el = get_test_element()\n",
" \n",
" for i=1:10\n",
" Logging.debug(\"Starting iteration $i\")\n",
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())\n",
" ass = JuliaFEM.Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], Dict{Int64,Array{Int64,1}}())\n",
" ass.gdofs[el.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n",
" Logging.debug(\"Assembling\")\n",
" assemble_element!(ass, el)\n",
"\n",
" # Boundary conditions\n",
" F = [0 0; 0 -2; 0 0; 0 0]'\n",
" F = reshape(F, prod(size(F)))\n",
" free_dofs = [1, 2, 3, 4]\n",
" F = [0.0 0.0; 0.0 0.0; 0.0 -2.0; 0.0 0.0]'\n",
" F = F[:]\n",
" free_dofs = [3, 4, 5, 6]\n",
"\n",
" # solution\n",
" K = sparse(ass.I, ass.J, ass.A)\n",
" R = full(sparsevec(ass.i, ass.b))\n",
" R = R - F\n",
" du = zeros(8) # must be determined from ass\n",
" du = zeros(8)\n",
" du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n",
"\n",
" Logging.debug(\"Solution norm = $(norm(du))\")\n",
@@ -462,7 +440,7 @@
" end\n",
" disp = el.attributes[\"displacement\"]\n",
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
" @fact norm(disp) => roughly(3.1292483947150043)\n",
" @fact norm(disp) --> roughly(3.1292483947150043)\n",
"end"
]
},
@@ -479,7 +457,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 7,
"metadata": {
"collapsed": false
},
@@ -493,83 +471,18 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 8,
"metadata": {
"collapsed": false,
"scrolled": false
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solve two element problem\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"15-Aug 17:03:59:DEBUG:root:Creating elements\n",
"15-Aug 17:03:59:INFO:root:create_ldof2gdofmap: dofs per node: 2\n",
"15-Aug 17:04:00:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],5=>[9,10],6=>[11,12],1=>[1,2])\n",
"15-Aug 17:04:04:INFO:root:solve!: dofs per node: 2\n",
"15-Aug 17:04:04:DEBUG:root:Problem size = 12\n",
"15-Aug 17:04:04:DEBUG:root:Starting iteration 1\n",
"15-Aug 17:04:04:DEBUG:root:Assembling\n",
"15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
"15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n",
"15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.6015838690633517\n",
"15-Aug 17:04:05:DEBUG:root:Starting iteration 2\n",
"15-Aug 17:04:05:DEBUG:root:Assembling\n",
"15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
"15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n",
"15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.013420417380980414\n",
"15-Aug 17:04:05:DEBUG:root:Starting iteration 3\n",
"15-Aug 17:04:05:DEBUG:root:Assembling\n",
"15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
"15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
"15-Aug 17:04:06:DEBUG:root:Solution norm du = 0.00032202957936854873\n",
"15-Aug 17:04:06:DEBUG:root:Starting iteration 4\n",
"15-Aug 17:04:06:DEBUG:root:Assembling\n",
"15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n",
"15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
"15-Aug 17:04:06:DEBUG:root:Solution norm du = 9.906677094801476e-8\n",
"15-Aug 17:04:06:DEBUG:root:Starting iteration 5\n",
"15-Aug 17:04:06:DEBUG:root:Assembling\n",
"15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n",
"15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
"15-Aug 17:04:06:DEBUG:root:Solution norm du = 5.027624635820698e-15\n",
"15-Aug 17:04:06:DEBUG:root:Converged in 5 iterations.\n",
"15-Aug 17:04:06:DEBUG:root:Displacement of element = \n",
"[-0.02442313597467864 -0.039356000063335075 0.021097993207232584 0.020877031423993653\n",
" -0.12673626841485705 -0.40433021969759375 -0.40656320177872923 -0.1275776940913048]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"0 facts verified.\n"
]
},
{
"data": {
"text/plain": [
"delayed_handler (generic function with 4 methods)"
"create_ldof2gdofmap (generic function with 1 method)"
]
},
"execution_count": 12,
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
@@ -578,18 +491,17 @@
"\"\"\"\n",
"Create local dof to global dof mapping for given elements\n",
"\"\"\"\n",
"function create_ldof2gdofmap(elements; ndofs=2)\n",
" Logging.info(\"create_ldof2gdofmap: dofs per node: $ndofs\")\n",
"function create_ldof2gdofmap(elements, field)\n",
"\n",
" ndofs = size(elements[1].attributes[field], 1)\n",
" \n",
" all_node_ids = Int64[]\n",
" for el in elements\n",
" eldim, elnodes = size(el.coordinates)\n",
" for nid in el.node_ids\n",
" push!(all_node_ids, nid)\n",
" end\n",
" end\n",
" all_node_ids = unique(all_node_ids)\n",
" sort!(all_node_ids)\n",
"\n",
" # Assign global dof for each node\n",
" pdim = 1\n",
@@ -600,24 +512,104 @@
" end\n",
"\n",
" return ngdofs\n",
"end\n",
"\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"solve one element problem\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"19-Aug 18:49:12:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n",
"19-Aug 18:49:13:INFO:root:solve!: dofs per node: 2\n",
"19-Aug 18:49:13:DEBUG:root:Problem size = 8\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 1\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 3.0900221367289444\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 2\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.32121316021534796\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 3\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.040431781940014504\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 4\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 0.0009291101052065917\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 5\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n",
"19-Aug 18:49:13:DEBUG:root:Starting iteration 6\n",
"19-Aug 18:49:13:DEBUG:root:Assembling\n",
"19-Aug 18:49:13:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
"19-Aug 18:49:13:DEBUG:root:Added 5 Lagrange multipliers to model\n",
"19-Aug 18:49:13:DEBUG:root:Adding Neumann boundary conditions\n",
"19-Aug 18:49:13:DEBUG:root:Solving system of equations. Total size = 12\n",
"19-Aug 18:49:13:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n",
"19-Aug 18:49:13:DEBUG:root:Converged in 6 iterations.\n",
"19-Aug 18:49:13:DEBUG:root:Displacement of element = \n",
"[0.0 -0.39914506095474345 -0.07228582695592467 0.0\n",
" 0.0 -2.177989231707351 -2.222244754401765 0.0]\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"1 fact verified.\n"
]
},
{
"data": {
"text/plain": [
"delayed_handler (generic function with 4 methods)"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function solve!(elements, dofmap, neumann_bcs, dirichlet_bcs; ndofs=2, max_iterations=10)\n",
"\n",
" dbc = \"lagrange\"\n",
" \n",
" Logging.info(\"solve!: dofs per node: $ndofs\")\n",
" pdim = length(dofmap)*ndofs\n",
" Logging.debug(\"Problem size = $pdim\")\n",
"\n",
" # Initialize elements ready for solution\n",
" for el in elements\n",
" eldim, elnodes = size(el.coordinates)\n",
" el.attributes[\"displacement\"] = zeros(ndofs, elnodes)\n",
" el.attributes[\"displacement nodal force\"] = zeros(ndofs, elnodes)\n",
" el.attributes[\"displacement tangent stiffness\"] = zeros(ndofs*elnodes, ndofs*elnodes)\n",
" end\n",
"\n",
" # Assign global dofs for elements\n",
" gdofs = Dict{Int64, Array{Int64,1}}()\n",
" for el in elements\n",
@@ -629,7 +621,7 @@
" end\n",
" end\n",
"\n",
" ass = Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n",
" ass = JuliaFEM.Assembly(Int64[], Int64[], Float64[], Int64[], Float64[], gdofs)\n",
"\n",
" for iter=1:max_iterations\n",
" Logging.debug(\"Starting iteration $iter\")\n",
@@ -642,30 +634,26 @@
" Logging.debug(\"Assembling\")\n",
" for el in elements\n",
" assemble_element!(ass, el)\n",
" #println(\"Element stiffness matrix\")\n",
" #dump(round(el.attributes[\"displacement tangent stiffness\"], 2))\n",
" end\n",
"\n",
" i = 1\n",
" if dbc == \"lagrange\"\n",
" Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n",
" # Dirichlet boundary conditions\n",
" for bc in dirichlet_bcs\n",
" for (dof, val) in zip(bc.dofs, bc.values)\n",
" #Logging.debug(\"dof $dof => $val\")\n",
" push!(ass.I, dof)\n",
" push!(ass.J, pdim+i)\n",
" push!(ass.A, 1)\n",
" push!(ass.I, pdim+i)\n",
" push!(ass.J, dof)\n",
" push!(ass.A, 1)\n",
" push!(ass.i, pdim+i)\n",
" push!(ass.b, 0)\n",
" i += 1\n",
" end\n",
" Logging.debug(\"Adding Dirichlet boundary conditions using Lagrange multipliers\")\n",
" # Dirichlet boundary conditions\n",
" for bc in dirichlet_bcs\n",
" for (dof, val) in zip(bc.dofs, bc.values)\n",
" #Logging.debug(\"dof $dof => $val\")\n",
" push!(ass.I, dof)\n",
" push!(ass.J, pdim+i)\n",
" push!(ass.A, 1)\n",
" push!(ass.I, pdim+i)\n",
" push!(ass.J, dof)\n",
" push!(ass.A, 1)\n",
" push!(ass.i, pdim+i)\n",
" push!(ass.b, 0)\n",
" i += 1\n",
" end\n",
" Logging.debug(\"Added $i Lagrange multipliers to model\")\n",
" end\n",
" Logging.debug(\"Added $i Lagrange multipliers to model\")\n",
" i -= 1\n",
"\n",
" Logging.debug(\"Adding Neumann boundary conditions\")\n",
@@ -682,20 +670,8 @@
" K = sparse(ass.I, ass.J, ass.A)\n",
" R = full(sparsevec(ass.i, ass.b))\n",
" R = R - F\n",
" #Logging.debug(dump(round(full(K), 2)))\n",
" #print_matrix(full(K))\n",
" #Logging.debug(dump(round(R', 1)))\n",
"\n",
" if dbc == \"eliminate\"\n",
" Logging.debug(\"Eliminating Dirichlet boundary conditions\")\n",
" throw(\"Implement this properly\")\n",
" free_dofs = [1, 2, 3, 4, 5, 6, 7, 8]\n",
" du = zeros(12)\n",
" Logging.debug(\"K norm = $(norm(full(K[free_dofs, free_dofs])))\")\n",
" du[free_dofs] = K[free_dofs, free_dofs] \\ -R[free_dofs]\n",
" else\n",
" du = K \\ -R\n",
" end\n",
" du = K \\ -R\n",
"\n",
" #du = reshape(du, 2, 6)\n",
" solnorm = norm(du[1:pdim])\n",
@@ -725,55 +701,27 @@
"\n",
"ENV[\"COLUMNS\"] = 160\n",
"\n",
"function test1():\n",
" facts(\"solve one element problem\") do\n",
" # Create model\n",
" Logging.debug(\"Creating nodes\")\n",
" node_ids = [1, 2, 3, 4]\n",
" coordinates = [10.0 0.0; 10.0 1.0; 0.0 1.0; 0.0 0.0]'\n",
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
" Logging.debug(\"Creating elements\")\n",
" el = Element(1, node_ids, coordinates, attributes)\n",
" elements = [el]\n",
" dofmap = create_ldof2gdofmap(elements)\n",
" Logging.debug(dofmap)\n",
" # Boundary conditions\n",
" # here we want to create nodal force for third dof, that is, node id 2, second dof\n",
" bc1 = BC([dofmap[2][2]], [-2.0])\n",
" # dirichlet bc, set dx=dy=0 on support\n",
" bc2 = BC([dofmap[3][1], dofmap[3][2], dofmap[4][1], dofmap[4][2]], [0.0, 0.0, 0.0, 0.0])\n",
" solve!(elements, dofmap, [bc1], [bc2]; max_iterations=7)\n",
" disp = elements[1].attributes[\"displacement\"]\n",
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
" @fact norm(disp) => roughly(3.1292483947150043)\n",
" end\n",
"end\n",
"\n",
"facts(\"solve two element problem\") do\n",
"facts(\"solve one element problem\") do\n",
" # Create model\n",
" attributes = Dict(\"Young\" => 90, \"Poisson\" => 0.25)\n",
"\n",
" Logging.debug(\"Creating elements\")\n",
" nids1 = [5, 1, 3, 6]\n",
" coords1 = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'\n",
" el1 = Element(1, nids1, coords1, copy(attributes))\n",
" nids2 = [1, 2, 4, 3]\n",
" coords2 = [5.0 0.0; 10.0 0.0; 10.0 1.0; 5.0 1.0]'\n",
" el2 = Element(2, nids2, coords2, copy(attributes))\n",
" elements = [el1, el2]\n",
"\n",
" dofmap = create_ldof2gdofmap(elements)\n",
" el = get_test_element()\n",
" elements = [el]\n",
" # Initialize elements ready for solution\n",
" for el in elements\n",
" eldim, elnodes = size(el.attributes[\"coordinates\"])\n",
" el.attributes[\"displacement\"] = zeros(2, elnodes)\n",
" end\n",
" dofmap = create_ldof2gdofmap(elements, \"displacement\")\n",
" Logging.debug(dofmap)\n",
" # Boundary conditions\n",
" # here we want to create nodal force for third dof, that is, node id 2, second dof\n",
" bc1 = BC([dofmap[4][2]], [-0.1])\n",
" bc1 = BC([dofmap[3][2]], [-2.0])\n",
" # dirichlet bc, set dx=dy=0 on support\n",
" bc2 = BC([dofmap[5][1], dofmap[5][2], dofmap[6][1], dofmap[6][2]], [0.0, 0.0, 0.0, 0.0])\n",
" bc2 = BC([dofmap[1][1], dofmap[1][2], dofmap[4][1], dofmap[4][2]], [0.0, 0.0, 0.0, 0.0])\n",
" solve!(elements, dofmap, [bc1], [bc2]; max_iterations=7)\n",
" disp = elements[2].attributes[\"displacement\"]\n",
" disp = elements[1].attributes[\"displacement\"]\n",
" Logging.debug(\"Displacement of element = \\n$disp\")\n",
" #@fact norm(disp) => roughly(3.1292483947150043)\n",
"end"
" @fact norm(disp) --> roughly(3.1292483947150043)\n",
"end\n"
]
},
{
@@ -787,7 +735,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 17,
"metadata": {
"collapsed": false
},
@@ -796,68 +744,44 @@
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
"WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:11.\n",
"Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
"\n",
"WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:32.\n",
"Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
"15-Aug 17:04:20:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n",
"WARNING: beginswith is deprecated, use startswith instead.\n",
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in beginswith at deprecated.jl:30\n",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n",
" in include_string at loading.jl:99\n",
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
" in anonymous at task.jl:365\n",
"while loading In[13], in expression starting on line 3\n",
"WARNING: beginswith is deprecated, use startswith instead.\n",
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in beginswith at deprecated.jl:30\n",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:113\n",
" in include_string at loading.jl:99\n",
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
" in anonymous at task.jl:365\n",
"while loading In[13], in expression starting on line 3\n",
"15-Aug 17:04:21:DEBUG:root:Found NODE section\n",
"15-Aug 17:04:22:DEBUG:root:Found ELEMENT section\n",
"WARNING: integer(s::AbstractString) is deprecated, use parse(Int,s) instead.\n",
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in integer at deprecated.jl:49\n",
" in map at abstractarray.jl:1251\n",
" in parse_element_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:59\n",
" in process_section at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:108\n",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:117\n",
" in include_string at loading.jl:99\n",
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
" in anonymous at task.jl:365\n",
"while loading In[13], in expression starting on line 3\n",
"15-Aug 17:04:27:DEBUG:root:120 elements found\n",
"15-Aug 17:04:27:INFO:root:Creating ELSET Body1\n",
"15-Aug 17:04:27:DEBUG:root:Found NSET section\n",
"15-Aug 17:04:28:DEBUG:root:Creating node set SUPPORT\n",
"15-Aug 17:04:28:DEBUG:root:Found NSET section\n",
"15-Aug 17:04:28:DEBUG:root:Creating node set LOAD\n",
"15-Aug 17:04:28:DEBUG:root:Found NSET section\n",
"15-Aug 17:04:28:DEBUG:root:Creating node set TOP\n"
"19-Aug 18:49:21:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n",
"WARNING: beginswith is deprecated, use startswith instead.\n"
]
},
{
"data": {
"text/plain": [
"Dict{Any,Any} with 4 entries:\n",
" \"nodes\" => Dict{Any,Any}(288=>[97.5,7.5,10.0],11=>[92.5,2.5,5.0],134=>[45.0,10.0,0.0],158=>[2.5,2.5,0.0],160=>[7.5,7.5,0.0],215=>[60.0,0.0,5.0],29=>[2.5,7…\n",
" \"elements\" => Dict{Any,Any}(68=>[71,144,149,198,51,150,57,43,50,214],2=>[204,199,175,130,207,208,209,3,4,176],89=>[95,78,104,52,127,126,106,60,68,67],11=>[15…\n",
" \"elsets\" => Dict{Any,Any}(\"Body1\"=>[1,2,3,4,5,6,7,8,9,10 … 111,112,113,114,115,116,117,118,119,120])\n",
" \"nsets\" => Dict{Any,Any}(\"LOAD\"=>[82,84,87,179,197,246,249,256,257],\"SUPPORT\"=>[108,109,111,155,162,216,225,281,298],\"TOP\"=>[70,75,76,84,88,90,95,96,98,10…"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
"ename": "LoadError",
"evalue": "LoadError: MethodError: `convert` has no method matching convert(::Type{SubString{ASCIIString}}, ::Dict{Any,Any})\nThis may have arisen from a call to the constructor SubString{ASCIIString}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T<:AbstractString}(::Type{T<:AbstractString}, !Matched::AbstractArray{Char,1})\n convert{T<:AbstractString}(::Type{SubString{T<:AbstractString}}, !Matched::T<:AbstractString)\n ...\nwhile loading In[17], in expression starting on line 3",
"output_type": "error",
"traceback": [
"LoadError: MethodError: `convert` has no method matching convert(::Type{SubString{ASCIIString}}, ::Dict{Any,Any})\nThis may have arisen from a call to the constructor SubString{ASCIIString}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n call{T}(::Type{T}, ::Any)\n convert{T<:AbstractString}(::Type{T<:AbstractString}, !Matched::AbstractArray{Char,1})\n convert{T<:AbstractString}(::Type{SubString{T<:AbstractString}}, !Matched::T<:AbstractString)\n ...\nwhile loading In[17], in expression starting on line 3",
"",
" in setindex! at dict.jl:615",
" in parse_header at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:40",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:120"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in beginswith at deprecated.jl:30\n",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:115\n",
" in include_string at loading.jl:99\n",
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
" in anonymous at task.jl:365\n",
"while loading In[17], in expression starting on line 3\n",
"WARNING: beginswith is deprecated, use startswith instead.\n",
" in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
" in beginswith at deprecated.jl:30\n",
" in parse_abaqus at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:115\n",
" in include_string at loading.jl:99\n",
" in execute_request_0x535c5df2 at /Users/jukka/.julia/v0.4/IJulia/src/execute_request.jl:157\n",
" in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
" in anonymous at task.jl:365\n",
"while loading In[17], in expression starting on line 3\n"
]
}
],
"source": [
@@ -1280,6 +1204,53 @@
"JuliaFEM.xdmf.xdmf_save_model(xdoc, \"/tmp/3d_solid_model.xmf\")"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[0.0,0.0,0.0,0.0,0.0])"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = zeros(10)\n",
"u1 = u[1:5]\n",
"u, u1"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[1.0,1.0,1.0,1.0,1.0])"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u1[:] = 1\n",
"u, u1"
]
},
{
"cell_type": "code",
"execution_count": null,
src/abaqus_reader.jl
+4 -2
View File
@@ -8,7 +8,9 @@ using Logging
VERSION < v"0.4-" && using Docile
eldims = Dict({"C3D10" => 10})
eldims = Dict(
"C3D10" => 10,
"C3D4" => 4)
global handlers = Dict()
@@ -29,7 +31,7 @@ end
function parse_header(header_line)
args = map(s -> strip(s), split(header_line, ","))
args[1] = strip(args[1], '*')
d = Dict({"section" => args[1]})
d = Dict("section" => args[1])
options = Dict()
for k in args[2:end]
args2 = split(k, "=")
src/math.jl
+12 -18
View File
@@ -73,9 +73,9 @@ end
"""
Return partial derivatives of shape functions w.r.t X using chain rule.
"""
function get_dbasisdX(el::Element, xi)
J = interpolate(el, "coordinates", xi; derivative=true)
dbasisdX = el.dbasis(xi)*inv(J')
function get_dbasisdX(el::Element, ip)
J = interpolate(el, "coordinates", ip.xi; derivative=true)
dbasisdX = el.dbasis(ip.xi)*inv(J')
return dbasisdX
end
@@ -156,11 +156,9 @@ f::Function
"""
function integrate(f::Function, el::JuliaFEM.Element)
target = []
for m = 1:length(el.iweights)
w = el.iweights[m]
xi = el.ipoints[:, m]
J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true)
push!(target, w*f(el, xi)*det(J))
for ip in el.integration_points
J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
push!(target, ip.weight*f(el, ip)*det(J))
end
return sum(target)
end
@@ -171,11 +169,9 @@ This version returns a function which must be operated with element e
function integrate(f::Function)
function integrate(el::JuliaFEM.Element)
target = []
for m = 1:length(el.iweights)
w = el.iweights[m]
xi = el.ipoints[:, m]
J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true)
push!(target, w*f(el, xi)*det(J))
for ip in el.integration_points
J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
push!(target, ip.weight*f(el, ip)*det(J))
end
return sum(target)
end
@@ -188,10 +184,8 @@ This version saves results inplace to target, garbage collection free
function integrate!(f::Function, el::JuliaFEM.Element, target)
# set target to zero
el.attributes[target][:] = 0.0
for m = 1:length(el.iweights)
w = el.iweights[m]
xi = el.ipoints[:, m]
J = JuliaFEM.interpolate(el, "coordinates", xi; derivative=true)
el.attributes[target][:,:] += w*f(el, xi)*det(J)
for ip in el.integration_points
J = JuliaFEM.interpolate(el, "coordinates", ip.xi; derivative=true)
el.attributes[target][:,:] += ip.weight*f(el, ip)*det(J)
end
end
src/types.jl
+19 -2
View File
@@ -1,6 +1,22 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
Integration point
xi :: Array{Float64, 1}
(dimensionless) coordinates of integration point
weight :: Float64
Integration weight
attributes :: Dict{ASCIIString, Any}
This is used to save internal variables of IP needed e.g. for incremental
material models.
"""
type IntegrationPoint
xi :: Array{Float64, 1}
weight :: Float64
attributes :: Dict{ASCIIString, Any}
end
type Element
id :: Int
@@ -8,9 +24,10 @@ type Element
node_ids :: Array{Int, 1}
basis :: Function
dbasis :: Function
integration_points :: Array{IntegrationPoint, 1}
attributes :: Dict{ASCIIString, Any}
ipoints :: Array{Float64, 2}
iweights :: Array{Float64, 1}
# ipoints :: Array{Float64, 2}
# iweights :: Array{Float64, 1}
end