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

SCM Repository

[diderot] Diff of /examples/iso2d-spatial/iso2d-glk3.diderot
ViewVC logotype

Diff of /examples/iso2d-spatial/iso2d-glk3.diderot

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3018, Mon Mar 9 21:51:21 2015 UTC revision 3019, Mon Mar 9 21:54:11 2015 UTC
# Line 1  Line 1 
1  input int isoIterMax = 10;  // limit on isocontour search  input real radius = 0.1;    // particle interaction radius
 input int forceIterMax = 0; // HEY get rid of this with global convergence criterion  
 input real radius = 0.2;    // particle interaction radius  
2  input real epsilon = 0.001; // convergence criterion  input real epsilon = 0.001; // convergence criterion
3  input int res = 12;         // initialization sampling resolution  input int res = 13;         // initialization sampling resolution
4  input real isoval = 1;      // isovalue  input real isoval = 1;      // isovalue
5    real motion = ∞;            // mean particle motion in last iter
6    int iter = 0;               // global iteration count
7    input int iterLimit = 0;    // upper limit on # iterations
8    
9  // scalar field in which we sample isocontour at isoval  // scalar field in which we sample isocontour at isoval
10  field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");  field#1(2)[] F = bspln3 ⊛ image("data/hex.nrrd");
11    
12  real motion = ∞;            // mean normalized particle motion  // inter-particle energy, and its derivative
13    function real  phi(real x)  =   (1 - |x|)^4 if |x| < 1 else 0.0;
14  function real phi(real x) = (1 - |x|)^4;  function real phi'(real x) = -4*(1 - |x|)^3 if |x| < 1 else 0.0;
 function real phi'(real x) = -4*(1 - |x|)^3;  
15    
16  strand Particle (int ID, vec2 pos0) {  strand Particle (vec2 pos0) {
17    vec2 pos = pos0;       // particle position    output vec2 pos = pos0;   // particle position
   // GLK asks: why can't we just say "output vec2 pos = pos0;"?  
   // this leads to "Error: pos needs to be a strand state variable"  
   // As is, have to make a separate outPos variable, set via stabilize  
18    vec2 dpos = [0,0];     // change in position    vec2 dpos = [0,0];     // change in position
19    bool foundIso = false; // initial isocontour search done    bool foundIso = false; // initial isocontour search done
   int isoIter = 0;       // iteration of isocontour search  
   int forceIter = 0;  // HEY get rid of this with global convergence criterion  
20    real tt = radius/2;    // line search step size    real tt = radius/2;    // line search step size
21    
   output vec2 outPos = [0,0];  
   stabilize {  
     outPos = pos;  
   }  
   
22    update {    update {
23      if (!foundIso) {      if (!foundIso) {
24        if (!inside(pos, F) || isoIter == isoIterMax) {        if (!inside(pos, F) || iter == 10) {
25          die; // quit if outside field or took too many steps          die; // quit if outside field or took too many steps
26        }        }
27        // Newton-Raphson step        // Newton-Raphson step
28        dpos = -normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;        dpos = -normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
29        pos += dpos;        pos += dpos;
       isoIter += 1;  
30        if (|dpos|/radius < epsilon) {        if (|dpos|/radius < epsilon) {
31          foundIso = true;          foundIso = true;
32        }        }
33      } else {      } else {
34        // if (0 == ID) print(ID, " tt=", tt, ", motion=", motion, "\n");        if (motion < epsilon)
       if (motion < epsilon) {  
35          stabilize;          stabilize;
36        }        if (iter == iterLimit)
       forceIter += 1;  
       if (forceIter == forceIterMax) {  
37          stabilize;          stabilize;
       }  
38        real energy=0;        real energy=0;
39        vec2 force=[0,0];        vec2 force=[0,0];
40        foreach (Particle P in sphere(radius)) {        foreach (Particle P in sphere(radius)) {
# Line 65  Line 51 
51          // take Newton-Raphson step back to surface          // take Newton-Raphson step back to surface
52          pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;          pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
53          pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;          pos -= normalize(∇F(pos)) * (F(pos) - isoval)/|∇F(pos)|;
54          real newenergy = 0;          real energyNew = 0;
         // GLK asks: is there a way to write this summation with the  
         // same syntax as used in the global reductions below? e.g.:  
         // newenergy = sum{phi(|pos - P.pos|/radius) | Particle P in sphere(radius)}  
55          foreach (Particle P in sphere(radius))          foreach (Particle P in sphere(radius))
56            newenergy += phi(|pos - P.pos|/radius);            energyNew += phi(|pos - P.pos|/radius);
57          if (newenergy > energy - 0.5*tt*|force|) {          if (energyNew > energy - 0.5*tt*|force|) {
58            tt *= 0.5; // backtrack            tt *= 0.5; // backtrack
59            pos = posLast;            pos = posLast;
60          } else {          } else {
# Line 83  Line 66 
66  }  }
67  global{  global{
68    motion = mean{ |P.dpos|/radius | P in Particle.all };    motion = mean{ |P.dpos|/radius | P in Particle.all };
69      print(iter, ": motion=", motion, "\n");
70      iter+=1;
71  }  }
72    
73  initially { Particle(ui + res*vi,  initially { Particle([lerp(-1.5, 1.5, 0, ui, res-1),
                      [lerp(-1.5, 1.5, 0, ui, res-1),  
74                        lerp(-1.5, 1.5, 0, vi, res-1)])                        lerp(-1.5, 1.5, 0, vi, res-1)])
75               | vi in 0..(res-1), ui in 0..(res-1) };               | vi in 0..(res-1), ui in 0..(res-1) };

Legend:
Removed from v.3018  
changed lines
  Added in v.3019

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