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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3282 - (download) (annotate)
Tue Oct 13 19:46:34 2015 UTC (4 years, 3 months ago) by lamonts
File size: 4864 byte(s)
Fixed bug with dead strands appearing in query lists
// unit-circle
//
// Demo of distributing particles on the unit circle
//
real{} initPosns = load("posns.nrrd");
int numOfParticles = length(initPosns)/2;
input int iterMax = 500;  // maximum number of steps to run 
input real rr = 0.2;      // actual particle radius
real d = rr + 0.5;   // desired distance for a new particle. 
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
real stdev = ∞;                     // used to test for convergence 
real oldStdev = ∞; 
real threshold = 0.24;       // covnergence threshold  
input real cellSize; 
int MAX_PARTICLES = 400; 
input real eps = 0.01;
bool newNeigbor = true; 
// ============================ 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);

real xSamples = floor((xDom[1] - xDom[0])/cellSize); // Lamont verify logic floor vs ceil here
real ySamples = floor((yDom[1] - yDom[0])/cellSize);

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 = [cellSize,cellSize];      // 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;
    int id = ii;
    bool hasNeighbors = false;  
    output vec2 outPos = pos;
    real energy = ii;     // HEY: can do convergence test based on variance of energies
    vec2 force = [0,0];   // or can test convergence based on sum of |force|
    update {
		hasNeighbors = false; 							
        if(!newNeigbor)  {
            stabilize;
        }
        
		int samePosCount = 0; 
        // compute energy and forces on us
        energy = 0;
        force = [0,0];
        foreach (Particle p_j in sphere(RR)) {
            real isZero = |[pos[0],pos[1],0] × [p_j.pos[0],p_j.pos[1],0]|; 
 			if(isZero > 0) {
            	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;
            	}
            }else { 
              samePosCount+=1; 
            }
        }
                
        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);
            real travel = |pos - posOld1| + |posOld1 - posOld2|;
            if (travel > 0) {
                // if we've moved in the past two steps, but we've moved back
                // to where we were two steps ago, we're oscillating ==>
                // okay = 0. Two steps in the same direction ==> okay = 1.
                real okay = |pos - posOld2|/travel;
                // slow down if oscillating, speed up a little if not
                hh *= lerp(0.8, 1.01, 0, okay, 1);
            }
        }else {  // we do have neighbors so create new neighbors
           hasNeighbors = true; 
           if(samePosCount > 0) { 
        		pos = normalize(pos + [d,d]); 
        	}else { 
           	   new Particle(id + 1,pos[0] + d, pos[1] + d); 
           	}
        }
        outPos = pos;
    }
}

// can add statements in here to do some global computation
global{
  iter+=1;
//  real var = variance{P.energy | P in Particle.all}; 
  newNeigbor = exists{P.hasNeighbors | P in Particle.all}; 
 // oldStdev = stdev; 
 // stdev = sqrt(var); 
  //print("stdev = ", stdev, "\n"); 
}

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


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