mesh based reaction diffusion

grasshopper csharp only

using System;
using System.Collections;
using System.Collections.Generic;

using Rhino;
using Rhino.Geometry;

using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Types;

using System.IO;
using System.Linq;
using System.Data;
using System.Drawing;
using System.Reflection;
using System.Windows.Forms;
using System.Xml;
using System.Xml.Linq;
using System.Runtime.InteropServices;

using Rhino.DocObjects;
using Rhino.Collections;
using GH_IO;
using GH_IO.Serialization;

/// <summary>
/// This class will be instantiated on demand by the Script component.
/// </summary>
public class Script_Instance : GH_ScriptInstance
{
#region Utility functions
  /// <summary>Print a String to the [Out] Parameter of the Script component.</summary>
  /// <param name="text">String to print.</param>
  private void Print(string text) { /* Implementation hidden. */ }
  /// <summary>Print a formatted String to the [Out] Parameter of the Script component.</summary>
  /// <param name="format">String format.</param>
  /// <param name="args">Formatting parameters.</param>
  private void Print(string format, params object[] args) { /* Implementation hidden. */ }
  /// <summary>Print useful information about an object instance to the [Out] Parameter of the Script component. </summary>
  /// <param name="obj">Object instance to parse.</param>
  private void Reflect(object obj) { /* Implementation hidden. */ }
  /// <summary>Print the signatures of all the overloads of a specific method to the [Out] Parameter of the Script component. </summary>
  /// <param name="obj">Object instance to parse.</param>
  private void Reflect(object obj, string method_name) { /* Implementation hidden. */ }
#endregion

#region Members
  /// <summary>Gets the current Rhino document.</summary>
  private readonly RhinoDoc RhinoDocument;
  /// <summary>Gets the Grasshopper document that owns this script.</summary>
  private readonly GH_Document GrasshopperDocument;
  /// <summary>Gets the Grasshopper script component that owns this script.</summary>
  private readonly IGH_Component Component;
  /// <summary>
  /// Gets the current iteration count. The first call to RunScript() is associated with Iteration==0.
  /// Any subsequent call within the same solution will increment the Iteration count.
  /// </summary>
  private readonly int Iteration;
#endregion

  /// <summary>
  /// This procedure contains the user code. Input parameters are provided as regular arguments,
  /// Output parameters as ref arguments. You don't have to assign output parameters,
  /// they will have a default value.
  /// </summary>
  private void RunScript(Mesh mesh, List<Vector3d> dir, List<double> dirF, List<double> v, List<double> dB, List<double> feed, List<double> kill, List<int> seedA, List<int> seedB, int iter, ref object A, ref object B)
  {
    //Adapted by Bathsheba Grossman, February 2019 to take more general inputs.
    //Written by Vicente Soler and changed by Laurent Delrieu
    //http://www.grasshopper3d.com/forum/topics/reaction-diffusion-on-triangular-mesh
    //29th of August 2015

    if (mesh == null) {
      Print("No mesh");
      return;
    }

    //Validate inputs: supply reasonable defaults for nulls and check # of values.  No range checks.
    int size = mesh.Vertices.Count;

    if (dir.Count == 0 || dirF.Count == 0) {        //direction vector and weight
      dirF = new List<double> {1.0};
      dir = new List<Vector3d> { new Vector3d(1, 0, 0) };
    } else {
      if (dir.Count != 1 && dir.Count != size) {
        Print("dir wrong number of values");
        return;
      }
    }

    if (v.Count == 0)     v = new List<double> {1.0};   //slows diffusion => shrinks scale 0-1
    if (dB.Count == 0)    dB = new List<double> {0.5};  //diffusion of B relative to A, always < 1

    if (feed.Count == 0)  feed = new List<double> {0.055};
    if (kill.Count == 0)  kill = new List<double> {0.062};


    if (seedA.Count == 0) seedA = new List<int> {1};    //default A is a field of 1

    if (seedB.Count == 0) {                             //seed B values with 0 or 1
      Random random = new Random(2);
      for (int i = 0; i < size; i++)
        seedB.Add((random.NextDouble() < 0.05) ? 1 : 0);
    }

    if (!CheckLen<Vector3d>(dir, "dir", size)) return;
    if (!CheckLen<double>(dirF, "dirF", size)) return;
    if (!CheckLen<double>(v, "v", size)) return;
    if (!CheckLen<double>(dB, "dB", size)) return;
    if (!CheckLen<double>(feed, "feed", size)) return;
    if (!CheckLen<double>(kill, "kill", size)) return;
    if (!CheckLen<int>(seedA, "seedA", size)) return;
    if (!CheckLen(seedB, "seedB", size)) return;
    //Done validating

    DateTime t;
    TimeSpan s;

    t = System.DateTime.Now;
    var reaction = new ReactionDiffusion(mesh, dir, dirF, v, dB, feed, kill, seedA, seedB, this);

    s = System.DateTime.Now.Subtract(t);
    Print("initialize: " + s.ToString());

    t = System.DateTime.Now;
    reaction.Run(iter);

    s = System.DateTime.Now.Subtract(t);
    Print("run: " + s.ToString());

    A = reaction.listA;
    B = reaction.listB;
  }

