Changeset 415 for trunk/matml/src

Timestamp:

02/22/2009 11:02:56 AM (3 years ago)

Author:

powell

Message:

New facet_type enum and function for better calculation of facet type.

Location:

trunk/matml/src/ternary

Files:

• qhull.c (modified) (12 diffs)
• ternary.h (modified) (4 diffs)

trunk/matml/src/ternary/qhull.c

@@ -76 +76 @@
 
   return 0;
+}
+
+
+/*++++++++++++++++++++++++++++++++++++++
+  This determines the facet type from the minimum energy at the midpoint of
+  each edge.
+
+  facet_type facettype It returns the facet type.
+
+  coordT *corners Coordinates C2, C3, G of the three facet corners.
+
+  energy_params *eparams Collection of energy parameter structures.
+
+  int nparams Number of energy parameter structures in eparams.
+  ++++++++++++++++++++++++++++++++++++++*/
+
+static facet_type facettype (coordT *corners, energy_params *eparams,
+                             int nparams, double T, double P)
+{
+  int i, cparam, ret=0;
+  double edgenergy, edgecenter[3];
+
+  edgecenter[0] = (corners[0] + corners[3]) / 2.;
+  edgecenter[1] = (corners[1] + corners[4]) / 2.;
+  edgecenter[2] = (corners[2] + corners[5]) / 2.;
+  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
+  cparam=0;
+  for (i=1; i<nparams; i++)
+    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
+      {
+        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
+        cparam = i;
+      }
+  if (edgenergy > edgecenter[2])
+    ret++;
+      
+  edgecenter[0] = (corners[0] + corners[6]) / 2.;
+  edgecenter[1] = (corners[1] + corners[7]) / 2.;
+  edgecenter[2] = (corners[2] + corners[8]) / 2.;
+  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
+  cparam=0;
+  for (i=1; i<nparams; i++)
+    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
+      {
+        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
+        cparam = i;
+      }
+  if (edgenergy > edgecenter[2])
+    ret++;
+      
+  edgecenter[0] = (corners[3] + corners[6]) / 2.;
+  edgecenter[1] = (corners[4] + corners[7]) / 2.;
+  edgecenter[2] = (corners[5] + corners[8]) / 2.;
+  edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams);
+  cparam=0;
+  for (i=1; i<nparams; i++)
+    if (edgenergy > free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i))
+      {
+        edgenergy = free_energy (edgecenter[0],edgecenter[1], T,P, eparams+i);
+        cparam = i;
+      }
+  if (edgenergy > edgecenter[2])
+    ret++;
+
+  return (facet_type) ret;
 }
 
