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

# SCM Repository

[diderot] Annotation of /tests/vis15-bugs/circle.diderot
 [diderot] / tests / vis15-bugs / circle.diderot # Annotation of /tests/vis15-bugs/circle.diderot

Revision 4729 - (view) (download)

 1 : glk 4729 /* 2 : running this under valgrind produces: 3 : Syscall param write(buf) points to uninitialised byte(s)" that seems to be about the memory in the nrrd not being set 4 : due to saving out the first (post-initialization, before first iteration) snapshot 5 : 6 : */ 7 : 8 : input vec2{} ipos ("initial positions for all particles") = load("vec2.nrrd"); 9 : input real rad ("radius of particle's potential energy support") = 0.25; 10 : input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05; 11 : 12 : real hhInit = 1; // initial step size 13 : real stepMax = rad; // limit on distance to travel per iter 14 : int iter = 0; // which iteration we're on 15 : 16 : // phi(0) and phi'(0) are bounded 17 : function real phi(real r) = (1 - r)^4; 18 : function real phi'(real r) = -4*(1 - r)^3; 19 : 20 : function real enr(vec2 x) = phi(|x|/rad); 21 : function vec2 frc(vec2 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule 22 : 23 : /* The particle is initialized at position pos0 */ 24 : strand particle (vec2 pos0) { 25 : /* These strand variables are visible to the global update and to other 26 : strands making spatial queries. Any variable inside the scope of 27 : update{} will not be visible in this way. NOTE: "pos" is a special 28 : variable name that *must* be used to enable spatial queries of 29 : neighboring particles via sphere(). */ 30 : output vec2 pos = pos0/|pos0|; 31 : real hh = hhInit; 32 : vec2 step = [0,0]; // step along force 33 : real closest = 0; // distance to closest neighbor 34 : int ncount = 0; // how many neighbors did we have 35 : 36 : update { 37 : // Compute energy and forces on us from neighbors 38 : real energyLast = 0; 39 : vec2 force = [0,0]; 40 : ncount = 0; 41 : foreach (particle P in sphere(rad)) { 42 : energyLast += enr(P.pos - pos); 43 : force += frc(P.pos - pos); 44 : ncount += 1; 45 : } 46 : if (0 == ncount) { 47 : /* no neighbors; do nothing to do but set closest to big value */ 48 : closest = rad; 49 : continue; 50 : } 51 : 52 : // Else we have interacting neighbors 53 : vec2 norm = normalize(pos); // surface normal for unit circle 54 : // project force onto tangent surface 55 : force = (identity - norm⊗norm)•force; 56 : 57 : /* Limiting the step size (even before testing for sufficient decrease, 58 : below) helps keep particles near the surface they are supposed to be 59 : sampling; this precaution matters more with a non-trivial 60 : (data-driven) constraint surface */ 61 : step = hh*force; // compute step along 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 (of length stepMax) 66 : step = hh*force; 67 : } 68 : 69 : // take step and re-apply constraint 70 : vec2 posLast = pos; 71 : pos = normalize(pos + step); 72 : // find energy at new location, and distance to closest neighbor 73 : real energy = 0; 74 : closest = rad; 75 : foreach (particle P in sphere(rad)) { 76 : energy += enr(P.pos - pos); 77 : closest = min(closest, |P.pos - pos|); 78 : } 79 : // save actual step taken 80 : step = pos - posLast; 81 : 82 : if (energy - energyLast > 0.5*(pos - posLast)•(-force)) { 83 : hh *= 0.5; // backtrack 84 : pos = posLast; 85 : } else { 86 : hh *= 1.02; // opportunistically increase step size 87 : } 88 : } 89 : } 90 : 91 : global { 92 : real meancl = mean { P.closest | P in particle.all}; 93 : real varicl = mean { (P.closest - meancl)^2 | P in particle.all}; 94 : real covcl = sqrt(varicl) / meancl; 95 : print("(iter ", iter, ") mean(hh)=", mean { P.hh | P in particle.all}, 96 : "; COV(closest) = ", covcl, "\n"); 97 : if (covcl < eps) { 98 : print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n"); 99 : stabilize; 100 : } 101 : iter += 1; 102 : } 103 : 104 : initially { particle(ipos[ii]) | ii in 0 .. length(ipos)-1 };

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