examples/twisted.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Twisted quad"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch as tn\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the the curves defining the edges along $x_3$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "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,w = tn.pi,v = 3):\n",
    "    phi = 0\n",
    "    r = 0.5\n",
    "    x = r*tn.cos(w*t+phi)\n",
    "    y = r*tn.sin(w*t+phi)\n",
    "    z = v * t\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))\n",
    "\n",
    "def curve1(t,w = tn.pi, v = 3):\n",
    "    phi = -tn.pi/2\n",
    "    r = 0.5\n",
    "    x = r*tn.cos(w*t+phi)\n",
    "    y = r*tn.sin(w*t+phi)\n",
    "    z = v * t\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))\n",
    "  \n",
    "def curve3(t,w = tn.pi, v = 3):\n",
    "    phi = tn.pi/2\n",
    "    r = 0.5\n",
    "    x = r*tn.cos(w*t+phi)\n",
    "    y = r*tn.sin(w*t+phi)\n",
    "    z = v * t\n",
    "    return tn.hstack((tn.reshape(x,[-1,1]),tn.reshape(y,[-1,1]),tn.reshape(z,[-1,1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the B-spline spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the parametrization and create the geometry object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scale_mult = 1\n",
    "xparam = lambda t : plane_spanner(curve1(t[:,2]),curve2(t[:,2]),curve3(t[:,2]),t[:,0],t[:,1])[0]\n",
    "yparam = lambda t : plane_spanner(curve1(t[:,2]),curve2(t[:,2]),curve3(t[:,2]),t[:,0],t[:,1])[1]\n",
    "zparam = lambda t : plane_spanner(curve1(t[:,2]),curve2(t[:,2]),curve3(t[:,2]),t[:,0],t[:,1])[2]\n",
    "\n",
    "# interpolate the geometry parametrization\n",
    "geom = tt_iga.Geometry(Basis)\n",
    "geom.interpolate([xparam, yparam, zparam])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the interpoalted TT control points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(geom.Xs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = geom.plot_domain([],[(0,1),(0,1),(0.0,1)],surface_color='blue', wireframe = True,alpha=0.1,n=64)\n",
    "fig.gca().set_xlabel(r'$x_1$')\n",
    "fig.gca().set_ylabel(r'$x_2$')\n",
    "fig.gca().set_zlabel(r'$x_3$')\n",
    "fig.gca().view_init(25, -60)\n",
    "fig.gca().set_box_aspect(aspect = (1,1,3))\n",
    "plt.xticks([-0.5,0,0.5])\n",
    "plt.yticks([-0.5,0,0.5])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute stiffness and mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Mass_tt = geom.mass_interp(eps=1e-11)\n",
    "Stiff_tt = geom.stiffness_interp( func=None, func_reference = None, qtt = False, verb=False)\n",
    "print(Mass_tt)\n",
    "print(Stiff_tt)"
   ]
  }
 ],
 "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
}