{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f389f465-1007-443d-873b-d27bfa74cbc7",
   "metadata": {},
   "source": [
    "# Composite vector space tutorial\n",
    "\n",
    "In romtools, a composite vector space refers to multiple vector spaces that are stichted together into one vector space. \n",
    "The use case of a composite vector space is if you want to have independent basis functions for different DOFs (e.g., mass, momentum)\n",
    "In this tutorial you will learn:\n",
    "- How to construct a composite vector space that comprises several vector spaces stiched together\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "358fb709-9c81-45a3-a871-2549de121580",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "module 'mpi4py' is not installed\n"
     ]
    }
   ],
   "source": [
    "#First, let's import the relavant modules:\n",
    "import romtools\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from romtools import vector_space, composite_vector_space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b1ffa4be-246c-44a9-b486-f517b9f7acb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now, we will load in snapshots from a FOM. Here, we use pre-computed snapshots of the 1D Euler equations obtained using pressio-demo-apps\n",
    "snapshots = np.load('snapshots.npz')['snapshots']\n",
    "\n",
    "## The snapshots are in tensor form:\n",
    "n_vars, nx, nt = snapshots.shape\n",
    "\n",
    "# Along the first axis we have snapshots of mass (rho), momentum (rho U), and energy (rho E), respectively\n",
    "\n",
    "## Note that romtools works with tensor forms (https://pressio.github.io/rom-tools-and-workflows/romtools/vector_space.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "217759fb-6ec9-4488-897b-e5b4db9f4b2c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# In this example, we make a composite vector space where each state variable in our snapshot matrix has each own basis\n",
    "# (This is often times referred to as \"vector POD\" in the literature\n",
    "\n",
    "#Like the last tuorial, let's use a truncater that controls for how we want to truncate our basis.\n",
    "my_truncater = vector_space.utils.EnergyBasedTruncater(0.999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "8dddecd2-01f7-4b7c-955e-456a1974fefb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The vector space for rho is of shape:  (1, 500, 4)\n"
     ]
    }
   ],
   "source": [
    "#Now, let's construct a vector space using POD for the DENSITY VARIABLE ONLY\n",
    "my_vector_space_rho = vector_space.VectorSpaceFromPOD(snapshots[0:1],truncater=my_truncater)\n",
    "\n",
    "#Let's look a the shape of our vector space:\n",
    "print('The vector space for rho is of shape: ', my_vector_space_rho.extents())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "90e36c9f-9587-479c-b943-a4857d30129d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The vector space for rhoU is of shape:  (1, 500, 31)\n",
      "The vector space for rhoE is of shape:  (1, 500, 6)\n"
     ]
    }
   ],
   "source": [
    "# Now, we will make vector spaces for the next two variables (momentum and energy)\n",
    "my_vector_space_rhoU = vector_space.VectorSpaceFromPOD(snapshots[1:2],truncater=my_truncater)\n",
    "print('The vector space for rhoU is of shape: ', my_vector_space_rhoU.extents())\n",
    "my_vector_space_rhoE = vector_space.VectorSpaceFromPOD(snapshots[2:3],truncater=my_truncater)\n",
    "print('The vector space for rhoE is of shape: ', my_vector_space_rhoE.extents())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9555169a-4579-489b-a9b8-21e2c3a67545",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The composite vector space is of shape:  [  3 500  41]\n"
     ]
    }
   ],
   "source": [
    "# Now we will create a composite vector space that stiches these three vector spaces together:\n",
    "my_composite_vector_space = composite_vector_space.CompositeVectorSpace([my_vector_space_rho,my_vector_space_rhoU,my_vector_space_rhoE])\n",
    "\n",
    "# Let's look at the size of this vector space:\n",
    "print('The composite vector space is of shape: ' , my_composite_vector_space.extents())\n",
    "#We see that the number of basis vectors is the sum of the individual vector spaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "81b37134-431e-4f77-b214-31fea9df8199",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdkAAAGiCAYAAACrorwAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAn2UlEQVR4nO3df3RU9Z3/8deEJENIk1sSTIaRQONuqmgi64Yugq5QgayuAbuerSg0pUdWyxbQLGCVdXvAnkooZ4u2hxW17Sl2rRvPnhLW7drUuMUgBxAMZA2w/jpmIWBirBvuJPyYhOTz/aPL/XYSfuTHfJzJ5Pk453N07n3P5fOeD8wrd+bOxGeMMQIAAFGXFOsJAACQqAhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsGVEh+/TTTys/P1+jR49WcXGx3njjjVhPqd/WrVsnn88XMQKBQKyndUk7d+7UvHnzFAwG5fP5tH379oj9xhitW7dOwWBQaWlpmjVrlg4fPhybyV7C5fr4xje+0WdtbrzxxthM9gIqKir0pS99SRkZGcrJydFXvvIVvfvuuxE1w2Et+tNHvK/Fli1bdP311yszM1OZmZmaPn26fv3rX3v7h8M6XK6HeF+DC6moqJDP51N5ebm3LVprMWJC9qWXXlJ5ebkee+wxHTx4UH/+53+u22+/XceOHYv11PrtuuuuU3NzszcaGhpiPaVLOnXqlKZMmaLNmzdfcP/GjRu1adMmbd68Wfv371cgENDcuXPV3t7+Gc/00i7XhyTddtttEWvzyiuvfIYzvLTa2lotW7ZMe/fuVU1Njc6dO6eSkhKdOnXKqxkOa9GfPqT4XosJEyZow4YNeuutt/TWW2/p1ltv1Z133uk9eQ+HdbhcD1J8r0Fv+/fv13PPPafrr78+YnvU1sKMEH/2Z39mli5dGrHtmmuuMY8++miMZjQwa9euNVOmTIn1NAZNkqmqqvJu9/T0mEAgYDZs2OBtO3v2rHEcxzzzzDMxmGH/9O7DGGMWL15s7rzzzpjMZzBaW1uNJFNbW2uMGb5r0bsPY4bfWhhjzNixY81PfvKTYbsOxvz/HowZXmvQ3t5uCgoKTE1NjZk5c6Z56KGHjDHR/TcxIs5kOzs7VVdXp5KSkojtJSUl2r17d4xmNXDvv/++gsGg8vPzdc899+jDDz+M9ZQGrbGxUS0tLRFr4vf7NXPmzGG1Jue9/vrrysnJ0Re/+EXdf//9am1tjfWULsp1XUlSVlaWpOG7Fr37OG+4rEV3d7cqKyt16tQpTZ8+fViuQ+8ezhsua7Bs2TLdcccdmjNnTsT2aK5FclRmGud+97vfqbu7W7m5uRHbc3Nz1dLSEqNZDcy0adP085//XF/84hf18ccf63vf+55mzJihw4cPKzs7O9bTG7Dzj/uF1uTo0aOxmNKg3X777frqV7+qSZMmqbGxUd/5znd06623qq6uTn6/P9bTi2CM0cqVK3XzzTersLBQ0vBciwv1IQ2PtWhoaND06dN19uxZfe5zn1NVVZWuvfZa78l7OKzDxXqQhscaSFJlZaUOHDig/fv399kXzX8TIyJkz/P5fBG3jTF9tsWr22+/3fv/oqIiTZ8+XX/0R3+k559/XitXrozhzIZmOK/JeQsWLPD+v7CwUFOnTtWkSZP0H//xH7rrrrtiOLO+li9frrffflu7du3qs284rcXF+hgOa3H11Vervr5eJ0+e1C9/+UstXrxYtbW13v7hsA4X6+Haa68dFmvQ1NSkhx56SK+++qpGjx590bporMWIeLl43LhxGjVqVJ+z1tbW1j4/qQwX6enpKioq0vvvvx/rqQzK+SujE2lNzhs/frwmTZoUd2uzYsUKvfzyy9qxY4cmTJjgbR9ua3GxPi4kHtciNTVVf/zHf6ypU6eqoqJCU6ZM0Q9/+MNhtQ4X6+FC4nEN6urq1NraquLiYiUnJys5OVm1tbX60Y9+pOTkZO/xjsZajIiQTU1NVXFxsWpqaiK219TUaMaMGTGa1dCEw2H993//t8aPHx/rqQxKfn6+AoFAxJp0dnaqtrZ22K7JeZ9++qmampriZm2MMVq+fLm2bdum3/72t8rPz4/YP1zW4nJ9XEi8rcWFGGMUDoeHzTpcyPkeLiQe12D27NlqaGhQfX29N6ZOnapFixapvr5eV111VfTWYsiXZw0TlZWVJiUlxfz0pz81R44cMeXl5SY9Pd38z//8T6yn1i+rVq0yr7/+uvnwww/N3r17TWlpqcnIyIjr+be3t5uDBw+agwcPGklm06ZN5uDBg+bo0aPGGGM2bNhgHMcx27ZtMw0NDebee+8148ePN6FQKMYzj3SpPtrb282qVavM7t27TWNjo9mxY4eZPn26ufLKK+Omj7/92781juOY119/3TQ3N3vj9OnTXs1wWIvL9TEc1mLNmjVm586dprGx0bz99tvm7//+701SUpJ59dVXjTHDYx0u1cNwWIOL+cOri42J3lqMmJA1xph/+qd/MpMmTTKpqanmT//0TyMu/Y93CxYsMOPHjzcpKSkmGAyau+66yxw+fDjW07qkHTt2GEl9xuLFi40xv79Mfu3atSYQCBi/329uueUW09DQENtJX8Cl+jh9+rQpKSkxV1xxhUlJSTETJ040ixcvNseOHYv1tD0Xmrsk87Of/cyrGQ5rcbk+hsNa3Hfffd5z0BVXXGFmz57tBawxw2MdLtXDcFiDi+kdstFaC58xxgzyjBsAAFzCiHhPFgCAWCBkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMCSuA/Zp59+Wvn5+Ro9erSKi4v1xhtvDPpY4XBY69atu+jXfw0HidCDlBh90EP8SIQ+6CE+RLuHuP4yipdeekllZWV6+umnddNNN+nZZ5/VT37yEx05ckQTJ04c8PFCoZAcx5HrusrMzLQwY/sSoQcpMfqgh/iRCH3QQ3yIdg9xfSa7adMmLVmyRH/zN3+jyZMn66mnnlJeXp62bNkS66kBAHBZcRuynZ2dqquri/jN9JJUUlIy4N9MDwBALMTtL23/3e9+p+7u7gv+Zvrev+NP+v3r6H/4GnpPT4/+93//V9nZ2d4v2Q2FQhH/HY4SoQcpMfqgh/iRCH3QQ3y4UA/GGLW3tysYDCopaWDnpnEbsuf19zfTV1RU6PHHH+/XMfPy8qIyt1hKhB6kxOiDHuJHIvRBD/HhQj00NTVpwoQJAzpO3IbsuHHjNGrUqH7/Zvo1a9Zo5cqV3m3XdQd1cdRnxXXdWE8BANAPoVBIeXl5ysjIGPB94zZkU1NTVVxcrJqaGv3VX/2Vt72mpkZ33nlnn3q/3y+/3/9ZTnFIhuuVdwAwUl3oVdTLiduQlaSVK1eqrKxMU6dO1fTp0/Xcc8/p2LFjWrp0aaynBgDAZcV1yC5YsECffvqpvvvd76q5uVmFhYV65ZVXNGnSpFhPDQCAy4rrL6MYivMfKI5XCfqwA0DCGcoXVMTt52QBABjuCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsCTqIVtRUaEvfelLysjIUE5Ojr7yla/o3XffjagxxmjdunUKBoNKS0vTrFmzdPjw4YiacDisFStWaNy4cUpPT9f8+fN1/PjxaE8XAABroh6ytbW1WrZsmfbu3auamhqdO3dOJSUlOnXqlFezceNGbdq0SZs3b9b+/fsVCAQ0d+5ctbe3ezXl5eWqqqpSZWWldu3apY6ODpWWlqq7uzvaUwYAwAqfMcbY/AM++eQT5eTkqLa2VrfccouMMQoGgyovL9cjjzwi6fdnrbm5ufr+97+vb37zm3JdV1dccYX++Z//WQsWLJAkffTRR8rLy9Mrr7yiv/iLv7jsnxsKheQ4js3WhsTyww4AiJLzeeK6rjIzMwd0X+vvybquK0nKysqSJDU2NqqlpUUlJSVejd/v18yZM7V7925JUl1dnbq6uiJqgsGgCgsLvZrewuGwQqFQxAAAIJashqwxRitXrtTNN9+swsJCSVJLS4skKTc3N6I2NzfX29fS0qLU1FSNHTv2ojW9VVRUyHEcb+Tl5UW7HQAABsRqyC5fvlxvv/22/uVf/qXPPp/PF3HbGNNnW2+XqlmzZo1c1/VGU1PT4CcOAEAUWAvZFStW6OWXX9aOHTs0YcIEb3sgEJCkPmekra2t3tltIBBQZ2en2traLlrTm9/vV2ZmZsQAACCWoh6yxhgtX75c27Zt029/+1vl5+dH7M/Pz1cgEFBNTY23rbOzU7W1tZoxY4Ykqbi4WCkpKRE1zc3NOnTokFcDAEC8S472AZctW6YXX3xR//Zv/6aMjAzvjNVxHKWlpcnn86m8vFzr169XQUGBCgoKtH79eo0ZM0YLFy70apcsWaJVq1YpOztbWVlZWr16tYqKijRnzpxoTxkAACuiHrJbtmyRJM2aNSti+89+9jN94xvfkCR9+9vf1pkzZ/Stb31LbW1tmjZtml599VVlZGR49U8++aSSk5N1991368yZM5o9e7a2bt2qUaNGRXvKAABYYf1zsrHC52QBANEQ15+TBQBgpCJkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwxHrIVlRUyOfzqby83NtmjNG6desUDAaVlpamWbNm6fDhwxH3C4fDWrFihcaNG6f09HTNnz9fx48ftz1dAACixmrI7t+/X88995yuv/76iO0bN27Upk2btHnzZu3fv1+BQEBz585Ve3u7V1NeXq6qqipVVlZq165d6ujoUGlpqbq7u21OGQCAqLEWsh0dHVq0aJF+/OMfa+zYsd52Y4yeeuopPfbYY7rrrrtUWFio559/XqdPn9aLL74oSXJdVz/96U/1gx/8QHPmzNENN9ygF154QQ0NDXrttddsTRkAgKiyFrLLli3THXfcoTlz5kRsb2xsVEtLi0pKSrxtfr9fM2fO1O7duyVJdXV16urqiqgJBoMqLCz0anoLh8MKhUIRAwCAWEq2cdDKykodOHBA+/fv77OvpaVFkpSbmxuxPTc3V0ePHvVqUlNTI86Az9ecv39vFRUVevzxx6MxfQAAoiLqZ7JNTU166KGH9MILL2j06NEXrfP5fBG3jTF9tvV2qZo1a9bIdV1vNDU1DXzyAABEUdRDtq6uTq2trSouLlZycrKSk5NVW1urH/3oR0pOTvbOYHufkba2tnr7AoGAOjs71dbWdtGa3vx+vzIzMyMGAACxFPWQnT17thoaGlRfX++NqVOnatGiRaqvr9dVV12lQCCgmpoa7z6dnZ2qra3VjBkzJEnFxcVKSUmJqGlubtahQ4e8GgAA4l3U35PNyMhQYWFhxLb09HRlZ2d728vLy7V+/XoVFBSooKBA69ev15gxY7Rw4UJJkuM4WrJkiVatWqXs7GxlZWVp9erVKioq6nMhFQAA8crKhU+X8+1vf1tnzpzRt771LbW1tWnatGl69dVXlZGR4dU8+eSTSk5O1t13360zZ85o9uzZ2rp1q0aNGhWLKQMAMGA+Y4yJ9SRsCIVCchwn1tO4qAR92AEg4ZzPE9d1B3y9D99dDACAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJcmxnoBtg/lN9p8Fn88X6ykAwJAZY2I9hbiW8CEbr6L9F5PQBoD4Q8gmiHj/aZIfAgCMRIQsPhPR/CFgpAU2r3oAwxchi2Gnv6GTKGGSKH0AIxEhi2GH0ElcI+0HKCQ+QhYX1N8nsVi8lMkT8dDE80v3rBkSDSGLIYl2GEczAGL1XmYi/BAwknoFbCJkMSSxeJKN1ZXU0f5z4/2KcIxs/P2MDkIWQxKLM5R4PyuK9yenaD5+8b4W6CvaPxjH+9/3WCNkcUH8w0lcrO3Ixg9Gny1CFsCIlwjBw9sZ8YmQBRJEPL+8F+8hlggXcMXqYjXC+NIIWWCEicVFaLF6Yo92KCZCoBCeny1CFkgQsQieeA+xeA6KWIUdZ7KfLUIWSBDx/FJmfyVCD9HGYzK88UvbAQCwhDNZAH0kwoVAUuL0EQu8DBwdhCwwwvTnyTPe32vtr1iEZ7xflNXfP5f3ZKODkI0RfnJGrPAtXXbFqteR9BgPJ7wnCwCAJVZC9sSJE/ra176m7OxsjRkzRn/yJ3+iuro6b78xRuvWrVMwGFRaWppmzZqlw4cPRxwjHA5rxYoVGjdunNLT0zV//nwdP37cxnQBALAi6iHb1tamm266SSkpKfr1r3+tI0eO6Ac/+IE+//nPezUbN27Upk2btHnzZu3fv1+BQEBz585Ve3u7V1NeXq6qqipVVlZq165d6ujoUGlpqbq7u6M9ZQAA7DBR9sgjj5ibb775ovt7enpMIBAwGzZs8LadPXvWOI5jnnnmGWOMMSdPnjQpKSmmsrLSqzlx4oRJSkoy1dXV/ZqH67pGknFdd5Cd2CWJwWAwhv0YCYaSJ1E/k3355Zc1depUffWrX1VOTo5uuOEG/fjHP/b2NzY2qqWlRSUlJd42v9+vmTNnavfu3ZKkuro6dXV1RdQEg0EVFhZ6NQAAxLuoh+yHH36oLVu2qKCgQL/5zW+0dOlSPfjgg/r5z38uSWppaZEk5ebmRtwvNzfX29fS0qLU1FSNHTv2ojW9hcNhhUKhiAEAQCxF/SM8PT09mjp1qtavXy9JuuGGG3T48GFt2bJFX//617263pebG2Muewn6pWoqKir0+OOPD3H2AABET9TPZMePH69rr702YtvkyZN17NgxSVIgEJCkPmekra2t3tltIBBQZ2en2traLlrT25o1a+S6rjeampqi0g8AAIMV9ZC96aab9O6770Zse++99zRp0iRJUn5+vgKBgGpqarz9nZ2dqq2t1YwZMyRJxcXFSklJiahpbm7WoUOHvJre/H6/MjMzIwYAADEV7auw9u3bZ5KTk80TTzxh3n//ffOLX/zCjBkzxrzwwgtezYYNG4zjOGbbtm2moaHB3HvvvWb8+PEmFAp5NUuXLjUTJkwwr732mjlw4IC59dZbzZQpU8y5c+f6NQ+uLmYwGAz7YyQYSp5YeYT+/d//3RQWFhq/32+uueYa89xzz0Xs7+npMWvXrjWBQMD4/X5zyy23mIaGhoiaM2fOmOXLl5usrCyTlpZmSktLzbFjx/o9B0KWwWAw7I+RYCh54jMmMb/dORQKyXEcua4bly8d8z2jABJBgkZIhKHkCd9dDACAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWRD1kz507p3/4h39Qfn6+0tLSdNVVV+m73/2uenp6vBpjjNatW6dgMKi0tDTNmjVLhw8fjjhOOBzWihUrNG7cOKWnp2v+/Pk6fvx4tKcLAIA9Jsq+973vmezsbPOrX/3KNDY2mn/91381n/vc58xTTz3l1WzYsMFkZGSYX/7yl6ahocEsWLDAjB8/3oRCIa9m6dKl5sorrzQ1NTXmwIED5stf/rKZMmWKOXfuXL/m4bqukWRc1412i1EhicFgMIb9GAmGkidRf4TuuOMOc99990Vsu+uuu8zXvvY1Y4wxPT09JhAImA0bNnj7z549axzHMc8884wxxpiTJ0+alJQUU1lZ6dWcOHHCJCUlmerq6n7Ng5BlMBgM+2MkGEqeRP3l4ptvvln/+Z//qffee0+S9F//9V/atWuX/vIv/1KS1NjYqJaWFpWUlHj38fv9mjlzpnbv3i1JqqurU1dXV0RNMBhUYWGhVwMAQLxLjvYBH3nkEbmuq2uuuUajRo1Sd3e3nnjiCd17772SpJaWFklSbm5uxP1yc3N19OhRryY1NVVjx47tU3P+/r2Fw2GFw2HvdigUilpPAAAMRtTPZF966SW98MILevHFF3XgwAE9//zz+sd//Ec9//zzEXU+ny/itjGmz7beLlVTUVEhx3G8kZeXN7RGAAAYoqiH7MMPP6xHH31U99xzj4qKilRWVqa/+7u/U0VFhSQpEAhIUp8z0tbWVu/sNhAIqLOzU21tbRet6W3NmjVyXdcbTU1N0W4NAIABiXrInj59WklJkYcdNWqU9xGe/Px8BQIB1dTUePs7OztVW1urGTNmSJKKi4uVkpISUdPc3KxDhw55Nb35/X5lZmZGDAAAYinq78nOmzdPTzzxhCZOnKjrrrtOBw8e1KZNm3TfffdJ+v3LxOXl5Vq/fr0KCgpUUFCg9evXa8yYMVq4cKEkyXEcLVmyRKtWrVJ2draysrK0evVqFRUVac6cOdGeMgAAdkT7UudQKGQeeughM3HiRDN69Ghz1VVXmccee8yEw2Gvpqenx6xdu9YEAgHj9/vNLbfcYhoaGiKOc+bMGbN8+XKTlZVl0tLSTGlpqTl27Fi/58FHeBgMBsP+GAmGkic+Y4xRAgqFQnIcR67rxuVLx5e7yAsAhoMEjZAIQ8kTvrsYAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAEkIWAABLCFkAACwhZAEAsISQBQDAkuRYT2CkMsbEegoAAMs4kwUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwYcsjt37tS8efMUDAbl8/m0ffv2iP3GGK1bt07BYFBpaWmaNWuWDh8+HFETDoe1YsUKjRs3Tunp6Zo/f76OHz8eUdPW1qaysjI5jiPHcVRWVqaTJ08OuEEAAGJlwCF76tQpTZkyRZs3b77g/o0bN2rTpk3avHmz9u/fr0AgoLlz56q9vd2rKS8vV1VVlSorK7Vr1y51dHSotLRU3d3dXs3ChQtVX1+v6upqVVdXq76+XmVlZYNoEQCAGDFDIMlUVVV5t3t6ekwgEDAbNmzwtp09e9Y4jmOeeeYZY4wxJ0+eNCkpKaaystKrOXHihElKSjLV1dXGGGOOHDliJJm9e/d6NXv27DGSzDvvvNOvubmuayQZ13WH0iIAYIQbSp5E9T3ZxsZGtbS0qKSkxNvm9/s1c+ZM7d69W5JUV1enrq6uiJpgMKjCwkKvZs+ePXIcR9OmTfNqbrzxRjmO49X0Fg6HFQqFIgYAALEU1a9VbGlpkSTl5uZGbM/NzdXRo0e9mtTUVI0dO7ZPzfn7t7S0KCcnp8/xc3JyvJreKioq9Pjjjw+5BwBA//l8vlhPIa5Zubq494NujLnsQvSuuVD9pY6zZs0aua7rjaampkHMHACA6IlqyAYCAUnqc7bZ2trqnd0GAgF1dnaqra3tkjUff/xxn+N/8sknfc6Sz/P7/crMzIwYAADEUlRDNj8/X4FAQDU1Nd62zs5O1dbWasaMGZKk4uJipaSkRNQ0Nzfr0KFDXs306dPluq727dvn1bz55ptyXderAQAg3g34PdmOjg598MEH3u3GxkbV19crKytLEydOVHl5udavX6+CggIVFBRo/fr1GjNmjBYuXChJchxHS5Ys0apVq5Sdna2srCytXr1aRUVFmjNnjiRp8uTJuu2223T//ffr2WeflSQ98MADKi0t1dVXXx2NvgEAsG+glyPv2LHDSOozFi9ebIz5/cd41q5dawKBgPH7/eaWW24xDQ0NEcc4c+aMWb58ucnKyjJpaWmmtLTUHDt2LKLm008/NYsWLTIZGRkmIyPDLFq0yLS1tfV7nnyEBwDsu1AeJOoYTJ74/u9BSjihUEiO48h1Xd6fBQBLRtLVxYPJE767GAAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsIWQBALCEkAUAwBJCFgAASwhZAAAsGXDI7ty5U/PmzVMwGJTP59P27du9fV1dXXrkkUdUVFSk9PR0BYNBff3rX9dHH30UcYxwOKwVK1Zo3LhxSk9P1/z583X8+PGImra2NpWVlclxHDmOo7KyMp08eXJQTQIAEAsDDtlTp05pypQp2rx5c599p0+f1oEDB/Sd73xHBw4c0LZt2/Tee+9p/vz5EXXl5eWqqqpSZWWldu3apY6ODpWWlqq7u9urWbhwoerr61VdXa3q6mrV19errKxsEC0CABAjZggkmaqqqkvW7Nu3z0gyR48eNcYYc/LkSZOSkmIqKyu9mhMnTpikpCRTXV1tjDHmyJEjRpLZu3evV7Nnzx4jybzzzjv9mpvrukaScV13gF0BAPpL0ogZg8kT6+/Juq4rn8+nz3/+85Kkuro6dXV1qaSkxKsJBoMqLCzU7t27JUl79uyR4ziaNm2aV3PjjTfKcRyvprdwOKxQKBQxAACIJashe/bsWT366KNauHChMjMzJUktLS1KTU3V2LFjI2pzc3PV0tLi1eTk5PQ5Xk5OjlfTW0VFhff+reM4ysvLi3I3AAAMjLWQ7erq0j333KOenh49/fTTl603xsjn83m3//D/L1bzh9asWSPXdb3R1NQ0+MkDABAFVkK2q6tLd999txobG1VTU+OdxUpSIBBQZ2en2traIu7T2tqq3Nxcr+bjjz/uc9xPPvnEq+nN7/crMzMzYgAAEEtRD9nzAfv+++/rtddeU3Z2dsT+4uJipaSkqKamxtvW3NysQ4cOacaMGZKk6dOny3Vd7du3z6t588035bquVwMAQLxLHugdOjo69MEHH3i3GxsbVV9fr6ysLAWDQf31X/+1Dhw4oF/96lfq7u723kPNyspSamqqHMfRkiVLtGrVKmVnZysrK0urV69WUVGR5syZI0maPHmybrvtNt1///169tlnJUkPPPCASktLdfXVV0ejbwAA7Bvo5cg7duy44KXNixcvNo2NjRe99HnHjh3eMc6cOWOWL19usrKyTFpamiktLTXHjh2L+HM+/fRTs2jRIpORkWEyMjLMokWLTFtbW7/nyUd4AMC+iz3nJ+IYTJ74/u9BSjihUEiO48h1Xd6fBQBLLnYxaiIaTJ7w3cUAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFhCyAIAYMmAQ3bnzp2aN2+egsGgfD6ftm/fftHab37zm/L5fHrqqacitofDYa1YsULjxo1Tenq65s+fr+PHj0fUtLW1qaysTI7jyHEclZWV6eTJkwOdLgAAMTPgkD116pSmTJmizZs3X7Ju+/btevPNNxUMBvvsKy8vV1VVlSorK7Vr1y51dHSotLRU3d3dXs3ChQtVX1+v6upqVVdXq76+XmVlZQOdLgAAsWOGQJKpqqrqs/348ePmyiuvNIcOHTKTJk0yTz75pLfv5MmTJiUlxVRWVnrbTpw4YZKSkkx1dbUxxpgjR44YSWbv3r1ezZ49e4wk88477/Rrbq7rGknGdd3BNQcAuCxJI2YMJk+i/p5sT0+PysrK9PDDD+u6667rs7+urk5dXV0qKSnxtgWDQRUWFmr37t2SpD179shxHE2bNs2rufHGG+U4jlcDAEC8S472Ab///e8rOTlZDz744AX3t7S0KDU1VWPHjo3Ynpubq5aWFq8mJyenz31zcnK8mt7C4bDC4bB3OxQKDbYFAACiIqpnsnV1dfrhD3+orVu3yufzDei+xpiI+1zo/r1r/lBFRYV3kZTjOMrLyxvY5AEAiLKohuwbb7yh1tZWTZw4UcnJyUpOTtbRo0e1atUqfeELX5AkBQIBdXZ2qq2tLeK+ra2tys3N9Wo+/vjjPsf/5JNPvJre1qxZI9d1vdHU1BTN1gAAGLCohmxZWZnefvtt1dfXeyMYDOrhhx/Wb37zG0lScXGxUlJSVFNT492vublZhw4d0owZMyRJ06dPl+u62rdvn1fz5ptvynVdr6Y3v9+vzMzMiAEAQCwN+D3Zjo4OffDBB97txsZG1dfXKysrSxMnTlR2dnZEfUpKigKBgK6++mpJkuM4WrJkiVatWqXs7GxlZWVp9erVKioq0pw5cyRJkydP1m233ab7779fzz77rCTpgQceUGlpqXccAADi3YBD9q233tKXv/xl7/bKlSslSYsXL9bWrVv7dYwnn3xSycnJuvvuu3XmzBnNnj1bW7du1ahRo7yaX/ziF3rwwQe9q5Dnz59/2c/mAgAQT3z/9zmnhBMKheQ4jlzX5aVjALBkoBe5DmeDyRO+uxgAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEsIWQAALCFkAQCwhJAFAMASQhYAAEuSYz0BW4wxkqRQKBTjmQAAEsH5XBmIhA3Z9vZ2SVJeXl6MZwIASATt7e1yHGdA9/GZwUTzMNDT06OPPvpIGRkZ8vl8kn5/VpuXl6empiZlZmbGeIaDkwg9SInRBz3Ej0Togx7iw4V6MMaovb1dwWBQSUkDe5c1Yc9kk5KSNGHChAvuy8zMHLZ/Ac5LhB6kxOiDHuJHIvRBD/Ghdw8DPYM9jwufAACwhJAFAMCSERWyfr9fa9euld/vj/VUBi0RepASow96iB+J0Ac9xIdo95CwFz4BABBrI+pMFgCAzxIhCwCAJYQsAACWELIAAFhCyAIAYAkhCwCAJYQsAACWELIAAFjy/wBxSoB8GYKGKQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Let's view the sparsity pattern if we view this as a tradititional basis matrix\n",
    "#First, let's grab the basis\n",
    "composite_basis = my_composite_vector_space.get_basis()\n",
    "#Now we will reshape it into a matrix. Note all of romtools uses order 'C' for reshaping\n",
    "basis_matrix = np.reshape(composite_basis,(composite_basis.shape[0]*composite_basis.shape[1],composite_basis.shape[2]),'C')\n",
    "\n",
    "plt.spy(basis_matrix)\n",
    "plt.gca().set_aspect(0.025)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb18a534-c204-4290-8957-e3830d399d2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We see the basis has a block structure. This means that each state variable will effectively have their own reduced coordinates."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}