Changeset 421

Timestamp:

02/23/2009 06:38:33 PM (3 years ago)

Author:

powell

Message:

Refactored hullCalculate interface, adding facet type, breaking refine.

Location:

trunk/matml/src/ternary

Files:

• ChangeLog (modified) (1 diff)
• qhull.c (modified) (4 diffs)
• ternary.c (modified) (5 diffs)
• ternary.h (modified) (2 diffs)
• ternary.tex.in (modified) (1 diff)

trunk/matml/src/ternary/ChangeLog

@@ -4 +4 @@
     implementations.
   * Added phase_boundary typedef struct.
-  * Changed hull calculation interface to permit (multiple) refinement between
-    original calculation and return.
+  * Changed hull calculation interface to only hold the qhull variables during
+    hull calculation, and permit (multiple) phase boundary refinement.
   * Added WITH_DERIVATIVES and NO_INFINITY options to free energies().
   * Added spinodal.c with calc_spinodal function.

trunk/matml/src/ternary/qhull.c

@@ -12 +12 @@
 #include <stdlib.h>
 #include <math.h>
-
-
-static coordT *qpoints; /*+ Point coordinate storage for qhull +*/
-
-
-/*++++++++++++++++++++++++++++++++++++++
-  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 does not return anything, but leaves all of the facets in
-  the qhull structures.
-
-  int hullCalculate It returns zero or an error code.
-
-  int dim Number of dimensions (should always be 3 for ternary).
-
-  ternary_point *points Point information.
-
-  int numpoints Number of points of which to calculate the convex hull.
-  ++++++++++++++++++++++++++++++++++++++*/
-
-int hullCalculate (int dim, ternary_point *points, int numpoints)
-{
-  int i, corner1, corner2, corner3;
-
-  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; }
-  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 corners
-      if (points[i].C2 == 0. && points[i].C3 == 0.)
-        corner1 = i;
-      if (points[i].C2 == 0. && points[i].C3 == 1.)
-        corner2 = i;
-      if (points[i].C2 == 1. && points[i].C3 == 0.)
-        corner3 = i;
-    }
-
-  // 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
-
-  return 0;
-}
 
 
@@ -145 +80 @@
 
 /*++++++++++++++++++++++++++++++++++++++
+  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).
+
+  ternary_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, ternary_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; }
+  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.) &&
+          (corner2 == -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");
+      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++;
+        }
+    }
+
+#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
+
+  free (qpoints);
+  qh_freeqhull (True);
+
+  return 0;
+}
+
+
+/*++++++++++++++++++++++++++++++++++++++
   This uses Newton iteration to refine a single point within a triangle.
 
@@ -534 +650 @@
     return 1;
           
-  qh_freeqhull (True);
-  free (qpoints);
+  //qh_freeqhull (True);
+  //free (qpoints);
 
   for (i=0; i<numnewpoints; i++)
@@ -565 +681 @@
   free(newpoints);
 
-  return hullCalculate (3, *points, *numpoints);
-}
-
-
-/*++++++++++++++++++++++++++++++++++++++
-  This returns the facets of the convex hull comprising the minimum energy
-  surface.  It removes the sides and top, currently by getting rid of
-  non-simplical facets (sides) and a facet whose projection area is the entire
-  triangle (top).
-
-  int hullReturnFacets It returns zero or an error code.
-
-  energy_params *eparams Collection of energy parameter structures.
-
-  int nparams Number of energy parameter structures in eparams.
-
-  double T Temperature.
-
-  double P Pressure.
-
-  int *numfacets Pointer to an integer which will contain the number of facets
-  in the convex hull on return.
-
-  int **facetverts Pointer which contains the array of facet vertices on
-  return.  This is changed by realloc() each time the function is called, so it
-  should be NULL or a valid array.
-
-  facet_type **facettypes Pointer which contains the array of facet types on
-  return.  This is changed by realloc() each time the function is called, so it
-  should be NULL or a valid array.
-  ++++++++++++++++++++++++++++++++++++++*/
-
-int hullReturnFacets
-(energy_params *eparams, int nparams, double T,double P,
- int *numfacets, int **facetverts, facet_type **facettypes)
-{
-  int i=0;
-  facetT *facet;
-
-  /*+ After qhull runs, this reallocates the array pointed to by facetverts,
-    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 ("facetverts: pre-realloc 0x%lx, ", *facetverts);
-#endif
-  if (!(*facetverts = realloc (*facetverts, *numfacets * 3 * sizeof (int))))
-    { printf ("hullReturnFacets: could not reallocate memory for vertices\n");
-      return -1; }
-#ifdef DEBUG
-  printf ("  post-realloc 0x%lx\n", *facetverts);
-  printf ("facettypes: pre-realloc 0x%lx, ", *facettypes);
-#endif
-  if (!(*facettypes = realloc (*facettypes, *numfacets * sizeof (facet_type))))
-    { printf ("hullReturnFacets: could not reallocate memory for facet types\n");
-      return -1; }
-#ifdef DEBUG
-  printf ("  post-realloc 0x%lx\n", *facettypes);
-#endif
-
-  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)
-            {
-              (*facetverts) [3*i+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
-        {
-          (*facettypes) [i] = facettype (corners, eparams, nparams, T, P);
-          printf ("  Type %d\n", (*facettypes) [i]);
-          i++;
-        }
-    }
-
-#ifdef DEBUG
-  printf ("After loop: numfacets=%d, i=%d\n", *numfacets, i);
-#endif
-
-  /* free (qpoints);
-     qh_freeqhull (True); */
-
-  return 0;
+  //return hullCalculate (3, *points, *numpoints);
 }
 

