Changeset 369

Timestamp:

05/09/08 12:58:25 (6 months ago)

Author:

powell

Message:

New code for doing the Newton refining; completely untested.

Files:

• trunk/matml/src/ternary/qhull.c (modified) (4 diffs)

trunk/matml/src/ternary/qhull.c

@@ -84 +84 @@
 
   double *corners Coordinates (including energies) of the corners of this
-  facet, used for the facet normal.
+  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.
+
+  double P Pressure.
 
   int globaldims Dimensionality of the global space.
@@ -103 +107 @@
 static inline void NewtonRefine
 (coordT *point, coordT *corners, energy_params *eparams, int nparams,
- int globaldims, int localdims, int side, double newton_tolerance)
-{
-  double distance = 0.;
+ double T, double P, int globaldims, int localdims, int side,
+ double newton_tolerance)
+{
+  int i, minfunc=0;
+  double distance = 1., G, slope2, slope3;
+  ternary_point current;
+
+  if (globaldims != 3)
+    { printf ("qhullCalcHull: Non-ternaries not yet supported\n"); return; }
+  current.C2 = point[0];
+  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 +*/
+  if (localdims==3)
+    {
+      double x2=corners[3]-corners[0], y2=corners[6]-corners[0],
+        x3=corners[4]-corners[1], y3=corners[7]-corners[1],
+        xG=corners[5]-corners[2], yG=corners[8]-corners[2];
+      double vG=x2*y3-x3*y2;
+      slope2 = (x3*yG-xG*y3)/vG;
+      slope3 = (xG*y2-x2*yG)/vG;
+    }
+  else
+    {
+      /*+ Slope for the binary case +*/
+      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]);
+    }
 
   while (distance >= newton_tolerance)
     {
+      double dC2, dC3;
+
       /*+ Calculate the energy derivatives +*/
+      free_energy_derivatives (&current, 1, T, P, eparams+minfunc);
 
       /*+ Subtract the facet slopes from the energy derivatives +*/
-
-      /*+ Solve the linear system to estimate the distance to the minimum +*/
-
-      /*+ Add that distance to the point +*/
-    }
+      current.G2 -= slope2;
+      current.G3 -= slope3;
+
+      /*+ Solve the linear system to estimate the vector to the minimum +*/
+      dC3 = current.G22*current.G33 - current.G23*current.G23;
+      dC2 = (current.G3*current.G23 - current.G2*current.G33) / dC3;
+      dC3 = (current.G2*current.G23 - current.G3*current.G22) / dC3;
+
+      /*+ Add that vector to the point +*/
+      current.C2 += dC2;
+      current.C3 += dC3;
+      distance = dC2*dC2 + dC3*dC3;
+    }
+
+  /* Return the current point */
+  point[0] = current.C2;
+  point[1] = current.C3;
 }
 
@@ -132 +185 @@
 
   int nparams Number of energy parameter structures in eparams.
+
+  double T Temperature.
+
+  double P Pressure.
 
   double newton_tolerance Stop Newton iterations when the vertex motion falls
@@ -142 +199 @@
   ++++++++++++++++++++++++++++++++++++++*/
 
-int hullRefine (energy_params *eparams, int nparams, double newton_tolerance,
-                double vertex_tolerance)
+int hullRefine (energy_params *eparams, int nparams, double T, double P,
+                double newton_tolerance, double vertex_tolerance)
 {
   facetT *facet;