Changeset 355

Timestamp:

03/27/08 18:02:07 (8 months ago)

Author:

powell

Message:

Better check of whether hull facet is in one-phase region.

Location:

trunk/matml/src/ternary

Files:

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

trunk/matml/src/ternary/ChangeLog

@@ -6 +6 @@
   * Moved many O(N) functions from main() out to new book.c file.
   * Multiple free energy functions work.
-  * TODO: better single-phase check in ternary.c; 0.3 refine qhull concave
-    facets.
+  * Better check to determine whether hull facets are in one-phase region.
 
  --

trunk/matml/src/ternary/ternary.c

@@ -31 +31 @@
   char gv_version[100], *qh_version;
   FILE *pfd = NULL;
-  double T=1.5;
+  double T=0.75;
   ternary_point *points;
   /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */
@@ -76 +76 @@
   /* Scale all free energy values to (0->1) */
   if (i=scale_energy_array ((loop_max+1)*(loop_max+2), points,
-                              0., -1.5, 1., 1., NULL, NULL))
+                              0., -.2, 1., .9, NULL, NULL))
     { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
 
@@ -94 +94 @@
   for (i=0; i<hullnumverts; i++)
     {
-      double x0 = points [hullverts [3*i]].C2,
-        y0 = points [hullverts [3*i]].C3,
-        x1 = points [hullverts [3*i+1]].C2,
-        y1 = points [hullverts [3*i+1]].C3,
-        x2 = points [hullverts [3*i+2]].C2,
-        y2 = points [hullverts [3*i+2]].C3;
+      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
-      /* XY-projection area test (only works for uniform triangles) */
-      /* Better test: is free energy at centroid above or below the facet? */
-      if (fabs ((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)) <= 1.1/loop_max/loop_max)
+
+      /* 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, &eparams1),
+                free_energy (C2av, C3av, T, &eparams2)) < Gav)
         {
           for (j=i; j<hullnumverts-1; j++)