{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Triangle Mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from lapy import TriaMesh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we create and initialize a TriaMesh of a cube. Note, this mesh is not oriented (some normals flip), so some functions below will complain until we fix the situation with .orient_() later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = [[0.0,0.0,0.0], [0,1,0], [1,1,0], [1,0,0],\n",
    "          [0,0,1], [0,1,1], [1,1,1], [1,0,1]]\n",
    "trias = [[0,1,2], [2,3,0], [4,5,6], [6,7,4], [0,4,7], [7,3,0],\n",
    "         [0,4,5], [5,1,0], [1,5,6], [6,2,1], [3,7,6], [6,2,3]]\n",
    "\n",
    "T = TriaMesh(points,trias)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we check a few properties of the mesh. Is it closed (no boundary)? Is it manifold (no edge has more than two triangles)? And is it consistently oriented (it won't be)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.is_closed()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.is_manifold()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.is_oriented()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that orientation is not OK also by looking at the max entry of the directed adjacency matrix. It is a sparse matrix (mainly zeros) adding a 1 for an edge from i to j. If any position has a 2 it means that two half edges are oriented in the same direction, which can only happen if the mesh is not oriented consistently or if it is non-manifold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.0\n"
     ]
    }
   ],
   "source": [
    "print(np.max(T.adj_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many functions work on non-oriented meshes. For example, we can compute the Euler characteristics (2 for spherical topology), the areas of the triangles and the total mesh area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.euler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.tria_areas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.999999999999998"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.area()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The volume inside the closed mesh, however, requires the mesh to be correctly oriented."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "Error: Can only compute volume for oriented triangle meshes!",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[11], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mT\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mvolume\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m  \u001b[38;5;66;03m# will fail for non-oriented meshes\u001b[39;00m\n",
      "File \u001b[0;32m~/Work/software/lapy/LaPy/examples/lapy/tria_mesh.py:304\u001b[0m, in \u001b[0;36mTriaMesh.volume\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    302\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;241m0.0\u001b[39m\n\u001b[1;32m    303\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mis_oriented():\n\u001b[0;32m--> 304\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m    305\u001b[0m         \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mError: Can only compute volume for oriented triangle meshes!\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m    306\u001b[0m     )\n\u001b[1;32m    307\u001b[0m v0 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mv[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mt[:, \u001b[38;5;241m0\u001b[39m], :]\n\u001b[1;32m    308\u001b[0m v1 \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mv[\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mt[:, \u001b[38;5;241m1\u001b[39m], :]\n",
      "\u001b[0;31mValueError\u001b[0m: Error: Can only compute volume for oriented triangle meshes!"
     ]
    }
   ],
   "source": [
    "T.volume()  # will fail for non-oriented meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But we can get the vertex degrees (number of edges at each vertex)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6, 4, 4, 4, 4, 4, 6, 4])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.vertex_degrees()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The area that a vertex represents is computed by summing the areas of all triangles at that vertex and dividing it by 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.        , 0.66666667, 0.66666667, 0.66666667, 0.66666667,\n",
       "       0.66666667, 1.        , 0.66666667])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.vertex_areas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course their sum should be equal to the total surface area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(T.vertex_areas())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also the average edge length can be computed on non-oriented meshes easily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1380711874576983"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.avg_edge_length()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the triangle normals, but some will point outwards and some inwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  0., -1.],\n",
       "       [ 0., -0., -1.],\n",
       "       [ 0.,  0., -1.],\n",
       "       [ 0., -0., -1.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [-1.,  0.,  0.],\n",
       "       [-1.,  0., -0.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [ 0.,  1.,  0.],\n",
       "       [-1.,  0.,  0.],\n",
       "       [-1.,  0., -0.]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.tria_normals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Triangle qualities measure how close a triangle is to the equilateral triangle. Triangles with a very short edge (sliver) are usually numerically unstable when solving equations with FEM and therefore get a low quality. Here all triangles are similar and have the same quality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.8660254, 0.8660254, 0.8660254, 0.8660254, 0.8660254, 0.8660254,\n",
       "       0.8660254, 0.8660254, 0.8660254, 0.8660254, 0.8660254, 0.8660254])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.tria_qualities()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since all vertices are used in at least one triangle, we do not have free vertices. Usually free vertices exist if the mesh is obtained as the boundary of a tetrahedral mesh and the inner vertices are not removed. This keeps vertex indices fixed and allows mapping results from the surface back to the solid. Anyway, removing free vertices in our case will not change anything."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.has_free_vertices()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 1, 2, 3, 4, 5, 6, 7]), [])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.rm_free_vertices_()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now orient the mesh consistently so that all triangle normals point outwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Searched mesh after 3 flooding iterations (0.0007233619689941406 sec).\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.orient_()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T.is_oriented()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can also compute the enclosed volume."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute Volume (works only for oriented meshes)\n",
    "T.volume() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once triangle normals are oriented we can also compute vertex normals (by aggregating the triangle normals at a vertex weighted by the angle)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.57735027, -0.57735027, -0.57735027],\n",
       "       [-0.40824829,  0.81649658, -0.40824829],\n",
       "       [ 0.40824829,  0.40824829, -0.81649658],\n",
       "       [ 0.81649658, -0.40824829, -0.40824829],\n",
       "       [-0.40824829, -0.40824829,  0.81649658],\n",
       "       [-0.81649658,  0.40824829,  0.40824829],\n",
       "       [ 0.57735027,  0.57735027,  0.57735027],\n",
       "       [ 0.40824829, -0.81649658,  0.40824829]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Vertex normals averages tria normals (only oriented meshes)\n",
    "T.vertex_normals()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vertex normals are important for rendering the mesh, or they can be used to shift vertices along the normals by an offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "T.normal_offset_(0.2*T.avg_edge_length())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far our mesh was closed, here we create one with a boundary by dropping triangles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tbdr = TriaMesh(T.v,T.t[2:,:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Meshes can be refined by subdividing each triangle into 4 similar smaller ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "Tbdr.refine_()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we print the boundary loops, these are the vertices oredered along the boundary. Here we only have one loop, but meshes can have many."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 8, 1, 13, 2, 16, 3, 9]]\n"
     ]
    }
   ],
   "source": [
    "print(Tbdr.boundary_loops())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python3",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}