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

SCM Repository

[diderot] Diff of /branches/lamont/test/implicit-surface/unit-circle.diderot
ViewVC logotype

Diff of /branches/lamont/test/implicit-surface/unit-circle.diderot

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

revision 2223, Thu Feb 28 04:37:47 2013 UTC revision 2224, Thu Feb 28 05:06:21 2013 UTC
# Line 2  Line 2 
2  //  //
3  // Demo of distributing particles on the unit circle  // Demo of distributing particles on the unit circle
4  //  //
5  real{} initPosns = load("../../data/implicit-surface/circlePosns.nrrd");  real{} initPosns = load("../../data/implicit-surface/circlePosns.nrrd"); // file has 300 particle positions
6  int numOfParticles = 100;  int numOfParticles = 100; // But I'm only look at 100 for right now.
7  input real r = 0.2;  // particle radius  input real r = 0.08;  // particle radius
8  input real d = 0.1;  // desired difference of radius (d < r && d > 0)  input real d = 0.16;  // desired difference of radius (d < r && d > 0)
9  input real R = 0.9;  // neighbor radius  input real R = 0.15;  // neighbor radius
10  input real h;  // integration stepsize scaling ∇E update lenght  real MAX_FORCE = 2;
11    input real h = 0.009;  // integration step size to control how quickly iterations change position.
12  input int maxi = 1000; // maximum number of steps in the program  input int maxi = 1000; // maximum number of steps in the program
13    vec4 qWinDim = [-1,1,-1,1]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)
14    vec2 qGridDim = [12,12];  // how many grid cells you want in each direction for the uniform grid (required for the query function)
15    vec2 qCellDim = [r,r];    // the length in each direction for a cell (required for the query function)
16    int steps = 1;
17    
18  strand Particle (real posx, real posy) {  strand Particle (int i, real posx, real posy) {
19      vec2 pos = [posx,posy];      vec2 pos = [posx,posy];
     real isoval = 1;  
     int steps = 0;  
20      output vec2 outPos  = pos;      output vec2 outPos  = pos;
21        int strandIdent = i;
22      update {      update {
23          real energy = 0.0;  
24          real force = 0.0;          vec2 energy =[0.0,0.0];
25            vec2 force = [0.0,0.0];
26    
27          foreach(Particle neighbor in sphere(R)){          foreach(Particle neighbor in sphere(R)){
28              real d = dist (pos,neighbor.pos);              real neighborDist = dist (pos,neighbor.pos);
29              if(d >= 0 && d <= 1) {              vec2 unitDirection = normalize(pos-neighbor.pos);
30                     force += -4.0 * ((1 - d)^3);  
31                     energy += ((1 - d)^4);              if(neighborDist >= 0 && neighborDist <= 1) {
32                       force += ((-4.0 * ((1 - neighborDist)^3)) * unitDirection);
33                       energy += (((1 - neighborDist)^4) * unitDirection);
34               }               }
35          }          }
36          printf("f : ", force, "\n");          real angle = 0.0;
37          pos *= force;  
38          pos = pos / |pos|;          //Looking at how strand #0 changes for 30 iterations.
39            /*if(strandIdent == 0) {
40               print("pos before = ", pos, "\n");
41               print("force = ", force, "\n");
42               angle = atan2(pos[1],pos[0])* 180/π;
43               angle = (angle + 360) if angle < 0 else angle; // atan2 rangs from (π,-π] so need to range from [0,2π)
44               print("Angle = ", angle, "\n");
45            }*/
46    
47            force *= h;
48    
49            real isZero = force • force;
50    
51            // limit the force
52                if ( isZero != 0 && |force| > MAX_FORCE) {
53                     force = normalize(force);
54                 force *= MAX_FORCE;
55                 }
56    
57            pos += force;
58            pos = normalize(pos);
59    
60          /*  if(strandIdent == 0) {
61               print("new pos = ", pos, "\n");
62               angle = atan2(pos[1],pos[0])* 180/π;
63               angle = (angle + 360) if angle < 0 else angle;
64               print("new angle = ", angle, "\n");
65               print("========= iteration done ==================\n");
66            }*/
67          outPos = pos;          outPos = pos;
68          steps+=1;  
69          if(steps >= maxi)          if(steps >= 30) {
70              stabilize;              stabilize;
71      }      }
72  }  }
73  initially {Particle(initPosns{i*2},initPosns{i*2+1}) | i in 0 .. numOfParticles-1 };  }
74    
75    // can add statements in here to do some global computation
76    global{
77      steps+=1;
78    }
79    
80    initially {Particle(i,initPosns{i*2},initPosns{i*2+1}) | i in 0 .. numOfParticles-1 };
81    

Legend:
Removed from v.2223  
changed lines
  Added in v.2224

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