Changeset 383

Timestamp:

06/29/08 22:28:08 (3 months ago)

Author:

powell

Message:

This sort of works, but I suspect an error in the freenergy.c derivatives.

Files:

• trunk/matml/src/ternary/qhull.c (modified) (14 diffs)
• trunk/matml/src/ternary/ternary.c (modified) (6 diffs)
• trunk/matml/src/ternary/ternary.h (modified) (1 diff)

trunk/matml/src/ternary/qhull.c

@@ -14 +14 @@
-
-
+static coordT *qpoints; /*+ Point coordinate storage for qhull +*/
+
+
-/*++++++++++++++++++++++++++++++++++++++
-  This initializes qhull and uses it to calculate the convex hull of a set of
@@ -32 +35 @@
-{
-  int i, corner1, corner2, corner3;
-  coordT *qpoints;
-
-  if (dim != 3)
@@ -141 +143 @@
-    }
-
-
+
+
-    {
-      double dC2, dC3;
-
-      /*+ Calculate the energy derivatives +*/
--1
+0
-
-      /*+ Subtract the facet slopes from the energy derivatives +*/
@@ -161 +164 @@
-      current.C3 += dC3;
-      distance = dC2*dC2 + dC3*dC3;
-
-       
+//
+      //
-    }
-
@@ -178 +181 @@
-  int hullRefine It returns zero or an error code.
-
-points Point information
+*points Pointer to point list
-
-  int numpoints Number of points of which to calculate the convex hull.
@@ -202 +205 @@
-
-int hullRefine
-points, int n
+*points, int *num
- double T,double P, double newton_tolerance, double vertex_tolerance,
- int one_phase_refine)
-{
-  facetT *facet;
+  ternary_point *newpoints;
+  int newpointsize=100, numnewpoints=0, i;
+
+  if (!(newpoints = malloc (100*sizeof(ternary_point))))
+    return 1;
-
-  /*+ Steps in the algorithm:
@@ -239 +247 @@
-              corners [3*j+2] = vertex->point[2];
-
-points
+(*points)
-            }
-          j++;
@@ -250 +258 @@
-        {
-          int eparamc=0;
-
-
-
-
+
+
-                  corners[0], corners[1], corners[3],
-                  corners[4], corners[6], corners[7]);
@@ -290 +296 @@
-              if (one_phase_refine || eparam[0] != eparamc)
-                {
-Will refine %s-phase triangle
+Refining %s-phase triangle at centroid
-                          (eparam[0]==eparamc) ? "one" : "multi");
-
+
+
+
+
+
+
+
+
-                }
-
@@ -328 +341 @@
-                      corners[3*j] + corners[3*j+1] == 1.)
-                    {
- centroid[0];
-      //
-      //
-      //
- centroid[1];
-      //
-      //
-      //
- centroid[2];
-      //
-      //
-      //
+
+       
+       
+       
+
+       
+       
+       
+
+       
+       
+       
-                    }
-
-
+
+
+
+
-                  printf ("    Refining point at %g, %g function %d\n",
-                          newcorners[3*j], newcorners[3*j+1], eparam[j]);
-2
+3
-                                T,P, 3,3, newton_tolerance);
-                  printf ("    -> %g, %g\n", newcorners[3*j],
@@ -352 +367 @@
-                  newcorners[3*j+2] = free_energy
-                    (newcorners[3*j],newcorners[3*j+1], T,P, eparams+eparam[j]);
-                    }
-                }
-
@@ -364 +378 @@
-
-              // Don't add the point, reallocate the points array instead!
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-              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]) >
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-              if ((newcorners[6]-newcorners[0])*(newcorners[6]-newcorners[0]) +
@@ -387 +419 @@
-                  (newcorners[7]-newcorners[4])*(newcorners[7]-newcorners[4]) +
-                  (newcorners[8]-newcorners[5])*(newcorners[8]-newcorners[5]) >
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            }
-        }
@@ -400 +442 @@
-    +*/
-
+  if (!(*points = realloc
+        (*points, (*numpoints + numnewpoints) * sizeof (ternary_point))))
+    return 1;
+          
+  printf ("Adding %d new points to hull\n", numnewpoints);
+  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];
+      // In case I try to implement this using qh_addpoint someday...
+      /*printf ("  Finding best facet for coords %g, %g, %g\n",
+              (*points) [(*numpoints)+i].C2,(*points) [(*numpoints)+i].C3,
+              (*points) [(*numpoints)+i].G);
+      facet = qh_findbestfacet (newcorners, True, &bestdist, &isoutside);
+      printf ("  Adding to %d facet 0x%lx: %g, %g, %g\n", isoutside, facet, newcorners[0], newcorners[1],
+              newcorners[2]);
+      if (isoutside)
+        printf ("  Result: %d\n", qh_addpoint (newcorners, facet, True));
+      else
+      (*numpoints)--;*/
+    }
+  *numpoints += numnewpoints;
+
+  qh_freeqhull (True);
+  free (qpoints);
+
-  printf ("Done, returning\n");
-
-
+
+
+
-}
-

trunk/matml/src/ternary/ternary.c

@@ -28 +28 @@
-int main (int argc, char *argv[])
-{
- j, index, loop_max=50, *verts, hullnumverts, 
+j, index, loop_max=50, numpoints, *verts, hullnumverts,
-  char gv_version[100];
-  FILE *pfd = NULL;
@@ -41 +41 @@
-  if (argc>1)
-    sscanf (argv[1], "%d", &loop_max);
+  numpoints=(loop_max+1)*(loop_max+2);
-
-  /* Allocate array pairs for ternary points and triangle vertices */
-calloc ((loop_max+1)*(loop_max+2),
+malloc (numpoints *
-    { printf ("Cannot allocate point coordinate array\n"); exit (1); }
-calloc (loop_max*loop_max * 6,
+malloc (loop_max*loop_max * 6 *
-    { printf ("Cannot allocate triangle vertex array\n"); exit (1); }
-
@@ -56 +57 @@
-  if (i=init_triangle_array (loop_max, points))
-    { printf ("main: Error %d in init_triangle_array\n", i); exit (i); }
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in init_triangle_array\n", i); exit (i); }
-
@@ -63 +64 @@
-    { printf ("main: Error %d in init_triangle_vertices\n", i); exit (i); }
-  if (i=init_triangle_vertices (loop_max, verts + loop_max*loop_max*3,
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in init_triangle_vertices\n", i); exit (i); }
-
-  /* Calculate free energies */
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in free_energies\n", i); exit (i); }
-
-
+
-    { printf ("main: Error %d in free_energies\n", i); exit (i); }
-
-  /* Scale all free energy values */
-
-
+
-    { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
-
@@ -81 +80 @@
-  if (i=GeomviewDisplayTriangleCOFF
-      (pfd, "Ternary Free Energy", "tfe", "shading smooth",
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
-
-  /* Calculate and refine the lower convex hull */
-  printf ("qhull version: %s\n", hullQHullVersion());
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
-points, (loop_max+1)*(loop_max+2)
+&points, &numpoints
-                    1e-10, 1e-7, 0))
-    { printf ("main: error %d in hullRefine\n", i); exit (i); }
@@ -130 +129 @@
-  printf ("Binodal region has %d facets\n", hullnumverts);
-
+  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",
-(loop_max+1)*(loop_max+2)
+numpoints
-    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
-

trunk/matml/src/ternary/ternary.h

@@ -121 +121 @@
-int hullCalculate (int dim, ternary_point *points, int numpoints);
-int hullRefine
-points, int 
+*points, int *
- double T,double P, double newton_tolerance, double vertex_tolerance,
- int one_phase_refine);