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

SCM Repository

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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4773 - (download) (annotate)
Mon Oct 17 19:08:42 2016 UTC (2 years, 8 months ago) by glk
File size: 4524 byte(s)
doesn't compile
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[0]| < |v[1]|) {
      c = 1;
   }
   // not v[c] because tensors can only be indexed by constants
   if (|v[1] if 1==c else v[0]| < |v[2]|) {
      c = 2;
   }
   // now c is index of v component with largest absolute value
   vec3 ret = ([v[1] - v[2], -v[0], v[0]] if (c == 0) else
               [-v[1], v[0] - v[2], v[1]] if (c == 1) else
               [-v[2], v[2], v[0] - v[1]]);
   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[3] - 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;
      closest = rad;
      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[0]*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