Changeset 378

Timestamp:

06/26/08 10:52:01 (4 months ago)

Author:

powell

Message:

Add pressure to ternary.c; try to use hullRefine (segfaults for now).

Files:

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

trunk/matml/src/ternary/qhull.c

@@ -80 +80 @@
-  This uses Newton iteration to refine a single 
-
-double
+coordT
-  result will be returned here.
-
-double
+coordT
-  facet, used for the facet normal: C21,C31,G1, C22,C32,G2, C23,C33,G3.
-
-  energy_params *eparams Collection of energy parameter structures.
-
-  int nparams Number of energy parameter structures in eparams.
-
-  double T Temperature.
-
@@ -97 +95 @@
-
-  int localdims Dimensionality over which to do the refinement.
-
-  int side Side of the phase diagram on which to refine, if
-  localdims<globaldims.
-
-  double newton_tolerance Stop Newton iterations when the vertex motion falls
@@ -106 +101 @@
-
-static void NewtonRefine
-
-
-
+
+
-{
-, minfunc=0
+
-  double distance = 1., G, slope2, slope3;
-  ternary_point current;
@@ -119 +113 @@
-  current.C3 = point[1];
-
-  /*+ Determine which is the lowest energy function at this point +*/
-  G = free_energy (current.C2, current.C3, T, P, eparams);
-  for (i=1; i<nparams; i++)
-    if (slope2 = free_energy (current.C2, current.C3, T, P, eparams+i) < G)
-      {
-        G = slope2;
-        minfunc = i;
-      }
-
-  /*+ Determine the slope(s) of the facet +*/
@@ -135 +122 @@
-        xG=corners[5]-corners[2], yG=corners[8]-corners[2];
-      double vG=x2*y3-x3*y2;
-3*yG-xG*y3
+G*y3-x3*yG
-      slope3 = (xG*y2-x2*yG)/vG;
+
+      printf ("      corners: %g %g %g %g %g %g %g %g %g\n      slopes: %g %g\n",
+              corners[0], corners[1], corners[2],
+              corners[3], corners[4], corners[5],
+              corners[6], corners[7], corners[8], slope2, slope3);
-    }
-  else
-    {
-
+
+
-      if (side==0 || side==2) // C3=0 or C2+C3=1 sides: dG/dC2
-        slope2 = (corners[5]-corners[2])/(corners[3]-corners[0]);
-      else if (side==1) // C2=0 side: dG/dC3
-        slope2 = (corners[5]-corners[2])/(corners[4]-corners[1]);
+      */
-    }
-
-
+ * newton_tolerance
-    {
-      double dC2, dC3;
-
-      /*+ Calculate the energy derivatives +*/
-+minfunc
+
-
-      /*+ Subtract the facet slopes from the energy derivatives +*/
--
--
++
++
-
-      /*+ Solve the linear system to estimate the vector to the minimum +*/
@@ -167 +161 @@
-      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);
-    }
-
@@ -218 +214 @@
-  FORALLfacets
-    {
-6], energies [3
+9
-      vertexT *vertex, **vertexp;
-
+, eparam [3] = {0,0,0}
-
-      /*+
-        +latex+\item
-        +html+ <li>
-
+
+
-        +html+ </li>
-        +*/
@@ -234 +231 @@
-          if (j<3)
-            {
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            }
-          j++;
-        }
-
-
+
+
+
+
-        {
+          int eparamc=0;
+          realT bestdist;
+          boolT isoutside;
+
+          printf ("Triangle: %g,%g %g,%g %g,%g\n",
+                  corners[0], corners[1], corners[3],
+                  corners[4], corners[6], corners[7]);
+          printf ("  Curves: %d %d %d, Area: %g\n",
+                  eparam[0], eparam[1], eparam[2],
+                  fabs((corners[3]-corners[0])*(corners[7]-corners[1]) -
+                       (corners[6]-corners[0])*(corners[4]-corners[1])));
+
-          /*+
-            +latex+\item
-            +html+ <li>
-            If the lowest free energy function at the centroid has lower energy
-
+
+
-            +html+ </li>
-            +*/
-2] + corners[4
-3] + corners[5
-energies[0] + energies[1] + energies[2
-/*** CHECK ALL EPARAMS OF CORNERS ***/
+3] + corners[6
+4] + corners[7
+corners[2] + corners[5] + corners[8
+
-          Gcenter = free_energy (centroid[0], centroid[1], T, P, eparams);
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
-          /*+
-            +latex+\item
-            +html+ <li>
-
+
+
-            +html+ </li>
-            +*/
-          else
-            {
+              coordT newcorners [9] = { corners[0], corners[1], corners[2],
+                                        corners[3], corners[4], corners[5],
+                                        corners[6], corners[7], corners[8]};
+
+              printf ("  Refining multi-phase triangle\n");
+              for (j=0; j<3; j++)
+                {
+                  /*+
+                    +latex+\item
+                    +html+ <li>
+                    If the vertex is on a side of the phase diagram (one or
+                    more composition variables is zero or they sum to one),
+                    refine generally starting a bit closer to the others.
+                    +html+ </li>
+                    +*/
+                  if (corners[3*j] == 0. || corners[3*j+1] == 0 ||
+                      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];
+                    }
+                  else
+                    {
+                  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],
+                                T,P, 3,3, newton_tolerance);
+                  printf ("    -> %g, %g\n", newcorners[3*j],
+                          newcorners[3*j+1]);
+
+                  newcorners[3*j+2] = free_energy
+                    (newcorners[3*j],newcorners[3*j+1], T,P, eparams+eparam[j]);
+                    }
+                }
+
-              /*+
-                +latex+\item
-                +html+ <li>
-
+
+
-                +html+ </li>
-                +*/
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-            }
-
-          /*+
-            +latex+\item
-            +html+ <li>
-            If any two new candidate vertices are in the same place (closer than
-            vertex_tolerance), eliminate one.
-            +html+ </li>
-            +*/
-        }
-    }
@@ -306 +411 @@
-    +html+ </ol>
-    +*/
+
+  printf ("Done, returning\n");
-
-  return 0;