trunk/matml/src/ternary/ternary.c

@@ -28 +28 @@
 int main (int argc, char *argv[])
 {
-  int i,j, index, loop_max=20, numpoints, *verts,
-    hullnumverts,*hullverts=NULL;
-  facet_type *hulltypes=NULL;
+  int i,j, index, loop_max=20, numpoints, *verts, hullnumfacets;
+  hull_facet *hullfacets = NULL;
   char gv_version[100];
   FILE *pfd = NULL;
@@ -82 +81 @@
     { printf ("main: error %d in spinodal calculation\n", i); exit (i); }
 
-  printf ("First edge: %d %d\n", solid_spinodal.edges [0], solid_spinodal.edges [1]);
-
   /* Scale all free energy values */
   if (i=scale_energy_array (numpoints, points, 0., -.2, 1., .9, NULL, NULL))
@@ -101 +98 @@
   /* Calculate and refine the lower convex hull */
   printf ("qhull version: %s\n", hullQHullVersion());
-  if (i=hullCalculate (3, points, numpoints))
+  if (i=hullCalculate (3, points, numpoints, allparams, 2, T, P,
+                       &hullfacets, &hullnumfacets))
     { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
-  if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
+  /*if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
                     1e-10, 1e-7, 0))
     { printf ("main: error %d in hullRefine\n", i); exit (i); }
@@ -111 +109 @@
   if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
                     1e-10, 1e-7, 0))
-    { printf ("main: error %d in hullRefine\n", i); exit (i); }
-  if (i=hullReturnFacets (allparams, 2, T, P, &hullnumverts, &hullverts,
-                          &hulltypes))
-    { printf ("main: Error %d in hullReturnFacets\n", i); exit (i); }
+    { printf ("main: error %d in hullRefine\n", i); exit (i); }*/
 
   /* Take out hull triangles which are on the main function (in one phase) */
-  for (i=0; i<hullnumverts; i++)
+  /*for (i=0; i<hullnumfacets; i++)
     {
-      /* double Gav = (points [hullverts [3*i]].G +
-                    points [hullverts [3*i+1]].G +
-                    points [hullverts [3*i+2]].G)/3.,
-        C2av = (points [hullverts [3*i]].C2 +
-                points [hullverts [3*i+1]].C2 +
-                points [hullverts [3*i+2]].C2)/3.,
-        C3av = (points [hullverts [3*i]].C3 +
-                points [hullverts [3*i+1]].C3 +
-                points [hullverts [3*i+2]].C3)/3.;
 
-#ifdef DEBUG
-      printf ("Facet %d vertices: %d %d %d\n", i, hullverts[3*i],
-              hullverts[3*i+1], hullverts[3*i+2]);
-              #endif */
-
-      /* Test whether hull facet is in one-phase region: free energy at
-         centroid is below facet; if so, remove it */
-      /*if (fmin (free_energy (C2av, C3av, T, P, allparams),
-        free_energy (C2av, C3av, T, P, allparams+1)) < Gav)*/
-      if (hulltypes [i] == ONE_PHASE)
+      // Remove one-phase facets to display just binodal and ternary
+      if (hullfacets [i].type == ONE_PHASE)
         {
-          for (j=i; j<hullnumverts-1; j++)
-            {
-              hullverts [3*j]   = hullverts [3*j+3];
-              hullverts [3*j+1] = hullverts [3*j+4];
-              hullverts [3*j+2] = hullverts [3*j+5];
-              hulltypes [j] = hulltypes [j+1];
-            }
-          hullnumverts--;
+          for (j=i; j<hullnumfacets-1; j++)
+            hullfacets [j] = hullfacets [j+1];
+          hullnumfacets--;
           i--;
         }
-    }
+    } */
 
