Changeset 432 for trunk/matml

Timestamp:

03/06/2009 04:49:33 PM (3 years ago)

Author:

powell

Message:

Major overhaul of phase boundary return code including miscibility gap.

Location:

trunk/matml/src/ternary

Files:

• qhull.c (modified) (12 diffs)
• ternary.c (modified) (2 diffs)

trunk/matml/src/ternary/qhull.c

@@ -705 +705 @@
   int numfacets Number of facets.
 
+  energy_params *eparams Collection of energy parameter structures.
+
+  double T Temperature.
+
+  double P Pressure.
+
   phase_boundary **boundaries Pointer to array of phase_boundary objects.  This
   is changed by realloc() each time the function is called, so it should be
@@ -714 +720 @@
 int hullReturnPhaseBoundaries
 (ternary_point *points, int n_points, hull_facet *facets, int n_facets,
+ energy_params *eparams, double T, double P,
  phase_boundary **boundaries, int *n_bounds)
 {
@@ -731 +738 @@
         v2 = facets[f].vertex[2];
 
-      printf ("Last facet %d vertices %d, %d, %d\n", n_facets-1,
-              facets[n_facets-1].vertex[0], facets[n_facets-1].vertex[1],
-              facets[n_facets-1].vertex[2]);
-
       printf ("Facet %d: type %d, vertices %d ",f, facets [f].type, v0);
-      fflush (stdout);
       printf ("(%d), %d (%d), %d(%d)\n",
               points[v0].efunc, v1,
@@ -745 +747 @@
         {
         case ONE_PHASE:
-          {
-            printf ("Type 0 facet\n"); fflush (stdout);
-            // Totally uninteresting
-            break;
-          }
+          // Totally uninteresting
+          break;
 
         case THREE_PHASE:
@@ -795 +794 @@
                 (*boundaries) [bound].edges [2] = MAX3 (v0, v1, v2);
                 (*n_bounds) ++;
+                printf ("  Corners for bound %d edge 0: %d, %d, %d\n", bound,
+                        (*boundaries) [bound].edges [0],
+                        (*boundaries) [bound].edges [1],
+                        (*boundaries) [bound].edges [2]);
               }
             else
@@ -811 +814 @@
                 (*boundaries) [bound].edges [3*e+2] = MAX3 (v0, v1, v2);
                 (*boundaries) [bound].n_edges += 1;
+                printf ("  Corners for bound %d edge %d: %d, %d, %d\n", bound,e,
+                        (*boundaries) [bound].edges [3*e],
+                        (*boundaries) [bound].edges [3*e+1],
+                        (*boundaries) [bound].edges [3*e+2]);
               }
 
@@ -820 +827 @@
         case TWO_PHASE:
           {
-            int b=0, bound=-1, p0, p1, p2,
+            int b=0, bound=-1, p0, p1, p2, nphases=2, e=0,
               phase0 = MIN3 (points [v0].efunc, points [v1].efunc,
                              points [v2].efunc),
@@ -833 +840 @@
                 p1 = MAX (v0, v1);
                 p2 = v2;
-                // If all three are on the same phase, need to check midpoint
-                // energies to determine which are one- and two-phase segments.
-                // Punt for the default case later!
+
+                // If all three are on the same phase, this is a miscibility
+                // gap facet, set nphases to 1 and figure out which vertex is
+                // across from the other two
                 if (points [v0].efunc == points [v2].efunc)
-                  { printf ("Miscibility gap!\n"); break; }
+                  {
+                    double C02 = points [v0].C2, C03 = points [v0].C3,
+                      C12 = points [v1].C2, C13 = points [v1].C3,
+                      C22 = points [v2].C2, C23 = points [v2].C3,
+                      G0 = points [v0].G, G1 = points [v1].G,
+                      G2 = points [v2].G;
+
+                    printf ("Misc gap: ");
+                    nphases=1;
+
+                    if (free_energy ((C02+C12)/2,(C03,C13)/2,T,P,eparams+phase0) <
+                        G0+G1/2)
+                      {
+                        printf ("v2 alone, Gav=%g, Gmid=%g\n",
+                                G0+G1/2, free_energy ((C02+C12)/2,(C03,C13)/2,T,P,eparams+phase0));
+                        p0 = MIN (v0, v1);
+                        p1 = MAX (v0, v1);
+                        p2 = v2;
+                      }
+                    else if (free_energy ((C02+C22)/2,(C03,C23)/2,T,P,eparams+phase0) <
+                        G0+G2/2)
+                      {
+                        printf ("v1 alone, Gav=%g, Gmid=%g\n",
+                                G0+G2/2, free_energy ((C02+C22)/2,(C03,C23)/2,T,P,eparams+phase0));
+                        p0 = MIN (v0, v2);
+                        p1 = MAX (v0, v2);
+                        p2 = v1;
+                      }
+                    else
+                      {
+                        printf ("v0 alone, Gav=%g, Gmid=%g\n",
+                                G1+G2/2, free_energy ((C22+C12)/2,(C23,C13)/2,T,P,eparams+phase0));
+                        p0 = MIN (v1, v2);
+                        p1 = MAX (v1, v2);
+                        p2 = v0;
+                      }
+                  }
               }
             else if (points [v0].efunc == points [v2].efunc)
@@ -854 +898 @@
             // Check for a boundary with tie lines between these two phases
             for (b=0; b<*n_bounds; b++)
-              if ((*boundaries) [b].n_phases == 2 &&
+              if ((*boundaries) [b].n_phases == nphases &&
                   (*boundaries) [b].flags == TIE_SIMPLICES &&
                   (*boundaries) [b].phases[0] == phase0 &&
-                  (*boundaries) [b].phases[1] == phase1)
+                  (*boundaries) [b].phases[nphases-1] == phase1)
                 bound = b;
 
@@ -870 +914 @@
 
                 (*boundaries) [bound].compos = 2;
-                (*boundaries) [bound].n_phases = 2;
+                (*boundaries) [bound].n_phases = nphases;
                 if (!((*boundaries) [bound].phases =
-                      malloc (2*sizeof (int))))
+                      malloc (nphases*sizeof (int))))
                   return 1;
                 (*boundaries) [bound].phases [0] = phase0;
