Metadata-Version: 2.4
Name: eikonalspiral
Version: 0.1.0
Summary: Eikonal-diffusion solver for periodic propagation
Author-email: Vincent Jacquemet <vincent.jacquemet@umontreal.ca>
License-Expression: MIT
Project-URL: Homepage, https://github.com/jacquemv/eikonalspiral
Classifier: Programming Language :: Python :: 3
Classifier: Operating System :: OS Independent
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Dynamic: license-file

![Illustration of streamlines](https://github.com/jacquemv/eikonalspiral/blob/main/image.png?raw=true)


### Objective

Phase fields may exhibit phase singularities. Our objective is to generate a phase field on a triangulated surface (a two-dimensional manifold embedded in three dimensions, possibly including boundaries and holes) that has phase singularities at predetermined locations with desired chirality (direction of rotation) [1]. The manifold is assumed to be a globally parameterizable surface, which is notably the case when it is homeomorphic to a sphere with holes. 

The approach is based on an initial phase field analytically calculated on the 2D parameter domain, then mapped onto the surface and regularized using the eikonal-diffusion equation, which creates wave fronts spiralling around phase singularities. We developped that approach to automatically create complex initial conditions for reaction-diffusion systems [2]. 

### Simple example

Here is a 2D example of phase map generation on a triangular mesh of a rectangle. The result is the image on the top of the page.
```python
import numpy as np
import matplotlib.pyplot as plt
import eikonalspiral as eik
# create a 500 x 150 rectangular geometry
ver, tri = eik.rectangular_mesh(500, 150, 0.1)
# principal direction of anisotropy
orient = np.full((tri.shape[0], 2), 1/np.sqrt(2))
# distribute 10 phase singularities randomly
location, chirality = eik.random_phase_singularities(ver, tri, n=10)
# generate the phase map
theta = eik.generate_phase_map(ver, tri, None, location, chirality,
                               fiber=orient, wavelength=10, aniso_ratio=4)
# visualize the phase map
theta_interp = np.angle(np.exp(1j*theta[tri]).mean(1))
plt.tripcolor(ver[:, 0], ver[:, 1], triangles=tri, edgecolors='none',
              facecolors=theta_interp, cmap='hsv')
plt.show()
```
The output phase field (theta in radians) is defined on the vertices of the triangular mesh.



### Syntax

```python
theta: np.ndarray = generate_phase_map(
    vertices: np.ndarray,
    triangles: np.ndarray,
    parametrization: np.ndarray|None, 
    location: np.ndarray,
    chirality: np.ndarray, *, 
    fiber: np.ndarray=None,
    aniso_ratio: float=4,
    wavelength: float=16, 
    c: np.ndarray|float=0.1,
    D: np.ndarray|float=0.01, 
    verbose=True
)
```

### Positional arguments

The arguments **vertices** and **triangles** define the triangulated surface (with nv vertices and nt triangles), and **orientation** the vector field on each triangle of that surface; **parametrization** continuously maps the surface onto the plane (two coordinates must be given for each vertex):
- **vertices** (nv-by-3 float array): vertex positions (in cm)
- **triangles** (nt-by-3 int array): vertex indices for each triangle
- **parametrization** (nv-by-2 float array): 2D parameterization of the
    surface (continuous mapping from the vertices in 3D space to the
    2D plane); can be None if the geometry is already in 2D
- **location** (int array): list of vertex indices where spirals should be
    located
- **chirality** (int array): direction of rotation for each spiral (+1 or -1)


### Optional keyword arguments

- **fiber** (nt-by-3 array): fiber orientation in each triangle (default: isotropic conduction).
- **aniso_ratio** (float): conductivity anisotropy ratio (default: 4)
- **wavelength** (float): target median wavelength in cm = conduction velocity * cycle length (default: 16 cm). A shorter wavelength makes the spiral more "curved".
- **c** (float): conduction velocity in cm/rad. For a small circuit around a phase singularity (length L ~ 0.6 cm), c should be set to L/2pi ~ 0.1 cm/rad (which is the default value). 
- **D** (float): diffusion coefficient that is used to stabilize the system. The value should of the order of dx**2 and is set to 0.01 cm2 by default. If the output phase field seems noisy, increase this value.
- **verbose** (bool): if True (default), print iterations.

### Outputs

The function outputs an array of size nv providing the phase in the range (-pi, pi) for each vertex of the mesh.


### Implementation

The eikonal-diffusion solver for periodic propagation [3] uses finite element methods for which implementation details are provided in [4]. It requires a solver for sparse complex linear equation (scipy.sparse.linalg.spsolve is used by default, but it could be replaced by another one if needed).

### Installation

The package can be installed using the command ``pip install eikonalspiral``. The code was tested using Anaconda 2025.12 (python 3.13) on Linux.

### Acknowledgements

This work was supported by the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC grant RGPIN-2020-05252).

### References

1. [Reconstruction of phase maps from the configuration of phase singularities in two-dimensional manifolds](https://hdl.handle.net/1866/32915), A. Herlin, V. Jacquemet, *Phys. Rev. E*, vol. 85, no. 5, pp. 051916, May 2012.

2. [Eikonal-based initiation of fibrillatory activity in thin-walled cardiac propagation models](https://hdl.handle.net/1866/32910), A. Herlin, V. Jacquemet, *Chaos*, vol. 21, pp. 043136, Dec. 2011.

3. [An eikonal approach for the initiation of reentrant cardiac propagation in reaction-diffusion models](https://hdl.handle.net/1866/32912), V. Jacquemet, *IEEE Trans. Biomed. Eng.*, vol. 57, no. 9, pp. 2090-2098, Sep. 2010.

4. [An eikonal-diffusion solver and its application to the interpolation and the simulation of reentrant cardiac activations](https://hdl.handle.net/1866/32901), V. Jacquemet, *Comput. Methods Programs Biomed.*, vol. 108, no. 2, pp. 548-558, Nov. 2012.
