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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2262 - (view) (download)

1 : lamonts 2131 // unit-circle
2 :     //
3 : glk 2230 // Demo of distributing particles on the unit circle
4 : lamonts 2131 //
5 : lamonts 2262 real{} initPosns = load("-");
6 : glk 2230 int numOfParticles = length(initPosns)/2;
7 : glk 2244 input int iterMax = 500; // maximum number of steps to run
8 : glk 2241 input real rr = 0.2; // actual particle radius
9 : glk 2240
10 :     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 : glk 2230 vec2 xDom = [-1,1];
17 :     vec2 yDom = [-1,1];
18 : lamonts 2236 real xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here
19 :     real ySamples = floor((yDom[1] - yDom[0])/RR);
20 : glk 2230 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)
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)
23 : glk 2240 // ============================ end setting up gridding of domain
24 : lamonts 2131
25 : glk 2230 strand Particle (int ii, real posx, real posy) {
26 :     vec2 pos = normalize([posx,posy]);
27 : glk 2240 real hh = hhInit;
28 :     vec2 posOld1 = pos; // remember last TWO positions
29 :     vec2 posOld2 = pos;
30 :     output vec2 outPos = pos;
31 : glk 2244 real energy = ii; // HEY: can do convergence test based on variance of energies
32 :     vec2 force = [0,0]; // or can test convergence based on sum of |force|
33 : lamonts 2131 update {
34 :    
35 : glk 2240 // print positions
36 : glk 2230 print(atan2(pos[1],pos[0]), " ");
37 : glk 2240 if (ii == numOfParticles-1) {
38 : glk 2245 print("\n");
39 : glk 2230 }
40 : lamonts 2224
41 : glk 2240 // compute energy and forces on us
42 :     energy = 0;
43 : glk 2244 force = [0,0];
44 : glk 2230 foreach (Particle p_j in sphere(RR)) {
45 : glk 2235 vec2 r_ij = (pos - p_j.pos)/rr;
46 : glk 2230 vec2 d_ij = normalize(r_ij);
47 :     if (|r_ij| < 1) {
48 :     energy += (1 - |r_ij|)^4;
49 :     force += - (-4*(1 - |r_ij|)^3) * d_ij;
50 :     }
51 : lamonts 2132 }
52 : glk 2240 force /= rr; // smaller particles make larger forces
53 : lamonts 2224
54 : glk 2240 // update position based on force
55 :     posOld2 = posOld1; // shuffle saved positions down
56 :     posOld1 = pos;
57 :     if (energy > 0.0) { // we have neighbors
58 :     tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);
59 :     force = pten•force; // project force onto tangent surface
60 :     vec2 step = hh*force;
61 :     if (|step| > stepMax) {
62 :     // decrease hh by factor by which step was too big
63 :     hh *= stepMax/|step|;
64 :     // and find smaller step
65 :     step = hh*force;
66 :     }
67 :     // take step and re-find implicit surface
68 :     pos = normalize(pos + step);
69 : glk 2245 real travel = |pos - posOld1| + |posOld1 - posOld2|;
70 :     if (travel > 0) {
71 :     // if we've moved in the past two steps, but we've moved back
72 :     // to where we were two steps ago, we're oscillating ==>
73 :     // okay = 0. Two steps in the same direction ==> okay = 1.
74 :     real okay = |pos - posOld2|/travel;
75 :     // slow down if oscillating, speed up a little if not
76 :     hh *= lerp(0.8, 1.01, 0, okay, 1);
77 : glk 2240 }
78 :     }
79 : lamonts 2224
80 : lamonts 2237 outPos = pos;
81 : glk 2240 if (iter >= iterMax) {
82 : lamonts 2224 stabilize;
83 :     }
84 : glk 2230
85 : lamonts 2131 }
86 :     }
87 :    
88 : lamonts 2224 // can add statements in here to do some global computation
89 :     global{
90 : glk 2230 iter+=1;
91 :     }
92 : lamonts 2224
93 : glk 2230 initially {Particle(ii, initPosns{ii*2}, initPosns{ii*2+1})
94 :     | ii in 0 .. numOfParticles-1 };
95 : lamonts 2224

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