/****************************************************************************\
* This file will do calculations of length and area and derivatives of these *
*									     *
* Look at a site s, and consider moving it.  If the sites around it are      *
* s0,s1,... and the voronoi corners are v0=incenter(s,s0,s1),v1=in(s,s1,s2)..*
* Then the lengths of these segments are l1=sqrt(norm^2(v1-v0))) etc and the *
* area of the cell is sum(li*dist(s,si)/4), as computed in voronoi.c         *
*									     *
* As s moves, the corners vi will move.  In fact, dvi = ((vi-s).ds) a(vi,s)  *
* where a(vi,s) is the vector *(s_i+1 - si)/area(s,si,si+1).  Here * is the  *
* Hodge *, which rotates the vector 90deg, and area is the area of the tri., *
* given by (s x si + si x si+i + si+1 x s).  So the length li changes by     *
* dli = (u.a(vi,s) (vi-s) + u.a(vi-1,s) (vi-1 - s)) . ds = grad(li,s) . ds   *
* Here u is the unit vector in the direction vi - vi-1.			     *
*									     *
* We must also consider changes in the length of the segments out from the vi*
* If Li is the length of the segment from vi (between si and si-1), then     *
* dLi = +-length(dvi) = length(a(vi,s)) (s - vi) . ds = grad(li,s) . ds,     *
* since vi moves along the segment whose length we are measuring.	     *
*									     *
*	   \             /						     *
*         L0\    s1.    /L1						     *
*            \   l1    / 						     *
*	s0.  v0-------v1   s2.						     *
*	     /         \						     *
*	  l0/    s.     \l2						     *
*	   /             \						     *
*									     *
*									     *
* Thus we will sum up these derivatives for each edge on the interface.      *
* We also must look at gradients of volumes (areas) of the voronoi cells.    *
* When site s moves, the volume of its cell changes just by                  *
*  dv(s) = sum(dli dist(s,si)),  since the other term you might expect is 0. *
* But the volumes of other cells will change.  Any adjacent cell (cell of    *
* a site si) will change its volume also.  The important terms in 4v(s1) are *
* dist(s1,s)l1 + dist(s0,s1)L0 + dist(s2,s1)L1				     *
* So the derivative is the dot product of ds with			     *
* dst(s1,s)grad(l1,s) + *(v1-v0) +dst(s0,s1)grad(L0,s)+dst(s2,s1)grad(L1,s)  *
* If we sum this over all si, we should get the negative of the grad(v(s),s) *
* since the total volume of space doesn't change.  The second term sums to 0 *
* around the circle (as we noted calculating dv(s)).  The first term looks   *
* like it might cancel dv(s), but actually the signs are the same.  It is the*
* combination of terms 1,3,4 that cancels dv(s).  The sum becomes	     *
* 2 sum (s-vi)/area(s,si+1,si) d(si,si+1) *				     *
* 	( dist(si,si+1) - dist(si,s)cos(th_i) - dist(si+1,s)cos(th_i+1))     *
* which is 0 as required.						     *
* We are interested in the volume gradient for the sum of all solid volumes. *
* If s and all its neighbors are the same type, we can ignore it.  Perhaps   *
* the easiest solution is to sum over all triangles where the three sites    *
* have different types.  Then we add to the Sol_agrad of each of these sites *
* the contribution due to this triangle's corner moving.		     *
*									     *
* Things to consider include pre-computing some parts of this calculation.   *
* Perhaps what I called a(v,s) above should be computed and stored.  Perhaps *
* the vectors (s1-s0) or (v1-v0) should be stored instead of just their      *
* lengths.  Or perhaps what we really need is their dot product, or the cos  *
* of the angle between them.  Someone should work out the formulas again, and*
* see if they have a simpler formulation.				     *
\****************************************************************************/

#include "defines.h"
struct VEC2 V2_times(),V2_sum(),V2_sub(),V2_cross();
double V2_cprod();

/* we assume that the voronoi program has calculated Area(t) for each
triangle and vor_length and del_length for each edge, as well as
Volume(s) for each site */

/* we use the field Type(s) to store 0 for liquid, 1 for solid */

