diff --git a/notebooks/2015-12-13-3d-tie-contact.ipynb b/notebooks/2015-12-13-3d-tie-contact.ipynb new file mode 100644 index 0000000..7b979ee --- /dev/null +++ b/notebooks/2015-12-13-3d-tie-contact.ipynb @@ -0,0 +1,931 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Tie contact 3d\n", + "\n", + "Author(s): Jukka Aho" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solid block\n", + "\n", + "\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM\n", + "using JuliaFEM.Preprocess: parse_aster_med_file\n", + "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Quad4, Hex8, LinearSolver" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: Found 4 element sets: SYM23, SYM12, SYM13, LOAD\n" + ] + }, + { + "data": { + "text/plain": [ + "Dict{ASCIIString,Any} with 2 entries:\n", + " \"nodes\" => Dict(2=>[0.0,0.0,0.0],11=>[0.0,0.3333333333333333,1.0],39=>…\n", + " \"connectivity\" => Dict(68=>(:QU4,:OTHER,[45,47,48,46]),2=>(:SE2,:OTHER,[9,10]…" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/mesh.med\")\n", + "#mesh = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/twoelem_box.med\")" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: created 36 elements.\n" + ] + } + ], + "source": [ + "# interior elements are of type HE8\n", + "field_problem = LinearElasticityProblem()\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " eltype == :HE8 || continue\n", + " element = Hex8(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"youngs modulus\"] = 900.0\n", + " element[\"poissons ratio\"] = 0.25\n", + " push!(field_problem, element)\n", + "end\n", + "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == :LOAD) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"displacement traction force\"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]\n", + " push!(field_problem, element)\n", + "end\n", + "info(\"created $(length(field_problem.elements)) elements.\")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: created 27 boundary elements.\n" + ] + } + ], + "source": [ + "# boundary conditions\n", + "boundary_problem = DirichletProblem(\"displacement\", 3)\n", + "for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " (eltype == :QU4) && (elset in [:SYM23, :SYM12, :SYM13]) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh[\"nodes\"][i] for i in get_connectivity(element)]\n", + " if elset == :SYM23\n", + " element[\"displacement 1\"] = 0.0\n", + " elseif elset == :SYM12\n", + " element[\"displacement 3\"] = 0.0\n", + " elseif elset == :SYM13\n", + " element[\"displacement 2\"] = 0.0\n", + " end\n", + " push!(boundary_problem, element)\n", + "end\n", + "info(\"created $(length(boundary_problem.elements)) boundary elements.\")" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: solving displacement problem, 3 dofs / nodes\n", + "INFO: solved problem in 2.11 seconds.\n" + ] + }, + { + "data": { + "text/plain": [ + "0.587944735792132" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ls = LinearSolver(\"block\")\n", + "push!(ls, field_problem)\n", + "push!(ls, boundary_problem)\n", + "norm = call(ls, 0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: nid near corner = 7\n" + ] + }, + { + "data": { + "text/plain": [ + "7" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "nid = 0\n", + "for (nid, coords) in mesh[\"nodes\"]\n", + " if isapprox(coords, [1.0, 1.0, 1.0])\n", + " info(\"nid near corner = $nid\")\n", + " break\n", + " end\n", + "end\n", + "nid" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: displacement X = [1.0,1.0,1.0], u = [-0.02777777777777779,-0.027777777777777773,0.11111111111111113]\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [-0.02777777777777779,-0.027777777777777773,0.11111111111111113]\n" + ] + } + ], + "source": [ + "using JuliaFEM.Test\n", + "known_value = [1/36, 1/36, -1/9]\n", + "for element in field_problem.elements\n", + " i = indexin([nid], get_connectivity(element))[1]\n", + " i != 0 || continue\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " info(\"displacement X = $(X[i]), u = $(u[i])\")\n", + " @test isapprox(-u[i], known_value)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: XDFM: ndim = 192\n" + ] + } + ], + "source": [ + "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n", + "coll = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n", + "grid = JuliaFEM.Postprocess.xdmf_new_grid(coll; time=0.0)\n", + "\n", + "Xg = Dict{Int64, Vector{Float64}}()\n", + "ug = Dict{Int64, Vector{Float64}}()\n", + "for element in field_problem.elements\n", + " conn = get_connectivity(element)\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " for (i, c) in enumerate(conn)\n", + " Xg[c] = X[i]\n", + " ug[c] = u[i]\n", + " end\n", + "end\n", + "perm = sort(collect(keys(Xg)))\n", + "nodes = Vector{Float64}[Xg[i] for i in perm]\n", + "disp = Vector{Float64}[ug[i] for i in perm]\n", + "elements = []\n", + "for el in field_problem.elements\n", + " isa(el, JuliaFEM.Core.Element{JuliaFEM.Core.Hex8}) || continue\n", + " push!(elements, (:Hex8, get_connectivity(el)))\n", + "end\n", + "#elements\n", + "JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements)\n", + "JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, \"displacement\", disp)\n", + "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/foobar2.xmf\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3x56 Array{Float64,2}:\n", + " 0.0 0.0 0.0 0.0 0.0 0.0 … 0.0 0.0 0.0 \n", + " 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 \n", + " -2.77778 0.0 -2.77778 0.0 -2.77778 0.0 -11.1111 -11.1111 -11.1111" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "reshape(full(JuliaFEM.Core.assemble(field_problem, 0.0).force_vector), 3, 56)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Block divived to two parts\n", + "\n", + "\n", + "" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: Found 4 element sets: LSYM23, LSYM12, LSYM13, LOWER_BLOCK_TO_UPPER_BLOCK\n", + "INFO: Found 4 element sets: USYM23, UPPER_BLOCK_TO_LOWER_BLOCK, USYM13, LOAD\n" + ] + }, + { + "data": { + "text/plain": [ + "Dict{ASCIIString,Any} with 2 entries:\n", + " \"nodes\" => Dict(2=>[0.0,0.0,1.0],11=>[0.0,0.6666666666666666,1.0],39=>…\n", + " \"connectivity\" => Dict(68=>(:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[42,14,4,26]),2=…" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using JuliaFEM\n", + "using JuliaFEM.Preprocess: parse_aster_med_file\n", + "using JuliaFEM.Core: LinearElasticityProblem, DirichletProblem, get_connectivity, Quad4, Hex8, LinearSolver\n", + "mesh_lower = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/LOWER_BLOCK.med\")\n", + "mesh_upper = parse_aster_med_file(Pkg.dir(\"JuliaFEM\")*\"/geometry/unit_box/UPPER_BLOCK.med\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x75 Array{Int64,2}:\n", + " 1 2 3 4 5 6 7 8 9 10 11 12 … 67 68 69 70 71 72 73 74 75" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sort(collect(keys(mesh_lower[\"nodes\"])))'" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x48 Array{Int64,2}:\n", + " 1 2 3 4 5 6 7 8 9 10 11 12 … 40 41 42 43 44 45 46 47 48" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sort(collect(keys(mesh_upper[\"nodes\"])))'" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Need to renumber nodes." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Dict{Int64,Int64} with 48 entries:\n", + " 2 => 77\n", + " 11 => 86\n", + " 39 => 114\n", + " 46 => 121\n", + " 25 => 100\n", + " 42 => 117\n", + " 29 => 104\n", + " 8 => 83\n", + " 20 => 95\n", + " 14 => 89\n", + " 31 => 106\n", + " 33 => 108\n", + " 18 => 93\n", + " 26 => 101\n", + " 35 => 110\n", + " 17 => 92\n", + " 44 => 119\n", + " 4 => 79\n", + " 37 => 112\n", + " 45 => 120\n", + " 13 => 88\n", + " 30 => 105\n", + " 1 => 76\n", + " 47 => 122\n", + " 32 => 107\n", + " ⋮ => ⋮" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "new_node_numbering = Dict{Int64, Int64}()\n", + "upper_node_ids = sort(collect(keys(mesh_upper[\"nodes\"])))\n", + "for (k, old_node_id) in enumerate(upper_node_ids)\n", + " new_node_numbering[old_node_id] = 75+k\n", + "end\n", + "new_node_numbering" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1x48 Array{Int64,2}:\n", + " 76 77 78 79 80 81 82 83 84 85 … 117 118 119 120 121 122 123" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function renumber!(mesh, node_numbering)\n", + " old_nodes = mesh[\"nodes\"]\n", + " new_nodes = typeof(old_nodes)()\n", + " for (node_id, node_coords) in old_nodes\n", + " new_node_id = node_numbering[node_id]\n", + " new_nodes[new_node_id] = node_coords\n", + " end\n", + " mesh[\"nodes\"] = new_nodes\n", + " for (elid, (eltype, elset, elcon)) in mesh[\"connectivity\"]\n", + " new_elcon = [node_numbering[node_id] for node_id in elcon]\n", + " mesh[\"connectivity\"][elid] = (eltype, elset, new_elcon)\n", + " end\n", + " return mesh\n", + "end\n", + "\n", + "renumber!(mesh_upper, new_node_numbering)\n", + "sort(collect(keys(mesh_upper[\"nodes\"])))'" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Dict{Int64,Tuple{Symbol,Symbol,Array{Int64,1}}} with 92 entries:\n", + " 68 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[117,89,79,101])\n", + " 2 => (:SE2,:OTHER,[84,77])\n", + " 89 => (:HE8,:OTHER,[121,123,115,114,118,119,103,102])\n", + " 11 => (:SE2,:OTHER,[80,90])\n", + " 39 => (:QU4,:USYM13,[80,90,106,93])\n", + " 46 => (:QU4,:LOAD,[95,108,109,96])\n", + " 85 => (:HE8,:OTHER,[107,120,116,94,106,121,114,93])\n", + " 25 => (:SE2,:OTHER,[83,99])\n", + " 55 => (:QU4,:OTHER,[96,112,113,95])\n", + " 42 => (:QU4,:USYM13,[90,81,92,106])\n", + " 66 => (:QU4,:UPPER_BLOCK_TO_LOWER_BLOCK,[94,76,88,116])\n", + " 58 => (:QU4,:OTHER,[112,100,101,113])\n", + " 29 => (:SE2,:OTHER,[101,79])\n", + " 59 => (:QU4,:OTHER,[113,101,79,87])\n", + " 8 => (:SE2,:OTHER,[76,88])\n", + " 74 => (:QU4,:OTHER,[119,103,83,99])\n", + " 90 => (:HE8,:OTHER,[105,87,79,89,122,113,101,117])\n", + " 57 => (:QU4,:OTHER,[99,83,100,112])\n", + " 20 => (:SE2,:OTHER,[95,96])\n", + " 78 => (:HE8,:OTHER,[85,86,105,104,110,108,122,120])\n", + " 14 => (:SE2,:OTHER,[91,92])\n", + " 31 => (:SE2,:OTHER,[102,103])\n", + " 70 => (:QU4,:OTHER,[97,118,119,98])\n", + " 33 => (:QU4,:USYM23,[76,84,104,88])\n", + " 52 => (:QU4,:LOAD,[110,91,92,111])\n", + " ⋮ => ⋮" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mesh_upper[\"connectivity\"]" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: lower: add element with connectivity [40,69,57,23,17,64,34,5]\n", + "INFO: lower: add element with connectivity [73,54,31,61,74,53,32,62]\n", + "INFO: lower: add element with connectivity [75,52,33,63,66,30,8,36]\n", + "INFO: lower: add element with connectivity [71,74,62,59,72,75,63,60]\n", + "INFO: lower: add element with connectivity [46,47,71,68,49,50,72,69]\n", + "INFO: lower: add element with connectivity [49,50,72,69,24,25,65,64]\n", + "INFO: lower: add element with connectivity [20,43,67,42,19,46,68,41]\n", + "INFO: lower: add element with connectivity [9,37,14,1,42,67,55,21]\n", + "INFO: lower: add element with connectivity [44,45,73,70,47,48,74,71]\n", + "INFO: lower: add element with connectivity [10,11,38,37,43,44,70,67]\n", + "INFO: lower: add element with connectivity [39,13,4,16,73,54,31,61]\n", + "INFO: lower: add element with connectivity [74,53,32,62,75,52,33,63]\n", + "INFO: lower: add element with connectivity [37,38,15,14,67,70,58,55]\n", + "INFO: lower: add element with connectivity [70,73,61,58,71,74,62,59]\n", + "INFO: lower: add element with connectivity [43,44,70,67,46,47,71,68]\n", + "INFO: lower: add element with connectivity [68,71,59,56,69,72,60,57]\n", + "INFO: lower: add element with connectivity [12,3,13,39,45,29,54,73]\n", + "INFO: lower: add element with connectivity [67,70,58,55,68,71,59,56]\n", + "INFO: lower: add element with connectivity [47,48,74,71,50,51,75,72]\n", + "INFO: lower: add element with connectivity [2,10,37,9,20,43,67,42]\n", + "INFO: lower: add element with connectivity [69,72,60,57,64,65,35,34]\n", + "INFO: lower: add element with connectivity [72,75,63,60,65,66,36,35]\n", + "INFO: lower: add element with connectivity [45,29,54,73,48,28,53,74]\n", + "INFO: lower: add element with connectivity [51,27,52,75,26,7,30,66]\n", + "INFO: lower: add element with connectivity [11,12,39,38,44,45,73,70]\n", + "INFO: lower: add element with connectivity [42,67,55,21,41,68,56,22]\n", + "INFO: lower: add element with connectivity [19,46,68,41,18,49,69,40]\n", + "INFO: lower: add element with connectivity [38,39,16,15,70,73,61,58]\n", + "INFO: lower: add element with connectivity [48,28,53,74,51,27,52,75]\n", + "INFO: lower: add element with connectivity [50,51,75,72,25,26,66,65]\n", + "INFO: lower: add element with connectivity [41,68,56,22,40,69,57,23]\n", + "INFO: lower: add element with connectivity [18,49,69,40,6,24,64,17]\n", + "INFO: upper: add element with connectivity [121,123,115,114,118,119,103,102]\n", + "INFO: upper: add element with connectivity [107,120,116,94,106,121,114,93]\n", + "INFO: upper: add element with connectivity [105,87,79,89,122,113,101,117]\n", + "INFO: upper: add element with connectivity [85,86,105,104,110,108,122,120]\n", + "INFO: upper: add element with connectivity [109,96,112,123,98,82,99,119]\n", + "INFO: upper: add element with connectivity [84,104,88,76,107,120,116,94]\n", + "INFO: upper: add element with connectivity [106,121,114,93,90,118,102,80]\n", + "INFO: upper: add element with connectivity [92,111,121,106,81,97,118,90]\n", + "INFO: upper: add element with connectivity [111,109,123,121,97,98,119,118]\n", + "INFO: upper: add element with connectivity [108,95,113,122,109,96,112,123]\n", + "INFO: upper: add element with connectivity [122,113,101,117,123,112,100,115]\n", + "INFO: upper: add element with connectivity [77,85,104,84,91,110,120,107]\n", + "INFO: upper: add element with connectivity [104,105,89,88,120,122,117,116]\n", + "INFO: upper: add element with connectivity [110,108,122,120,111,109,123,121]\n", + "INFO: upper: add element with connectivity [91,110,120,107,92,111,121,106]\n", + "INFO: upper: add element with connectivity [86,78,87,105,108,95,113,122]\n", + "INFO: upper: add element with connectivity [120,122,117,116,121,123,115,114]\n", + "INFO: upper: add element with connectivity [123,112,100,115,119,99,83,103]\n", + "INFO: created 59 elements.\n" + ] + } + ], + "source": [ + "# interior elements are of type HE8\n", + "field_problem = LinearElasticityProblem()\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n", + " eltype == :HE8 || continue\n", + " element = Hex8(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"youngs modulus\"] = 900.0\n", + " element[\"poissons ratio\"] = 0.25\n", + " info(\"lower: add element with connectivity $elcon\")\n", + " push!(field_problem, element)\n", + "end\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", + " eltype == :HE8 || continue\n", + " element = Hex8(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"youngs modulus\"] = 900.0\n", + " element[\"poissons ratio\"] = 0.25\n", + " info(\"upper: add element with connectivity $elcon\")\n", + " push!(field_problem, element)\n", + "end\n", + "\n", + "# Neumann boundary condition, traction force -100 on Z direction for element set LOAD\n", + "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == :LOAD) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"displacement traction force\"] = Vector{Float64}[[0.0, 0.0, 100.0] for i=1:4]\n", + " push!(field_problem, element)\n", + "end\n", + "info(\"created $(length(field_problem.elements)) elements.\")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: created 44 boundary elements.\n" + ] + } + ], + "source": [ + "# boundary conditions\n", + "boundary_problem = DirichletProblem(\"displacement\", 3)\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh_lower[\"connectivity\"]\n", + " (eltype == :QU4) && (elset in [:LSYM23, :LSYM12, :LSYM13]) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh_lower[\"nodes\"][i] for i in get_connectivity(element)]\n", + " if elset == :LSYM23\n", + " element[\"displacement 1\"] = 0.0\n", + " elseif elset == :LSYM12\n", + " element[\"displacement 3\"] = 0.0\n", + " elseif elset == :LSYM13\n", + " element[\"displacement 2\"] = 0.0\n", + " end\n", + " push!(boundary_problem, element)\n", + "end\n", + "\n", + "for (elid, (eltype, elset, elcon)) in mesh_upper[\"connectivity\"]\n", + " (eltype == :QU4) && (elset in [:USYM23, :USYM12, :USYM13]) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mesh_upper[\"nodes\"][i] for i in get_connectivity(element)]\n", + " if elset == :USYM23\n", + " element[\"displacement 1\"] = 0.0\n", + " elseif elset == :USYM12\n", + " element[\"displacement 3\"] = 0.0\n", + " elseif elset == :USYM13\n", + " element[\"displacement 2\"] = 0.0\n", + " end\n", + " push!(boundary_problem, element)\n", + "end\n", + "\n", + "info(\"created $(length(boundary_problem.elements)) boundary elements.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Creating tie contact\n", + "- define slave element surface (the one where integration happend)\n", + "- define potential master elements for slave elements" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!\n", + "\n", + "mortar_side = mesh_upper\n", + "nonmortar_side = mesh_lower\n", + "\n", + "master_elements = JuliaFEM.Core.Element[]\n", + "for (elid, (eltype, elset, elcon)) in mortar_side[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == :UPPER_BLOCK_TO_LOWER_BLOCK) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[mortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n", + " push!(master_elements, element)\n", + "end\n", + "\n", + "contact_problem = MortarProblem(\"displacement\", 3)\n", + "for (elid, (eltype, elset, elcon)) in nonmortar_side[\"connectivity\"]\n", + " (eltype == :QU4) && (elset == :LOWER_BLOCK_TO_UPPER_BLOCK) || continue\n", + " element = Quad4(elcon)\n", + " element[\"geometry\"] = Vector{Float64}[nonmortar_side[\"nodes\"][i] for i in get_connectivity(element)]\n", + " element[\"master elements\"] = master_elements\n", + " calculate_normal_tangential_coordinates!(element, 0.0)\n", + " push!(contact_problem, element)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: # of field problems: 1\n", + "INFO: # of boundary problems: 2\n", + "INFO: Starting iteration 1\n", + "INFO: Assembling field problems...\n", + "INFO: Assembling body 1...\n", + "INFO: dim = 369\n", + "INFO: Assembling boundary problems...\n", + "INFO: Assembling boundary 1...\n", + "INFO: Assembling boundary 2...\n", + "INFO: dumping matrices to disk, file = matrices_tie_contact_3d_host_1_iteration_1.jld\n", + "INFO: Solving system\n", + "INFO: UMFPACK: solved in 0.1845378875732422 seconds. norm = 0.37585220926688956\n", + "INFO: timing info for non-linear iteration:\n" + ] + }, + { + "data": { + "text/plain": [ + "(1,false)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: boundary assembly : 1.6115069389343262\n", + "INFO: field assembly : 1.110853910446167\n", + "INFO: dump matrices to disk : 1.0786750316619873\n", + "INFO: solve problem : 0.3312990665435791\n", + "INFO: update element data : 0.022845029830932617\n", + "INFO: non-linear iteration : 4.155207872390747\n", + "INFO: Warning: did not coverge in 1 iterations!\n" + ] + } + ], + "source": [ + "using JuliaFEM.Core: DirectSolver\n", + "solver = DirectSolver()\n", + "solver.name = \"tie_contact_3d\"\n", + "solver.method = :UMFPACK\n", + "solver.max_iterations = 1\n", + "solver.dump_matrices = true\n", + "push!(solver, field_problem)\n", + "push!(solver, boundary_problem)\n", + "push!(solver, contact_problem)\n", + "call(solver, 0.0)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: nid near corner = 82\n", + "INFO: displacement X = [1.0,1.0,1.0], u = [0.022782767229558416,0.022002026379764242,-0.08215452308106556]\n" + ] + }, + { + "ename": "LoadError", + "evalue": "LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14", + "output_type": "error", + "traceback": [ + "LoadError: test failed: isapprox(-(u[i]),known_value)\n in expression: isapprox(-(u[i]),known_value)\nwhile loading In[11], in expression starting on line 14", + "", + " in error at ./error.jl:21", + " in default_handler at test.jl:30", + " in do_test at test.jl:53", + " [inlined code] from In[11]:20", + " in anonymous at no file:0" + ] + } + ], + "source": [ + "using JuliaFEM.Test\n", + "\n", + "nid = 0\n", + "for (nid, coords) in mesh_upper[\"nodes\"]\n", + " if isapprox(coords, [1.0, 1.0, 1.0])\n", + " info(\"nid near corner = $nid\")\n", + " break\n", + " end\n", + "end\n", + "nid\n", + "\n", + "known_value = [1/36, 1/36, -1/9]\n", + "\n", + "for element in field_problem.elements\n", + " i = indexin([nid], get_connectivity(element))[1]\n", + " i != 0 || continue\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " info(\"displacement X = $(X[i]), u = $(u[i])\")\n", + " @test isapprox(-u[i], known_value)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "INFO: XDFM: ndim = 369\n" + ] + } + ], + "source": [ + "xdoc, xmodel = JuliaFEM.Postprocess.xdmf_new_model()\n", + "coll = JuliaFEM.Postprocess.xdmf_new_temporal_collection(xmodel)\n", + "grid = JuliaFEM.Postprocess.xdmf_new_grid(coll; time=0.0)\n", + "\n", + "Xg = Dict{Int64, Vector{Float64}}()\n", + "ug = Dict{Int64, Vector{Float64}}()\n", + "for element in field_problem.elements\n", + " conn = get_connectivity(element)\n", + " X = element(\"geometry\", 0.0)\n", + " u = element(\"displacement\", 0.0)\n", + " for (i, c) in enumerate(conn)\n", + " Xg[c] = X[i]\n", + " ug[c] = u[i]\n", + " end\n", + "end\n", + "perm = sort(collect(keys(Xg)))\n", + "nodes = Vector{Float64}[Xg[i] for i in perm]\n", + "disp = Vector{Float64}[ug[i] for i in perm]\n", + "elements = []\n", + "for el in field_problem.elements\n", + " isa(el, JuliaFEM.Core.Element{JuliaFEM.Core.Hex8}) || continue\n", + " push!(elements, (:Hex8, get_connectivity(el)))\n", + "end\n", + "#elements\n", + "JuliaFEM.Postprocess.xdmf_new_mesh!(grid, nodes, elements)\n", + "JuliaFEM.Postprocess.xdmf_new_nodal_field!(grid, \"displacement\", disp)\n", + "JuliaFEM.Postprocess.xdmf_save_model(xdoc, \"/tmp/foobar3.xmf\");" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.2", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "0.4.2" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index a55fc49..1ad400a 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -20,7 +20,7 @@ end module Preprocess include("abaqus_reader.jl") -include("aster_reader.jl") +include("preprocess_aster_reader.jl") end module Postprocess diff --git a/src/aster_reader.jl b/src/aster_reader.jl deleted file mode 100644 index 7121ab7..0000000 --- a/src/aster_reader.jl +++ /dev/null @@ -1,53 +0,0 @@ -# This file is a part of JuliaFEM. -# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md - - -function aster_parse_nodes(section::ASCIIString; strip_characters=true) - nodes = Dict{Any, Vector{Float64}}() - has_started = false - for line in split(section, '\n') - m = matchall(r"[\w.-]+", line) - if (length(m) != 1) && (!has_started) - continue - end - if length(m) == 1 - if (m[1] == "COOR_2D") || (m[1] == "COOR_3D") - has_started = true - continue - end - if m[1] == "FINSF" - break - end - end - if length(m) == 4 - nid = m[1] - if strip_characters - nid = matchall(r"\d", nid) - nid = parse(Int, nid[1]) - end - nodes[nid] = float(m[2:end]) - end - end - return nodes -end - -function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}}) - model = Dict{ASCIIString, Any}() - header = nothing - data = ASCIIString[] - for line in split(mesh, '\n') - length(line) != 0 || continue - info("line: $line") - if is_aster_mail_keyword(strip(line)) - header = parse_aster_header(line) - empty!(data) - continue - end - if line == "FINSF" - info(data) - header = nothing - process_aster_section!(model, join(data, ""), header, Val{header[1]}) - end - end - return model -end diff --git a/src/directsolver.jl b/src/directsolver.jl index 3912118..6f463a7 100644 --- a/src/directsolver.jl +++ b/src/directsolver.jl @@ -207,6 +207,7 @@ function call(solver::DirectSolver, time::Number=0.0) end K = sparse(field_assembly.stiffness_matrix) dim = size(K, 1) + info("dim = $dim") f = sparse(field_assembly.force_vector, dim, 1) field_assembly = nothing gc() diff --git a/src/dirichlet.jl b/src/dirichlet.jl index b7c473b..82b560c 100644 --- a/src/dirichlet.jl +++ b/src/dirichlet.jl @@ -15,7 +15,14 @@ function assemble!(assembly::Assembly, problem::BoundaryProblem{DirichletProblem gdofs = get_gdofs(element, field_dim) for ip in get_integration_points(element, Val{2}) - w = ip.weight * det(element, ip, time) + w = ip.weight + J = get_jacobian(element, ip, time) + JT = transpose(J) + if size(JT, 2) == 1 # plane problem + w *= norm(JT) + else + w *= norm(cross(JT[:,1], JT[:,2])) + end N = element(ip, time) A = w*N'*N diff --git a/src/elements.jl b/src/elements.jl index 1510562..58639f7 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -8,6 +8,10 @@ type Element{E<:AbstractElement} fields :: Dict{ASCIIString, Field} end +function Base.size{E}(::Element{E}) + return size(E) +end + function convert{E}(::Type{Element{E}}, connectivity::Vector{Int}) # return Element{E}(connectivity, get_integration_points(E), Dict()) return Element{E}(connectivity, Dict()) @@ -211,15 +215,20 @@ end """ Return the determinant of jacobian. """ function LinAlg.det{E<:AbstractElement}(element::Element{E}, xi::Vector{Float64}, time::Real) - warn("det(element, ip, time) is ambiguous: use J = get_jacobian(element, ip, time); det(J) instead.") J = get_jacobian(element, xi, time) n, m = size(J) if n == m + warn("det(element, ip, time) is ambiguous: use J = get_jacobian(element, ip, time); det(J) instead.") return det(J) end JT = transpose(J) - s = size(JT, 2) == 1 ? norm(JT) : norm(cross(JT[:,1], JT[:,2])) - return s + if size(JT, 2) == 1 + warn("det(element, ip, time) is ambiguous: use J = get_jacobian(element, ip, time); norm(J) instead.") + return norm(JT) + else + warn("det(element, ip, time) is ambiguous: use J = get_jacobian(element, ip, time); norm(cross(...)) instead.") + return norm(cross(JT[:,1], JT[:,2])) + end end function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Real) return det(element, ip.xi, time) diff --git a/src/mortar.jl b/src/mortar.jl index 934cb60..1206bad 100644 --- a/src/mortar.jl +++ b/src/mortar.jl @@ -376,6 +376,29 @@ function get_points_inside_triangle(Y::Matrix, X::Matrix) return P end +""" +Determine is point P inside or on boudary of polygon X. + +http://paulbourke.net/geometry/polygonmesh/#insidepoly +""" +function is_point_inside_convex_polygon(P, X) + x, y = P + for i=1:length(X) + x0, y0 = X[i] + x1, y1 = X[mod(i, length(X))+1] + if (y-y0)*(x1-x0) - (x-x0)*(y1-y0) < 0 + return false + end + end + return true +end + +function get_points_inside_convex_polygon(pts, X) + # TODO: Make more readable + X2 = [X[:,i] for i=1:size(X,2)] + c = filter(P->is_point_inside_convex_polygon(P, X2), [pts[:,i] for i=1:size(pts, 2)]) + return length(c) == 0 ? zeros(2, 0) : hcat(c...) +end """ Return unique objects with some given tolerance. This is used in next function because traditional unique() command returns row vectors as non-unique if they @@ -434,9 +457,18 @@ julia> n """ function clip_polygon(S::Matrix, M::Matrix) P1, neighbours = get_edge_intersections(M, S) - P2 = get_points_inside_triangle(M, S) - P3 = get_points_inside_triangle(S, M) + #P2 = get_points_inside_triangle(M, S) + #P3 = get_points_inside_triangle(S, M) + P2 = get_points_inside_convex_polygon(M, S) + P3 = get_points_inside_convex_polygon(S, M) +# info("polygon clipping: P1 = $P1") +# info("polygon clipping: P2 = $P2") +# info("polygon clipping: P3 = $P3") +# info("hcat P = $P") P = hcat(P1, P2, P3) + if length(P) == 0 + return nothing, nothing + end P = uniquetol(P, 2) meanval = mean(P, 2) tmp = P .- meanval @@ -654,7 +686,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::BoundaryPro end -typealias MortarElements3D Union{Tri3} +typealias MortarElements3D Union{Tri3, Quad4} function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryProblem{MortarProblem}, slave_element::Element{E}, time::Real) field_dim = problem.parent_field_dim @@ -666,13 +698,14 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro # create auxiliary plane and project slave nodes to it # x0 = origo, Q = local basis x0, Q = create_auxiliary_plane(slave_element, time) - S = Vector{Float64}[] + Sl = Vector{Float64}[] for p in slave_element("geometry", time) - push!(S, project_point_to_auxiliary_plane(p, x0, Q)) + push!(Sl, project_point_to_auxiliary_plane(p, x0, Q)) end - S = reshape([S...;], 2, 3) + #S = reshape([S...;], 2, size(slave_element)[2]) + S = hcat(Sl...) - integration_points = get_integration_points(E, Val{5}) + integration_points = get_integration_points(Tri3, Val{5}) for master_element in slave_element["master elements"] master_dofs = get_gdofs(master_element, field_dim) @@ -681,25 +714,48 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro for p in master_element("geometry", time) push!(M, project_point_to_auxiliary_plane(p, x0, Q)) end - M = reshape([M...;], 2, 3) - P, neighbours = clip_polygon(S, M) + #M = reshape([M...;], 2, size(master_element)[2]) + M = hcat(M...) + P = nothing + neighbours = nothing + try + P, neighbours = clip_polygon(S, M) + catch + info("polygon clipping failed") + info("S = ") + dump(S) + info("M = ") + dump(M) + info("original Sl = ") + info(Sl) + error("cannot continue") + end + isa(P, Void) && continue # no clipping +# info("polygon on auxilyary plane: ") +# dump(round(P, 3)) C = calculate_polygon_centerpoint(P) +# info("center point = $C") npts = size(P, 2) # number of vertices in polygon +# info("number of vectices in polygon: $npts") # S = zeros(3, 3) # M = zeros(3, 3) for i=1:npts # loop vertices and create temporary integrate cells xvec = [C[1], P[1, i], P[1, mod(i, npts)+1]] yvec = [C[2], P[2, i], P[2, mod(i, npts)+1]] X = hcat(xvec, yvec)' +# info("cell $i, coords = ") +# dump(round(X, 3)) geom = Field(Vector{Float64}[X[:,j] for j=1:size(X,2)]) for ip in integration_points # calculate determiant of jacobian - dN = get_dbasis(E, ip.xi) + #dN = get_dbasis(E, ip.xi) + dN = get_dbasis(Tri3, ip.xi) J = sum([kron(dN[:,j], geom[j]') for j=1:length(geom)]) w = ip.weight*det(J) # gauss point in auxiliary plane - N = get_basis(E, ip.xi) + #N = get_basis(E, ip.xi) + N = get_basis(Tri3, ip.xi) x = vec(N*geom) # find projection of gauss point to master and slave elements theta1 = project_point_from_plane_to_surface(x, x0, Q, slave_element, time) @@ -707,13 +763,13 @@ function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::BoundaryPro # evaluate shape functions values in gauss point and add contribution to matrices N1 = slave_element(theta1[2:3], time) N2 = master_element(theta2[2:3], time) - S = w*N1'*N1 - M = w*N1'*N2 + Sm = w*N1'*N1 + Mm = w*N1'*N2 for k=1:field_dim sd = slave_dofs[k:field_dim:end] md = master_dofs[k:field_dim:end] - add!(assembly.stiffness_matrix, sd, sd, S) - add!(assembly.stiffness_matrix, sd, md, -M) + add!(assembly.stiffness_matrix, sd, sd, Sm) + add!(assembly.stiffness_matrix, sd, md, -Mm) # info("sd = $sd") # info("md = $md") end diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl new file mode 100644 index 0000000..c833dad --- /dev/null +++ b/src/preprocess_aster_reader.jl @@ -0,0 +1,157 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using HDF5 + +function aster_parse_nodes(section::ASCIIString; strip_characters=true) + nodes = Dict{Any, Vector{Float64}}() + has_started = false + for line in split(section, '\n') + m = matchall(r"[\w.-]+", line) + if (length(m) != 1) && (!has_started) + continue + end + if length(m) == 1 + if (m[1] == "COOR_2D") || (m[1] == "COOR_3D") + has_started = true + continue + end + if m[1] == "FINSF" + break + end + end + if length(m) == 4 + nid = m[1] + if strip_characters + nid = matchall(r"\d", nid) + nid = parse(Int, nid[1]) + end + nodes[nid] = float(m[2:end]) + end + end + return nodes +end + +function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}}) + model = Dict{ASCIIString, Any}() + header = nothing + data = ASCIIString[] + for line in split(mesh, '\n') + length(line) != 0 || continue + info("line: $line") + if is_aster_mail_keyword(strip(line)) + header = parse_aster_header(line) + empty!(data) + continue + end + if line == "FINSF" + info(data) + header = nothing + process_aster_section!(model, join(data, ""), header, Val{header[1]}) + end + end + return model +end + +""" +Code Aster binary file (.med), which is exported from SALOME. +""" +type MEDFile + data :: Dict +end + +function MEDFile(fn::ASCIIString) + MEDFile(h5read(fn, "/")) +end + +function get_mesh_names(med::MEDFile) + return collect(keys(med.data["FAS"])) +end + +function get_nodes(med::MEDFile, mesh_name) + increments = keys(med.data["ENS_MAA"][mesh_name]) + @assert length(increments) == 1 + increment = first(increments) + nodes = med.data["ENS_MAA"][mesh_name][increment]["NOE"] + node_ids = nodes["NUM"] + nnodes = length(node_ids) + node_coords = nodes["COO"] + dim = round(Int, length(node_coords)/nnodes) + node_coords = reshape(node_coords, nnodes, dim)' + d = Dict{Int64}{Vector{Float64}}() + for i=1:nnodes + d[node_ids[i]] = node_coords[:, i] + end + return d +end + +function get_element_sets(med::MEDFile, mesh_name) + es = Dict{Int64, Symbol}() + if !haskey(med.data["FAS"][mesh_name], "ELEME") + return es + end + elsets = med.data["FAS"][mesh_name]["ELEME"] + for elset in keys(elsets) + k = split(elset, '_') + elset_id = parse(Int, k[2]) + elset_name = ascii(pointer(convert(Vector{UInt8}, elsets[elset]["GRO"]["NOM"][1]))) + es[elset_id] = Symbol(elset_name) + end + return es +end + +function get_connectivity(med::MEDFile, elsets, mesh_name) + elsets[0] = :OTHER + increments = keys(med.data["ENS_MAA"][mesh_name]) + @assert length(increments) == 1 + increment = first(increments) + all_elements = med.data["ENS_MAA"][mesh_name][increment]["MAI"] + d = Dict{Int64, Tuple{Symbol, Symbol, Vector{Int64}}}() + for eltype in keys(all_elements) + elements = all_elements[eltype] + elset_ids = elements["FAM"] + element_ids = elements["NUM"] + nelements = length(element_ids) + element_connectivity = elements["NOD"] + element_dim = round(Int, length(element_connectivity)/nelements) + element_connectivity = reshape(element_connectivity, nelements, element_dim)' + for i=1:nelements + d[element_ids[i]] = (Symbol(eltype), Symbol(elsets[elset_ids[i]]), element_connectivity[:, i]) + end + end + return d +end + +""" Parse code aster .med file. + +Paramters +--------- +fn :: ASCIIString + file name to parse +mesh_name :: ASCIIString, optional + mesh name, if several meshes in one file + +Returns +------- +Dict containing fields "nodes" and "connectivity". + +""" +function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing) + med = MEDFile(fn) + if isa(mesh_name, Void) + mesh_names = get_mesh_names(med::MEDFile) + all_meshes = join(mesh_names, ", ") + length(mesh_names) == 1 || error("several meshes found from med, pick one: $all_meshes") + mesh_name = mesh_names[1] + end + elsets = get_element_sets(med, mesh_name) + elset_names = join(values(elsets), ", ") + info("Found $(length(elsets)) element sets: $elset_names") + nodes = get_nodes(med, mesh_name) + conn = get_connectivity(med, elsets, mesh_name) + result = Dict{ASCIIString, Any}() + result["nodes"] = nodes + result["connectivity"] = conn + return result +end + diff --git a/src/solvers.jl b/src/solvers.jl index dc7ada6..da3491a 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -87,6 +87,7 @@ common situation, i.e., some main field problem and it's Dirichlet boundary. """ function call(solver::LinearSolver, time::Float64) + t0 = Base.time() field_name = get_unknown_field_name(solver.field_problems[1]) field_dim = get_unknown_field_dimension(solver.field_problems[1]) info("solving $field_name problem, $field_dim dofs / nodes") @@ -94,7 +95,7 @@ function call(solver::LinearSolver, time::Float64) field_assembly = assemble(solver.field_problems[1], time) boundary_assembly = assemble(solver.boundary_problems[1], time) - info("Creating sparse matrices") + #info("Creating sparse matrices") K = sparse(field_assembly.stiffness_matrix) dim = size(K, 1) f = sparse(field_assembly.force_vector, dim, 1) @@ -132,6 +133,8 @@ function call(solver::LinearSolver, time::Float64) end end + t1 = round(Base.time()-t0, 2) + info("solved problem in $t1 seconds.") return norm(u) end diff --git a/src/xdmf.jl b/src/xdmf.jl index 4f87610..170aa5e 100644 --- a/src/xdmf.jl +++ b/src/xdmf.jl @@ -37,8 +37,9 @@ using LightXML # > #define XDMF_3DCORECTMESH 0x1102 global eltypes = Dict{Symbol, Int}( - :Tet4 => 0x6, :Quad4 => 0x5, + :Tet4 => 0x6, + :Hex8 => 0x9, :Tet10 => 0x0026) function xdmf_new_model(xdmf_version="2.1") diff --git a/test/test_linear_elasticity.jl b/test/test_linear_elasticity.jl index aaa5e2c..3eb13f6 100644 --- a/test/test_linear_elasticity.jl +++ b/test/test_linear_elasticity.jl @@ -7,7 +7,8 @@ using JuliaFEM using JuliaFEM.Test using JuliaFEM.Core: Seg2, Quad4, Hex8, LinearElasticityProblem, get_connectivity, - assemble, PlaneStressLinearElasticityProblem + assemble, PlaneStressLinearElasticityProblem, DirichletProblem, + LinearSolver using JuliaFEM.Preprocess: aster_parse_nodes @@ -85,7 +86,7 @@ function test_continuum_elasticity_with_surface_load() problem = LinearElasticityProblem() push!(problem, element1) push!(problem, element2) - +#= free_dofs = zeros(Bool, 8, 3) x = 1 y = 2 @@ -107,15 +108,38 @@ function test_continuum_elasticity_with_surface_load() # dump(reshape(f, 3, 8)) # info("initial stiffness matrix") # dump(round(Int, K)[free_dofs, free_dofs]) - u = zeros(3, 8) u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs] - info("result vector") dump(u) +=# + + dx = Quad4([1, 4, 8, 5]) + dx["displacement 1"] = 0.0 + dy = Quad4([1, 5, 6, 2]) + dy["displacement 2"] = 0.0 + dz = Quad4([1, 2, 3, 4]) + dz["displacement 3"] = 0.0 + bc = DirichletProblem("displacement", 3) + for el in [dx, dy, dz] + set_geometry!(el, nodes) + push!(bc, el) + end + solver = LinearSolver() + push!(solver, problem) + push!(solver, bc) +# solver.dump_matrices = true +# solver.name = "3d_hex8" + solver(0.0) + + X = element1("geometry", [1.0, 1.0, 1.0], 0.0) + u = element1("displacement", [1.0, 1.0, 1.0], 0.0) + info("displacement at $X = $u") + # verified using Code Aster. # 2015-12-12-continuum-elasticity/c3d_linear.* - @test isapprox(u[:,7], [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01]) + # [1/36, 1/36, -1/9] + @test isapprox(u, [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01]) end #test_continuum_elasticity_with_surface_load() diff --git a/test/test_mortar.jl b/test/test_mortar.jl index 4c9f493..a856083 100644 --- a/test/test_mortar.jl +++ b/test/test_mortar.jl @@ -16,7 +16,8 @@ using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane, get_edge_intersections, get_points_inside_triangle, clip_polygon, calculate_polygon_centerpoint, project_point_from_plane_to_surface, assemble, - calculate_normal_tangential_coordinates! + calculate_normal_tangential_coordinates!, + is_point_inside_convex_polygon function get_test_2d_model() @@ -455,14 +456,38 @@ end #test_get_points_inside_triangle() -function test_polygon_clipping() +function test_is_point_inside_convex_polygon() + X = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]] + @test is_point_inside_convex_polygon([0.5, 0.5], X) == true + @test is_point_inside_convex_polygon([1.0, 0.5], X) == true + @test is_point_inside_convex_polygon([1.1, 0.5], X) == false + @test is_point_inside_convex_polygon([1.0, 1.0], X) == true + @test is_point_inside_convex_polygon([0.0, 0.3], X) == true + @test is_point_inside_convex_polygon([0.0, -0.000001], X) == false +end + + +function test_polygon_clipping_easy() S = [0 0; 3 0; 0 3]' M = [-1 1; 2 -1/2; 2 2]' P, n = clip_polygon(S, M) @test isapprox(P, [0.0 0.5; 1.0 0.0; 2.0 0.0; 2.0 1.0; 1.25 1.75; 0.0 4/3]') @test isapprox(n, [1 0 1; 1 1 0; 0 1 1]) end -#test_polygon_clipping() + +function test_polygon_clipping_no_clip() + # no clipping at all + S = [-0.125 0.125 0.125 -0.125 + -0.125 -0.125 0.125 0.125] + M = [-0.291667 -0.625 -0.625 -0.291667 + -0.208333 -0.208333 0.125 0.125 ] + P, n = clip_polygon(S, M) + # FIXME: check better. + @test isa(P, Void) + @test isa(n, Void) + +end +#test_polygon_clipping_no_clip() function test_calculate_polygon_centerpoint() @@ -477,7 +502,7 @@ end -function test_assemble_3d_problem() +function test_assemble_3d_problem_tri3() nodes = Vector{Float64}[ [0.0, 0.0, 0.0], [1.0, 0.0, 0.0], @@ -533,7 +558,47 @@ function test_assemble_3d_problem() @test isapprox(stiffness_matrix, B) end -test_assemble_3d_problem() +#test_assemble_3d_problem_tri3() +function test_assemble_3d_problem_quad4() + nodes = Vector{Float64}[ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.1], + [1.0, 0.0, 0.1], + [1.0, 1.0, 0.1], + [0.0, 1.0, 0.1]] +#= + nodes = Vector{Float64}[ + [-1.0, -1.0, 0.0], + [+1.0, -1.0, 0.0], + [+1.0, +1.0, 0.0], + [-1.0, +1.0, 0.0], + [-1.0, -1.0, 0.1], + [+1.0, -1.0, 0.1], + [+1.0, +1.0, 0.1], + [-1.0, +1.0, 0.1]] +=# + mel = Quad4([5, 6, 7, 8]) + mel["geometry"] = Vector{Float64}[nodes[5], nodes[6], nodes[7], nodes[8]] + sel = Quad4([1, 2, 3, 4]) + sel["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3], nodes[4]] + calculate_normal_tangential_coordinates!(sel, 0.0) + sel["master elements"] = Element[mel] + prob = MortarProblem("temperature", 1) + + push!(prob, sel) + stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix) + info("stiffness matrix for this problem:\n$stiffness_matrix") + M = D = 1/36*[4 2 1 2; 2 4 2 1; 1 2 4 2; 2 1 2 4] + B = [D -M] # slave dofs are first in this. + info("expected matrix for this problem:\n$B") + @test isapprox(stiffness_matrix, B) + +end +#test_assemble_3d_problem_quad4() + end diff --git a/test/test_aster_reader.jl b/test/test_preprocess_aster_reader.jl similarity index 100% rename from test/test_aster_reader.jl rename to test/test_preprocess_aster_reader.jl