-                (*boundaries) [bound].phases [1] = phase1;
-                (*boundaries) [bound].n_edges = 2;
-                if (!((*boundaries) [bound].edges =
-                      malloc (4*sizeof (int))))
-                  return 1;
-                (*boundaries) [bound].edges [0] = MIN (p0, p2);
-                (*boundaries) [bound].edges [1] = MAX (p0, p2);
-                (*boundaries) [bound].edges [2] = MIN (p1, p2);
-                (*boundaries) [bound].edges [3] = MAX (p1, p2);
+                (*boundaries) [bound].phases [nphases-1] = phase1;
+                (*boundaries) [bound].n_edges = 0;
+                (*boundaries) [bound].edges = NULL;
                 (*boundaries) [bound].flags = TIE_SIMPLICES;
                 (*n_bounds) ++;
+                e=0;
               }
             else
-              {
-                int e = (*boundaries) [bound].n_edges;
-                printf ("Tie-line boundary %d adding edges %d,%d\n",bound,e,e+1);
-                // All of these reallocs are pretty expensive...
-                // Later use compos to store size of array, sort out later
-                if (!((*boundaries) [bound].edges = realloc
-                      ((*boundaries) [bound].edges, (2*e+4)*sizeof (int))))
-                  return 1;
-                (*boundaries) [bound].edges [2*e]   = MIN (p0, p2);
-                (*boundaries) [bound].edges [2*e+1] = MAX (p0, p2);
-                (*boundaries) [bound].edges [2*e+2] = MIN (p1, p2);
-                (*boundaries) [bound].edges [2*e+3] = MAX (p1, p2);
-                (*boundaries) [bound].n_edges += 2;
-              }
+              e = (*boundaries) [bound].n_edges;
+
+            // All of these reallocs are pretty expensive...
+            // Later use compos to store size of array, sort out later
+            if (!((*boundaries) [bound].edges = realloc
+                  ((*boundaries) [bound].edges, (2*e+4) * sizeof (int))))
+              return 1;
+            (*boundaries) [bound].n_edges += 2;
+
+            (*boundaries) [bound].edges [2*e]   = MIN (p0, p2);
+            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p2);
+            (*boundaries) [bound].edges [2*e+2] = MIN (p1, p2);
+            (*boundaries) [bound].edges [2*e+3] = MAX (p1, p2);
+            printf ("  Corners for bound %d edges %d,%d: %d,%d, %d,%d\n",
+                    bound, e, e+1, (*boundaries) [bound].edges [2*e],
+                    (*boundaries) [bound].edges [2*e+1],
+                    (*boundaries) [bound].edges [2*e+2],
+                    (*boundaries) [bound].edges [2*e+3]);
 
             // Check for a boundary with an edge on the "majority" phase
             for (b=0, bound=-1; b<*n_bounds; b++)
