root/trunk/matml/src/ternary/qhull.c

/***************************************
  $Header$

  This file includes the function(s) which interact with qhull to compute the
  convex hull of the free energy function points.
  ***************************************/


#include "freenergy.h"
#include "gibbs.h"
#include <qhull/qset.h>
#include <qhull/qhull.h>
#include <stdlib.h>
#include <math.h>


static int qhull_lock=0; /*+ Lock for qhull memory and structures +*/


/*++++++++++++++++++++++++++++++++++++++
  This determines the facet type from the minimum energy at the midpoint of
  each edge.

  facet_type facettype It returns the facet type.

  coordT *corners Coordinates C2, C3, G of the three facet corners.

  energy_params *eparams Collection of energy parameter structures.

  int nparams Number of energy parameter structures in eparams.
  ++++++++++++++++++++++++++++++++++++++*/

static facet_type facettype (coordT *corners, energy_params *eparams,
                             int nparams, double T, double P)
{
  int i, cparam, ret=0;
  double edgenergy, edgecenter[3];

  edgecenter[0] = (corners[0] + corners[3]) / 2.;
  edgecenter[1] = (corners[1] + corners[4]) / 2.;
  edgecenter[2] = (corners[2] + corners[5]) / 2.;
  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
  cparam=0;
  for (i=1; i<nparams; i++)
    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
      {
        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
        cparam = i;
      }
  if (edgenergy > edgecenter[2])
    ret++;
#ifdef DEBUG
  printf ("   Type: side 0-1 edge %g mid %g, ", edgenergy, edgecenter [2]);
#endif
      
  edgecenter[0] = (corners[0] + corners[6]) / 2.;
  edgecenter[1] = (corners[1] + corners[7]) / 2.;
  edgecenter[2] = (corners[2] + corners[8]) / 2.;
  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
  cparam=0;
  for (i=1; i<nparams; i++)
    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
      {
        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
        cparam = i;
      }
  if (edgenergy > edgecenter[2])
    ret++;
#ifdef DEBUG
  printf ("0-2 edge %g mid %g, ", edgenergy, edgecenter [2]);
#endif
      
  edgecenter[0] = (corners[3] + corners[6]) / 2.;
  edgecenter[1] = (corners[4] + corners[7]) / 2.;
  edgecenter[2] = (corners[5] + corners[8]) / 2.;
  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
  cparam=0;
  for (i=1; i<nparams; i++)
    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
      {
        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
        cparam = i;
      }
  if (edgenergy > edgecenter[2])
    ret++;
#ifdef DEBUG
  printf ("1-2 edge %g mid %g", edgenergy, edgecenter [2]);
#endif

  return (facet_type) ret;
}


/*++++++++++++++++++++++++++++++++++++++
  This initializes qhull and uses it to calculate the convex hull of a set of
  ternary compositions, clipping everything above the energy defined by the
  three corners.  It then builds a set of hull_facet structures, and returns
  those, along with the number of calculated facets.  It frees the qhull
  structures before returning.

  int hullCalculate It returns zero or an error code.

  int dim Number of dimensions (should always be 3 for ternary).

  energy_point *points Point information.

  int numpoints Number of points of which to calculate the convex hull.

  energy_params *eparams Collection of energy parameter structures.

  int nparams Number of energy parameter structures in eparams.

  double T Temperature.

  double P Pressure.

  hull_facet **tfacets Array of returned facets.  This is changed by realloc()
  each time the function is called, so it should be NULL or a valid array.

  int *numfacets Number of facets returned.
++++++++++++++++++++++++++++++++++++++*/

