epfl-archive/cs457-gc/assignment_2_4/notebook/inverse_design_dino.ipynb

475 lines
34 KiB
Plaintext
Raw Permalink Normal View History

2022-04-07 18:46:57 +02:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import torch\n",
"import igl\n",
"\n",
"import meshplot as mp\n",
"import sys as _sys\n",
"_sys.path.append(\"../src\")\n",
"from elasticenergy import *\n",
"from elasticsolid import *\n",
"from adjoint_sensitivity import *\n",
"from vis_utils import *\n",
"from objectives import *\n",
"from harmonic_interpolator import *\n",
"from shape_optimizer import *\n",
"\n",
"from utils import *\n",
"\n",
"shadingOptions = {\n",
" \"flat\":True,\n",
" \"wireframe\":False, \n",
"}\n",
"\n",
"rot = np.array(\n",
" [[1, 0, 0 ],\n",
" [0, 0, 1],\n",
" [0, -1, 0 ]]\n",
")\n",
"\n",
"torch.set_default_dtype(torch.float64)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Create the deformed object\n",
"\n",
"## Load the mesh"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "35cf42b4eded4cb99f3811a16dcb4d67",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-1.987469…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"<meshplot.Viewer.Viewer at 0x7fa5b879cc10>"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"vNP, _, _, tNP, _, _ = igl.read_obj(\"../data/dinosaur.obj\")\n",
"# vNP, _, _, tNP, _, _ = igl.read_obj(\"../data/beam.obj\")\n",
"\n",
"aabb = np.max(vNP, axis=0) - np.min(vNP, axis=0)\n",
"length_scale = np.mean(aabb)\n",
"\n",
"\n",
"v, t = torch.tensor(vNP), torch.tensor(tNP)\n",
"eNP = igl.edges(tNP)\n",
"beNP = igl.edges(igl.boundary_facets(tNP))\n",
"\n",
"bvNP, ivNP = get_boundary_and_interior(v.shape[0], tNP)\n",
"\n",
"mp.plot(vNP @ rot.T, np.array(tNP), shading=shadingOptions)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Add some physical characteristics"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pinned vertices: [ 46 47 50 59 60 62 64 65 88 89 91 98 99 102 103 104]\n"
]
}
],
"source": [
"rho = 131 # [kg.m-3], if aabb[0] ~ 14m, and m_tot = 6000kg\n",
"young = 3e8 # [Pa] \n",
"poisson = 0.2\n",
"\n",
"# Find some of the lowest vertices and pin them\n",
"minZ = torch.min(v[:, 2])\n",
"pin_idx = torch.arange(v.shape[0])[v[:, 2] < minZ + 0.01*aabb[2]]\n",
"vIdx = np.arange(v.shape[0])\n",
"pin_idx = vIdx[np.in1d(vIdx, bvNP) & np.in1d(vIdx, pin_idx)]\n",
"print(\"Pinned vertices: {}\".format(pin_idx))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Initial guess\n",
"\n",
"The idea is that we start deforming the mesh by inverting gravity."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "5d62f192ccc940388718972e4a447f2b",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-2.468080…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Inverted gravity\n",
"force_mass = torch.zeros(size=(3,))\n",
"force_mass[2] = + rho * 9.81\n",
"\n",
"# Gravity going in the wrong direction\n",
"\n",
"ee = NeoHookeanElasticEnergy(young, poisson)\n",
"\n",
"v = HarmonicInterpolator(v, t, ivNP).interpolate(v[bvNP])\n",
"solid_init = ElasticSolid(v, t, ee, rho=rho, pin_idx=pin_idx, f_mass=force_mass)\n",
"\n",
"solid_init.find_equilibrium()\n",
"plot_torch_solid(solid_init, beNP, rot, length_scale)\n",
"\n",
"# Use these as initial guesses\n",
"v_init_rest = solid_init.v_def.clone().detach()\n",
"v_init_def = solid_init.v_rest.clone().detach()\n",
"\n",
"# v_init_rest = solid_init.v_rest.clone().detach()\n",
"# v_init_def = solid_init.v_def.clone().detach()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Inverse design\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Initial objective: 6.6351e+02\n",
"\n"
]
}
],
"source": [
"force_mass = torch.zeros(size=(3,))\n",
"force_mass[2] = - rho * 9.81\n",
"use_linear = False\n",
"\n",
"# The target is the initial raw mesh\n",
"vt_surf = torch.tensor(vNP[bvNP, :])\n",
"\n",
"# Create solid\n",
"if use_linear:\n",
" ee = LinearElasticEnergy(young, poisson)\n",
"else:\n",
" ee = NeoHookeanElasticEnergy(young, poisson)\n",
"solid_ = ElasticSolid(v_init_rest, t, ee, rho=rho, pin_idx=pin_idx, f_mass=force_mass)\n",
"solid_.update_def_shape(v_init_def)\n",
"\n",
"optimizer = ShapeOptimizer(solid_, vt_surf, weight_reg=0.)\n",
"\n",
"v_eq_init = optimizer.solid.v_def.clone().detach() #bookkeeping"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Objective after 19 optimization step(s): 5.8586e+02\n",
" Line search Iters: 8\n",
"Elapsed time: 683.7s. \n",
"Estimated remaining time: 755.7s\n",
"\n"
]
},
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fa1a13b5101d4629a3d07e84205409dd",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.497634…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"An exception occured: . \n",
"Halving the step.\n",
"Line search can't find a step to take\n"
]
}
],
"source": [
"optimizer.optimize(step_size_init=1e-4, max_l_iter=10, n_optim_steps=40)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAGHCAYAAAD1HvUOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABA2klEQVR4nO3deZxcVZ338c+3t+okXZ0ASTWQhD2EVRACiiB2jCPq6OCCI86Mj6Ij48gw+DzjxiyOOsO4LzjIKIIrYERcwAVGECIuLAIi+xIgkoVsQNLpLL3+nj/u7aZS6e70Ult3fd+vV72q6t5z655zqzv9zTn33qOIwMzMzMyqQ12lK2BmZmZmz3M4MzMzM6siDmdmZmZmVcThzMzMzKyKOJyZmZmZVRGHMzMzM7Mq4nBmVuUkLZN00UTLFKkuIemMUu9nspDUnh6T2RP8nG9K+mmx6jXMPg5I67qolPsxs4lzODOrEElzJV0iaZWkbkmrJX1N0rxxfNwbgfOLWLfhwsI+wE+KtZ/JRNIKSe8vWPw7kmPyzAQ//jzgbyb4GYOGCesrSep6T7H2Y2al4XBmVgGSDgTuBI4C3g4cQvLH+Ujg95IOGMvnRcSzEbGl2PUcYj9rI6Kr1PuZLCKiOz0mE7qbd0RsjohNRarWcPvoS+vaW8r9mNnEOZyZVcaXgX7gFRHxy4h4KiJuBl6RLv9yQfkGSRdKei59fEbS4O9vYU+JpCZJn0p75bZK+r2k0/I/UNJhkq6VtFlSp6RbJR0t6aMkgfHP02GwkNSebjM4rJmW/1zBZ7ZK2i7pDaOtRyFJr5L067Sdz0r6X0mHF5T5iKQ/SeqStFbSt3fzmadKul3SDknrJH1BUlPB8fvKcMdY0jJgf+AzA8ckXb7TsKakd6TH8tWSHpa0LT3GMyWdIemx9Hh/R9K0vP0P9lTmfWbhY1m6fi9J302P6XZJD0g6K/+zgJcB5+Rte8BQw5qjPC4XS/ovSRslrZf02fyfvWGO9zslPZW2/yeS3jtwzPLK/J2k5Up6jZdLenfB+plKepbXS9oi6VcFdZ+ZHsf1af2fkPS+keplNlk4nJmVmaQ9gVcBX46Ibfnr0vcXA6+WtEfeqr8m+X09Cfg74GzgfSPs5hskf6D/Cjga+BbwE0nHpHXYF/gNEMCfAceRBMJ64LPAVcCNJMNg+5AM3xW6HDiz4A/1m4DtwM9GU49hzAC+CJwItAOb022a0rq/CXg/8F5gAfBa4I7hPkzSXOA64A/AC4F3AW8FPlFQdKRj/EZgFfBxnj8mw8kA/5R+3hJgEXA1SeB9E/D6tM7vHWb7gaHSgcciYBOwLF3fDNydfsaRwIXAVyUtSdefB9xKcuwHPmNl4U7GeFx6gZcA/0ByTN4yXOMlnQRcSvLzdCxwLfCxgjJvAC4i+Z6PSttwsaTXpetF8jM0N23nC4FbgJskDRz7/yT5mXotcBjwTmD1cPUym1Qiwg8//CjjA3gRSSh6wzDr35CuPzF9vwx4FFBemX8FVuW9XwZclL4+mKT3bb+Cz/0xcHH6+gLgT0DTMHX4JvDTIZYHcEb6ei+gG1iSt/5G4Kujrccoj9cMoA84JX3//4BHgMZRbn8BsByoy1v2DqALmD6GY7wCeH/BZ7enx2R23ucGsDCvzGfT+s8e7viOcLynkQx//zC/bkOUWwpcOtTPQ96yA9K6LRrjcbm14HNuyN/XEHX5LnB9wbJLgMh7/1vg60P8zP0mff1yoBOYVlDmHuCD6etrgW8U43fSDz+q7eGeM7PKGe48JQ2x/raIyH9/KzBXUusQ2x+XfsaD6RBbp6RO4M9JAhMkPRG/iYjucVc+4hngf0l6Vkh7NBaT9KiNth67kHSwpCslPS6pA1hH0qO1X1rk+yS9R09KukzSmyVlRqjq4SQBoz9v2W+AJpJz/QaM5RiPpCsiHsl7vw5YGxEbC5blRvqQtPfomyS9mW8bqJukekn/IuleSc+kx/SNPH98Rmu0x+Xegu3W7Kbuh7FrT+btQ+z7twXLfgMckb4+HpgObCj42TmK5392/gf4S0l/TIdaXzZCncwmlYZKV8CsBj1GEryOJOlFKnR4uv7xcX5+Xbr9CUBPwbrt6bMojsuBSyS9l2RIbCXJH9nR1mMoPyEZnvq79LkXeJAkNBARKyUtJBkyfAXwOeDfJb0oIrYO8Xli+CA8oRP5h1F4wn2wa/uD3Z9W8hHgVOCEgna9n2TY9DzgPpIepv9iN2FvCKM9LmOt+0ifO9w+CpfVkQTYlw5RpgMgIq6TtD/wapKfhZ9J+n5EnDXENmaTinvOzMosIp4l6XF6r6Tp+evS9+cA16XlBrwo7UkZ8GJgTUR0DLGLP5D8gdw7IpYXPAbOybkbOCX/5O8C3SQ9NrtzTfr8WpIetCvyep9GU4+dSNqLJJz+V0TcGBEPAVkK/iMZETsi4mcR8X9Jwt+RwMnD1PFB4KSCc+NOSduYH4B3d4xHe0wmTMlFFx8ETo+IVQWrTwF+EhHfiYh7SNpwaEGZ0dR1tMdlrB4iOV8wX+H7h9J95TslrRMkP59tQP8QPzvrBzaIiI3pcXgHyTlzb99NL6rZpOBwZlYZ/0ASOG6U9HJJ85VcEXkDSaD5h4Ly+wJflLQw/cP9AeALQ31wRDwKXAF8U8kVggdJWiTp/ZLemBa7GGgBrpJ0gqRDJL1V0rHp+hXAUen+ZktqHGZfO0jOh/pXkmHMy/PWjaYehZ4DNgLvTuv0MuAr5PVGKbki8m+VXFl6IHAWSe/OY8N85sXp8btY0uGS/hz4JMk5WfkXZOzuGK8AXqrk/nQTuunsSCQdRXLhxD8DT0naO33smRZ5FFgi6RRJh5GcWH9gwcesAE5UcoXmbA19deVoj8tYfQl4paQPSFog6V0k51Hm+wzwNknnpGXOJQn3n07X30gy7HmNkitfD5R0kqSPSXopgKSPS3p9uv3hJEO7T4Rv9WJTgMOZWQVExOMkV+E9AHwHeAK4kqRH4YSIeLJgkytIekJuB74GXMYw4Sx1FsnVep8GHgZ+SjJE9qd0/6vT903AzSS9XOfyfAj6WlqXO4ENDN8rRVr/Y4C7056uUdejUHr+01uAFwD3k1zx928kJ6kP2ETSS/LrtMybgDcOccwGPnM1ydDXC0lOKP86yUnr/1xQdHfH+CPAfJJepQ1D7atIFpGcb/VF4Om8xw/T9f9Jck7XdSRXMG5N657vsyQ9YA+mdd3lfLQxHJcxiYhbgXcD/0hyvtrrgU8BO/LK/Jjk5+3/pnU8D3hvRPwkXR/Aa4CbSL6LR0iuIF5Ics4bJD8TFwB/JAlyWeB1E6m7WbXQzue/mtlkJOlW4FcR8eFK12UyUnIPsfsjorDH0opA0hdI7ul3dKXrYjYZuOfMbBKTlFFyY84jSXqRzCouHdI8Nh2afg/wHpKeOTMbBV+taTa5vRr4NskVjt+rcF3MBiwiuap0JvAkybyvF1a0RmaTiIc1zczMzKqIhzXNzMzMqojDmZmZmVkVmTLnnM2ePTsOOOCAku9n69atzJgxo+T7qUa13Hao7fa77bXZdqjt9tdy26G221+Ott91110bI2LOUOumTDg74IADuPPOO0u+n2XLltHe3l7y/VSjWm471Hb73fb2SlejYmq5/bXcdqjt9pej7ZKGvN8jeFjTzMzMrKo4nJmZmZlVEYczMzMzsyricGZmZmZWRRzOzMzMzKqIw5mZmZlZFXE4MzMzM6siDmdmZmZmVcThzMzMzKyKOJyZmZmZVRGHMzMzM7MqMmXm1iy1iGB7Tx/be4PN23uICPr6g/6A/gj60/cRpMsHHoy4buB9/rpZ05s4bO8szY31lW62mZmZlZnD2Sh1bO/lmI//Inlz4y9Kvr+GOnFoW5aj587k6HkzOXruTA7bJ0umwYHNzMxsKnM4G6XmpjrOf/VhPPHEExy64BDqBHUSdXWiTlAv7fy+TkhKl5Mu3/l14br6OpDE+o4d3Ld6M/eu2swvHlzL9+5cCSSBbeHeOwe2hXs7sJmZmU0lDmejlGmo5+9edjDLYiXtpxxY8v296qh9gGQ4ddVz27l/9WbuSx/
"text/plain": [
"<Figure size 720x432 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"plt.figure(figsize=(10, 6))\n",
"plt.plot(to_numpy(optimizer.objectives[optimizer.objectives > 0]))\n",
"plt.title(\"Objective as optimization goes\", fontsize=14)\n",
"plt.xlabel(\"Optimization steps\", fontsize=12)\n",
"plt.ylabel(\"Objective\", fontsize=12)\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Green (Initial guess for rest state) deploys to Black\n",
"\n",
"Blue (Optimized rest state) deploys to Yellow\n",
"\n",
"Red is the Target Shape\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "70afc138cf264b47b4fddfd6b9ae8624",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"Renderer(camera=PerspectiveCamera(children=(DirectionalLight(color='white', intensity=0.6, position=(-0.497634…"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
"4"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = mp.plot(np.array(optimizer.solid.v_def) @ rot.T, tNP, shading=shadingOptions)\n",
"# p.add_points(np.array(optimizer.solid.v_def)[pin_idx, :] @ rot.T, shading={\"point_color\":\"black\", \"point_size\": 0.2})\n",
"p.add_edges(np.array(v_init_rest) @ rot.T, beNP, shading={\"line_color\": \"green\"})\n",
"p.add_edges(vNP @ rot.T, beNP, shading={\"line_color\": \"red\"})\n",
"p.add_edges(np.array(v_eq_init) @ rot.T, beNP, shading={\"line_color\": \"black\"})\n",
"p.add_edges(np.array(optimizer.solid.v_rest) @ rot.T, beNP, shading={\"line_color\": \"blue\"})\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"v_rest_optim_g = optimizer.solid.v_rest.clone().detach() #bookkeeping"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Add point load to the right most vertices\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"maxX = torch.min(v[:, 0])\n",
"f_point_idx = torch.arange(v.shape[0])[v[:, 0] > maxX - 0.01*aabb[0]]\n",
"\n",
"f_point = torch.zeros(size=(f_point_idx.shape[0], 3))\n",
"f_point[:, 2] = -5e3\n",
"\n",
"optimizer.solid.add_point_load(f_point_idx, f_point)\n",
"optimizer.set_params(optimizer.params)\n",
"v_def_optim_g_under_point = optimizer.solid.v_def.clone().detach() #bookkeeping"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"optimizer.reset_BFGS()\n",
"optimizer.optimize(step_size_init=1e-4, max_l_iter=10, n_optim_steps=20)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Green (Optimum rest state under gravity) deploys to Black with the additional point load\n",
"\n",
"Blue (Optimized rest state) deploys to Yellow\n",
"\n",
"Red is the Target Shape\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p = mp.plot(np.array(optimizer.solid.v_def) @ rot.T, tNP, shading=shadingOptions)\n",
"# p.add_points(np.array(optimizer.solid.v_def)[pin_idx, :] @ rot.T, shading={\"point_color\":\"black\", \"point_size\": 0.2})\n",
"p.add_edges(np.array(v_rest_optim_g) @ rot.T, beNP, shading={\"line_color\": \"green\"})\n",
"p.add_edges(vNP @ rot.T, beNP, shading={\"line_color\": \"red\"})\n",
"p.add_edges(np.array(v_def_optim_g_under_point) @ rot.T, beNP, shading={\"line_color\": \"black\"})\n",
"p.add_edges(np.array(optimizer.solid.v_rest) @ rot.T, beNP, shading={\"line_color\": \"blue\"})\n"
]
}
],
"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.9.7"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}