<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Elasticity solver examples\n",
    "\n",
    "Author(s): <NAME> <<EMAIL>>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# These are internal module functions and not intended to use like this.\n",
    "using JuliaFEM.elasticity_solver\n",
    "using JuliaFEM.xdmf\n",
    "using JuliaFEM.abaqus_reader\n",
    "using Logging\n",
    "Logging.configure(level=DEBUG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using LightXML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2d beam with linear elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2x4 Array{Float64,2}:\n",
       " 0.0  -0.399145  -0.0722858  0.0\n",
       " 0.0  -2.17799   -2.22224    0.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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",
    "  nodalloads = [0 0; 0 0; 0 -2; 0 0]'\n",
    "  dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0]'\n",
    "\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",
    "  la = la*ones(1, 4)\n",
    "  mu = mu*ones(1, 4)\n",
    "  u = zeros(2, 4)\n",
    "  du = zeros(2, 4)\n",
    "\n",
    "  N(xi) = [\n",
    "      (1-xi[1])*(1-xi[2])/4\n",
    "      (1+xi[1])*(1-xi[2])/4\n",
    "      (1+xi[1])*(1+xi[2])/4\n",
    "      (1-xi[1])*(1+xi[2])/4\n",
    "    ]\n",
    "\n",
    "  dNdξ(ξ) = [-(1-ξ[2])/4.0    -(1-ξ[1])/4.0\n",
    "              (1-ξ[2])/4.0    -(1+ξ[1])/4.0\n",
    "              (1+ξ[2])/4.0     (1+ξ[1])/4.0\n",
    "             -(1+ξ[2])/4.0     (1-ξ[1])/4.0]\n",
    "\n",
    "  ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
    "  iweights = [1 1 1 1]\n",
    "\n",
    "    return (X, u, du, elmap, nodalloads, dirichletbc,\n",
    "     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",
    "    \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",
    "    end\n",
    "end\n",
    "\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "u3d = [u; 0 0 0 0] # extend to 3d vector field\n",
    "X3d = [X; 0 0 0 0]\n",
    "elmap2 = [0x5; elmap]'';"
   ]<commit_msg>unfinished, some ideas<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Elasticity solver examples\n",
    "\n",
    "Author(s): <NAME> <<EMAIL>>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Docile: upgrading cache to 0.0.2.\n"
     ]
    },
    {
     "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",
      "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"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# These are internal module functions and not intended to use like this.\n",
    "using JuliaFEM.elasticity_solver\n",
    "using JuliaFEM.xdmf\n",
    "using JuliaFEM.abaqus_reader\n",
    "using Logging\n",
    "Logging.configure(level=DEBUG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "using LightXML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converged\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2x4 Array{Float64,2}:\n",
       " 0.0  -0.399145  -0.0722858  0.0\n",
       " 0.0  -2.17799   -2.22224    0.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "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",
    "  nodalloads = [0 0; 0 0; 0 -2; 0 0]'\n",
    "  dirichletbc = [0 0; NaN NaN; NaN NaN; 0 0]'\n",
    "\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",
    "  la = la*ones(1, 4)\n",
    "  mu = mu*ones(1, 4)\n",
    "  u = zeros(2, 4)\n",
    "  du = zeros(2, 4)\n",
    "\n",
    "  N(xi) = [\n",
    "      (1-xi[1])*(1-xi[2])/4\n",
    "      (1+xi[1])*(1-xi[2])/4\n",
    "      (1+xi[1])*(1+xi[2])/4\n",
    "      (1-xi[1])*(1+xi[2])/4\n",
    "    ]\n",
    "\n",
    "  dNdξ(ξ) = [-(1-ξ[2])/4.0    -(1-ξ[1])/4.0\n",
    "              (1-ξ[2])/4.0    -(1+ξ[1])/4.0\n",
    "              (1+ξ[2])/4.0     (1+ξ[1])/4.0\n",
    "             -(1+ξ[2])/4.0     (1-ξ[1])/4.0]\n",
    "\n",
    "  ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
    "  iweights = [1 1 1 1]\n",
    "\n",
    "    return (X, u, du, elmap, nodalloads, dirichletbc,\n",
    "     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",
    "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",
    "    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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "u3d = [u; 0 0 0 0] # extend to 3d vector field\n",
    "X3d = [X; 0 0 0 0]\n",
    "elmap2 = [0x5; elmap]'';"
   ]
{'fraction_non_alphanumeric': 0.22116283417461208, 'fraction_numerical': 0.034835171683180656, 'mean_word_length': 1.881149012567325, 'pattern_counts': {'":': 108, '<': 11, '<?xml version=': 0, '>': 23, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3627108', 'n_tokens_mistral': 6934, 'n_tokens_neox': 6587, 'n_words': 1342}