examples/waveguide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helmholz equation on a parameter dependent domain\n",
    "\n",
    "Solve the scalar Helmholz on a domain (depending on 3 parameters) resembeling a rectangular waveguide. \n",
    "\n",
    "Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "import torch as tn\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors\n",
    "import datetime\n",
    "try:\n",
    "    import torchtt as tntt\n",
    "except:\n",
    "    print('Installing torchTT...')\n",
    "    %pip install git+https://github.com/ion-g-ion/torchtt\n",
    "    import torchtt as tntt\n",
    "try:\n",
    "    import tt_iga\n",
    "except:\n",
    "    print('Installing TT-IGA...')\n",
    "    %pip install git+https://github.com/ion-g-ion/code-paper-tt-iga\n",
    "    import tt_iga\n",
    "\n",
    "tn.set_default_dtype(tn.float64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Curves and surfaces parametrizations for the shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "xc = lambda u,v: u*tn.sqrt(1-v**2/2)\n",
    "yc = lambda u,v: v*tn.sqrt(1-u**2/2)\n",
    "\n",
    "line = lambda t,a,b: t*(b-a)+a\n",
    "\n",
    "def plane_spanner(P1,P2,P3,t1,t2):\n",
    "    x = (P1[:,0]-P2[:,0])*t1+(P3[:,0]-P2[:,0])*t2+P2[:,0]\n",
    "    y = (P1[:,1]-P2[:,1])*t1+(P3[:,1]-P2[:,1])*t2+P2[:,1]\n",
    "    z = (P1[:,2]-P2[:,2])*t1+(P3[:,2]-P2[:,2])*t2+P2[:,2]\n",
    "    return x,y,z\n",
    "\n",
    "def curve2(t,a,b,c,d):\n",
    "    w2 = 1.5\n",
    "    ry = 1.5\n",
    "    h = 0.5\n",
    "    rz = 1.5+c\n",
    "    x = -1*(t<0.25)+tn.logical_and(t>=0.25,t<0.5)*line((t-0.25)/0.25,-1,-w2)+tn.logical_and(t>=0.5,t<=1)*(-w2)\n",
    "    y = 0*(t<=0.75) + (t>0.75) * ( -ry*tn.cos((t-0.75)/0.25*np.pi/2) + ry )\n",
    "    z = line(t/0.25,-3,-2+a)*(t<=0.25) + tn.logical_and(t>0.25,t<=0.5)*line((t-0.25)/0.25,-2+a,-1+b) + tn.logical_and(t>0.5,t<=0.75)*line((t-0.5)/0.25,-1+b,-0) + (t>0.75) * ( rz*tn.sin((t-0.75)/0.25*np.pi/2) )\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))\n",
    "\n",
    "def curve1(t,a,b,c,d):\n",
    "    w2 = 1.5\n",
    "    ry = 1.5\n",
    "    h = 0.5\n",
    "    rz = 1.5+c\n",
    "    x = 1*(t<0.25)+tn.logical_and(t>=0.25,t<0.5)*line((t-0.25)/0.25,1,w2)+tn.logical_and(t>=0.5,t<=1)*(w2)\n",
    "    y = 0*(t<=0.75) + (t>0.75) * ( -ry*tn.cos((t-0.75)/0.25*np.pi/2) + ry )\n",
    "    z = line(t/0.25,-3,-2+a)*(t<=0.25) + tn.logical_and(t>0.25,t<=0.5)*line((t-0.25)/0.25,-2+a,-1+b) + tn.logical_and(t>0.5,t<=0.75)*line((t-0.5)/0.25,-1+b,-0) + (t>0.75) * ( rz*tn.sin((t-0.75)/0.25*np.pi/2) )\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))\n",
    "  \n",
    "def curve3(t,a,b,c,d):\n",
    "    w2 = 1.5\n",
    "    ry = 1.5\n",
    "    h = 0.5\n",
    "    rz = 1.5+c\n",
    "    x = -1*(t<0.25)+tn.logical_and(t>=0.25,t<0.5)*line((t-0.25)/0.25,-1,-w2)+tn.logical_and(t>=0.5,t<=1)*(-w2)\n",
    "    y = h*(t<=0.75) + (t>0.75) * ( -(ry-h)*tn.cos((t-0.75)/0.25*np.pi/2) + ry )\n",
    "    z = line(t/0.25,-3,-2+a)*(t<=0.25) + tn.logical_and(t>0.25,t<=0.5)*line((t-0.25)/0.25,-2+a,-1+b) + tn.logical_and(t>0.5,t<=0.75)*line((t-0.5)/0.25,-1+b,-0) + (t>0.75) * ( (rz-h)*tn.sin((t-0.75)/0.25*np.pi/2) )\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))\n",
    "      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the B-spline bases and the parameter space bases. \n",
    "For the parameter space, Gauss-Legendre nodes are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "deg = 2\n",
    "Ns = np.array([60,60,120])-deg+1\n",
    "Ns = np.array([40,40,82])-deg+1\n",
    "baza1 = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Ns[0]),deg)\n",
    "baza2 = tt_iga.bspline.BSplineBasis(np.linspace(0,1,Ns[1]),deg)\n",
    "baza3 = tt_iga.bspline.BSplineBasis(np.concatenate((np.linspace(0,0.25,Ns[2]//4),np.linspace(0.25,0.5,Ns[2]//4),np.linspace(0.5,0.75,Ns[2]//4),np.linspace(0.75,1,Ns[2]//4-1))),deg)\n",
    "\n",
    "Basis = [baza1,baza2,baza3]\n",
    "N = [baza1.N,baza2.N,baza3.N]\n",
    "\n",
    "nl = 8\n",
    "Basis_param = [tt_iga.lagrange.LagrangeLeg(nl,[-0.2,0.2])]*2+[tt_iga.lagrange.LagrangeLeg(nl,[-0.3,0.3])]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the parametrization and create the geometry object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "scale_mult = 1\n",
    "xparam = lambda t : plane_spanner(curve1(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve2(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve3(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),t[:,0],t[:,1])[0]\n",
    "yparam = lambda t : plane_spanner(curve1(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve2(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve3(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),t[:,0],t[:,1])[1]\n",
    "zparam = lambda t : plane_spanner(curve1(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve2(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),curve3(t[:,2],scale_mult*t[:,3],scale_mult*t[:,4],scale_mult*t[:,5],0),t[:,0],t[:,1])[2]\n",
    "\n",
    "# interpolate the geometry parametrization\n",
    "geom = tt_iga.Geometry(Basis+Basis_param)\n",
    "geom.interpolate([xparam, yparam, zparam])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Firstly, plot the  geometry for the parameter $\\theta_1 = \\theta_2 = \\theta_3 = 0$.\n",
    "\n",
    "Secondly, plot the geometries corrrecponding to the edges of the parameter space $[-0.2,0.2] \\times [-0.2,0.2] \\times [-0.3,0.3]$, i.e. $(\\theta_1,\\theta_2,\\theta_3) \\in \\{-0.2,0.2\\} \\times \\{-0.2,0.2\\} \\times \\{-0.3,0.3\\}$.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "\n",
    "fig = geom.plot_domain([tn.tensor([0.0]),tn.tensor([-0.0]),tn.tensor([0.0]),tn.tensor([0.1])],[(0,1),(0,1),(0.0,1)],surface_color='blue', wireframe = False,alpha=0.1,n=64)\n",
    "fig.gca().set_xlabel(r'$x_1$', fontsize=14)\n",
    "fig.gca().set_ylabel(r'$x_2$', fontsize=14)\n",
    "fig.gca().set_zlabel(r'$x_3$', fontsize=14)\n",
    "fig.gca().view_init(15, -60)\n",
    "fig.gca().zaxis.set_rotate_label(False)\n",
    "fig.gca().set_xticks([-1.5, 0, 1.5])\n",
    "fig.gca().set_yticks([0,1])\n",
    "fig.gca().set_zticks([-3,-1.5,0,1.5])\n",
    "fig.gca().tick_params(axis='both', labelsize=14)\n",
    "fig.gca().set_box_aspect(aspect = (3,1.5,4.5))\n",
    "\n",
    "fig = geom.plot_domain([tn.tensor([0.0]), tn.tensor([-0.0]), tn.tensor([0.0])],[(0,1),(0,1),(0.0,1)],surface_color=None, wireframe = False, alpha = 0.1, n = 64, frame_color = 'k')\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[0]]),tn.tensor([Basis_param[1].interval[0]]), tn.tensor([Basis_param[2].interval[0]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[0]]),tn.tensor([Basis_param[1].interval[0]]), tn.tensor([Basis_param[2].interval[1]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[0]]),tn.tensor([Basis_param[1].interval[1]]), tn.tensor([Basis_param[2].interval[0]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[0]]),tn.tensor([Basis_param[1].interval[1]]), tn.tensor([Basis_param[2].interval[1]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[1]]),tn.tensor([Basis_param[1].interval[0]]), tn.tensor([Basis_param[2].interval[0]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[1]]),tn.tensor([Basis_param[1].interval[0]]), tn.tensor([Basis_param[2].interval[1]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[1]]),tn.tensor([Basis_param[1].interval[1]]), tn.tensor([Basis_param[2].interval[0]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "geom.plot_domain([tn.tensor([Basis_param[0].interval[1]]),tn.tensor([Basis_param[1].interval[1]]), tn.tensor([Basis_param[2].interval[1]]) ],[(0,1),(0,1),(0.0,1)],fig = fig,surface_color='blue', wireframe = False,alpha=0.01,n=64)\n",
    "fig.gca().set_xlabel(r'$x_1$', fontsize=14)\n",
    "fig.gca().set_ylabel(r'$x_2$', fontsize=14)\n",
    "fig.gca().set_zlabel(r'$x_3$', fontsize=14)\n",
    "fig.gca().view_init(15, -60)\n",
    "fig.gca().zaxis.set_rotate_label(False)\n",
    "fig.gca().set_xticks([-1.5, 0, 1.5])\n",
    "fig.gca().set_yticks([0,1])\n",
    "fig.gca().set_zticks([-3,-1.5,0,1.5])\n",
    "fig.gca().tick_params(axis='both', labelsize=14)\n",
    "fig.gca().set_box_aspect(aspect = (3,1.5,4.5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the stiffness and mass matrix for the entire parameter grid in the TT format."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "tme = datetime.datetime.now() \n",
    "Mass_tt = geom.mass_interp(eps=1e-11)\n",
    "tme = datetime.datetime.now() -tme\n",
    "print('Time mass matrix ',tme.total_seconds())\n",
    "\n",
    "tme = datetime.datetime.now() \n",
    "Stt = geom.stiffness_interp( func=None, func_reference = None, qtt = False, verb=False)\n",
    "tme = datetime.datetime.now() -tme\n",
    "print('Time stiffness matrix ',tme.total_seconds())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get the boundary/interior projection operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "Pin_tt,Pbd_tt = tt_iga.projectors.get_projectors(N,[[0,0],[0,0],[0,0]])\n",
    "# Pbd_tt = (1/N[0]) * Pbd_tt\n",
    "\n",
    "Pin_tt = Pin_tt ** tntt.eye([nl]*3)\n",
    "Pbd_tt = Pbd_tt ** tntt.eye([nl]*3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the excitation on the bottom port in TT format. \n",
    "The tensor is extended using Kronecker product since the excitation does not change with the parameters.\n",
    "The boundary projectpor is applied to let only the desired dofs nonzero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "f_tt = tntt.zeros(Stt.N)\n",
    "\n",
    "excitation_dofs = tt_iga.Function(Basis).interpolate(lambda t: tn.sin(t[:,0]*np.pi)*tn.sin(t[:,1]*np.pi))\n",
    "tmp = tn.zeros(N)\n",
    "tmp[:,:,0] = excitation_dofs[:,:,0].full()\n",
    "g_tt = Pbd_tt@ (tntt.TT(tmp) ** tntt.ones([nl]*3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the system matrix and include the boundary adaptation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "k = 49\n",
    "\n",
    "eps_solver = 1e-7\n",
    "M_tt = (Pin_tt@(Stt-k*Mass_tt)+Pbd_tt).round(1e-12)\n",
    "rhs_tt = (Pbd_tt @ g_tt).round(1e-12)\n",
    "\n",
    "M_tt = M_tt.round(1e-11)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve the problem (can solve on GPU if one is available) and report ranks, times, memory consumption, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "cuda = True\n",
    "tme_amen = datetime.datetime.now() \n",
    "if cuda and tn.cuda.is_available():\n",
    "    print('Solving in TT (GPU)...')\n",
    "    dofs_tt = tntt.solvers.amen_solve(M_tt.cuda(), rhs_tt.cuda(), x0 = tntt.ones(rhs_tt.N).cuda(), eps = eps_solver, nswp=40, kickrank=2, verbose=False, preconditioner = 'c', local_iterations=24, resets=10).cpu()\n",
    "else:\n",
    "    print('Solving in TT (CPU)...')\n",
    "    dofs_tt = tntt.solvers.amen_solve(M_tt, rhs_tt, x0 = tntt.ones(rhs_tt.N), eps = eps_solver, nswp=40, kickrank=2, verbose=False, preconditioner = 'c', local_iterations=24, resets=10)\n",
    "tme_amen = (datetime.datetime.now() -tme_amen).total_seconds() \n",
    "\n",
    "print('',flush=True)\n",
    "print('Time system solve in TT :',tme_amen)\n",
    "print('Relative residual       :', (M_tt@dofs_tt-rhs_tt).norm()/rhs_tt.norm())\n",
    "print('N                       :',N)\n",
    "print('Rank Mtt                :',M_tt.R)\n",
    "print('Rank rhstt              :',rhs_tt.R)\n",
    "print('Rank solution           :',dofs_tt.R)\n",
    "print('size stiff              :',tntt.numel(Stt)*8/1e6,' MB')\n",
    "print('size mass               :',tntt.numel(Mass_tt)*8/1e6,' MB')\n",
    "print('size system mat         :',tntt.numel(M_tt)*8/1e6,' MB')\n",
    "print('size rhstt              :',tntt.numel(rhs_tt)*8/1e6,' MB')\n",
    "print('size solution           :',tntt.numel(dofs_tt)*8/1e6,' MB, one full solution: ',np.prod(N)*8/1e6,' MB')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the solution in the reference domain for $(y_1,y_2,y_3)\\in[0,1]\\times\\{0.5\\}\\times[0,1]$ and $\\theta_1=\\theta_2=\\theta_3=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "fspace = tt_iga.Function(Basis+Basis_param)\n",
    "fspace.dofs = dofs_tt\n",
    "\n",
    "fval = fspace([tn.linspace(0,1,128),tn.tensor([0.5]),tn.linspace(0,1,128),tn.tensor([-0.2]),tn.tensor([-0.2]),tn.tensor([-0.3])])\n",
    "y1, y3 = tn.meshgrid(tn.linspace(0,1,128),tn.linspace(0,1,128))\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(y1, y3, fval.numpy().squeeze(), levels = 128)\n",
    "plt.colorbar()\n",
    "plt.xlabel(r'$y_1$')\n",
    "plt.ylabel(r'$y_3$')\n",
    "\n",
    "fval = fspace([tn.linspace(0,1,128),tn.tensor([0.5]),tn.linspace(0,1,128),tn.tensor([0.2]),tn.tensor([0.2]),tn.tensor([0.3])])\n",
    "\n",
    "plt.figure()\n",
    "plt.contourf(y1, y3, fval.numpy().squeeze(), levels = 128)\n",
    "plt.colorbar()\n",
    "plt.xlabel(r'$y_1$')\n",
    "plt.ylabel(r'$y_3$')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the domain and the solution in the physical domain for $(y_1,y_2,y_3)\\in[0,1]\\times\\{0.5\\}\\times[0,1]$ and $\\theta_1=\\theta_2=\\theta_3=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "ename": "",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31mFailed to start the Kernel. \n",
      "request to http://localhost:8888/api/sessions?1650304554540 failed, reason: getaddrinfo EEXIST localhost. \n",
      "View Jupyter <a href='command:jupyter.viewOutput'>log</a> for further details."
     ]
    }
   ],
   "source": [
    "fval = fspace([tn.linspace(0,1,128),tn.tensor([0.5]),tn.linspace(0,1,128),tn.tensor([0.0]),tn.tensor([0.0]),tn.tensor([0.0]),tn.tensor([0.0])])\n",
    "x,y,z =  geom([tn.linspace(0,1,128),tn.tensor([0.5]),tn.linspace(0,1,128),tn.tensor([0.0]),tn.tensor([0.0]),tn.tensor([0.0]),tn.tensor([0.0])])\n",
    "\n",
    "fig = geom.plot_domain([tn.tensor([0.0])]*4,[(0,1),(0,1),(0.0,1)],surface_color=None, wireframe = False,frame_color='k',n = 64)\n",
    "geom.plot_domain([tn.tensor([0.0])]*4,[(0,1),(0,1),(0.25,0.5)],fig=fig,surface_color=None, wireframe = False,frame_color='k',n = 64)\n",
    "geom.plot_domain([tn.tensor([0.0])]*4,[(0,1),(0,1),(0.75,1)],fig=fig,surface_color=None, wireframe = False,frame_color='k',n = 64)\n",
    "\n",
    "ax = fig.gca()\n",
    "C = fval.numpy().squeeze()\n",
    "norm = matplotlib.colors.Normalize(vmin=C.min(),vmax=C.max())\n",
    "C = plt.cm.jet(norm(C))\n",
    "C[:,:,-1] = 1\n",
    "ax.plot_surface(x.numpy().squeeze(), y.numpy().squeeze(), z.numpy().squeeze(), edgecolors=None, linewidth=0, facecolors = C, antialiased=True, rcount=256, ccount=256, alpha=0.5)\n",
    "fig.gca().set_xlabel(r'$x_1$', fontsize=14)\n",
    "fig.gca().set_ylabel(r'$x_2$', fontsize=14)\n",
    "fig.gca().set_zlabel(r'$x_3$', fontsize=14)\n",
    "fig.gca().view_init(15, -60)\n",
    "fig.gca().zaxis.set_rotate_label(False)\n",
    "fig.gca().set_xticks([-1.5, 0, 1.5])\n",
    "fig.gca().set_yticks([0,1])\n",
    "fig.gca().set_zticks([-3,-1.5,0,1.5])\n",
    "fig.gca().tick_params(axis='both', labelsize=14)\n",
    "fig.gca().set_box_aspect(aspect = (3,1.5,4.5))\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "df6fc3a9b7a9c6f4b0308ab6eb361a4cabbf6b5db181383d07014ff4304e5cb3"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 64-bit ('pytorch': conda)",
   "language": "python",
   "name": "python3"
  },
  "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}