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 2985 - (download) (annotate)
Sat Mar 7 00:22:32 2015 UTC (4 years, 2 months ago) by glk
File size: 4185 byte(s)
more comments, another bug
int gridSize = 10;
real isoval = 1;
field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");
input int newtonIterMax = 50;
input real stepScale = 1.0;
real epsilon = 0.0001;
//////////////// Global Variables for Spacing Particles //////////////////////
input real rr = 0.4;      // actual particle radius
real stepMax = rr/2;
real evariance = ∞;
real hhInit = 10.0;  // initial integration step size; can err too big, will be trimmed down during iterations
int iterSpacing = 1;             // which iteration we're on
input int forceIterMax = 20;
/////////////////////////////////////////////////////////////////////////////

strand Particle (int ID, vec2 pos0) {
    vec2 pos = [pos0[0] - pos0[1]/gridSize, pos0[1] + pos0[0]/gridSize];
    output vec2 outPos = [0,0];

    bool foundContour = false;
    int forceIter = 1;
    int newtonIter = 0;
    real energy = 0;     // HEY: can do convergence test based on variance of energies
    vec2 force = [0,0];   // or can test convergence based on sum of |force|

    vec2 posOld = pos;
    real hh = hhInit;

    stabilize {
      outPos = pos;
    }

    update {
       if (!foundContour) {
            // We bail if we're no longer inside or taken too many steps.
            if (!inside(pos, F) || newtonIter > newtonIterMax) {
    	         die;
            }
            if (|∇F(pos)| == 0.0) {  // can't compute step if |∇F|, so have to bail
    	         die;
            }
            vec2 delta = -((F(pos) - isoval)/|∇F(pos)|)*normalize(∇F(pos));  // Newton-Raphson step
            if (|delta| < epsilon) {    // we've converged if step is small enough
                foundContour = true;
            } else {
                pos += delta;
                newtonIter += 1;
            }
       } else { // if (!foundContour)
          if (forceIter > forceIterMax) {
            stabilize;
          }
          energy = 0;
          force = [0,0];
          foreach (Particle p_j in sphere(rr)) {
            if (p_j.foundContour) {
              vec2 r_ij = (pos - p_j.pos)/rr;
              if (|r_ij| < 1) {
                energy += (1 - |r_ij|)^4;
                force += - (-4*(1 - |r_ij|)^3) * normalize(r_ij);
              }
            }
          }
          force /= rr;     // smaller particles make larger forces
          // update position based on force
          posOld = pos;
          if (energy > 0.0) {  // we have neighbors
              tensor[2,2] pten = identity[2] - normalize(∇F(pos))⊗normalize(∇F(pos));
              force = pten•force;  // project force onto tangent surface
              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
              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_j in sphere(rr)) {
                if (p_j.foundContour) {
                  vec2 r_ij = (pos - p_j.pos)/rr;
                  if (|r_ij| < 1) {
                    newenergy += (1 - |r_ij|)^4;
                  }
                }
              }
              if (newenergy > energy) {
                // bad step: energy went up;
                // try again next time with smaller step
                hh *= 0.5;
                pos = posOld;
              }
          } // if (energy > 0.0)
          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