Changeset 366

Show
Ignore:
Timestamp:
04/18/08 10:37:21 (7 months ago)
Author:
powell
Message:

New hull calculation interface and small implementation changes.

Location:
trunk/matml/src/ternary
Files:
4 modified

Legend:

Unmodified
Added
Removed

  • trunk/matml/src/ternary/ChangeLog

    r364 r366  
    33  * Overhauled the free energy interface and functions, adding derivative 
    44    implementations. 
     5  * Changed hull calculation interface to permit (multiple) refinement between 
     6    original calculation and return. 
    57 
    68 -- 

  • trunk/matml/src/ternary/qhull.c

    r338 r366  
    1515 
    1616/*++++++++++++++++++++++++++++++++++++++ 
    17   This calculates the convex hull of a set of ternary compositions, clipping 
    18   everything above the plane with the three corners, and removing all 
    19   non-simplical facets.  It should return only the facets making up the minimum 
    20   energy surface. 
     17  This initializes qhull and uses it to calculate the convex hull of a set of 
     18  ternary compositions, clipping everything above the energy defined by the 
     19  three corners.  It does not return anything, but leaves all of the facets in 
     20  the qhull structures. 
    2121 
    22   int qhullCalcHull It returns zero or an error code. 
     22  int hullCalculate It returns zero or an error code. 
    2323 
    2424  int dim Number of dimensions (should always be 3 for ternary). 
    … …  
    2727 
    2828  ternary_point *points Point information. 
    29  
    30   int *numfacets Pointer which contains the number of facets in the convex hull 
    31   on return. 
    32  
    33   int **facetverts Pointer which contains the array of facet vertices on 
    34   return.  This is changed by realloc() each time the function is called, so it 
    35   should be NULL when first called. 
    36  
    37   char **version Pointer which contains the pointer to qh_version on return, or 
    38   NULL to ignore this. 
    3929  ++++++++++++++++++++++++++++++++++++++*/ 
    4030 
    41 int qhullCalcHull (int dim, int numpoints, ternary_point *points, 
    42                    int *numfacets, int **facetverts, char **version) 
     31int hullCalculate (int dim, int numpoints, ternary_point *points) 
    4332{ 
    4433  int i, corner1, corner2, corner3; 
    4534  coordT *qpoints; 
    46   facetT *facet; 
    4735 
    4836  if (dim != 3) 
    49     { printf ("qhullCalcHull: Non-3-D spaces not supported\n"); return -1; } 
     37    { printf ("qhullCalcHull: Non-3-D spaces not yet supported\n"); return -1; } 
    5038 
    51   /*+ This first re-creates the points array for three reasons: the coordT type 
    52     may or may not be the same as double, the qhull array needs to have just 
    53     three entries per point vs. 10+ for ternary_point, and this provides the 
    54     opportunity to clip or otherwise transform the array. +*/ 
     39  /*+ This first copies the points array for qhull, and clips the array at the 
     40    plane defined by the corners. +*/ 
    5541  if (!(qpoints = malloc (dim*numpoints*sizeof(coordT)))) 
    5642    { printf ("qhullCalcHull: could not allocate memory for points\n"); 
    … …  
    5844  for (i=0; i<numpoints; i++) 
    5945    { 
     46      // Copy compositions into the new array 
    6047      qpoints[3*i]   = (coordT) points[i].C2; 
    6148      qpoints[3*i+1] = (coordT) points[i].C3; 
    6249 
     50      // Find the corners 
    6351      if (points[i].C2 == 0. && points[i].C3 == 0.) 
    6452        corner1 = i; 
    … …  
    6957    } 
    7058 
    71   /* Clip the free energy at the plane given by the corners */ 
     59  // Copy free energy into the new array, clipping it at the plane given by the 
     60  // corners 
    7261  for (i=0; i<numpoints; i++) 
    7362    qpoints[3*i+2] = (coordT) 
    … …  
    7665            (points[corner3].G-points[corner1].G) * points[i].C2); 
    7766 
     67  // Initialize qhull and calculate the hull 
    7868  qh_init_A (stdin, stdout, stderr, 0, NULL); 
    7969  qh_init_B (qpoints, numpoints, dim, False); 
    … …  
    8272  qh_printsummary (stdout); 
    8373#endif 
     74 
     75  return 0; 
     76} 
     77 
     78 
     79/*++++++++++++++++++++++++++++++++++++++ 
     80  This returns the facets of the convex hull comprising the minimum energy 
     81  surface.  It removes the sides and top, currently by getting rid of 
     82  non-simplical facets (sides) and a facet whose projection area is the entire 
     83  triangle (top). 
     84 
     85  int hullReturnFacets It returns zero or an error code. 
     86 
     87  int *numfacets Pointer to an integer which will contain the number of facets 
     88  in the convex hull on return. 
     89 
     90  int **facetverts Pointer which contains the array of facet vertices on 
     91  return.  This is changed by realloc() each time the function is called, so it 
     92  should be NULL or a valid array. 
     93  ++++++++++++++++++++++++++++++++++++++*/ 
     94 
     95int hullReturnFacets (int *numfacets, int **facetverts) 
     96{ 
     97  int i=0; 
     98  facetT *facet; 
    8499 
    85100  /*+ After qhull runs, this reallocates the array pointed to by facetverts, 
    … …  
    95110      return -1; } 
    96111#ifdef DEBUG 
    97   printf ("post-realloc 0x%lx\n", *facetverts); 
     112  printf ("  post-realloc 0x%lx\n", *facetverts); 
    98113#endif 
    99114 
    100   i=0; 
    101115  FORALLfacets 
    102116    { 
    103117      int j=0; 
     118      coordT corners [6], projarea; 
    104119      vertexT *vertex, **vertexp; 
    105120 
    … …  
    113128 
    114129          if (j<3) 
    115             (*facetverts) [3*i+j] = thepoint; 
     130            { 
     131              (*facetverts) [3*i+j] = thepoint; 
     132              corners [2*j]   = vertex->point[0]; 
     133              corners [2*j+1] = vertex->point[1]; 
    116134#ifdef DEBUG 
    117           printf ("  vertex %d, coords %g %g %g\n", thepoint, 
    118                   vertex->point[0], vertex->point[1], vertex->point[2]); 
     135              printf ("  vertex %d, coords %g %g %g\n", thepoint, 
     136                      vertex->point[0], vertex->point[1], vertex->point[2]); 
    119137#endif 
    120           j++; 
     138              j++; 
     139            } 
    121140        } 
     141           
     142      /*+ This calculates the area of the projection of each facet. 
     143        + That area is given by the determinant of the matrix formed by the 
     144        + vectors making up the edges of the facet divided by 
     145        + (dimensionality-1)!. 
     146        +latex+For a ternary system, that's 
     147        +latex+$\frac{(x_1-x_0)(y_2-y_0) - (x_2-x_0)(y_1-y_0)}{2}$. 
     148        +*/ 
     149      projarea = 
     150        (corners[2]-corners[0]) * (corners[5]-corners[1]) - 
     151        (corners[4]-corners[0]) * (corners[3]-corners[1]); 
     152#ifdef DEBUG 
     153      printf ("  xyarea %g\n", projarea); 
     154#endif 
    122155 
    123       /* Remove non-simplical facets */ 
    124       if (j>3) 
    125         { 
    126           (*numfacets)--; 
    127         } 
    128  
    129       /* Remove facet with the three corners using an area test */ 
     156      // Remove side facets (zero area) and facet with the three corners 
     157      if (projarea == 0. || fabs (projarea) == 1.) 
     158        (*numfacets)--; 
    130159      else 
    131         { 
    132           coordT x0 = qpoints [3* (*facetverts) [3*i]], 
    133             y0 = qpoints [3* (*facetverts) [3*i]+1], 
    134             x1 = qpoints [3* (*facetverts) [3*i+1]], 
    135             y1 = qpoints [3* (*facetverts) [3*i+1]+1], 
    136             x2 = qpoints [3* (*facetverts) [3*i+2]], 
    137             y2 = qpoints [3* (*facetverts) [3*i+2]+1]; 
    138 #ifdef DEBUG 
    139           printf ("  xyarea %g\n", 
    140                   fabs ((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0))); 
    141 #endif 
    142           if (fabs ((x1-x0)*(y2-y0) - (x2-x0)*(y1-y0)) == 1.) 
    143             (*numfacets)--; 
    144  
    145           else 
    146             i++; 
    147         } 
     160        i++; 
    148161    } 
    149162 
    … …  
    152165#endif 
    153166 
    154   if (version) 
    155     *version = qh_version; 
    156  
    157   free (qpoints); 
    158   qh_freeqhull (True); 
     167  /* free (qpoints); 
     168     qh_freeqhull (True); */ 
    159169 
    160170  return 0; 
    161171} 
     172 
     173 
     174/*++++++++++++++++++++++++++++++++++++++ 
     175  This is what it sounds like. 
     176 
     177  char *hullQHullVersion It returns the QHull version string. 
     178  ++++++++++++++++++++++++++++++++++++++*/ 
     179 
     180char *hullQHullVersion () 
     181{ 
     182  return qh_version; 
     183} 

  • trunk/matml/src/ternary/ternary.c

    r361 r366  
    2929{ 
    3030  int i, j, index, loop_max=50, *verts, hullnumverts, *hullverts=NULL; 
    31   char gv_version[100], *qh_version; 
     31  char gv_version[100]; 
    3232  FILE *pfd = NULL; 
    3333  double T=0.75; 
    … …  
    8484 
    8585  /* Calculate the lower convex hull */ 
    86   if (i=qhullCalcHull (3, (loop_max+1)*(loop_max+2), points, &hullnumverts, 
    87                        &hullverts, &qh_version)) 
    88     { printf ("main: qhullCalcHull returned %d\n", i); exit (i); } 
    89   printf ("qhull version: %s\n", qh_version); 
     86  printf ("qhull version: %s\n", hullQHullVersion()); 
     87  if (i=hullCalculate (3, (loop_max+1)*(loop_max+2), points)) 
     88    { printf ("main: Error %d in hullCalculate\n", i); exit (i); } 
     89  if (i=hullReturnFacets (&hullnumverts, &hullverts)) 
     90    { printf ("main: Error %d in hullReturnFacets\n", i); exit (i); } 
    9091 
    9192  /* Take out hull triangles which are on the main function (in one phase) */ 

  • trunk/matml/src/ternary/ternary.h

    r360 r366  
    104104 
    105105/* From qhull.c */ 
    106 int qhullCalcHull (int dim, int numpoints, ternary_point *points, int *facets, 
    107                    int **verts, char **version); 
     106int hullCalculate (int dim, int numpoints, ternary_point *points); 
     107int hullReturnFacets (int *facets, int **verts); 
     108char *hullQHullVersion (); 
    108109 
    109110/* From book.c */