2D-Wave-DamageDG-Fsource.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b0b318bc",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<frozen importlib._bootstrap>:228: RuntimeWarning: scipy._lib.messagestream.MessageStream size changed, may indicate binary incompatibility. Expected 56 from C header, got 64 from PyObject\n"
     ]
    }
   ],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "from dolfin import *\n",
    "parameters[\"form_compiler\"][\"representation\"] = 'quadrature'\n",
    "set_log_level(30)\n",
    "import warnings\n",
    "from ffc.quadrature.deprecation import QuadratureRepresentationDeprecationWarning\n",
    "warnings.simplefilter(\"once\", QuadratureRepresentationDeprecationWarning)\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from scipy.optimize import curve_fit\n",
    "import matplotlib.pyplot as plt\n",
    "import obspy\n",
    "from obspy.imaging import spectrogram as spect\n",
    "from obspy.signal.invsim import cosine_taper\n",
    "\n",
    "import math\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2c7d8f0",
   "metadata": {},
   "source": [
    "# With slow dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a045be7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plotting a RectangleMesh\n"
     ]
    },
    {
     "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",
      "Time step: 0\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 1\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 2\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 3\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 4\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 5\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 6\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 7\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 8\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 9\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 10\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 11\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 12\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 13\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 14\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 15\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 16\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 17\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 18\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 19\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 20\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 21\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 22\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 23\n",
      "Strain: 0.0\n",
      "===============\n",
      "Time step: 24\n",
      "Strain: -1.12750060078e-28\n",
      "===============\n",
      "Time step: 25\n",
      "Strain: 4.28808628186e-28\n",
      "===============\n",
      "Time step: 26\n",
      "Strain: 1.11840898585e-26\n",
      "===============\n",
      "Time step: 27\n",
      "Strain: 4.02893053607e-26\n",
      "===============\n",
      "Time step: 28\n",
      "Strain: -3.66226612197e-26\n",
      "===============\n",
      "Time step: 29\n",
      "Strain: -6.36972801422e-25\n",
      "===============\n",
      "Time step: 30\n",
      "Strain: -1.71245686895e-24\n",
      "===============\n",
      "Time step: 31\n",
      "Strain: -1.30772972545e-25\n",
      "===============\n",
      "Time step: 32\n",
      "Strain: 1.0717629146e-23\n",
      "===============\n",
      "Time step: 33\n",
      "Strain: 2.67301601824e-23\n",
      "===============\n",
      "Time step: 34\n",
      "Strain: 1.56160876015e-24\n",
      "===============\n",
      "Time step: 35\n",
      "Strain: -1.52284609158e-22\n",
      "===============\n",
      "Time step: 36\n",
      "Strain: -4.2315133429e-22\n",
      "===============\n",
      "Time step: 37\n",
      "Strain: -4.06390523879e-22\n",
      "===============\n",
      "Time step: 38\n",
      "Strain: 7.87598896654e-22\n",
      "===============\n",
      "Time step: 39\n",
      "Strain: 3.67406443321e-21\n",
      "===============\n",
      "Time step: 40\n",
      "Strain: 6.23935484289e-21\n",
      "===============\n",
      "Time step: 41\n",
      "Strain: 2.10716192591e-21\n",
      "===============\n",
      "Time step: 42\n",
      "Strain: -1.67061330034e-20\n",
      "===============\n",
      "Time step: 43\n",
      "Strain: -4.83821991832e-20\n",
      "===============\n",
      "Time step: 44\n",
      "Strain: -6.46351441074e-20\n",
      "===============\n",
      "Time step: 45\n",
      "Strain: -6.48879059904e-21\n",
      "===============\n",
      "Time step: 46\n",
      "Strain: 1.7989694485e-19\n",
      "===============\n",
      "Time step: 47\n",
      "Strain: 4.54647923889e-19\n",
      "===============\n",
      "Time step: 48\n",
      "Strain: 5.83693232147e-19\n",
      "===============\n",
      "Time step: 49\n",
      "Strain: 1.4117598175e-19\n",
      "===============\n",
      "Time step: 50\n",
      "Strain: -1.24132672455e-18\n",
      "===============\n",
      "Time step: 51\n",
      "Strain: -3.3393527589e-18\n",
      "===============\n",
      "Time step: 52\n",
      "Strain: -4.7514946403e-18\n",
      "===============\n",
      "Time step: 53\n",
      "Strain: -2.84995091341e-18\n",
      "===============\n",
      "Time step: 54\n",
      "Strain: 5.02051382563e-18\n",
      "===============\n",
      "Time step: 55\n",
      "Strain: 1.88239374953e-17\n",
      "===============\n",
      "Time step: 56\n",
      "Strain: 3.2502460004e-17\n",
      "===============\n",
      "Time step: 57\n",
      "Strain: 3.24981112115e-17\n",
      "===============\n",
      "Time step: 58\n",
      "Strain: 1.50158022707e-18\n",
      "===============\n",
      "Time step: 59\n",
      "Strain: -7.03270591441e-17\n",
      "===============\n",
      "Time step: 60\n",
      "Strain: -1.68510755821e-16\n",
      "===============\n",
      "Time step: 61\n",
      "Strain: -2.40010434865e-16\n",
      "===============\n",
      "Time step: 62\n",
      "Strain: -1.95494237001e-16\n",
      "===============\n",
      "Time step: 63\n",
      "Strain: 5.78728119229e-17\n",
      "===============\n",
      "Time step: 64\n",
      "Strain: 5.49266936556e-16\n",
      "===============\n",
      "Time step: 65\n",
      "Strain: 1.16022457008e-15\n",
      "===============\n",
      "Time step: 66\n",
      "Strain: 1.56809752661e-15\n",
      "===============\n",
      "Time step: 67\n",
      "Strain: 1.27937962149e-15\n",
      "===============\n",
      "Time step: 68\n",
      "Strain: -1.93171416505e-16\n",
      "===============\n",
      "Time step: 69\n",
      "Strain: -3.00496822971e-15\n",
      "===============\n",
      "Time step: 70\n",
      "Strain: -6.58975805148e-15\n",
      "===============\n",
      "Time step: 71\n",
      "Strain: -9.392571326e-15\n",
      "===============\n",
      "Time step: 72\n",
      "Strain: -8.97926503054e-15\n",
      "===============\n",
      "Time step: 73\n",
      "Strain: -2.74528171927e-15\n",
      "===============\n",
      "Time step: 74\n",
      "Strain: 1.072962322e-14\n",
      "===============\n",
      "Time step: 75\n",
      "Strain: 2.99526137549e-14\n",
      "===============\n",
      "Time step: 76\n",
      "Strain: 4.90712706098e-14\n",
      "===============\n",
      "Time step: 77\n",
      "Strain: 5.76603262002e-14\n",
      "===============\n",
      "Time step: 78\n",
      "Strain: 4.26974449964e-14\n",
      "===============\n",
      "Time step: 79\n",
      "Strain: -6.77410328402e-15\n",
      "===============\n",
      "Time step: 80\n",
      "Strain: -9.27697242631e-14\n",
      "===============\n",
      "Time step: 81\n",
      "Strain: -2.01077270139e-13\n",
      "===============\n",
      "Time step: 82\n",
      "Strain: -2.96649551233e-13\n",
      "===============\n",
      "Time step: 83\n",
      "Strain: -3.25422232665e-13\n",
      "===============\n",
      "Time step: 84\n",
      "Strain: -2.25713953204e-13\n",
      "===============\n",
      "Time step: 85\n",
      "Strain: 4.98789526072e-14\n",
      "===============\n",
      "Time step: 86\n",
      "Strain: 5.05014791693e-13\n",
      "===============\n",
      "Time step: 87\n",
      "Strain: 1.07018446699e-12\n",
      "===============\n",
      "Time step: 88\n",
      "Strain: 1.58510993984e-12\n",
      "===============\n",
      "Time step: 89\n",
      "Strain: 1.80650633164e-12\n",
      "===============\n",
      "Time step: 90\n",
      "Strain: 1.45232891622e-12\n",
      "===============\n",
      "Time step: 91\n",
      "Strain: 2.85925625019e-13\n",
      "===============\n",
      "Time step: 92\n",
      "Strain: -1.76916068277e-12\n",
      "===============\n",
      "Time step: 93\n",
      "Strain: -4.50849282794e-12\n",
      "===============\n",
      "Time step: 94\n",
      "Strain: -7.35919314078e-12\n",
      "===============\n",
      "Time step: 95\n",
      "Strain: -9.37240994792e-12\n",
      "===============\n",
      "Time step: 96\n",
      "Strain: -9.33863467233e-12\n",
      "===============\n",
      "Time step: 97\n",
      "Strain: -6.04912586345e-12\n",
      "===============\n",
      "Time step: 98\n",
      "Strain: 1.30881447102e-12\n",
      "===============\n",
      "Time step: 99\n",
      "Strain: 1.26820826905e-11\n",
      "===============\n",
      "Time step: 100\n",
      "Strain: 2.67162101183e-11\n",
      "===============\n",
      "Time step: 101\n",
      "Strain: 4.05006043167e-11\n",
      "===============\n",
      "Time step: 102\n",
      "Strain: 4.96524704568e-11\n",
      "===============\n",
      "Time step: 103\n",
      "Strain: 4.88685196607e-11\n",
      "===============\n",
      "Time step: 104\n",
      "Strain: 3.29976047368e-11\n",
      "===============\n",
      "Time step: 105\n",
      "Strain: -1.43277326369e-12\n",
      "===============\n",
      "Time step: 106\n",
      "Strain: -5.44520047355e-11\n",
      "===============\n",
      "Time step: 107\n",
      "Strain: -1.21009249581e-10\n",
      "===============\n",
      "Time step: 108\n",
      "Strain: -1.89941267585e-10\n",
      "===============\n",
      "Time step: 109\n",
      "Strain: -2.44057455572e-10\n",
      "===============\n",
      "Time step: 110\n",
      "Strain: -2.61776467738e-10\n",
      "===============\n",
      "Time step: 111\n",
      "Strain: -2.2053924416e-10\n",
      "===============\n",
      "Time step: 112\n",
      "Strain: -1.01887488985e-10\n",
      "===============\n",
      "Time step: 113\n",
      "Strain: 1.0232965078e-10\n",
      "===============\n",
      "Time step: 114\n",
      "Strain: 3.83608300383e-10\n",
      "===============\n",
      "Time step: 115\n",
      "Strain: 7.11672891236e-10\n",
      "===============\n",
      "Time step: 116\n",
      "Strain: 1.03217391253e-09\n",
      "===============\n",
      "Time step: 117\n",
      "Strain: 1.26873217888e-09\n",
      "===============\n",
      "Time step: 118\n",
      "Strain: 1.33047452264e-09\n",
      "===============\n",
      "Time step: 119\n",
      "Strain: 1.12546513852e-09\n",
      "===============\n",
      "Time step: 120\n",
      "Strain: 5.79387358612e-10\n",
      "===============\n",
      "Time step: 121\n",
      "Strain: -3.42411159641e-10\n",
      "===============\n",
      "Time step: 122\n",
      "Strain: -1.61264831895e-09\n",
      "===============\n",
      "Time step: 123\n",
      "Strain: -3.12378177975e-09\n",
      "===============\n",
      "Time step: 124\n",
      "Strain: -4.67835970265e-09\n",
      "===============\n",
      "Time step: 125\n",
      "Strain: -5.99265214913e-09\n",
      "===============\n",
      "Time step: 126\n",
      "Strain: -6.71719252673e-09\n",
      "===============\n",
      "Time step: 127\n",
      "Strain: -6.47599055447e-09\n",
      "===============\n",
      "Time step: 128\n",
      "Strain: -4.92350991905e-09\n",
      "===============\n",
      "Time step: 129\n",
      "Strain: -1.81528025809e-09\n",
      "===============\n",
      "Time step: 130\n",
      "Strain: 2.91537747759e-09\n",
      "===============\n",
      "Time step: 131\n",
      "Strain: 9.08509341237e-09\n",
      "===============\n",
      "Time step: 132\n",
      "Strain: 1.62054195938e-08\n",
      "===============\n",
      "Time step: 133\n",
      "Strain: 2.34585600063e-08\n",
      "===============\n",
      "Time step: 134\n",
      "Strain: 2.97163258098e-08\n",
      "===============\n",
      "Time step: 135\n",
      "Strain: 3.36116418282e-08\n",
      "===============\n",
      "Time step: 136\n",
      "Strain: 3.36667245024e-08\n",
      "===============\n",
      "Time step: 137\n",
      "Strain: 2.84750803581e-08\n",
      "===============\n",
      "Time step: 138\n",
      "Strain: 1.69263346152e-08\n",
      "===============\n",
      "Time step: 139\n",
      "Strain: -1.5455479374e-09\n",
      "===============\n",
      "Time step: 140\n",
      "Strain: -2.67177933488e-08\n",
      "===============\n",
      "Time step: 141\n",
      "Strain: -5.73710756666e-08\n",
      "===============\n",
      "Time step: 142\n",
      "Strain: -9.11543526668e-08\n",
      "===============\n",
      "Time step: 143\n",
      "Strain: -1.24557163013e-07\n",
      "===============\n",
      "Time step: 144\n",
      "Strain: -1.53018891101e-07\n",
      "===============\n",
      "Time step: 145\n",
      "Strain: -1.71194597822e-07\n",
      "===============\n",
      "Time step: 146\n",
      "Strain: -1.73382962003e-07\n",
      "===============\n",
      "Time step: 147\n",
      "Strain: -1.54104593434e-07\n",
      "===============\n",
      "Time step: 148\n",
      "Strain: -1.0879986845e-07\n",
      "===============\n",
      "Time step: 149\n",
      "Strain: -3.45962974097e-08\n",
      "===============\n",
      "Time step: 150\n",
      "Strain: 6.89217073347e-08\n",
      "===============\n",
      "Time step: 151\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Strain: 1.99020856144e-07\n",
      "===============\n",
      "Time step: 152\n",
      "Strain: 3.49292065426e-07\n",
      "===============\n",
      "Time step: 153\n",
      "Strain: 5.09375422103e-07\n",
      "===============\n",
      "Time step: 154\n",
      "Strain: 6.65014938545e-07\n",
      "===============\n",
      "Time step: 155\n",
      "Strain: 7.98503391412e-07\n",
      "===============\n",
      "Time step: 156\n",
      "Strain: 8.89551790707e-07\n",
      "===============\n",
      "Time step: 157\n",
      "Strain: 9.16585532075e-07\n",
      "===============\n",
      "Time step: 158\n",
      "Strain: 8.58432174842e-07\n",
      "===============\n",
      "Time step: 159\n",
      "Strain: 6.96326787346e-07\n",
      "===============\n",
      "Time step: 160\n",
      "Strain: 4.16123024773e-07\n",
      "===============\n",
      "Time step: 161\n",
      "Strain: 1.05648015723e-08\n",
      "===============\n",
      "Time step: 162\n",
      "Strain: -5.1855230458e-07\n",
      "===============\n",
      "Time step: 163\n",
      "Strain: -1.1585162376e-06\n",
      "===============\n",
      "Time step: 164\n",
      "Strain: -1.88431971372e-06\n",
      "===============\n",
      "Time step: 165\n",
      "Strain: -2.65794668341e-06\n",
      "===============\n",
      "Time step: 166\n",
      "Strain: -3.42847910318e-06\n",
      "===============\n",
      "Time step: 167\n",
      "Strain: -4.13314212592e-06\n",
      "===============\n",
      "Time step: 168\n",
      "Strain: -4.69935760803e-06\n",
      "===============\n",
      "Time step: 169\n",
      "Strain: -5.04781895207e-06\n",
      "===============\n",
      "Time step: 170\n",
      "Strain: -5.09653794871e-06\n",
      "===============\n",
      "Time step: 171\n",
      "Strain: -4.76575005929e-06\n",
      "===============\n",
      "Time step: 172\n",
      "Strain: -3.98350225981e-06\n",
      "===============\n",
      "Time step: 173\n",
      "Strain: -2.69169083095e-06\n",
      "===============\n",
      "Time step: 174\n",
      "Strain: -8.52268691215e-07\n",
      "===============\n",
      "Time step: 175\n",
      "Strain: 1.54669411174e-06\n",
      "===============\n",
      "Time step: 176\n",
      "Strain: 4.48543510931e-06\n",
      "===============\n",
      "Time step: 177\n",
      "Strain: 7.90775417976e-06\n",
      "===============\n",
      "Time step: 178\n",
      "Strain: 1.17174335846e-05\n",
      "===============\n",
      "Time step: 179\n",
      "Strain: 1.57761327689e-05\n",
      "===============\n",
      "Time step: 180\n",
      "Strain: 1.9903026459e-05\n",
      "===============\n",
      "Time step: 181\n",
      "Strain: 2.38763998443e-05\n",
      "===============\n",
      "Time step: 182\n",
      "Strain: 2.74373481313e-05\n",
      "===============\n",
      "Time step: 183\n",
      "Strain: 3.02956519808e-05\n",
      "===============\n",
      "Time step: 184\n",
      "Strain: 3.21378181485e-05\n",
      "===============\n",
      "Time step: 185\n",
      "Strain: 3.26371891296e-05\n",
      "===============\n",
      "Time step: 186\n",
      "Strain: 3.14659400249e-05\n",
      "===============\n",
      "Time step: 187\n",
      "Strain: 2.83086985147e-05\n",
      "===============\n",
      "Time step: 188\n",
      "Strain: 2.2877447984e-05\n",
      "===============\n",
      "Time step: 189\n",
      "Strain: 1.49273075902e-05\n",
      "===============\n",
      "Time step: 190\n",
      "Strain: 4.27272918925e-06\n",
      "===============\n",
      "Time step: 191\n",
      "Strain: -9.19638807509e-06\n",
      "===============\n",
      "Time step: 192\n",
      "Strain: -2.54991871937e-05\n",
      "===============\n",
      "Time step: 193\n",
      "Strain: -4.45494730454e-05\n",
      "===============\n",
      "Time step: 194\n",
      "Strain: -6.61429751032e-05\n",
      "===============\n",
      "Time step: 195\n",
      "Strain: -8.99467339817e-05\n",
      "===============\n",
      "Time step: 196\n",
      "Strain: -0.000115491044041\n",
      "===============\n",
      "Time step: 197\n",
      "Strain: -0.000142164320518\n",
      "===============\n",
      "Time step: 198\n",
      "Strain: -0.000169211183406\n",
      "===============\n",
      "Time step: 199\n",
      "Strain: -0.000195733964673\n",
      "===============\n",
      "Time step: 200\n",
      "Strain: -0.000220697753906\n",
      "===============\n",
      "Time step: 201\n",
      "Strain: -0.000242939003758\n",
      "===============\n",
      "Time step: 202\n",
      "Strain: -0.000261177624336\n",
      "===============\n",
      "Time step: 203\n",
      "Strain: -0.000274032408229\n",
      "===============\n",
      "Time step: 204\n",
      "Strain: -0.000280039548181\n",
      "===============\n",
      "Time step: 205\n",
      "Strain: -0.000277673939694\n",
      "===============\n",
      "Time step: 206\n",
      "Strain: -0.000265372902698\n",
      "===============\n",
      "Time step: 207\n",
      "Strain: -0.000241561910554\n",
      "===============\n",
      "Time step: 208\n",
      "Strain: -0.000204681881037\n",
      "===============\n",
      "Time step: 209\n",
      "Strain: -0.000153217561971\n",
      "===============\n",
      "Time step: 210\n",
      "Strain: -8.57265324543e-05\n",
      "===============\n",
      "Time step: 211\n",
      "Strain: -8.68337685391e-07\n",
      "===============\n",
      "Time step: 212\n",
      "Strain: 0.000102566720655\n",
      "===============\n",
      "Time step: 213\n",
      "Strain: 0.000225629607368\n",
      "===============\n",
      "Time step: 214\n",
      "Strain: 0.000369185846154\n",
      "===============\n",
      "Time step: 215\n",
      "Strain: 0.000533890522295\n",
      "===============\n",
      "Time step: 216\n",
      "Strain: 0.000720165419365\n",
      "===============\n",
      "Time step: 217\n",
      "Strain: 0.00092817869497\n",
      "===============\n",
      "Time step: 218\n",
      "Strain: 0.00115782747982\n",
      "===============\n",
      "Time step: 219\n",
      "Strain: 0.00140872376096\n",
      "===============\n",
      "Time step: 220\n",
      "Strain: 0.00168018388267\n",
      "===============\n",
      "Time step: 221\n",
      "Strain: 0.00197122196597\n",
      "===============\n",
      "Time step: 222\n",
      "Strain: 0.00228054750887\n",
      "===============\n",
      "Time step: 223\n",
      "Strain: 0.00260656738316\n",
      "===============\n",
      "Time step: 224\n",
      "Strain: 0.0029473923889\n",
      "===============\n",
      "Time step: 225\n",
      "Strain: 0.00330084846498\n",
      "===============\n",
      "Time step: 226\n",
      "Strain: 0.0036644925827\n",
      "===============\n",
      "Time step: 227\n",
      "Strain: 0.00403563327136\n",
      "===============\n",
      "Time step: 228\n",
      "Strain: 0.00441135564048\n",
      "===============\n",
      "Time step: 229\n",
      "Strain: 0.00478855067558\n",
      "===============\n",
      "Time step: 230\n",
      "Strain: 0.00516394849571\n",
      "===============\n",
      "Time step: 231\n",
      "Strain: 0.00553415517297\n",
      "===============\n",
      "Time step: 232\n",
      "Strain: 0.00589569263156\n",
      "===============\n",
      "Time step: 233\n",
      "Strain: 0.00624504106734\n",
      "===============\n",
      "Time step: 234\n",
      "Strain: 0.00657868326314\n",
      "===============\n",
      "Time step: 235\n",
      "Strain: 0.00689315012026\n",
      "===============\n",
      "Time step: 236\n",
      "Strain: 0.00718506668637\n",
      "===============\n",
      "Time step: 237\n",
      "Strain: 0.00745119793448\n",
      "===============\n",
      "Time step: 238\n",
      "Strain: 0.00768849353809\n",
      "===============\n",
      "Time step: 239\n",
      "Strain: 0.00789413089383\n",
      "===============\n",
      "Time step: 240\n",
      "Strain: 0.00806555566504\n",
      "===============\n",
      "Time step: 241\n",
      "Strain: 0.00820051915598\n",
      "===============\n",
      "Time step: 242\n",
      "Strain: 0.00829711187718\n",
      "===============\n",
      "Time step: 243\n",
      "Strain: 0.00835379272407\n",
      "===============\n",
      "Time step: 244\n",
      "Strain: 0.00836941326428\n",
      "===============\n",
      "Time step: 245\n"
     ]
    }
   ],
   "source": [
    "# Form compiler options\n",
    "parameters[\"form_compiler\"][\"cpp_optimize\"] = True\n",
    "parameters[\"form_compiler\"][\"optimize\"] = True\n",
    "\n",
    "# Define mesh\n",
    "lx = 1\n",
    "ly = 1*lx\n",
    "m_size = 50\n",
    "mesh = RectangleMesh(Point(0., 0.), Point(lx, ly), m_size, int(m_size*ly/lx))\n",
    "print(\"Plotting a RectangleMesh\")\n",
    "plt.figure()\n",
    "plot(mesh, title=\"RectangleMesh\")\n",
    "\n",
    "i_VelBC = True\n",
    "i_heter = False\n",
    "\n",
    "# Material properties\n",
    "mu = Constant(14e9)\n",
    "lamb = Constant(14e9)\n",
    "rho = Constant(2.2e3)\n",
    "gamma = Constant(20)\n",
    "beta_E = Constant(-2000e9) # m\n",
    "# delta_E = delta_E # l\n",
    "if i_heter == True:\n",
    "    gamma = Expression(\"(x[0] < lx0) ? gamma : gamma + (500-gamma)*(x[0]-lx0)/(lx-lx0)\"\\\n",
    "                       , degree=1, lx0=0.5*lx, lx = lx, gamma = 5)\n",
    "    beta_E = Expression(\"(x[0] < lx0) ? gamma : gamma + (0-gamma)*(x[0]-lx0)/(lx-lx0)\"\\\n",
    "                       , degree=1, lx0=0.1*lx, lx = lx, gamma = beta_E)\n",
    "    delta_E = Expression(\"(x[0] < lx0) ? gamma : gamma + (0-gamma)*(x[0]-lx0)/(lx-lx0)\"\\\n",
    "                       , degree=1, lx0=0.1*lx, lx = lx, gamma = delta_E)\n",
    "tau_min = 1e-3\n",
    "tau_max = 10e-3\n",
    "tau = tau_min*gamma # \\tau_1 in the literature\n",
    "\n",
    "delta_t = 0.5e-6\n",
    "t_total = 2e-4\n",
    "\n",
    "t_excitation = 0.25e-4\n",
    "\n",
    "fc = 20e3\n",
    "\n",
    "A_amp = 1e4\n",
    "\n",
    "g0 = 1e0\n",
    "\n",
    "deg_dg = 1 # DG degree\n",
    "\n",
    "def strain_from_q(q_f,i_case):\n",
    "    (e11,e12,e21,e22,v1,v2,g_int) = split(q_f)\n",
    "    if i_case == 0:\n",
    "        strain_tensor = as_tensor([[e11,e12],\\\n",
    "                                   [e21,e22]])\n",
    "    elif i_case == 1:\n",
    "        strain_tensor = as_tensor([[e11('+'),e12('+')],\\\n",
    "                                   [e21('+'),e22('+')]])\n",
    "    elif i_case == 2:\n",
    "        strain_tensor = as_tensor([[e11('-'),e12('-')],\\\n",
    "                                   [e21('-'),e22('-')]])\n",
    "    return strain_tensor\n",
    "\n",
    "def mu_fun(eps):\n",
    "    return mu\n",
    "\n",
    "def stress_fun(eps):\n",
    "    sig11 = lamb*(eps[0][0]+eps[1][1]) + (mu)*(eps[0][0]+eps[0][0])\n",
    "    sig22 = lamb*(eps[0][0]+eps[1][1]) + (mu)*(eps[1][1]+eps[1][1])\n",
    "    sig21 = (mu)*(eps[0][1]+eps[1][0])\n",
    "    sig12 = (mu)*(eps[0][1]+eps[1][0])\n",
    "    return as_tensor([[sig11,sig12],\\\n",
    "                      [sig21,sig22]])\n",
    "\n",
    "def energy_fun(eps):\n",
    "    return (lamb+2*mu)*(eps[0][0]+eps[1][1])*(eps[0][0]+eps[1][1]) \\\n",
    "          -2*mu*(eps[0][0]*eps[1][1] - (eps[0][1]+eps[1][0])*(eps[0][1]+eps[1][0])/4)\n",
    "\n",
    "EPS_SMALL = 1e-3\n",
    "def phig_fun(g_inti):\n",
    "    return g_inti/(tau_max*(exp(-g_inti/g0*math.log(tau_max/tau_min)))\\\n",
    "                            *(g_inti/g0*math.log(tau_max/tau_min)* (1-ln(g_inti/g0+EPS_SMALL)) + 1) \\\n",
    "                    + tau_min*(g_inti)/g0*math.log(tau_max/tau_min) )\n",
    "\n",
    "# Berjamin's model\n",
    "# def phig_fun(g_inti):\n",
    "#     return 1/tau_min*g_inti/(1-g_inti*g_inti)\n",
    "\n",
    "def Max(a, b): return (a+b+abs(a-b))/Constant(2)\n",
    "def Min(a, b): return (a+b-abs(a-b))/Constant(2)\n",
    "\n",
    "# Function space\n",
    "Ee11 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ee12 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ee21 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ee22 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ve1 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ve2 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Ve2 = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "Gint = FiniteElement(\"DG\", mesh.ufl_cell(), deg_dg)\n",
    "V_dg = FunctionSpace(mesh, MixedElement([Ee11, Ee12, Ee21, Ee22, Ve1, Ve2, Gint]))\n",
    "V_out = VectorFunctionSpace(mesh,\"DG\", deg_dg) \n",
    "\n",
    "q_vout = Function(V_out)\n",
    "\n",
    "phi, v = TrialFunction(V_dg), TestFunction(V_dg)\n",
    "q_new, q_old = Function(V_dg), Function(V_dg, name='q_field')\n",
    "q_rk = Function(V_dg)\n",
    "\n",
    "q0 = Expression((\"0\",\"0\",\"0\",\"0\",\"0\",\"0\",\"0\"), degree=1, t=0, name='q0') # \n",
    "q_old = project(q0, V_dg)\n",
    "\n",
    "v_exact = Expression(\"(t<tc1 ) ? A*(lx-x[0])/lx*sin(2*pi*(x[0] + c*t)) : 0 \", \\\n",
    "                    degree=1, t=0, tc1=t_excitation, c=fc, A=A_amp, lx = lx) # \n",
    "\n",
    "# v_exact2 = Expression(\"(t<tc1 ) ? A*(lx-x[0])/lx*sin(2*pi*(x[0] + c*t)) : 0 \", \\\n",
    "#                     degree=1, t=0, tc1=t_excitation, c=fc, A=A_amp2, lx = lx) # \n",
    "\n",
    "v_exactr = Expression(\" 0 \", \\\n",
    "                    degree=1, t=0, tc1=t_excitation, c=fc, A=A_amp, lx = lx) # \n",
    "\n",
    "# f_ext = Expression((\"(t<tc1 ) ? -A1*(lx-x[0])/lx*sin(2*pi*(x[0] + c*t)) : 0 \",\\\n",
    "#                     \"(t<tc1 ) ? -A2*(lx-x[0])/lx*sin(2*pi*(x[0] + c*t)) : 0 \",\\\n",
    "#                     \"0\"), \\\n",
    "#                 degree=1, t=0, tc1=t_excitation, c=fc, A1=A_amp2, A2=A_amp1, lx = lx) #\n",
    "\n",
    "DOFIN_EPS = 1e-10 # + lx/m_size/2\n",
    "if deg_dg == 0: DOFIN_EPS = 1e-10 + lx/m_size/2\n",
    "def left(x, on_boundary):\n",
    "    return (abs(x[0]-0.) < DOFIN_EPS and abs(x[1]-ly/2) < 0.03)\n",
    "def right(x, on_boundary):\n",
    "    return abs(x[0]-lx) < DOFIN_EPS\n",
    "\n",
    "boundary_markers = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)\n",
    "\n",
    "class BoundaryNR(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        DOFIN_EPS = 1e-10 # + lx/m_size/2\n",
    "        if deg_dg == 0: DOFIN_EPS = 1e-10 + lx/m_size/2\n",
    "        return on_boundary and near(x[0], lx, DOFIN_EPS)\n",
    "\n",
    "class BoundaryNL(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        DOFIN_EPS = 1e-10 # + lx/m_size/2\n",
    "        if deg_dg == 0: DOFIN_EPS = 1e-10 + lx/m_size/2\n",
    "        return on_boundary and near(x[0], 0, DOFIN_EPS)\n",
    "\n",
    "class BoundaryNT(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        DOFIN_EPS = 1e-10 # + lx/m_size/2\n",
    "        if deg_dg == 0: DOFIN_EPS = 1e-10 + ly/m_size/2\n",
    "        return on_boundary and near(x[1], ly, DOFIN_EPS)\n",
    "    \n",
    "class BoundaryNB(SubDomain):\n",
    "    def inside(self, x, on_boundary):\n",
    "        DOFIN_EPS = 1e-10 # + lx/m_size/2\n",
    "        if deg_dg == 0: DOFIN_EPS = 1e-10 + ly/m_size/2\n",
    "        return on_boundary and near(x[1], 0, DOFIN_EPS)\n",
    "\n",
    "bx0 = BoundaryNR()\n",
    "bx0.mark(boundary_markers, 0)\n",
    "bx0 = BoundaryNL()\n",
    "bx0.mark(boundary_markers, 1)\n",
    "bx0 = BoundaryNT()\n",
    "bx0.mark(boundary_markers, 0)\n",
    "bx0 = BoundaryNB()\n",
    "bx0.mark(boundary_markers, 1)\n",
    "\n",
    "dss = Measure('ds', domain=mesh, subdomain_data=boundary_markers)\n",
    "\n",
    "if i_VelBC == False:\n",
    "    bc =  [DirichletBC(V_dg.sub(0), v_exact, left, method=\"geometric\")]\n",
    "else:\n",
    "    bc =  [DirichletBC(V_dg.sub(4), v_exact, left, method=\"geometric\")]\n",
    "    \n",
    "# Weak form\n",
    "n = FacetNormal(mesh)\n",
    "\n",
    "## First stage\n",
    "(e11_1,e12_1,e21_1,e22_1,v1_1,v2_1,g_int1) = split(q_old)\n",
    "strainp1 = strain_from_q(q_old,1)\n",
    "stressp1 = stress_fun(strainp1)\n",
    "strainm1 = strain_from_q(q_old,2)\n",
    "stressm1 = stress_fun(strainm1)\n",
    "stre_avg1 = ((1-g_int1('+'))*stressp1+(1-g_int1('-'))*stressm1)/2\n",
    "\n",
    "strain1 = strain_from_q(q_old,0)\n",
    "stress1 = stress_fun(strain1)\n",
    "W_energy1 = energy_fun(strain1)\n",
    "\n",
    "spec_rad1 = sqrt((lamb + 2*mu)/rho)\n",
    "\n",
    "F = dot( (phi-q_old)/delta_t, v) * dx \n",
    "\n",
    "flux_innerf1 = as_vector([-v1_1, 0, -v2_1, 0, -(1-g_int1)*stress1[0][0]/rho, -(1-g_int1)*stress1[1][0]/rho, 0])\n",
    "flux_innerg1 = as_vector([0, -v1_1, 0, -v2_1, -(1-g_int1)*stress1[0][1]/rho, -(1-g_int1)*stress1[1][1]/rho, 0])\n",
    "\n",
    "F1 = F - dot(flux_innerf1, grad(v)[:,0]) * dx - dot(flux_innerg1, grad(v)[:,1]) * dx\n",
    "\n",
    "f_ext = Expression((\"(t<tc1 && sqrt((x[1]-ly/2)*(x[1]-ly/2) + (x[0]-lx/2)*(x[0]-lx/2)) < 0.03) ? 0*sin(2*pi*c*t) : 0 \"\\\n",
    "                    ,\"0\",\"0\",\"0\"\\\n",
    "                    ,\"(t<tc1 && sqrt((x[1]-ly/2)*(x[1]-ly/2) + (x[0]-lx/2)*(x[0]-lx/2)) < 0.03) ? A*sin(2*pi*c*t) : 0 \"\\\n",
    "                    ,\"(t<tc1 && sqrt((x[1]-ly/2)*(x[1]-ly/2) + (x[0]-lx/2)*(x[0]-lx/2)) < 0.03) ? 0*sin(2*pi*c*t) : 0 \"\\\n",
    "                    ,\"0\")\\\n",
    "                    ,degree=1, t=0, tc1=t_excitation, c=fc, A=A_amp, lx = lx, ly=ly)\n",
    "\n",
    "F1 += - dot(f_ext, v) * dx\n",
    "\n",
    "r1 = as_vector([0, 0, 0, 0, 0, 0, 1/tau*W_energy1 - phig_fun(g_int1)])\n",
    "F1 += - dot(r1, v) * dx\n",
    "\n",
    "# Rusanov solver\n",
    "flux1 = 1*n('+')[0]\\\n",
    "                *as_vector([-(v1_1('+')+v1_1('-'))/2, 0, -(v2_1('+')+v2_1('-'))/2, 0 \\\n",
    "                            , -stre_avg1[0,0]/rho, -stre_avg1[1,0]/rho, 0]) \\\n",
    "            + 1*n('+')[1]\\\n",
    "                *as_vector([0, -(v1_1('+')+v1_1('-'))/2, 0, -(v2_1('+')+v2_1('-'))/2 \\\n",
    "                            , -stre_avg1[0,1]/rho, -stre_avg1[1,1]/rho, 0]) \\\n",
    "            + 1.0*spec_rad1/2*jump(q_old)\n",
    "\n",
    "F1 += dot(flux1, jump(v)) * dS + dot(flux_innerf1,v)*n[0] * dss(1) + dot(flux_innerg1,v)*n[1] * dss(1)\\\n",
    "                                + dot(flux_innerf1,v)*n[0] * dss(0) + dot(flux_innerg1,v)*n[1] * dss(0) \n",
    "# + dot(flux_inner1,v)*n[0] * dss(0) \n",
    "\n",
    "## Second stage\n",
    "(e11_2,e12_2,e21_2,e22_2,v1_2,v2_2,g_int2) = split(q_rk)\n",
    "strainp2 = strain_from_q(q_rk,1)\n",
    "stressp2 = stress_fun(strainp2)\n",
    "strainm2 = strain_from_q(q_rk,2)\n",
    "stressm2 = stress_fun(strainm2)\n",
    "stre_avg2 = ((1-g_int2('+'))*stressp2+(1-g_int2('-'))*stressm2)/2\n",
    "\n",
    "strain2 = strain_from_q(q_rk,0)\n",
    "stress2 = stress_fun(strain2)\n",
    "W_energy2 = energy_fun(strain2)\n",
    "\n",
    "spec_rad2 = sqrt((lamb + 2*mu)/rho)\n",
    "\n",
    "flux_innerf2 = as_vector([-v1_2, 0, -v2_2, 0, -(1-g_int2)*stress2[0][0]/rho, -(1-g_int2)*stress2[1][0]/rho, 0])\n",
    "flux_innerg2 = as_vector([0, -v1_2, 0, -v2_2, -(1-g_int2)*stress2[0][1]/rho, -(1-g_int2)*stress2[1][1]/rho, 0])\n",
    "\n",
    "# F2 = - dot(dot(A_adv2,q_rk), grad(v)[:,0]) * dx\n",
    "F2 = - dot(flux_innerf2, grad(v)[:,0]) * dx - dot(flux_innerg2, grad(v)[:,1]) * dx\n",
    "\n",
    "# f_ext = Expression((\"0\",\"0\",\"0\",\"0\"\\\n",
    "#                     ,\"(t<tc1 && sqrt((x[1]-ly/2)*(x[1]-ly/2) + (x[0]-lx/2)*(x[0]-lx/2)) < 0.03) ? A*sin(2*pi*c*t) : 0 \"\\\n",
    "#                     ,\"(t<tc1 && sqrt((x[1]-ly/2)*(x[1]-ly/2) + (x[0]-lx/2)*(x[0]-lx/2)) < 0.03) ? A*sin(2*pi*c*t) : 0 \"\\\n",
    "#                     ,\"0\")\\\n",
    "#                     ,degree=1, t=0, tc1=t_excitation, c=fc, A=A_amp, lx = lx, ly=ly)\n",
    "\n",
    "F2 += - dot(f_ext, v) * dx\n",
    "\n",
    "r2 = as_vector([0, 0, 0, 0, 0, 0, 1/tau*W_energy1 - phig_fun(g_int1)])\n",
    "\n",
    "F2 += - dot(r2, v) * dx\n",
    "\n",
    "# Rusanov solver\n",
    "flux2 = 1*n('+')[0]\\\n",
    "                *as_vector([-(v1_2('+')+v1_2('-'))/2, 0, -(v2_2('+')+v2_2('-'))/2, 0 \\\n",
    "                            , -stre_avg2[0,0]/rho, -stre_avg2[1,0]/rho, 0]) \\\n",
    "            + 1*n('+')[1]\\\n",
    "                *as_vector([0, -(v1_2('+')+v1_2('-'))/2, 0, -(v2_2('+')+v2_2('-'))/2 \\\n",
    "                            , -stre_avg2[0,1]/rho, -stre_avg2[1,1]/rho, 0]) \\\n",
    "            + 1.0*spec_rad2/2*jump(q_rk)\n",
    "\n",
    "F2 += dot(flux2, jump(v)) * dS + dot(flux_innerf2,v)*n[0] * dss(1) + dot(flux_innerg2,v)*n[1] * dss(1)\\\n",
    "                               + dot(flux_innerf2,v)*n[0] * dss(0) + dot(flux_innerg2,v)*n[1] * dss(0) \n",
    "F += 0.5*((F1-F) + F2)\n",
    "\n",
    "#===========\n",
    "Sta1 = time.time()\n",
    "T_sol = 0\n",
    "i_print = True\n",
    "i_plot = False\n",
    "i_file = True\n",
    "#===========\n",
    "\n",
    "xdmf_file = XDMFFile(\"2D-waveDG_F.xdmf\")\n",
    "xdmf_file.parameters[\"flush_output\"] = True\n",
    "xdmf_file.parameters[\"functions_share_mesh\"] = True\n",
    "xdmf_file.parameters[\"rewrite_function_mesh\"] = False\n",
    "\n",
    "t = 0\n",
    "time_steps = int(t_total/delta_t)\n",
    "uuu = np.zeros((time_steps,m_size,7))\n",
    "\n",
    "a1,L1 = lhs(F1),rhs(F1)\n",
    "a,L = lhs(F),rhs(F)\n",
    "\n",
    "for n in range(time_steps):\n",
    "    if i_print:\n",
    "        print('===============')\n",
    "        print('Time step:',n)\n",
    "\n",
    "    for i in range(m_size):\n",
    "        uuu[n,i,:] = q_old(i/(m_size-1)*lx,ly/2)[:]\n",
    "\n",
    "\n",
    "    T_incre1 = time.time()\n",
    "    solve(a1==L1, q_new)\n",
    "    T_incre2 = time.time()\n",
    "    T_sol += T_incre2-T_incre1\n",
    "    q_rk.assign(q_new)\n",
    "\n",
    "\n",
    "    T_incre1 = time.time()\n",
    "    solve(a==L, q_new)\n",
    "    T_incre2 = time.time()\n",
    "    T_sol += T_incre2-T_incre1\n",
    "    q_old.assign(q_new)\n",
    "    if i_print:\n",
    "        print('Strain:',q_old(0,ly/2)[4])\n",
    "\n",
    "    t += delta_t\n",
    "    v_exact.t = t\n",
    "    f_ext.t = t\n",
    "    \n",
    "    if i_file: \n",
    "        v1_out = q_old.split()[4]\n",
    "        v1_out.rename(\"Vx\", \"tmp\")\n",
    "        v2_out = q_old.split()[5]\n",
    "        v2_out.rename(\"Vy\", \"tmp\")\n",
    "        xdmf_file.write(v1_out, n)\n",
    "        xdmf_file.write(v2_out, n)\n",
    "        q_vout=project(as_vector([v1_out,v2_out]))\n",
    "        q_vout.rename('V','tmp')\n",
    "        xdmf_file.write(q_vout, n)\n",
    "        g_out = q_old.split()[6]\n",
    "        g_out.rename(\"Damage\", \"tmp\")\n",
    "        xdmf_file.write(g_out, n)\n",
    "\n",
    "T_total = time.time()-Sta1\n",
    "if i_plot:\n",
    "    plt.subplots()\n",
    "    plot(q_old[2])\n",
    "    plt.show()\n",
    "\n",
    "print(\"Total time:\",T_total)\n",
    "print(\"Solving time:\",T_sol)\n",
    "# return uuu"
   ]
  }
 ],
 "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
}