Changeset 430 for trunk

Timestamp:

02/24/2009 09:35:26 PM (3 years ago)

Author:

powell

Message:

Construct all phase boundaries in the convex hull (almost works).

Location:

trunk/matml/src/ternary

Files:

• ChangeLog (modified) (1 diff)
• qhull.c (modified) (1 diff)
• ternary.c (modified) (3 diffs)
• ternary.h (modified) (2 diffs)

trunk/matml/src/ternary/ChangeLog

@@ -3 +3 @@
   * Overhauled the free energy interface and functions, adding derivative
     implementations.
-  * Added phase_boundary typedef struct.
+  * Added phase_boundary object.
   * Changed hull calculation interface to only hold the qhull variables during
     hull calculation (with a lock to prevent attempts at multiple usage).
   * Phase boundary refinement using a cool, and I think novel, Newton energy
     minimization algorithm.
+  * Report all phase boundaries in the hull as phase_boundary objects.
   * Added WITH_DERIVATIVES and NO_INFINITY options to free energies().
-  * Added spinodal.c with calc_spinodal function.
+  * Added spinodal.c with calc_spinodal function returning a phase_boundary.
 
  --

trunk/matml/src/ternary/qhull.c

@@ -677 +677 @@
 
 
+#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.
+
+  ternary_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.
+
+  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
+(ternary_point *points, int n_points, hull_facet *facets, int n_facets,
+ phase_boundary **boundaries, int *n_bounds)
+{
+  int f;
+  *n_bounds = 0;
+
+  printf ("Starting ReturnPhaseBoundaries with %d facets, %d points\n",
+          n_facets, n_points);
+
+  printf ("Last facet %d vertices %d, %d, %d\n", n_facets-1,
+          facets[n_facets-1].vertex[0], facets[n_facets-1].vertex[1],
+          facets[n_facets-1].vertex[2]);
+
+  for (f=0; f<n_facets; f++)
+    {
+      int v0 = facets[f].vertex[0], v1 = facets[f].vertex[1],
+        v2 = facets[f].vertex[2];
+
+      printf ("Last facet %d vertices %d, %d, %d\n", n_facets-1,
+              facets[n_facets-1].vertex[0], facets[n_facets-1].vertex[1],
+              facets[n_facets-1].vertex[2]);
+
+      printf ("Facet %d: type %d, vertices %d ",f, facets [f].type, v0);
+      fflush (stdout);
+      printf ("(%d), %d (%d), %d(%d)\n",
+              points[v0].efunc, v1,
+              points[v1].efunc, v2, points[v2].efunc);
+      fflush (stdout);
+
+      switch (facets [f].type)
+        {
+        case ONE_PHASE:
+          {
+            printf ("Type 0 facet\n"); fflush (stdout);
+            // 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].flags == 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;
+                printf ("Making new three-phase boundary %d\n", *n_bounds);
+                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].flags = 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) ++;
+              }
+            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;
+                printf ("Ternary boundary %d adding edge %d\n",bound,e);
+                // 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;
+              }
+
+            // Need to handle the corner "boundaries" later...
+
+            break;
+          }
+
+        case TWO_PHASE:
+          {
+            int b=0, bound=-1, p0, p1, p2,
+              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, need to check midpoint
+                // energies to determine which are one- and two-phase segments.
+                // Punt for the default case later!
+                if (points [v0].efunc == points [v2].efunc)
+                  { printf ("Miscibility gap!\n"); break; }
+              }
+            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 == 2 &&
+                  (*boundaries) [b].flags == TIE_SIMPLICES &&
+                  (*boundaries) [b].phases[0] == phase0 &&
+                  (*boundaries) [b].phases[1] == phase1)
+                bound = b;
+
+            if (bound == -1)
+              {
+                bound = *n_bounds;
+                // No boundary with these two phases, make a new one
+                printf ("Making new tie-line boundary %d\n", *n_bounds);
+                if (!(*boundaries = realloc
+                      (*boundaries, (*n_bounds+1) * sizeof (phase_boundary))))
+                  return 1;
+
+                (*boundaries) [bound].compos = 2;
+                (*boundaries) [bound].n_phases = 2;
+                if (!((*boundaries) [bound].phases =
+                      malloc (2*sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].phases [0] = phase0;
+                (*boundaries) [bound].phases [1] = phase1;
+                (*boundaries) [bound].n_edges = 2;
+                if (!((*boundaries) [bound].edges =
+                      malloc (4*sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].edges [0] = MIN (p0, p2);
+                (*boundaries) [bound].edges [1] = MAX (p0, p2);
+                (*boundaries) [bound].edges [2] = MIN (p1, p2);
+                (*boundaries) [bound].edges [3] = MAX (p1, p2);
+                (*boundaries) [bound].flags = TIE_SIMPLICES;
+                (*n_bounds) ++;
+              }
+            else
+              {
+                int e = (*boundaries) [bound].n_edges;
+                printf ("Tie-line boundary %d adding edges %d,%d\n",bound,e,e+1);
+                // 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].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;
+              }
+
+            // 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 == 2 &&
+                  (*boundaries) [b].flags == ONE_PHASE_EDGE &&
+                  (*boundaries) [b].phases[0] == points[p0].efunc &&
+                  (*boundaries) [b].phases[1] == points[p2].efunc)
+                bound = b;
+
+            if (bound == -1)
+              {
+                // No boundary with these two phases, make a new one
+                bound = *n_bounds;
+                printf ("Making new phase edge boundary %d\n", *n_bounds);
+                if (!(*boundaries = realloc
+                      (*boundaries, (bound+1) * sizeof (phase_boundary))))
+                  return 1;
+
+                (*boundaries) [bound].compos = 2;
+                (*boundaries) [bound].n_phases = 2;
+                if (!((*boundaries) [bound].phases =
+                      malloc (2*sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].phases [0] = points[p0].efunc;
+                (*boundaries) [bound].phases [1] = points[p2].efunc;
+                (*boundaries) [bound].n_edges = 1;
+                if (!((*boundaries) [bound].edges =
+                      malloc (2*sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].edges [0] = MIN (p0, p1);
+                (*boundaries) [bound].edges [1] = MAX (p0, p1);
+                (*boundaries) [bound].flags = ONE_PHASE_EDGE;
+                (*n_bounds) ++;
+              }
+            else
+              {
+                int e = (*boundaries) [bound].n_edges;
+                printf ("Phase edge boundary %d adding edge %d\n",bound,e);
+                // 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;
+              }
+
+            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 ("Done!\n");
+  return 0;
+}
+
+
 /*++++++++++++++++++++++++++++++++++++++
   This is what it sounds like.

trunk/matml/src/ternary/ternary.c

@@ -28 +28 @@
 int main (int argc, char *argv[])
 {
-  int i,j, index, loop_max=20, numpoints, *verts, hullnumfacets;
+  int i,j, index, loop_max=20, numpoints, *verts, hullnumfacets, numbounds;
   hull_facet *hullfacets = NULL;
   char gv_version[100];
@@ -34 +34 @@
   double T=0.7, P=1.;
   ternary_point *points;
-  phase_boundary solid_spinodal;
+  phase_boundary solid_spinodal, *allbounds=NULL;
   /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */
   energy_params /* Solid, liquid */
@@ -120 +120 @@
     { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
 
+  if (i=hullReturnPhaseBoundaries (points,numpoints, hullfacets,hullnumfacets,
+                                   &allbounds, &numbounds))
+    { printf ("main: error %d in hullReturnPhaseBoundaries\n", i); exit (i); }
+
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 0", "bd0", "-face +edge",
+       numpoints, points, allbounds))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 1", "bd1", "-face +edge",
+       numpoints, points, allbounds+1))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 2", "bd2", "-face +edge",
+       numpoints, points, allbounds+2))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 3", "bd3", "-face +edge",
+       numpoints, points, allbounds+3))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
   printf ("Press <return> to close up... ");
   fflush (stdout);

trunk/matml/src/ternary/ternary.h

@@ -105 +105 @@
 typedef enum {
   ONE_PHASE_EDGE=0, /*+ One phase's boundary edges +*/
-  TIE_SIMPLICES,    /*+ Tie lines (ternary) or other simplices between two
-                      phases +*/
+  TIE_SIMPLICES,    /*+ Tie lines between two phases, triangles between three,
+                      etc. +*/
   SPINODAL          /*+ Spinodal boundary +*/
 } phase_boundary_flags;
 
 typedef struct {
-  int compos;   /*+ Number of components of the space (ternary binary etc.) +*/
+  int compos;   /*+ Number of points in a boundary edge +*/
   int *edges;   /*+ Edges with indices of ternary_points comprising a phase
                   boundary; each edge has compos points +*/
@@ -174 +174 @@
  energy_params *eparams, int nparams, double T,double P,
  double newton_tolerance, double vertex_tolerance, int one_phase_refine);
+int hullReturnPhaseBoundaries
+(ternary_point *points, int n_points, hull_facet *facets, int n_facets,
+ phase_boundary **boundaries, int *n_bounds);
 char *hullQHullVersion ();