int hullCalculate
(int dim, energy_point *points, int numpoints, energy_params *eparams,
 int nparams, double T, double P, hull_facet **facets, int *numfacets)
{
  int i, corner1=-1, corner2=-1, corner3=-1;
  double G1, G2, G3;
  coordT *qpoints; /*+ Point coordinate storage for qhull +*/
  facetT *facet; /*+ For the FORALLfacets loop +*/

  if (dim != 3)
    { printf ("qhullCalcHull: Non-3-D spaces not yet supported\n"); return -1; }

  /*+ This first copies the points array for qhull, and clips the array at the
    plane defined by the corners. +*/
  if (!(qpoints = malloc (dim*numpoints*sizeof(coordT))))
    { printf ("qhullCalcHull: could not allocate memory for points\n");
      return -1; }

  if (qhull_lock)
    { printf ("Qhull use is currently locked\n"); return 1; }
  qhull_lock=1;

  for (i=0; i<numpoints; i++)
    {
      // Copy compositions into the new array
      qpoints[3*i]   = (coordT) points[i].C2;
      qpoints[3*i+1] = (coordT) points[i].C3;

      // Find the top corners
      if ((points[i].C2 == 0. && points[i].C3 == 0.) &&
          (corner1 == -1 || points[i].G > G1))
        {
          corner1 = i;
          G1 = points [corner1].G;
        }
      if ((points[i].C2 == 0. && points[i].C3 == 1.) &&
          (corner2 == -1 || points[i].G > G2))
        {
          corner2 = i;
          G2 = points [corner2].G;
        }
      if ((points[i].C2 == 1. && points[i].C3 == 0.) &&
          (corner3 == -1 || points[i].G > G2))
        {
          corner3 = i;
          G2 = points [corner3].G;
        }
    }
#ifdef DEBUG
  printf ("Corners: %d, %d, %d\n", corner1, corner2, corner3);
#endif

  // Copy free energy into the new array, clipping it at the plane given by the
  // corners
  for (i=0; i<numpoints; i++)
    qpoints[3*i+2] = (coordT)
      fmin (points[i].G, points[corner1].G +
            (points[corner2].G-points[corner1].G) * points[i].C3 +
            (points[corner3].G-points[corner1].G) * points[i].C2);

  // Initialize qhull and calculate the hull
  qh_init_A (stdin, stdout, stderr, 0, NULL);
  qh_init_B (qpoints, numpoints, dim, False);
  qh_qhull ();
#ifdef DEBUG
  qh_printsummary (stdout);
#endif

  /*+ After qhull runs, this reallocates the array pointed to by facets,
    then fills it, eliminating non-simplical facets along the way (which seem
    to typically be the binaries). +*/
  *numfacets = qh num_facets;
#ifdef DEBUG
  printf ("%d total facets\n", *numfacets);
  printf ("facets: pre-realloc 0x%lx, ", *facets);
#endif
  if (!(*facets = realloc (*facets, *numfacets * sizeof (hull_facet))))
    {
      printf ("hullReturnFacets: could not reallocate memory for vertices\n");
      free (qpoints);
      qh_freeqhull (True);
      qhull_lock=0;
      return -1;
    }
#ifdef DEBUG
  printf ("  post-realloc 0x%lx\n", *facets);
#endif

  i=0;
  FORALLfacets
    {
      int j=0;
      coordT corners [9], projarea;
      vertexT *vertex, **vertexp;

#ifdef DEBUG
      printf ("facet %d:\n", i);
#endif

      FOREACHvertex_(facet->vertices)
        {
          int thepoint = (vertex->point-qh first_point)/3;

          if (j<3)
            {
              (*facets) [i].vertex[j] = thepoint;
              corners [3*j]   = vertex->point[0];
              corners [3*j+1] = vertex->point[1];
              corners [3*j+2] = vertex->point[2];
#ifdef DEBUG
              printf ("  vertex %d, coords %g %g %g\n", thepoint,
                      vertex->point[0], vertex->point[1], vertex->point[2]);
#endif
              j++;
            }
        }
          
      /*+ This calculates the area of the projection of each facet.
        That area is given by the determinant of the matrix formed by the
        vectors making up the edges of the facet divided by
        (dimensionality-1)!.
        +latex+For a ternary system, that's
        +latex+$\frac12[(x_1-x_0)(y_2-y_0) - (x_2-x_0)(y_1-y_0)]$.
        +*/
      projarea =
        (corners[3]-corners[0]) * (corners[7]-corners[1]) -
        (corners[6]-corners[0]) * (corners[4]-corners[1]);
#ifdef DEBUG
      printf ("  xyarea %g\n", projarea);
#endif

      // Remove side facets (zero area) and facet with the three corners
      if (projarea == 0. || fabs (projarea) == 1.)
        (*numfacets)--;
      else
        {
          (*facets) [i].type = facettype (corners, eparams, nparams, T, P);
#ifdef DEBUG
          printf ("  Type %d\n", (*facets) [i].type);
#endif
          i++;
        }
    }

  free (qpoints);
  qh_freeqhull (True);
  qhull_lock=0;

#ifdef DEBUG
  printf ("After trimming loop: %d total facets, i=%d\n", *numfacets, i);
  printf ("facets: pre-realloc 0x%lx, ", *facets);
#endif
  if (!(*facets = realloc (*facets, *numfacets * sizeof (hull_facet))))
    { printf ("hullReturnFacets: could not reallocate memory for vertices\n");
      return -1; }
#ifdef DEBUG
  printf ("  post-realloc 0x%lx\n", *facets);
#endif

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This uses Newton iteration to refine a single point within a triangle.

  coordT *point Composition coordinates of the starting point for refining; the
  result will be returned here.

  coordT *corners Coordinates (including energies) of the corners of this
  facet, used for the facet normal: C21,C31,G1, C22,C32,G2, C23,C33,G3.

  energy_params *eparams Collection of energy parameter structures.

  double T Temperature.

  double P Pressure.

  int globaldims Dimensionality of the global space.

  int localdims Dimensionality over which to do the refinement.

  double newton_tolerance Stop Newton iterations when the vertex motion falls
  below this value.
++++++++++++++++++++++++++++++++++++++*/

static void NewtonRefine
(coordT *point, coordT *corners, energy_params *eparams, double T, double P,
 int globaldims, int localdims, double newton_tolerance)
{
  int i;
  double distance = 1., G, slope2, slope3;
  energy_point current;

  if (globaldims != 3)
    { printf ("qhullCalcHull: Non-ternaries not yet supported\n"); return; }
  current.C2 = point[0];
  current.C3 = point[1];

  G = free_energy (current.C2, current.C3, T, P, eparams);

  /*+ Determine the slope(s) of the facet +*/
  if (localdims==3)
    {
      double x2=corners[3]-corners[0], y2=corners[6]-corners[0],
        x3=corners[4]-corners[1], y3=corners[7]-corners[1],
        xG=corners[5]-corners[2], yG=corners[8]-corners[2];
      double vG=x2*y3-x3*y2;
      slope2 = (xG*y3-x3*yG)/vG;
      slope3 = (x2*yG-xG*y2)/vG;

#ifdef DEBUG
      printf ("      corners: %g %g %g %g %g %g %g %g %g\n      slopes: %g %g\n",
              corners[0], corners[1], corners[2],
              corners[3], corners[4], corners[5],
              corners[6], corners[7], corners[8], slope2, slope3);
#endif
    }
  else
    {
      /*
      // Slope for the binary case
      if (side==0 || side==2) // C3=0 or C2+C3=1 sides: dG/dC2
        slope2 = (corners[5]-corners[2])/(corners[3]-corners[0]);
      else if (side==1) // C2=0 side: dG/dC3
        slope2 = (corners[5]-corners[2])/(corners[4]-corners[1]);
      */
    }

  while (distance >= newton_tolerance * newton_tolerance &&
         current.C2>=0. && current.C3>=0. &&
         current.C2<=1. && current.C3<=1.)
    {
      double dC2, dC3, det;

      /*+ Calculate the energy derivatives +*/
      free_energies (&current, 1, T, P, eparams, 0, WITH_DERIVATIVES);

      /*+ Subtract the facet slopes from the energy derivatives +*/
      current.G2 -= slope2;
      current.G3 -= slope3;

      /*+ Solve the linear system to estimate the vector to the minimum +*/
      det = current.G22*current.G33 - current.G23*current.G23;
      dC2 = (current.G3*current.G23 - current.G2*current.G33) / det;
      dC3 = (current.G2*current.G23 - current.G3*current.G22) / det;

      /*+ Add that vector to the point +*/
      current.C2 += dC2;
      current.C3 += dC3;
      distance = dC2*dC2 + dC3*dC3;
#ifdef DEBUG
      printf ("        Added %g,%g (dist %g) to go to %g,%g\n", dC2,dC3,
              sqrt(distance), current.C2,current.C3);
#endif
    }

  /* Return the current point */
  point[0] = current.C2;
  point[1] = current.C3;
}


/*++++++++++++++++++++++++++++++++++++++
  This refines each of the facets in the hull, using Newton-Raphson iteration
  to refine down the vertices of multi-phase facets, and (optionally) adding a
  new vertex at its centroid if in a one-phase region.

  int hullRefine It returns zero or an error code.

  energy_point **points Pointer to point list, which will be modified to
  accommodate the new points created in the refining process.

  int *numpoints Pointer to number of points of which to calculate the convex
  hull, which will also be modified.

  hull_facet **facets Pointer to the array of facets to be refined, which will
  be modified to accommodate new facets created -- and some removed -- in the
  refining process.

  int *numfacets Pointer to the number of facets comprising the convex hull,
  which will also be modified.

  energy_params *eparams Collection of energy parameter structures.

  int nparams Number of energy parameter structures in eparams.

  double T Temperature.

  double P Pressure.

  double newton_tolerance Stop Newton iterations when the vertex motion falls
  below this value.

  double vertex_tolerance Minimum distance between two vertices.  When refining
  a two-phase facet, two of the vertices should go to the same local minimum;
  choose this parameter to indicate how close to consider them the "same"
  point.

  int one_phase_refine Non-zero if the one-phase regions should be refined.
  ++++++++++++++++++++++++++++++++++++++*/

int hullRefine
(energy_point **points, int *numpoints, hull_facet **facets, int *numfacets,
 energy_params *eparams, int nparams, double T,double P,
 double newton_tolerance, double vertex_tolerance, int one_phase_refine)
{
  energy_point *newpoints;
  int newpointsize=100, numnewpoints=0, i;

  if (!(newpoints = (energy_point *) malloc (100*sizeof(energy_point))))
    return 1;

  /*+ Steps in the algorithm:
    +latex+\begin{enumerate}
    +latex+\item
    +html+ <ol><li>
    Loop over the facets
    +html+ </li>
    +*/

  for (i=0; i<*numfacets; i++)
    {
      double corners [9], centroid [3], Gcenter;
      int j=0, eparam [3] = {0,0,0}, thepoint [3];

      /*+
        +latex+\item
        +html+ <li>
        Copy vertex information into local variables: composition, energy,
        closest energy function.
        +html+ </li>
        +*/
      for (j=0; j<3; j++)
        {
          thepoint [j] = (*facets) [i].vertex[j];

          corners [3*j]   = (*points) [thepoint[j]].C2;
          corners [3*j+1] = (*points) [thepoint[j]].C3;
          corners [3*j+2] = (*points) [thepoint[j]].G;
          eparam [j] = (*points)[thepoint[j]].efunc;
        }

#ifdef DEBUG
      printf ("Facet %d: %g,%g %g,%g %g,%g\n", i,
              corners[0], corners[1], corners[3],
              corners[4], corners[6], corners[7]);
      printf ("  Curves: %d %d %d, points %d %d %d, Area: %g\n",
              eparam[0], eparam[1], eparam[2],
              /*thepoint[0], thepoint[1], thepoint[2],I dare you to uncomment!*/
              fabs((corners[3]-corners[0])*(corners[7]-corners[1]) -
                   (corners[6]-corners[0])*(corners[4]-corners[1])));
#endif

      /*+
        +latex+\item
        +html+ <li>
        If this is a one-phase facet, optionally add the centroid as a new
        vertex.
        +html+ </li>
        +*/
      centroid[0] = (corners[0] + corners[3] + corners[6]) / 3.;
      centroid[1] = (corners[1] + corners[4] + corners[7]) / 3.;
      centroid[2] = (corners[2] + corners[5] + corners[8]) / 3.;

      if ((*facets) [i].type == ONE_PHASE && one_phase_refine)
        {
#ifdef DEBUG
          printf ("  Refining one-phase triangle at centroid\n");
#endif
          if (numnewpoints >= newpointsize)
            newpoints = (energy_point *) realloc
              (newpoints, (newpointsize+=100)*sizeof(energy_point));
          newpoints [numnewpoints].C2 = centroid[0];
          newpoints [numnewpoints].C3 = centroid[1];
          newpoints [numnewpoints].G  = free_energy
            (centroid[0], centroid[1], T, P, eparams+eparam [0]);
          newpoints [numnewpoints].efunc = eparam [0];
          numnewpoints++;
        }

      /*+
        +latex+\item
        +html+ <li>
        Otherwise, use Newton-Raphson iteration to refine each vertex,
        which consists of the remaining steps.
        +html+ </li>
        +*/
      else if ((*facets) [i].type != ONE_PHASE)
        {
          coordT newcorners [9] = { corners[0], corners[1], corners[2],
                                    corners[3], corners[4], corners[5],
                                    corners[6], corners[7], corners[8]};

#ifdef DEBUG
          printf ("  Refining multi-phase triangle\n");
#endif
          /*+
            +latex+\item
            +html+ <li>
            If the vertex is on a side of the phase diagram (one or more
            composition variables is zero or they sum to one), refine generally
            starting a bit closer to the others.
            +html+ </li>
            +*/
          for (j=0; j<3; j++)
            {
              if (corners[3*j] == 0. || corners[3*j+1] == 0. ||
                  corners[3*j] == 1. || corners[3*j+1] == 1. ||
                  corners[3*j] + corners[3*j+1] == 1.)
                {
                  newcorners[3*j] =
                    0.9999999998 * corners[3*j] +
                    0.0000000001 * corners[3*((j+1)%3)] +
                    0.0000000001 * corners[3*((j+2)%3)];
                  newcorners[3*j+1] =
                    0.9999999998 * corners[3*j+1] +
                    0.0000000001 * corners[3*((j+1)%3)+1] +
                    0.0000000001 * corners[3*((j+2)%3)+1];
                  newcorners[3*j+2] =
                    0.9999999998 * corners[3*j+2] +
                    0.0000000001 * corners[3*((j+1)%3)+2] +
                    0.0000000001 * corners[3*((j+2)%3)+2];
                }
            }

          /*+
            +latex+\item
            +html+ <li>
            Refine the vertices of the facet using
            +latex+{\tt NewtonRefine()} (section
            +latex+\ref{func_NewtonRefine_qhull.c} page
            +latex+\pageref{func_NewtonRefine_qhull.c}).
            +html+ <a href="#func-NewtonRefine"><tt>NewtonRefine()</tt></a>.</li>
            +*/
          for (j=0; j<3; j++)
            {
#ifdef DEBUG
              printf ("    Refining point at %g, %g function %d\n",
                      newcorners[3*j], newcorners[3*j+1], eparam[j]);
#endif
              NewtonRefine (newcorners+3*j, corners, eparams+eparam[j],
                            T,P, 3,3, newton_tolerance);
#ifdef DEBUG
              printf ("    -> %g, %g\n", newcorners[3*j],
                      newcorners[3*j+1]);
#endif

              newcorners[3*j+2] = free_energy
                (newcorners[3*j],newcorners[3*j+1], T,P, eparams+eparam[j]);
            }

          /*+
            +latex+\item
            +html+ <li>
            Add new points, leaving out duplicates in same place (closer
            than vertex_tolerance).
            +html+ </li>
            +*/

          // Don't add the point, reallocate the points array instead!
          if (newcorners[0]>=0. && newcorners[1]>=0. &&
              newcorners[0]<=1. && newcorners[1]<=1.)
            {
#ifdef DEBUG
              printf ("    Adding new point %g,%g,%g\n",
                      newcorners[0], newcorners[1], newcorners[2]);
#endif
              if (numnewpoints >= newpointsize)
                newpoints = (energy_point *) realloc
                  (newpoints, (newpointsize+=100)*sizeof(energy_point));
              newpoints [numnewpoints].C2 = newcorners[0];
              newpoints [numnewpoints].C3 = newcorners[1];
              newpoints [numnewpoints].G  = newcorners[2];
              newpoints [numnewpoints].efunc = eparam[0];
              numnewpoints++;
            }

          if ((newcorners[3]-newcorners[0])*(newcorners[3]-newcorners[0]) +
              (newcorners[4]-newcorners[1])*(newcorners[4]-newcorners[1]) +
              (newcorners[5]-newcorners[2])*(newcorners[5]-newcorners[2]) >
              vertex_tolerance * vertex_tolerance &&
              newcorners[3]>=0. && newcorners[4]>=0. &&
              newcorners[3]<=1. && newcorners[4]<=1.)
            {
#ifdef DEBUG
              printf ("    Adding new point %g,%g,%g\n",
                      newcorners[3], newcorners[4], newcorners[5]);
#endif
              if (numnewpoints >= newpointsize)
                newpoints = (energy_point *) realloc
                  (newpoints, (newpointsize+=100)*sizeof(energy_point));
              newpoints [numnewpoints].C2 = newcorners[3];
              newpoints [numnewpoints].C3 = newcorners[4];
              newpoints [numnewpoints].G  = newcorners[5];
              newpoints [numnewpoints].efunc = eparam[1];
              numnewpoints++;
            }

          if ((newcorners[6]-newcorners[0])*(newcorners[6]-newcorners[0]) +
              (newcorners[7]-newcorners[1])*(newcorners[7]-newcorners[1]) +
              (newcorners[8]-newcorners[2])*(newcorners[8]-newcorners[2]) >
              vertex_tolerance * vertex_tolerance &&
              (newcorners[6]-newcorners[3])*(newcorners[6]-newcorners[3]) +
              (newcorners[7]-newcorners[4])*(newcorners[7]-newcorners[4]) +
              (newcorners[8]-newcorners[5])*(newcorners[8]-newcorners[5]) >
              vertex_tolerance * vertex_tolerance &&
              newcorners[6]>=0. && newcorners[7]>=0. &&
              newcorners[6]<=1. && newcorners[7]<=1.)
            {
#ifdef DEBUG
              printf ("    Adding new point %g,%g,%g\n",
                      newcorners[6], newcorners[7], newcorners[8]);
#endif
              if (numnewpoints >= newpointsize)
                newpoints = (energy_point *) realloc
                  (newpoints, (newpointsize+=100)*sizeof(energy_point));
              newpoints [numnewpoints].C2 = newcorners[6];
              newpoints [numnewpoints].C3 = newcorners[7];
              newpoints [numnewpoints].G  = newcorners[8];
              newpoints [numnewpoints].efunc = eparam[2];
              numnewpoints++;
            }
        }
    }
  /*+
    +latex+\end{enumerate}
    +html+ </ol>
    +*/

#ifdef DEBUG
  printf ("Adding %d new points to hull\n", numnewpoints);
#endif

  if (!(*points = (energy_point *) realloc
        (*points, (*numpoints + numnewpoints) * sizeof (energy_point))))
    return 1;
          
  for (i=0; i<numnewpoints; i++)
    {
      pointT newcorners[3] =
        { newpoints[i].C2, newpoints[i].C3, newpoints[i].G };
      realT bestdist;
      boolT isoutside;
      facetT *facet;

      (*points) [(*numpoints)+i] = newpoints [i];
    }
  *numpoints += numnewpoints;
  free(newpoints);

#ifdef DEBUG
  printf ("Re-calculating hull using hullCalculate()\n");
#endif

  /*+ Finally, after generating the new points and adding them to the hull, it
    re-calculates the convex hull using
    +latex+{\tt hullCalculate()} (section
    +latex+\ref{func_hullCalculate_qhull.c} page
    +latex+\pageref{func_hullCalculate_qhull.c}).
    +html+ <a href="#func-hullCalculate"><tt>hullCalculate()</tt></a>.</li>
    +*/
  
  return hullCalculate (3, *points,*numpoints, eparams,nparams, T,P, facets,
                        numfacets);
}


#ifndef MIN
#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#endif

#ifndef MAX
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#endif

// For more than three points, qsort() is just easier!
#define MIN3(a,b,c) (MIN (MIN (a,b), c))
#define MAX3(a,b,c) (MAX (MAX (a,b), c))
#define MID3(a,b,c) (((a) < (b)) ? \
                     (((a) < (c)) ? MIN (b,c) : (a)) : \
                     (((a) < (c)) ? (a) : MAX (b,c)))

