Changeset 383

Timestamp:

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

Author:

powell

Message:

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

Location:

trunk/matml/src/ternary

Files:

• qhull.c (modified) (14 diffs)
• ternary.c (modified) (6 diffs)
• 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 @@
     }
 
-  while (distance >= newton_tolerance * newton_tolerance)
+  while (distance >= newton_tolerance * newton_tolerance &&
+         current.C2>=0. && current.C3>=0. && current.C2+current.C3<=1.)
     {
       double dC2, dC3;
 
       /*+ Calculate the energy derivatives +*/
-      free_energies (&current, 1, T, P, eparams, -1, 1);
+      free_energies (&current, 1, T, P, eparams, 0, 1);
 
       /*+ Subtract the facet slopes from the energy derivatives +*/
@@ -161 +164 @@
       current.C3 += dC3;
       distance = dC2*dC2 + dC3*dC3;
-      printf ("        Added %g,%g (dist %g) to go to %g,%g\n", dC2,dC3,
-              sqrt(distance), current.C2,current.C3);
+      //printf ("        Added %g,%g (dist %g) to go to %g,%g\n", dC2,dC3,
+      //      sqrt(distance), current.C2,current.C3);
     }
 
@@ -178 +181 @@
   int hullRefine It returns zero or an error code.
 
-  ternary_point *points Point information.
+  ternary_point **points Pointer to point list.
 
   int numpoints Number of points of which to calculate the convex hull.
@@ -202 +205 @@
 
 int hullRefine
-(ternary_point *points, int npoints, energy_params *eparams, int nparams,
+(ternary_point **points, int *numpoints, energy_params *eparams, int nparams,
  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];
 
-              eparam [j] = points[thepoint[j]].efunc;
+              eparam [j] = (*points)[thepoint[j]].efunc;
             }
           j++;
@@ -250 +258 @@
         {
           int eparamc=0;
-          realT bestdist;
-          boolT isoutside;
-
-          printf ("Triangle: %g,%g %g,%g %g,%g\n",
+
+          printf ("Facet 0x%lx: %g,%g %g,%g %g,%g\n", facet,
                   corners[0], corners[1], corners[3],
                   corners[4], corners[6], corners[7]);
@@ -290 +296 @@
               if (one_phase_refine || eparam[0] != eparamc)
                 {
-                  printf ("  Will refine %s-phase triangle\n",
+                  printf ("  Refining %s-phase triangle at centroid\n",
                           (eparam[0]==eparamc) ? "one" : "multi");
-                  // Change that when we add this as a new vertex
+                  if (numnewpoints >= newpointsize)
+                    newpoints = realloc
+                      (newpoints, (newpointsize+=100)*sizeof(ternary_point));
+                  newpoints [numnewpoints].C2 = centroid[0];
+                  newpoints [numnewpoints].C3 = centroid[1];
+                  newpoints [numnewpoints].G  = Gcenter;
+                  newpoints [numnewpoints].efunc = eparamc;
+                  numnewpoints++;
                 }
 
@@ -328 +341 @@
                       corners[3*j] + corners[3*j+1] == 1.)
                     {
-                      newcorners[3*j] = centroid[0];
-                      //0.9999999998 * corners[3*j] +
-                      //0.0000000001 * corners[3*((j+1)%3)] +
-                      //0.0000000001 * corners[3*((j+2)%3)];
-                      newcorners[3*j+1] = centroid[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] = centroid[2];
-                      //0.9999999998 * corners[3*j+2] +
-                      //0.0000000001 * corners[3*((j+1)%3)+2] +
-                      //0.0000000001 * corners[3*((j+2)%3)+2];
+                      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];
                     }
-                  else
-                    {
+                }
+
+              for (j=0; j<3; j++)
+                {
                   printf ("    Refining point at %g, %g function %d\n",
                           newcorners[3*j], newcorners[3*j+1], eparam[j]);
-                  NewtonRefine (newcorners+2*j, corners, eparams+eparam[j],
+                  NewtonRefine (newcorners+3*j, corners, eparams+eparam[j],
                                 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!
-              printf ("    Adding new point %g,%g,%g; result %d\n",
-                      newcorners[0], newcorners[1], newcorners[2],
-                      qh_addpoint ((pointT *) newcorners, 
-                                   qh_findbestfacet ((pointT *)newcorners,
-                                                     True, &bestdist, &isoutside),
-                                   True));
-
-              /*
+              if (newcorners[0]>=0. && newcorners[1]>=0. &&
+                  newcorners[0]+newcorners[1]<=1.)
+                {
+                  printf ("    Adding new point %g,%g,%g\n",
+                          newcorners[0], newcorners[1], newcorners[2]);
+                  if (numnewpoints >= newpointsize)
+                    newpoints = realloc
+                      (newpoints, (newpointsize+=100)*sizeof(ternary_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)
-                printf ("    Adding new point %g,%g,%g; result %d\n",
-                        newcorners[3], newcorners[4], newcorners[5],
-                        qh_addpoint ((pointT *) newcorners+3, facet, True));
+                  vertex_tolerance * vertex_tolerance &&
+                  newcorners[3]>=0. && newcorners[4]>=0. &&
+                  newcorners[3]+newcorners[4]<=1.)
+                {
+                  printf ("    Adding new point %g,%g,%g\n",
+                          newcorners[3], newcorners[4], newcorners[5]);
+                  if (numnewpoints >= newpointsize)
+                    newpoints = realloc
+                      (newpoints, (newpointsize+=100)*sizeof(ternary_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]) +
@@ -387 +419 @@
                   (newcorners[7]-newcorners[4])*(newcorners[7]-newcorners[4]) +
                   (newcorners[8]-newcorners[5])*(newcorners[8]-newcorners[5]) >
-                  vertex_tolerance * vertex_tolerance)
-                printf ("    Adding new point %g,%g,%g; result %d\n",
-                        newcorners[6], newcorners[7], newcorners[8],
-                        qh_addpoint ((pointT *) newcorners+6, facet, True));
-              */
+                  vertex_tolerance * vertex_tolerance &&
+                  newcorners[6]>=0. && newcorners[7]>=0. &&
+                  newcorners[6]+newcorners[7]<=1.)
+                {
+                  printf ("    Adding new point %g,%g,%g\n",
+                          newcorners[6], newcorners[7], newcorners[8]);
+                  if (numnewpoints >= newpointsize)
+                    newpoints = realloc
+                      (newpoints, (newpointsize+=100)*sizeof(ternary_point));
+                  newpoints [numnewpoints].C2 = newcorners[6];
+                  newpoints [numnewpoints].C3 = newcorners[7];
+                  newpoints [numnewpoints].G  = newcorners[8];
+                  newpoints [numnewpoints].efunc = eparam[1];
+                  numnewpoints++;
+                }
             }
         }
@@ -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");
-
-  return 0;
+  free(newpoints);
+
+  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=50, *verts, hullnumverts, *hullverts=NULL;
+  int i,j, index, loop_max=50, numpoints, *verts, hullnumverts,*hullverts=NULL;
   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 */
-  if (!(points = calloc ((loop_max+1)*(loop_max+2), sizeof (ternary_point))))
+  if (!(points = malloc (numpoints * sizeof (ternary_point))))
     { printf ("Cannot allocate point coordinate array\n"); exit (1); }
-  if (!(verts = calloc (loop_max*loop_max * 6, sizeof (int))))
+  if (!(verts = malloc (loop_max*loop_max * 6 * sizeof (int))))
     { 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); }