@@ -214 +279 @@
 {
   facetT *facet;
+  facet_type type;
   ternary_point *newpoints;
   int newpointsize=100, numnewpoints=0, i;
@@ -262 +328 @@
         {
           int eparamc=0;
+          facet_type facetype;
 
 #ifdef DEBUG
@@ -277 +344 @@
             +latex+\item
             +html+ <li>
-            If the lowest free energy function at the centroid has lower energy
-            than the average energy of the vertices, and the three vertices are
-            on the same free energy curve, add it as a new vertex.
+            If this is a one-phase facet, optionally add the centroid as a new
+            vertex.
             +html+ </li>
             +*/
@@ -285 +351 @@
           centroid[1] = (corners[1] + corners[4] + corners[7]) / 3.;
           centroid[2] = (corners[2] + corners[5] + corners[8]) / 3.;
-          
-          Gcenter = free_energy (centroid[0], centroid[1], T, P, eparams);
-          for (j=1; j<nparams; j++)
-            if (Gcenter > free_energy (centroid[0],centroid[1], T,P, eparams+j))
-              {
-                Gcenter = free_energy (centroid[0],centroid[1], T,P, eparams+j);
-                eparamc = j;
-              }
-
-#ifdef DEBUG
-          printf ("  Centroid %g, free energy %g (func %d)\n",
-            centroid[2], Gcenter, eparamc);
-#endif
-
-          if (Gcenter < centroid[2] &&
+
+          facetype = facettype (corners, eparams, nparams, T, P);
+
+          if (facetype == ONE_PHASE &&
               eparam[0] == eparam[1] && eparam[1] == eparam[2])
             {
-              if (one_phase_refine || eparam[0] != eparamc)
+              if (one_phase_refine)
                 {
 #ifdef DEBUG
@@ -522 +577 @@
   int hullReturnFacets It returns zero or an error code.
 
+  energy_params *eparams Collection of energy parameter structures.
+
+  int nparams Number of energy parameter structures in eparams.
+
+  double T Temperature.
+
+  double P Pressure.
+
   int *numfacets Pointer to an integer which will contain the number of facets
   in the convex hull on return.
@@ -528 +591 @@
   return.  This is changed by realloc() each time the function is called, so it
   should be NULL or a valid array.
+
+  facet_type **facettypes Pointer which contains the array of facet types on
+  return.  This is changed by realloc() each time the function is called, so it
+  should be NULL or a valid array.
   ++++++++++++++++++++++++++++++++++++++*/
 
-int hullReturnFacets (int *numfacets, int **facetverts)
+int hullReturnFacets
+(energy_params *eparams, int nparams, double T,double P,
+ int *numfacets, int **facetverts, facet_type **facettypes)
 {
   int i=0;
@@ -544 +613 @@
 #endif
   if (!(*facetverts = realloc (*facetverts, *numfacets * 3 * sizeof (int))))
-    { printf ("qhullCalcHull: could not reallocate memory for vertices\n");
+    { printf ("hullReturnFacets: could not reallocate memory for vertices\n");
       return -1; }
 #ifdef DEBUG
   printf ("  post-realloc 0x%lx\n", *facetverts);
+  printf ("facettypes: pre-realloc 0x%lx, ", *facettypes);
+#endif
+  if (!(*facettypes = realloc (*facettypes, *numfacets * sizeof (facet_type))))
+    { printf ("hullReturnFacets: could not reallocate memory for facet types\n");
+      return -1; }
+#ifdef DEBUG
+  printf ("  post-realloc 0x%lx\n", *facettypes);
 #endif
 
@@ -553 +629 @@
     {
       int j=0;
-      coordT corners [6], projarea;
+      coordT corners [9], projarea;
       vertexT *vertex, **vertexp;
 
@@ -567 +643 @@
             {
               (*facetverts) [3*i+j] = thepoint;
-              corners [2*j]   = vertex->point[0];
-              corners [2*j+1] = vertex->point[1];
+              corners [3*j]   = vertex->point[0];
+              corners [3*j+1] = vertex->point[1];
+              corners [3*j+2] = vertex->point[2];
 #ifdef DEBUG
               printf ("  vertex %d, coords %g %g %g\n", thepoint,
@@ -585 +662 @@
         +*/
       projarea =
-        (corners[2]-corners[0]) * (corners[5]-corners[1]) -
-        (corners[4]-corners[0]) * (corners[3]-corners[1]);
+        (corners[3]-corners[0]) * (corners[7]-corners[1]) -
+        (corners[6]-corners[0]) * (corners[4]-corners[1]);
 #ifdef DEBUG
       printf ("  xyarea %g\n", projarea);
@@ -595 +672 @@
         (*numfacets)--;
       else
-        i++;
+        {
+          (*facettypes) [i] = facettype (corners, eparams, nparams, T, P);
+          printf ("  Type %d\n", (*facettypes) [i]);
+          i++;
+        }
     }
 

trunk/matml/src/ternary/ternary.h

@@ -97 +97 @@
 
 typedef enum {
+  STANDARD=0,
+  WITH_DERIVATIVES=1,
+  NO_INFINITY=2,
+  WITH_DERIVATIVES_NO_INFINITY=3
+} free_energy_flags;
+
+typedef enum {
   ONE_PHASE_EDGE=0, /*+ One phase's boundary edges +*/
   TIE_SIMPLICES,    /*+ Tie lines (ternary) or other simplices between two
@@ -103 +110 @@
 } phase_boundary_flags;
 
-typedef enum {
-  STANDARD=0,
-  WITH_DERIVATIVES=1,
-  NO_INFINITY=2,
-  WITH_DERIVATIVES_NO_INFINITY=3
-} free_energy_flags;
-
 typedef struct {
   int compos;   /*+ Number of components of the space (ternary binary etc.) +*/
   int *edges;   /*+ Edges with indices of ternary_points comprising a phase
-                  boundary; each edge has n_phases points +*/
+                  boundary; each edge has compos points +*/
   int n_edges;  /*+ Number of edges comprising this boundary +*/
   int *phases;  /*+ Phases involved (energy_params indices) starting with the
@@ -125 +125 @@
 #define TERNARY_INFINITY HUGE_VALF
 #define TERNARY_NEGATIVE_INFINITY (-HUGE_VALF)
+
+typedef enum {
+  ONE_PHASE=0,
+  CRITICAL_POINT,
+  TWO_PHASE,
+  THREE_PHASE
+} facet_type;
 
 /* From freenergy.c */
@@ -152 +159 @@
  double T,double P, double newton_tolerance, double vertex_tolerance,
  int one_phase_refine);
-int hullReturnFacets (int *facets, int **verts);
+int hullReturnFacets
+(energy_params *eparams, int nparams, double T,double P,
+ int *numfacets, int **facetverts, facet_type **facettypes);
 char *hullQHullVersion ();