Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] View of /examples/iso2d-spatial/iso2d-glk2.diderot
ViewVC logotype

View of /examples/iso2d-spatial/iso2d-glk2.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2993 - (download) (annotate)
Sat Mar 7 13:42:29 2015 UTC (4 years, 5 months ago) by glk
File size: 2787 byte(s)
mostly white-space change
int gridSize = 10;
real isoval = 1;
input real rad = 0.2;      // actual particle radius
real stepMax = rad/2;
real evariance = ∞;
real hhInit = 10;  // initial integration step size; can err too big, will be trimmed down during iterations

input int newtonIterMax = 10;
input int forceIterMax = 20;

field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");

strand Particle (int ID, vec2 pos0) {
  vec2 pos = pos0;
  output vec2 outPos = [0,0];

  bool foundContour = false;
  int iter = 0;
  int forceIter = 0;
  real energy = 0;    // has to be a member variable so that global reduction can see it
  vec2 force = [0,0];
  real hh = hhInit;

  stabilize {
    outPos = pos;
  }

  update {
    if (!foundContour) {
      if (!inside(pos, F) || iter == newtonIterMax) {
        die; // quit if outside field or took too many steps
      }
      // Newton-Raphson step
      vec2 step = -normalize(∇F(pos))*(F(pos) - isoval)/|∇F(pos)|;
      pos += step;
      iter += 1;
      if (|step| < 0.001) {
        foundContour = true;
      }
    } else {
      if (forceIter == forceIterMax) {
        stabilize;
      }
      energy = 0; force = [0,0];
      foreach (Particle p_j in sphere(rad)) {
        vec2 r_ij = (pos - p_j.pos)/rad;
        energy += (1 - |r_ij|)^4;
        force -= normalize(r_ij)*(-4*(1 - |r_ij|)^3)/rad;
      }
      // update position based on force
      force -= normalize(∇F(pos))⊗normalize(∇F(pos))•force;
      vec2 step = hh*force;
      if (|step| > stepMax) {
        // decrease hh by factor by which step was too big
        hh *= stepMax/|step|;
        // and find smaller step
        step = hh*force;
      }
      // take step and re-find implicit surface
      vec2 posLast = pos;
      pos += step;
      pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));  // Newton-Raphson step
      pos += -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));  // Newton-Raphson step
      real newenergy = 0;
      foreach (Particle p_j2 in sphere(rad)) {
        vec2 r_ij = (pos - p_j2.pos)/rad;
        newenergy += (1 - |r_ij|)^4;
      }
      if (newenergy > energy) {
        // bad step; try next time with smaller step
        hh *= 0.2;
        pos = posLast;
      } else {   // this reduces oscillication, but we still converge
        hh *= 1.2;
      }
      forceIter +=1;
    } // if (!foundContour)
  }
}
global{
  real energyMean = mean{P.energy | P in Particle.all};
  evariance = mean{ (P.energy - energyMean) * (P.energy - energyMean) | P in Particle.all};
}

initially { Particle(ui + gridSize*vi,
                     [lerp(-1.5, 1.5, -0.5, real(ui), real(gridSize)-0.5),
                      lerp(-1.5, 1.5, -0.5, real(vi), real(gridSize)-0.5)])
             | vi in 0..(gridSize-1), ui in 0..(gridSize-1) };

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0