-  if (i=init_triangle_array (loop_max, points+(loop_max+1)*(loop_max+2)/2))
+  if (i=init_triangle_array (loop_max, points+numpoints/2))
     { 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)/2))
+                                numpoints/2))
     { printf ("main: Error %d in init_triangle_vertices\n", i); exit (i); }
 
   /* Calculate free energies */
-  if (i=free_energies (points, (loop_max+1)*(loop_max+2)/2, T,P, allparams,0,0))
+  if (i=free_energies (points, numpoints/2, T,P, allparams,0,0))
     { printf ("main: Error %d in free_energies\n", i); exit (i); }
-  if (i=free_energies (points+(loop_max+1)*(loop_max+2)/2,
-                       (loop_max+1)*(loop_max+2)/2, T,P, allparams,1,0))
+  if (i=free_energies (points+numpoints/2, numpoints/2, T,P, allparams,1,0))
     { printf ("main: Error %d in free_energies\n", i); exit (i); }
 
   /* Scale all free energy values */
-  if (i=scale_energy_array ((loop_max+1)*(loop_max+2), points,
-                            0., -.2, 1., .9, NULL, NULL))
+  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); }
 
@@ -81 +80 @@
   if (i=GeomviewDisplayTriangleCOFF
       (pfd, "Ternary Free Energy", "tfe", "shading smooth",
-       (loop_max+1)*(loop_max+2), points, loop_max*loop_max*2, verts))
+       numpoints, points, loop_max*loop_max*2, verts))
     { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
 
   /* Calculate and refine the lower convex hull */
   printf ("qhull version: %s\n", hullQHullVersion());
-  if (i=hullCalculate (3, points, (loop_max+1)*(loop_max+2)))
+  if (i=hullCalculate (3, points, numpoints))
     { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
-  if (i=hullRefine (points, (loop_max+1)*(loop_max+2), 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); }
@@ -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), points, hullnumverts, hullverts))
+       numpoints, points, hullnumverts, hullverts))
     { 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
-(ternary_point *points, int numpoints, energy_params *eparams, int nparams,
+(ternary_point **points, int *numpoints, energy_params *eparams, int nparams,
  double T,double P, double newton_tolerance, double vertex_tolerance,
  int one_phase_refine);