-              if ((*boundaries) [b].n_phases == 2 &&
+              if ((*boundaries) [b].n_phases == nphases &&
                   (*boundaries) [b].flags == ONE_PHASE_EDGE &&
                   (*boundaries) [b].phases[0] == points[p0].efunc &&
-                  (*boundaries) [b].phases[1] == points[p2].efunc)
+                  (*boundaries) [b].phases[nphases-1] == points[p2].efunc)
                 bound = b;
 
@@ -921 +964 @@
 
                 (*boundaries) [bound].compos = 2;
-                (*boundaries) [bound].n_phases = 2;
+                (*boundaries) [bound].n_phases = nphases;
                 if (!((*boundaries) [bound].phases =
-                      malloc (2*sizeof (int))))
+                      malloc (nphases*sizeof (int))))
                   return 1;
                 (*boundaries) [bound].phases [0] = points[p0].efunc;
-                (*boundaries) [bound].phases [1] = points[p2].efunc;
-                (*boundaries) [bound].n_edges = 1;
-                if (!((*boundaries) [bound].edges =
-                      malloc (2*sizeof (int))))
-                  return 1;
-                (*boundaries) [bound].edges [0] = MIN (p0, p1);
-                (*boundaries) [bound].edges [1] = MAX (p0, p1);
+                (*boundaries) [bound].phases [nphases-1] = points[p2].efunc;
+                (*boundaries) [bound].n_edges = 0;
+                (*boundaries) [bound].edges = NULL;
                 (*boundaries) [bound].flags = ONE_PHASE_EDGE;
                 (*n_bounds) ++;
+                e=0;
               }
             else
-              {
-                int e = (*boundaries) [bound].n_edges;
-                printf ("Phase edge boundary %d adding edge %d\n",bound,e);
-                // All of these reallocs are pretty expensive...
-                if (!((*boundaries) [bound].edges = realloc
-                      ((*boundaries) [bound].edges, (2*e+2)*sizeof (int))))
-                  return 1;
-                (*boundaries) [bound].edges [2*e]   = MIN (p0, p1);
-                (*boundaries) [bound].edges [2*e+1] = MAX (p0, p1);
-                (*boundaries) [bound].n_edges += 1;
-              }
-
+              e = (*boundaries) [bound].n_edges;
+
+            // All of these reallocs are pretty expensive...
+            if (!((*boundaries) [bound].edges = realloc
+                  ((*boundaries) [bound].edges, (2*e+2)*sizeof (int))))
+              return 1;
+            (*boundaries) [bound].edges [2*e]   = MIN (p0, p1);
+            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p1);
+            (*boundaries) [bound].n_edges += 1;
+            printf ("  Corners for bound %d edge %d: %d,%d\n",
+                    bound, e, (*boundaries) [bound].edges [2*e],
+                    (*boundaries) [bound].edges [2*e+1]);
+
+            break;
+          }
+
+        case CRITICAL_POINT:
+          {
+            printf ("Critical point case\n");
             break;
           }
@@ -962 +1009 @@
                     points [facets[f].vertex[1]].efunc,
                     points [facets[f].vertex[2]].efunc);
+            printf ("  Coords: %d %g,%g, %d %g,%g, %d %g,%g\n",
+                    facets[f].vertex[0], points [facets[f].vertex[0]].C2,
+                    points [facets[f].vertex[0]].C3,
+                    facets[f].vertex[1], points [facets[f].vertex[1]].C2,
+                    points [facets[f].vertex[1]].C3,
+                    facets[f].vertex[2], points [facets[f].vertex[2]].C2,
+                    points [facets[f].vertex[2]].C3);
           }
         }

trunk/matml/src/ternary/ternary.c

@@ -121 +121 @@
 
   if (i=hullReturnPhaseBoundaries (points,numpoints, hullfacets,hullnumfacets,
-                                   &allbounds, &numbounds))
+                                   allparams, T, P, &allbounds, &numbounds))
     { printf ("main: error %d in hullReturnPhaseBoundaries\n", i); exit (i); }
 
@@ -144 +144 @@
     { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
 
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 4", "bd4", "-face +edge",
+       numpoints, points, allbounds+4))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
+  if (i=GeomviewDisplayPhaseBoundary
+      (pfd, "Boundary 5", "bd5", "-face +edge",
+       numpoints, points, allbounds+5))
+    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }  
+
   printf ("Press <return> to close up... ");
   fflush (stdout);