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

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

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