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 2237, Sat Mar 2 02:00:56 2013 UTC revision 2240, Sun Mar 3 00:58:00 2013 UTC
# Line 4  Line 4
4  //  //
6  int numOfParticles = length(initPosns)/2;  int numOfParticles = length(initPosns)/2;
7  input real rr = 0.08;     // actual particle radius  input int iterMax = 300;  // maximum number of steps to run
8  input real RR = 0.15;     // neighbor query radius (MUST be bigger than rr, but should get same results for any RR > rr)  input real rr = 2.0;      // actual particle radius
9  input real hh = 0.001;    // integration step size to control how quickly iterations change position.
10  input int MAX_ITER = 200; // maximum number of steps in the program  real RR = rr+0.0;         // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
11    real hhInit = 10.0;       // initial integration step size; can err too big, will be trimmed down during iterations
12    int iter = 1;             // which iteration we're on
13    real stepMax = 2*π/numOfParticles;  // limit on distance to travel per iter
14
15    // ============================ begin setting up gridding of domain
16  vec2 xDom = [-1,1];  vec2 xDom = [-1,1];
17  vec2 yDom = [-1,1];  vec2 yDom = [-1,1];
18  real xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here  real xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here
# Line 15  Line 20
20  vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)  vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)
21  vec2 qGridDim = [xSamples,ySamples];   // how many grid cells you want in each direction for the uniform grid (required for the query function)  vec2 qGridDim = [xSamples,ySamples];   // how many grid cells you want in each direction for the uniform grid (required for the query function)
22  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)  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)
23  int iter = 1;  // ============================ end setting up gridding of domain
24
25  strand Particle (int ii, real posx, real posy) {  strand Particle (int ii, real posx, real posy) {
26      vec2 pos = normalize([posx,posy]);      vec2 pos = normalize([posx,posy]);
27        real hh = hhInit;
28        vec2 posOld1 = pos;   // remember last TWO positions
29        vec2 posOld2 = pos;
30      output vec2 outPos  = pos;      output vec2 outPos  = pos;
31      int id = ii;      real energy = ii;     // should eventually do convergence test based on variance of energies
32      update {      update {
33
34  /*          // print positions
35          print(atan2(pos[1],pos[0]), " ");          print(atan2(pos[1],pos[0]), " ");
36          if (i == numOfParticles-1) {          if (ii == numOfParticles-1) {
37            print("\n");            print("\n");
38          }          }
*/
39
40          real energy = 0;          // compute energy and forces on us
41            energy = 0;
42          vec2 force = [0,0];          vec2 force = [0,0];
43          foreach (Particle p_j in sphere(RR)) {          foreach (Particle p_j in sphere(RR)) {
44              vec2 r_ij = (pos - p_j.pos)/rr;              vec2 r_ij = (pos - p_j.pos)/rr;
45              vec2 d_ij = normalize(r_ij);              vec2 d_ij = normalize(r_ij);
print("===== ", id, ": with ", p_j.id, "; |r_ij| = ", |r_ij|, "\n");

46              if (|r_ij| < 1) {              if (|r_ij| < 1) {
47                  energy += (1 - |r_ij|)^4;                  energy += (1 - |r_ij|)^4;
48                  force += - (-4*(1 - |r_ij|)^3) * d_ij;                  force += - (-4*(1 - |r_ij|)^3) * d_ij;
49              }              }
50          }          }
51            force /= rr;     // smaller particles make larger forces
52
53          vec2 step = hh*force/rr;          // update position based on force
54            posOld2 = posOld1;   // shuffle saved positions down
55            posOld1 = pos;
56            if (energy > 0.0) {  // we have neighbors
57                tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);
58                force = pten•force;  // project force onto tangent surface
59                vec2 step = hh*force;
60                if (|step| > stepMax) {
61                    // decrease hh by factor by which step was too big
62                    hh *= stepMax/|step|;
63                    // and find smaller step
64                    step = hh*force;
65                }
66                // take step and re-find implicit surface
67          pos = normalize(pos + step);          pos = normalize(pos + step);
68                // oscillating if moved close to where we were two steps ago
69                real osci = |pos - posOld2|/|pos - posOld1|;
70                if (osci < 0.1) {
71                    // altering step size slightly disrupts oscillation
72                    hh *= 0.9;
73                }
74            }
75
76          outPos = pos;          outPos = pos;
77            if (iter >= iterMax) {
if (iter >= MAX_ITER) {
78              stabilize;              stabilize;
79          }          }
80

Legend:
 Removed from v.2237 changed lines Added in v.2240