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 2240 - (view) (download)

1 : lamonts 2131 // unit-circle
2 :     //
3 : glk 2230 // Demo of distributing particles on the unit circle
4 : lamonts 2131 //
5 : glk 2230 real{} initPosns = load("-");
6 :     int numOfParticles = length(initPosns)/2;
7 : glk 2240 input int iterMax = 300; // maximum number of steps to run
8 :     input real rr = 2.0; // actual particle radius
9 :    
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 :     real energy = ii; // should eventually do convergence test based on variance of energies
32 : lamonts 2131 update {
33 :    
34 : glk 2240 // print positions
35 : glk 2230 print(atan2(pos[1],pos[0]), " ");
36 : glk 2240 if (ii == numOfParticles-1) {
37 : glk 2230 print("\n");
38 :     }
39 : lamonts 2224
40 : glk 2240 // compute energy and forces on us
41 :     energy = 0;
42 : glk 2230 vec2 force = [0,0];
43 :     foreach (Particle p_j in sphere(RR)) {
44 : glk 2235 vec2 r_ij = (pos - p_j.pos)/rr;
45 : glk 2230 vec2 d_ij = normalize(r_ij);
46 :     if (|r_ij| < 1) {
47 :     energy += (1 - |r_ij|)^4;
48 :     force += - (-4*(1 - |r_ij|)^3) * d_ij;
49 :     }
50 : lamonts 2132 }
51 : glk 2240 force /= rr; // smaller particles make larger forces
52 : lamonts 2224
53 : glk 2240 // 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);
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 : lamonts 2224
76 : lamonts 2237 outPos = pos;
77 : glk 2240 if (iter >= iterMax) {
78 : lamonts 2224 stabilize;
79 :     }
80 : glk 2230
81 : lamonts 2131 }
82 :     }
83 :    
84 : lamonts 2224 // can add statements in here to do some global computation
85 :     global{
86 : glk 2230 iter+=1;
87 :     }
88 : lamonts 2224
89 : glk 2230 initially {Particle(ii, initPosns{ii*2}, initPosns{ii*2+1})
90 :     | ii in 0 .. numOfParticles-1 };
91 : lamonts 2224

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