  // <Custom additional code> 

  class ReactionDiffusion
  {
    protected List<Particle> particles = new List<Particle>();
    static protected Script_Instance GHSI;       //this to make Print function accessible in class methods.  so much hate.

    public List<GH_Number> listA {get {return particles.Select(particle => new GH_Number(particle.A)).ToList();}}
    public List<GH_Number> listB {get {return particles.Select(particle => new GH_Number(particle.B)).ToList();}}

    public ReactionDiffusion(Mesh mesh, List<Vector3d> dir, List<double > dirF, List<double> v, List<double > dB, List <double> feed, List<double > kill, List<int> seedA, List<int> seedB, Script_Instance SI)
    {
      int size = mesh.Vertices.Count;
      GHSI = SI;

      for (int i = 0; i < size; i++)
      {
        double a = (seedA.Count == 1 ? seedA[0] : seedA[i]);   //a and b starts
        double b = seedB[i];
        double f = (feed.Count == 1 ? feed[0] : feed[i]);     //f and k
        double k = (kill.Count == 1 ? kill[0] : kill[i]);
        Vector3d d = (dir.Count == 1 ? dir[0] : dir[i]);      //direction & strength
        double dFac = (dirF.Count == 1 ? dirF[0] : dirF[i]);
        double vs = (v.Count == 1 ? v[0] : v[i]);             //speed of reaction relative to time
        double diffB = (dB.Count == 1 ? dB[0] : dB[i]);       //diffusion speed of b

        particles.Add(new Particle(a, b, f, k, d, dFac, vs, diffB, mesh.Vertices[i], this, GHSI));
      }

      System.Threading.Tasks.Parallel.For(0, size, i => particles[i].SetNeighbours(i, mesh));
    }

    public void Run(int iterations)
    {
      while(iterations-- > 0)
      {
        System.Threading.Tasks.Parallel.ForEach(particles, particle => particle.Laplacian());

        System.Threading.Tasks.Parallel.ForEach(particles, particle => particle.ReactionDiffusion());
      }
    }

    protected class Particle
    {
      public double A {get; set;}
      public double B {get; set;}
      static protected Script_Instance GHSI;       //this makes Print function accessible in class methods.  so much hate.

      public double dxA, dxB;
      Point3d point;
      double f, k, v, dB, dirF;
      Vector3d dir;
      List<Particle> neighbours = new List<Particle>();
      public List<double> weights = new List<double>();
      ReactionDiffusion reaction;

      public Particle(double A, double B, double f, double k, Vector3d dir, double dirF, double v, double dB, Point3d point, ReactionDiffusion reacc, Script_Instance SI)
      {
        this.A = A;
        this.B = B;
        this.f = f;
        this.k = k;
        this.dir = dir;
        this.dirF = dirF;
        this.v = v;
        this.dB = dB;
        this.point = point;
        this.reaction = reacc;
        GHSI = SI;
      }

      public void SetNeighbours(int i, Mesh mesh)
      {
        double s,c;
        double weightTotal = 0;

        foreach(int j in mesh.Vertices.GetConnectedVertices(i).Where(x => x != i))
        {
          neighbours.Add(reaction.particles[j]);

          double angle = Vector3d.VectorAngle(reaction.particles[j].point - point, dir);
          //if dirF == 1 weight is 1. > 1, higher weight if parallel with curve.  <1, higher weight if perpendicular to curve.
          s = Math.Sin(angle);
          c = Math.Cos(angle);
          double weight = Math.Sqrt(s * s + dirF * dirF * c * c);
          weights.Add(weight);
          weightTotal += weight;
        }
        for (int j = 0; j < weights.Count; j++)
          weights[j] /= weightTotal;
      }

      public void Laplacian()
      {
        double nA = 0, nB = 0;
        for(int i = 0;i < neighbours.Count;i++)
        {
          nA += neighbours[i].A;
          nB += neighbours[i].B * weights[i];
        }
        nA /= neighbours.Count;

        dxA = nA - A;
        dxB = nB - B;
      }