trunk/matml/src/ternary/ternary.c

@@ -31 +31 @@
-  char gv_version[100];
-  FILE *pfd = NULL;
-
+, P=1.
-  ternary_point *points;
-  /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */
-  energy_params /* Solid, liquid */
-    eparams1 = {1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 1.,2.,5., -.5,2.,.2, 0.},
-
+
+
-
-  if (argc>1)
@@ -66 +67 @@
-
-  /* Calculate free energies */
-1.
+P
-    { printf ("main: Error %d in free_energies\n", i); exit (i); }
-  if (i=free_energies (points+(loop_max+1)*(loop_max+2)/2,
-1.
+P
-    { 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,
-  
+
-    { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
-
@@ -87 +88 @@
-  if (i=hullCalculate (3, (loop_max+1)*(loop_max+2), points))
-    { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
+  if (i=hullRefine (allparams, 2, T, P, 1e-15, 1e-10, 0))
+    { printf ("main: error %d in hullRefine\n", i); exit (i); }
-  if (i=hullReturnFacets (&hullnumverts, &hullverts))
-    { printf ("main: Error %d in hullReturnFacets\n", i); exit (i); }
@@ -110 +113 @@
-      /* Test whether hull facet is in one-phase region: free energy at
-         centroid is below facet; if so, remove it */
-1.
-1.
+P
+P
-        {
-          for (j=i; j<hullnumverts-1; j++)

trunk/matml/src/ternary/ternary.h

@@ -119 +119 @@
-/* From qhull.c */
-int hullCalculate (int dim, int numpoints, ternary_point *points);
+int hullRefine (energy_params *eparams, int nparams, double T, double P,
+                double newton_tolerance, double vertex_tolerance,
+                int one_phase_refine);
-int hullReturnFacets (int *facets, int **verts);
-char *hullQHullVersion ();