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

SCM Repository

[diderot] View of /tests/examples/sphere/sphere.diderot
ViewVC logotype

View of /tests/examples/sphere/sphere.diderot

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4675 - (download) (annotate)
Sat Oct 1 12:36:25 2016 UTC (2 years, 11 months ago) by glk
File size: 4920 byte(s)
another particle example, starting with minor modification of circle
/* ==========================================
Mutually-repelling particles on a unit sphere

This is based on the unit circle example; see that program for more
detailed and explanatory comments.
========================================== */

input vec3{} ipos ("initial positions for all particles") = load("vec3.nrrd");
input real rad ("radius of particle's energy support") = 0.1;
input real eps ("system convergence threshold, computed as the coefficient-of-variation of distances to nearest neighbors") = 0.05;

real hhInit = 1;        // initial step size
real stepMax = rad;     // limit on distance to travel per iter
int iter = 0;           // which iteration we're on

/* phi defines the potential energy around a particle, as a function of
distance from particle r > 0. By design, the potential function smoothly
falls to 0 at r == 1, so that the potential has compact support, which
simplifies computing the inter-particle interactions. phi'(r) is the
derivative of phi(r). The ' in phi' is merely suggestive; ' is not acting
as an operator. Below are some different possibilities for phi and phi';
only one should be uncommented.

In a conservative force field, force is the *negative* gradient of some
potential. phi(r) is potential at r > 0, due to a particle at 0. It is also
the potential at 0, due to a particle at r. Moving closer to the particle
at r by a small dr > 0, the change in potential is dphi = -phi'(r)*dr,
which is positive for a monotonically decreasing phi(r). The force, or
negative gradient of the potential, at 0 is then phi'(r), which justifies
how no negation is used in the frc() definition below. */

// 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);
*/

/* The enr and frc functions use phi and phi' to define the potential
   due to, and force from, a particle at 2-D offset x with a potential
   field extending to radius rad. */
function real enr(vec3 x) = phi(|x|/rad);
function vec3 frc(vec3 x) = phi'(|x|/rad) * (1/rad) * x/|x|; // chain rule

/* The particle is initialized at position pos0, with stepsize hh0 */
strand particle (vec3 pos0, real hh0) {
   /* "pos" is a special variable name that must be used to enable
      spatial queries of neighboring particles via sphere() */
   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;
      }
      if (0 == ncount) {
         /* no neighbors; do nothing to do but set closest to big value */
         closest = rad;
         continue;
      }

      /* Else we have interacting neighbors */
      vec3 norm = normalize(pos); // surface normal for unit circle
      // projection onto tangent surface
      force = (identity[3] - norm⊗norm)•force;

      // Limit step size
      step = hh*force;
      if (|step| > stepMax) {
         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;

      /* The Armijo-Goldstein sufficient decrease condition */
      if (energy - energyLast > 0.5*(pos - posLast)•(-force)) {
         hh *= 0.5; // backtrack
         pos = posLast;
      } else {
         hh *= 1.02; // opportunistically increase step size
      }
   }
}

global {
   //print("(iter ", iter, ") hello from global\n");
   /* 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("   mean hh = ", mean { P.hh | P in particle.all}, "; covcl = ", covcl, "\n");
   if (covcl < eps) {
      print("Stabilizing with COV(closest) ", covcl, " < ", eps, " (iter ", iter, ")\n");
      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