Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

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

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

Sun Mar 3 01:02:44 2013 UTC (7 years, 1 month ago) by glk
File size: 3509 byte(s)
`a sane default radius`
```// unit-circle
//
// Demo of distributing particles on the unit circle
//
int numOfParticles = length(initPosns)/2;
input int iterMax = 300;  // maximum number of steps to run
input real rr = 0.2;      // actual particle radius

real RR = rr+0.0;         // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr)
real hhInit = 10.0;       // initial integration step size; can err too big, will be trimmed down during iterations
int iter = 1;             // which iteration we're on
real stepMax = 2*π/numOfParticles;  // limit on distance to travel per iter

// ============================ begin setting up gridding of domain
vec2 xDom = [-1,1];
vec2 yDom = [-1,1];
real xSamples = floor((xDom[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here
real ySamples = floor((yDom[1] - yDom[0])/RR);
vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // i.e. [XMIN, XMAX, YMIN, YMAX]  (required for the query function)
vec2 qGridDim = [xSamples,ySamples];   // how many grid cells you want in each direction for the uniform grid (required for the query function)
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)
// ============================ end setting up gridding of domain

strand Particle (int ii, real posx, real posy) {
vec2 pos = normalize([posx,posy]);
real hh = hhInit;
vec2 posOld1 = pos;   // remember last TWO positions
vec2 posOld2 = pos;
output vec2 outPos = pos;
real energy = ii;     // should eventually do convergence test based on variance of energies
update {

// print positions
print(atan2(pos[1],pos[0]), " ");
if (ii == numOfParticles-1) {
print("\n");
}

// compute energy and forces on us
energy = 0;
vec2 force = [0,0];
foreach (Particle p_j in sphere(RR)) {
vec2 r_ij = (pos - p_j.pos)/rr;
vec2 d_ij = normalize(r_ij);
if (|r_ij| < 1) {
energy += (1 - |r_ij|)^4;
force += - (-4*(1 - |r_ij|)^3) * d_ij;
}
}
force /= rr;     // smaller particles make larger forces

// update position based on force
posOld2 = posOld1;   // shuffle saved positions down
posOld1 = pos;
if (energy > 0.0) {  // we have neighbors
tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos);
force = pten•force;  // project force onto tangent surface
vec2 step = hh*force;
if (|step| > stepMax) {
// decrease hh by factor by which step was too big
hh *= stepMax/|step|;
// and find smaller step
step = hh*force;
}
// take step and re-find implicit surface
pos = normalize(pos + step);
// oscillating if moved close to where we were two steps ago
real osci = |pos - posOld2|/|pos - posOld1|;
if (osci < 0.1) {
// altering step size slightly disrupts oscillation
hh *= 0.9;
}
}

outPos = pos;
if (iter >= iterMax) {
stabilize;
}

}
}

// can add statements in here to do some global computation
global{
iter+=1;
}

initially {Particle(ii, initPosns{ii*2}, initPosns{ii*2+1})
| ii in 0 .. numOfParticles-1 };

```