Changeset 424 for trunk

Timestamp:

02/23/2009 07:27:26 PM (3 years ago)

Author:

powell

Message:

Fixed hullRefine.

Location:

trunk/matml/src/ternary

Files:

• qhull.c (modified) (13 diffs)
• ternary.c (modified) (1 diff)
• ternary.h (modified) (1 diff)

trunk/matml/src/ternary/qhull.c

@@ -366 +366 @@
   int hullRefine It returns zero or an error code.
 
-  ternary_point **points Pointer to point list.
-
-  int *numpoints Number of points of which to calculate the convex hull.
+  ternary_point **points Pointer to point list, which will be modified to
+  accommodate the new points created in the refining process.
+
+  int *numpoints Pointer to number of points of which to calculate the convex
+  hull, which will also be modified.
+
+  hull_facet **facets Pointer to the array of facets to be refined, which will
+  be modified to accommodate new facets created -- and some removed -- in the
+  refining process.
+
+  int *numfacets Pointer to the number of facets comprising the convex hull,
+  which will also be modified.
 
   energy_params *eparams Collection of energy parameter structures.
@@ -390 +399 @@
 
 int hullRefine
-(ternary_point **points, int *numpoints, energy_params *eparams, int nparams,
- double T,double P, double newton_tolerance, double vertex_tolerance,
- int one_phase_refine)
+(ternary_point **points, int *numpoints, hull_facet **facets, int *numfacets,
+ energy_params *eparams, int nparams, double T,double P,
+ double newton_tolerance, double vertex_tolerance, int one_phase_refine)
 {
-  facetT *facet;
-  facet_type type;
   ternary_point *newpoints;
   int newpointsize=100, numnewpoints=0, i;
@@ -410 +417 @@
     +*/
 
-  FORALLfacets
+  for (i=0; i<*numfacets; i++)
     {
       double corners [9], centroid [3], Gcenter;
-      vertexT *vertex, **vertexp;
       int j=0, eparam [3] = {0,0,0}, thepoint [3];
 
@@ -423 +429 @@
         +html+ </li>
         +*/
-      FOREACHvertex_(facet->vertices)
+      for (j=0; j<3; j++)
         {
-          if (j<3)
-            {
-              thepoint [j] = (vertex->point-qh first_point)/3;
-
-              corners [3*j]   = vertex->point[0];
-              corners [3*j+1] = vertex->point[1];
-              corners [3*j+2] = vertex->point[2];
-
-              eparam [j] = (*points)[thepoint[j]].efunc;
-            }
-          j++;
+          thepoint [j] = (*facets) [i].vertex[j];
+
+          corners [3*j]   = (*points) [thepoint[j]].C2;
+          corners [3*j+1] = (*points) [thepoint[j]].C3;
+          corners [3*j+2] = (*points) [thepoint[j]].G;
+          eparam [j] = (*points)[thepoint[j]].efunc;
         }
 
-      // Ignore non-simplical side facets and top facet
-      if (fabs((corners[3]-corners[0])*(corners[7]-corners[1]) -
-               (corners[6]-corners[0])*(corners[4]-corners[1])) < 0.999999
-          && j==3)
-        {
           int eparamc=0;
-          facet_type facetype;
-
-#ifdef DEBUG
-          printf ("Facet 0x%lx: %g,%g %g,%g %g,%g\n", facet,
+
+#ifdef DEBUG
+          printf ("Facet %d: %g,%g %g,%g %g,%g\n", i,
                   corners[0], corners[1], corners[3],
                   corners[4], corners[6], corners[7]);
@@ -468 +463 @@
           centroid[2] = (corners[2] + corners[5] + corners[8]) / 3.;
 
-          facetype = facettype (corners, eparams, nparams, T, P);
-
-          if (facetype == ONE_PHASE &&
-              eparam[0] == eparam[1] && eparam[1] == eparam[2])
+          if ((*facets) [i].type == ONE_PHASE && one_phase_refine)
             {
-              if (one_phase_refine)
-                {
-#ifdef DEBUG
-                  printf ("  Refining %s-phase triangle at centroid\n",
-                    (eparam[0]==eparamc) ? "one" : "multi");
+#ifdef DEBUG
+                  printf ("  Refining one-phase triangle at centroid\n");
 #endif
                   if (numnewpoints >= newpointsize)
@@ -484 +473 @@
                   newpoints [numnewpoints].C2 = centroid[0];
                   newpoints [numnewpoints].C3 = centroid[1];
-                  newpoints [numnewpoints].G  = Gcenter;
-                  newpoints [numnewpoints].efunc = eparamc;
+                  newpoints [numnewpoints].G  = free_energy
+                    (centroid[0], centroid[1], T, P, eparams+eparam [0]);
+                  newpoints [numnewpoints].efunc = eparam [0];
                   numnewpoints++;
-                }
-
-              else
-                {
-#ifdef DEBUG
-                  printf ("  Not refining one-phase triangle\n");
-#endif
-                }
             }
 
@@ -504 +486 @@
             +html+ </li>
             +*/
-          else
+          else if ((*facets) [i].type != ONE_PHASE)
             {
               coordT newcorners [9] = { corners[0], corners[1], corners[2],
@@ -541 +523 @@
                 }
 
+              /*+
+                +latex+\item
+                +html+ <li>
+                Refine the vertices of the facet using
+                +latex+{\tt NewtonRefine()} (section
+                +latex+\ref{func_NewtonRefine_qhull.c} page
+                +latex+\pageref{func_NewtonRefine_qhull.c}).
+                +html+ <a href="#func-NewtonRefine"><tt>NewtonRefine()</tt></a>.</li>
+                +*/
               for (j=0; j<3; j++)
                 {
@@ -561 +552 @@
                 +latex+\item
                 +html+ <li>
-                If any two new candidate vertices are in the same place (closer
-                than vertex_tolerance), eliminate one.
+                Add new points, leaving out duplicates in same place (closer
+                than vertex_tolerance).
                 +html+ </li>
                 +*/
@@ -630 +621 @@
                 }
             }
-        }
     }
   /*+
@@ -637 +627 @@
     +*/
 
-  /*+ Unfortunately, adding new points can create dissonance between qhull
-    indices and Ternary indices.  So this frees the old hull and creates a new
-    one around the old points and new refined ones.  This adds an N log N
-    process to an otherwise order-N refining algorithm, but can't be
-    helped. +*/
 #ifdef DEBUG
   printf ("Adding %d new points to hull\n", numnewpoints);
@@ -650 +635 @@
     return 1;
           
-  //qh_freeqhull (True);
-  //free (qpoints);
-
   for (i=0; i<numnewpoints; i++)
     {
@@ -662 +644 @@
 
       (*points) [(*numpoints)+i] = newpoints [i];
-      // In case I try to implement this using qh_addpoint someday...
-      /*printf ("  Finding best facet for coords %g, %g, %g\n",
-              (*points) [(*numpoints)+i].C2,(*points) [(*numpoints)+i].C3,
-              (*points) [(*numpoints)+i].G);
-      facet = qh_findbestfacet (newcorners, True, &bestdist, &isoutside);
-      printf ("  Adding to %d facet 0x%lx: %g, %g, %g\n", isoutside, facet, newcorners[0], newcorners[1],
-              newcorners[2]);
-      if (isoutside)
-        printf ("  Result: %d\n", qh_addpoint (newcorners, facet, True));
-      else
-      (*numpoints)--;*/
     }
   *numpoints += numnewpoints;
-
-#ifdef DEBUG
-  printf ("Done, returning\n");
-#endif
   free(newpoints);
 
-  //return hullCalculate (3, *points, *numpoints);
+#ifdef DEBUG
+  printf ("Re-calculating hull using hullCalculate()\n");
+#endif
+
+  /*+ Finally, after generating the new points and adding them to the hull, it
+    re-calculates the convex hull using
+    +latex+{\tt hullCalculate()} (section
+    +latex+\ref{func_hullCalculate_qhull.c} page
+    +latex+\pageref{func_hullCalculate_qhull.c}).
+    +html+ <a href="#func-hullCalculate"><tt>hullCalculate()</tt></a>.</li>
+    +*/
+  
+  return hullCalculate (3, *points,*numpoints, eparams,nparams, T,P, facets,
+                        numfacets);
 }
 

trunk/matml/src/ternary/ternary.c

@@ -101 +101 @@
                        &hullfacets, &hullnumfacets))
     { printf ("main: Error %d in hullCalculate\n", i); exit (i); }
