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

SCM Repository

[diderot] Annotation of /tests/vis15-bugs/sphere-snapbug2.diderot
ViewVC logotype

Annotation of /tests/vis15-bugs/sphere-snapbug2.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4773 - (view) (download)

1 : glk 4773 vec3{} ipos = {
2 :     [0.15629026,-0.56885761,-0.47557172],
3 :     [-0.70781976,-0.56275266,-0.072711878],
4 :     [0.087024987,-0.4896704,0.68623573],
5 :     [-0.48833442,0.38590688,-0.40620103],
6 :     [0.11184707,-1.0173737,1.712598],
7 :     [-0.66177833,0.6881938,-0.21358061],
8 :     [0.93325704,0.37870207,-0.85690439],
9 :     [-0.82517642,0.30594504,0.090624742],
10 :     [0.47916487,-0.2020634,0.44536155],
11 :     [0.76371771,0.47714719,0.36225933]
12 :     };
13 :     real rad = 0.4;
14 :     input int pcp ("periodicity of population control (if non-zero)") = 2;
15 :    
16 :     real hhInit = 1; // initial step size
17 :     real newDist = 0.5*rad; // how far away to birth new particles
18 :     real stepMax = rad; // limit on distance to travel per iter
19 :     int iter = 0; // which iteration we're on
20 :    
21 :     /* energy functions */
22 :    
23 :     // phi(0) and phi'(0) are bounded
24 :     function real phi(real r) = (1 - r)^4;
25 :     function real phi'(real r) = -4*(1 - r)^3;
26 :    
27 :     /*
28 :     // electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1
29 :     function real phi(real r) = (1/r)*(1 - r)^3;
30 :     function real phi'(real r) = 3 - 1/(r^2) - 2*r;
31 :     */
32 :     /*
33 :     // Cotangent-based potential from Meyer et al. SMI'05
34 :     function real phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2;
35 :     function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
36 :     */
37 :    
38 :     // energy and force from particle at vec3 x
39 :     function real enr(vec3 x) = phi(|x|/rad);
40 :     function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule
41 :    
42 :     // returns a non-zero vector perpendicular to given non-zero vector v
43 :     function vec3 perp3(vec3 v) {
44 :     int c = 0;
45 :     if (|v[0]| < |v[1]|) {
46 :     c = 1;
47 :     }
48 :     // not v[c] because tensors can only be indexed by constants
49 :     if (|v[1] if 1==c else v[0]| < |v[2]|) {
50 :     c = 2;
51 :     }
52 :     // now c is index of v component with largest absolute value
53 :     vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
54 :     [-v[1], v[0] - v[2], v[1]] if (c == 1) else
55 :     [-v[2], v[2], v[0] - v[1]]);
56 :     return ret;
57 :     }
58 :    
59 :     /* The particle is initialized at position pos0, with stepsize hh0 */
60 :     strand particle (vec3 pos0, real hh0) {
61 :     output vec3 pos = pos0/|pos0|;
62 :     real hh = hh0;
63 :     vec3 step = [0,0,0]; // step along force
64 :     real closest = 0; // distance to closest neighbor
65 :     update {
66 :     // compute energy and forces on us from neighbors
67 :     real energyLast = 0;
68 :     vec3 force = [0,0,0];
69 :     int ncount = 0;
70 :     foreach (particle P in sphere(rad)) {
71 :     energyLast += enr(P.pos - pos);
72 :     force += frc(P.pos - pos);
73 :     ncount += 1;
74 :     }
75 :     vec3 norm = normalize(pos); // surface normal for unit circle
76 :     if (0 == ncount && (pcp > 0 && 0 == iter % pcp)) {
77 :     /* no neighbors, so let's make one */
78 :     closest = rad; // hamper convergence
79 :     //new particle(pos + newDist*normalize(perp3(norm)), hh);
80 :     continue;
81 :     }
82 :    
83 :     /* Else we have interacting neighbors */
84 :     // projection onto tangent surface
85 :     force = (identity[3] - norm⊗norm)•force;
86 :     step = hh*force;
87 :     if (|step| > stepMax) { // Limit step size
88 :     hh *= stepMax/|step|;
89 :     step = hh*force;
90 :     }
91 :    
92 :     // take step, re-apply constraint, find new energy
93 :     vec3 posLast = pos;
94 :     pos = normalize(pos + step);
95 :     real energy = 0;
96 :     closest = rad;
97 :     foreach (particle P in sphere(rad)) {
98 :     energy += enr(P.pos - pos);
99 :     closest = min(closest, |P.pos - pos|);
100 :     }
101 :     step = pos - posLast;
102 :    
103 :     /* Armijo-Goldstein sufficient decrease condition */
104 :     if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
105 :     hh *= 0.5; // backtrack
106 :     pos = posLast;
107 :     } else {
108 :     hh *= 1.01; // opportunistically increase step size
109 :     if (pcp > 0 && 0 == iter % pcp) {
110 :     if (ncount < 5) {
111 :     //new particle(pos + newDist*normalize(force), hh);
112 :     } else if (ncount > 8) {
113 :     if (|fmod(pos[0]*1000, 1)| < 1.0/ncount) {
114 :     //die;
115 :     }
116 :     }
117 :     }
118 :     }
119 :     }
120 :     }
121 :    
122 :     global {
123 :     /* could test convergence based on movement */
124 :     //real mvmt = max { |P.step|/rad | P in particle.all};
125 :     /* testing convergence based on distance to closest neighbor */
126 :     real meancl = mean { P.closest | P in particle.all};
127 :     real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
128 :     real covcl = sqrt(varicl) / meancl;
129 :     print("(iter ", iter, ") mean hh = ", mean { P.hh | P in particle.all}, "; covcl = ", covcl, "\n");
130 :     if (iter > 15) {
131 :     stabilize;
132 :     }
133 :     iter += 1;
134 :     }
135 :    
136 :     initially { particle(ipos[ii], hhInit) | ii in 0 .. length(ipos)-1 };

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