-  printf ("Binodal region has %d facets\n", hullnumverts);
+  printf ("Binodal region has %d facets\n", hullnumfacets);
 
   if (i=scale_energy_array (numpoints, points, 0., -.2, 1., .9, NULL, NULL))
     { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
 
-  if (i=GeomviewDisplayTriangleCOFF
-      (pfd, "Binodal and Tie Lines", "ech", "-face +edge",
-       numpoints, points, hullnumverts, hullverts))
-    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
-
   printf ("Press <return> to close up... ");
+  fflush (stdout);
   fgets (gv_version, 99, stdin);
 
@@ -168 +137 @@
   free (points);
   free (verts);
-  free (hullverts);
+  free (hullfacets);
 
   return 0;

trunk/matml/src/ternary/ternary.h

@@ -130 +130 @@
   CRITICAL_POINT,
   TWO_PHASE,
-  THREE_PHASE
+  THREE_PHASE,
+  ALL_TYPES=127
 } facet_type;
+
+typedef struct {
+  int vertex [3];
+  facet_type type;
+} hull_facet;
 
 /* From freenergy.c */
@@ -157 +163 @@
 
 /* From qhull.c */
-int hullCalculate (int dim, ternary_point *points, int numpoints);
+int hullCalculate
+(int dim, ternary_point *points, int numpoints, energy_params *eparams,
+ int nparams, double T, double P, hull_facet **facets, int *numfacets);
 int hullRefine
 (ternary_point **points, int *numpoints, energy_params *eparams, int nparams,
  double T,double P, double newton_tolerance, double vertex_tolerance,
  int one_phase_refine);
-int hullReturnFacets
-(energy_params *eparams, int nparams, double T,double P,
- int *numfacets, int **facetverts, facet_type **facettypes);
 char *hullQHullVersion ();
 

trunk/matml/src/ternary/ternary.tex.in

@@ -41 +41 @@
 {\tt Ternary} is nearly complete, with the one un-implemented feature being to
 build and neatly visualize how free energy functions create a ternary phase
-diagram.  Also, {\tt qhull.c} needs re-architecting with its own ``facet''
-structure, such that we can do the hull analysis, store what we need, and free
-{\tt qhull}'s structure arrays and not need to rely on its loops beyond that.
-(This will also make {\tt qhull.c} functions thread-safe.)
+diagram.
 
 \section{How it works}