Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

 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 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