/*++++++++++++++++++++++++++++++++++++++
  This distills a set of facets to a list of phase boundaries.

  int hullReturnPhaseBoundaries It returns zero or an error code.

  energy_point *points Array of points in the ternary system.

  int n_points Number of points.

  hull_facet *facets Facets from which to create phase boundaries.

  int numfacets Number of facets.

  energy_params *eparams Collection of energy parameter structures.

  double T Temperature.

  double P Pressure.

  phase_boundary **boundaries Pointer to array of phase_boundary objects.  This
  is changed by realloc() each time the function is called, so it should be
  NULL or a valid array.

  int *n_bounds Number of phase boundaries returned in boundaries.
  ++++++++++++++++++++++++++++++++++++++*/

int hullReturnPhaseBoundaries
(energy_point *points, int n_points, hull_facet *facets, int n_facets,
 energy_params *eparams, double T, double P,
 phase_boundary **boundaries, int *n_bounds)
{
  int f;
  *n_bounds = 0;

#ifdef DEBUG
  printf ("Starting ReturnPhaseBoundaries with %d facets, %d points\n",
          n_facets, n_points);
#endif

  for (f=0; f<n_facets; f++)
    {
      int v0 = facets[f].vertex[0], v1 = facets[f].vertex[1],
        v2 = facets[f].vertex[2];

#ifdef DEBUG
      printf ("Facet %d: type %d, vertices %d ",f, facets [f].type, v0);
      printf ("(%d), %d (%d), %d(%d)\n",
              points[v0].efunc, v1,
              points[v1].efunc, v2, points[v2].efunc);
      fflush (stdout);
#endif

      switch (facets [f].type)
        {
        case ONE_PHASE:
          // Totally uninteresting
          break;

        case THREE_PHASE:
          {
            int b=0, bound=-1,
              phase0 = MIN3 (points [v0].efunc, points [v1].efunc,
                             points [v2].efunc),
              phase1 = MID3 (points [v0].efunc, points [v1].efunc,
                             points [v2].efunc),
              phase2 = MAX3 (points [v0].efunc, points [v1].efunc,
                             points [v2].efunc);

            // Check for a boundary with the same three phases
            for (b=0; b<*n_bounds; b++)
              if ((*boundaries) [b].n_phases == 3)
                if ((*boundaries) [b].compos == 3 &&
                    (*boundaries) [b].type == TIE_SIMPLICES &&
                    (*boundaries) [b].phases[0] == phase0 &&
                    (*boundaries) [b].phases[1] == phase1 &&
                    (*boundaries) [b].phases[2] == phase2)
                  bound = b;

            if (bound == -1)
              {
                // No boundary with these three phases, make a new one
                bound = *n_bounds;
#ifdef DEBUG
                printf ("Making new three-phase boundary %d\n", *n_bounds);
#endif
                if (!(*boundaries = realloc
                      (*boundaries, (*n_bounds+1) * sizeof (phase_boundary))))
                  return 1;

                (*boundaries) [bound].compos = 3;
                (*boundaries) [bound].n_phases = 3;
                if (!((*boundaries) [bound].phases = malloc (3*sizeof (int))))
                  return 1;
                (*boundaries) [bound].phases [0] = phase0;
                (*boundaries) [bound].phases [1] = phase1;
                (*boundaries) [bound].phases [2] = phase2;
                (*boundaries) [bound].type = TIE_SIMPLICES;
                (*boundaries) [bound].n_edges = 1;
                if (!((*boundaries) [bound].edges = malloc (3*sizeof (int))))
                  return 1;
                (*boundaries) [bound].edges [0] = MIN3 (v0, v1, v2);
                (*boundaries) [bound].edges [1] = MID3 (v0, v1, v2);
                (*boundaries) [bound].edges [2] = MAX3 (v0, v1, v2);
                (*n_bounds) ++;
#ifdef DEBUG
                printf ("  Corners for bound %d edge 0: %d, %d, %d\n", bound,
                        (*boundaries) [bound].edges [0],
                        (*boundaries) [bound].edges [1],
                        (*boundaries) [bound].edges [2]);
#endif
              }
            else
              {
                // Boundary found, add to it; it's likely to be rare enough
                // that we need not worry about the realloc overhead
                int e = (*boundaries) [bound].n_edges;
#ifdef DEBUG
                printf ("Ternary boundary %d adding edge %d\n",bound,e);
#endif
                // All of these reallocs are pretty expensive...
                // Later use compos to store size of array, sort out later
                if (!((*boundaries) [bound].edges = realloc
                      ((*boundaries) [bound].edges, (3*e+3)*sizeof (int))))
                  return 1;
                (*boundaries) [bound].edges [3*e]   = MIN3 (v0, v1, v2);
                (*boundaries) [bound].edges [3*e+1] = MID3 (v0, v1, v2);
                (*boundaries) [bound].edges [3*e+2] = MAX3 (v0, v1, v2);
                (*boundaries) [bound].n_edges += 1;
#ifdef DEBUG
                printf ("  Corners for bound %d edge %d: %d, %d, %d\n", bound,e,
                        (*boundaries) [bound].edges [3*e],
                        (*boundaries) [bound].edges [3*e+1],
                        (*boundaries) [bound].edges [3*e+2]);
#endif
              }

            // Need to handle the corner "boundaries" later...

            break;
          }

        case TWO_PHASE:
          {
            int b=0, bound=-1, p0, p1, p2, nphases=2, e=0,
              phase0 = MIN3 (points [v0].efunc, points [v1].efunc,
                             points [v2].efunc),
              phase1 = MAX3 (points [v0].efunc, points [v1].efunc,
                             points [v2].efunc);

            // Set p0 and p1 to index first and second vertices in same phase,
            // p2 to lone vertex in the other phase.
            if (points [v0].efunc == points [v1].efunc)
              {
                p0 = MIN (v0, v1);
                p1 = MAX (v0, v1);
                p2 = v2;

                // If all three are on the same phase, this is a miscibility
                // gap facet, set nphases to 1 and figure out which vertex is
                // across from the other two
                if (points [v0].efunc == points [v2].efunc)
                  {
                    double C02 = points [v0].C2, C03 = points [v0].C3,
                      C12 = points [v1].C2, C13 = points [v1].C3,
                      C22 = points [v2].C2, C23 = points [v2].C3,
                      G0 = points [v0].G, G1 = points [v1].G,
                      G2 = points [v2].G;

#ifdef DEBUG
                    printf ("Misc gap: ");
#endif
                    nphases=1;

                    if (free_energy ((C02+C12)/2,(C03+C13)/2,T,P,eparams+phase0) <
                        (G0+G1)/2)
                      {
#ifdef DEBUG
                        printf ("v2 alone, Gav=%g, Gmid=%g\n",
                                (G0+G1)/2, free_energy ((C02+C12)/2,(C03+C13)/2,T,P,eparams+phase0));
#endif
                        p0 = MIN (v0, v1);
                        p1 = MAX (v0, v1);
                        p2 = v2;
                      }
                    else if (free_energy ((C02+C22)/2,(C03+C23)/2,T,P,eparams+phase0) <
                             (G0+G2)/2)
                      {
#ifdef DEBUG
                        printf ("v1 alone, Gav=%g, Gmid=%g\n",
                                (G0+G2)/2, free_energy ((C02+C22)/2,(C03+C23)/2,T,P,eparams+phase0));
#endif
                        p0 = MIN (v0, v2);
                        p1 = MAX (v0, v2);
                        p2 = v1;
                      }
                    else
                      {
#ifdef DEBUG
                        printf ("v0 alone, Gav=%g, Gmid=%g\n",
                                (G1+G2)/2, free_energy ((C22+C12)/2,(C23+C13)/2,T,P,eparams+phase0));
#endif
                        p0 = MIN (v1, v2);
                        p1 = MAX (v1, v2);
                        p2 = v0;
                      }
                  }
              }
            else if (points [v0].efunc == points [v2].efunc)
              {
                p0 = MIN (v0, v2);
                p1 = MAX (v0, v2);
                p2 = v1;
              }
            else
              {
                p0 = MIN (v1, v2);
                p1 = MAX (v1, v2);
                p2 = v0;
              }

            // Check for a boundary with tie lines between these two phases
            for (b=0; b<*n_bounds; b++)
              if ((*boundaries) [b].n_phases == nphases &&
                  (*boundaries) [b].type == TIE_SIMPLICES &&
                  (*boundaries) [b].phases[0] == phase0 &&
                  (*boundaries) [b].phases[nphases-1] == phase1)
                bound = b;

            if (bound == -1)
              {
                bound = *n_bounds;
                // No boundary with these two phases, make a new one
#ifdef DEBUG
                printf ("Making new tie-line boundary %d\n", *n_bounds);
#endif
                if (!(*boundaries = realloc
                      (*boundaries, (*n_bounds+1) * sizeof (phase_boundary))))
                  return 1;

                (*boundaries) [bound].compos = 2;
                (*boundaries) [bound].n_phases = nphases;
                if (!((*boundaries) [bound].phases =
                      malloc (nphases*sizeof (int))))
                  return 1;
                (*boundaries) [bound].phases [0] = phase0;
                (*boundaries) [bound].phases [nphases-1] = phase1;
                (*boundaries) [bound].n_edges = 0;
                (*boundaries) [bound].edges = NULL;
                (*boundaries) [bound].type = TIE_SIMPLICES;
                (*n_bounds) ++;
                e=0;
              }
            else
              e = (*boundaries) [bound].n_edges;

            // All of these reallocs are pretty expensive...
            // Later use compos to store size of array, sort out later
            if (!((*boundaries) [bound].edges = realloc
                  ((*boundaries) [bound].edges, (2*e+4) * sizeof (int))))
              return 1;
            (*boundaries) [bound].n_edges += 2;

            (*boundaries) [bound].edges [2*e]   = MIN (p0, p2);
            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p2);
            (*boundaries) [bound].edges [2*e+2] = MIN (p1, p2);
            (*boundaries) [bound].edges [2*e+3] = MAX (p1, p2);
#ifdef DEBUG
            printf ("  Corners for bound %d edges %d,%d: %d,%d, %d,%d\n",
                    bound, e, e+1, (*boundaries) [bound].edges [2*e],
                    (*boundaries) [bound].edges [2*e+1],
                    (*boundaries) [bound].edges [2*e+2],
                    (*boundaries) [bound].edges [2*e+3]);
#endif

            // Check for a boundary with an edge on the "majority" phase
            for (b=0, bound=-1; b<*n_bounds; b++)
              if ((*boundaries) [b].n_phases == nphases &&
                  (*boundaries) [b].type == ONE_PHASE_EDGE &&
                  (*boundaries) [b].phases[0] == points[p0].efunc &&
                  (*boundaries) [b].phases[nphases-1] == points[p2].efunc)
                bound = b;

            if (bound == -1)
              {
                // No boundary with these two phases, make a new one
                bound = *n_bounds;
#ifdef DEBUG
                printf ("Making new phase edge boundary %d\n", *n_bounds);
#endif
                if (!(*boundaries = realloc
                      (*boundaries, (bound+1) * sizeof (phase_boundary))))
                  return 1;

                (*boundaries) [bound].compos = 2;
                (*boundaries) [bound].n_phases = nphases;
                if (!((*boundaries) [bound].phases =
                      malloc (nphases*sizeof (int))))
                  return 1;
                (*boundaries) [bound].phases [0] = points[p0].efunc;
                (*boundaries) [bound].phases [nphases-1] = points[p2].efunc;
                (*boundaries) [bound].n_edges = 0;
                (*boundaries) [bound].edges = NULL;
                (*boundaries) [bound].type = ONE_PHASE_EDGE;
                (*n_bounds) ++;
                e=0;
              }
            else
              e = (*boundaries) [bound].n_edges;

            // All of these reallocs are pretty expensive...
            if (!((*boundaries) [bound].edges = realloc
                  ((*boundaries) [bound].edges, (2*e+2)*sizeof (int))))
              return 1;
            (*boundaries) [bound].edges [2*e]   = MIN (p0, p1);
            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p1);
            (*boundaries) [bound].n_edges += 1;
#ifdef DEBUG
            printf ("  Corners for bound %d edge %d: %d,%d\n",
                    bound, e, (*boundaries) [bound].edges [2*e],
                    (*boundaries) [bound].edges [2*e+1]);
#endif

            break;
          }

        case CRITICAL_POINT:
          {
#ifdef DEBUG
            printf ("Critical point case\n");
#endif
            int b=0, bound=-1, p0, p1, p2, e=0, phase = points [v0].efunc;
            double C02 = points [v0].C2, C03 = points [v0].C3,
              C12 = points [v1].C2, C13 = points [v1].C3,
              C22 = points [v2].C2, C23 = points [v2].C3,
              G0 = points [v0].G, G1 = points [v1].G,
              G2 = points [v2].G;

            if (free_energy ((C02+C12)/2,(C03+C13)/2,T,P,eparams+phase) >
                (G0+G1)/2)
              {
#ifdef DEBUG
                printf ("v2 alone, Gav=%g, Gmid=%g\n",
                        (G0+G1)/2, free_energy ((C02+C12)/2,(C03+C13)/2,T,P,eparams+phase));
#endif
                p0 = MIN (v0, v1);
                p1 = MAX (v0, v1);
                p2 = v2;
              }
            else if (free_energy ((C02+C22)/2,(C03+C23)/2,T,P,eparams+phase) >
                     (G0+G2)/2)
              {
#ifdef DEBUG
                printf ("v1 alone, Gav=%g, Gmid=%g\n",
                        (G0+G2)/2, free_energy ((C02+C22)/2,(C03+C23)/2,T,P,eparams+phase));
#endif
                p0 = MIN (v0, v2);
                p1 = MAX (v0, v2);
                p2 = v1;
              }
            else
              {
#ifdef DEBUG
                printf ("v0 alone, Gav=%g, Gmid=%g\n",
                        (G1+G2)/2, free_energy ((C22+C12)/2,(C23+C13)/2,T,P,eparams+phase));
#endif
                p0 = MIN (v1, v2);
                p1 = MAX (v1, v2);
                p2 = v0;
              }

            // Check for a boundary with tie lines across this miscibility gap
            for (b=0; b<*n_bounds; b++)
              if ((*boundaries) [b].n_phases == 1 &&
                  (*boundaries) [b].type == TIE_SIMPLICES &&
                  (*boundaries) [b].phases[0] == phase)
                bound = b;

            if (bound == -1)
              {
                bound = *n_bounds;
                // No tie-line boundary in this phase, make a new one
#ifdef DEBUG
                printf ("Making new tie-line boundary %d\n", *n_bounds);
#endif
                if (!(*boundaries = realloc
                      (*boundaries, (*n_bounds+1) * sizeof (phase_boundary))))
                  return 1;

                (*boundaries) [bound].compos = 2;
                (*boundaries) [bound].n_phases = 1;
                if (!((*boundaries) [bound].phases = malloc (sizeof (int))))
                  return 1;
                (*boundaries) [bound].phases [0] = phase;
                (*boundaries) [bound].n_edges = 0;
                (*boundaries) [bound].edges = NULL;
                (*boundaries) [bound].type = TIE_SIMPLICES;
                (*n_bounds) ++;
                e=0;
              }
            else
              e = (*boundaries) [bound].n_edges;

            // All of these reallocs are pretty expensive...
            // Later use compos to store size of array, sort out later
            if (!((*boundaries) [bound].edges = realloc
                  ((*boundaries) [bound].edges, (2*e+2) * sizeof (int))))
              return 1;
            (*boundaries) [bound].n_edges += 1;

            (*boundaries) [bound].edges [2*e]   = MIN (p0, p1);
            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p1);
#ifdef DEBUG
            printf ("  Corners for bound %d edges %d,%d: %d,%d\n",
                    bound, e, e+1, (*boundaries) [bound].edges [2*e],
                    (*boundaries) [bound].edges [2*e+1]);
#endif

            // Check for a boundary with an edge on this miscibility gap
            for (b=0, bound=-1; b<*n_bounds; b++)
              if ((*boundaries) [b].n_phases == 1 &&
                  (*boundaries) [b].type == ONE_PHASE_EDGE &&
                  (*boundaries) [b].phases[0] == phase)
                bound = b;

            if (bound == -1)
              {
                // No edge boundary with this miscibility gap, make a new one
                bound = *n_bounds;
#ifdef DEBUG
                printf ("Making new phase edge boundary %d\n", *n_bounds);
#endif
                if (!(*boundaries = realloc
                      (*boundaries, (bound+1) * sizeof (phase_boundary))))
                  return 1;

                (*boundaries) [bound].compos = 2;
                (*boundaries) [bound].n_phases = 1;
                if (!((*boundaries) [bound].phases = malloc (sizeof (int))))
                  return 1;
                (*boundaries) [bound].phases [0] = phase;
                (*boundaries) [bound].n_edges = 0;
                (*boundaries) [bound].edges = NULL;
                (*boundaries) [bound].type = ONE_PHASE_EDGE;
                (*n_bounds) ++;
                e=0;
              }
            else
              e = (*boundaries) [bound].n_edges;

            // All of these reallocs are pretty expensive...
            if (!((*boundaries) [bound].edges = realloc
                  ((*boundaries) [bound].edges, (2*e+4)*sizeof (int))))
              return 1;
            (*boundaries) [bound].edges [2*e]   = MIN (p0, p2);
            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p2);
            (*boundaries) [bound].edges [2*e+2] = MIN (p1, p2);
            (*boundaries) [bound].edges [2*e+3] = MAX (p1, p2);
            (*boundaries) [bound].n_edges += 2;
