Changeset 437 for trunk/matml

Timestamp:

03/16/2009 03:40:00 PM (3 years ago)

Author:

powell

Message:

Critical point implementation, completes this function.

Files:

• trunk/matml/src/ternary/qhull.c (modified) (1 diff)

trunk/matml/src/ternary/qhull.c

@@ -996 +996 @@
           {
             printf ("Critical point case\n");
+            int b=0, bound=-1, p0, p1, p2, e=0, phase = points [v0].efunc;
+            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;
+
+            if (free_energy ((C02+C12)/2,(C03+C13)/2,T,P,eparams+phase) >
+                (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+phase));
+                p0 = MIN (v0, v1);
+                p1 = MAX (v0, v1);
+                p2 = v2;
+              }
+            else if (free_energy ((C02+C22)/2,(C03+C23)/2,T,P,eparams+phase) >
+                     (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+phase));
+                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+phase));
+                p0 = MIN (v1, v2);
+                p1 = MAX (v1, v2);
+                p2 = v0;
+              }
+
+            // Check for a boundary with tie lines across this miscibility gap
+            for (b=0; b<*n_bounds; b++)
+              if ((*boundaries) [b].n_phases == 1 &&
+                  (*boundaries) [b].flags == TIE_SIMPLICES &&
+                  (*boundaries) [b].phases[0] == phase)
+                bound = b;
+
+            if (bound == -1)
+              {
+                bound = *n_bounds;
+                // No tie-line boundary in this phase, make a new one
+                printf ("Making new tie-line boundary %d\n", *n_bounds);
+                if (!(*boundaries = realloc
+                      (*boundaries, (*n_bounds+1) * sizeof (phase_boundary))))
+                  return 1;
+
+                (*boundaries) [bound].compos = 2;
+                (*boundaries) [bound].n_phases = 1;
+                if (!((*boundaries) [bound].phases = malloc (sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].phases [0] = phase;
+                (*boundaries) [bound].n_edges = 0;
+                (*boundaries) [bound].edges = NULL;
+                (*boundaries) [bound].flags = TIE_SIMPLICES;
+                (*n_bounds) ++;
+                e=0;
+              }
+            else
+              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+2) * sizeof (int))))
+              return 1;
+            (*boundaries) [bound].n_edges += 1;
+
+            (*boundaries) [bound].edges [2*e]   = MIN (p0, p1);
+            (*boundaries) [bound].edges [2*e+1] = MAX (p0, p1);
+            printf ("  Corners for bound %d edges %d,%d: %d,%d\n",
+                    bound, e, e+1, (*boundaries) [bound].edges [2*e],
+                    (*boundaries) [bound].edges [2*e+1]);
+
+            // Check for a boundary with an edge on this miscibility gap
+            for (b=0, bound=-1; b<*n_bounds; b++)
+              if ((*boundaries) [b].n_phases == 1 &&
+                  (*boundaries) [b].flags == ONE_PHASE_EDGE &&
+                  (*boundaries) [b].phases[0] == phase)
+                bound = b;
+
+            if (bound == -1)
+              {
+                // No edge boundary with this miscibility gap, make a new one
+                bound = *n_bounds;
+                printf ("Making new phase edge boundary %d\n", *n_bounds);
+                if (!(*boundaries = realloc
+                      (*boundaries, (bound+1) * sizeof (phase_boundary))))
+                  return 1;
+
+                (*boundaries) [bound].compos = 2;
+                (*boundaries) [bound].n_phases = 1;
+                if (!((*boundaries) [bound].phases = malloc (sizeof (int))))
+                  return 1;
+                (*boundaries) [bound].phases [0] = phase;
+                (*boundaries) [bound].n_edges = 0;
+                (*boundaries) [bound].edges = NULL;
+                (*boundaries) [bound].flags = ONE_PHASE_EDGE;
+                (*n_bounds) ++;
+                e=0;
+              }
+            else
+              e = (*boundaries) [bound].n_edges;
+
+            // All of these reallocs are pretty expensive...
+            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;
+            printf ("  Corners for bound %d edge %d: %d,%d, %d,%d\n",
+                    bound, e, (*boundaries) [bound].edges [2*e],
+                    (*boundaries) [bound].edges [2*e+1],
+                    (*boundaries) [bound].edges [2*e+2],
+                    (*boundaries) [bound].edges [2*e+3]);
+
             break;
           }