WaveProp1D-Veri2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "8e060377",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/aangniu/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/jitcompiler.py:234: QuadratureRepresentationDeprecationWarning: \n",
      "*** ===================================================== ***\n",
      "*** FFC: quadrature representation is deprecated! It will ***\n",
      "*** likely be removed in 2018.2.0 release. Use uflacs     ***\n",
      "*** representation instead.                               ***\n",
      "*** ===================================================== ***\n",
      "  issue_deprecation_warning()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------------\n",
      "Calling FFC just-in-time (JIT) compiler, this may take some time.\n",
      "  Ignoring precision in integral metadata compiled using quadrature representation. Not implemented.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/aangniu/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/ffc/quadrature/quadraturerepresentation.py:58: QuadratureRepresentationDeprecationWarning: \n",
      "*** ===================================================== ***\n",
      "*** FFC: quadrature representation is deprecated! It will ***\n",
      "*** likely be removed in 2018.2.0 release. Use uflacs     ***\n",
      "*** representation instead.                               ***\n",
      "*** ===================================================== ***\n",
      "  issue_deprecation_warning()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution misfit: 0.011882819221281066 16999.11449720282  dt: 0.05\n",
      "---------------------------------\n",
      "Solution misfit: 0.005831806396599088 8758.676657696735  dt: 0.025\n",
      "---------------------------------\n",
      "Solution misfit: 0.0028845804608949597 4431.406443301782  dt: 0.0125\n",
      "---------------------------------\n",
      "Solution misfit: 0.0014339345012653112 2231.069521820808  dt: 0.00625\n",
      "---------------------------------\n",
      "Solution misfit: 0.0007148112227730335 1118.029826373461  dt: 0.003125\n",
      "---------------------------------\n",
      "Solution misfit: 0.0003568580618774018 559.8039444706573  dt: 0.0015625\n"
     ]
    }
   ],
   "source": [
    "from dolfin import *\n",
    "import numpy as np\n",
    "parameters[\"form_compiler\"][\"representation\"] = 'quadrature'\n",
    "import warnings\n",
    "from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning\n",
    "warnings.simplefilter(\"once\", QuadratureRepresentationDeprecationWarning)\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Form compiler options\n",
    "parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
    "parameters[\"form_compiler\"][\"optimize\"] = True\n",
    "\n",
    "\n",
    "\n",
    "N_refines = np.array([10,20,40,80,160,320])\n",
    "dts = np.zeros((N_refines.size,1))\n",
    "errors_result = np.zeros((N_refines.size,1))\n",
    "i_ref = 0\n",
    "\n",
    "for N_refine in N_refines:\n",
    "    # Define mesh\n",
    "    N_mesh = N_refine\n",
    "    mesh = UnitIntervalMesh(N_mesh)\n",
    "    #mesh = RectangleMesh(Point(0., 0.), Point(1., 0.1), 50, 5)\n",
    "\n",
    "    tol = 1E-12\n",
    "    # Sub domain for clamp at left end\n",
    "    def left(x, on_boundary):\n",
    "        return near(x[0], 0.,tol) and on_boundary\n",
    "\n",
    "    # Sub domain for rotation at right end\n",
    "    def right(x, on_boundary):\n",
    "        return near(x[0], 1.,tol) and on_boundary\n",
    "\n",
    "    def left_and_right(x, on_boundary):\n",
    "        return on_boundary\n",
    "\n",
    "    # Create mesh function over the cell facets\n",
    "    boundary_subdomains = MeshFunction(\"size_t\", mesh, mesh.topology().dim() - 1)\n",
    "    boundary_subdomains.set_all(0)\n",
    "    force_boundary = AutoSubDomain(right)\n",
    "    force_boundary.mark(boundary_subdomains, 3)\n",
    "\n",
    "    # Define measure for boundary condition integral\n",
    "    dss = ds(subdomain_data=boundary_subdomains)\n",
    "\n",
    "    E = Constant(70e3)\n",
    "    nu = Constant(0.3)\n",
    "    lmbda = E*nu/(1+nu)/(1-2*nu)\n",
    "    mu = E/2./(1+nu)\n",
    "    sig0 = Constant((1e0))  # yield strength\n",
    "    Et = E/100.  # tangent modulus\n",
    "    H = E*Et/(E-Et)  # hardening modulus\n",
    "    rho = Constant(10e4)\n",
    "\n",
    "    Nitermax, tol = 200, 1e-8  # parameters of the Newton-Raphson procedure\n",
    "    Nincr = N_refine*2\n",
    "    # Nincr = 100\n",
    "    load_steps = np.linspace(0, 1, Nincr+1)[1:]\n",
    "\n",
    "    dt = Constant(load_steps[1]-load_steps[0])\n",
    "    dt_ = float(dt)\n",
    "\n",
    "    # Newmark-beta method\n",
    "    betaN = Constant(0.25)\n",
    "    gammaN = Constant(0.5)\n",
    "\n",
    "    deg_u = 2\n",
    "    deg_stress = 2\n",
    "    V = FunctionSpace(mesh, \"CG\", deg_u)\n",
    "    We = VectorElement(\"Quadrature\", mesh.ufl_cell(), degree=deg_stress, dim=1, quad_scheme='default')\n",
    "    W = FunctionSpace(mesh, We)\n",
    "    W0e = FiniteElement(\"Quadrature\", mesh.ufl_cell(), degree=deg_stress, quad_scheme='default')\n",
    "    W0 = FunctionSpace(mesh, W0e)\n",
    "    # Define function space for stresses\n",
    "    Vsig = VectorFunctionSpace(mesh, \"DG\", 0)\n",
    "\n",
    "    sig = Function(W)\n",
    "    sig_old = Function(W)\n",
    "    beta = Function(W0)\n",
    "    beta1 = Function(V)\n",
    "    gamma0 = Function(W0)\n",
    "    gamma0.interpolate(Constant(1e-6))\n",
    "    u = Function(V, name=\"Total displacement\")\n",
    "    u_pre = Function(V)\n",
    "    u_old = Function(V)\n",
    "    u_d_old = Function(V)\n",
    "    u_dd_old = Function(V)\n",
    "    du = Function(V, name=\"Iteration correction\")\n",
    "    Du = Function(V, name=\"Current increment\")\n",
    "    v = TrialFunction(V)\n",
    "    u_ = TestFunction(V)\n",
    "\n",
    "    u_theo = Function(V)\n",
    "\n",
    "    sig_out = Function(Vsig, name=\"sigma\")\n",
    "    eps_out = Function(Vsig, name=\"epsilon\")\n",
    "\n",
    "    # Set up boundary condition at left end\n",
    "    # bc_expression = Expression( (\"t <= tc ? x[0]*1e-6*pow(t,1) : x[0]*0.5e-6\"), tc=0.5, t=0, degree=2) # Constant(0.0)\n",
    "    # t*x[0]*(1-t*x[0]/L)/L -BC1\n",
    "    # t*x[0]*(1-t*x[0]*x[0]/L)/L -BC2\n",
    "    # t*t*x[0]*(1-x[0]*x[0]*t/L)/L - BC3\n",
    "    # \"t*t*x[0]*(1-x[0]*x[0]/L)/L\" - 0 BC\n",
    "    bc_expression = Expression((\"t*t*x[0]*(1-x[0]*x[0]*t/L)/L\"), L=1, t=0, degree=2)\n",
    "    bc = DirichletBC(V, bc_expression, left_and_right)    \n",
    "\n",
    "    zero = Constant(0.0)\n",
    "    bc_newton = DirichletBC(V, zero, left_and_right)\n",
    "\n",
    "    # bc = [DirichletBC(V.sub(1), 0, facets, 1), DirichletBC(V.sub(0), 0, facets, 3)]\n",
    "    n = FacetNormal(mesh)\n",
    "\n",
    "    q_lim = float(2/sqrt(3)*sig0*mu*2)\n",
    "    # loading = Expression(\"q*t\", q=q_lim, t=0, degree=2)\n",
    "    loading = Expression((\"t <= tc ? q*t*(1-t)*4 : 0\"), q=0.2*q_lim, t=0, tc=0.7, degree=2)\n",
    "\n",
    "    # DEFINING SOURCE TERM\n",
    "    # source = Expression((\"t <= tc && x[0]<0.55 && x[0]>0.45  ? q*t*(1-t)*4 : 0\"), q=0.2*q_lim, t=0, tc=1.2, degree=2)\n",
    "    # source = Constant(0.0)\n",
    "\n",
    "    # MMS - Nonlinear\n",
    "    # \"(2*eps0*eps0*(G0*t*t-rho*x[0]*x[0])*pow(L,4)-4*pow(L,3)*eps0*rho*t*x[0]*x[0]+8*t*t*(eps0*x[0]-1/4)*rho*pow(x[0],2)*L*L+8*L*rho*pow(t,3)*pow(x[0],3)-8*rho*pow(t,4)*pow(x[0],4))/(L*L*pow(L*L*eps0-2*t*t*x[0]+L*t,2))\"\n",
    "    # BC3 - Nonlinear\n",
    "    # (2*(pow(L,5)*eps0*eps0*rho+3*t*eps0*eps0*(G0*t*t-rho*x[0]*x[0])*pow(L,4)+9*G0*L*L*pow(t,7)-54*G0*L*pow(t,8)*x[0]*x[0]+81*G0*pow(t,9)*pow(x[0],4)))*x[0]/(pow(eps0,2)*pow(L,6))\n",
    "\n",
    "    # MMS - linear\n",
    "    # (2*G0*t*t-2*rho*x[0]*x[0])/L/L - BC1\n",
    "    # 6*x[0]*(G0*t*t-(1/3)*rho*x[0]*x[0])/L/L - BC2\n",
    "    # 2*x[0]*(rho*(-3*t*x[0]*x[0]+L)+3*G0*pow(t,3))/L/L - BC3\n",
    "    # (2*(rho*(-x[0]*x[0]+L)+3*G0*t*t))*x[0]/L/L - 0 BC case\n",
    "\n",
    "    source = Expression(\\\n",
    "        (\"(2*(pow(L,5)*eps0*eps0*rho+3*t*eps0*eps0*(G0*t*t-rho*x[0]*x[0])*pow(L,4)+9*G0*L*L*pow(t,7)-54*G0*L*pow(t,8)*x[0]*x[0]+81*G0*pow(t,9)*pow(x[0],4)))*x[0]/(pow(eps0,2)*pow(L,6))\")\\\n",
    "        ,G0=float(mu), eps0=float(sig0), rho=float(rho), t=0, L=1, degree=2)\n",
    "\n",
    "    # (2*eps0*eps0*(G0*t*t-rho*x*x)*pow(L,4)\n",
    "    #  -4*pow(L,3)*eps0*rho*t*x*x\n",
    "    #  +8*x*x*t*t*rho*(x*eps0-1/4)*L*L\n",
    "    #  +8*L*rho*pow(t,3)*pow(x,3)\n",
    "    #  -8*rho*pow(t,4)*pow(x,4)\n",
    "    # )\n",
    "    # /(L*L*pow(L*L*eps0-2*t*2*x+L*t,2))\n",
    "\n",
    "\n",
    "    def F_ext(v):\n",
    "        return dot(loading,v)*dss(3)\n",
    "        # return loading*dot(n, v)*dss(3)\n",
    "\n",
    "    # Update formula for acceleration\n",
    "    # a = 1/(2*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) - (1-2*beta)*a0)\n",
    "    def update_a(u, u_old, v_old, a_old, ufl=True):\n",
    "        if ufl:\n",
    "            dt_ = dt\n",
    "            betaN_ = betaN\n",
    "        else:\n",
    "            dt_ = float(dt)\n",
    "            betaN_ = float(betaN)\n",
    "        return (u-u_old-dt_*v_old)/betaN_/dt_**2 - (1-2*betaN_)/2/betaN_*a_old\n",
    "\n",
    "    # Update formula for velocity\n",
    "    # v = dt * ((1-gamma)*a0 + gamma*a) + v0\n",
    "    def update_v(a, u_old, v_old, a_old, ufl=True):\n",
    "        if ufl:\n",
    "            dt_ = dt\n",
    "            gamma_ = gammaN\n",
    "        else:\n",
    "            dt_ = float(dt)\n",
    "            gamma_ = float(gammaN)\n",
    "        return v_old + dt_*((1-gamma_)*a_old + gamma_*a)\n",
    "\n",
    "    def update_fields(u, u_old, v_old, a_old):\n",
    "        \"\"\"Update fields at the end of each time step.\"\"\"\n",
    "\n",
    "        # Get vectors (references)\n",
    "        u_vec, u0_vec  = u.vector(), u_old.vector()\n",
    "        v0_vec, a0_vec = v_old.vector(), a_old.vector()\n",
    "\n",
    "        # use update functions using vector arguments\n",
    "        a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)\n",
    "        v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)\n",
    "\n",
    "        # Update (u_old <- u)\n",
    "        v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec\n",
    "        u_old.vector()[:] = u.vector()\n",
    "\n",
    "    def eps(v):\n",
    "        e = 0.5*grad(v)\n",
    "        return e\n",
    "\n",
    "    def sigma(eps_el):\n",
    "        return 2*mu*eps_el\n",
    "\n",
    "    def sigma_out(r):\n",
    "        return 2.0*mu*sym(grad(r))\n",
    "\n",
    "    def proj_sig(old_u,old_sig,old_Du):\n",
    "        eps_old = eps(old_u)\n",
    "        beta = sig0/(abs(eps_old[0])+sig0)\n",
    "        # beta = 1\n",
    "        beta = 1+3*(2*eps_old[0]/sig0)*(2*eps_old[0]/sig0)\n",
    "        # beta = gamma0\n",
    "        # beta /= (eps_old[0]+as_tensor(gamma0[0]))\n",
    "        return 2*mu*beta*eps(old_Du)+old_sig, beta\n",
    "\n",
    "    def sigma_tang(e):\n",
    "        return 2*mu*beta*e\n",
    "\n",
    "    metadata = {\"quadrature_degree\": deg_stress, \"quadrature_scheme\": \"default\"}\n",
    "    dxm = dx(metadata=metadata)\n",
    "\n",
    "    # a_Newton = inner(eps(v), sigma_tang(eps(u_)))*dxm\n",
    "\n",
    "    # a_Newton = inner(eps(u_), sigma_tang(eps(v)))*dxm\n",
    "\n",
    "    a_Newton = rho*inner(u_,v)/(betaN*dt*dt)*dxm + inner(2*eps(u_), sigma_tang(eps(v)) )*dxm\n",
    "    # F_ext(u_) + \\\n",
    "    res = (-rho)*inner(u_,u)/(betaN*dt*dt)*dxm + \\\n",
    "          -inner(2*eps(u_), sig)*dxm + \\\n",
    "          F_ext(u_) + \\\n",
    "          inner(u_,source)*dxm + \\\n",
    "          rho*inner(u_old+dt*u_d_old,u_)*1./(betaN*dt*dt)*dxm + \\\n",
    "          (1.-2*betaN)/(2*betaN)*rho*inner(u_dd_old,u_)*dxm\n",
    "\n",
    "    def local_project(v, V, u=None):\n",
    "        dv = TrialFunction(V)\n",
    "        v_ = TestFunction(V)\n",
    "        a_proj = inner(dv, v_)*dxm\n",
    "        b_proj = inner(v, v_)*dxm\n",
    "        solver = LocalSolver(a_proj, b_proj)\n",
    "        solver.factorize()\n",
    "        if u is None:\n",
    "            u = Function(V)\n",
    "            solver.solve_local_rhs(u)\n",
    "            return u\n",
    "        else:\n",
    "            solver.solve_local_rhs(u)\n",
    "            return\n",
    "\n",
    "    file_results = XDMFFile(\"nonlinearElas1DVeri__results.xdmf\")\n",
    "    file_results.parameters[\"flush_output\"] = True\n",
    "    file_results.parameters[\"functions_share_mesh\"] = True\n",
    "\n",
    "    results = np.zeros((Nincr+1, 5))\n",
    "    results_1d = np.zeros((Nincr+1, N_mesh))\n",
    "\n",
    "    # Initial conditions\n",
    "    u_old.interpolate( Expression((\"0\"), degree=2) )\n",
    "    u_d_old.interpolate( Expression((\"0\"), degree=2) )\n",
    "    # u.assign(u_old)\n",
    "\n",
    "    for (i, t) in enumerate(load_steps):\n",
    "        loading.t = t\n",
    "        bc_expression.t = t\n",
    "        source.t = t\n",
    "\n",
    "        u_pre.assign(u)\n",
    "\n",
    "        # Prepare for U_n+1^0 from U_n\n",
    "        bc.apply(u.vector())\n",
    "        # print('U at the left and right:',u(0), \" \", u(1))\n",
    "        # print('Source at the left and right:',source(0), \" \", source(1))\n",
    "        # -- Compute components in F(U_n+1^0)\n",
    "        Du.assign(u-u_pre)\n",
    "        sig_, beta_ = proj_sig(u_pre,sig_old,Du)\n",
    "        local_project(sig_, W, sig)\n",
    "        local_project(beta_, W0, beta)\n",
    "        A, Res = assemble_system(a_Newton, res, bc_newton)\n",
    "\n",
    "        nRes0 = Res.norm(\"l2\")\n",
    "        nRes = nRes0\n",
    "        # print(\"Increment:\", str(i+1),\"   Initial Residual:\", nRes0)\n",
    "        niter = 0\n",
    "\n",
    "        while (nRes/nRes0 > tol and nRes0>1e-10) and niter < Nitermax:\n",
    "            solve(A, du.vector(), Res, \"mumps\")\n",
    "            Du.assign(Du+du)\n",
    "            u.assign(u+du) # u^{k+1}_{n+1}\n",
    "            sig_, beta_ = proj_sig(u_pre,sig_old,Du)\n",
    "            local_project(sig_, W, sig)\n",
    "            local_project(beta_, W0, beta)\n",
    "            A, Res = assemble_system(a_Newton, res, bc_newton)\n",
    "            nRes = Res.norm(\"l2\")\n",
    "            # print(\"    Residual:\", nRes,\" \",niter)\n",
    "            niter += 1\n",
    "        u.assign(u_pre+Du)\n",
    "\n",
    "        update_fields(u, u_old, u_d_old, u_dd_old)\n",
    "\n",
    "        sig_old.assign(sig)\n",
    "\n",
    "    # ----------------\n",
    "    # Post-processing\n",
    "    # ----------------\n",
    "    #\n",
    "    # Inside the incremental loop, the displacement and plastic strains are exported\n",
    "    # at each time increment, the plastic strain must first be projected onto the\n",
    "    # previously defined DG FunctionSpace. We also monitor the value of the cylinder\n",
    "    # displacement on the inner boundary. The load-displacement curve is then plotted::\n",
    "\n",
    "        local_project(sig, Vsig, sig_out)\n",
    "        local_project(eps(u), Vsig, eps_out)\n",
    "\n",
    "        file_results.write(u, t)\n",
    "        for i_x in range(N_mesh):\n",
    "            results_1d[i+1,i_x] = u(i_x/(N_mesh-1))\n",
    "        tc = 0.7;\n",
    "        if t<tc:\n",
    "            # results[i+1, :] = (u(1),t,4*t*(1-t))\n",
    "            results[i+1, :] = (u(1), sig_out(1),t,4*t*(1-t),eps_out(1))\n",
    "        else:\n",
    "            results[i+1, :] = (u(1), sig_out(1),t,0,eps_out(1))\n",
    "\n",
    "    print(\"---------------------------------\")\n",
    "    \n",
    "    sig_expression = Expression((\"G0*t*t*(-3*t*x[0]*x[0]+L)*(9*pow(t,6)*pow(x[0],4)-6*L*pow(t,5)*x[0]*x[0]+pow(L,4)*eps0*eps0+L*L*pow(t,4))/(eps0*eps0*pow(L,6))\",)\\\n",
    "                                ,G0=float(mu), eps0=float(sig0), rho=float(rho), L=1, t=1, degree=2)\n",
    "    sig_theo = Function(Vsig)\n",
    "    sig_theo.interpolate(sig_expression)\n",
    "    \n",
    "    u_theo.interpolate(bc_expression)\n",
    "    error = (u_theo - u)**2*dxm\n",
    "    Errors_theo = sqrt(abs(assemble(error)))\n",
    "    error = (sig_theo - sig)**2*dxm\n",
    "    Errors_theo1 = sqrt(abs(assemble(error)))\n",
    "    errors_result[i_ref] = Errors_theo\n",
    "    dts[i_ref] = 1/Nincr\n",
    "    i_ref += 1\n",
    "\n",
    "    print(\"Solution misfit:\", Errors_theo,Errors_theo1,\" dt:\",1/Nincr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "7e561f80",
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'scipy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_52623/2961353623.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m# Convergence study\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moptimize\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mcurve_fit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'scipy'"
     ]
    }
   ],
   "source": [
    "# Convergence study\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "def func(x, a, b):\n",
    "    return a*x + b\n",
    "\n",
    "popt, pcov = curve_fit(func,np.log(dts),np.log(errors_result))\n",
    "\n",
    "plt.loglog(dts,errors_result,'*')\n",
    "plt.axis('equal')\n",
    "plt.xlabel('dt')\n",
    "plt.ylabel('RMS Errors')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f2e49b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "\n",
    "from matplotlib.animation import FuncAnimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37e95632",
   "metadata": {},
   "outputs": [],
   "source": [
    "# MPII\n",
    "# 2*float(mu)*float(sig0)*np.log(1+abs(results[:N_pt, 4])/float(sig0))\\\n",
    "#           *abs(results[:N_pt, 4])/results[:N_pt, 4]\n",
    "\n",
    "# Simple case\n",
    "# 2*float(mu)*results[:N_pt, 4]*(1+ (2*results[:N_pt, 4]/float(sig0))**2 )\n",
    "\n",
    "N_pt = 500;\n",
    "stress_theo = 2*float(mu)*results[:N_pt, 4]*(1+ (2*results[:N_pt, 4]/float(sig0))**2 )\n",
    "stress_lin = 2*float(mu)*results[:N_pt, 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "263a0c08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.subplot(1,4,1)\n",
    "plt.plot(results[:, 0], results[:, 1], \"-*\")\n",
    "plt.plot(results[:, 0], results[:, 3]*0.2*q_lim, \".\")\n",
    "plt.xlabel(\"Displacement\")\n",
    "plt.ylabel(\"$\\sigma_{xx,yy}$ at the right end\")\n",
    "plt.legend([\"$\\sigma_{yy}$\",\"$\\sigma_{xx}$\",\"Loads\"])\n",
    "plt.subplot(1,4,2)\n",
    "plt.plot(results[:, 2], results[:, 1], \"-.\")\n",
    "plt.plot(results[:, 2], results[:, 3]*0.2*q_lim, \".\")\n",
    "plt.xlabel(\"Loading steps\")\n",
    "plt.ylabel(\"$\\sigma_{xx}$ at the right end\")\n",
    "plt.legend([\"$\\sigma_{yy}$\",\"Loads\"])\n",
    "plt.subplot(1,4,3)\n",
    "plt.plot(results[:, 2], results[:, 0], \"-.\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Displacement at the right end\")\n",
    "plt.subplot(1,4,4)\n",
    "\n",
    "plt.plot(results[:N_pt, 4], results[:N_pt, 1], \".-\")\n",
    "plt.plot(results[:N_pt, 4], stress_theo, \"-\")\n",
    "plt.plot(results[:N_pt, 4], stress_lin, \"-\")\n",
    "plt.xlabel(\"Strain and the right end\")\n",
    "plt.ylabel(\"Stress at the right end\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d02305e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_grid = np.linspace(0,1,N_mesh)\n",
    "u_theo = x_grid/1*(1-x_grid**2/1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "341bd088",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "line1, = ax.plot([])\n",
    "line2, = ax.plot([])\n",
    "\n",
    "ax.set_xlim(0,1)\n",
    "ax.set_ylim(-1e0,1e0)\n",
    "\n",
    "ax.set_xlabel('x')\n",
    "ax.set_ylabel('Displacement')\n",
    "\n",
    "def animate(frame_num):\n",
    "    # ax.clear()\n",
    "    line1.set_data((np.linspace(0,1,N_mesh), results_1d[frame_num, :]))\n",
    "    line2.set_data(x_grid,u_theo)\n",
    "    return line1, line2\n",
    "\n",
    "ani = FuncAnimation(fig, animate, frames=Nincr, interval=10)\n",
    "\n",
    "\n",
    "plt.show()\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fenicsproject",
   "language": "python",
   "name": "fenicsproject"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}