mesh_tissue/numba_functions.py

   1 # ##### BEGIN GPL LICENSE BLOCK #####
   2 #
   3 #  This program is free software; you can redistribute it and/or
   4 #  modify it under the terms of the GNU General Public License
   5 #  as published by the Free Software Foundation; either version 2
   6 #  of the License, or (at your option) any later version.
   7 #
   8 #  This program is distributed in the hope that it will be useful,
   9 #  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 #  GNU General Public License for more details.
  12 #
  13 #  You should have received a copy of the GNU General Public License
  14 #  along with this program; if not, write to the Free Software Foundation,
  15 #  Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  16 #
  17 # ##### END GPL LICENSE BLOCK #####
  18
  19 import numpy as np
  20 try:
  21     from numba import jit
  22
  23     @jit
  24     def numba_reaction_diffusion(n_verts, n_edges, edge_verts, a, b, diff_a, diff_b, f, k, dt, time_steps):
  25         arr = np.arange(n_edges)*2
  26         id0 = edge_verts[arr]     # first vertex indices for each edge
  27         id1 = edge_verts[arr+1]   # second vertex indices for each edge
  28         for i in range(time_steps):
  29             lap_a = np.zeros(n_verts)
  30             lap_b = np.zeros(n_verts)
  31             lap_a0 =  a[id1] -  a[id0]   # laplacian increment for first vertex of each edge
  32             lap_b0 =  b[id1] -  b[id0]   # laplacian increment for first vertex of each edge
  33
  34             for i, j, la0, lb0 in zip(id0,id1,lap_a0,lap_b0):
  35                 lap_a[i] += la0
  36                 lap_b[i] += lb0
  37                 lap_a[j] -= la0
  38                 lap_b[j] -= lb0
  39             ab2 = a*b**2
  40             #a += eval("(diff_a*lap_a - ab2 + f*(1-a))*dt")
  41             #b += eval("(diff_b*lap_b + ab2 - (k+f)*b)*dt")
  42             a += (diff_a*lap_a - ab2 + f*(1-a))*dt
  43             b += (diff_b*lap_b + ab2 - (k+f)*b)*dt
  44         return a, b
  45 except:
  46     pass