# Nonlinear Discrete Shell Energies - The Easy Way
This package provides a stand-alone way to use the nonlinear shell deformation energies from [GOAST](https://gitlab.com/numod/goast) with triangle meshes described by matrices as used by libigl. 
Furthermore, it provides convenient interfaces for MATLAB ([gptoolbox](https://github.com/alecjacobson/gptoolbox)) and Python ([libigl](https://github.com/libigl/libigl-python-bindings)).

For any questions, feel free to contact [Josua Sassen](https://josuasassen.com).

## Setup instructions

## Compiling
Independent of your OS, you will first need to install a C++ compiler, CMake, and Boost.
Then, to compile, you have to follow the usual steps of compiling a CMake-based project in a terminal (PowerShell on Windows):
```bash
cd <SHELL_DIR>
mkdir build
cd build
cmake ..
cmake --build . --config Release
```
Here, `<SHELL_DIR>` needs to be replaced with the directory in which you cloned the git repository.
This will compile the C++ library that is directly compatible with libigl as demonstrated by the included example.

On Windows, don't worry if a range of warnings appears during compilation (last step), that is normal.

### MATLAB interface
To use this package with MATLAB you will need to compile it, similar to the mex modules of [gptoolbox](https://github.com/alecjacobson/gptoolbox).
To this end, you have to explicitly tell CMake that you want the MATLAB interface by setting the `BUILD_MEX` option to `ON`.
This means, in the above instructions, you have to replace the last two steps with
```bash
cmake -DBUILD_MEX=ON ..
cmake --build . --config Release
```
This compiles the MATLAB interface together with the rest of the package.

On Linux and macOS, this will automatically put the mex-files for MATLAB in the `mex` directory.
On Windows, to move the mex/files in the right directory, you need to additionally execute the following command
```bash
cp ../mex/Release/*.mexw64 ../mex
```

Sometimes CMake has a troubles finding MATLAB.
In this case, you need to give the path to MATLAB explicitly by adding `-DMatlab_ROOT_DIR=/path/to/matlab-version/` to the `cmake ..` command.
For more helpful tips, also check out [this](https://github.com/odedstein/sgi-introduction-course/blob/main/205_mex/compilation_instructions.md) tutorial on compiling gptoolbox, which is very similar.


### Python interface
Similarly to the MATLAB interface, you also have to explicitly tell CMake that you want the Python interface by setting the `BUILD_PYTHON` option to `ON`.
This means, in the above instructions, you have to replace the last two steps with
```bash
cmake -DBUILD_PYTHON=ON ..
cmake --build . --config Release
```
Of course, you can also combine the two options!

If you are using another version of Python than is provided by your system (e.g. if you are using Anaconda), then you need to provide this information to CMake via the option `-DPython_EXECUTABLE=/path/to/python-exectuable`. 
For Anaconda, this path will typically look like `<ANACONDA_DIRECTORY>/bin/python`.

On Linux and macOS, the Python interface is provided as `pyshell.<SOMETHING>.so` in the `python`-subdirectory of your `build`-directory.
A nice installation script for this will come in the future.

## Using in MATLAB
First, add the mex directory to your MATLAB paths: (starting from here, everything is done in MATLAB!)
```matlab
addpath("<SHELL_DIR>/mex");
```

Then create two spheres (using gptoolbox) of different radii and deform the second one slightly
```matlab
[V,F] = subdivided_sphere(1);
V1 = V;
V2 = V * 1.1;
V2(3:10,:) = V2(3:10,:)*1.1;
```
for the bending energy, we furthermore need to construct edge flaps (again using gptoolbox)
```matlab
[EF,EI,uE,EMAP] = edge_flaps(F);
```
Now, we can evaluate the complete shell energy as
```matlab
en = shell_energy(V1, V2, F, uE, EMAP, EF, EI, 0.001)
```
where `0.001` is a chosen bending weight. Furthermore, we can evaluate its gradient and Hessian:
```matlab
g = shell_deformed_gradient(V1, V2, F, uE, EMAP, EF, EI, 0.001);
H = shell_deformed_hessian(V1, V2, F, uE, EMAP, EF, EI, 0.001);
```
Note that the gradient is a dense matrix of shape `#V x 3` while the Hessian is a sparse matrix with shape `3*#V x 3*#V` for which the coordinates `V2` are assumed to be serialized in column-major fashion, i.e. via `V2(:)`.

## Using in Python
To use the Python interface, you first need to either add the path were the resulting the library is stored to your `PYTHONPATH` environment variable, or copy the library to the path of your proeject/where you are executing python.
Then the following Python scripts demonstrates the 'usage' of all included functions:
```python
import igl
import pyshell 
import numpy as np
import scipy

# load example data, remember to adjust the paths!
v_undef, f = igl.read_triangle_mesh("../../data/finger/finger0.ply")
v_def, _ = igl.read_triangle_mesh("../../data/finger/finger1.ply")  # connectivity of second mesh is not needed, since meshes are assumed to be in dense correspondence
bending_weight = 0.001 # often times a good starting point

# change type of integer of the connecitivy, since (for now) pyshell only supports 32-bit integers
f = f.astype(np.int32)

# edge flaps
E, EMAP, EF, EI = igl.edge_flaps(f)

# Energies
print("Membrane energy:", pyshell.membrane_energy(v_undef, v_def, f))
print("Bending energy:", pyshell.bending_energy(v_undef, v_def, f, E, EMAP, EF, EI))
print("Shell energy:", pyshell.shell_energy(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)) 

# Gradients
grad = pyshell.membrane_deformed_gradient(v_undef, v_def, f)
print("Membrane deformed gradient norm:", np.linalg.norm(grad, 'fro'))

grad = pyshell.membrane_undeformed_gradient(v_undef, v_def, f)
print("Membrane undeformed gradient norm:", np.linalg.norm(grad, 'fro'))

grad = pyshell.bending_deformed_gradient(v_undef, v_def, f, E, EMAP, EF, EI)
print("Bending deformed gradient norm:", np.linalg.norm(grad, 'fro'))

grad = pyshell.bending_undeformed_gradient(v_undef, v_def, f, E, EMAP, EF, EI)
print("Bending undeformed gradient norm:", np.linalg.norm(grad, 'fro'))

grad = pyshell.shell_deformed_gradient(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)
print("Shell deformed gradient norm:", np.linalg.norm(grad, 'fro'))

grad = pyshell.shell_undeformed_gradient(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)
print("Shell undeformed gradient norm:", np.linalg.norm(grad, 'fro'))

# Hessians
hess = pyshell.membrane_deformed_hessian(v_undef, v_def, f)
print("Membrane deformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.membrane_undeformed_hessian(v_undef, v_def, f)
print("Membrane undeformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.membrane_mixed_hessian(v_undef, v_def, f)
print("Membrane mixed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

print("Symmetry of mixed membrane Hessian:", scipy.sparse.linalg.norm(hess.T - pyshell.membrane_mixed_hessian(v_undef, v_def, f, False), 'fro'))


hess = pyshell.bending_deformed_hessian(v_undef, v_def, f, E, EMAP, EF, EI)
print("Bending deformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.bending_undeformed_hessian(v_undef, v_def, f, E, EMAP, EF, EI)
print("Bending undeformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.bending_mixed_hessian(v_undef, v_def, f, E, EMAP, EF, EI)
print("Bending mixed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

print("Symmetry of mixed bending Hessian:", scipy.sparse.linalg.norm(hess.T - pyshell.bending_mixed_hessian(v_undef, v_def, f, E, EMAP, EF, EI, False), 'fro'))

hess = pyshell.shell_deformed_hessian(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)
print("Shell deformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.shell_undeformed_hessian(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)
print("Shell undeformed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

hess = pyshell.shell_mixed_hessian(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight)
print("Shell mixed Hessian norm:", scipy.sparse.linalg.norm(hess, 'fro'))

print("Symmetry of mixed shell Hessian:", scipy.sparse.linalg.norm(hess.T - pyshell.shell_mixed_hessian(v_undef, v_def, f, E, EMAP, EF, EI, bending_weight, False), 'fro'))
```