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

# SCM Repository

[diderot] View of /tests/vis15-bugs/sphere-snapbug.diderot
 [diderot] / tests / vis15-bugs / sphere-snapbug.diderot # View of /tests/vis15-bugs/sphere-snapbug.diderot

Mon Oct 24 11:31:07 2016 UTC (2 years, 4 months ago) by glk
File size: 4524 byte(s)
`does not compile with r4793`
```vec3{} ipos = {
[0.15629026,-0.56885761,-0.47557172],
[-0.70781976,-0.56275266,-0.072711878],
[0.087024987,-0.4896704,0.68623573],
[-0.48833442,0.38590688,-0.40620103],
[0.11184707,-1.0173737,1.712598],
[-0.66177833,0.6881938,-0.21358061],
[0.93325704,0.37870207,-0.85690439],
[-0.82517642,0.30594504,0.090624742],
[0.47916487,-0.2020634,0.44536155],
[0.76371771,0.47714719,0.36225933]
};
real rad = 0.4;
input int pcp ("periodicity of population control (if non-zero)") = 2;

real hhInit = 1;        // initial step size
real newDist = 0.5*rad; // how far away to birth new particles
real stepMax = rad;     // limit on distance to travel per iter
int iter = 0;           // which iteration we're on

/* energy functions */

// phi(0) and phi'(0) are bounded
function real  phi(real r) = (1 - r)^4;
function real phi'(real r) = -4*(1 - r)^3;

/*
// electrostatic potential 1/r, scaled to be C^2 continuous with 0 at r==1
function real  phi(real r) =  (1/r)*(1 - r)^3;
function real phi'(real r) = 3 - 1/(r^2) - 2*r;
*/
/*
// Cotangent-based potential from Meyer et al. SMI'05
function real  phi(real r) = 1/tan(r*π/2) + r*π/2 - π/2;
function real phi'(real r) = (π/2)*(1 - (1/sin(r*π/2))^2);
*/

// energy and force from particle at vec3 x
function real enr(vec3 x) = phi(|x|/rad);
function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule

// returns a non-zero vector perpendicular to given non-zero vector v
function vec3 perp3(vec3 v) {
int c = 0;
if (|v| < |v|) {
c = 1;
}
// not v[c] because tensors can only be indexed by constants
if (|v if 1==c else v| < |v|) {
c = 2;
}
// now c is index of v component with largest absolute value
vec3 ret = ([v - v, -v, v] if (c == 0) else
[-v, v - v, v] if (c == 1) else
[-v, v, v - v]);
return ret;
}

/* The particle is initialized at position pos0, with stepsize hh0 */
strand particle (vec3 pos0, real hh0) {
output vec3 pos = pos0/|pos0|;
real hh = hh0;
vec3 step = [0,0,0]; // step along force
real closest = 0;   // distance to closest neighbor
update {
// compute energy and forces on us from neighbors
real energyLast = 0;
vec3 force = [0,0,0];
int ncount = 0;
foreach (particle P in sphere(rad)) {
energyLast += enr(P.pos - pos);
force += frc(P.pos - pos);
ncount += 1;
}
vec3 norm = normalize(pos); // surface normal for unit circle
if (0 == ncount && (pcp > 0 && 0 == iter % pcp)) {
/* no neighbors, so let's make one */
closest = rad; // hamper convergence
//new particle(pos + newDist*normalize(perp3(norm)), hh);
continue;
}

/* Else we have interacting neighbors */
// projection onto tangent surface
force = (identity - norm⊗norm)•force;
step = hh*force;
if (|step| > stepMax) {  // Limit step size
hh *= stepMax/|step|;
step = hh*force;
}

// take step, re-apply constraint, find new energy
vec3 posLast = pos;
pos = normalize(pos + step);
real energy = 0;
foreach (particle P in sphere(rad)) {
energy += enr(P.pos - pos);
closest = min(closest, |P.pos - pos|);
}
step = pos - posLast;

/* Armijo-Goldstein sufficient decrease condition */
if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
hh *= 0.5; // backtrack
pos = posLast;
} else {
hh *= 1.01; // opportunistically increase step size
if (pcp > 0 && 0 == iter % pcp) {
if (ncount < 5) {
//new particle(pos + newDist*normalize(force), hh);
} else if (ncount > 8) {
if (|fmod(pos*1000, 1)| < 1.0/ncount) {
//die;
}
}
}
}
}
}

global {
/* could test convergence based on movement */
//real mvmt = max { |P.step|/rad | P in particle.all};
/* testing convergence based on distance to closest neighbor */
real meancl = mean { P.closest | P in particle.all};
real varicl = mean { (P.closest - meancl)^2 | P in particle.all};
real covcl = sqrt(varicl) / meancl;
print("(iter ", iter, ") mean hh = ", mean { P.hh | P in particle.all}, "; covcl = ", covcl, "\n");
if (iter > 15) {
stabilize;
}
iter += 1;
}

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