WaveProp1D-Veri3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
      "Solution misfit: 0.009577915614926758 6134.652199368199  dx: 0.1 1.0\n",
      "---------------------------------\n",
      "Solution misfit: 0.00524160939762612 3489.7327798633314  dx: 0.05 1.0\n",
      "---------------------------------\n",
      "Solution misfit: 0.0027363138020314753 1853.5177577223155  dx: 0.025 1.0\n",
      "---------------------------------\n",
      "Solution misfit: 0.0013968591960122623 954.5476537100484  dx: 0.0125 1.0\n",
      "---------------------------------\n",
      "Solution misfit: 0.0007055448938828217 484.2725548864858  dx: 0.00625 1.0\n",
      "---------------------------------\n",
      "Solution misfit: 0.0003545420000861747 243.88870492481797  dx: 0.003125 1.0\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_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,2))\n",
    "i_ref = 0\n",
    "\n",
    "for N_refine in N_refines:\n",
    "    \n",
    "    # Define mesh\n",
    "    N_mesh = N_refine\n",
    "    mesh = UnitIntervalMesh(N_mesh)\n",
    "\n",
    "    # =========================================================\n",
    "    # DEFINE BOUNDARY ELEMENTS\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",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE MAT. PROPERTIES\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",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE ELEMENT FUNCTION SPACE AND VARIABLES\n",
    "\n",
    "    deg_u = 1\n",
    "    deg_stress = 0\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",
    "    metadata = {\"quadrature_degree\": deg_stress, \"quadrature_scheme\": \"default\"}\n",
    "    dxm = dx(metadata=metadata)\n",
    "\n",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE NEUMANN BCs \n",
    "    \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",
    "    n = FacetNormal(mesh)\n",
    "\n",
    "    q_lim = float(2/sqrt(3)*sig0*mu*2)\n",
    "    \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",
    "    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",
    "\n",
    "    def F_ext(v):\n",
    "        return dot(loading,v)*dss(3)\n",
    "\n",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE FUNCTIONS FOR FUNCTIONAL CONSTRUCTION\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",
    "\n",
    "    def sigma_tang(e):\n",
    "        return 2*mu*beta*e\n",
    "\n",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE LHS MATRIX AND RHS RESIDUAL\n",
    "    a_Newton = rho*inner(u_,v)/(betaN*dt*dt)*dxm + inner(2*eps(u_), sigma_tang(eps(v)) )*dxm\n",
    "\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",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE FUNCTIONS FOR FUNCTION PROJECTION\n",
    "    \n",
    "    def proj_sig(old_u,old_sig,old_Du):\n",
    "        eps_old = eps(old_u)\n",
    "        beta = 1+3*(2*eps_old[0]/sig0)*(2*eps_old[0]/sig0)\n",
    "        return 2*mu*beta*eps(old_Du)+old_sig, beta\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",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # DEFINE FUNCTIONS FOR U, V, A UPDATE\n",
    "    \n",
    "    # Update formula for acceleration\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",
    "    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",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # FILE OUTPUT\n",
    "    file_results = XDMFFile(\"nonlinearElas1DVeri__results.xdmf\")\n",
    "    file_results.parameters[\"flush_output\"] = True\n",
    "    file_results.parameters[\"functions_share_mesh\"] = True\n",
    "    # =========================================================\n",
    "    \n",
    "    # =========================================================\n",
    "    # ITERATIVE SOLVER\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[0:]):\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",
    "        # -- 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",
    "        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=t, 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(assemble(error))\n",
    "    errors_result[i_ref,0] = Errors_theo\n",
    "    error = (sig_theo - sig)**2*dxm\n",
    "    Errors_theo1 = sqrt(assemble(error))\n",
    "    errors_result[i_ref,1] = Errors_theo1\n",
    "    dts[i_ref] = 1/Nincr\n",
    "    i_ref += 1\n",
    "\n",
    "    print(\"Solution misfit:\", Errors_theo,Errors_theo1,\" dx:\",1/N_mesh,t)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37e95632",
   "metadata": {},
   "outputs": [],
   "source": [
    "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
}