Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

revision 2229, Fri Mar 1 20:01:56 2013 UTC revision 2230, Fri Mar 1 22:56:17 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"); // file has 300 particle positions  real{} initPosns = load("-");
6  input int numOfParticles = 100;  // But I'm only look at 100 for right now.  int numOfParticles = length(initPosns)/2;
7  input real r = 0.08;       // particle radius  input real rr = 0.08;       // particle radius
8  input real d = 0.16;       // desired difference of radius (d < r && d > 0)  input real dd = 0.16;       // desired difference of radius (d < r && d > 0)
9  input real R = 0.15;       // neighbor radius  input real RR = 0.15;       // neighbor radius
10  real MAX_FORCE = 2;  input real hh = 0.001;        // integration step size to control how quickly iterations change position.
11  input real h = 0.09;        // integration step size to control how quickly iterations change position.  input int MAX_ITER = 200;   // maximum number of steps in the program
12  input int MAX_ITER = 10;   // maximum number of steps in the program  vec2 xDom = [-1,1];
13  vec4 qWinDim = [-1,1,-1,1]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)  vec2 yDom = [-1,1];
14  vec2 qGridDim = [12,12];    // how many grid cells you want in each direction for the uniform grid (required for the query function)  // HEY we need floor
15  vec2 qCellDim = [r,r];      // the length in each direction for a cell (required for the query function)  int xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here
16  int steps = 1;  int ySamples = floor((yDom[1] - yDom[0])/RR);
17    vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)
18    vec2 qGridDim = [xSamples,ySamples];   // how many grid cells you want in each direction for the uniform grid (required for the query function)
19    vec2 qCellDim = [(xDom[1] - xDom[0])/xSamples,(yDom[1] - yDom[0])/ySamples];      // the length in each direction for a cell (required for the query function)
20    int iter = 1;
21
22  strand Particle (int i, real posx, real posy) {  strand Particle (int ii, real posx, real posy) {
23      vec2 pos = [posx,posy];      vec2 pos = normalize([posx,posy]);
24      output vec2 outPos  = pos;      output vec2 outPos  = pos;
25      int id = i;      int id = ii;
26      update {      update {
27
28          vec2 energy =[0.0,0.0];  /*
29          vec2 force = [0.0,0.0];          print(atan2(pos[1],pos[0]), " ");
30            if (i == numOfParticles-1) {
31          foreach(Particle neighbor in sphere(R)){            print("\n");
real neighborDist = dist (pos,neighbor.pos);
vec2 unitDirection = normalize(pos-neighbor.pos);

if(neighborDist >= 0 && neighborDist <= 1) {
force += ((-4.0 * ((1 - neighborDist)^3)) * unitDirection);
energy += (((1 - neighborDist)^4) * unitDirection);
}
32          }          }
33          real angle = 0.0;  */
34
35          force *= h;          real energy = 0;
36            vec2 force = [0,0];
37            foreach (Particle p_j in sphere(RR)) {
38                vec2 r_ij = pos - p_j.pos;
39                vec2 d_ij = normalize(r_ij);
40                print("===== ", id, ": with ", p_j.id, "; |r_ij| = ", |r_ij|, "\n");
41
42          real isZero = force • force;              if (|r_ij| < 1) {
43                    energy += (1 - |r_ij|)^4;
44          // limit the force                  force += - (-4*(1 - |r_ij|)^3) * d_ij;
45              if ( isZero != 0 && |force| > MAX_FORCE) {              }
force = normalize(force);
force *= MAX_FORCE;
46               }               }
47
48          pos += force;          vec2 step = hh*force;
49          pos = normalize(pos);          pos = normalize(pos + step);

angle = atan2(pos[1],pos[0]);
print(angle, " ");
50
51          outPos = pos;          outPos = pos;
52
53          if(steps >= MAX_ITER) {          if (iter >= MAX_ITER) {
54              stabilize;              stabilize;
55          }          }
56
# Line 62  Line 59
59
60  // can add statements in here to do some global computation  // can add statements in here to do some global computation
61  global{  global{
62    steps+=1;    iter+=1;
63  }  }
64
65  initially {Particle(i,initPosns{i*2},initPosns{i*2+1}) | i in 0 .. numOfParticles-1 };  initially {Particle(ii, initPosns{ii*2}, initPosns{ii*2+1})
66               | ii in 0 .. numOfParticles-1 };
67

Legend:
 Removed from v.2229 changed lines Added in v.2230