mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-22 18:52:16 +00:00
unfinished, some ideas
This commit is contained in:
Jukka Aho
@@ -16,6 +16,13 @@
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Docile: upgrading cache to 0.0.2.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
@@ -25,7 +32,10 @@
|
||||
"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"
|
||||
"Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
|
||||
"Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n",
|
||||
"Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n",
|
||||
"Warning: requiring \"JuliaFEM\" did not define a corresponding module.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -50,7 +60,7 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"execution_count": 2,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
@@ -66,6 +76,28 @@
|
||||
"## 2d beam with linear elements"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"Array{Float64,2}"
|
||||
]
|
||||
},
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"typeof(Float64[1 2; 2 3])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
@@ -94,6 +126,64 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"\"\"\"\n",
|
||||
"#### Parameters\n",
|
||||
"X : Array{Float64, 2}\n",
|
||||
" \n",
|
||||
"\"\"\"\n",
|
||||
"function solve_elasticity_increment!(X, u, du, elmap, nodalloads,\n",
|
||||
" dirichletbc, lambda, mu, N, dNdchi, ipoints,\n",
|
||||
" iweights)\n",
|
||||
" if length(size(elmap)) == 1\n",
|
||||
" # quick hack for just one element\n",
|
||||
" elmap = elmap''\n",
|
||||
" end\n",
|
||||
" nelnodes, nelements = size(elmap)\n",
|
||||
" dim, nnodes = size(u)\n",
|
||||
" dofs = dim*nelnodes\n",
|
||||
"\n",
|
||||
" Imat = Int64[]\n",
|
||||
" Jmat = Int64[]\n",
|
||||
" Vmat = Float64[]\n",
|
||||
" Ivec = Int64[]\n",
|
||||
" Vvec = Float64[]\n",
|
||||
"\n",
|
||||
" # FIXME: different number of nodes/element\n",
|
||||
" R = zeros(dim, nelnodes)\n",
|
||||
" Kt = zeros(dofs, dofs)\n",
|
||||
"\n",
|
||||
" # this can be parallelized\n",
|
||||
" for i in 1:nelements\n",
|
||||
" eldofs = elmap[:,i]\n",
|
||||
" calc_local_matrices!(X[:, eldofs], u[:, eldofs], R, Kt, N, dNdchi,\n",
|
||||
" lambda[eldofs], mu[eldofs], ipoints, iweights)\n",
|
||||
" assemble!(Kt, eldofs, Imat, Jmat, Vmat)\n",
|
||||
" assemble!(R, eldofs, Ivec, Vvec)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # add additional neumann boundary conditions to force vector\n",
|
||||
" for (i, nodal_load) in enumerate(nodalloads)\n",
|
||||
" if nodal_load == 0\n",
|
||||
" continue\n",
|
||||
" end\n",
|
||||
" push!(Ivec, i)\n",
|
||||
" push!(Vvec, -nodal_load)\n",
|
||||
" end\n",
|
||||
"\n",
|
||||
" # Create sparse matrix and vector\n",
|
||||
" A = sparse(Imat, Jmat, Vmat)\n",
|
||||
" b = sparsevec(Ivec, Vvec)\n",
|
||||
"\n",
|
||||
" # Remove dirichlet boundary conditions\n",
|
||||
" free_dofs = find(isnan(dirichletbc))\n",
|
||||
" #Imat, Jmat, Vmat = eliminate_boundary_conditions(dirichletbc, Imat, Jmat, Vmat)\n",
|
||||
" b = b[free_dofs]\n",
|
||||
" A = A[free_dofs, free_dofs]\n",
|
||||
"\n",
|
||||
" # solution\n",
|
||||
" du[free_dofs] = lufact(A) \\ -full(b)\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"function one_elem_fixture()\n",
|
||||
" X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
|
||||
" elmap = [1; 2; 3; 4]\n",
|
||||
@@ -130,17 +220,137 @@
|
||||
" la, mu, N, dNdξ, ipoints, iweights)\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
" la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()\n",
|
||||
"#(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
"# la, mu, N, dNdξ, ipoints, iweights) = one_elem_fixture()\n",
|
||||
"function solve_one_element()\n",
|
||||
"\n",
|
||||
" m = new_model()\n",
|
||||
"\n",
|
||||
" # create nodes separately ...\n",
|
||||
" n1 = new_node()\n",
|
||||
" set_node_id(n1, 1)\n",
|
||||
" set_node_coords(n1, [0.0, 0.0])\n",
|
||||
"\n",
|
||||
" n2 = new_node()\n",
|
||||
" set_node_id(n2, 2)\n",
|
||||
" set_node_coords(n2, [10.0, 0.0])\n",
|
||||
" \n",
|
||||
"for i=1:10\n",
|
||||
" JuliaFEM.elasticity_solver.solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
" la, mu, N, dNdξ, ipoints, iweights)\n",
|
||||
" u += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged\")\n",
|
||||
" break\n",
|
||||
" n3 = new_node()\n",
|
||||
" set_node_id(n3, 3)\n",
|
||||
" set_node_coords(n3, [10.0, 1.0])\n",
|
||||
"\n",
|
||||
" n4 = new_node()\n",
|
||||
" set_node_id(n4, 4)\n",
|
||||
" set_node_coords(n4, [0.0, 1.0])\n",
|
||||
"\n",
|
||||
" nodes = [n1, n2, n3, n4]\n",
|
||||
" add_nodes(m, \"NALL\", nodes)\n",
|
||||
"\n",
|
||||
" # .. or create somewhat simpler syntax\n",
|
||||
" # nodes = Dict(1 => [0.0, 0.0], 2 => [10.0, 0.0], 3 => [10.0, 1.0], 4 => [0.0, 1.0])\n",
|
||||
" # add_nodes(m, \"NALL\", nodes)\n",
|
||||
"\n",
|
||||
" # create element (hard way)\n",
|
||||
" e = new_element()\n",
|
||||
" set_element_id(1)\n",
|
||||
" set_node_ids(e, [1, 2, 3, 4])\n",
|
||||
"\n",
|
||||
" # we can set function spaces and integration schema for element-wise ...\n",
|
||||
" set_function_space(e, \"Lagrange(1)\") # use linear Lagrange function space to approximate unknown field\n",
|
||||
" set_integration_schema(e, \"FPG4\") # use four integration points\n",
|
||||
" # or for model as a \"default value\"\n",
|
||||
" # set_function_space(m, \"Lagrange\", 1)\n",
|
||||
" # set_integration_schema(m, \"FPG4\")\n",
|
||||
"\n",
|
||||
" # set necessary material parameters\n",
|
||||
" E = 90\n",
|
||||
" nu = 0.25\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",
|
||||
" # we can add attribute for each element node ...\n",
|
||||
" #add_attribute(e, \"lambda\", [la, la, la, la])\n",
|
||||
" #add_attribute(e, \"mu\", [mu, mu, mu, mu])\n",
|
||||
" # or just one for each element ...\n",
|
||||
" add_attribute(e, \"lambda\", lambda)\n",
|
||||
" add_attribute(e, \"mu\", mu)\n",
|
||||
" # .. or just for model as a default value\n",
|
||||
" # add_attribute(m, \"lambda\", lambda)\n",
|
||||
" # add_attribute(m, \"mu\", mu)\n",
|
||||
" # note that we don't assign attributes to nodes here so we can describe discontinous fields\n",
|
||||
" # in attributes\n",
|
||||
"\n",
|
||||
" # or we can use convenient syntax\n",
|
||||
" # e = new_element(element_id=1, node_ids=[1, 2, 3, 4], function_space=\"Lagrange(1)\",\n",
|
||||
" # integration_schema=\"FPG4\", attributes=Dict(\"lambda\" => lambda, \"mu\" => mu))\n",
|
||||
"\n",
|
||||
" elements = [e]\n",
|
||||
" add_elements(m, \"EALL\", elements)\n",
|
||||
"\n",
|
||||
" # boundary conditions are always assignet to sets\n",
|
||||
"\n",
|
||||
" # element boundary conditions\n",
|
||||
" \n",
|
||||
" # add load for element surface S1 in normal-tangential coordinate system\n",
|
||||
" # add_attribute(e, \"displacement S1 normal load\", 1)\n",
|
||||
" # add load for element surface S1 from -1 to -2\n",
|
||||
" # add_attribute(e, \"displacement S1 load\", [-1, -2])\n",
|
||||
"\n",
|
||||
" # nodal boundary conditions\n",
|
||||
" # add neumann boundary condition to node 3, in y direction\n",
|
||||
" add_attribute(n3, \"displacement 2 load\", -2)\n",
|
||||
" # add dirichlet boundary condition to nodes 1 and 2 (encastre)\n",
|
||||
" add_attribute(n1, \"displacement\", [0.0, 0.0])\n",
|
||||
" # or\n",
|
||||
" # add_attribute(n1, \"displacement 1\", 0.0)\n",
|
||||
" # add_attribute(n1, \"displacement 2\", 0.0)\n",
|
||||
" add_attribute(n2, \"displacement\", [0.0, 0.0])\n",
|
||||
" # or, for example, add dirichlet boundary conditions in directions 1 and 3 for node\n",
|
||||
" #add_attribute(n2, \"displacement 1,3\", [0.0, 1.0])\n",
|
||||
" # or fix all dofs of a nodes in nodeset \"SUPPORT\"\n",
|
||||
" # support_bc = add_nodeset(m, \"SUPPORT\", [n1, n2])\n",
|
||||
" # add_attribute(support_bc, \"displacement\", 0)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" # Everything is very general so far. Datamodel is well defined and in this point\n",
|
||||
" # we can save or load it to disk\n",
|
||||
"\n",
|
||||
" # save_model(m, \"mymodel\") # save model to xml/h5 (Xdmf)\n",
|
||||
" # m = load_model(\"mymodel\") # load model from file\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # END OF MODEL DEFINITON\n",
|
||||
" \n",
|
||||
" # Now we kick in elasticity iterations and solve displacement\n",
|
||||
" # field but of course it could be something else too\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" # m1, m2 = make_domain_decomposition(parts=2, keep_in_one_domain=[\"CONTACT_BOUNDARY\"])\n",
|
||||
"\n",
|
||||
"\n",
|
||||
" \n",
|
||||
" # In newton iteration, there might be situations where we just want to update RHS\n",
|
||||
" # like in radiation problems, it makes no sense to update stiffness matrix in that\n",
|
||||
" # case. For this reason local element matrix and force vector can be updated both\n",
|
||||
" # or just one of them. This time we don't have anything nonlinear in RHS so we can\n",
|
||||
" # calculate it outside of iteration loop\n",
|
||||
" p = new_problem(\"displacement\")\n",
|
||||
" \n",
|
||||
" RHS = zeros(length(nodes)*dofs)\n",
|
||||
" \n",
|
||||
" for i=1:10\n",
|
||||
" \n",
|
||||
" solve_elasticity_increment!(X, u, du, elmap, nodalloads, dirichletbc,\n",
|
||||
" la, mu, N, dNdξ, ipoints, iweights)\n",
|
||||
" u += du\n",
|
||||
" if norm(du) < 1.0e-9\n",
|
||||
" println(\"Converged\")\n",
|
||||
" break\n",
|
||||
" end\n",
|
||||
" end\n",
|
||||
" return u\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"u"
|
||||
|
||||
Reference in New Issue
Block a user