-  /*if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
-                    1e-10, 1e-7, 0))
+  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
+                    allparams, 2, T, P, 1e-10, 1e-7, 0))
     { printf ("main: error %d in hullRefine\n", i); exit (i); }
-  if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
-                    1e-10, 1e-7, 0))
+  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
+                    allparams, 2, T, P, 1e-10, 1e-7, 0))
     { printf ("main: error %d in hullRefine\n", i); exit (i); }
-  if (i=hullRefine (&points, &numpoints, allparams, 2, T, P,
-                    1e-10, 1e-7, 0))
-    { printf ("main: error %d in hullRefine\n", i); exit (i); }*/
+  if (i=hullRefine (&points, &numpoints, &hullfacets, &hullnumfacets,
+                    allparams, 2, T, P, 1e-10, 1e-7, 0))
+    { printf ("main: error %d in hullRefine\n", i); exit (i); }
 
+  /* Scale and display everything */
   if (i=scale_energy_array (numpoints, points, 0., -.2, 1., .9, NULL, NULL))
     { printf ("main: Error %d in scale_triangle_array\n", i); exit (i); }

trunk/matml/src/ternary/ternary.h

@@ -171 +171 @@
  int nparams, double T, double P, hull_facet **facets, int *numfacets);
 int hullRefine
-(ternary_point **points, int *numpoints, energy_params *eparams, int nparams,
- double T,double P, double newton_tolerance, double vertex_tolerance,
- int one_phase_refine);
+(ternary_point **points, int *numpoints, hull_facet **facets, int *numfacets,
+ energy_params *eparams, int nparams, double T,double P,
+ double newton_tolerance, double vertex_tolerance, int one_phase_refine);
 char *hullQHullVersion ();