Changeset 373

Timestamp:

05/23/08 05:59:49 (6 months ago)

Author:

powell

Message:

New phase boundary, doesn't quite build.

Location:

trunk/matml/src/ternary

Files:

• ChangeLog (modified) (1 diff)
• qhull.c (modified) (5 diffs)
• ternary.h (modified) (1 diff)

trunk/matml/src/ternary/ChangeLog

@@ -3 +3 @@
   * Overhauled the free energy interface and functions, adding derivative
     implementations.
+  * Added phase_boundary typedef struct.
   * Changed hull calculation interface to permit (multiple) refinement between
     original calculation and return.

trunk/matml/src/ternary/qhull.c

@@ -105 +105 @@
 ++++++++++++++++++++++++++++++++++++++*/
 
-static inline void NewtonRefine
+static void NewtonRefine
 (coordT *point, coordT *corners, energy_params *eparams, int nparams,
  double T, double P, int globaldims, int localdims, int side,
@@ -140 +140 @@
   else
     {
-      /*+ Slope for the binary case +*/
+      /* 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]);
@@ -176 +176 @@
 
 /*++++++++++++++++++++++++++++++++++++++
-  This refines each of the facets in the hull, adding a new vertex at its
-  centroid if in a one-phase region, and using Newton-Raphson iteration to
-  refine down the vertices of multi-phase facets.
+  This refines each of the facets in the hull, using Newton-Raphson iteration
+  to refine down the vertices of multi-phase facets, and (optionally) adding a
+  new vertex at its centroid if in a one-phase region.
 
   int hullRefine It returns zero or an error code.
@@ -197 +197 @@
   choose this parameter to indicate how close to consider them the "same"
   point.
+
+  int one_phase_refine Non-zero if the one-phase regions should be refined by
+  adding a new vertex at the centroid.
   ++++++++++++++++++++++++++++++++++++++*/
 
 int hullRefine (energy_params *eparams, int nparams, double T, double P,
-                double newton_tolerance, double vertex_tolerance)
+                double newton_tolerance, double vertex_tolerance,
+                int one_phase_refine)
 {
   facetT *facet;
@@ -214 +218 @@
   FORALLfacets
     {
-      coordT corners [6], centroid [2];
+      double corners [6], energies [3], centroid [3], Gcenter;
       vertexT *vertex, **vertexp;
+      int j=0;
 
       /*+
         +latex+\item
         +html+ <li>
-        If the lowest free energy function at the centroid has lower energy
-        than the average energy of the vertices, add it as a new vertex.
+        Copy vertex information into local variables.
         +html+ </li>
         +*/
-      if (centroid[0] < corners[0])
-        ;
-
-      /*+
-        +latex+\item
-        +html+ <li>
-        Otherwise, use Newton-Raphson iteration to refine each vertex.
-        +html+ </li>
-        +*/
-      else
+      FOREACHvertex_(facet->vertices)
         {
-
+          int thepoint = (vertex->point-qh first_point)/3;
+
+          if (j<3)
+            {
+              (*facetverts) [3*i+j] = thepoint;
+              corners [2*j]   = vertex->point[0];
+              corners [2*j+1] = vertex->point[1];
+              energies [j]    = vertex->point[2];
+            }
+          j++;
+        }
+
+      if (j<3) /* Ignore non-simplical facets */
+        {
           /*+
             +latex+\item
             +html+ <li>
-            Find the lowest free energy function at this vertex.
+            If the lowest free energy function at the centroid has lower energy
+            than the average energy of the vertices, add it as a new vertex.
             +html+ </li>
             +*/
+          centroid[0] = (corners[0] + corners[2] + corners[4]) / 3.;
+          centroid[1] = (corners[1] + corners[3] + corners[5]) / 3.;
+          centroid[2] = (energies[0] + energies[1] + energies[2]) / 3.;
+          /*** CHECK ALL EPARAMS OF CORNERS ***/
+          Gcenter = free_energy (centroid[0], centroid[1], T, P, eparams);
+
+          if (Gcenter < centroid[2] && one_phase_refine)
+              ;
 
           /*+
             +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 along
-            each side that it's on, to avoid the "infinite slope at the edge"
-            problem.  If it's on more than one side, add a new candidate vertex
-            for each side it's on.
+            Otherwise, use Newton-Raphson iteration to refine each vertex.
             +html+ </li>
             +*/
+          else
+            {
+              /*+
+                +latex+\item
+                +html+ <li>
+                Find the lowest free energy function at this vertex.
+                +html+ </li>
+                +*/
+
+              /*+
+                +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 along
+                each side that it's on, to avoid the "infinite slope at the
+                edge" problem.  If it's on more than one side, add a new
+                candidate vertex for each side it's on.
+                +html+ </li>
+                +*/
+
+              /*+
+                +latex+\item
+                +html+ <li>
+                Otherwise refine starting at this vertex.
+                +html+ </li>
+                +*/
+            }
 
           /*+
             +latex+\item
             +html+ <li>
-            Otherwise refine starting at this vertex.
+            If any two new candidate vertices are in the same place (closer than
+            vertex_tolerance), eliminate one.
             +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>
-        +*/
     }
   /*+

trunk/matml/src/ternary/ternary.h

@@ -88 +88 @@
 } energy_params;   /*+ Free energy parameters structure +*/
 
+typedef struct {
+  int compos;   /*+ Number of components of the space +*/
+  int *edges;   /*+ Edges with indices of ternary_points comprising a phase
+                  boundary; each edge has compos-n_phases+1 points +*/
+  int n_edges;  /*+ Number of points comprising this boundary +*/
+  int *phases;  /*+ Phases involved (energy_params indices) starting with the
+                  "current" phase +*/
+  int n_phases; /*+ Number of phases involved +*/
+  int ties;     /*+ Zero if this is one phase's boundary edges (the "current"
+                  phase), non-zero if this is a set of tie simplicies +*/
+} phase_boundary; /*+ Phase boundary structure which can store either the
+                    boundary between a single phase and a multi-phase region,
+                    or the tie lines/simplices across a multi-phase region +*/
+
 /* From freenergy.c */
 double free_energy