root/trunk/matml/src/ternary/ternary.c

Revision 475, 6.3 kB (checked in by powell, 3 years ago)

Polynomial coefficients are now temperature-dependent.
Line 
1/***************************************
2  $Header$
3
4  This little program uses Geomview to display the free energy function on a
5  ternary system in three dimensions.
6  ***************************************/
7
8
9#include "gibbs.h" /*+ Ternary prototypes, typedefs, etc. +*/
10#include <stdlib.h>  /*+ For calloc() +*/
11#include <math.h>    /*+ For fabs() +*/
12
13
14/*++++++++++++++++++++++++++++++++++++++
15  This program starts Geomview, displays a single ternary phase diagram, waits
16  for user interaction, and closes.  Its only argument right now is an optional
17  "resolution" of the discretization: "ternary 30" will use a 30x30/2 triangle
18  array to display free energy and calculate the phase diagram.
19
20  int main It returns zero or an error code.
21
22  int argc Number of arguments.
23
24  char *argv[] Arguments: right now only the "resolution" of the
25  discretization.
26  ++++++++++++++++++++++++++++++++++++++*/
27
28int main (int argc, char *argv[])
29{
30  int i,j, index, loop_max=20, numpoints, *verts, hullnumfacets, numbounds;
31  hull_facet *hullfacets = NULL;
32  char gv_version[100];
33  FILE *pfd = NULL;
34  double T=0.7, P=1.;
35  energy_point *points;
36  visual_point *tpoints=NULL;
37  phase_boundary solid_spinodal, *allbounds=NULL;
38  energy_polynomial squarer = { -2., 0.5, 1., 3., 2. };
39  /* eparams: name, R,T0, G1@T0,G2,G3, C1,C2,C3, S1,S2,S3,S4,S5,
40     O12,O13,O23, O123,O234,O235,O245,O345,O2345 gauss ptr, gauss number */
41  energy_params /* Solid, liquid */
42    eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 1.,.5,.2,0.,0.,
43                -.5,2.,.2, 0.,0.,0.,0.,0.,0., &squarer,1, NULL,0},
44    eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 1.,1.,.67,0.,0.,
45                -.3,-.4,-.1, -.5,0.,0.,0.,0.,0., NULL,0, NULL,0},
46      allparams[2] = {eparams0, eparams1};
47
48  if (argc>1)
49    sscanf (argv[1], "%d", &loop_max);
50  numpoints=(loop_max+1)*(loop_max+2);
51
52  /* Allocate array pairs for ternary points and triangle vertices */
53  if (!(points = (energy_point *) malloc (numpoints * sizeof (energy_point))))
54    { printf ("Cannot allocate point coordinate array\n"); exit (1); }
55  if (!(verts = (int *) malloc (loop_max*loop_max * 6 * sizeof (int))))
56    { printf ("Cannot allocate triangle vertex array\n"); exit (1); }
57
58  /* Start Geomview process */
59  if (i=GeomviewBegin (&pfd, gv_version))
60    { printf ("main: Error %d in GeomviewBegin\n", i); exit (i); }
61  printf("Geomview version: %s\n", gv_version);
62
63  /* Initialize coordinates of both triangle array halves */
64  if (i=init_triangle_array (loop_max, points, 0, 0.,0., 1.,0., 0.,1.))
65    { printf ("main: Error %d in init_triangle_array\n", i); exit (i); }
66  if (i=init_triangle_array (loop_max, points+numpoints/2, 1, 0.,0., 1.,0., 0.,1.))
67    { printf ("main: Error %d in init_triangle_array\n", i); exit (i); }
68
69  /* Calculate triangle vertex indices for both arrays */
70  if (i=init_triangle_vertices (loop_max, verts, 0))
71    { printf ("main: Error %d in init_triangle_vertices\n", i); exit (i); }
72  if (i=init_triangle_vertices (loop_max, verts + loop_max*loop_max*3,
73                                numpoints/2))
74    { printf ("main: Error %d in init_triangle_vertices\n", i); exit (i); }
75
76  /* Calculate free energies */
77  if (i=free_energies (points, numpoints/2, T,P, allparams,1,WITH_DERIVATIVES_NO_INFINITY))
78    { printf ("main: Error %d in free_energies\n", i); exit (i); }
79  if (i=free_energies (points+numpoints/2, numpoints/2, T,P, allparams,0,
80                       WITH_DERIVATIVES_NO_INFINITY))
81    { printf ("main: Error %d in free_energies\n", i); exit (i); }
82
83  // Calculate the spinodal: first with no return, then more accurate with one
84  if (i=calc_spinodal (&points, &numpoints, verts, loop_max*loop_max*2, T, P,
85                       allparams, 0, &solid_spinodal))
86    { printf ("main: error %d in spinodal calculation\n", i); exit (i); }
87
88  /* Scale all free energy values */
89  if (i=energy_to_visual_array (numpoints, points, &tpoints, 0., -.2, 1., .9,
90                                NULL, NULL, 0))
91    { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
92
93  /* Send free energy surface data to Geomview */
94  if (i=GeomviewDisplayTriangleCOFF
95      (pfd, "Ternary Free Energy", "tfe", "shading smooth",
96       numpoints, tpoints, loop_max*loop_max*2, verts))
97    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
98
99  if (i=GeomviewDisplayPhaseBoundary
100      (pfd, "Spinodal", "spin", "-face +edge",
101       numpoints, tpoints, &solid_spinodal))
102    { printf ("main: Error %d in Geomview Display\n", i); exit (i); }
103
104  /* Calculate and refine the lower convex hull */
105  printf ("qhull version: %s\n", hullQHullVersion());
106  if (i=hullCalculate (3, points, numpoints, allparams, 2, T, P,
107                       &hullfacets, &hullnumfacets))
108    { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
109  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
110                    allparams, 2, T, P, 1e-10, 1e-7, 0))
111    { printf ("main: error %d in hullRefine\n", i); exit (i); }
112  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
113                    allparams, 2, T, P, 1e-10, 1e-7, 0))
114    { printf ("main: error %d in hullRefine\n", i); exit (i); }
115  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
116                    allparams, 2, T, P, 1e-10, 1e-7, 0))
117    { printf ("main: error %d in hullRefine\n", i); exit (i); }
118
119  /* Scale and display everything */
120  if (i=energy_to_visual_array (numpoints, points, &tpoints, 0., -.2, 1., .9,
121                                NULL, NULL, 0))
122    { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
123
124  if (i=hullReturnPhaseBoundaries (points,numpoints, hullfacets,hullnumfacets,
125                                   allparams, T, P, &allbounds, &numbounds))
126    { printf ("main: error %d in hullReturnPhaseBoundaries\n", i); exit (i); }
127
128  for (i=0; i<numbounds; i++)
129    {
130      char name [20], label [10];
131      int ret;
132
133      snprintf (name, 18, "Boundary %d", i);
134      snprintf (label, 8, "bd%d", i);
135
136      printf ("Boundary %d: %d-phase %s, phases:", i, (allbounds+i)->n_phases,
137              ((allbounds+i)->type==TIE_SIMPLICES) ? "tie-simplices":"edge");
138      for (ret=0; ret<(allbounds+i)->n_phases; ret++)
139        printf (" %d", (allbounds+i)->phases[ret]);
140      printf ("\n");
141
142      if (ret=GeomviewDisplayPhaseBoundary
143          (pfd, name, label, "-face +edge", numpoints, tpoints, allbounds+i))
144        { printf ("main: Error %d in Geomview Display\n", ret); exit (ret); }
145    }
146
147  printf ("Press <return> to close up... ");
148  fflush (stdout);
149  fgets (gv_version, 99, stdin);
150
151  GeomviewEnd (&pfd);
152  free (points);
153  free (verts);
154  free (hullfacets);
155
156  return 0;
157}
Note: See TracBrowser for help on using the browser.