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

# SCM Repository

[diderot] View of /tests/lamont-tests/unit-circle-pc3.diderot
 [diderot] / tests / lamont-tests / unit-circle-pc3.diderot

# View of /tests/lamont-tests/unit-circle-pc3.diderot

Tue Sep 27 20:54:47 2016 UTC (2 years, 9 months ago) by glk
File size: 3855 byte(s)
`initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/`
```#version 2
// unit-circle
//
// Demo of distributing particles on the unit circle
//
input vec2[] initPos = load_sequence("vec2.nrrd");
//input vec2[] initPos = {[0.15606558, -0.56803983]}; // BUG
input int iterMax ("maximum number of iterations to run") = 1;
input real rad ("radius of particle's energy support") = 0.2;
input int pcp ("periodicity of population control") = 5;
input real mvmtEps ("system convergence criterion on particle movement ") = 0.001;

real hhInit = 10;       // initial step size, can be big
real newDist = 0.5*rad; // how far away to birth new particles
real stepMax = 2*rad;   // limit on distance to travel per iter

int iter = 0;           // which iteration we're on

// inter-particle energy, and its derivative
function real  φ(real r) = (1 - r)^4;
function real φ'(real r) = -4*(1 - r)^3;
//function real  φ(real r)  =  (1/r)*(1-r)^3;
//function real φ'(real r) = 3 - 1/(r^2) - 2*r;
function real enr(vec2 x) = φ(|x|/rad);
function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule

strand Particle (int parent_, int ii, vec2 pos0, real hh0) {
int parent = parent_;
vec2 pos = pos0/|pos0|;
real hh = hh0;
output vec2 outPos = pos;
vec2 delta = [0,0];
int id = ii;
int born = 0;
real energyLast = ii;     // HEY: can do convergence test based on variance of energies
vec2 force = [0,0];   // or can test convergence based on sum of |force|
bool first = true;
update {
if (first) {
//print("(iter ", iter, ") hello from ", ii, "\n");
first = false; // HEY remove me
}
// compute energy and forces on us
energyLast = 0;
force = [0,0];
int ncount = 0;
foreach (Particle p_j in sphere(rad)) {
ncount += 1;
vec2 r_ij = p_j.pos - pos;
energyLast += enr(r_ij);
force += frc(r_ij);
}
vec2 norm = normalize(pos); // surface normal (for unit circle)
if (ncount == 0 && (pcp > 0 && 0 == iter % pcp)) {
vec2 npos = pos + newDist*rad*[norm[1],-norm[0]];
new Particle(ii, born, npos, hh);
born += 1;
continue;
}
// else update position based on force from neighbors.
// first, project force onto tangent surface
force = (identity[2] - norm⊗norm)•force;
delta = hh*force;
if (|delta| > stepMax) {
// decrease hh by factor by which step was too big
hh *= stepMax/|delta|;
// and find smaller step
delta = hh*force; // == stepMax*normalize(force);
}
vec2 posLast = pos;
// take step and re-apply constraint
pos = normalize(pos + delta);
real energy = 0;
ncount = 0;
foreach (Particle p_j in sphere(rad)) {
energy += enr(p_j.pos - pos);
ncount += 1;
}
if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
hh *= 0.5; // backtrack
pos = posLast;
} else {
hh *= 1.1;
/* putting new particle creation here avoids possibility of
starting multiple particles at same position while not
moving due to backtracking */
if (ncount < 2 && (pcp > 0 && 0 == iter % pcp)) {
vec2 npos = pos + newDist*rad*normalize(force);
new Particle(ii, born, npos, hh);
born += 1;
}
}
outPos = pos;
}
}

update {
//print("(iter ", iter, ") hello from global\n");
real mvmt = max { |P.delta|/rad | P in Particle.all};
if (mvmt < mvmtEps) {
print("stabilizing with movment ", mvmt, " at iter ", iter, "\n");
stabilize;
}
iter += 1;
if (iter >= iterMax)  {
print("stopping with movment ", mvmt, " at iter ", iter, "\n");
die;
}
}

create_collection {Particle(-1, ii, initPos[ii], hhInit)
| ii in 0 .. length(initPos)-1 }
```

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