      public void ReactionDiffusion()
      {
        double AB2 = A * B * B;
        A += 1.0 * dxA * v - AB2 + f * (1.0 - A);
        B += dB * dxB * v + AB2 - (k + f) * B;

        A = (A < 0 ? 0 : (A > 1 ? 1 : A));
        B = (B < 0 ? 0 : (B > 1 ? 1 : B));
      }
    }
  }
  private bool CheckLen<T>(List<T> stuff , string name, int size)
  {
  if (stuff.Count != 1 && stuff.Count != size)
  {
  Print(name + " - wrong number of values, should be 1 or {0}",size);
  return false;
}
  return true;
}


  // </Custom additional code> 
}

grasshopper cuda

using System;
using System.Collections;
using System.Collections.Generic;

using Rhino;
using Rhino.Geometry;

using Grasshopper;
using Grasshopper.Kernel;
using Grasshopper.Kernel.Data;
using Grasshopper.Kernel.Types;

using System.IO;
using System.Linq;
using System.Data;
using System.Drawing;
using System.Reflection;
using System.Windows.Forms;
using System.Xml;
using System.Xml.Linq;
using System.Runtime.InteropServices;

using Rhino.DocObjects;
using Rhino.Collections;
using GH_IO;
using GH_IO.Serialization;

using Alea;
using Alea.CSharp;
using Alea.Parallel;
using Alea.Fody;


/// <summary>
/// This class will be instantiated on demand by the Script component.
/// </summary>
public class Script_Instance : GH_ScriptInstance
{
#region Utility functions
  /// <summary>Print a String to the [Out] Parameter of the Script component.</summary>
  /// <param name="text">String to print.</param>
  private void Print(string text) { /* Implementation hidden. */ }
  /// <summary>Print a formatted String to the [Out] Parameter of the Script component.</summary>
  /// <param name="format">String format.</param>
  /// <param name="args">Formatting parameters.</param>
  private void Print(string format, params object[] args) { /* Implementation hidden. */ }
  /// <summary>Print useful information about an object instance to the [Out] Parameter of the Script component. </summary>
  /// <param name="obj">Object instance to parse.</param>
  private void Reflect(object obj) { /* Implementation hidden. */ }
  /// <summary>Print the signatures of all the overloads of a specific method to the [Out] Parameter of the Script component. </summary>
  /// <param name="obj">Object instance to parse.</param>
  private void Reflect(object obj, string method_name) { /* Implementation hidden. */ }
#endregion

#region Members
  /// <summary>Gets the current Rhino document.</summary>
  private readonly RhinoDoc RhinoDocument;
  /// <summary>Gets the Grasshopper document that owns this script.</summary>
  private readonly GH_Document GrasshopperDocument;
  /// <summary>Gets the Grasshopper script component that owns this script.</summary>
  private readonly IGH_Component Component;
  /// <summary>
  /// Gets the current iteration count. The first call to RunScript() is associated with Iteration==0.
  /// Any subsequent call within the same solution will increment the Iteration count.
  /// </summary>
  private readonly int Iteration;
#endregion