gradients()
{
    SITE_PTR s;
    EDGE_PTR e;

    struct VEC2 Zero;
    Zero.x = Zero.y = 0.;

    for_sites(s)
	Len_grad(s) = Sol_agrad(s) = Zero;
    end_sites(s)

    for_int_edges(e)
	if (Type(orig(e))==Type(dest(e)))
	    continue;
	/* now compute gradient of length of e wrt change in each of the */
	/* four sites determining e */
	orig_lgrad(e);
	orig_lgrad(sym(e));
	right_lgrad(e);
	right_lgrad(sym(e));
    end_int_edges(e)
    /* Now all Len_grad fields should be set correctly; go on to volumes */

    for_tris(e)
	if (Type(orig(e))==Type(dest(e)) && Type(orig(e))==Type(dest(onext(e))))
	    continue;
	tri_agrad(e,orig(lprev(e)));
	tri_agrad(lnext(e),orig(e));
	tri_agrad(lprev(e),dest(e));
    end_tris(e)
}

orig_lgrad(e)	/* gradient of length of e if orig(e) moves */
EDGE_PTR e;	/* pretend e is s->s1 in picture above */
{
    struct VEC2 uvor,aright,aleft,lengrad;

    lengrad = V2_sub(rightv(e),leftv(e)); /* use this for temp storage */
    if (!n_wulff)
	uvor = V2_times(1/info(e).vor_length,lengrad);
	    /* unit vector from v0 to v1 */
    else
    {
	uvor = Wulffv(e);
	    /* the vector we are currently using to compute energy */
	if (V2_dotq(uvor,lengrad)<0)
	    uvor = V2_times(-1.,uvor); /* get orientation right */
    }
    aright = V2_times(1/Area(right(e)),
			V2_cross(V2_sub(destv(e),destv(oprev(e)))));
	/* points out along L1 */
    aleft = V2_times(1/Area(left(e)),
			V2_cross(V2_sub(destv(onext(e)),destv(e))));
	/* points out along L0 */
    lengrad = V2_sub(V2_times(V2_dotq(uvor,aright),V2_sub(rightv(e),origv(e))),
		     V2_times(V2_dotq(uvor,aleft),V2_sub(leftv(e),origv(e))));
    Len_grad(orig(e)) = V2_sum(Len_grad(orig(e)),lengrad);
}

right_lgrad(e)	/* gradient of length of e if right(e) moves */
EDGE_PTR e;	/* pretend e is s0->s1 in picture above */
{
    double factor;
    double fabs();
    struct VEC2 lengrad;
    SITE_PTR s;

    s = dest(rprev(e));
    if (!n_wulff)
	factor = info(e).del_length;
    else /* should be a better way to express this */
    {
	lengrad = V2_sub(destv(e),origv(e)); /* just use this as temp var */
	factor = fabs(V2_cprodq(lengrad, Wulffv(e)));
    }
    factor = factor/Area(right(e)),
    lengrad = V2_times(factor, V2_sub(Site(s),rightv(e)));
    Len_grad(s) = V2_sum(Len_grad(s),lengrad);
}

tri_agrad(e,s)
EDGE_PTR e;	/* pretend e is s2->s1 in picture above */
SITE_PTR s;	/* then s is s in picture above (always called this way) */
{
    struct VEC2 avs;

/* if v1 alone moves, along L1, (due to motion of s) how do volumes change? */
/* we don't use formulas above, but just draw a simple picture, where dv1   */
/* creates two thin triangles with base v2-v1 and v1-v0, which move from    */
/* the cells of sites s2 and s1 to the cell of site s			    */
/* Their areas are  (dv1 x (v2-v1)) and (dv1 x (v1 - v0))		    */

    avs = V2_times(1/Area(left(e)),
			V2_cross(V2_sub(destv(e),Site(s))));
    Sol_agrad(s) = V2_sum(Sol_agrad(s),V2_times(
	Type(s)-Type(orig(e)) * V2_cprod(avs,V2_sub(leftv(onext(e)),leftv(e))),
	V2_sub(leftv(e),Site(s))));
    Sol_agrad(s) = V2_sum(Sol_agrad(s),V2_times(
	Type(s)-Type(dest(e)) * V2_cprod(avs,V2_sub(leftv(dprev(e)),leftv(e))),
	V2_sub(leftv(e),Site(s))));
}