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

Revision 471, 6.2 kB (checked in by powell, 3 years ago)

Opening 0.5 with new arbitrary polynomials in free energy function.
  • Property svn:eol-style set to native
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, res2=20,res3=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  /* eparams: R,T0, G1@T0,G2,G3, C1,C2,C3, M1,M2,M3, O12,O13,O23,O123 */
39  energy_gaussian species = { "Species 0.8 0.8", .8, .8, -.2, .1, 2. };
40  energy_params /* Solid, liquid */
41    eparams0 = {"Solid", 1.,1., 0.,-.1,-.2, -1.,-1.1,-1.2, 0.,1.,1.,1.,1.,
42                -.5,2.,.2, 0.,-.1,-.2,.1,.2,-.1, NULL,0, &species,1},
43    eparams1 = {"Liquid", 1.,1., .3,.1,-.1, -2.2,-2.,-1.8, 0.,1.,.67,1.,2.,
44                -.3,-.4,-.1, 0.,-.1,-.15,-.05,0.,.3, NULL,0, NULL,0},
45      allparams[2] = {eparams0, eparams1};
46
47  if (argc>1)
48    sscanf (argv[1], "%d", &res2);
49  if (argc>2)
50    sscanf (argv[2], "%d", &res3);
51  else
52    res3 = res2;
53  numpoints= 2 * (res2+1) * (res3+1);
54
55  /* Allocate array pairs for ternary points and triangle vertices */
56  if (!(points = (energy_point *) malloc (numpoints * sizeof (energy_point))))
57    { printf ("Cannot allocate point coordinate array\n"); exit (1); }
58  if (!(verts = (int *) malloc (res2*res3 * 12 * sizeof (int))))
59    { printf ("Cannot allocate triangle vertex array\n"); exit (1); }
60
61  /* Start Geomview process */
62  if (i=GeomviewBegin (&pfd, gv_version))
63    { printf ("main: Error %d in GeomviewBegin\n", i); exit (i); }
64  printf("Geomview version: %s\n", gv_version);
65
66  /* Initialize coordinates of both rectangle array halves */
67  if (i=init_rectangle_array (res2,res3, points, 0, 0.,0., 1.,1.))
68    { printf ("main: Error %d in init_rectangle_array\n", i); exit (i); }
69  if (i=init_rectangle_array (res2,res3, points+numpoints/2, 1, 0.,0., 1.,1.))
70    { printf ("main: Error %d in init_rectangle_array\n", i); exit (i); }
71
72  /* Calculate triangle vertex indices for both arrays */
73  if (i=init_rectangle_vertices (res2,res3, verts, 0))
74    { printf ("main: Error %d in init_rectangle_vertices\n", i); exit (i); }
75  if (i=init_rectangle_vertices (res2,res3, verts + res2*res3*6,
76                                 numpoints/2))
77    { printf ("main: Error %d in init_rectangle_vertices\n", i); exit (i); }
78
79  /* Calculate free energies */
80  if (i=free_energies (points, numpoints/2, T,P, allparams,1,WITH_DERIVATIVES_NO_INFINITY))
81    { printf ("main: Error %d in free_energies\n", i); exit (i); }
82  if (i=free_energies (points+numpoints/2, numpoints/2, T,P, allparams,0,
83                       WITH_DERIVATIVES_NO_INFINITY))
84    { printf ("main: Error %d in free_energies\n", i); exit (i); }
85
86  /* Calculate the spinodal: first with no return, then more accurate with one
87     Problem: what's the spinodal criterion in a square array?
88  if (i=calc_spinodal (&points, &numpoints, verts, res2*res3*4, T, P,
89                       allparams, 0, &solid_spinodal))
90    { printf ("main: error %d in spinodal calculation\n", i); exit (i); }
91  */
92
93  /* Scale all free energy values */
94  if (i=energy_to_visual_array (numpoints, points, &tpoints, 0., -.2, 1., .9,
95                                NULL, NULL, 1))
96    { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }
97
98  /* Send free energy surface data to Geomview */
99  if (i=GeomviewDisplayTriangleCOFF
100      (pfd, "Ternary Free Energy", "tfe", "shading smooth",
101       numpoints, tpoints, res2*res3*4, verts))
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, 1))
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.