  /// <summary>
  /// This procedure contains the user code. Input parameters are provided as regular arguments,
  /// Output parameters as ref arguments. You don't have to assign output parameters,
  /// they will have a default value.
  /// </summary>
  private void RunScript(Mesh mesh, List<Vector3d> dir, List<double> dirF, List<double> v, List<double> dB, List<double> feed, List<double> kill, List<int> seedA, List<int> seedB, int iter, ref object A, ref object B)
  {
    //Written by Vicente Soler and changed by Laurent Delrieu
    //http://www.grasshopper3d.com/forum/topics/reaction-diffusion-on-triangular-mesh
    //29th of August 2015
    //Adapted by Bathsheba Grossman February 2019 to take more general inputs and use Alea to compute on the GPU.

    if (mesh == null) {
      Print("No mesh");
      return;
    }

    //Validate inputs: supply reasonable defaults for nulls and check # of values.  No range checks.
    int size = mesh.Vertices.Count;

    if (dir.Count == 0 || dirF.Count == 0) {        //direction vector and weight
      dirF = new List<double> {1.0};
      dir = new List<Vector3d> { new Vector3d(1, 0, 0) };
    } else {
      if (!CheckLen<Vector3d>(dir, "dir", size)) return;
    }

    if (v.Count == 0 || v == null)    v = new List<double> {1.0};   //slows diffusion => shrinks scale 0-1
    if (dB.Count == 0 || dB == null)  dB = new List<double> {0.5};  //diffusion of B relative to A, always < 1

    if (feed.Count == 0 || feed == null)  feed = new List<double> {0.055};
    if (kill.Count == 0 || feed == null)  kill = new List<double> {0.062};

    if (seedA.Count == 0 || seedA == null) seedA = new List<int> {1};    //default A seeds all 1

    if (seedB.Count == 0 || seedB == null) {                             //default B seeds 1/20 1 rest 0
      Random random = new Random(2);
      for (int i = 0; i < size; i++)
        seedB.Add((random.NextDouble() < 0.05) ? 1 : 0);
    }

    if (!CheckLen<double>(dirF, "dirF", size)) return;
    if (!CheckLen<double>(v, "v", size)) return;
    if (!CheckLen<double>(dB, "dB", size)) return;
    if (!CheckLen<double>(feed, "feed", size)) return;
    if (!CheckLen<double>(kill, "kill", size)) return;
    if (!CheckLen<int>(seedA, "seedA", size)) return;
    if (!CheckLen(seedB, "seedB", size)) return;
    //Done validating

    var devices = Device.Devices;       //hello Alea, what's my video card?
    Print("Alea: {0} GPU(s) available", devices.Length);
    Print("Alea: {0}", devices[0].ToString());

    DateTime t;
    TimeSpan s;

    t = System.DateTime.Now;
    var reaction = new ReactionDiffusion(mesh, dir, dirF, v, dB, feed, kill, seedA, seedB);

    s = System.DateTime.Now.Subtract(t);
    Print("initialize: " + s.ToString());

    t = System.DateTime.Now;
    reaction.Run(iter);

    s = System.DateTime.Now.Subtract(t);
    Print("run: " + s.ToString());

    A = reaction.listA;
    B = reaction.listB;
  }

  // <Custom additional code> 

  class ReactionDiffusion
  {
    protected List<Particle> particles = new List<Particle>();

    public List<GH_Number> listA {get {return particles.Select(p => new GH_Number(p.A)).ToList();}}
    public List<GH_Number> listB {get {return particles.Select(p => new GH_Number(p.B)).ToList();}}

    public ReactionDiffusion(Mesh mesh, List<Vector3d> dir, List<double > dirF, List<double> v, List<double > dB, List <double> feed, List<double > kill, List<int> seedA, List<int> seedB)
    {
      int size = mesh.Vertices.Count;

      for (int i = 0; i < size; i++)                          //make particles
      {
        double a = (seedA.Count == 1 ? seedA[0] : seedA[i]);  //a and b starts
        double b = seedB[i];
        double f = (feed.Count == 1 ? feed[0] : feed[i]);     //f and k
        double k = (kill.Count == 1 ? kill[0] : kill[i]);
        Vector3d d = (dir.Count == 1 ? dir[0] : dir[i]);      //direction & strength
        double dFac = (dirF.Count == 1 ? dirF[0] : dirF[i]);
        double vs = (v.Count == 1 ? v[0] : v[i]);             //speed of reaction relative to time
        double diffB = (dB.Count == 1 ? dB[0] : dB[i]);       //diffusion speed of b

        particles.Add(new Particle(a, b, f, k, d, dFac, vs, diffB, mesh.Vertices[i], this));

      }
      //find particle neighbors and weights
      System.Threading.Tasks.Parallel.For(0, size, i => particles[i].SetNeighbours(i, mesh));
    }