#ifdef DEBUG
            printf ("  Corners for bound %d edge %d: %d,%d, %d,%d\n",
                    bound, e, (*boundaries) [bound].edges [2*e],
                    (*boundaries) [bound].edges [2*e+1],
                    (*boundaries) [bound].edges [2*e+2],
                    (*boundaries) [bound].edges [2*e+3]);
#endif

            break;
          }

        default:
          {
            printf ("Default case: type %s, phases %d, %d, %d\n",
                    (facets [f].type == ONE_PHASE) ? "1-phase" :
                    ((facets [f].type == CRITICAL_POINT) ? "critical" :
                     ((facets [f].type == TWO_PHASE) ? "2-phase" :
                      ((facets [f].type == THREE_PHASE) ? "3-phase" :
                       "unknown"))), points [facets[f].vertex[0]].efunc,
                    points [facets[f].vertex[1]].efunc,
                    points [facets[f].vertex[2]].efunc);
            printf ("  Coords: %d %g,%g, %d %g,%g, %d %g,%g\n",
                    facets[f].vertex[0], points [facets[f].vertex[0]].C2,
                    points [facets[f].vertex[0]].C3,
                    facets[f].vertex[1], points [facets[f].vertex[1]].C2,
                    points [facets[f].vertex[1]].C3,
                    facets[f].vertex[2], points [facets[f].vertex[2]].C2,
                    points [facets[f].vertex[2]].C3);
          }
        }
    }

  return 0;
}


/*++++++++++++++++++++++++++++++++++++++
  This is what it sounds like.

  char *hullQHullVersion It returns the QHull version string.
  ++++++++++++++++++++++++++++++++++++++*/

char *hullQHullVersion ()
{
  return qh_version;
}