    public void Run(int iterations)
    {
      /*//how it used to be
      while(iterations-- > 0)
      {
        System.Threading.Tasks.Parallel.ForEach(particles, particle => particle.Laplacian());
        System.Threading.Tasks.Parallel.ForEach(particles, particle => particle.ReactionDiffusion());
      }*/

      int size = particles.Count;

      //make neighbours and weights into non-ragged 2d arrays
      int[] nblen = particles.Select(particle => particle.nbris.Count).ToArray();
      int maxnbr = nblen.Max();

      int[,] nbrs = new int[size, maxnbr];
      double[,] weights = new double[size, maxnbr];

      Particle p;
      for (int i = 0; i < size; i++)
      {
        p = particles[i];
        for (int j = 0; j < nblen[i]; j++) {  //oh c#, no faster way?
          nbrs[i, j] = p.nbris[j];
          weights[i, j] = p.weights[j];
        }
      }

      var gpu = Alea.Gpu.Default;             //put particle values into arrays and copy them to the GPU

      var A = gpu.Allocate<double>(particles.Select(particle => particle.A).ToArray());
      var B = gpu.Allocate<double>(particles.Select(particle => particle.B).ToArray());
      var f = gpu.Allocate<double>(particles.Select(particle => particle.f).ToArray());
      var k = gpu.Allocate<double>(particles.Select(particle => particle.k).ToArray());
      var v = gpu.Allocate<double>(particles.Select(particle => particle.v).ToArray());
      var dB = gpu.Allocate<double>(particles.Select(particle => particle.dB).ToArray());

      var gnblen = gpu.Allocate(nblen);
      int[,] gnbrs = gpu.Allocate(nbrs);      //because these arrays aren't ragged they can be copied en bloc
      double[,] gweights = gpu.Allocate(weights);

      var dxA = gpu.Allocate<double>(new double[size]);    //these start at 0
      var dxB = gpu.Allocate<double>(new double[size]);

      while(iterations-- > 0)                 //do it
      {
        //Laplacians
        gpu.For(0, size, i =>
          {
          double nA = 0, nB = 0;
          for(int n = 0; n < gnblen[i]; n++)
          {
            nA += A[gnbrs[i, n]];
            nB += B[gnbrs[i, n]] * gweights[i, n];
          }
          nA /= gnblen[i];

          dxA[i] = nA - A[i];
          dxB[i] = nB - B[i];
          });

        //reaction diffusion
        gpu.For(0, size, i =>
          {
          double AB2 = A[i] * B[i] * B[i];
          A[i] += 1.0 * dxA[i] * v[i] - AB2 + f[i] * (1.0 - A[i]);
          B[i] += dB[i] * dxB[i] * v[i] + AB2 - (k[i] + f[i]) * B[i];

          A[i] = (A[i] < 0 ? 0 : (A[i] > 1 ? 1 : A[i]));
          B[i] = (B[i] < 0 ? 0 : (B[i] > 1 ? 1 : B[i]));
          });
      }

      double[] resA = new double[size];
      double[] resB = new double[size];
      Gpu.CopyToHost(A).CopyTo(resA, 0);    //pull results out of GPU
      Gpu.CopyToHost(B).CopyTo(resB, 0);

      for (int i = 0; i < size; i++)        //put them back into particles. can't linq do this?
      {
        particles[i].A = resA[i];
        particles[i].B = resB[i];
      }

      Gpu.Free(A);                          //because we're nice people
      Gpu.Free(B);
      Gpu.Free(dxA);
      Gpu.Free(dxB);
      Gpu.Free(f);
      Gpu.Free(k);
      Gpu.Free(v);
      Gpu.Free(dB);
      Gpu.Free(gnblen);
      Gpu.Free(gnbrs);
      Gpu.Free(gweights);
    }

    protected class Particle
    {
      public double A {get; set;}
      public double B {get; set;}

      Point3d point;
      public double f, k, v, dB, dirF;
      Vector3d dir;

      public List<int>nbris = new List<int>();          //neighbour indices
      public List<double> weights = new List<double>();
      public double weightTotal;
      ReactionDiffusion reaction;

      public Particle(double A, double B, double f, double k, Vector3d dir, double dirF, double v, double dB, Point3d point, ReactionDiffusion reacc)
      {
        this.A = A;
        this.B = B;
        this.f = f;
        this.k = k;
        this.dir = dir;
        this.dirF = dirF;
        this.v = v;
        this.dB = dB;
        this.point = point;
        this.reaction = reacc;
      }

      public void SetNeighbours(int i, Mesh mesh)
      {
        double s,c;

        foreach(int j in mesh.Vertices.GetConnectedVertices(i).Where(x => x != i))
        {
          nbris.Add(j);

          double angle = Vector3d.VectorAngle(reaction.particles[j].point - point, dir);
          //if dirF == 1 weight is 1. > 1, higher weight if parallel with curve.  <1, higher weight if perpendicular to curve.
          s = Math.Sin(angle);
          c = Math.Cos(angle);
          double weight = Math.Sqrt(s * s + dirF * dirF * c * c);
          weights.Add(weight);
          weightTotal += weight;
        }
        for (int j = 0; j < weights.Count; j++)
          weights[j] /= weightTotal;
      }

    } //end Particle
  }  //end ReactionDiffusion

  private bool CheckLen<T>(List<T> stuff , string name, int size)
  {
  if (stuff.Count != 1 && stuff.Count != size)
  {
  Print(name + " - wrong number of values, should be 1 or {0}",size);
  return false;
}
    return true;
  }

  // </Custom additional code> 
}