Index: ThirdParty/mpqc_open/src/lib/math/isosurf/LIBS.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/LIBS.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,5 @@
+libSCisosurf.LIBSUF
+#include <math/scmat/LIBS.h>
+#include <util/container/LIBS.h>
+#include <math/optimize/LIBS.h>
+#include <util/render/LIBS.h>
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/Makefile
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/Makefile	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,73 @@
+#
+# Makefile
+#
+# Copyright (C) 1996 Limit Point Systems, Inc.
+#
+# Author: Curtis Janssen <cljanss@ca.sandia.gov>
+# Maintainer: LPS
+#
+# This file is part of the SC Toolkit.
+#
+# The SC Toolkit is free software; you can redistribute it and/or modify
+# it under the terms of the GNU Library General Public License as published by
+# the Free Software Foundation; either version 2, or (at your option)
+# any later version.
+#
+# The SC Toolkit is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Library General Public License for more details.
+#
+# You should have received a copy of the GNU Library General Public License
+# along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+# the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+#
+# The U.S. Government is granted a limited license as per AL 91-7.
+#
+
+TOPDIR=../../../..
+ifndef SRCDIR
+  SRCDIR=$(shell pwd)
+endif
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalMakefile
+
+LD=$(CXX)
+
+CXXSRC = volume.cc shape.cc volume_i.cc \
+         surf.cc surfse.cc surfor.cc surfst.cc \
+         vertex.cc edge.cc triangle.cc isosurf.cc \
+         tricoef.cc tmplinst.cc
+CSRC = implicit.c
+
+LIBOBJ = $(CXXSRC:%.cc=%.$(OBJSUF)) $(CSRC:%.c=%.$(OBJSUF))
+
+INC = volume.h shape.h volume_i.h implicit.h tricoef.h
+
+DEPENDINCLUDE = $(INC)
+
+BIN_OR_LIB = LIB
+TARGET_TO_MAKE = libSCisosurf
+
+TESTPROGS = isotest
+
+DISTFILES = isotest.cc $(CXXSRC) $(CSRC) $(INC) Makefile LIBS.h
+
+default:: $(DEPENDINCLUDE)
+
+tmplinst.$(OBJSUF): tmplinst.cc
+	$(LTCOMP) $(CXX) $(CXXFLAGS:-fno-implicit-templates=) $(CPPFLAGS) -c $<
+
+isotest.$(OBJSUF): isotest.cc
+	$(LTCOMP) $(CXX) $(CXXFLAGS) $(CPPFLAGS) -DSRCDIR=\"$(SRCDIR)\" -c $<
+
+isotest: isotest.$(OBJSUF) \
+	$(shell $(LISTLIBS) $(INCLUDE) $(SRCDIR)/LIBS.h)
+	$(LTLINK) $(LD) $(LDFLAGS) -o isotest $^ $(SYSLIBS) $(LTLINKBINOPTS)
+
+include $(SRCDIR)/$(TOPDIR)/lib/GlobalRules
+
+$(LIBOBJ:.$(OBJSUF)=.d) isotest.d: $(DEPENDINCLUDE)
+ifneq ($(DODEPEND),no)
+include $(LIBOBJ:.$(OBJSUF)=.d) isotest.d
+endif
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/edge.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/edge.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/edge.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,175 @@
+//
+// edge.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <math/isosurf/triangle.h>
+#include <math/isosurf/edge.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+// Edge
+
+Edge::Edge(const Ref<Vertex> &p1,
+     const Ref<Vertex> &p2,
+     const Ref<Vertex> &p3)
+{
+  _order = 2;
+  _vertices = new Ref<Vertex>[3];
+  _vertices[0]=p1; _vertices[1]=p2; _vertices[2]=p3;
+}
+
+Edge::Edge(const Ref<Vertex> &p1,
+     const Ref<Vertex> &p2,
+     const Ref<Vertex> &p3,
+     const Ref<Vertex> &p4)
+{
+  _order = 3;
+  _vertices = new Ref<Vertex>[4];
+  _vertices[0]=p1; _vertices[1]=p2; _vertices[2]=p3; _vertices[3]=p4;
+}
+
+Edge::~Edge()
+{
+  delete[] _vertices;
+}
+
+double Edge::straight_length()
+{
+  SCVector3 BA = vertex(1)->point() - vertex(0)->point();
+  return BA.norm();
+}
+
+void Edge::add_vertices(std::set<Ref<Vertex> >&set)
+{
+  set.insert(_vertices[0]);
+  set.insert(_vertices[_order]);
+}
+
+void
+Edge::set_order(int order, const Ref<Volume>&vol,double isovalue)
+{
+  Ref<Vertex> *newvertices = new Ref<Vertex>[order+1];
+  newvertices[0] = vertex(0);
+  newvertices[order] = vertex(1);
+  delete[] _vertices;
+  _vertices = newvertices;
+  _order = order;
+
+  SCVector3 pv[2];
+  SCVector3 norm[2];
+
+  int i;
+  for (i=0; i<2; i++) {
+      pv[i] = vertex(i)->point();
+      norm[i] = vertex(i)->normal();
+    }
+
+  for (i=1; i<_order; i++) {
+      double I = (1.0*i)/order;
+      double J = (1.0*(_order - i))/order;
+      SCVector3 interpv = I * pv[0] + J * pv[1];
+      SCVector3 start(interpv);
+      SCVector3 interpnorm = I * norm[0] + J * norm[1];
+      SCVector3 newpoint;
+      vol->solve(start,interpnorm,isovalue,newpoint);
+      vol->set_x(newpoint);
+      if (vol->gradient_implemented()) {
+          vol->get_gradient(interpnorm);
+        }
+      interpnorm.normalize();
+      _vertices[i] = new Vertex(newpoint,interpnorm);
+    }
+
+}
+
+int
+Edge::interpolate(double r, SCVector3&point, SCVector3&norm)
+{
+  int i;
+  double s = 1.0 - r;
+  double rcoef[Triangle::max_order+1];
+  double scoef[Triangle::max_order+1];
+  double spacing = 1.0/_order;
+  rcoef[0] = 1.0;
+  scoef[0] = 1.0;
+  for (i=1; i<=_order; i++) {
+      rcoef[i] = rcoef[i-1]*(r-(i-1)*spacing)/(i*spacing);
+      scoef[i] = scoef[i-1]*(s-(i-1)*spacing)/(i*spacing);
+    }
+  int has_norm = 1;
+  for (i=0; i<=_order; i++) {
+      if ((has_norm = (has_norm && _vertices[i]->has_normal()))) break;
+    }
+  point = 0.0;
+  norm = 0.0;
+  for (i=0; i<=_order; i++) {
+      int j=_order-i;
+      double coef = rcoef[i]*scoef[j];
+      point += coef * _vertices[i]->point();
+      if (has_norm) norm += coef * _vertices[i]->normal();
+    }
+  if (has_norm) norm.normalize();
+  return has_norm;
+}
+
+int
+Edge::interpolate(double r, SCVector3&point, SCVector3&norm,
+                  const Ref<Volume> &vol, double isovalue)
+{
+  // first guess
+  int has_norm = interpolate(r,point,norm);
+
+  if (!has_norm) {
+      ExEnv::errn() << "Edge::interpolate with volume requires norm"
+           << endl;
+      abort();
+    }
+
+  // refine guess
+  SCVector3 newpoint;
+  vol->solve(point,norm,isovalue,newpoint);
+  // compute the true normal
+  vol->set_x(newpoint);
+  if (vol->gradient_implemented()) {
+      vol->get_gradient(norm);
+      norm.normalize();
+    }
+
+  return has_norm || vol->gradient_implemented();
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/edge.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/edge.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/edge.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,90 @@
+//
+// edge.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_edge_h
+#define _math_isosurf_edge_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <set>
+
+#include <math/isosurf/vertex.h>
+
+namespace sc {
+
+class Edge: public RefCount {
+  private:
+    int _order;
+    Ref<Vertex> *_vertices; // nvertices = _order + 1
+  public:
+    Edge(const Ref<Vertex> &p1,
+         const Ref<Vertex> &p2)
+    {
+      _order = 1;
+      _vertices = new Ref<Vertex>[2];
+      _vertices[0]=p1; _vertices[1]=p2;
+    }
+    Edge(const Ref<Vertex> &p1,
+         const Ref<Vertex> &p2,
+         const Ref<Vertex> &p3);
+    Edge(const Ref<Vertex> &p1,
+         const Ref<Vertex> &p2,
+         const Ref<Vertex> &p3,
+         const Ref<Vertex> &p4);
+    ~Edge();
+    int order() const { return _order; }
+    double straight_length();
+    // return the endpoints
+    Ref<Vertex> vertex(int i) const
+    {
+      return i?_vertices[_order]:_vertices[0];
+    }
+    // returns endpoints or interior vertex 0 <= i <= order
+    Ref<Vertex> interior_vertex(int i) const
+    {
+      return _vertices[i];
+    }
+    // add the endpoints to the set
+    void add_vertices(std::set<Ref<Vertex> >&);
+    void set_order(int order, const Ref<Volume>&vol,double isovalue);
+    // find the position of a point on the edge
+    int interpolate(double location, SCVector3&point, SCVector3&norm);
+    // find the true position of a point using the isosurface
+    int interpolate(double location, SCVector3&point, SCVector3&norm,
+                     const Ref<Volume> &vol, double isovalue);
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.c
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.c	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.c	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,1040 @@
+/*
+ * C code from the article
+ * "An Implicit Surface Polygonizer"
+ * by Jules Bloomenthal, jbloom@beauty.gmu.edu
+ * in "Graphics Gems IV", Academic Press, 1994
+ */
+
+/* Modified by Curtis Janssen:
+ *  1. Eliminate memory leaks.
+ *  2. Make main routine optional (-DMAIN to compile a main routine).
+ */
+
+/* implicit.c
+ *     an implicit surface polygonizer, translated from Mesa
+ *     applications should call polygonize()
+ *
+ * to compile a test program for ASCII output:
+ *     cc -DMAIN implicit.c -o implicit -lm
+ *
+ * to compile a test program for display on an SGI workstation:
+ *     cc -DMAIN -DSGIGFX implicit.c -o implicit -lgl_s -lm
+ *
+ * Authored by Jules Bloomenthal, Xerox PARC.
+ * Copyright (c) Xerox Corporation, 1991.  All rights reserved.
+ * Permission is granted to reproduce, use and distribute this code for
+ * any and all purposes, provided that this notice appears in all copies. */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <stdio.h>
+#include <sys/types.h>
+
+#define TET	0  /* use tetrahedral decomposition */
+#define NOTET	1  /* no tetrahedral decomposition  */
+
+#define RES	10 /* # converge iterations    */
+
+#define L	0  /* left direction:	-x, -i */
+#define R	1  /* right direction:	+x, +i */
+#define B	2  /* bottom direction: -y, -j */
+#define T	3  /* top direction:	+y, +j */
+#define N	4  /* near direction:	-z, -k */
+#define F	5  /* far direction:	+z, +k */
+#define LBN	0  /* left bottom near corner  */
+#define LBF	1  /* left bottom far corner   */
+#define LTN	2  /* left top near corner     */
+#define LTF	3  /* left top far corner      */
+#define RBN	4  /* right bottom near corner */
+#define RBF	5  /* right bottom far corner  */
+#define RTN	6  /* right top near corner    */
+#define RTF	7  /* right top far corner     */
+
+/* the LBN corner of cube (i, j, k), corresponds with location
+ * (start.x+(i-.5)*size, start.y+(j-.5)*size, start.z+(k-.5)*size) */
+
+#define RAND()	    ((rand()&32767)/32767.)    /* random number between 0 and 1 */
+#define HASHBIT	    (5)
+#define HASHSIZE    (size_t)(1<<(3*HASHBIT))   /* hash table size (32768) */
+#define MASK	    ((1<<HASHBIT)-1)
+#define HASH(i,j,k) ((((((i)&MASK)<<HASHBIT)|((j)&MASK))<<HASHBIT)|((k)&MASK))
+#define BIT(i, bit) (((i)>>(bit))&1)
+#define FLIP(i,bit) ((i)^1<<(bit)) /* flip the given bit of i */
+
+typedef struct point {		   /* a three-dimensional point */
+    double x, y, z;		   /* its coordinates */
+} POINT;
+
+typedef struct test {		   /* test the function for a signed value */
+    POINT p;			   /* location of test */
+    double value;		   /* function value at p */
+    int ok;			   /* if value is of correct sign */
+} TEST;
+
+typedef struct vertex {		   /* surface vertex */
+    POINT position, normal;	   /* position and surface normal */
+} VERTEX;
+
+typedef struct vertices {	   /* list of vertices in polygonization */
+    int count, max;		   /* # vertices, max # allowed */
+    VERTEX *ptr;		   /* dynamically allocated */
+} VERTICES;
+
+typedef struct corner {		   /* corner of a cube */
+    int i, j, k;		   /* (i, j, k) is index within lattice */
+    double x, y, z, value;	   /* location and function value */
+} CORNER;
+
+typedef struct cube {		   /* partitioning cell (cube) */
+    int i, j, k;		   /* lattice location of cube */
+    CORNER *corners[8];		   /* eight corners */
+} CUBE;
+
+typedef struct cubes {		   /* linked list of cubes acting as stack */
+    CUBE cube;			   /* a single cube */
+    struct cubes *next;		   /* remaining elements */
+} CUBES;
+
+typedef struct centerlist {	   /* list of cube locations */
+    int i, j, k;		   /* cube location */
+    struct centerlist *next;	   /* remaining elements */
+} CENTERLIST;
+
+typedef struct cornerlist {	   /* list of corners */
+    int i, j, k;		   /* corner id */
+    double value;		   /* corner value */
+    struct cornerlist *next;	   /* remaining elements */
+} CORNERLIST;
+
+typedef struct edgelist {	   /* list of edges */
+    int i1, j1, k1, i2, j2, k2;	   /* edge corner ids */
+    int vid;			   /* vertex id */
+    struct edgelist *next;	   /* remaining elements */
+} EDGELIST;
+
+typedef struct intlist {	   /* list of integers */
+    int i;			   /* an integer */
+    struct intlist *next;	   /* remaining elements */
+} INTLIST;
+
+typedef struct intlists {	   /* list of list of integers */
+    INTLIST *list;		   /* a list of integers */
+    struct intlists *next;	   /* remaining elements */
+} INTLISTS;
+
+typedef struct process {	   /* parameters, function, storage */
+    double (*function)();	   /* implicit surface function */
+    int (*triproc)();		   /* triangle output function */
+    double size, delta;		   /* cube size, normal delta */
+    int bounds;			   /* cube range within lattice */
+    POINT start;		   /* start point on surface */
+    CUBES *cubes;		   /* active cubes */
+    VERTICES vertices;		   /* surface vertices */
+    CENTERLIST **centers;	   /* cube center hash table */
+    CORNERLIST **corners;	   /* corner value hash table */
+    EDGELIST **edges;		   /* edge and vertex id hash table */
+} PROCESS;
+
+void *calloc();
+
+#define mycalloc(n,nbyte) _mycalloc(n,nbyte,__LINE__)
+#define myfree(ptr) _myfree(ptr,__LINE__)
+
+static void makecubetable ();
+static void free_cubetable();
+static void converge(POINT*,POINT*,double,double(*f)(),POINT*);
+static CORNER *setcorner (PROCESS*, int, int, int);
+static int setcenter(CENTERLIST *table[], int, int, int);
+static int dotet (CUBE*, int, int, int, int, PROCESS*);
+static int docube(CUBE*,PROCESS*);
+static void testface (int,int,int,CUBE*,int,int,int,int,int,PROCESS*);
+static TEST find (int,PROCESS*,double,double,double);
+static void vnormal (POINT*,PROCESS*,POINT*);
+static void addtovertices (VERTICES*, VERTEX);
+static int vertid (CORNER*,CORNER*,PROCESS*);
+static void free_process_data(PROCESS *);
+static void clean_malloc();
+static char *_mycalloc (int nitems, int nbytes, int line);
+static void _myfree(void*ptr, int lineno);
+
+#ifdef MAIN
+
+/**** A Test Program ****/
+
+
+/* ffunction: a piece of an atomic f function */
+
+double ffunction (x, y, z)
+double x, y, z;
+{
+  return x*y*z;
+}
+
+/* torus: a torus with major, minor radii = 0.5, 0.1, try size = .05 */
+
+double torus (x, y, z)
+double x, y, z;
+{
+    double x2 = x*x, y2 = y*y, z2 = z*z;
+    double a = x2+y2+z2+(0.5*0.5)-(0.1*0.1);
+    return a*a-4.0*(0.5*0.5)*(y2+z2);
+}
+
+
+/* sphere: an inverse square function (always positive) */
+
+double sphere (x, y, z)
+double x, y, z;
+{
+    double rsq = x*x+y*y+z*z;
+    return 1.0/(rsq < 0.00001? 0.00001 : rsq);
+}
+
+
+/* blob: a three-pole blend function, try size = .1 */
+
+double blob (x, y, z)
+double x, y, z;
+{
+    return 4.0-sphere(x+1.0,y,z)-sphere(x,y+1.0,z)-sphere(x,y,z+1.0);
+}
+
+#ifdef SGIGFX /**************************************************************/
+
+#include <math/isosurf/gl.h>
+
+/* triangle: called by polygonize() for each triangle; set SGI lines */
+
+triangle (i1, i2, i3, vertices)
+int i1, i2, i3;
+VERTICES vertices;
+{
+    float v[3];
+    int i, ids[3];
+    ids[0] = i1;
+    ids[1] = i2;
+    ids[2] = i3;
+    bgnclosedline();
+    for (i = 0; i < 3; i++) {
+	POINT *p = &vertices.ptr[ids[i]].position;
+	v[0] = p->x; v[1] = p->y; v[2] = p->z;
+	v3f(v);
+    }
+    endclosedline();
+    return 1;
+}
+
+
+/* main: call polygonize() with torus function
+ * display lines on SGI */
+
+main ()
+{
+    char *err, *polygonize();
+
+    keepaspect(1, 1);
+    winopen("implicit");
+    doublebuffer();
+    gconfig();
+    perspective(450, 1.0/1.0, 0.1, 10.0);
+    color(7);
+    clear();
+    swapbuffers();
+    makeobj(1);
+    if ((err = polygonize(torus, .1, 20, 0.,0.,0., triangle, TET)) != NULL) {
+	fprintf(stderr, "%s\n", err);
+	exit(1);
+    }
+    closeobj();
+    translate(0.0, 0.0, -2.0);
+    pushmatrix();
+    while(1) { /* spin the object */
+	reshapeviewport();
+	color(7);
+	clear();
+	color(0);
+	callobj(1);
+	rot(0.8, 'x');
+	rot(0.3, 'y');
+	rot(0.1, 'z');
+	swapbuffers();
+
+    }
+}
+
+#else /***********************************************************************/
+
+int gntris;	     /* global needed by application */
+VERTICES gvertices;  /* global needed by application */
+
+
+/* triangle: called by polygonize() for each triangle; write to stdout */
+
+triangle (i1, i2, i3, vertices)
+int i1, i2, i3;
+VERTICES vertices;
+{
+    gvertices = vertices;
+    gntris++;
+    fprintf(stdout, "%d %d %d\n", i1, i2, i3);
+    return 1;
+}
+
+
+/* main: call polygonize() with torus function
+ * write points-polygon formatted data to stdout */
+
+main ()
+    {
+    int i;
+    char *err, *polygonize();
+    gntris = 0;
+    fprintf(stdout, "triangles\n\n");
+    if ((err = polygonize(torus, .05, 20, 0.,0.,0., triangle, TET)) != NULL) {
+	fprintf(stdout, "%s\n", err);
+	exit(1);
+	}
+    fprintf(stdout, "\n%d triangles, %d vertices\n", gntris, gvertices.count);
+    fprintf(stdout, "\nvertices\n\n");
+    for (i = 0; i < gvertices.count; i++) {
+	VERTEX v;
+	v = gvertices.ptr[i];
+	fprintf(stdout, "%f  %f	 %f\t%f	 %f  %f\n",
+	    v.position.x, v.position.y,	 v.position.z,
+	    v.normal.x,	  v.normal.y,	 v.normal.z);
+    }
+    fprintf(stderr, "%d triangles, %d vertices\n", gntris, gvertices.count);
+    exit(0);
+}
+
+#endif /**********************************************************************/
+
+#endif /* MAIN */
+
+
+/**** An Implicit Surface Polygonizer ****/
+
+
+/* polygonize: polygonize the implicit surface function
+ *   arguments are:
+ *	 double function (x, y, z)
+ *		 double x, y, z (an arbitrary 3D point)
+ *	     the implicit surface function
+ *	     return negative for inside, positive for outside
+ *	 double size
+ *	     width of the partitioning cube
+ *	 int bounds
+ *	     max. range of cubes (+/- on the three axes) from first cube
+ *	 double x, y, z
+ *	     coordinates of a starting point on or near the surface
+ *	     may be defaulted to 0., 0., 0.
+ *	 int triproc (i1, i2, i3, vertices)
+ *		 int i1, i2, i3 (indices into the vertex array)
+ *		 VERTICES vertices (the vertex array, indexed from 0)
+ *	     called for each triangle
+ *	     the triangle coordinates are (for i = i1, i2, i3):
+ *		 vertices.ptr[i].position.x, .y, and .z
+ *	     vertices are ccw when viewed from the out (positive) side
+ *		 in a left-handed coordinate system
+ *	     vertex normals point outwards
+ *	     return 1 to continue, 0 to abort
+ *	 int mode
+ *	     TET: decompose cube and polygonize six tetrahedra
+ *	     NOTET: polygonize cube directly
+ *   returns error or NULL
+ */
+
+char *polygonize (function, size, bounds, x, y, z, triproc, mode)
+double (*function)(), size, x, y, z;
+int bounds, (*triproc)(), mode;
+{
+    PROCESS p;
+    int n, noabort;
+    CORNER *setcorner();
+    TEST in, out, find();
+
+    p.function = function;
+    p.triproc = triproc;
+    p.size = size;
+    p.bounds = bounds;
+    p.delta = size/(double)(RES*RES);
+
+    /* allocate hash tables and build cube polygon table: */
+    p.centers = (CENTERLIST **) mycalloc(HASHSIZE,sizeof(CENTERLIST *));
+    p.corners = (CORNERLIST **) mycalloc(HASHSIZE,sizeof(CORNERLIST *));
+    p.edges =	(EDGELIST   **) mycalloc(2*HASHSIZE,sizeof(EDGELIST *));
+    makecubetable();
+
+    /* find point on surface, beginning search at (x, y, z): */
+    srand(1);
+    in = find(1, &p, x, y, z);
+    out = find(0, &p, x, y, z);
+    if (!in.ok || !out.ok) {
+        free_cubetable();
+        free_process_data(&p);
+        clean_malloc();
+        return "can't find starting point";
+      }
+    converge(&in.p, &out.p, in.value, p.function, &p.start);
+
+    /* push initial cube on stack: */
+    p.cubes = (CUBES *) mycalloc(1, sizeof(CUBES)); /* list of 1 */
+    p.cubes->cube.i = p.cubes->cube.j = p.cubes->cube.k = 0;
+    p.cubes->next = NULL;
+
+    /* set corners of initial cube: */
+    for (n = 0; n < 8; n++)
+	p.cubes->cube.corners[n] = setcorner(&p, BIT(n,2), BIT(n,1), BIT(n,0));
+
+    p.vertices.count = p.vertices.max = 0; /* no vertices yet */
+    p.vertices.ptr = NULL;
+
+    setcenter(p.centers, 0, 0, 0);
+
+    while (p.cubes != NULL) { /* process active cubes till none left */
+        int i;
+	CUBE c;
+	CUBES *temp = p.cubes;
+	c = p.cubes->cube;
+
+	noabort = mode == TET?
+	       /* either decompose into tetrahedra and polygonize: */
+	       dotet(&c, LBN, LTN, RBN, LBF, &p) &&
+	       dotet(&c, RTN, LTN, LBF, RBN, &p) &&
+	       dotet(&c, RTN, LTN, LTF, LBF, &p) &&
+	       dotet(&c, RTN, RBN, LBF, RBF, &p) &&
+	       dotet(&c, RTN, LBF, LTF, RBF, &p) &&
+	       dotet(&c, RTN, LTF, RTF, RBF, &p)
+	       :
+	       /* or polygonize the cube directly: */
+	       docube(&c, &p);
+	if (! noabort) {
+            free_cubetable();
+            free_process_data(&p);
+            clean_malloc();
+            return "aborted";
+          }
+
+	/* pop current cube from stack */
+	p.cubes = p.cubes->next;
+
+	/* test six face directions, maybe add to stack: */
+	testface(c.i-1, c.j, c.k, &c, L, LBN, LBF, LTN, LTF, &p);
+	testface(c.i+1, c.j, c.k, &c, R, RBN, RBF, RTN, RTF, &p);
+	testface(c.i, c.j-1, c.k, &c, B, LBN, LBF, RBN, RBF, &p);
+	testface(c.i, c.j+1, c.k, &c, T, LTN, LTF, RTN, RTF, &p);
+	testface(c.i, c.j, c.k-1, &c, N, LBN, LTN, RBN, RTN, &p);
+	testface(c.i, c.j, c.k+1, &c, F, LBF, LTF, RBF, RTF, &p);
+
+        /* get rid of the current cube */
+        for (i=0; i<8; i++) {
+            myfree(temp->cube.corners[i]);
+            temp->cube.corners[i]=0;
+          }
+	myfree(temp);
+    }
+    free_cubetable();
+    free_process_data(&p);
+    clean_malloc();
+    return NULL;
+}
+
+static void
+free_process_data(p)
+    PROCESS *p;
+{
+  int i;
+  CUBES *cubes,*nextcubes;
+
+  if (p->vertices.ptr) myfree(p->vertices.ptr);
+
+  for (i=0; i<HASHSIZE; i++) {
+      CENTERLIST *l,*next;
+      for (l=p->centers[i]; l; l=next) {
+          next = l->next;
+          myfree(l);
+        }
+    }
+
+  for (i=0; i<HASHSIZE; i++) {
+      CORNERLIST *l,*next;
+      for (l=p->corners[i]; l; l=next) {
+          next = l->next;
+          myfree(l);
+        }
+    }
+
+  for (i=0; i<2*HASHSIZE; i++) {
+      EDGELIST *l,*next;
+      for (l=p->edges[i]; l; l=next) {
+          next = l->next;
+          myfree(l);
+        }
+    }
+
+  for (cubes=p->cubes; cubes; cubes=nextcubes) {
+      nextcubes = cubes->next;
+      for (i=0; i<8; i++) {
+          myfree(cubes->cube.corners[i]);
+        }
+      myfree(cubes);
+    }
+  
+  myfree(p->centers);
+  myfree(p->corners);
+  myfree(p->edges);
+}
+
+
+/* testface: given cube at lattice (i, j, k), and four corners of face,
+ * if surface crosses face, compute other four corners of adjacent cube
+ * and add new cube to cube stack */
+
+static void
+testface (i, j, k, old, face, c1, c2, c3, c4, p)
+CUBE *old;
+PROCESS *p;
+int i, j, k, face, c1, c2, c3, c4;
+{
+    CUBE new;
+    CUBES *oldcubes = p->cubes;
+    CORNER *setcorner();
+    int n, pos = old->corners[c1]->value > 0.0 ? 1 : 0;
+    /* static int facebit[6] = {2, 2, 1, 1, 0, 0}; */
+    /* int bit = facebit[face]; */
+
+    /* test if no surface crossing, cube out of bounds, or already visited: */
+    if ((old->corners[c2]->value > 0) == pos &&
+	(old->corners[c3]->value > 0) == pos &&
+	(old->corners[c4]->value > 0) == pos) return;
+    if (abs(i) > p->bounds || abs(j) > p->bounds || abs(k) > p->bounds) {
+        static int have_been_warned = 0;
+        if (!have_been_warned) {
+            fprintf(stderr,"WARNING: testface: cube out of bounds\n");
+            have_been_warned = 1;
+          }
+        /* abort(); */
+        return;
+      }
+    if (setcenter(p->centers, i, j, k)) return;
+
+    /* create new cube: */
+    new.i = i;
+    new.j = j;
+    new.k = k;
+    /* CLJ: changed this to make memory management possible. */
+/*     for (n = 0; n < 8; n++) new.corners[n] = NULL; */
+/*     new.corners[FLIP(c1, bit)] = old->corners[c1]; */
+/*     new.corners[FLIP(c2, bit)] = old->corners[c2]; */
+/*     new.corners[FLIP(c3, bit)] = old->corners[c3]; */
+/*     new.corners[FLIP(c4, bit)] = old->corners[c4]; */
+/*     for (n = 0; n < 8; n++) */
+/* 	if (new.corners[n] == NULL) */
+/* 	    new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0)); */
+    for (n = 0; n < 8; n++)
+        new.corners[n] = setcorner(p, i+BIT(n,2), j+BIT(n,1), k+BIT(n,0));
+
+    /*add cube to top of stack: */
+    p->cubes = (CUBES *) mycalloc(1, sizeof(CUBES));
+    p->cubes->cube = new;
+    p->cubes->next = oldcubes;
+}
+
+
+/* setcorner: return corner with the given lattice location
+   set (and cache) its function value */
+
+static CORNER *setcorner (p, i, j, k)
+int i, j, k;
+PROCESS *p;
+{
+    /* for speed, do corner value caching here */
+    CORNER *c = (CORNER *) mycalloc(1, sizeof(CORNER));
+    int index = HASH(i, j, k);
+    CORNERLIST *l = p->corners[index];
+    c->i = i; c->x = p->start.x+((double)i-.5)*p->size;
+    c->j = j; c->y = p->start.y+((double)j-.5)*p->size;
+    c->k = k; c->z = p->start.z+((double)k-.5)*p->size;
+    for (; l != NULL; l = l->next)
+	if (l->i == i && l->j == j && l->k == k) {
+	    c->value = l->value;
+	    return c;
+	    }
+    l = (CORNERLIST *) mycalloc(1, sizeof(CORNERLIST));
+    l->i = i; l->j = j; l->k = k;
+    l->value = c->value = p->function(c->x, c->y, c->z);
+    if (c->value > 100.0 || c->value < -100.0) {
+        fprintf(stderr,"suspicious\n");
+        abort();
+      }
+    l->next = p->corners[index];
+    p->corners[index] = l;
+    return c;
+}
+
+
+/* find: search for point with value of given sign (0: neg, 1: pos) */
+
+static TEST find (sign, p, x, y, z)
+int sign;
+PROCESS *p;
+double x, y, z;
+{
+    int i;
+    TEST test;
+    double range = p->size;
+    test.ok = 1;
+    for (i = 0; i < 10000; i++) {
+	test.p.x = x+range*(RAND()-0.5);
+	test.p.y = y+range*(RAND()-0.5);
+	test.p.z = z+range*(RAND()-0.5);
+	test.value = p->function(test.p.x, test.p.y, test.p.z);
+	if (sign == (test.value > 0.0)) return test;
+	range = range*1.0005; /* slowly expand search outwards */
+    }
+    test.ok = 0;
+    return test;
+}
+
+
+/**** Tetrahedral Polygonization ****/
+
+
+/* dotet: triangulate the tetrahedron
+ * b, c, d should appear clockwise when viewed from a
+ * return 0 if client aborts, 1 otherwise */
+
+static int dotet (cube, c1, c2, c3, c4, p)
+CUBE *cube;
+int c1, c2, c3, c4;
+PROCESS *p;
+{
+    CORNER *a = cube->corners[c1];
+    CORNER *b = cube->corners[c2];
+    CORNER *c = cube->corners[c3];
+    CORNER *d = cube->corners[c4];
+    int index = 0, apos, bpos, cpos, dpos, e1=0, e2=0, e3=0, e4=0, e5=0, e6=0;
+    if ((apos = (a->value > 0.0))) index += 8;
+    if ((bpos = (b->value > 0.0))) index += 4;
+    if ((cpos = (c->value > 0.0))) index += 2;
+    if ((dpos = (d->value > 0.0))) index += 1;
+    /* index is now 4-bit number representing one of the 16 possible cases */
+    if (apos != bpos) e1 = vertid(a, b, p);
+    if (apos != cpos) e2 = vertid(a, c, p);
+    if (apos != dpos) e3 = vertid(a, d, p);
+    if (bpos != cpos) e4 = vertid(b, c, p);
+    if (bpos != dpos) e5 = vertid(b, d, p);
+    if (cpos != dpos) e6 = vertid(c, d, p);
+    /* 14 productive tetrahedral cases (0000 and 1111 do not yield polygons */
+    switch (index) {
+	case 1:	 return p->triproc(e5, e6, e3, p->vertices);
+	case 2:	 return p->triproc(e2, e6, e4, p->vertices);
+	case 3:	 return p->triproc(e3, e5, e4, p->vertices) &&
+			p->triproc(e3, e4, e2, p->vertices);
+	case 4:	 return p->triproc(e1, e4, e5, p->vertices);
+	case 5:	 return p->triproc(e3, e1, e4, p->vertices) &&
+			p->triproc(e3, e4, e6, p->vertices);
+	case 6:	 return p->triproc(e1, e2, e6, p->vertices) &&
+			p->triproc(e1, e6, e5, p->vertices);
+	case 7:	 return p->triproc(e1, e2, e3, p->vertices);
+	case 8:	 return p->triproc(e1, e3, e2, p->vertices);
+	case 9:	 return p->triproc(e1, e5, e6, p->vertices) &&
+			p->triproc(e1, e6, e2, p->vertices);
+	case 10: return p->triproc(e1, e3, e6, p->vertices) &&
+			p->triproc(e1, e6, e4, p->vertices);
+	case 11: return p->triproc(e1, e5, e4, p->vertices);
+	case 12: return p->triproc(e3, e2, e4, p->vertices) &&
+			p->triproc(e3, e4, e5, p->vertices);
+	case 13: return p->triproc(e6, e2, e4, p->vertices);
+	case 14: return p->triproc(e5, e3, e6, p->vertices);
+    }
+    return 1;
+}
+
+
+/**** Cubical Polygonization (optional) ****/
+
+
+#define LB	0  /* left bottom edge	*/
+#define LT	1  /* left top edge	*/
+#define LN	2  /* left near edge	*/
+#define LF	3  /* left far edge	*/
+#define RB	4  /* right bottom edge */
+#define RT	5  /* right top edge	*/
+#define RN	6  /* right near edge	*/
+#define RF	7  /* right far edge	*/
+#define BN	8  /* bottom near edge	*/
+#define BF	9  /* bottom far edge	*/
+#define TN	10 /* top near edge	*/
+#define TF	11 /* top far edge	*/
+
+static INTLISTS *cubetable[256];
+
+/*			edge: LB, LT, LN, LF, RB, RT, RN, RF, BN, BF, TN, TF */
+static int corner1[12]	   = {LBN,LTN,LBN,LBF,RBN,RTN,RBN,RBF,LBN,LBF,LTN,LTF};
+static int corner2[12]	   = {LBF,LTF,LTN,LTF,RBF,RTF,RTN,RTF,RBN,RBF,RTN,RTF};
+static int leftface[12]	   = {B,  L,  L,  F,  R,  T,  N,  R,  N,  B,  T,  F};
+			     /* face on left when going corner1 to corner2 */
+static int rightface[12]   = {L,  T,  N,  L,  B,  R,  R,  F,  B,  F,  N,  T};
+			     /* face on right when going corner1 to corner2 */
+
+
+/* docube: triangulate the cube directly, without decomposition */
+
+static int docube (cube, p)
+CUBE *cube;
+PROCESS *p;
+{
+    INTLISTS *polys;
+    int i, index = 0;
+    for (i = 0; i < 8; i++) if (cube->corners[i]->value > 0.0) index += (1<<i);
+    for (polys = cubetable[index]; polys; polys = polys->next) {
+	INTLIST *edges;
+	int a = -1, b = -1, count = 0;
+	for (edges = polys->list; edges; edges = edges->next) {
+	    CORNER *c1 = cube->corners[corner1[edges->i]];
+	    CORNER *c2 = cube->corners[corner2[edges->i]];
+	    int c = vertid(c1, c2, p);
+	    if (++count > 2 && ! p->triproc(a, b, c, p->vertices)) return 0;
+	    if (count < 3) a = b;
+	    b = c;
+	}
+    }
+    return 1;
+}
+
+
+/* nextcwedge: return next clockwise edge from given edge around given face */
+
+static int nextcwedge (edge, face)
+int edge, face;
+{
+    switch (edge) {
+	case LB: return (face == L)? LF : BN;
+	case LT: return (face == L)? LN : TF;
+	case LN: return (face == L)? LB : TN;
+	case LF: return (face == L)? LT : BF;
+	case RB: return (face == R)? RN : BF;
+	case RT: return (face == R)? RF : TN;
+	case RN: return (face == R)? RT : BN;
+	case RF: return (face == R)? RB : TF;
+	case BN: return (face == B)? RB : LN;
+	case BF: return (face == B)? LB : RF;
+	case TN: return (face == T)? LT : RN;
+	case TF: return (face == T)? RT : LF;
+    }
+
+    return -1;
+}
+
+
+/* otherface: return face adjoining edge that is not the given face */
+
+static int otherface (edge, face)
+int edge, face;
+{
+    int other = leftface[edge];
+    return face == other? rightface[edge] : other;
+}
+
+
+/* makecubetable: create the 256 entry table for cubical polygonization */
+
+static void makecubetable ()
+{
+    int i, e, c, done[12], pos[8];
+    memset(cubetable, 0, sizeof(cubetable));
+    for (i = 0; i < 256; i++) {
+	for (e = 0; e < 12; e++) done[e] = 0;
+	for (c = 0; c < 8; c++) pos[c] = BIT(i, c);
+	for (e = 0; e < 12; e++)
+	    if (!done[e] && (pos[corner1[e]] != pos[corner2[e]])) {
+		INTLIST *ints = 0;
+		INTLISTS *lists = (INTLISTS *) mycalloc(1, sizeof(INTLISTS));
+		int start = e, edge = e;
+		/* get face that is to right of edge from pos to neg corner: */
+		int face = pos[corner1[e]]? rightface[e] : leftface[e];
+		while (1) {
+		    edge = nextcwedge(edge, face);
+		    done[edge] = 1;
+		    if (pos[corner1[edge]] != pos[corner2[edge]]) {
+			INTLIST *tmp = ints;
+			ints = (INTLIST *) mycalloc(1, sizeof(INTLIST));
+			ints->i = edge;
+			ints->next = tmp; /* add edge to head of list */
+			if (edge == start) break;
+			face = otherface(edge, face);
+		    }
+		}
+		lists->list = ints; /* add ints to head of table entry */
+		lists->next = cubetable[i];
+		cubetable[i] = lists;
+	    }
+    }
+}
+
+static void
+free_cubetable()
+{
+  int i;
+  for (i=0; i<256; i++) {
+      INTLISTS *l,*nextl;
+      for (l=cubetable[i]; l; l=nextl) {
+          INTLIST *m, *nextm;
+          for (m=l->list; m; m=nextm) {
+              nextm = m->next;
+              myfree(m);
+            }
+          nextl = l->next;
+          myfree(l);
+        }
+    }
+}
+
+/**** Storage ****/
+
+#undef CHECK_MALLOC
+
+#ifdef CHECK_MALLOC
+static char allocwarn[10000];
+static char delwarn[10000];
+#endif
+
+/* mycalloc: return successful calloc or exit program */
+
+typedef struct mallocdata {
+    int lineno;
+    char* ptr;
+    size_t size;
+    struct mallocdata* next;
+} MALLOCDATA;
+
+#ifdef CHECK_MALLOC
+static MALLOCDATA *malloc_list;
+static void add_mallocdata(char* ptr, int lineno, size_t size)
+{
+  MALLOCDATA * old = malloc_list;
+  malloc_list = (MALLOCDATA*) malloc(sizeof(MALLOCDATA));
+  malloc_list->next = old;
+  malloc_list->ptr = ptr;
+  malloc_list->size = size;
+  malloc_list->lineno = lineno;
+}
+
+static size_t del_mallocdata(char* ptr,int lineno)
+{
+  MALLOCDATA *i, *ilast = 0;
+  int size;
+  for (i=malloc_list; i; ilast=i,i=i->next) {
+      if (i->ptr == ptr) {
+          if (ilast) {
+              MALLOCDATA * tmp = i->next;
+              ilast->next = i->next;
+            }
+          else {
+              malloc_list = i->next;
+            }
+          size = i->size;
+          free(i);
+          return size;
+        }
+    }
+  if (!delwarn[lineno]) {
+      fprintf(stderr,"tried to delete unknown data at line %d\n",lineno);
+      delwarn[lineno] = 1;
+    }
+  return 0;
+}
+#endif
+
+static void clean_malloc()
+{
+#ifdef CHECK_MALLOC
+  MALLOCDATA*i;
+  int count=0;
+  for (i=malloc_list; i; i=i->next) {
+      if (!allocwarn[i->lineno]) {
+          fprintf(stderr,"have memory allocated from line %d\n",i->lineno);
+          allocwarn[i->lineno] = 1;
+        }
+      count++;
+    }
+  fprintf(stderr,"%d allocated pieces of memory remain\n",count);
+#endif
+}
+
+static char *_mycalloc (nitems, nbytes, line)
+int nitems, nbytes, line;
+{
+   char *ptr = calloc(nitems, nbytes);
+#ifdef CHECK_MALLOC
+   add_mallocdata(ptr,line,nitems*nbytes);
+#endif
+   if (ptr != NULL) return ptr;
+   fprintf(stderr, "can't calloc %d bytes\n", nitems*nbytes);
+   abort();
+   return 0;
+}
+
+static void _myfree(ptr, lineno)
+    void* ptr;
+    int lineno;
+{
+#ifdef CHECK_MALLOC
+  size_t size = del_mallocdata(ptr,lineno);
+  char*tmp = ptr;
+  for (int i=0; i<size; i++) {
+      *tmp++ = 0x00;
+    }
+#endif
+
+  free(ptr);
+}
+
+
+/* setcenter: set (i,j,k) entry of table[]
+ * return 1 if already set; otherwise, set and return 0 */
+
+static int setcenter(table, i, j, k)
+CENTERLIST *table[];
+int i, j, k;
+{
+    int index = HASH(i, j, k);
+    CENTERLIST *new, *l, *q = table[index];
+    for (l = q; l != NULL; l = l->next)
+	if (l->i == i && l->j == j && l->k == k) return 1;
+    new = (CENTERLIST *) mycalloc(1, sizeof(CENTERLIST));
+    new->i = i; new->j = j; new->k = k; new->next = q;
+    table[index] = new;
+    return 0;
+}
+
+
+/* setedge: set vertex id for edge */
+
+static void setedge (table, i1, j1, k1, i2, j2, k2, vid)
+EDGELIST *table[];
+int i1, j1, k1, i2, j2, k2, vid;
+{
+    unsigned int index;
+    EDGELIST *new;
+    if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
+	int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
+    }
+    index = HASH(i1, j1, k1) + HASH(i2, j2, k2);
+    new = (EDGELIST *) mycalloc(1, sizeof(EDGELIST));
+    new->i1 = i1; new->j1 = j1; new->k1 = k1;
+    new->i2 = i2; new->j2 = j2; new->k2 = k2;
+    new->vid = vid;
+    new->next = table[index];
+    table[index] = new;
+}
+
+
+/* getedge: return vertex id for edge; return -1 if not set */
+
+static int getedge (table, i1, j1, k1, i2, j2, k2)
+EDGELIST *table[];
+int i1, j1, k1, i2, j2, k2;
+{
+    EDGELIST *q;
+    if (i1>i2 || (i1==i2 && (j1>j2 || (j1==j2 && k1>k2)))) {
+	int t=i1; i1=i2; i2=t; t=j1; j1=j2; j2=t; t=k1; k1=k2; k2=t;
+    };
+    q = table[HASH(i1, j1, k1)+HASH(i2, j2, k2)];
+    for (; q != NULL; q = q->next)
+	if (q->i1 == i1 && q->j1 == j1 && q->k1 == k1 &&
+	    q->i2 == i2 && q->j2 == j2 && q->k2 == k2)
+	    return q->vid;
+    return -1;
+}
+
+
+/**** Vertices ****/
+
+
+/* vertid: return index for vertex on edge:
+ * c1->value and c2->value are presumed of different sign
+ * return saved index if any; else compute vertex and save */
+
+static int vertid (c1, c2, p)
+CORNER *c1, *c2;
+PROCESS *p;
+{
+    VERTEX v;
+    POINT a, b;
+    int vid = getedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k);
+    if (vid != -1) return vid;			     /* previously computed */
+    a.x = c1->x; a.y = c1->y; a.z = c1->z;
+    b.x = c2->x; b.y = c2->y; b.z = c2->z;
+    converge(&a, &b, c1->value, p->function, &v.position); /* position */
+    vnormal(&v.position, p, &v.normal);			   /* normal */
+    addtovertices(&p->vertices, v);			   /* save vertex */
+    vid = p->vertices.count-1;
+    setedge(p->edges, c1->i, c1->j, c1->k, c2->i, c2->j, c2->k, vid);
+    return vid;
+}
+
+
+/* addtovertices: add v to sequence of vertices */
+
+static void addtovertices (vertices, v)
+VERTICES *vertices;
+VERTEX v;
+{
+    if (vertices->count == vertices->max) {
+	int i;
+	VERTEX *new;
+	vertices->max = vertices->count == 0 ? 10 : 2*vertices->count;
+	new = (VERTEX *) mycalloc(vertices->max, sizeof(VERTEX));
+	for (i = 0; i < vertices->count; i++) new[i] = vertices->ptr[i];
+	if (vertices->ptr != NULL) myfree(vertices->ptr);
+	vertices->ptr = new;
+    }
+    vertices->ptr[vertices->count++] = v;
+}
+
+
+/* vnormal: compute unit length surface normal at point */
+
+static void vnormal (point, p, v)
+POINT *point, *v;
+PROCESS *p;
+{
+    double f = p->function(point->x, point->y, point->z);
+    v->x = p->function(point->x+p->delta, point->y, point->z)-f;
+    v->y = p->function(point->x, point->y+p->delta, point->z)-f;
+    v->z = p->function(point->x, point->y, point->z+p->delta)-f;
+    f = sqrt(v->x*v->x + v->y*v->y + v->z*v->z);
+    if (f != 0.0) {v->x /= f; v->y /= f; v->z /= f;}
+}
+
+
+/* converge: from two points of differing sign, converge to zero crossing */
+
+static void converge (p1, p2, v, function, p)
+double v;
+double (*function)();
+POINT *p1, *p2, *p;
+{
+    int i = 0;
+    POINT pos, neg;
+    if (v < 0) {
+	pos.x = p2->x; pos.y = p2->y; pos.z = p2->z;
+	neg.x = p1->x; neg.y = p1->y; neg.z = p1->z;
+    }
+    else {
+	pos.x = p1->x; pos.y = p1->y; pos.z = p1->z;
+	neg.x = p2->x; neg.y = p2->y; neg.z = p2->z;
+    }
+    while (1) {
+	p->x = 0.5*(pos.x + neg.x);
+	p->y = 0.5*(pos.y + neg.y);
+	p->z = 0.5*(pos.z + neg.z);
+	if (i++ == RES) return;
+	if ((function(p->x, p->y, p->z)) > 0.0)
+	     {pos.x = p->x; pos.y = p->y; pos.z = p->z;}
+	else {neg.x = p->x; neg.y = p->y; neg.z = p->z;}
+    }
+}
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/implicit.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,28 @@
+
+#ifndef _math_isosurf_implicit_h
+#define _math_isosurf_implicit_h
+
+typedef struct point {		   /* a three-dimensional point */
+    double x, y, z;		   /* its coordinates */
+} POINT;
+
+typedef struct vertex {		   /* surface vertex */
+    POINT position, normal;	   /* position and surface normal */
+} VERTEX;
+
+typedef struct vertices {	   /* list of vertices in polygonization */
+    int count, max;		   /* # vertices, max # allowed */
+    VERTEX *ptr;		   /* dynamically allocated */
+} VERTICES;
+
+#define TET	0  /* use tetrahedral decomposition */
+#define NOTET	1  /* no tetrahedral decomposition  */
+
+extern "C" {
+    char * polygonize(double(*function)(double,double,double),
+                      double size, int bounds,
+                      double x, double y, double z,
+                      int(*triproc)(int,int,int,VERTICES), int mode);
+}
+
+#endif
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,255 @@
+//
+// isosurf.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <iostream>
+
+#include <math/isosurf/isosurf.h>
+#include <math/isosurf/implicit.h>
+#include <math/scmat/vector3.h>
+
+using namespace std;
+using namespace sc;
+
+// This static datum is used to interface to the implicit.c routine
+// provided in Graphics Gems IV.
+ImplicitSurfacePolygonizer* ImplicitSurfacePolygonizer::current = 0;
+
+// These functions are used to interface to the implicit.c routine provided
+// in Graphics Gems IV.
+extern "C" int
+ImplicitSurfacePolygonizer_add_triangle_to_current(int i1, int i2, int i3,
+                                                   VERTICES v)
+{
+  return ImplicitSurfacePolygonizer::add_triangle_to_current(i1, i2, i3, v);
+}
+
+extern "C" double
+ImplicitSurfacePolygonizer_value_of_current(double x,double y,double z)
+{
+  return ImplicitSurfacePolygonizer::value_of_current(x,y,z);
+}
+
+////////////////////////////////////////////////////////////////////////////
+// IsosurfaceGen members
+
+IsosurfaceGen::IsosurfaceGen():
+  _resolution(0.25)
+{
+}
+
+IsosurfaceGen::~IsosurfaceGen()
+{
+}
+
+void
+IsosurfaceGen::set_resolution(double r)
+{
+  _resolution = r;
+}
+
+////////////////////////////////////////////////////////////////////////////
+// ImplicitSurfacePolygonizer members
+
+ImplicitSurfacePolygonizer::ImplicitSurfacePolygonizer(const Ref<Volume>&vol):
+  _volume(vol)
+  ,_tmp_vertices(0)
+{
+}
+
+ImplicitSurfacePolygonizer::~ImplicitSurfacePolygonizer()
+{
+}
+
+static SCVector3 current_x;
+void
+ImplicitSurfacePolygonizer::isosurface(double value,
+                                       TriangulatedSurface& surf)
+{
+  surf.clear();
+
+  // Find the bounding box.
+  SCVector3 p0;
+  SCVector3 p1;
+  _volume->boundingbox(value - 0.001, value + 0.001, p0, p1);
+  SCVector3 diag = p1 - p0;
+  SCVector3 midpoint = 0.5*diag + p0;
+  double biggest_width = diag.maxabs();
+  int bounds = (int)(0.5*biggest_width/_resolution) + 2;
+
+  // polygonize will find a starting point and do bounds checking
+  // from that point.  To make sure the bounding box is big enough
+  // its size must be doubled.  Since polygonization is implicit
+  // there is no performance penalty.
+  bounds *= 2;
+
+  // Initialize the static pointer to this, so the C polygonizer can find us.
+  current = this;
+  _surf = &surf;
+  _value = value;
+  // Find the polygons.
+  char *msg = polygonize(ImplicitSurfacePolygonizer_value_of_current, _resolution, bounds,
+                         midpoint[0], midpoint[1], midpoint[2],
+                         ImplicitSurfacePolygonizer_add_triangle_to_current,
+                         NOTET);
+  current = 0;
+  _surf = 0;
+  if (msg) {
+      ExEnv::errn() << "ImplicitSurfacePolygonizer::isosurface: failed: "
+           << msg << endl;
+      abort();
+    }
+
+  // Clean up temporaries.
+  _tmp_vertices.clear();
+
+  ExEnv::out0() << "about to complete the surface" << endl;
+
+  // finish the surface
+  surf.complete_surface();
+
+  ExEnv::out0() << "completed the surface" << endl;
+  ExEnv::out0() << "flat area = " << surf.flat_area() << endl;
+  ExEnv::out0() << "  ntri = " << setw(10) << surf.ntriangle()
+       << " bytes = "
+       << setw(10) << surf.ntriangle() * sizeof(Triangle)
+       << endl;
+  ExEnv::out0() << "  nedg = " << setw(10) << surf.nedge()
+       << " bytes = "
+       << setw(10) << surf.nedge() * sizeof(Edge)
+       << endl;
+  ExEnv::out0() << "  nver = " << setw(10) << surf.nvertex()
+       << " bytes = "
+       << setw(10) << surf.nvertex() * sizeof(Vertex)
+       << endl;
+
+  // compute normals if they weren't computed from the gradients
+  if (!_volume->gradient_implemented()) {
+      int i,j;
+      // compute the normals as the average of the normals of
+      // all the connected triangles
+      for (i=0; i<surf.ntriangle(); i++) {
+          Ref<Triangle> t = surf.triangle(i);
+          SCVector3 tmp;
+          SCVector3 BA = t->vertex(1)->point() - t->vertex(0)->point();
+          SCVector3 CA = t->vertex(2)->point() - t->vertex(0)->point();
+          SCVector3 N = BA.cross(CA);
+          double n = N.norm();
+          if (n < 1.0e-8) {
+              tmp = 0.0;
+            }
+          else {
+              n = 1.0/n;
+              for (int j=0; j<3; j++) {
+                  tmp[j] = - N[j]*n;
+                }
+            }
+          for (j=0; j<3; j++) {
+              int iv = surf.vertex_index(t->vertex(j));
+              if (iv>=0) {
+                  Ref<Vertex> v = surf.vertex(iv);
+                  if (v->has_normal()) {
+                      tmp += v->normal();
+                    }
+                  tmp.normalize();
+                  v->set_normal(tmp);
+                }
+            }
+        }
+      // normalize all the normals
+      for (i=0; i<surf.nvertex(); i++) {
+          Ref<Vertex> v = surf.vertex(i);
+          if (v->has_normal()) {
+              SCVector3 n = v->normal();
+              n.normalize();
+              v->set_normal(n);
+            }
+          else {
+              ExEnv::outn() << "ERROR: isosurf has a vertex without a triangle" << endl;
+              abort();
+            }
+        }
+    }
+}
+
+double
+ImplicitSurfacePolygonizer::value_of_current(double x,double y,double z)
+{
+  current_x[0] = x; current_x[1] = y; current_x[2] = z;
+  current->_volume->set_x(current_x);
+  return current->_volume->value() - current->_value;
+}
+
+int
+ImplicitSurfacePolygonizer::add_triangle_to_current(int i1, int i2, int i3,
+                                                    VERTICES v)
+{
+  int oldlength = current->_tmp_vertices.size();
+  current->_tmp_vertices.resize(v.count);
+  for (int i=oldlength; i<v.count; i++) {
+      SCVector3 newpoint;
+      newpoint[0] = v.ptr[i].position.x;
+      newpoint[1] = v.ptr[i].position.y;
+      newpoint[2] = v.ptr[i].position.z;
+      current->_volume->set_x(newpoint);
+      SCVector3 normal;
+      if (current->_volume->gradient_implemented()) {
+          current->_volume->get_gradient(normal);
+          normal.normalize();
+          current->_tmp_vertices[i] = new Vertex(newpoint, normal);
+        }
+      else {
+          current->_tmp_vertices[i] = new Vertex(newpoint);
+        }
+    }
+
+  Ref<Vertex> v1 = current->_tmp_vertices[i1];
+  Ref<Vertex> v2 = current->_tmp_vertices[i2];
+  Ref<Vertex> v3 = current->_tmp_vertices[i3];
+  
+  static int tricnt = 0;
+  if (++tricnt%100 == 0) {
+      ExEnv::out0() << "adding triangle " << tricnt << endl;
+      ExEnv::out0() << "  ntri = " << setw(10) << current->_surf->ntriangle()
+           << " bytes = "
+           << setw(10) << current->_surf->ntriangle() * sizeof(Triangle)
+           << endl;
+    }
+  current->_surf->add_triangle(v1,v2,v3);
+
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/isosurf.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,85 @@
+//
+// isosurf.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_isosurf_h
+#define _math_isosurf_isosurf_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <vector>
+
+#include <math/isosurf/surf.h>
+
+namespace sc {
+
+class IsosurfaceGen {
+  protected:
+    double _resolution;
+  public:
+    IsosurfaceGen();
+    virtual ~IsosurfaceGen();
+    virtual void isosurface(double value,
+                            TriangulatedSurface& surf) = 0;
+    virtual void set_resolution(double);
+};
+
+class ImplicitSurfacePolygonizer: public IsosurfaceGen {
+    // These static data and members are used to interface to the
+    // implicit.c routine provided in Graphics Gems IV.
+    static ImplicitSurfacePolygonizer* current;
+    // The following should not really be used publically.
+    // they are public to permit access through internal
+    // C-language functions.
+  public:
+    /// For internal use only.
+    static int add_triangle_to_current(int,int,int,VERTICES);
+    /// For internal use only.
+    static double value_of_current(double x, double y, double z);
+  protected:
+    Ref<Volume> _volume;
+
+    std::vector<Ref<Vertex> >  _tmp_vertices;
+
+    TriangulatedSurface* _surf;
+    double _value;
+  public:
+    ImplicitSurfacePolygonizer(const Ref<Volume>&);
+    virtual ~ImplicitSurfacePolygonizer();
+    virtual void isosurface(double value,
+                            TriangulatedSurface& surf);
+};  
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,106 @@
+//
+// isotest.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+extern "C" {
+#include <math.h>
+#ifdef SGI
+#include <sigfpe.h>
+#endif // SGI
+}
+
+#include <util/misc/ieee.h>
+#include <util/misc/formio.h>
+
+#include <util/keyval/keyval.h>
+#include <math/isosurf/volume.h>
+#include <math/isosurf/shape.h>
+#include <math/isosurf/surf.h>
+#include <math/isosurf/isosurf.h>
+
+using namespace std;
+using namespace sc;
+
+// Force linkages:
+#ifndef __PIC__
+static ForceLink<SphereShape> fl0;
+static ForceLink<TriangulatedImplicitSurface> fl1;
+#endif
+
+main(int argc,char** argv)
+{
+  ieee_trap_errors();
+
+  // The first argument is the input filename.  If it doesn't
+  // exist or is a '.', then the default ioput file is used.
+  const char *defaultinput = SRCDIR "/isotest.in";
+  const char *input;
+  if (argc == 1 || !strcmp(argv[1], ".")) input = defaultinput;
+  else input = argv[1];
+
+  const char *keyword = (argc > 2)? argv[2] : "surf";
+
+  // open keyval input
+  Ref<KeyVal> rpkv(new ParsedKeyVal(input));
+
+  Ref<TriangulatedSurface> surf;
+  surf << rpkv->describedclassvalue(keyword);
+
+  cout << scprintf("surf->flat_area() = %f\n", surf->flat_area());
+  cout << scprintf("surf->flat_volume() = %f\n", surf->flat_volume());
+
+  cout << scprintf("surf->area() = %f\n", surf->area());
+  cout << scprintf("surf->volume() = %f\n", surf->volume());
+
+  for (int order=1; order<4; order++) { 
+      cout << scprintf("order = %d", order) << endl;
+      int nir = 4;
+      for (int ir=0; ir<nir; ir++) {
+          double r = double(ir)/(nir-1);
+          double s = 1.0 - r;
+          cout << scprintf("  r = %6.4f, s = %6.4f", r, s) << endl;
+          TriInterpCoefKey key(order, r, s);
+          Ref<TriInterpCoef> coef = new TriInterpCoef(key);
+          for (int i=0; i<=order; i++) {
+              for (int j=0; j+i<=order; j++) {
+                  int k = order - i - j;
+                  cout << scprintf("    coef(%d,%d,%d) = %12.8f",
+                                   i,j,k,coef->coef(i,j,k))
+                       << endl;
+                }
+            }
+        }
+    }
+
+  return 0;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.in
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.in	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/isotest.in	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,22 @@
+% Emacs should use -*- KeyVal -*- mode.
+
+surf<TriangulatedImplicitSurface>: (
+    % TriangulatedSurface
+    verbose = yes
+    debug = no
+    integrator<GaussTriangleIntegrator>: (
+         n = 7
+      )
+    % TriangulatedImplicitSurface
+    volume<SphereShape>: (
+        origin = [ 0.0 0.0 0.0 ]
+        radius = 1.0
+      )
+    resolution = 0.25
+    remove_short_edges = yes
+    remove_slender_triangles = yes
+    short_edge_factor = 0.4
+    slender_triangle_factor = 0.2
+    order = 3
+    initialize = yes
+  )
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/shape.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/shape.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/shape.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,1180 @@
+//
+// shape.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <stdio.h>
+#include <util/misc/math.h>
+
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/isosurf/shape.h>
+
+using namespace std;
+using namespace sc;
+
+static const double shape_infinity = 1.0e23;
+
+// given a vector X find which of the points in the vector of
+// vectors, A, is closest to it and return the distance
+static double
+closest_distance(SCVector3& X,SCVector3*A,int n,SCVector3*grad)
+{
+  SCVector3 T = X-A[0];
+  double min = T.dot(T);
+  int imin = 0;
+  for (int i=1; i<n; i++) {
+      T = X-A[i];
+      double tmp = T.dot(T);
+      if (tmp < min) {min = tmp; imin = i;}
+    }
+  if (grad) {
+      T = X - A[imin];
+      T.normalize();
+      *grad = T;
+    }
+  return sqrt(min);
+}
+
+//////////////////////////////////////////////////////////////////////
+// Shape
+
+static ClassDesc Shape_cd(
+  typeid(Shape),"Shape",1,"public Volume",
+  0, 0, 0);
+
+Shape::Shape():
+  Volume()
+{
+}
+
+Shape::Shape(const Ref<KeyVal>& keyval):
+  Volume(keyval)
+{
+}
+
+Shape::~Shape()
+{
+}
+
+void
+Shape::compute()
+{
+  SCVector3 r;
+  get_x(r);
+  if (gradient_needed()) {
+      if (!gradient_implemented()) {
+          ExEnv::errn() << "Shape::compute: gradient not implemented" << endl;
+          abort();
+        }
+      SCVector3 v;
+      set_value(distance_to_surface(r,&v));
+      set_actual_value_accuracy(desired_value_accuracy());
+      set_gradient(v);
+      set_actual_gradient_accuracy(desired_gradient_accuracy());
+    }
+  else if (value_needed()) {
+      set_value(distance_to_surface(r));
+      set_actual_value_accuracy(desired_value_accuracy());
+    }
+  if (hessian_needed()) {
+      ExEnv::errn() << "Shape::compute(): can't do hessian yet" << endl;
+      abort();
+    }
+}
+
+int
+Shape::is_outside(const SCVector3&r) const
+{
+  if (distance_to_surface(r)>0.0) return 1;
+  return 0;
+}
+
+// Shape overrides volume's interpolate so it always gets points on
+// the outside of the shape are always returned.
+
+// interpolate using the bisection algorithm
+void
+Shape::interpolate(const SCVector3& A,
+                   const SCVector3& B,
+                   double val,
+                   SCVector3& result)
+{
+  if (val < 0.0) {
+      failure("Shape::interpolate(): val is < 0.0");
+    }
+
+  set_x(A);
+  double value0 = value() - val;
+
+  set_x(B);
+  double value1 = value() - val;
+
+  if (value0*value1 > 0.0) {
+      failure("Shape::interpolate(): values at endpoints don't bracket val");
+    }
+  else if (value0 == 0.0) {
+      result = A;
+      return;
+    }
+  else if (value1 == 0.0) {
+      result = B;
+      return;
+    }
+
+  SCVector3 BA = B - A;
+
+  double length = BA.norm();
+  int niter = (int) (log(length/interpolation_accuracy())/M_LN2);
+
+  double f0 = 0.0;
+  double f1 = 1.0;
+  double fnext = 0.5;
+
+  SCVector3 X = A + fnext*BA;
+  set_x(X);
+  double valuenext = value() - val;
+
+  do {
+      for (int i=0; i<niter; i++) {
+          if (valuenext*value0 <= 0.0) {
+              value1 = valuenext;
+              f1 = fnext;
+              fnext = (f0 + fnext)*0.5;
+            }
+          else {
+              value0 = valuenext;
+              f0 = fnext;
+              fnext = (fnext + f1)*0.5;
+            }
+          X = A + fnext*BA;
+          set_x(X);
+          valuenext = value() - val;
+        }
+      niter = 1;
+    } while (valuenext < 0.0);
+
+  result = X;
+}
+
+int
+Shape::value_implemented() const
+{
+  return 1;
+}
+
+//////////////////////////////////////////////////////////////////////
+// SphereShape
+
+static ClassDesc SphereShape_cd(
+  typeid(SphereShape),"SphereShape",1,"public Shape",
+  0, create<SphereShape>, 0);
+
+SphereShape::SphereShape(const SCVector3&o,double r):
+  _origin(o),
+  _radius(r)
+{
+}
+
+SphereShape::SphereShape(const SphereShape&s):
+  _origin(s._origin),
+  _radius(s._radius)
+{
+}
+
+SphereShape::SphereShape(const Ref<KeyVal>& keyval):
+  _origin(new PrefixKeyVal(keyval,"origin")),
+  _radius(keyval->doublevalue("radius"))
+{
+}
+
+SphereShape::~SphereShape()
+{
+}
+
+double
+SphereShape::distance_to_surface(const SCVector3&p,SCVector3*grad) const
+{
+  int i;
+  double r2 = 0.0;
+  for (i=0; i<3; i++) {
+      double tmp = p[i] - _origin[i];
+      r2 += tmp*tmp;
+    }
+  double r = sqrt(r2);
+  double d = r - _radius;
+  if (grad) {
+      SCVector3 pv(p);
+      SCVector3 o(_origin);
+      SCVector3 unit = pv - o;
+      unit.normalize();
+      for (i=0; i<3; i++) grad->elem(i) = unit[i];
+    }
+  return d;
+}
+
+void SphereShape::print(ostream&o) const
+{
+  o << indent
+    << scprintf("SphereShape: r = %8.4f o = (%8.4f %8.4f %8.4f)",
+                radius(),origin()[0],origin()[1],origin()[2])
+    << endl;
+}
+
+void
+SphereShape::boundingbox(double valuemin, double valuemax,
+                         SCVector3& p1,
+                         SCVector3& p2)
+{
+  if (valuemax < 0.0) valuemax = 0.0;
+
+  int i;
+  for (i=0; i<3; i++) {
+      p1[i] = _origin[i] - _radius - valuemax;
+      p2[i] = _origin[i] + _radius + valuemax;
+    }
+}
+
+int
+SphereShape::gradient_implemented() const
+{
+  return 1;
+}
+
+////////////////////////////////////////////////////////////////////////
+// UncappedTorusHoleShape
+
+static ClassDesc UncappedTorusHoleShape_cd(
+  typeid(UncappedTorusHoleShape),"UncappedTorusHoleShape",1,"public Shape",
+  0, 0, 0);
+
+UncappedTorusHoleShape::UncappedTorusHoleShape(double r,
+                               const SphereShape& s1,
+                               const SphereShape& s2):
+_s1(s1),
+_s2(s2),
+_r(r)
+{
+}
+
+UncappedTorusHoleShape*
+UncappedTorusHoleShape::newUncappedTorusHoleShape(double r,
+                                                  const SphereShape&s1,
+                                                  const SphereShape&s2)
+{
+  // if the probe sphere fits between the two spheres, then there
+  // is no need to include this shape
+  SCVector3 A(s1.origin());
+  SCVector3 B(s2.origin());
+  SCVector3 BA = B - A;
+  if (2.0*r <  BA.norm() - s1.radius() - s2.radius()) return 0;
+
+  // check to see if the surface is reentrant
+  double rrs1 = r+s1.radius();
+  double rrs2 = r+s2.radius();
+  SCVector3 R12 = ((SCVector3)s1.origin()) - ((SCVector3)s2.origin());
+  double r12 = sqrt(R12.dot(R12));
+  if (sqrt(rrs1*rrs1-r*r) + sqrt(rrs2*rrs2-r*r) < r12)
+    return new ReentrantUncappedTorusHoleShape(r,s1,s2);
+
+  // otherwise create an ordinary torus hole
+  return new NonreentrantUncappedTorusHoleShape(r,s1,s2);
+}
+
+// Given a node, finds a sphere in the plane of n and the centers
+// of _s1 and _s2 that touches the UncappedTorusHole.  There are two
+// candidates, the one closest to n is chosen.
+void
+UncappedTorusHoleShape::in_plane_sphere(
+    const SCVector3& n,
+    SCVector3& P) const
+{
+  // the center of the sphere is given by the vector equation
+  // P = A + a R(AB) + b U(perp),
+  // where
+  // A is the center of _s1
+  // B is the center of _s2
+  // R(AB) is the vector from A to B, R(AB) = B - A
+  // U(perp) is a unit vect perp to R(AB) and lies in the plane of n, A, and B
+  // The unknown scalars, a and b are given by solving the following
+  // equations:
+  // | P - A | = r(A) + _r, and
+  // | P - B | = r(B) + _r,
+  // which give
+  // | a R(AB) + b U(perp) | = r(A) + _r, and
+  // | (a-1) R(AB) + b U(perp) | = r(B) + _r.
+  // These further simplify to
+  // a^2 r(AB)^2 + b^2 = (r(A)+_r)^2, and
+  // (a-1)^2 r(AB)^2 + b^2 = (r(B)+_r)^2.
+  // Thus,
+  // a = (((r(A)+_r)^2 - (r(B)+_r)^2 )/(2 r(AB)^2)) + 1/2
+  // b^2 = (r(A)+r)^2 - a^2 r(AB)^2
+
+  SCVector3 A = _s1.origin();
+  SCVector3 B = _s2.origin();
+  SCVector3 N = n;
+  SCVector3 R_AB = B - A;
+  SCVector3 R_AN = N - A;
+
+  // vector projection of R_AN onto R_AB
+  SCVector3 P_AN_AB = R_AB * (R_AN.dot(R_AB)/R_AB.dot(R_AB));
+  // the perpendicular vector
+  SCVector3 U_perp = R_AN - P_AN_AB;
+
+  // if |U| is tiny, then any vector perp to AB will do
+  double u = U_perp.dot(U_perp);
+  if (u<1.0e-23) {
+      SCVector3 vtry = R_AB;
+      vtry[0] += 1.0;
+      vtry = vtry - R_AB * (vtry.dot(R_AB)/R_AB.dot(R_AB));
+      if (vtry.dot(vtry) < 1.0e-23) {
+          vtry = R_AB;
+          vtry[1] += 1.0;
+          vtry = vtry - R_AB * (vtry.dot(R_AB)/R_AB.dot(R_AB));
+        }
+      U_perp = vtry;
+    }
+
+  U_perp.normalize();
+  //ExEnv::outn() << "A: "; A.print();
+  //ExEnv::outn() << "U_perp: "; U_perp.print();
+  //ExEnv::outn() << "R_AB: "; R_AB.print();
+
+  double r_AB = sqrt(R_AB.dot(R_AB));
+  double r_A = _s1.radius();
+  double r_B = _s2.radius();
+
+  double r_Ar = r_A + _r;
+  double r_Br = r_B + _r;
+  double a = ((r_Ar*r_Ar - r_Br*r_Br)/(2.0*r_AB*r_AB)) + 0.5;
+  double b = sqrt(r_Ar*r_Ar - a*a*r_AB*r_AB);
+
+  //ExEnv::outn() << scprintf("r_Ar = %f, r_AB = %f\n",r_Ar,r_AB);
+  //ExEnv::outn() << scprintf("a = %f, b = %f\n",a,b);
+
+  P = A + a * R_AB + b * U_perp;
+  //ExEnv::outn() << "a*R_AB: "; (a*R_AB).print();
+  //ExEnv::outn() << "b*U_perp: "; (b*U_perp).print();
+}
+
+void
+UncappedTorusHoleShape::print(ostream&o) const
+{
+  o << indent << "UncappedTorusHoleShape:" << endl;
+  o << incindent;
+  o << indent << "r = " << _r << endl;
+  o << indent << "s1 = ";
+  o << incindent << skipnextindent;
+  _s1.print(o);
+  o << decindent;
+  o << indent << "s2 = ";
+  o << incindent << skipnextindent;
+  _s2.print(o);
+  o << decindent;
+  o << decindent;
+}
+
+void
+UncappedTorusHoleShape::boundingbox(double valuemin, double valuemax,
+                                    SCVector3& p1,
+                                    SCVector3& p2)
+{
+  SCVector3 p11;
+  SCVector3 p12;
+  SCVector3 p21;
+  SCVector3 p22;
+
+  _s1.boundingbox(valuemin,valuemax,p11,p12);
+  _s2.boundingbox(valuemin,valuemax,p21,p22);
+
+  int i;
+  for (i=0; i<3; i++) {
+      if (p11[i] < p21[i]) p1[i] = p11[i];
+      else p1[i] = p21[i];
+      if (p12[i] > p22[i]) p2[i] = p12[i];
+      else p2[i] = p22[i];
+    }
+}
+
+int
+UncappedTorusHoleShape::gradient_implemented() const
+{
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////
+// is in triangle
+
+static int
+is_in_unbounded_triangle(const SCVector3&XP,const SCVector3&AP,const SCVector3&BP)
+{
+  SCVector3 axis = BP.cross(AP);
+
+  SCVector3 BP_perp = BP; BP_perp.rotate(M_PI_2,axis);
+  double u = BP_perp.dot(XP)/BP_perp.dot(AP);
+  if (u<0.0) return 0;
+
+  SCVector3 AP_perp = AP; AP_perp.rotate(M_PI_2,axis);
+  double w = AP_perp.dot(XP)/AP_perp.dot(BP);
+  if (w<0.0) return 0;
+
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////
+// ReentrantUncappedTorusHoleShape
+
+static ClassDesc ReentrantUncappedTorusHoleShape_cd(
+  typeid(ReentrantUncappedTorusHoleShape),"ReentrantUncappedTorusHoleShape",1,"public UncappedTorusHoleShape",
+  0, 0, 0);
+
+ReentrantUncappedTorusHoleShape::ReentrantUncappedTorusHoleShape(double r,
+                                                 const SphereShape& s1,
+                                                 const SphereShape& s2):
+  UncappedTorusHoleShape(r,s1,s2)
+{
+  rAP = r + s1.radius();
+  rBP = r + s2.radius();
+  BA = B() - A();
+
+  // Find the points at the ends of the two cone-like objects, I[0] and I[1].
+  // they are given by:
+  //   I = z BA, where BA = B-A and I is actually IA = I - A
+  //   r^2 = PI.PI, where PI = PA-I and P is the center of a probe sphere
+  // this gives
+  //  z^2 BA.BA - 2z PA.BA + PA.PA - r^2 = 0
+
+  SCVector3 arbitrary; 
+  arbitrary[0] = arbitrary[1] = arbitrary[2] = 0.0;
+  SCVector3 P;
+  in_plane_sphere(arbitrary,P);
+  SCVector3 PA = P - A();
+
+  double a = BA.dot(BA);
+  double minus_b = 2.0 * PA.dot(BA);
+  double c = PA.dot(PA) - r*r;
+  double b2m4ac = minus_b*minus_b - 4*a*c;
+  double sb2m4ac;
+  if (b2m4ac >= 0.0) {
+      sb2m4ac = sqrt(b2m4ac);
+    }
+  else if (b2m4ac > -1.0e-10) {
+      sb2m4ac = 0.0;
+    }
+  else {
+      ExEnv::errn() << "ReentrantUncappedTorusHoleShape:: imaginary point" << endl;
+      abort();
+    }
+  double zA = (minus_b - sb2m4ac)/(2.0*a);
+  double zB = (minus_b + sb2m4ac)/(2.0*a);
+  I[0] = BA * zA + A();
+  I[1] = BA * zB + A();
+}
+ReentrantUncappedTorusHoleShape::~ReentrantUncappedTorusHoleShape()
+{
+}
+int
+ReentrantUncappedTorusHoleShape::
+  is_outside(const SCVector3&X) const
+{
+  SCVector3 Xv(X);
+
+  SCVector3 P;
+  in_plane_sphere(X,P);
+  SCVector3 XP = Xv - P;
+  double rXP = XP.norm();
+  if (rXP > rAP || rXP > rBP) return 1;
+
+  SCVector3 AP = A() - P;
+  SCVector3 BP = B() - P;
+
+  if (!is_in_unbounded_triangle(XP,AP,BP)) return 1;
+
+  if (rXP < radius()) {
+      return 1;
+    }
+
+  return 0;
+}
+double
+ReentrantUncappedTorusHoleShape::
+  distance_to_surface(const SCVector3&X,SCVector3*grad) const
+{
+  SCVector3 Xv(X);
+
+  SCVector3 P;
+  in_plane_sphere(X,P);
+  SCVector3 XP = Xv - P;
+  double rXP = XP.norm();
+  if (rXP > rAP || rXP > rBP) return shape_infinity;
+
+  SCVector3 AP = A() - P;
+  SCVector3 BP = B() - P;
+
+  if (!is_in_unbounded_triangle(XP,AP,BP)) return shape_infinity;
+
+  SCVector3 I1P = I[0] - P;
+  SCVector3 I2P = I[1] - P;
+
+  if (!is_in_unbounded_triangle(XP,I1P,I2P)) {
+      if (rXP < radius()) {
+          if (grad) {
+              SCVector3 unit(XP);
+              unit.normalize();
+              *grad = unit;
+            }
+          return radius() - rXP;
+        }
+      else return -1.0;
+    }
+
+  return closest_distance(Xv,(SCVector3*)I,2,grad);
+}
+
+int
+ReentrantUncappedTorusHoleShape::gradient_implemented() const
+{
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////
+// NonreentrantUncappedTorusHoleShape
+
+static ClassDesc NonreentrantUncappedTorusHoleShape_cd(
+  typeid(NonreentrantUncappedTorusHoleShape),"NonreentrantUncappedTorusHoleShape",1,"public UncappedTorusHoleShape",
+  0, 0, 0);
+
+NonreentrantUncappedTorusHoleShape::
+  NonreentrantUncappedTorusHoleShape(double r,
+                                     const SphereShape& s1,
+                                     const SphereShape& s2):
+  UncappedTorusHoleShape(r,s1,s2)
+{
+  rAP = r + s1.radius();
+  rBP = r + s2.radius();
+  BA = B() - A();
+}
+NonreentrantUncappedTorusHoleShape::~NonreentrantUncappedTorusHoleShape()
+{
+}
+double NonreentrantUncappedTorusHoleShape::
+  distance_to_surface(const SCVector3&X,SCVector3* grad) const
+{
+  SCVector3 Xv(X);
+
+  SCVector3 P;
+  in_plane_sphere(X,P);
+  SCVector3 PX = P - Xv;
+  double rPX = PX.norm();
+  if (rPX > rAP || rPX > rBP) return shape_infinity;
+
+  SCVector3 PA = P - A();
+  SCVector3 XA = Xv - A();
+
+  SCVector3 axis = BA.cross(PA);
+
+  SCVector3 BA_perp = BA; BA_perp.rotate(M_PI_2,axis);
+  double u = BA_perp.dot(XA)/BA_perp.dot(PA);
+  if (u<0.0 || u>1.0) return shape_infinity;
+
+  SCVector3 PA_perp = PA; PA_perp.rotate(M_PI_2,axis);
+  double w = PA_perp.dot(XA)/PA_perp.dot(BA);
+  if (w<0.0 || w>1.0) return shape_infinity;
+
+  double uw = u+w;
+  if (uw<0.0 || uw>1.0) return shape_infinity;
+
+  if (rPX < radius()) {
+      if (grad) {
+          SCVector3 unit(PX);
+          unit.normalize();
+          *grad = unit;
+        }
+      return radius() - rPX;
+    }
+
+  return -1;
+}
+
+int
+NonreentrantUncappedTorusHoleShape::gradient_implemented() const
+{
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////
+// Uncapped5SphereExclusionShape
+
+static ClassDesc Uncapped5SphereExclusionShape_cd(
+  typeid(Uncapped5SphereExclusionShape),"Uncapped5SphereExclusionShape",1,"public Shape",
+  0, 0, 0);
+
+Uncapped5SphereExclusionShape*
+Uncapped5SphereExclusionShape::
+  newUncapped5SphereExclusionShape(double r,
+                                   const SphereShape& s1,
+                                   const SphereShape& s2,
+                                   const SphereShape& s3)
+{
+  Uncapped5SphereExclusionShape* s =
+    new Uncapped5SphereExclusionShape(r,s1,s2,s3);
+  if (s->solution_exists()) {
+      return s;
+    }
+  delete s;
+  return 0;
+}
+static int verbose = 0;
+Uncapped5SphereExclusionShape::
+  Uncapped5SphereExclusionShape(double radius,
+                                const SphereShape&s1,
+                                const SphereShape&s2,
+                                const SphereShape&s3):
+  _s1(s1),
+  _s2(s2),
+  _s3(s3),
+  _r(radius)
+{
+  double rAr = rA() + r();
+  double rAr2 = rAr*rAr;
+  double rBr = rB() + r();
+  double rBr2 = rBr*rBr;
+  double rCr = rC() + r();
+  double rCr2 = rCr*rCr;
+  double A2 = A().dot(A());
+  double B2 = B().dot(B());
+  double C2 = C().dot(C());
+  SCVector3 BA = B()-A();
+  double DdotBA = 0.5*(rAr2 - rBr2 + B2 - A2);
+  double DAdotBA = DdotBA - A().dot(BA);
+  double BA2 = BA.dot(BA);
+  SCVector3 CA = C() - A();
+  double CAdotBA = CA.dot(BA);
+  SCVector3 CA_perpBA = CA - (CAdotBA/BA2)*BA;
+  double CA_perpBA2 = CA_perpBA.dot(CA_perpBA);
+  if (CA_perpBA2 < 1.0e-23) {
+      _solution_exists = 0;
+      return;
+    }
+  double DdotCA_perpBA = 0.5*(rAr2 - rCr2 + C2 - A2)
+    - CAdotBA*DdotBA/BA2;
+  double DAdotCA_perpBA = DdotCA_perpBA - A().dot(CA_perpBA);
+  double rAt2 = DAdotBA*DAdotBA/BA2 + DAdotCA_perpBA*DAdotCA_perpBA/CA_perpBA2;
+  double h2 = rAr2 - rAt2;
+  if (h2 <= 0.0) {
+      _solution_exists = 0;
+      return;
+    }
+  _solution_exists = 1;
+
+  double h = sqrt(h2);
+  if (h<r()) {
+      _reentrant = 1;
+      //ExEnv::outn() << "WARNING: throwing out reentrant shape" << endl;
+      //_solution_exists = 0;
+      //return;
+    }
+  else {
+      _reentrant = 0;
+      //ExEnv::outn() << "WARNING: throwing out nonreentrant shape" << endl;
+      //_solution_exists = 0;
+      //return;
+    }
+
+  // The projection of D into the ABC plane
+  SCVector3 MA = (DAdotBA/BA2)*BA + (DAdotCA_perpBA/CA_perpBA2)*CA_perpBA;
+  M = MA + A();
+  SCVector3 BAxCA = BA.cross(CA);
+  D[0] = M + h * BAxCA * ( 1.0/BAxCA.norm() );
+  D[1] = M - h * BAxCA * ( 1.0/BAxCA.norm() );
+
+  // The projection of D into the ABC plane, M, does not lie in the
+  // ABC triangle, then this shape must be treated carefully and it
+  // will be designated as folded.
+  SCVector3 MC = M - C();
+  if (!(is_in_unbounded_triangle(MA, BA, CA)
+        &&is_in_unbounded_triangle(MC, B() - C(), A() - C()))) {
+      _folded = 1;
+      SCVector3 MB = M - B();
+      double MA2 = MA.dot(MA);
+      double MB2 = MB.dot(MB);
+      double MC2 = MC.dot(MC);
+      if (MA2 < MB2) {
+          F1 = A();
+          if (MB2 < MC2) F2 = B();
+          else F2 = C();
+        }
+      else {
+          F1 = B();
+          if (MA2 < MC2) F2 = A();
+          else F2 = C();
+        }
+    }
+  else _folded = 0;
+  
+  //ExEnv::outn() << scprintf("r = %14.8f, h = %14.8f\n",r(),h);
+  //M.print();
+  //D[0].print();
+  //D[1].print();
+  //A().print();
+  //B().print();
+  //C().print();
+
+  int i;
+  for (i=0; i<2; i++) {
+      SCVector3 AD = A() - D[i];
+      SCVector3 BD = B() - D[i];
+      SCVector3 CD = C() - D[i];
+      BDxCD[i] = BD.cross(CD);
+      BDxCDdotAD[i] = BDxCD[i].dot(AD);
+      CDxAD[i] = CD.cross(AD);
+      CDxADdotBD[i] = CDxAD[i].dot(BD);
+      ADxBD[i] = AD.cross(BD);
+      ADxBDdotCD[i] = ADxBD[i].dot(CD);
+    }
+
+  for (i=0; i<2; i++) MD[i] = M - D[i];
+
+  // reentrant surfaces need a whole bunch more to be able to compute
+  // the distance to the surface
+  if (_reentrant) {
+      int i;
+      double rMD = MD[0].norm(); // this is the same as rMD[1]
+      theta_intersect = M_PI_2 - asin(rMD/r());
+      r_intersect = r() * sin(theta_intersect);
+
+      {
+        // Does the circle of intersection intersect with BA?
+        SCVector3 MA = M - A();
+        SCVector3 uBA = B() - A(); uBA.normalize();
+        SCVector3 XA = uBA * MA.dot(uBA);
+        SCVector3 XM = XA - MA;
+        double rXM2 = XM.dot(XM);
+        double d_intersect_from_x2 = r_intersect*r_intersect - rXM2;
+        if (d_intersect_from_x2 > 0.0) {
+            _intersects_AB = 1;
+            double tmp = sqrt(d_intersect_from_x2);
+            double d_intersect_from_x[2];
+            d_intersect_from_x[0] = tmp;
+            d_intersect_from_x[1] = -tmp;
+            for (i=0; i<2; i++) {
+                for (int j=0; j<2; j++) {
+                    IABD[i][j] = XM + d_intersect_from_x[j]*uBA + MD[i];
+                  }
+              }
+          }
+        else _intersects_AB = 0;
+      }
+
+      {
+        // Does the circle of intersection intersect with BC?
+        SCVector3 MC = M - C();
+        SCVector3 uBC = B() - C(); uBC.normalize();
+        SCVector3 XC = uBC * MC.dot(uBC);
+        SCVector3 XM = XC - MC;
+        double rXM2 = XM.dot(XM);
+        double d_intersect_from_x2 = r_intersect*r_intersect - rXM2;
+        if (d_intersect_from_x2 > 0.0) {
+            _intersects_BC = 1;
+            double tmp = sqrt(d_intersect_from_x2);
+            double d_intersect_from_x[2];
+            d_intersect_from_x[0] = tmp;
+            d_intersect_from_x[1] = -tmp;
+            for (i=0; i<2; i++) {
+                for (int j=0; j<2; j++) {
+                    IBCD[i][j] = XM + d_intersect_from_x[j]*uBC + MD[i];
+                  }
+              }
+          }
+        else _intersects_BC = 0;
+      }
+
+      {
+        // Does the circle of intersection intersect with CA?
+        SCVector3 MA = M - A();
+        SCVector3 uCA = C() - A(); uCA.normalize();
+        SCVector3 XA = uCA * MA.dot(uCA);
+        SCVector3 XM = XA - MA;
+        double rXM2 = XM.dot(XM);
+        double d_intersect_from_x2 = r_intersect*r_intersect - rXM2;
+        if (d_intersect_from_x2 > 0.0) {
+            _intersects_CA = 1;
+            double tmp = sqrt(d_intersect_from_x2);
+            double d_intersect_from_x[2];
+            d_intersect_from_x[0] = tmp;
+            d_intersect_from_x[1] = -tmp;
+            for (i=0; i<2; i++) {
+                for (int j=0; j<2; j++) {
+                    ICAD[i][j] = XM + d_intersect_from_x[j]*uCA + MD[i];
+                  }
+              }
+          }
+        else _intersects_CA = 0;
+      }
+
+    }
+
+#if 0 // test code
+  ExEnv::outn() << "Uncapped5SphereExclusionShape: running some tests" << endl;
+  verbose = 1;
+
+  FILE* testout = fopen("testout.vect", "w");
+
+  const double scalefactor_inc = 0.05;
+  const double start = -10.0;
+  const double end = 10.0;
+
+  SCVector3 middle = (1.0/3.0)*(A()+B()+C());
+
+  int nlines = 1;
+  int nvert = (int) ( (end-start) / scalefactor_inc);
+  int ncolor = nvert;
+
+  fprintf(testout, "VECT\n%d %d %d\n", nlines, nvert, ncolor);
+
+  fprintf(testout, "%d\n", nvert);
+  fprintf(testout, "%d\n", ncolor);
+
+  double scalefactor = start;
+  for (int ii = 0; ii<nvert; ii++) {
+      SCVector3 position = (D[0] - middle) * scalefactor + middle;
+      double d = distance_to_surface(position);
+      fprintf(testout, "%f %f %f # value = %f\n",
+              position[0], position[1], position[2], d);
+      scalefactor += scalefactor_inc;
+    }
+  scalefactor = start;
+  for (ii = 0; ii<nvert; ii++) {
+      SCVector3 position = (D[0] - middle) * scalefactor + middle;
+      double d = distance_to_surface(position);
+      ExEnv::outn() << scprintf("d = %f\n", d);
+      if (d<0.0) fprintf(testout,"1.0 0.0 0.0 0.5\n");
+      else fprintf(testout,"0.0 0.0 1.0 0.5\n");
+      scalefactor += scalefactor_inc;
+    }
+
+  fclose(testout);
+  ExEnv::outn() << "testout.vect written" << endl;
+
+  verbose = 0;
+#endif // test code
+  
+}
+int
+Uncapped5SphereExclusionShape::is_outside(const SCVector3&X) const
+{
+  SCVector3 Xv(X);
+
+  if (verbose) ExEnv::outn() << scprintf("point %14.8f %14.8f %14.8f\n",X(0),X(1),X(2));
+
+  // The folded case isn't handled correctly here, so use
+  // the less efficient distance_to_surface routine.
+  if (_folded) {
+      return distance_to_surface(X) >= 0.0;
+    }
+
+  for (int i=0; i<2; i++) {
+      SCVector3 XD = Xv - D[i];
+      double rXD = XD.norm();
+      if (rXD <= r()) return 1;
+      double u = BDxCD[i].dot(XD)/BDxCDdotAD[i];
+      if (u <= 0.0) return 1;
+      double v = CDxAD[i].dot(XD)/CDxADdotBD[i];
+      if (v <= 0.0) return 1;
+      double w = ADxBD[i].dot(XD)/ADxBDdotCD[i];
+      if (w <= 0.0) return 1;
+    }
+
+  if (verbose) ExEnv::outn() << "is_inside" << endl;
+
+  return 0;
+}
+static int
+is_contained_in_unbounded_pyramid(SCVector3 XD,
+                                  SCVector3 AD,
+                                  SCVector3 BD,
+                                  SCVector3 CD)
+{
+  SCVector3 BDxCD = BD.cross(CD);
+  SCVector3 CDxAD = CD.cross(AD);
+  SCVector3 ADxBD = AD.cross(BD);
+  double u = BDxCD.dot(XD)/BDxCD.dot(AD);
+  if (u <= 0.0) return 0;
+  double v = CDxAD.dot(XD)/CDxAD.dot(BD);
+  if (v <= 0.0) return 0;
+  double w = ADxBD.dot(XD)/ADxBD.dot(CD);
+  if (w <= 0.0) return 0;
+  return 1;
+}
+double
+Uncapped5SphereExclusionShape::
+  distance_to_surface(const SCVector3&X,SCVector3*grad) const
+{
+  SCVector3 Xv(X);
+
+  // Find out if I'm on the D[0] side or the D[1] side of the ABC plane.
+  int side;
+  SCVector3 XM = Xv - M;
+  if (MD[0].dot(XM) > 0.0) side = 1;
+  else side = 0;
+
+  if (verbose) {
+      ExEnv::outn() << scprintf("distance_to_surface: folded = %d, side = %d\n",
+                       _folded, side);
+      ExEnv::outn() << "XM = "; XM.print();
+      ExEnv::outn() << "MD[0] = "; MD[0].print();
+      ExEnv::outn() << "MD[0].dot(XM) = " << MD[0].dot(XM) << endl;
+    }
+
+  SCVector3 XD = Xv - D[side];
+  double u = BDxCD[side].dot(XD)/BDxCDdotAD[side];
+  if (verbose) ExEnv::outn() << scprintf("u = %14.8f\n", u);
+  if (u <= 0.0) return shape_infinity;
+  double v = CDxAD[side].dot(XD)/CDxADdotBD[side];
+  if (verbose) ExEnv::outn() << scprintf("v = %14.8f\n", v);
+  if (v <= 0.0) return shape_infinity;
+  double w = ADxBD[side].dot(XD)/ADxBDdotCD[side];
+  if (verbose) ExEnv::outn() << scprintf("w = %14.8f\n", w);
+  if (w <= 0.0) return shape_infinity;
+  double rXD = XD.norm();
+  if (verbose) ExEnv::outn() << scprintf("r() - rXD = %14.8f\n", r() - rXD);
+  if (rXD <= r()) {
+      if (!_reentrant) return r() - rXD;
+      // this shape is reentrant
+      // make sure that we're on the right side
+      if ((side == 1) || (u + v + w <= 1.0)) {
+          // see if we're outside the cone that intersects
+          // the opposite sphere
+          double angle = acos(fabs(XD.dot(MD[side]))
+                              /(XD.norm()*MD[side].norm()));
+          if (angle >= theta_intersect) {
+              if (grad) {
+                  *grad = (-1.0/rXD)*XD;
+                }
+              return r() - rXD;
+            }
+          if (_intersects_AB
+              &&is_contained_in_unbounded_pyramid(XD,
+                                                  MD[side],
+                                                  IABD[side][0],
+                                                  IABD[side][1])) {
+              //ExEnv::outn() << scprintf("XD: "); XD.print();
+              //ExEnv::outn() << scprintf("MD[%d]: ",i); MD[i].print();
+              //ExEnv::outn() << scprintf("IABD[%d][0]: ",i); IABD[i][0].print();
+              //ExEnv::outn() << scprintf("IABD[%d][1]: ",i); IABD[i][1].print();
+              return closest_distance(XD,(SCVector3*)IABD[side],2,grad);
+            }
+          if (_intersects_BC
+              &&is_contained_in_unbounded_pyramid(XD,
+                                                  MD[side],
+                                                  IBCD[side][0],
+                                                  IBCD[side][1])) {
+              return closest_distance(XD,(SCVector3*)IBCD[side],2,grad);
+            }
+          if (_intersects_CA
+              &&is_contained_in_unbounded_pyramid(XD,
+                                                  MD[side],
+                                                  ICAD[side][0],
+                                                  ICAD[side][1])) {
+              return closest_distance(XD,(SCVector3*)ICAD[side],2,grad);
+            }
+          // at this point we are closest to the ring formed
+          // by the intersection of the two probe spheres
+          double distance_to_plane;
+          double distance_to_ring_in_plane;
+          double MDnorm = MD[side].norm();
+          SCVector3 XM = XD - MD[side];
+          SCVector3 XM_in_plane;
+          if (MDnorm<1.0e-16) {
+              distance_to_plane = 0.0;
+              XM_in_plane = XD;
+            }
+          else {
+              distance_to_plane = XM.dot(MD[side])/MDnorm;
+              XM_in_plane = XM - (distance_to_plane/MDnorm)*MD[side];
+            }
+          if (grad) {
+              double XM_in_plane_norm = XM_in_plane.norm();
+              if (XM_in_plane_norm < 1.e-8) {
+                  // the grad points along MD
+                  double scale = - distance_to_plane
+                         /(MDnorm*sqrt(r_intersect*r_intersect
+                                       + distance_to_plane*distance_to_plane));
+                  *grad = MD[side] * scale;
+                }
+              else {
+                  SCVector3 point_on_ring;
+                  point_on_ring = XM_in_plane
+                                * (r_intersect/XM_in_plane_norm) + M;
+                  SCVector3 gradv = Xv - point_on_ring;
+                  gradv.normalize();
+                  *grad = gradv;
+                }
+            }
+          distance_to_ring_in_plane =
+                         r_intersect - sqrt(XM_in_plane.dot(XM_in_plane));
+          return sqrt(distance_to_ring_in_plane*distance_to_ring_in_plane
+                      +distance_to_plane*distance_to_plane);
+        }
+    }
+
+  if (verbose) ExEnv::outn() << "returning -1.0" << endl;
+  return -1.0;
+}
+
+void
+Uncapped5SphereExclusionShape::boundingbox(double valuemin, double valuemax,
+                                           SCVector3& p1,
+                                           SCVector3& p2)
+{
+  SCVector3 p11;
+  SCVector3 p12;
+  SCVector3 p21;
+  SCVector3 p22;
+  SCVector3 p31;
+  SCVector3 p32;
+
+  _s1.boundingbox(valuemin,valuemax,p11,p12);
+  _s2.boundingbox(valuemin,valuemax,p21,p22);
+  _s3.boundingbox(valuemin,valuemax,p31,p32);
+
+  int i;
+  for (i=0; i<3; i++) {
+      if (p11[i] < p21[i]) p1[i] = p11[i];
+      else p1[i] = p21[i];
+      if (p31[i] < p1[i]) p1[i] = p31[i];
+      if (p12[i] > p22[i]) p2[i] = p12[i];
+      else p2[i] = p22[i];
+      if (p32[i] > p2[i]) p2[i] = p32[i];
+    }
+}
+
+int
+Uncapped5SphereExclusionShape::gradient_implemented() const
+{
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////
+// Unionshape
+
+static ClassDesc UnionShape_cd(
+  typeid(UnionShape),"UnionShape",1,"public Shape",
+  0, 0, 0);
+
+UnionShape::UnionShape()
+{
+}
+
+UnionShape::~UnionShape()
+{
+}
+
+void
+UnionShape::add_shape(Ref<Shape> s)
+{
+  _shapes.insert(s);
+}
+
+// NOTE: this underestimates the distance to the surface when
+//inside the surface
+double
+UnionShape::distance_to_surface(const SCVector3&p,SCVector3* grad) const
+{
+  std::set<Ref<Shape> >::const_iterator imin = _shapes.begin();
+  if (imin == _shapes.end()) return 0.0;
+  double min = (*imin)->distance_to_surface(p);
+  for (std::set<Ref<Shape> >::const_iterator i=imin; i!=_shapes.end(); i++) {
+      double d = (*i)->distance_to_surface(p);
+      if (min <= 0.0) {
+          if (d < 0.0 && d > min) { min = d; imin = i; }
+        }
+      else {
+          if (min > d) { min = d; imin = i; }
+        }
+    }
+
+  if (grad) {
+      (*imin)->distance_to_surface(p,grad);
+    }
+  return min;
+}
+
+int
+UnionShape::is_outside(const SCVector3&p) const
+{
+  for (std::set<Ref<Shape> >::const_iterator i=_shapes.begin();
+       i!=_shapes.end(); i++) {
+      if (!(*i)->is_outside(p)) return 0;
+    }
+
+  return 1;
+}
+
+void
+UnionShape::boundingbox(double valuemin, double valuemax,
+                        SCVector3& p1,
+                        SCVector3& p2)
+{
+  if (_shapes.begin() == _shapes.end()) {
+      for (int i=0; i<3; i++) p1[i] = p2[i] = 0.0;
+      return;
+    }
+  
+  SCVector3 pt1;
+  SCVector3 pt2;
+  
+  std::set<Ref<Shape> >::iterator j = _shapes.begin();
+  int i;
+  (*j)->boundingbox(valuemin,valuemax,p1,p2);
+  for (j++; j!=_shapes.end(); j++) {
+      (*j)->boundingbox(valuemin,valuemax,pt1,pt2);
+      for (i=0; i<3; i++) {
+          if (pt1[i] < p1[i]) p1[i] = pt1[i];
+          if (pt2[i] > p2[i]) p2[i] = pt2[i];
+        }
+    }
+}
+
+int
+UnionShape::gradient_implemented() const
+{
+  for (std::set<Ref<Shape> >::const_iterator j=_shapes.begin();
+       j!=_shapes.end(); j++) {
+      if (!(*j)->gradient_implemented()) return 0;
+    }
+  return 1;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/shape.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/shape.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/shape.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,252 @@
+//
+// shape.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_shape_h
+#define _math_isosurf_shape_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <set>
+
+#include <math/isosurf/volume.h>
+#include <math/scmat/matrix.h>
+#include <math/scmat/vector3.h>
+
+namespace sc {
+
+/** A Shape is a Volume represents an 3D solid.  The value of the Shape at
+each point in space is the distance to the surface.  The distance is
+negative if the point is inside the solid.  For Shape specializations that
+cannot compute the distance to the surface, the value will be 1.0 outside
+and -1.0 inside the solid. */
+class Shape: public Volume {
+  public:
+    Shape();
+    Shape(const Ref<KeyVal>&keyval);
+    virtual double distance_to_surface(const SCVector3&r,
+                                       SCVector3*grad=0) const = 0;
+    virtual int is_outside(const SCVector3&r) const;
+    virtual ~Shape();
+    void compute();
+    void interpolate(const SCVector3& p1,
+                     const SCVector3& p2,
+                     double val,
+                     SCVector3& result);
+
+    int value_implemented() const;
+};
+
+
+
+class SphereShape: public Shape {
+  private:
+    SCVector3 _origin;
+    double _radius;
+  public:
+    SphereShape(const SCVector3&,double);
+    SphereShape(const Ref<KeyVal>&);
+    SphereShape(const SphereShape&);
+    ~SphereShape();
+    void boundingbox(double minvalue, double maxvalue,
+                     SCVector3& p1, SCVector3&p2);
+    double radius() const { return _radius; }
+    const SCVector3& origin() const { return _origin; }
+    double distance_to_surface(const SCVector3&r,SCVector3*grad=0) const;
+    void print(std::ostream&o=ExEnv::out0()) const;
+
+    // these are used to update the parameters describing the sphere
+    double radius(double r);
+    const SCVector3& origin(const SCVector3& o);
+
+    int gradient_implemented() const;
+};
+
+inline double
+SphereShape::radius(double r)
+{
+  obsolete();
+  return _radius = r;
+}
+
+inline const SCVector3&
+SphereShape::origin(const SCVector3& o)
+{
+  obsolete();
+  _origin = o;
+  return _origin;
+}
+
+class UncappedTorusHoleShape: public Shape
+{
+  private:
+    SphereShape _s1;
+    SphereShape _s2;
+    double _r;
+  protected:
+    void in_plane_sphere(const SCVector3& point,
+                         SCVector3& origin) const;
+    UncappedTorusHoleShape(double r,const SphereShape&,const SphereShape&);
+  public:
+    static UncappedTorusHoleShape*
+    newUncappedTorusHoleShape(double r,
+                              const SphereShape&,
+                              const SphereShape&);
+    inline ~UncappedTorusHoleShape() {};
+    inline const SphereShape& sphere(int i) const { return (i?_s2:_s1); };
+    inline const SCVector3 A() const { SCVector3 v(_s1.origin()); return v; }
+    inline const SCVector3 B() const { SCVector3 v(_s2.origin()); return v; }
+    inline double radius() const { return _r; };
+    void print(std::ostream&o=ExEnv::out0()) const;
+    void boundingbox(double valuemin, double valuemax,
+                     SCVector3& p1, SCVector3&p2);
+
+    int gradient_implemented() const;
+};
+
+class NonreentrantUncappedTorusHoleShape: public UncappedTorusHoleShape
+{
+  private:
+    double rAP;
+    double rBP;
+    SCVector3 BA;
+  public:
+    NonreentrantUncappedTorusHoleShape(double r,
+                                       const SphereShape&,
+                                       const SphereShape&);
+    ~NonreentrantUncappedTorusHoleShape();
+    double distance_to_surface(const SCVector3&r,SCVector3*grad=0) const;
+
+    int gradient_implemented() const;
+};
+
+class ReentrantUncappedTorusHoleShape: public UncappedTorusHoleShape
+{
+  private:
+    double rAP;
+    double rBP;
+    SCVector3 BA;
+    SCVector3 I[2]; // the intersect points
+  public:
+    ReentrantUncappedTorusHoleShape(double r,
+                                    const SphereShape&,
+                                    const SphereShape&);
+    ~ReentrantUncappedTorusHoleShape();
+    int is_outside(const SCVector3&r) const;
+    double distance_to_surface(const SCVector3&r,SCVector3*grad=0) const;
+
+    int gradient_implemented() const;
+};
+
+class Uncapped5SphereExclusionShape: public Shape
+{
+  private:
+    int _solution_exists;
+    int _reentrant;
+    int _folded;
+    SphereShape _s1;
+    SphereShape _s2;
+    SphereShape _s3;
+    SCVector3 D[2];
+    double BDxCDdotAD[2];
+    SCVector3 BDxCD[2];
+    double CDxADdotBD[2];
+    SCVector3 CDxAD[2];
+    double ADxBDdotCD[2];
+    SCVector3 ADxBD[2];
+    double _r;
+
+    // these are needed for folded shapes
+    // F1 and F2 are the two points of A, B, and C that are closed to M
+    SCVector3 F1;
+    SCVector3 F2;
+    
+    // these are needed for reentrant surfaces to compute distances
+    SCVector3 M;   // projection of D onto ABC plane
+    SCVector3 MD[2];  // M - D 
+    double theta_intersect; // angle M - D - intersect_point
+    double r_intersect; // the radius of the intersect circle
+    int _intersects_AB;
+    SCVector3 IABD[2][2];
+    int _intersects_BC;
+    SCVector3 IBCD[2][2];
+    int _intersects_CA;
+    SCVector3 ICAD[2][2];
+    
+  protected:
+    Uncapped5SphereExclusionShape(double r,
+                                  const SphereShape&,
+                                  const SphereShape&,
+                                  const SphereShape&);
+  public:
+    static Uncapped5SphereExclusionShape*
+    newUncapped5SphereExclusionShape(double r,
+                                     const SphereShape&,
+                                     const SphereShape&,
+                                     const SphereShape&);
+    inline ~Uncapped5SphereExclusionShape() {};
+    inline const SCVector3 A() const { SCVector3 v(_s1.origin()); return v; }
+    inline const SCVector3 B() const { SCVector3 v(_s2.origin()); return v; }
+    inline const SCVector3 C() const { SCVector3 v(_s3.origin()); return v; }
+    inline double rA() const { return _s1.radius(); };
+    inline double rB() const { return _s2.radius(); };
+    inline double rC() const { return _s3.radius(); };
+    inline double r() const { return _r; };
+    inline int solution_exists() const { return _solution_exists; };
+    double distance_to_surface(const SCVector3&r,SCVector3*grad=0) const;
+    int is_outside(const SCVector3&) const;
+    void boundingbox(double valuemin, double valuemax,
+                     SCVector3& p1, SCVector3&p2);
+
+    int gradient_implemented() const;
+};
+
+/** A UnionShape is volume enclosed by a set of Shape's. */
+class UnionShape: public Shape {
+  protected:
+    std::set<Ref<Shape> > _shapes;
+  public:
+    void add_shape(Ref<Shape>);
+    UnionShape();
+    ~UnionShape();
+    double distance_to_surface(const SCVector3&r,SCVector3*grad=0) const;
+    int is_outside(const SCVector3&r) const;
+    void boundingbox(double valuemin, double valuemax,
+                     SCVector3& p1, SCVector3& p2);
+
+    int gradient_implemented() const;
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/surf.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/surf.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/surf.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,970 @@
+//
+// surf.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/scmat/matrix.h>
+#include <math/scmat/vector3.h>
+#include <math/isosurf/surf.h>
+#include <math/isosurf/isosurf.h>
+#include <util/render/polygons.h>
+
+using namespace std;
+using namespace sc;
+
+#ifndef WRITE_OOGL
+#define WRITE_OOGL 1
+#endif
+
+#if WRITE_OOGL
+#include <util/render/oogl.h>
+#endif
+
+/////////////////////////////////////////////////////////////////////////
+// TriangulatedSurface
+static ClassDesc TriangulatedSurface_cd(
+  typeid(TriangulatedSurface),"TriangulatedSurface",1,"public DescribedClass",
+  create<TriangulatedSurface>, create<TriangulatedSurface>, 0);
+
+TriangulatedSurface::TriangulatedSurface():
+  _verbose(0),
+  _debug(0),
+  _triangle_vertex(0),
+  _triangle_edge(0),
+  _edge_vertex(0),
+  _integrator(new GaussTriangleIntegrator(1))
+{
+  clear();
+}
+
+TriangulatedSurface::TriangulatedSurface(const Ref<KeyVal>& keyval):
+  _triangle_vertex(0),
+  _triangle_edge(0),
+  _edge_vertex(0)
+{
+  _verbose = keyval->booleanvalue("verbose");
+  _debug = keyval->booleanvalue("debug");
+  Ref<TriangleIntegrator> triint;
+  triint << keyval->describedclassvalue("integrator");
+  if (triint.null()) {
+      triint = new GaussTriangleIntegrator(1);
+    }
+  set_integrator(triint);
+  triint << keyval->describedclassvalue("fast_integrator");
+  set_fast_integrator(triint);
+  triint << keyval->describedclassvalue("accurate_integrator");
+  set_accurate_integrator(triint);
+  clear();
+}
+
+TriangulatedSurface::~TriangulatedSurface()
+{
+  clear();
+}
+
+void
+TriangulatedSurface::topology_info(ostream&o)
+{
+  topology_info(nvertex(), nedge(), ntriangle(), o);
+}
+
+void
+TriangulatedSurface::topology_info(int v, int e, int t, ostream&o)
+{
+  // Given v vertices i expect 2*v - 4*n_surface triangles
+  // and 3*v - 6*n_surface edges
+  o << indent
+    << scprintf("n_vertex = %d, n_edge = %d, n_triangle = %d:",
+                v, e, t)
+    << endl;
+  int nsurf_e = ((3*v - e)%6 == 0)? (3*v - e)/6 : -1;
+  int nsurf_t = ((2*v - t)%4 == 0)? (2*v - t)/4 : -1;
+  if ((nsurf_e!=-1) && (nsurf_e == nsurf_t)) {
+      o << indent
+        << scprintf("  this is consistent with n_closed_surface - n_hole = %d",
+                    nsurf_e)
+        << endl;
+    }
+  else {
+      o << indent
+        << scprintf("  this implies that some surfaces are not closed")
+        << endl;
+    }
+}
+
+void
+TriangulatedSurface::set_integrator(const Ref<TriangleIntegrator>& i)
+{
+  _integrator = i;
+}
+
+void
+TriangulatedSurface::set_fast_integrator(const Ref<TriangleIntegrator>& i)
+{
+  _fast_integrator = i;
+}
+
+void
+TriangulatedSurface::set_accurate_integrator(const Ref<TriangleIntegrator>& i)
+{
+  _accurate_integrator = i;
+}
+
+Ref<TriangleIntegrator>
+TriangulatedSurface::integrator(int)
+{
+  // currently the argument, the integer index of the triangle, is ignored
+  return _integrator;
+}
+
+Ref<TriangleIntegrator>
+TriangulatedSurface::fast_integrator(int)
+{
+  // currently the argument, the integer index of the triangle, is ignored
+  return _fast_integrator.null()?_integrator:_fast_integrator;
+}
+
+Ref<TriangleIntegrator>
+TriangulatedSurface::accurate_integrator(int)
+{
+  // currently the argument, the integer index of the triangle, is ignored
+  return _accurate_integrator.null()?_integrator:_accurate_integrator;
+}
+
+void
+TriangulatedSurface::clear_int_arrays()
+{
+  if (_triangle_vertex) {
+      for (int i=0; i<_triangles.size(); i++) {
+          delete[] _triangle_vertex[i];
+        }
+      delete[] _triangle_vertex;
+    }
+  _triangle_vertex = 0;
+
+  if (_triangle_edge) {
+      for (int i=0; i<_triangles.size(); i++) {
+          delete[] _triangle_edge[i];
+        }
+      delete[] _triangle_edge;
+    }
+  _triangle_edge = 0;
+
+  if (_edge_vertex) {
+      for (int i=0; i<_edges.size(); i++) {
+          delete[] _edge_vertex[i];
+        }
+      delete[] _edge_vertex;
+    }
+  _edge_vertex = 0;
+
+  _completed_surface = 0;
+}
+
+void
+TriangulatedSurface::clear()
+{
+  _completed_surface = 0;
+
+  clear_int_arrays();
+
+  _have_values = 0;
+  _values.clear();
+
+  _vertices.clear();
+  _edges.clear();
+  _triangles.clear();
+
+  _tmp_edges.clear();
+}
+
+void
+TriangulatedSurface::complete_surface()
+{
+  complete_ref_arrays();
+  complete_int_arrays();
+
+  _completed_surface = 1;
+}
+
+void
+TriangulatedSurface::complete_ref_arrays()
+{
+  _tmp_edges.clear();
+  _index_to_edge.clear();
+  _edge_to_index.clear();
+
+  int i;
+  int ntri = ntriangle();
+  _edges.clear();
+  for (i=0; i<ntri; i++) {
+      Ref<Triangle> tri = triangle(i);
+      add_edge(tri->edge(0));
+      add_edge(tri->edge(1));
+      add_edge(tri->edge(2));
+    }
+  int ne = nedge();
+  _vertices.clear();
+  _index_to_vertex.clear();
+  _vertex_to_index.clear();
+  for (i=0; i<ne; i++) {
+      Ref<Edge> e = edge(i);
+      add_vertex(e->vertex(0));
+      add_vertex(e->vertex(1));
+    }
+}
+
+void
+TriangulatedSurface::complete_int_arrays()
+{
+  clear_int_arrays();
+  
+  int i;
+  int ntri = ntriangle();
+  int ne = nedge();
+
+  // construct the array that converts the triangle number and vertex
+  // number within the triangle to the overall vertex number
+  _triangle_vertex = new int*[ntri];
+  for (i=0; i<ntri; i++) {
+      _triangle_vertex[i] = new int[3];
+      for (int j=0; j<3; j++) {
+          Ref<Vertex> v = triangle(i)->vertex(j);
+          _triangle_vertex[i][j] = _vertex_to_index[v];
+        }
+    }
+
+  // construct the array that converts the triangle number and edge number
+  // within the triangle to the overall edge number
+  _triangle_edge = new int*[ntri];
+  for (i=0; i<ntri; i++) {
+      _triangle_edge[i] = new int[3];
+      for (int j=0; j<3; j++) {
+          Ref<Edge> e = triangle(i)->edge(j);
+          _triangle_edge[i][j] = _edge_to_index[e];
+        }
+    }
+
+  // construct the array that converts the edge number and vertex number
+  // within the edge to the overall vertex number
+  _edge_vertex = new int*[ne];
+  for (i=0; i<ne; i++) {
+      _edge_vertex[i] = new int[2];
+      for (int j=0; j<2; j++) {
+          Ref<Vertex> v = edge(i)->vertex(j);
+          _edge_vertex[i][j] = _vertex_to_index[v];
+        }
+    }
+}
+
+void
+TriangulatedSurface::compute_values(Ref<Volume>&vol)
+{
+  int n = _vertices.size();
+  _values.resize(n);
+
+  for (int i=0; i<n; i++) {
+      vol->set_x(vertex(i)->point());
+      _values[i] = vol->value();
+    }
+  _have_values = 1;
+}
+
+double
+TriangulatedSurface::flat_area()
+{
+  double result = 0.0;
+  for (std::set<Ref<Triangle> >::iterator i=_triangles.begin();
+       i!=_triangles.end(); i++) {
+      result += (*i)->flat_area();
+    }
+  return result;
+}
+
+double
+TriangulatedSurface::flat_volume()
+{
+  double result = 0.0;
+  for (int i=0; i<_triangles.size(); i++) {
+
+      // get the vertices of the triangle
+      SCVector3 A(vertex(triangle_vertex(i,0))->point());
+      SCVector3 B(vertex(triangle_vertex(i,1))->point());
+      SCVector3 C(vertex(triangle_vertex(i,2))->point());
+
+      // project the vertices onto the xy plane
+      SCVector3 Axy(A); Axy[2] = 0.0;
+      SCVector3 Bxy(B); Bxy[2] = 0.0;
+      SCVector3 Cxy(C); Cxy[2] = 0.0;
+
+      // construct the legs of the triangle in the xy plane
+      SCVector3 BAxy = Bxy - Axy;
+      SCVector3 CAxy = Cxy - Axy;
+
+      // find the lengths of the legs of the triangle in the xy plane
+      double baxy = sqrt(BAxy.dot(BAxy));
+      double caxy = sqrt(CAxy.dot(CAxy));
+
+      // if one of the legs is of length zero, then there is
+      // no contribution from this triangle
+      if (baxy < 1.e-16 || caxy < 1.e-16) continue;
+
+      // find the sine of the angle between the legs of the triangle
+      // in the xy plane
+      double costheta = BAxy.dot(CAxy)/(baxy*caxy);
+      double sintheta = sqrt(1.0 - costheta*costheta);
+
+      // the area of the triangle in the xy plane
+      double areaxy = 0.5 * baxy * caxy * sintheta;
+
+      // the height of the three corners of the triangle
+      // (relative to the z plane)
+      double hA = A[2];
+      double hB = B[2];
+      double hC = C[2];
+
+      // the volume of the space under the triangle
+      double volume = areaxy * (hA + (hB + hC - 2.0*hA)/3.0);
+
+      // the orientation of the triangle along the projection axis (z)
+      SCVector3 BA(B-A);
+      SCVector3 CA(C-A);
+      double z_orientation = BA.cross(CA)[2];
+
+      if (z_orientation > 0.0) {
+          result += volume;
+        }
+      else {
+          result -= volume;
+        }
+
+    }
+
+  // If the volume is negative, then the surface gradients were
+  // opposite in sign to the direction assumed.  Flip the sign
+  // to fix.
+  return fabs(result);
+}
+
+double
+TriangulatedSurface::area()
+{
+  double area = 0.0;
+  TriangulatedSurfaceIntegrator triint(this);
+  for (triint = 0; triint.update(); triint++) {
+      area += triint.w();
+    }
+  return area;
+}
+
+double
+TriangulatedSurface::volume()
+{
+  double volume = 0.0;
+  TriangulatedSurfaceIntegrator triint(this);
+  for (triint = 0; triint.update(); triint++) {
+      volume += triint.weight()*triint.dA()[2]*triint.current()->point()[2];
+    }
+  return volume;
+}
+
+void
+TriangulatedSurface::add_vertex(const Ref<Vertex>&t)
+{
+  int i = _vertices.size();
+  _vertices.insert(t);
+  if (i != _vertices.size()) {
+      _index_to_vertex.push_back(t);
+      _vertex_to_index[t] = i;
+      if (_index_to_vertex.size() != _vertex_to_index.size()) {
+          ExEnv::errn() << "TriangulatedSurface::add_vertex: length mismatch" << endl;
+          abort();
+        }
+    }
+}
+
+void
+TriangulatedSurface::add_edge(const Ref<Edge>&t)
+{
+  int i = _edges.size();
+  _edges.insert(t);
+  if (i != _edges.size()) {
+      _index_to_edge.push_back(t);
+      _edge_to_index[t] = i;
+      if (_index_to_edge.size() != _edge_to_index.size()) {
+          ExEnv::errn() << "TriangulatedSurface::add_edge: length mismatch" << endl;
+          abort();
+        }
+    }
+}
+
+void
+TriangulatedSurface::add_triangle(const Ref<Triangle>&t)
+{
+  if (_completed_surface) clear();
+  int i = _triangles.size();
+  _triangles.insert(t);
+  if (i != _triangles.size()) {
+      _index_to_triangle.push_back(t);
+      _triangle_to_index[t] = i;
+      if (_index_to_triangle.size() != _triangle_to_index.size()) {
+          ExEnv::errn() << "TriangulatedSurface::add_triangle: length mismatch" << endl;
+          abort();
+        }
+    }
+}
+
+void
+TriangulatedSurface::add_triangle(const Ref<Vertex>& v1,
+                                  const Ref<Vertex>& v2,
+                                  const Ref<Vertex>& v3)
+{
+  // Find this triangle's edges if they have already be created
+  // for some other triangle.
+  Ref<Edge> e0, e1, e2;
+
+  const std::set<Ref<Edge> > &v1edges = _tmp_edges[v1];
+
+  const std::set<Ref<Edge> > &v2edges = _tmp_edges[v2];
+
+  std::set<Ref<Edge> >::const_iterator ix;
+  for (ix = v1edges.begin(); ix != v1edges.end(); ix++) {
+      const Ref<Edge>& e = *ix;
+      if (e->vertex(0) == v2 || e->vertex(1) == v2) {
+          e0 = e;
+        }
+      else if (e->vertex(0) == v3 || e->vertex(1) == v3) {
+          e2 = e;
+        }
+    }
+  for (ix = v2edges.begin(); ix != v2edges.end(); ix++) {
+      const Ref<Edge>& e = *ix;
+      if (e->vertex(0) == v3 || e->vertex(1) == v3) {
+          e1 = e;
+        }
+    }
+
+  if (e0.null()) {
+      e0 = newEdge(v1,v2);
+      _tmp_edges[v1].insert(e0);
+      _tmp_edges[v2].insert(e0);
+    }
+  if (e1.null()) {
+      e1 = newEdge(v2,v3);
+      _tmp_edges[v2].insert(e1);
+      _tmp_edges[v3].insert(e1);
+    }
+  if (e2.null()) {
+      e2 = newEdge(v3,v1);
+      _tmp_edges[v3].insert(e2);
+      _tmp_edges[v1].insert(e2);
+    }
+  
+  int orientation;
+  if (e0->vertex(0) == v1) {
+      orientation = 0;
+    }
+  else {
+      orientation = 1;
+    }
+  
+  add_triangle(newTriangle(e0,e1,e2,orientation));
+}
+
+// If a user isn't keeping track of edges while add_triangle is being
+// used to build the surface, then this can be called to see if an edge
+// already exists (at a great performance cost).
+Ref<Edge>
+TriangulatedSurface::find_edge(const Ref<Vertex>& v1, const Ref<Vertex>& v2)
+{
+  std::set<Ref<Triangle> >::iterator i;
+
+  for (i=_triangles.begin(); i!=_triangles.end(); i++) {
+      Ref<Triangle> t = *i;
+      Ref<Edge> e1 = t->edge(0);
+      Ref<Edge> e2 = t->edge(1);
+      Ref<Edge> e3 = t->edge(2);
+      if (e1->vertex(0) == v1 && e1->vertex(1) == v2) return e1;
+      if (e1->vertex(1) == v1 && e1->vertex(0) == v2) return e1;
+      if (e2->vertex(0) == v1 && e2->vertex(1) == v2) return e2;
+      if (e2->vertex(1) == v1 && e2->vertex(0) == v2) return e2;
+      if (e3->vertex(0) == v1 && e3->vertex(1) == v2) return e3;
+      if (e3->vertex(1) == v1 && e3->vertex(0) == v2) return e3;
+    }
+
+  return 0;
+}
+
+void
+TriangulatedSurface::print(ostream&o) const
+{
+  o << indent << "TriangulatedSurface:" << endl;
+  int i;
+
+  int np = nvertex();
+  o << indent << scprintf(" %3d Vertices:",np) << endl;
+  for (i=0; i<np; i++) {
+      Ref<Vertex> p = vertex(i);
+      o << indent << scprintf("  %3d:",i);
+      for (int j=0; j<3; j++) {
+          o << scprintf(" % 15.10f", p->point()[j]);
+        }
+      o << endl;
+    }
+
+  int ne = nedge();
+  o << indent << scprintf(" %3d Edges:",ne) << endl;
+  for (i=0; i<ne; i++) {
+      Ref<Edge> e = edge(i);
+      Ref<Vertex> v0 = e->vertex(0);
+      Ref<Vertex> v1 = e->vertex(1);
+      std::map<Ref<Vertex>,int>::const_iterator v0i=_vertex_to_index.find(v0);
+      std::map<Ref<Vertex>,int>::const_iterator v1i=_vertex_to_index.find(v1);
+      int v0int = v0i==_vertex_to_index.end()? -1: v0i->second;
+      int v1int = v1i==_vertex_to_index.end()? -1: v1i->second;
+      o << indent
+        << scprintf("  %3d: %3d %3d",i, v0int, v1int)
+        << endl;
+    }
+
+  int nt = ntriangle();
+  o << indent << scprintf(" %3d Triangles:",nt) << endl;
+  for (i=0; i<nt; i++) {
+      Ref<Triangle> tri = triangle(i);
+      Ref<Edge> e0 = tri->edge(0);
+      Ref<Edge> e1 = tri->edge(1);
+      Ref<Edge> e2 = tri->edge(2);
+      std::map<Ref<Edge>,int>::const_iterator e0i = _edge_to_index.find(e0);
+      std::map<Ref<Edge>,int>::const_iterator e1i = _edge_to_index.find(e1);
+      std::map<Ref<Edge>,int>::const_iterator e2i = _edge_to_index.find(e2);
+      int e0int = e0i==_edge_to_index.end()? -1: e0i->second;
+      int e1int = e1i==_edge_to_index.end()? -1: e1i->second;
+      int e2int = e2i==_edge_to_index.end()? -1: e2i->second;
+      o << indent
+        << scprintf("  %3d: %3d %3d %3d",i, e0int, e1int, e2int)
+        << endl;
+    }
+}
+
+void
+TriangulatedSurface::print_vertices_and_triangles(ostream&o) const
+{
+  o << indent << "TriangulatedSurface:" << endl;
+  int i;
+
+  int np = nvertex();
+  o << indent << scprintf(" %3d Vertices:",np) << endl;
+  for (i=0; i<np; i++) {
+      Ref<Vertex> p = vertex(i);
+      o << indent << scprintf("  %3d:",i);
+      for (int j=0; j<3; j++) {
+          o << scprintf(" % 15.10f", p->point()[j]);
+        }
+      o << endl;
+    }
+
+  int nt = ntriangle();
+  o << indent << scprintf(" %3d Triangles:",nt) << endl;
+  for (i=0; i<nt; i++) {
+      Ref<Triangle> tri = triangle(i);
+      o << indent
+        << scprintf("  %3d: %3d %3d %3d",i,
+                    _triangle_vertex[i][0],
+                    _triangle_vertex[i][1],
+                    _triangle_vertex[i][2])
+        << endl;
+    }
+}
+
+void
+TriangulatedSurface::render(const Ref<Render> &render)
+{
+  Ref<RenderedPolygons> poly = new RenderedPolygons;
+  poly->initialize(_vertices.size(), _triangles.size(),
+                   RenderedPolygons::Vertex);
+  std::set<Ref<Vertex> >::iterator iv;
+  std::set<Ref<Triangle> >::iterator it;
+  std::map<Ref<Vertex>, int> vertex_to_index;
+  int i = 0;
+  for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
+      Ref<Vertex> v = *iv;
+      vertex_to_index[v] = i;
+      poly->set_vertex(i,
+                       v->point()[0],
+                       v->point()[1],
+                       v->point()[2]);
+      poly->set_vertex_rgb(i, 0.3, 0.3, 0.3);
+    }
+  i = 0;
+  for (it = _triangles.begin(); it != _triangles.end(); it++, i++) {
+      Ref<Triangle> t = *it;
+      poly->set_face(i,
+                     vertex_to_index[t->vertex(0)],
+                     vertex_to_index[t->vertex(1)],
+                     vertex_to_index[t->vertex(2)]);
+    }
+  render->render(poly.pointer());
+}
+
+void
+TriangulatedSurface::print_geomview_format(ostream&o) const
+{
+  o << "OFF" << endl;
+
+  o << nvertex() << " " << ntriangle() << " " << nedge() << endl;
+  int i;
+
+  int np = nvertex();
+  for (i=0; i<np; i++) {
+      Ref<Vertex> p = vertex(i);
+      for (int j=0; j<3; j++) {
+          o << scprintf(" % 15.10f", p->point()[j]);
+        }
+      o << endl;
+    }
+
+  int nt = ntriangle();
+  for (i=0; i<nt; i++) {
+      Ref<Triangle> tri = triangle(i);
+      o << scprintf(" 3 %3d %3d %3d",
+                    _triangle_vertex[i][0],
+                    _triangle_vertex[i][1],
+                    _triangle_vertex[i][2])
+          << endl;
+    }
+}
+
+void
+TriangulatedSurface::recompute_index_maps()
+{
+  int i;
+  std::set<Ref<Vertex> >::iterator iv;
+  std::set<Ref<Edge> >::iterator ie;
+  std::set<Ref<Triangle> >::iterator it;
+
+  // fix the index maps
+  _vertex_to_index.clear();
+  _edge_to_index.clear();
+  _triangle_to_index.clear();
+
+  _index_to_vertex.clear();
+  _index_to_edge.clear();
+  _index_to_triangle.clear();
+
+  _index_to_vertex.resize(_vertices.size());
+  for (i=0, iv = _vertices.begin(); iv != _vertices.end(); i++, iv++) {
+      _vertex_to_index[*iv] = i;
+      _index_to_vertex[i] = *iv;
+    }
+
+  _index_to_edge.resize(_edges.size());
+  for (i=0, ie = _edges.begin(); ie != _edges.end(); i++, ie++) {
+      _edge_to_index[*ie] = i;
+      _index_to_edge[i] = *ie;
+    }
+
+  _index_to_triangle.resize(_triangles.size());
+  for (i=0, it = _triangles.begin(); it != _triangles.end(); i++, it++) {
+      _triangle_to_index[*it] = i;
+      _index_to_triangle[i] = *it;
+    }
+
+}
+
+Edge*
+TriangulatedSurface::newEdge(const Ref<Vertex>& v0, const Ref<Vertex>& v1) const
+{
+  return new Edge(v0,v1);
+}
+
+Triangle*
+TriangulatedSurface::newTriangle(const Ref<Edge>& e0,
+                                 const Ref<Edge>& e1,
+                                 const Ref<Edge>& e2,
+                                 int orientation) const
+{
+  return new Triangle(e0,e1,e2,orientation);
+}
+
+//////////////////////////////////////////////////////////////////////
+// TriangulatedSurfaceIntegrator
+
+TriangulatedSurfaceIntegrator::
+  TriangulatedSurfaceIntegrator(const Ref<TriangulatedSurface>&ts)
+{
+  set_surface(ts);
+  use_default_integrator();
+
+  _itri = 0;
+  _irs = 0;
+}
+
+TriangulatedSurfaceIntegrator::
+  TriangulatedSurfaceIntegrator()
+{
+  use_default_integrator();
+}
+
+void
+TriangulatedSurfaceIntegrator::
+  operator =(const TriangulatedSurfaceIntegrator&i)
+{
+  set_surface(i._ts);
+
+  _integrator = i._integrator;
+  _itri = i._itri;
+  _irs = i._irs;
+}
+
+TriangulatedSurfaceIntegrator::
+  ~TriangulatedSurfaceIntegrator()
+{
+}
+
+void
+TriangulatedSurfaceIntegrator::set_surface(const Ref<TriangulatedSurface>&s)
+{
+  _ts = s;
+  _current = new Vertex();
+}
+
+int
+TriangulatedSurfaceIntegrator::
+  vertex_number(int i)
+{
+  return _ts->triangle_vertex(_itri,i);
+}
+
+Ref<Vertex>
+TriangulatedSurfaceIntegrator::
+  current()
+{
+  return _current;
+}
+
+int
+TriangulatedSurfaceIntegrator::n()
+{
+  int result = 0;
+  int ntri = _ts->ntriangle();
+  for (int i=0; i<ntri; i++) {
+      result += (_ts.pointer()->*_integrator)(i)->n();
+    }
+  return result;
+}
+
+int
+TriangulatedSurfaceIntegrator::update()
+{
+  if (_itri < 0 || _itri >= _ts->ntriangle()) return 0;
+
+  TriangleIntegrator* i = (_ts.pointer()->*_integrator)(_itri).pointer();
+  _s = i->s(_irs);
+  _r = i->r(_irs);
+  _weight = i->w(_irs);
+  Ref<Triangle> t = _ts->triangle(_itri);
+  Ref<TriInterpCoef> coef = i->coef(t->order(),_irs);
+  t->interpolate(coef, _r, _s, _current, _dA);
+  _surface_element = _dA.norm();
+  static double cum;
+  if (_irs == 0) cum = 0.0;
+  cum += _surface_element * _weight;
+  //ExEnv::outn() << scprintf("%2d dA = %12.8f, w = %12.8f, Sum wdA = %12.8f",
+  //                 _irs, _surface_element, _weight, cum)
+  //     << endl;
+
+  return (int) 1;
+}
+
+void
+TriangulatedSurfaceIntegrator::
+  operator ++()
+{
+  int n = (_ts.pointer()->*_integrator)(_itri)->n();
+  if (_irs == n-1) {
+      _irs = 0;
+      if (_grp.null()) _itri++;
+      else _itri += _grp->n();
+    }
+  else {
+      _irs++;
+    }
+}
+
+void
+TriangulatedSurfaceIntegrator::distribute(const Ref<MessageGrp> &grp)
+{
+  _grp = grp;
+}
+
+int
+TriangulatedSurfaceIntegrator::
+  operator = (int i)
+{
+  _itri = i;
+  _irs = 0;
+  return i;
+}
+
+void
+TriangulatedSurfaceIntegrator::use_default_integrator()
+{
+  _integrator = &TriangulatedSurface::integrator;
+}
+
+void
+TriangulatedSurfaceIntegrator::use_fast_integrator()
+{
+  _integrator = &TriangulatedSurface::fast_integrator;
+}
+
+void
+TriangulatedSurfaceIntegrator::use_accurate_integrator()
+{
+  _integrator = &TriangulatedSurface::accurate_integrator;
+}
+
+/////////////////////////////////////////////////////////////////////////
+// TriangulatedImplicitSurface
+static ClassDesc TriangulatedImplicitSurface_cd(
+  typeid(TriangulatedImplicitSurface),"TriangulatedImplicitSurface",1,"public TriangulatedSurface",
+  0, create<TriangulatedImplicitSurface>, 0);
+
+TriangulatedImplicitSurface::
+TriangulatedImplicitSurface(const Ref<KeyVal>&keyval):
+  TriangulatedSurface(keyval)
+{
+  inited_ = 0;
+
+  vol_ << keyval->describedclassvalue("volume");
+  if (keyval->error() != KeyVal::OK) {
+      ExEnv::errn() << "TriangulatedImplicitSurface(const Ref<KeyVal>&keyval): "
+           << "requires \"volume\"" << endl;
+      abort();
+    }
+
+  isovalue_ = keyval->doublevalue("value");
+  if (keyval->error() != KeyVal::OK) isovalue_ = 0.0;
+
+  fix_orientation_ = keyval->booleanvalue("fix_orientation");
+  if (keyval->error() != KeyVal::OK) fix_orientation_ = 1;
+
+  remove_short_edges_ = keyval->booleanvalue("remove_short_edges");
+  if (keyval->error() != KeyVal::OK) remove_short_edges_ = 1;
+
+  remove_slender_triangles_ = keyval->booleanvalue("remove_slender_triangles");
+  if (keyval->error() != KeyVal::OK) remove_slender_triangles_ = 0;
+
+  remove_small_triangles_ = keyval->booleanvalue("remove_small_triangles");
+  if (keyval->error() != KeyVal::OK) remove_small_triangles_ = 0;
+
+  short_edge_factor_ = keyval->doublevalue("short_edge_factor");
+  if (keyval->error() != KeyVal::OK) short_edge_factor_ = 0.4;
+
+  slender_triangle_factor_ = keyval->doublevalue("slender_triangle_factor");
+  if (keyval->error() != KeyVal::OK) slender_triangle_factor_ = 0.2;
+
+  small_triangle_factor_ = keyval->doublevalue("small_triangle_factor");
+  if (keyval->error() != KeyVal::OK) small_triangle_factor_ = 0.2;
+
+  resolution_ = keyval->doublevalue("resolution");
+  if (keyval->error() != KeyVal::OK) resolution_ = 1.0;
+
+  order_ = keyval->intvalue("order");
+  if (keyval->error() != KeyVal::OK) order_ = 1;
+
+  int initialize = keyval->booleanvalue("initialize");
+  if (initialize) init();
+}
+
+void
+TriangulatedImplicitSurface::init()
+{
+  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: init start" << endl;
+  ImplicitSurfacePolygonizer isogen(vol_);
+  isogen.set_resolution(resolution_);
+
+  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: isosurface" << endl;
+  isogen.isosurface(isovalue_,*this);
+#if WRITE_OOGL
+  if (_debug) {
+      render(new OOGLRender("surfiso.oogl"));
+    }
+#endif
+  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
+  if (fix_orientation_) fix_orientation();
+#if WRITE_OOGL
+  if (_debug) {
+      render(new OOGLRender("surffix.oogl"));
+    }
+#endif
+  if (remove_short_edges_) {
+      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: short edges" << endl;
+      remove_short_edges(short_edge_factor_*resolution_,vol_,isovalue_);
+      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
+      if (fix_orientation_) fix_orientation();
+    }
+  if (remove_slender_triangles_ || remove_small_triangles_) {
+      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: slender" << endl;
+      double height_cutoff = slender_triangle_factor_ * resolution_;
+      double area_cutoff = small_triangle_factor_*resolution_*resolution_*0.5;
+      remove_slender_triangles(remove_slender_triangles_, height_cutoff,
+                               remove_small_triangles_, area_cutoff,
+                               vol_,isovalue_);
+      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: orientation" << endl;
+      if (fix_orientation_) fix_orientation();
+    }
+
+  // see if a higher order approximation to the surface is required
+  if (order_ > 1) {
+      if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: order" << endl;
+      int i;
+      for (i=0; i<nedge(); i++) {
+          edge(i)->set_order(order_, vol_, isovalue_);
+        }
+      for (i=0; i<ntriangle(); i++) {
+          triangle(i)->set_order(order_, vol_, isovalue_);
+        }
+    }
+  inited_ = 1;
+  if (_verbose) ExEnv::outn() << "TriangulatedImplicitSurface: init done" << endl;
+}
+
+TriangulatedImplicitSurface::~TriangulatedImplicitSurface()
+{
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/surf.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/surf.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/surf.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,302 @@
+//
+// surf.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_surf_h
+#define _math_isosurf_surf_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+
+#include <map>
+#include <set>
+#include <vector>
+
+#include <math/isosurf/triangle.h>
+#include <math/isosurf/volume.h>
+#include <util/render/render.h>
+
+namespace sc {
+
+template <class C, class I>
+inline void
+erase_elements_by_value(C &container, I begin, I end)
+{
+  for (I i=begin; i!=end; i++) {
+      container.erase(*i);
+    }
+}
+
+class TriangulatedSurface: public DescribedClass {
+  protected:
+    int _verbose;
+    int _debug;
+
+    int _completed_surface;
+
+    // sets of objects that make up the surface
+    std::set<Ref<Vertex> > _vertices;
+    std::set<Ref<Edge> > _edges;
+    std::set<Ref<Triangle> > _triangles;
+
+    // map objects to an integer index
+    std::map<Ref<Vertex>,int> _vertex_to_index;
+    std::map<Ref<Edge>,int> _edge_to_index;
+    std::map<Ref<Triangle>,int> _triangle_to_index;
+
+    // map integer indices to an object
+    std::vector<Ref<Vertex> > _index_to_vertex;
+    std::vector<Ref<Edge> > _index_to_edge;
+    std::vector<Ref<Triangle> > _index_to_triangle;
+
+    // mappings between array element numbers
+    int** _triangle_vertex;
+    int** _triangle_edge;
+    int** _edge_vertex;
+
+    // values for each of the vertices
+    int _have_values;
+    std::vector<double> _values;
+
+    // what to use to integrate over the surface, by default
+    Ref<TriangleIntegrator> _integrator;
+    // other integrators, in terms of time & accuracy:
+    // _fast_integrator <= _integrator <= _accurate_interator
+    Ref<TriangleIntegrator> _fast_integrator;
+    Ref<TriangleIntegrator> _accurate_integrator;
+
+    void clear_int_arrays();
+
+    void complete_ref_arrays();
+    void complete_int_arrays();
+
+    void recompute_index_maps();
+
+    void add_triangle(const Ref<Triangle>&);
+    void add_vertex(const Ref<Vertex>&);
+    void add_edge(const Ref<Edge>&);
+
+    // these members must be used to allocate new triangles and edges
+    // since specializations of TriangulatedSurface might need to
+    // override these to produce triangles and edges with interpolation
+    // data.
+    virtual Triangle* newTriangle(const Ref<Edge>&,
+                                  const Ref<Edge>&,
+                                  const Ref<Edge>&,
+                                  int orientation) const;
+    virtual Edge* newEdge(const Ref<Vertex>&,const Ref<Vertex>&) const;
+
+    // this map of edges to vertices is used to construct the surface
+    std::map<Ref<Vertex>,std::set<Ref<Edge> > > _tmp_edges;
+  public:
+    TriangulatedSurface();
+    TriangulatedSurface(const Ref<KeyVal>&);
+    virtual ~TriangulatedSurface();
+
+    // control printing
+    int verbose() const { return _verbose; }
+    void verbose(int v) { _verbose = v; }
+
+    // set up an integrator
+    void set_integrator(const Ref<TriangleIntegrator>&);
+    void set_fast_integrator(const Ref<TriangleIntegrator>&);
+    void set_accurate_integrator(const Ref<TriangleIntegrator>&);
+    virtual Ref<TriangleIntegrator> integrator(int itri);
+    virtual Ref<TriangleIntegrator> fast_integrator(int itri);
+    virtual Ref<TriangleIntegrator> accurate_integrator(int itri);
+
+    // construct the surface
+    void add_triangle(const Ref<Vertex>&,
+                      const Ref<Vertex>&,
+                      const Ref<Vertex>&);
+    Ref<Edge> find_edge(const Ref<Vertex>&, const Ref<Vertex>&);
+    virtual void complete_surface();
+
+    // clean up the surface
+    virtual void remove_short_edges(double cutoff_length = 1.0e-6,
+                                    const Ref<Volume> &vol=0, double isoval=0.0);
+    virtual void remove_slender_triangles(
+                                    int remove_slender, double height_cutoff,
+                                    int remove_small, double area_cutoff,
+                                    const Ref<Volume> &vol=0, double isoval=0.0);
+    virtual void fix_orientation();
+    virtual void clear();
+
+    // get information from the object sets
+    int nvertex() const { return _vertices.size(); };
+    Ref<Vertex> vertex(int i) const { return _index_to_vertex[i]; };
+    int vertex_index(const Ref<Vertex> &o) {
+      std::map<Ref<Vertex>,int>::iterator i = _vertex_to_index.find(o);
+      if (i != _vertex_to_index.end()) return i->second;
+      return -1;
+    }
+    int nedge() const { return _edges.size(); };
+    Ref<Edge> edge(int i) const { return _index_to_edge[i]; };
+    int edge_index(const Ref<Edge> &o) {
+      std::map<Ref<Edge>,int>::iterator i = _edge_to_index.find(o);
+      if (i != _edge_to_index.end()) return i->second;
+      return -1;
+    }
+    int ntriangle() const { return _triangles.size(); };
+    Ref<Triangle> triangle(int i) const { return _index_to_triangle[i]; }
+    int triangle_index(const Ref<Triangle> &o) {
+      std::map<Ref<Triangle>,int>::iterator i = _triangle_to_index.find(o);
+      if (i != _triangle_to_index.end()) return i->second;
+      return -1;
+    }
+
+    // information from the index mappings
+    int triangle_vertex(int i,int j) const { return _triangle_vertex[i][j]; };
+    int triangle_edge(int i,int j) const { return _triangle_edge[i][j]; };
+    int edge_vertex(int i,int j) const { return _edge_vertex[i][j]; };
+
+    // associate values with vertices
+    //void compute_colors(Volume&);
+    void compute_values(Ref<Volume>&);
+
+    // properties of the surface
+    virtual double flat_area(); // use flat triangles
+    virtual double flat_volume(); // use flat triangles
+    virtual double area();
+    virtual double volume();
+
+    // output of the surface
+    virtual void print(std::ostream&o=ExEnv::out0()) const;
+    virtual void print_vertices_and_triangles(std::ostream&o=ExEnv::out0()) const;
+    virtual void print_geomview_format(std::ostream&o=ExEnv::out0()) const;
+    virtual void render(const Ref<Render> &render);
+
+    // print information about the topology
+    void topology_info(std::ostream&o=ExEnv::out0());
+    void topology_info(int nvertex, int nedge, int ntri, std::ostream&o=ExEnv::out0());
+};
+
+
+class TriangulatedSurfaceIntegrator {
+  private:
+    Ref<TriangulatedSurface> _ts;
+    int _itri;
+    int _irs;
+    double _r;
+    double _s;
+    double _weight;
+    double _surface_element;
+    Ref<Vertex> _current;
+    SCVector3 _dA;
+    Ref<TriangleIntegrator> (TriangulatedSurface::*_integrator)(int itri);
+    Ref<MessageGrp> _grp;
+  public:
+    TriangulatedSurfaceIntegrator();
+    // the surface cannot be changed until this is destroyed
+    TriangulatedSurfaceIntegrator(const Ref<TriangulatedSurface>&);
+    ~TriangulatedSurfaceIntegrator();
+    // Objects initialized by these operators are not automatically
+    // updated.  This must be done with the update member.
+    // The _grp is not copied.
+    void operator = (const TriangulatedSurfaceIntegrator&);
+    TriangulatedSurfaceIntegrator(const TriangulatedSurfaceIntegrator&i) {
+        operator = (i);
+      }
+    // Return the number of integration points.
+    int n();
+    // Assign the surface.  Don't do this while iterating.
+    void set_surface(const Ref<TriangulatedSurface>&);
+    // returns the number of the vertex in the current triangle
+    int vertex_number(int i);
+    inline double r() const { return _r; }
+    inline double s() const { return _s; }
+    inline double w() const { return _weight*_surface_element; }
+    double surface_element() const { return _surface_element; }
+    double weight() const { return _weight; }
+    const SCVector3& dA() const { return _dA; }
+    Ref<Vertex> current();
+    // Tests to see if this point is valid, if it is then
+    // _r, _s, etc are computed and 1 is returned.
+    int update();
+    // This can be used to loop through unique pairs of points.
+    // The argument should be a TriangulatedSurfaceIntegrator for
+    // the same surface as this.
+    int operator < (TriangulatedSurfaceIntegrator&i) {
+        update();
+        return _itri<i._itri?1:(_itri>i._itri?0:(_irs<i._irs?1:0));
+      }
+    // Goes to the next point.  Does not update.
+    void operator++();
+    inline void operator++(int) { operator++(); }
+    // setting TSI = i sets TSI to begin at the triangle i
+    int operator = (int);
+    int itri() const { return _itri; }
+    int irs() const { return _irs; }
+    // the number of points in the current triangle
+    int n_in_tri() const { return (_ts.pointer()->*_integrator)(_itri)->n(); }
+    void distribute(const Ref<MessageGrp> &);
+    void use_fast_integrator();
+    void use_accurate_integrator();
+    void use_default_integrator();
+};
+
+class TriangulatedImplicitSurface: public TriangulatedSurface {
+  private:
+    // The surface is defined as an isosurface of the volume vol_.
+    Ref<Volume> vol_;
+    double isovalue_;
+
+    int fix_orientation_;
+    int remove_short_edges_;
+    double short_edge_factor_;
+    int remove_slender_triangles_;
+    double slender_triangle_factor_;
+    int remove_small_triangles_;
+    double small_triangle_factor_;
+    double resolution_;
+
+    int order_;
+
+    int inited_;
+  public:
+    TriangulatedImplicitSurface(const Ref<KeyVal>&);
+    ~TriangulatedImplicitSurface();
+
+    Ref<Volume> volume_object() const { return vol_; }
+    double isovalue() const { return isovalue_; }
+
+    void init();
+    int inited() const { return inited_; }
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/surfor.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/surfor.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/surfor.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,146 @@
+//
+// surfor.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <util/misc/formio.h>
+
+#include <math/scmat/matrix.h>
+#include <math/isosurf/surf.h>
+
+using namespace std;
+using namespace sc;
+
+void
+TriangulatedSurface::fix_orientation()
+{
+  int i,j;
+  std::set<Ref<Triangle> >::iterator I;
+  int nflip = 0;
+  
+  int ne = nedge();
+  int ntri = ntriangle();
+  Ref<Triangle> *edge_to_triangle0;
+  Ref<Triangle> *edge_to_triangle1;
+  edge_to_triangle0 = new Ref<Triangle>[ne];
+  edge_to_triangle1 = new Ref<Triangle>[ne];
+
+  for (I = _triangles.begin(); I != _triangles.end(); I++) {
+      Ref<Triangle> tri = *I;
+      for (j=0; j<3; j++) {
+          Ref<Edge> e = tri->edge(j);
+          int e_index = _edge_to_index[e];
+          if (edge_to_triangle0[e_index].null()) {
+              edge_to_triangle0[e_index] = tri;
+            }
+          else if (edge_to_triangle1[e_index].null()) {
+              edge_to_triangle1[e_index] = tri;
+            }
+          else {
+              ExEnv::errn() << "TriangulatedSurface::fix_orientation:"
+                   << " more than two triangles to an edge" << endl;
+              abort();
+            }
+        }
+    }
+
+  std::set<Ref<Triangle> > unfixed;
+  std::set<Ref<Triangle> > fixed;
+  std::set<Ref<Triangle> > finished;
+
+  unfixed.insert(_triangles.begin(), _triangles.end());
+
+  while (unfixed.size()) {
+      // define unfixed.first()'s orientation to be the fixed orientation
+      std::set<Ref<Triangle> >::iterator first = unfixed.begin();
+      fixed.insert(*first);
+      unfixed.erase(*first);
+      while (fixed.size()) {
+          Ref<Triangle> tri = *fixed.begin();
+          // make all neighbors of tri oriented the same as tri
+          for (i=0; i<3; i++) {
+              Ref<Edge> e = tri->edge(i);
+              int e_index = _edge_to_index[e];
+              Ref<Triangle> othertri;
+              if (edge_to_triangle0[e_index] == tri) {
+                  othertri = edge_to_triangle1[e_index];
+                }
+              else {
+                  othertri = edge_to_triangle0[e_index];
+                }
+              for (j=0; j<3; j++) {
+                  if (othertri->edge(j) == e) break;
+                }
+              if (j == 3) {
+                  ExEnv::errn() << "TriangulatedSurface::fix_orientation: "
+                       << " edge_to_triangle wrong" << endl;
+                  abort();
+                }
+              if (tri->orientation(i) == othertri->orientation(j)) {
+                  if (unfixed.find(othertri) != unfixed.end()) {
+                      unfixed.erase(othertri);
+                      fixed.insert(othertri);
+                      othertri->flip();
+                      nflip++;
+                    }
+                  else {
+                      ExEnv::errn() << "TriangulatedSurface::fix_orientation:"
+                           << " tried to flip a fixed triangle" << endl;
+                      abort();
+                    }
+                }
+              else if (unfixed.find(othertri) != unfixed.end()) {
+                  unfixed.erase(othertri);
+                  fixed.insert(othertri);
+                }
+            }
+          // by this point all of tri's neighbors have had their
+          // orientation fixed to match that of tri
+          fixed.erase(tri);
+          finished.insert(tri);
+        }
+    }
+
+  if (_verbose) {
+      ExEnv::outn()
+          << scprintf("%d out of %d triangles were flipped\n", nflip, ntri);
+    }
+
+  // just in case
+  for (i=0; i<ne; i++) {
+      edge_to_triangle0[i] = 0;
+      edge_to_triangle1[i] = 0;
+    }
+
+  delete[] edge_to_triangle0;
+  delete[] edge_to_triangle1;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/surfse.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/surfse.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/surfse.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,320 @@
+//
+// surfse.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <algorithm>
+
+#include <util/misc/formio.h>
+
+#include <math/scmat/matrix.h>
+#include <math/isosurf/surf.h>
+
+using namespace std;
+using namespace sc;
+
+void
+TriangulatedSurface::remove_short_edges(double length_cutoff,
+                                        const Ref<Volume> &vol, double isoval)
+{
+  int j,k;
+  std::set<Ref<Triangle> >::iterator it,jt,kt;
+  std::set<Ref<Edge> >::iterator ie,je;
+
+  int surface_was_completed = _completed_surface;
+
+  if (!_completed_surface) {
+      complete_ref_arrays();
+    }
+  else {
+      clear_int_arrays();
+    }
+
+  _have_values = 0;
+  _values.clear();
+  
+  if (_verbose) {
+      ExEnv::outn() << "TriangulatedSurface::remove_short_edges:" << endl
+           << "initial: ";
+      topology_info();
+    }
+
+  int deleted_edges_length;
+  do {
+      std::set<Ref<Triangle> > deleted_triangles;
+      std::set<Ref<Edge> > deleted_edges;
+      std::set<Ref<Vertex> > deleted_vertices;
+
+      std::set<Ref<Triangle> > new_triangles;
+      std::set<Ref<Edge> > new_edges;
+      std::set<Ref<Vertex> > new_vertices;
+
+      // a vertex to set-of-connected-triangles map
+      std::map<Ref<Vertex>,std::set<Ref<Triangle> > > connected_triangle_map;
+      for (it = _triangles.begin(); it != _triangles.end(); it++) {
+          Ref<Triangle> tri = *it;
+          for (j = 0; j<3; j++) {
+              Ref<Vertex> v  = tri->vertex(j);
+              connected_triangle_map[v].insert(tri);
+            }
+        }
+
+      // a vertex to set-of-connected-edges map
+      std::map<Ref<Vertex>,std::set<Ref<Edge> > > connected_edge_map;
+      for (ie = _edges.begin(); ie != _edges.end(); ie++) {
+          Ref<Edge> e = *ie;
+          for (j = 0; j<2; j++) {
+              Ref<Vertex> v = e->vertex(j);
+              connected_edge_map[v].insert(e);
+            }
+        }
+  
+      for (ie = _edges.begin(); ie != _edges.end(); ie++) {
+          Ref<Edge> edge = *ie;
+          double length = edge->straight_length();
+          if (length < length_cutoff) {
+              std::set<Ref<Triangle> > connected_triangles;
+              Ref<Vertex> v0 = edge->vertex(0);
+              Ref<Vertex> v1 = edge->vertex(1);
+              connected_triangles.insert(connected_triangle_map[v0].begin(),
+                                         connected_triangle_map[v0].end());
+              connected_triangles.insert(connected_triangle_map[v1].begin(),
+                                         connected_triangle_map[v1].end());
+              int skip = 0;
+              for (jt = connected_triangles.begin();
+                   jt != connected_triangles.end();
+                   jt++) {
+                  Ref<Triangle> tri = *jt;
+                  if (tri->edge(0)->straight_length() < length
+                      ||tri->edge(1)->straight_length() < length
+                      ||tri->edge(2)->straight_length() < length
+                      ||deleted_triangles.find(tri)!=deleted_triangles.end()) {
+                      skip = 1;
+                      break;
+                    }
+                }
+              if (skip) continue;
+              deleted_triangles.insert(connected_triangles.begin(),
+                                       connected_triangles.end());
+	      v0 = edge->vertex(0);
+	      v1 = edge->vertex(1);
+              deleted_vertices.insert(v0);
+              deleted_vertices.insert(v1);
+              // find all of the edges connected to the short edge
+              // (including the short edge)
+              std::set<Ref<Edge> > connected_edges;
+	      v0 = edge->vertex(0);
+	      v1 = edge->vertex(1);
+              connected_edges.insert(connected_edge_map[v0].begin(),
+                                     connected_edge_map[v0].end());
+              connected_edges.insert(connected_edge_map[v1].begin(),
+                                     connected_edge_map[v1].end());
+              deleted_edges.insert(connected_edges.begin(),
+                                   connected_edges.end());
+              // find the edges forming the perimeter of the deleted triangles
+              // (these are used to form the new triangles)
+              std::set<Ref<Edge> > perimeter_edges;
+              int embedded_triangle = 0;
+              for (jt = connected_triangles.begin();
+                   jt != connected_triangles.end();
+                   jt++) {
+                  Ref<Triangle> tri = *jt;
+                  for (j=0; j<3; j++) {
+                      Ref<Edge> e = tri->edge(j);
+                      if (connected_edges.find(e) == connected_edges.end()) {
+                          // check to see if another triangle has claimed
+                          // that this edge is a perimeter edge.  if so
+                          // then this isn't a perimeter edge after all
+                          // and it must be deleted.  this also implies
+                          // that there is at least one triangle that
+                          // has no perimeter edge.  these triangles and
+                          // their nonperimeter vertices must be
+                          // deleted.
+                          Ref<Edge> e = tri->edge(j);
+                          if (perimeter_edges.find(e)
+                              != perimeter_edges.end()) {
+                              perimeter_edges.erase(e);
+                              deleted_edges.insert(e);
+                              embedded_triangle = 1;
+                            }
+                          else {
+                              perimeter_edges.insert(e);
+                            }
+                        }
+                    }
+                }
+              // if a triangle is embedded make sure its vertices are
+              // all deleted.  (the triangle itself should already be
+              // connected and thus deleted).
+              if (embedded_triangle) {
+                  // make a list of vertices on the perimeter (so i
+                  // don't delete them
+                  std::set<Ref<Vertex> > perimeter_vertices;
+                  for (je = perimeter_edges.begin();
+                       je != perimeter_edges.end();
+                       je++) {
+                      Ref<Edge> e = *je;
+                      for (j=0; j<2; j++) {
+                          Ref<Vertex> v = e->vertex(j);
+                          perimeter_vertices.insert(v);
+                        }
+                    }
+                  // find the embedded_triangle
+                  for (jt = connected_triangles.begin();
+                       jt != connected_triangles.end();
+                       jt++) {
+                      Ref<Triangle> tri = *jt;
+                      // see if this triangle is embedded
+                      for (j=0; j<3; j++) {
+                          Ref<Edge> e = tri->edge(j);
+                          if (perimeter_edges.find(e) != perimeter_edges.end())
+                              break;
+                        }
+                      // if embedded then delete the triangle's vertices
+                      if (j==3) {
+                          for (j=0; j<3; j++) {
+                              Ref<Vertex> v = tri->vertex(j);
+                              if (perimeter_vertices.find(v) == perimeter_vertices.end())
+                                  deleted_vertices.insert(v);
+                            }
+                        }
+                    }
+                }
+              // find a new point that replaces the deleted edge
+              // (for now use one of the original, since it must lie on the
+              // analytic surface)
+              Ref<Vertex> replacement_vertex = edge->vertex(0);
+              // however, if we have a volume, find a new vertex on
+              // the analytic surface near the center of the edge
+              if (vol.nonnull()) {
+                  SCVector3 point, norm;
+                  int hn = edge->interpolate(0.5,point,norm,vol,isoval);
+                  replacement_vertex = new Vertex(point);
+                  if (hn) replacement_vertex->set_normal(norm);
+                }
+              new_vertices.insert(replacement_vertex);
+              // for each vertex on the perimeter form a new edge to the
+              // replacement vertex
+              std::map<Ref<Vertex>,Ref<Edge> > new_edge_map;
+              for (je = perimeter_edges.begin(); je!=perimeter_edges.end();
+                   je++) {
+                  Ref<Edge> e = *je;
+                  for (k = 0; k<2; k++) {
+                      Ref<Vertex> v = e->vertex(k);
+                      if (new_edge_map.find(v) == new_edge_map.end()) {
+                          Ref<Edge> new_e = newEdge(
+                              replacement_vertex,
+                              v
+                              );
+                          new_edge_map[v] = new_e;
+                          new_edges.insert(new_e);
+                        }
+                    }
+                }
+              // for each edge on the perimeter form a new triangle with the
+              // replacement vertex
+              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
+                   je++) {
+                  Ref<Edge> e1 = *je;
+                  Ref<Vertex> v0 = e1->vertex(0);
+                  Ref<Vertex> v1 = e1->vertex(1);
+                  Ref<Edge> e2 = new_edge_map[v0];
+                  Ref<Edge> e3 = new_edge_map[v1];
+                  if (e1.null() || e2.null() || e3.null()) {
+                      ExEnv::errn() << "TriangulatedSurface::remove_short_edges: "
+                           << "building new triangle but edges are null:"
+                           << endl;
+                      if (e1.null()) ExEnv::errn() << "  e1" << endl;
+                      if (e2.null()) ExEnv::errn() << "  e2" << endl;
+                      if (e3.null()) ExEnv::errn() << "  e3" << endl;
+                      abort();
+                    }
+                  // Compute the correct orientation of e1 within the new
+                  // triangle, by finding the orientation within the old
+                  // triangle.
+                  int orientation = 0;
+                  for (kt = connected_triangles.begin();
+                       kt != connected_triangles.end();
+                       kt++) {
+                      if ((*kt)->contains(e1)) {
+                          orientation
+                              = (*kt)->orientation(e1);
+                          break;
+                        }
+                    }
+                  Ref<Triangle> newtri(newTriangle(e1,e2,e3,orientation));
+                  new_triangles.insert(newtri);
+                }
+              //ExEnv::outn() << "WARNING: only one short edge removed" << endl;
+              //break;
+            }
+        }
+
+      erase_elements_by_value(_triangles,
+                              deleted_triangles.begin(),
+                              deleted_triangles.end());
+
+      erase_elements_by_value(_edges,
+                              deleted_edges.begin(),
+                              deleted_edges.end());
+
+      erase_elements_by_value(_vertices,
+                              deleted_vertices.begin(),
+                              deleted_vertices.end());
+
+      _triangles.insert(new_triangles.begin(), new_triangles.end());
+      _edges.insert(new_edges.begin(), new_edges.end());
+      _vertices.insert(new_vertices.begin(), new_vertices.end());
+
+      if (_verbose) {
+          topology_info();
+        }
+
+      deleted_edges_length = deleted_edges.size();
+      //ExEnv::outn() << "WARNING: one pass short edge removal" << endl;
+      //deleted_edges_length = 0; // do one pass
+    } while(deleted_edges_length != 0);
+
+  // fix the index maps
+  recompute_index_maps();
+
+  // restore the int arrays if they were there to begin with
+  if (surface_was_completed) {
+      complete_int_arrays();
+      _completed_surface = 1;
+    }
+
+  if (_verbose) {
+      ExEnv::outn() << "final: ";
+      topology_info();
+    }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/surfst.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/surfst.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/surfst.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,488 @@
+//
+// surfst.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#include <stdio.h>
+#include <math.h>
+
+#include <util/misc/formio.h>
+
+#include <math/scmat/matrix.h>
+#include <math/isosurf/surf.h>
+
+using namespace std;
+using namespace sc;
+
+#ifndef WRITE_OOGL // this is useful for debugging this routine
+#define WRITE_OOGL 1
+#endif
+
+#if WRITE_OOGL
+#include <util/render/oogl.h>
+#include <util/render/polygons.h>
+#endif
+
+void
+TriangulatedSurface::remove_slender_triangles(
+    int remove_slender, double height_cutoff,
+    int remove_small, double area_cutoff,
+    const Ref<Volume> &vol, double isoval)
+{
+  int i,j,k;
+  std::set<Ref<Triangle> >::iterator it,jt,kt;
+  std::set<Ref<Edge> >::iterator ie,je,ke;
+  std::set<Ref<Vertex> >::iterator iv,jv;
+
+  int surface_was_completed = _completed_surface;
+
+  if (!_completed_surface) {
+      complete_ref_arrays();
+    }
+  else {
+      clear_int_arrays();
+    }
+
+  _have_values = 0;
+  _values.clear();
+  
+  if (_verbose) {
+      ExEnv::outn() << "TriangulatedSurface::remove_slender_triangles:" << endl
+           << "initial: ";
+      topology_info();
+    }
+
+#if WRITE_OOGL
+  if (_debug) {
+      render(new OOGLRender("surfstinit.oogl"));
+    }
+#endif
+
+  int deleted_edges_length;
+  do {
+      std::set<Ref<Triangle> > deleted_triangles;
+      std::set<Ref<Edge> > deleted_edges;
+      std::set<Ref<Vertex> > deleted_vertices;
+
+      std::set<Ref<Vertex> > vertices_of_deleted_triangles;
+
+      std::set<Ref<Triangle> > new_triangles;
+      std::set<Ref<Edge> > new_edges;
+      std::set<Ref<Vertex> > new_vertices;
+
+      // a vertex to set-of-connected-triangles map
+      std::map<Ref<Vertex>,std::set<Ref<Triangle> > > connected_triangle_map;
+      for (it = _triangles.begin(); it != _triangles.end(); it++) {
+          Ref<Triangle> tri = *it;
+          for (j = 0; j<3; j++) {
+              Ref<Vertex> v  = tri->vertex(j);
+              connected_triangle_map[v].insert(tri);
+            }
+        }
+
+      // a vertex to set-of-connected-edges map
+      std::map<Ref<Vertex>,std::set<Ref<Edge> > > connected_edge_map;
+      for (ie = _edges.begin(); ie != _edges.end(); ie++) {
+          Ref<Edge> e = *ie;
+          for (j = 0; j<2; j++) {
+              Ref<Vertex> v = e->vertex(j);
+              connected_edge_map[v].insert(e);
+            }
+        }
+  
+      for (it = _triangles.begin(); it != _triangles.end(); it++) {
+          Ref<Triangle> tri = *it;
+          // find the heights of the vertices in tri
+          double l[3], l2[3], h[3];
+          for (j=0; j<3; j++) {
+              l[j] = tri->edge(j)->straight_length();
+              if (l[j] <= 0.0) {
+                  ExEnv::errn() << "TriangulatedSurface::"
+                       << "remove_slender_triangles: bad edge length"
+                       << endl;
+                  abort();
+                }
+              l2[j] = l[j]*l[j];
+            }
+          double y = 2.0*(l2[0]*l2[1]+l2[0]*l2[2]+l2[1]*l2[2])
+                     - l2[0]*l2[0] - l2[1]*l2[1] - l2[2]*l2[2];
+          if (y < 0.0) y = 0.0;
+          double x = 0.5*sqrt(y);
+          for (j=0; j<3; j++) h[j] = x/l[j];
+
+          // find the shortest height
+          int hmin;
+          if (h[0] < h[1]) hmin = 0;
+          else hmin = 1;
+          if (h[2] < h[hmin]) hmin = 2;
+
+          // see if the shortest height is below the cutoff
+          if (remove_slender && h[hmin] < height_cutoff) {
+              // find the vertex that gets eliminated
+              Ref<Vertex> vertex;
+              for (j=0; j<3; j++) {
+                  if (tri->vertex(j) != tri->edge(hmin)->vertex(0)
+                      &&tri->vertex(j) != tri->edge(hmin)->vertex(1)) {
+                      vertex = tri->vertex(j);
+                      break;
+                    }
+                }
+              std::set<Ref<Triangle> > connected_triangles;
+              connected_triangles = connected_triangle_map[vertex];
+
+              // if one of the connected triangles has a vertex
+              // in a deleted triangle, save this one until the
+              // next pass
+              int skip = 0;
+              for (jt = connected_triangles.begin();
+                   jt != connected_triangles.end();
+                   jt++) {
+                  Ref<Triangle> tri = *jt;
+                  for (j=0; j<3; j++) {
+                      Ref<Vertex> v = tri->vertex(j);
+                      if (vertices_of_deleted_triangles.find(v)
+                          != vertices_of_deleted_triangles.end()) {
+                          skip = 1;
+                          break;
+                        }
+                    }
+                  if (skip) break;
+                }
+              if (skip) continue;
+
+              // find all of the edges contained in the connected triangles
+              std::set<Ref<Edge> > all_edges;
+              for (jt = connected_triangles.begin();
+                   jt != connected_triangles.end();
+                   jt++) {
+                  Ref<Triangle> ctri = *jt;
+                  Ref<Edge> e0 = ctri->edge(0);
+                  Ref<Edge> e1 = ctri->edge(1);
+                  Ref<Edge> e2 = ctri->edge(2);
+                  all_edges.insert(e0);
+                  all_edges.insert(e1);
+                  all_edges.insert(e2);
+                }
+              // find all of the edges connected to the deleted vertex
+              // (including the short edge)
+              std::set<Ref<Edge> > connected_edges;
+              connected_edges = connected_edge_map[vertex];
+              // find the edges forming the perimeter of the deleted triangles
+              // (these are used to form the new triangles)
+              std::set<Ref<Edge> > perimeter_edges;
+              perimeter_edges = all_edges;
+              erase_elements_by_value(perimeter_edges,
+                                      connected_edges.begin(),
+                                      connected_edges.end());
+
+              // If deleting this point causes a flattened piece of
+              // surface, reject it.  This is tested by checking for
+              // triangles that have all vertices contained in the set
+              // of vertices connected to the deleted vertex.
+              std::set<Ref<Vertex> > connected_vertices;
+              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
+                   je++) {
+                  Ref<Edge> e = *je;
+                  Ref<Vertex> v0 = e->vertex(0);
+                  Ref<Vertex> v1 = e->vertex(1);
+                  connected_vertices.insert(v0);
+                  connected_vertices.insert(v1);
+                }
+              std::set<Ref<Triangle> > triangles_connected_to_perimeter;
+              for (jv = connected_vertices.begin();
+                   jv != connected_vertices.end();
+                   jv++) {
+                  triangles_connected_to_perimeter.insert(
+                      connected_triangle_map[*jv].begin(),
+                      connected_triangle_map[*jv].end());
+                }
+              for (jt = triangles_connected_to_perimeter.begin();
+                   jt != triangles_connected_to_perimeter.end();
+                   jt++) {
+                  Ref<Triangle> t = *jt;
+                  Ref<Vertex> v0 = t->vertex(0);
+                  Ref<Vertex> v1 = t->vertex(1);
+                  Ref<Vertex> v2 = t->vertex(2);
+                  if (connected_vertices.find(v0)!=connected_vertices.end()
+                      &&connected_vertices.find(v1)!=connected_vertices.end()
+                      &&connected_vertices.find(v2)!=connected_vertices.end())
+                      {
+                      skip = 1;
+                      break;
+                    }
+                }
+              if (skip) {
+                  continue;
+                }
+
+              deleted_triangles.insert(connected_triangles.begin(),
+                                       connected_triangles.end());
+              deleted_vertices.insert(vertex);
+              deleted_edges.insert(connected_edges.begin(),
+                                   connected_edges.end());
+
+              for (jt = deleted_triangles.begin();
+                   jt != deleted_triangles.end();
+                   jt++) {
+                  Ref<Triangle> t = *jt;
+                  for (j=0; j<2; j++) {
+                      Ref<Vertex> v = t->vertex(j);
+                      vertices_of_deleted_triangles.insert(v);
+                    }
+                }
+
+              // find a new point that replaces the deleted vertex
+              // (for now use one of the original, since it must lie on the
+              // analytic surface)
+              Ref<Vertex> replacement_vertex;
+              Ref<Edge> short_edge;
+              if (hmin==0) {
+                  if (l[1] < l[2]) short_edge = tri->edge(1);
+                  else short_edge = tri->edge(2);
+                }
+              else if (hmin==1) {
+                  if (l[0] < l[2]) short_edge = tri->edge(0);
+                  else short_edge = tri->edge(2);
+                }
+              else {
+                  if (l[0] < l[1]) short_edge = tri->edge(0);
+                  else short_edge = tri->edge(1);
+                }
+              if (short_edge->vertex(0) == vertex) {
+                  replacement_vertex = short_edge->vertex(1);
+                }
+              else {
+                  replacement_vertex = short_edge->vertex(0);
+                }
+              new_vertices.insert(replacement_vertex);
+              // for each vertex on the perimeter form a new edge to the
+              // replacement vertex (unless the replacement vertex
+              // is equal to the perimeter vertex)
+              std::map<Ref<Vertex>,Ref<Edge> > new_edge_map;
+              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
+                   je++) {
+                  Ref<Edge> e = *je;
+                  for (k = 0; k<2; k++) {
+                      Ref<Vertex> v = e->vertex(k);
+                      if (v == replacement_vertex) continue;
+                      if (new_edge_map.find(v) == new_edge_map.end()) {
+                          Ref<Edge> new_e;
+                          // if the edge already exists then use the
+                          // existing edge
+                          for (ke = perimeter_edges.begin();
+                               ke != perimeter_edges.end();
+                               ke++) {
+                              Ref<Edge> tmp = *ke;
+                              if ((tmp->vertex(0) == replacement_vertex
+                                   &&tmp->vertex(1) == v)
+                                  ||(tmp->vertex(1) == replacement_vertex
+                                     &&tmp->vertex(0) == v)) {
+                                  new_e = tmp;
+                                  break;
+                                }
+                            }
+                          if (ke == perimeter_edges.end()) {
+                              new_e = newEdge(replacement_vertex,
+                                              v);
+                            }
+                          new_edge_map[v] = new_e;
+                          new_edges.insert(new_e);
+                        }
+                    }
+                }
+              // for each edge on the perimeter form a new triangle with the
+              // replacement vertex (unless the edge contains the replacement
+              // vertex)
+              for (je = perimeter_edges.begin(); je != perimeter_edges.end();
+                   je++) {
+                  Ref<Edge> e1 = *je;
+                  Ref<Vertex> v0 = e1->vertex(0);
+                  Ref<Vertex> v1 = e1->vertex(1);
+                  Ref<Edge> e2 = new_edge_map[v0];
+                  Ref<Edge> e3 = new_edge_map[v1];
+                  if (v0 == replacement_vertex
+                      || v1 == replacement_vertex) continue;
+                  // Compute the correct orientation of e1 within the new
+                  // triangle, by finding the orientation within the old
+                  // triangle.
+                  int orientation = 0;
+                  for (kt = connected_triangles.begin();
+                       kt != connected_triangles.end();
+                       kt++) {
+                      if ((*kt)->contains(e1)) {
+                          orientation
+                              = (*kt)->orientation(e1);
+                          break;
+                        }
+                    }
+                  Ref<Triangle> newtri(newTriangle(e1,e2,e3,orientation));
+                  new_triangles.insert(newtri);
+                }
+            }
+        }
+
+#if WRITE_OOGL
+      if (_debug) {
+          char filename[100];
+          static int pass = 0;
+          sprintf(filename, "surfst%04d.oogl", pass);
+          ExEnv::outn() << scprintf("PASS = %04d\n", pass);
+          Ref<Render> render = new OOGLRender(filename);
+          Ref<RenderedPolygons> poly = new RenderedPolygons;
+          poly->initialize(_vertices.size(), _triangles.size(),
+                           RenderedPolygons::Vertex);
+          // the number of triangles and edges touching a vertex
+          int *n_triangle = new int[_vertices.size()];
+          int *n_edge = new int[_vertices.size()];
+          memset(n_triangle,0,sizeof(int)*_vertices.size());
+          memset(n_edge,0,sizeof(int)*_vertices.size());
+          std::set<Ref<Triangle> >::iterator it;
+          std::set<Ref<Edge> >::iterator ie;
+          std::set<Ref<Vertex> >::iterator iv;
+          std::map<Ref<Vertex>, int> vertex_to_index;
+          int i = 0;
+          for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
+              Ref<Vertex> v = *iv;
+              vertex_to_index[v] = i;
+              poly->set_vertex(i,
+                               v->point()[0],
+                               v->point()[1],
+                               v->point()[2]);
+              if (deleted_vertices.find(v) != deleted_vertices.end()) {
+                  poly->set_vertex_rgb(i, 1.0, 0.0, 0.0);
+                }
+              else {
+                  poly->set_vertex_rgb(i, 0.3, 0.3, 0.3);
+                }
+            }
+          i = 0;
+          for (it = _triangles.begin(); it != _triangles.end(); it++, i++) {
+              Ref<Triangle> t = *it;
+              int i0 = vertex_to_index[t->vertex(0)];
+              int i1 = vertex_to_index[t->vertex(1)];
+              int i2 = vertex_to_index[t->vertex(2)];
+              n_triangle[i0]++;
+              n_triangle[i1]++;
+              n_triangle[i2]++;
+              poly->set_face(i,i0,i1,i2);
+            }
+          for (ie = _edges.begin(); ie != _edges.end(); ie++, i++) {
+              Ref<Edge> e = *ie;
+              int i0 = vertex_to_index[e->vertex(0)];
+              int i1 = vertex_to_index[e->vertex(1)];
+              n_edge[i0]++;
+              n_edge[i1]++;
+            }
+          i = 0;
+          for (iv = _vertices.begin(); iv != _vertices.end(); iv++, i++) {
+              Ref<Vertex> v = *iv;
+              if (n_triangle[i] != n_edge[i]) {
+                  ExEnv::outn() << "found bad vertex"
+                       << " nedge = " << n_edge[i]
+                       << " ntriangle = " << n_triangle[i]
+                       << endl;
+                  if (deleted_vertices.find(v) != deleted_vertices.end()) {
+                      poly->set_vertex_rgb(i, 1.0, 1.0, 0.0);
+                    }
+                  else {
+                      poly->set_vertex_rgb(i, 0.0, 1.0, 0.0);
+                    }
+                }
+            }
+          render->render(poly.pointer());
+          pass++;
+          delete[] n_triangle;
+          delete[] n_edge;
+        }
+#endif
+
+      erase_elements_by_value(_triangles,
+                              deleted_triangles.begin(),
+                              deleted_triangles.end());
+      erase_elements_by_value(_edges,
+                              deleted_edges.begin(),
+                              deleted_edges.end());
+      erase_elements_by_value(_vertices,
+                              deleted_vertices.begin(),
+                              deleted_vertices.end());
+
+      _triangles.insert(new_triangles.begin(), new_triangles.end());
+      _edges.insert(new_edges.begin(), new_edges.end());
+      _vertices.insert(new_vertices.begin(), new_vertices.end());
+
+      if (_verbose) {
+          ExEnv::outn() << "intermediate: ";
+          topology_info();
+        }
+
+      deleted_edges_length = deleted_edges.size();
+    } while(deleted_edges_length != 0);
+
+  // fix the index maps
+  _vertex_to_index.clear();
+  _edge_to_index.clear();
+  _triangle_to_index.clear();
+
+  _index_to_vertex.clear();
+  _index_to_edge.clear();
+  _index_to_triangle.clear();
+
+  _index_to_vertex.resize(_vertices.size());
+  for (i=0, iv = _vertices.begin(); iv != _vertices.end(); i++, iv++) {
+      _vertex_to_index[*iv] = i;
+      _index_to_vertex[i] = *iv;
+    }
+
+  _index_to_edge.resize(_edges.size());
+  for (i=0, ie = _edges.begin(); ie != _edges.end(); i++, ie++) {
+      _edge_to_index[*ie] = i;
+      _index_to_edge[i] = *ie;
+    }
+
+  _index_to_triangle.resize(_triangles.size());
+  for (i=0, it = _triangles.begin(); it != _triangles.end(); i++, it++) {
+      _triangle_to_index[*it] = i;
+      _index_to_triangle[i] = *it;
+    }
+
+  // restore the int arrays if they were there to begin with
+  if (surface_was_completed) {
+      complete_int_arrays();
+      _completed_surface = 1;
+    }
+
+  if (_verbose) {
+      ExEnv::outn() << "final: ";
+      topology_info();
+    }
+}
+
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/tmplinst.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/tmplinst.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/tmplinst.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,71 @@
+//
+// tmplinst.cc - template instantations for isosurf
+//
+// Copyright (C) 1998 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef HAVE_CONFIG_H
+#include <scconfig.h>
+#endif
+#include <math/isosurf/shape.h>
+#include <math/isosurf/triangle.h>
+
+#include <vector>
+
+using namespace sc;
+
+#ifdef EXPLICIT_TEMPLATE_INSTANTIATION
+
+template class vector<Ref<Vertex> >;
+template class vector<Ref<Edge> >;
+template class vector<Ref<Triangle> >;
+
+// Vertex
+template class map<Ref<Vertex>, int>;
+template class set<Ref<Vertex> >;
+
+// Edge
+template class map<Ref<Edge>, int>;
+template class set<Ref<Edge> >;
+
+// Triangle
+template class map<Ref<Triangle>, int>;
+template class set<Ref<Triangle> >;
+
+// Shape
+template class map<Ref<Shape>, int>;
+template class set<Ref<Shape> >;
+
+// (mixed)
+template class map<Ref<Vertex>, set<Ref<Edge> > >;
+template class map<Ref<Vertex>, Ref<Edge> >;
+template class map<Ref<Vertex>, set<Ref<Triangle> > >;
+#endif
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,566 @@
+//
+// triangle.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/isosurf/triangle.h>
+#include <math/scmat/vector3.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+// Triangle
+///////////////////////////////////////////////////////////////////////////
+// Here is the layout of the vertices and edges in the triangles:         |
+//                                                                        |
+//                       vertex(1) (r=0, s=1)                             |
+//                          +                                             |
+//                         / \  _edges[1].vertex(_orientation1)           |
+//                        /   \                                           |
+//                       /     \                                          |
+//                      /       \                                         |
+//                     /         \                                        |
+//                    /           \                                       |
+//         _edges[0] /             \  _edges[1]                           |
+//          (r = 0) /               \   (1-r-s = 0)                       |
+//                 /                 \                                    |
+//                /                   \                                   |
+//               /                     \                                  |
+//              /                       \ _edges[1].vertex(!_orientation1)|
+//             /                         \                                |
+//   vertex(0)+---------------------------+                               |
+// (r=0, s=0)        _edges[2] (s = 0)      vertex(2) (r=1, s=0)          |
+//                                                                        |
+//  Zienkiewicz and Taylor, "The Finite Element Method", 4th Ed, Vol 1,   |
+//  use                                                                   |
+//      L1 = 1 - r - s                                                    |
+//      L2 = r,                                                           |
+//      L3 = s.                                                           |
+//  I also use these below.                                               |
+///////////////////////////////////////////////////////////////////////////
+
+Triangle::Triangle(const Ref<Edge>& v1, const Ref<Edge>& v2, const Ref<Edge>& v3,
+                   unsigned int orientation0)
+{
+  _orientation0 = orientation0;
+
+  _edges[0] = v1;
+  _edges[1] = v2;
+  _edges[2] = v3;
+
+  // edge 0 corresponds to r = 0
+  // edge 1 corresponds to (1-r-s) = 0
+  // edge 2 corresponds to s = 0
+  // edge 0 vertex _orientation0 is (r=0,s=0)
+  // edge 1 vertex _orientation1 is (r=0,s=1)
+  // edge 2 vertex _orientation2 is (r=1,s=0)
+  // edge 0 vertex (1-_orientation0) is edge 1 vertex _orientation1
+  // edge 1 vertex (1-_orientation1) is edge 2 vertex _orientation2
+  // edge 2 vertex (1-_orientation2) is edge 0 vertex _orientation0
+
+  Ref<Edge> *e = _edges;
+
+  // swap edges 1 and 2 if necessary
+  if (e[0]->vertex(1-_orientation0) != e[1]->vertex(0)) {
+      if (e[0]->vertex(1-_orientation0) != e[1]->vertex(1)) {
+          e[1] = v3;
+          e[2] = v2;
+        }
+    }
+
+  // compute the orientation of _edge[1]
+  if (e[0]->vertex(1-_orientation0) == e[1]->vertex(0)) {
+      _orientation1 = 0;
+    }
+  else {
+      _orientation1 = 1;
+    }
+
+  // compute the orientation of _edge[2]
+  if (e[1]->vertex(1-_orientation1) == e[2]->vertex(0)) {
+      _orientation2 = 0;
+    }
+  else {
+      _orientation2 = 1;
+    }
+
+  // make sure that the triangle is really a triangle
+  if ( e[0]->vertex(1-_orientation0) != e[1]->vertex(_orientation1)
+       || e[1]->vertex(1-_orientation1) != e[2]->vertex(_orientation2)
+       || e[2]->vertex(1-_orientation2) != e[0]->vertex(_orientation0))
+    {
+      ExEnv::errn() << "Triangle: given edges that don't form a triangle" << endl;
+      abort();
+    }
+
+  _order = 1;
+  _vertices = new Ref<Vertex>[3];
+
+  _vertices[TriInterpCoef::ijk_to_index(_order, 0, 0)] = vertex(0);
+  _vertices[TriInterpCoef::ijk_to_index(0, 0, _order)] = vertex(1);
+  _vertices[TriInterpCoef::ijk_to_index(0, _order, 0)] = vertex(2);
+}
+
+Triangle::~Triangle()
+{
+  if (_vertices) delete[] _vertices;
+}
+
+Ref<Vertex>
+Triangle::vertex(int i)
+{
+  return _edges[i]->vertex(orientation(i));
+}
+
+unsigned int
+Triangle::orientation(const Ref<Edge>& e) const
+{
+  if (e == _edges[0]) return orientation(0);
+  if (e == _edges[1]) return orientation(1);
+  return orientation(2);
+}
+
+int
+Triangle::contains(const Ref<Edge>& e) const
+{
+  if (_edges[0] == e) return 1;
+  if (_edges[1] == e) return 1;
+  if (_edges[2] == e) return 1;
+  return 0;
+}
+
+double Triangle::flat_area()
+{
+  double a = _edges[0]->straight_length();
+  double b = _edges[1]->straight_length();
+  double c = _edges[2]->straight_length();
+  double a2 = a*a;
+  double b2 = b*b;
+  double c2 = c*c;
+  return 0.25 * sqrt( 2.0 * (c2*b2 + c2*a2 + a2*b2)
+                      - c2*c2 - b2*b2 - a2*a2);
+}
+
+void Triangle::add_vertices(std::set<Ref<Vertex> >&set)
+{
+  for (int i=0; i<3; i++) set.insert(_edges[i]->vertex(orientation(i)));
+}
+
+void Triangle::add_edges(std::set<Ref<Edge> >&set)
+{
+  for (int i=0; i<3; i++) set.insert(_edges[i]);
+}
+
+void
+Triangle::interpolate(double r,double s,const Ref<Vertex>&result,SCVector3&dA)
+{
+  TriInterpCoefKey key(_order, r, s);
+  Ref<TriInterpCoef> coef = new TriInterpCoef(key);
+  interpolate(coef, r, s, result, dA);
+}
+
+void
+Triangle::interpolate(const Ref<TriInterpCoef>& coef,
+                      double r, double s,
+                      const Ref<Vertex>&result, SCVector3&dA)
+{
+  unsigned int i, j, k;
+
+  //double L1 = 1 - r - s;
+  //double L2 = r;
+  //double L3 = s;
+
+  SCVector3 tmp(0.0);
+  SCVector3 x_s(0.0);
+  SCVector3 x_r(0.0);
+  for (i=0; i<=_order; i++) {
+      for (j=0; j <= _order-i; j++) {
+          k = _order - i - j;
+          int index = TriInterpCoef::ijk_to_index(i,j,k);
+          tmp += coef->coef(i,j,k)*_vertices[index]->point();
+          x_s += coef->sderiv(i,j,k)*_vertices[index]->point();
+          x_r += coef->rderiv(i,j,k)*_vertices[index]->point();
+        }
+    }
+  result->set_point(tmp);
+
+  if (result->has_normal()) {
+      for (i=0; i<_order; i++) {
+          for (j=0; j <= _order-i; j++) {
+              k = _order - i - j;
+              int index = TriInterpCoef::ijk_to_index(i,j,k);
+              tmp += coef->coef(i,j,k)*_vertices[index]->point();
+            }
+        }
+      result->set_normal(tmp);
+    }
+
+  // Find the surface element
+  dA = x_r.cross(x_s);
+}
+
+void
+Triangle::interpolate(double r, double s,
+                      const Ref<Vertex>&result, SCVector3&dA,
+                      const Ref<Volume> &vol, double isoval)
+{
+  // set up an initial dummy norm
+  SCVector3 norm(0.0);
+  result->set_normal(norm);
+
+  // initial guess
+  interpolate(r,s,result,dA);
+
+  // now refine that guess
+  SCVector3 trialpoint = result->point();
+  SCVector3 trialnorm = result->normal();
+  SCVector3 newpoint;
+  vol->solve(trialpoint,trialnorm,isoval,newpoint);
+  // compute the true normal
+  vol->set_x(newpoint);
+  if (vol->gradient_implemented()) {
+      vol->get_gradient(trialnorm);
+    }
+  trialnorm.normalize();
+  result->set_point(newpoint);
+  result->set_normal(trialnorm);
+}
+
+void
+Triangle::flip()
+{
+  _orientation0 = _orientation0?0:1;
+  _orientation1 = _orientation1?0:1;
+  _orientation2 = _orientation2?0:1;
+}
+
+void
+Triangle::set_order(int order, const Ref<Volume>&vol, double isovalue)
+{
+  if (order > max_order) {
+      ExEnv::errn() << scprintf("Triangle::set_order: max_order = %d but order = %d\n",
+                       max_order, order);
+      abort();
+    }
+
+  unsigned int i, j, k;
+
+  if (edge(0)->order() != order
+      ||edge(1)->order() != order
+      ||edge(2)->order() != order) {
+      ExEnv::errn() << "Triangle::set_order: edge order doesn't match" << endl;
+      abort();
+    }
+
+  _order = order;
+  delete[] _vertices;
+
+  _vertices = new Ref<Vertex>[TriInterpCoef::order_to_nvertex(_order)];
+
+  // fill in the corner vertices
+  _vertices[TriInterpCoef::ijk_to_index(_order, 0, 0)] = vertex(0);
+  _vertices[TriInterpCoef::ijk_to_index(0, 0, _order)] = vertex(1);
+  _vertices[TriInterpCoef::ijk_to_index(0, _order, 0)] = vertex(2);
+
+  // fill in the interior edge vertices
+  for (i = 1; i < _order; i++) {
+      j = _order - i;
+      _vertices[TriInterpCoef::ijk_to_index(0, i, j)]
+          = _edges[1]->interior_vertex(_orientation1?i:j);
+      _vertices[TriInterpCoef::ijk_to_index(j, 0, i)]
+          = _edges[0]->interior_vertex(_orientation0?i:j);
+      _vertices[TriInterpCoef::ijk_to_index(i, j, 0)]
+          = _edges[2]->interior_vertex(_orientation2?i:j);
+    }
+
+  const SCVector3& p0 = vertex(0)->point();
+  const SCVector3& p1 = vertex(1)->point();
+  const SCVector3& p2 = vertex(2)->point();
+  const SCVector3& norm0 = vertex(0)->normal();
+  const SCVector3& norm1 = vertex(1)->normal();
+  const SCVector3& norm2 = vertex(2)->normal();
+
+  for (i=0; i<=_order; i++) {
+      double I = (1.0*i)/_order;
+      for (j=0; j<=_order-i; j++) {
+          SCVector3 trialpoint;
+          SCVector3 trialnorm;
+          SCVector3 newpoint;
+          double J = (1.0*j)/_order;
+          k = _order - i - j;
+          if (!i || !j || !k) continue; // interior point check
+          double K = (1.0*k)/_order;
+          int index = TriInterpCoef::ijk_to_index(i,j,k);
+          // this get approximate vertices and normals
+          trialpoint = I*p0 + J*p1 + K*p2;
+          trialnorm = I*norm0 + J*norm1 + K*norm2;
+          // now refine that guess
+          vol->solve(trialpoint,trialnorm,isovalue,newpoint);
+          // compute the true normal
+          vol->set_x(newpoint);
+          if (vol->gradient_implemented()) {
+              vol->get_gradient(trialnorm);
+            }
+          trialnorm.normalize();
+          _vertices[index] = new Vertex(newpoint,trialnorm);
+        }
+    }
+}
+
+/////////////////////////////////////////////////////////////////////////
+// TriangleIntegrator
+
+static ClassDesc TriangleIntegrator_cd(
+  typeid(TriangleIntegrator),"TriangleIntegrator",1,"public DescribedClass",
+  0, create<TriangleIntegrator>, 0);
+
+TriangleIntegrator::TriangleIntegrator(int order):
+  _n(order)
+{
+  _r = new double [_n];
+  _s = new double [_n];
+  _w = new double [_n];
+  coef_ = 0;
+}
+
+TriangleIntegrator::TriangleIntegrator(const Ref<KeyVal>& keyval)
+{
+  _n = keyval->intvalue("n");
+  if (keyval->error() != KeyVal::OK) {
+      _n = 7;
+    }
+  _r = new double [_n];
+  _s = new double [_n];
+  _w = new double [_n];
+  coef_ = 0;
+}
+
+TriangleIntegrator::~TriangleIntegrator()
+{
+  delete[] _r;
+  delete[] _s;
+  delete[] _w;
+  clear_coef();
+}
+
+void
+TriangleIntegrator::set_n(int n)
+{
+  delete[] _r;
+  delete[] _s;
+  delete[] _w;
+  _n = n;
+  _r = new double [_n];
+  _s = new double [_n];
+  _w = new double [_n];
+  clear_coef();
+}
+
+void
+TriangleIntegrator::set_w(int i,double w)
+{
+  _w[i] = w;
+}
+
+void
+TriangleIntegrator::set_r(int i,double r)
+{
+  _r[i] = r;
+}
+
+void
+TriangleIntegrator::set_s(int i,double s)
+{
+  _s[i] = s;
+}
+
+void
+TriangleIntegrator::init_coef()
+{
+  int i, j;
+
+  clear_coef();
+
+  coef_ = new Ref<TriInterpCoef>*[Triangle::max_order];
+  for (i=0; i<Triangle::max_order; i++) {
+      coef_[i] = new Ref<TriInterpCoef>[_n];
+      for (j=0; j<_n; j++) {
+          TriInterpCoefKey key(i+1,_r[j],_s[j]);
+          coef_[i][j] = new TriInterpCoef(key);
+        }
+    }
+}
+
+void
+TriangleIntegrator::clear_coef()
+{
+  int i, j;
+
+  if (coef_) {
+      for (i=0; i<Triangle::max_order; i++) {
+          for (j=0; j<_n; j++) {
+              coef_[i][j] = 0;
+            }
+          delete[] coef_[i];
+        }
+    }
+  delete[] coef_;
+  coef_ = 0;
+}
+
+/////////////////////////////////////////////////////////////////////////
+// GaussTriangleIntegrator
+
+static ClassDesc GaussTriangleIntegrator_cd(
+  typeid(GaussTriangleIntegrator),"GaussTriangleIntegrator",1,"public TriangleIntegrator",
+  0, create<GaussTriangleIntegrator>, 0);
+
+GaussTriangleIntegrator::GaussTriangleIntegrator(const Ref<KeyVal>& keyval):
+  TriangleIntegrator(keyval)
+{
+  ExEnv::out0() << "Created a GaussTriangleIntegrator with n = " << n() << endl;
+  init_rw(n());
+  init_coef();
+}
+
+GaussTriangleIntegrator::GaussTriangleIntegrator(int order):
+  TriangleIntegrator(order)
+{
+  init_rw(n());
+  init_coef();
+}
+
+void
+GaussTriangleIntegrator::set_n(int n)
+{
+  TriangleIntegrator::set_n(n);
+  init_rw(n);
+  init_coef();
+}
+
+void
+GaussTriangleIntegrator::init_rw(int order)
+{
+  if (order == 1) {
+      set_r(0, 1.0/3.0);
+      set_s(0, 1.0/3.0);
+      set_w(0, 1.0);
+    }
+  else if (order == 3) {
+      set_r(0, 1.0/6.0);
+      set_r(1, 2.0/3.0);
+      set_r(2, 1.0/6.0);
+      set_s(0, 1.0/6.0);
+      set_s(1, 1.0/6.0);
+      set_s(2, 2.0/3.0);
+      set_w(0, 1.0/3.0);
+      set_w(1, 1.0/3.0);
+      set_w(2, 1.0/3.0);
+    }
+  else if (order == 4) {
+      set_r(0, 1.0/3.0);
+      set_r(1, 1.0/5.0);
+      set_r(2, 3.0/5.0);
+      set_r(3, 1.0/5.0);
+      set_s(0, 1.0/3.0);
+      set_s(1, 1.0/5.0);
+      set_s(2, 1.0/5.0);
+      set_s(3, 3.0/5.0);
+      set_w(0, -27.0/48.0);
+      set_w(1, 25.0/48.0);
+      set_w(2, 25.0/48.0);
+      set_w(3, 25.0/48.0);
+    }
+  else if (order == 6) {
+      set_r(0, 0.091576213509771);
+      set_r(1, 0.816847572980459);
+      set_r(2, r(0));
+      set_r(3, 0.445948490915965);
+      set_r(4, 0.108103018168070);
+      set_r(5, r(3));
+      set_s(0, r(0));
+      set_s(1, r(0));
+      set_s(2, r(1));
+      set_s(3, r(3));
+      set_s(4, r(3));
+      set_s(5, r(4));
+      set_w(0, 0.109951743655322);
+      set_w(1, w(0));
+      set_w(2, w(0));
+      set_w(3, 0.223381589678011);
+      set_w(4, w(3));
+      set_w(5, w(3));
+    }
+  else if (order == 7) {
+      set_r(0, 1.0/3.0);
+      set_r(1, 0.101286507323456);
+      set_r(2, 0.797426985353087);
+      set_r(3, r(1));
+      set_r(4, 0.470142064105115);
+      set_r(5, 0.059715871789770);
+      set_r(6, r(4));
+      set_s(0, r(0));
+      set_s(1, r(1));
+      set_s(2, r(1));
+      set_s(3, r(2));
+      set_s(4, r(4));
+      set_s(5, r(4));
+      set_s(6, r(5));
+      set_w(0, 0.225);
+      set_w(1, 0.125939180544827);
+      set_w(2, w(1));
+      set_w(3, w(1));
+      set_w(4, 0.132394152788506);
+      set_w(5, w(4));
+      set_w(6, w(4));
+    }
+  else {
+      ExEnv::errn() << "GaussTriangleIntegrator: invalid order " << order << endl;
+      abort();
+    }
+
+  // scale the weights by the area of the unit triangle
+  for (int i=0; i<order; i++) {
+      set_w(i, w(i)*0.5);
+    }
+}
+
+GaussTriangleIntegrator::~GaussTriangleIntegrator()
+{
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/triangle.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,144 @@
+//
+// triangle.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_triangle_h
+#define _math_isosurf_triangle_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <math/isosurf/tricoef.h>
+#include <math/isosurf/edge.h>
+
+namespace sc {
+
+class Triangle: public RefCount {
+  protected:
+    // these break gcc 2.5.8
+    //unsigned int _order:5;
+    //unsigned int _orientation0:1;
+    //unsigned int _orientation1:1;
+    //unsigned int _orientation2:1;
+    //unsigned char _order;
+    //unsigned char _orientation0;
+    //unsigned char _orientation1;
+    //unsigned char _orientation2;
+    unsigned int _order;
+    unsigned int _orientation0;
+    unsigned int _orientation1;
+    unsigned int _orientation2;
+    Ref<Edge> _edges[3];
+    Ref<Vertex> *_vertices;
+  public:
+    enum {max_order = 10};
+
+    Triangle(const Ref<Edge>& v1, const Ref<Edge>& v2, const Ref<Edge>& v3,
+             unsigned int orient0 = 0);
+    Ref<Edge> edge(int i) { return _edges[i]; };
+    int contains(const Ref<Edge>&) const;
+    unsigned int orientation(int i) const
+    {
+      return i==0?_orientation0:i==1?_orientation1:_orientation2;
+    }
+    unsigned int orientation(const Ref<Edge>&) const;
+    ~Triangle();
+    void add_edges(std::set<Ref<Edge> >&);
+    void add_vertices(std::set<Ref<Vertex> >&);
+
+    // returns the surface area element
+    // 0<=r<=1, 0<=s<=1, 0<=r+s<=1
+    // Ref<Vertex> is the intepolated vertex (both point and normal)
+    void interpolate(const Ref<TriInterpCoef>&,
+                     double r,double s,const Ref<Vertex>&v, SCVector3& dA);
+    void interpolate(double r,double s,const Ref<Vertex>&v, SCVector3& dA);
+    void interpolate(double r,double s,const Ref<Vertex>&v, SCVector3& dA,
+                     const Ref<Volume> &vol, double isovalue);
+
+    // returns a corner vertex from the triangle
+    // i = 0 is the (0,0) vertex (or L1 = 1, L2 = 0, L3 = 0)
+    // i = 1 is the (r=1,s=0) vertex (or L1 = 0, L2 = 1, L3 = 0)
+    // i = 2 is the (r=0,s=1) vertex (or L1 = 0, L2 = 0, L3 = 1)
+    Ref<Vertex> vertex(int i);
+
+    double flat_area();
+
+    // flip the orientation
+    void flip();
+
+    unsigned int order() const { return _order; }
+
+    void set_order(int order, const Ref<Volume>&vol,double isovalue);
+};
+
+
+
+class TriangleIntegrator: public DescribedClass {
+  private:
+    int _n;
+    double* _r;
+    double* _s;
+    double* _w;
+    // precomputed interpolation coefficients for triangles of various orders
+    Ref<TriInterpCoef> **coef_; // max_order by _n
+  protected:
+    void set_r(int i,double r);
+    void set_s(int i,double s);
+    void set_w(int i,double w);
+    void init_coef();
+    void clear_coef();
+  public:
+    TriangleIntegrator(const Ref<KeyVal>&);
+    TriangleIntegrator(int n);
+    virtual ~TriangleIntegrator();
+    inline double w(int i) { return _w[i]; }
+    inline double r(int i) { return _r[i]; }
+    inline double s(int i) { return _s[i]; }
+    inline int n() { return _n; }
+    virtual void set_n(int n);
+    Ref<TriInterpCoef> coef(int order, int i) { return coef_[order-1][i]; }
+};
+
+
+class GaussTriangleIntegrator: public TriangleIntegrator {
+  private:
+    void init_rw(int order);
+  public:
+    GaussTriangleIntegrator(const Ref<KeyVal>&);
+    GaussTriangleIntegrator(int order);
+    ~GaussTriangleIntegrator();
+    void set_n(int n);
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,129 @@
+//
+// tricoef.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <math/isosurf/triangle.h>
+#include <math/isosurf/tricoef.h>
+
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+// Utility functions
+
+static inline void
+init_coef_deriv(double L, int order, double *Lcoef, double *Lcoefderiv)
+{
+  int i;
+  Lcoef[0] = 1.0;
+  Lcoefderiv[0] = 0.0;
+  double spacing = 1.0/order;
+  for (i=1; i<=order; i++) {
+      Lcoef[i] = Lcoef[i-1] * (L - (i-1)*spacing)/(i*spacing);
+      Lcoefderiv[i] = Lcoefderiv[i-1] * (L - (i-1)*spacing)/(i*spacing)
+                      + Lcoef[i-1]/(i*spacing);
+    }
+}
+
+
+/////////////////////////////////////////////////////////////////////////
+// The TriInterpCoef Utility Class
+
+TriInterpCoef::TriInterpCoef(const TriInterpCoefKey& key)
+{
+  int i,j,k;
+
+  int order = key.order();
+  double L1 = key.L1();
+  double L2 = key.L2();
+  double L3 = key.L3();
+  int n = order_to_nvertex(order);
+  coef_ = new double[n];
+  rderiv_ = new double[n];
+  sderiv_ = new double[n];
+
+  double L1coef[Triangle::max_order+1];
+  double L2coef[Triangle::max_order+1];
+  double L3coef[Triangle::max_order+1];
+
+  double L1coefderiv[Triangle::max_order+1];
+  double L2coefderiv[Triangle::max_order+1];
+  double L3coefderiv[Triangle::max_order+1];
+
+  init_coef_deriv(L1, order, L1coef, L1coefderiv);
+  init_coef_deriv(L2, order, L2coef, L2coefderiv);
+  init_coef_deriv(L3, order, L3coef, L3coefderiv);
+
+  // the r derivatives
+  double L1coef_r[Triangle::max_order+1];
+  double L2coef_r[Triangle::max_order+1];
+  double L3coef_r[Triangle::max_order+1];
+
+  // the s derivatives
+  double L1coef_s[Triangle::max_order+1];
+  double L2coef_s[Triangle::max_order+1];
+  double L3coef_s[Triangle::max_order+1];
+
+  // convert into r and s derivatives
+  for (i=0; i<=order; i++) {
+      L1coef_r[i] = -L1coefderiv[i];
+      L1coef_s[i] = -L1coefderiv[i];
+      L2coef_r[i] =  L2coefderiv[i];
+      L2coef_s[i] =  0.0;
+      L3coef_r[i] =  0.0;
+      L3coef_s[i] =  L3coefderiv[i];
+    }
+
+  for (i=0; i<=order; i++) {
+      for (j=0; j <= order-i; j++) {
+          k = order - i - j;
+          coef(i,j,k) = L1coef[i]*L2coef[j]*L3coef[k];
+          sderiv(i,j,k) = L1coef_s[i]*L2coef[j]*L3coef[k]
+                         +L1coef[i]*L2coef_s[j]*L3coef[k]
+                         +L1coef[i]*L2coef[j]*L3coef_s[k];
+          rderiv(i,j,k) = L1coef_r[i]*L2coef[j]*L3coef[k]
+                         +L1coef[i]*L2coef_r[j]*L3coef[k]
+                         +L1coef[i]*L2coef[j]*L3coef_r[k];
+        }
+    }
+}
+
+TriInterpCoef::~TriInterpCoef()
+{
+  delete[] coef_;
+  delete[] rderiv_;
+  delete[] sderiv_;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/tricoef.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,98 @@
+//
+// tricoef.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_tricoef_h
+#define _math_isosurf_tricoef_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <util/ref/ref.h>
+
+namespace sc {
+
+class TriInterpCoefKey {
+  private:
+    int order_;
+    double L2_;
+    double L3_;
+  public:
+    TriInterpCoefKey(int order, double L2, double L3):
+      order_(order), L2_(L2), L3_(L3) {}
+    int order() const { return order_; }
+    double L1() const { return 1.0 - L2_ - L3_; }
+    double L2() const { return L2_; }
+    double L3() const { return L3_; }
+    int cmp(const TriInterpCoefKey&t) const {
+        if (order_ < t.order_) return -1;
+        if (order_ > t.order_) return 1;
+        if (L2_ < t.L2_) return -1;
+        if (L2_ > t.L2_) return 1;
+        if (L3_ < t.L3_) return -1;
+        if (L3_ > t.L3_) return 1;
+        return 0;
+      }
+};
+
+#define TriInterpCoefKeyEQ(k1,k2) ((k1).cmp(k2)==0)
+#define TriInterpCoefKeyCMP(k1,k2) ((k1).cmp(k2))
+
+class TriInterpCoef: public RefCount {
+    double *coef_;
+    double *rderiv_;
+    double *sderiv_;
+  public:
+    TriInterpCoef(const TriInterpCoefKey& key);
+    ~TriInterpCoef();
+    double& coef(int i, int j, int k) {return coef_[ijk_to_index(i,j,k)];}
+    double& rderiv(int i, int j, int k) {return rderiv_[ijk_to_index(i,j,k)];}
+    double& sderiv(int i, int j, int k) {return sderiv_[ijk_to_index(i,j,k)];}
+
+    static int
+    ijk_to_index(int i, int j, int k)
+    {
+      int n = i + j + k;
+      int ir = n - i;
+      return (ir*(ir+1)>>1) + j;
+    }
+
+    static int
+    order_to_nvertex(int order)
+    {
+      return ((order+1)*(order+2)>>1);
+    }
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,116 @@
+//
+// vertex.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <util/misc/formio.h>
+#include <math/isosurf/vertex.h>
+
+using namespace std;
+using namespace sc;
+
+/////////////////////////////////////////////////////////////////////////
+// Vertex (a point and a gradient)
+
+Vertex::Vertex(const SCVector3&point):
+  _point(point),
+  _normal(0)
+{
+}
+
+Vertex::Vertex(const SCVector3&point,const SCVector3&normal):
+  _point(point),
+  _normal(new SCVector3(normal))
+{
+  double dot = _normal->dot(*_normal);
+  if (dot < 0.999999 || dot > 1.000001) {
+      ExEnv::outn() << "Vertex: ctor: bad normal\n" << endl;
+      abort();
+    }
+}
+
+Vertex::Vertex():
+  _normal(0)
+{
+}
+
+Vertex::~Vertex()
+{
+  if (_normal) delete _normal;
+}
+
+void
+Vertex::set_point(const SCVector3&p)
+{
+  _point = p;
+}
+
+void
+Vertex::set_normal(const SCVector3&p)
+{
+  if (_normal) {
+      *_normal = p;
+    }
+  else {
+      _normal = new SCVector3(p);
+    }
+  double dot = _normal->dot(*_normal);
+  if (dot < 0.999999 || dot > 1.000001) {
+      ExEnv::outn() << "Vertex::set_normal: bad normal\n" << endl;
+      abort();
+    }
+}
+
+Vertex::operator SCVector3&()
+{
+  return _point;
+}
+
+void
+Vertex::print(ostream&o)
+{
+  int i;
+  o << indent << "Vertex:";
+  for (i=0; i<3; i++)  {
+      o << scprintf(" %8.5f", _point[i]);
+    }
+  if (_normal) {
+      for (i=0; i<3; i++)  {
+          o << scprintf(" %8.5f", normal()[i]);
+        }
+    }
+  o << endl;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/vertex.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,70 @@
+//
+// vertex.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_vertex_h
+#define _math_isosurf_vertex_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <stdlib.h>
+#include <math.h>
+#include <util/ref/ref.h>
+#include <math/scmat/matrix.h>
+#include <math/isosurf/volume.h>
+#include <math/isosurf/implicit.h>
+
+namespace sc {
+
+class Vertex: public RefCount {
+  private:
+    SCVector3 _point;
+    SCVector3 *_normal; // _normal is optional
+  public:
+    Vertex();
+    Vertex(const SCVector3& point,const SCVector3& normal);
+    Vertex(const SCVector3& point);
+    ~Vertex();
+    const SCVector3& point() const { return _point; }
+    int has_normal() const { return _normal != 0; }
+    const SCVector3& normal() const { return *_normal; }
+    void set_point(const SCVector3&p);
+    void set_normal(const SCVector3&p);
+    operator SCVector3&();
+
+    void print(std::ostream&o=ExEnv::out0());
+};
+
+}
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/volume.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/volume.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/volume.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,211 @@
+//
+// volume.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#include <util/misc/math.h>
+#include <util/misc/formio.h>
+#include <util/keyval/keyval.h>
+#include <math/scmat/vector3.h>
+#include <math/scmat/local.h>
+#include <math/isosurf/volume.h>
+
+using namespace std;
+using namespace sc;
+
+static ClassDesc Volume_cd(
+  typeid(Volume),"Volume",1,"public Function",
+  0, 0, 0);
+
+Volume::Volume():
+  _interp_acc(1.0e-6)
+{
+  set_dimension(new SCDimension(3));
+}
+
+Volume::Volume(const Ref<KeyVal>&keyval):
+  Function(keyval)
+{
+  set_dimension(new SCDimension(3));
+  _interp_acc = keyval->doublevalue("interpolation_accuracy");
+  if (keyval->error() != KeyVal::OK) _interp_acc = 1.0e-6;
+}
+
+Volume::~Volume()
+{
+}
+
+// interpolate using the bisection algorithm
+void
+Volume::interpolate(const SCVector3& A,
+                    const SCVector3& B,
+                    double val,
+                    SCVector3& result)
+{
+  set_x(A);
+  double value0 = value() - val;
+
+  set_x(B);
+  double value1 = value() - val;
+
+  if (value0*value1 > 0.0) {
+      failure("interpolate(): values at endpoints don't bracket val");
+    }
+  else if (value0 == 0.0) {
+      result = A;
+      return;
+    }
+  else if (value1 == 0.0) {
+      result = B;
+      return;
+    }
+
+  SCVector3 BA = B - A;
+
+  double length = sqrt(BA.dot(BA));
+  int niter = (int) (log(length/_interp_acc)/M_LN2);
+
+  double f0 = 0.0;
+  double f1 = 1.0;
+  double fnext = 0.5;
+
+  for (int i=0; i<niter; i++) {
+      SCVector3 X = A + fnext*BA;
+      set_x(X);
+      double valuenext = value() - val;
+
+      if (valuenext*value0 <= 0.0) {
+          value1 = valuenext;
+          f1 = fnext;
+          fnext = (f0 + fnext)*0.5;
+        }
+      else {
+          value0 = valuenext;
+          f0 = fnext;
+          fnext = (fnext + f1)*0.5;
+        }
+    }
+
+  result = A + fnext*BA;
+}
+
+void
+Volume::solve(const SCVector3& start,
+              const SCVector3& grad,
+              double val,
+              SCVector3& result)
+{
+  double direction;
+  set_x(start);
+  double startvalue = value();
+  if (startvalue == val) {
+      result = start;
+      return;
+    }
+  else if (startvalue < val) direction = 1.0;
+  else direction = -1.0;
+  int i=0;
+  SCVector3 next;
+  double trialvalue;
+  do {
+      if (i>10) {
+          ExEnv::errn() << "Volume::solve: couldn't find end points" << endl;
+          abort();
+        }
+      i++;
+      next = start + (direction*i)*grad;
+      set_x(next);
+      trialvalue = value();
+    } while ((startvalue-val)*(trialvalue-val)>0.0);
+  interpolate(start,next,val,result);
+}
+
+void
+Volume::failure(const char * msg)
+{
+  ExEnv::errn() << scprintf("Volume::failure: \"%s\"\n",msg);
+  abort();
+}
+
+void
+Volume::set_gradient(const SCVector3& g)
+{
+  RefSCVector grad(dimension(), matrixkit());
+  grad.set_element(0, g[0]);
+  grad.set_element(1, g[1]);
+  grad.set_element(2, g[2]);
+  set_gradient(grad);
+}
+
+void
+Volume::set_gradient(RefSCVector& g)
+{
+  Function::set_gradient(g);
+}
+
+void
+Volume::get_gradient(SCVector3& g)
+{
+  const RefSCVector v = gradient();
+  g[0] = v.get_element(0);
+  g[1] = v.get_element(1);
+  g[2] = v.get_element(2);
+}
+
+void
+Volume::set_x(const SCVector3& x)
+{
+  RefSCVector xx(dimension(), matrixkit());
+  xx.set_element(0, x[0]);
+  xx.set_element(1, x[1]);
+  xx.set_element(2, x[2]);
+  set_x(xx);
+}
+
+void
+Volume::set_x(const RefSCVector& x)
+{
+  Function::set_x(x);
+}
+
+void
+Volume::get_x(SCVector3& x)
+{
+  const RefSCVector& v = get_x_no_copy();
+  x[0] = v.get_element(0);
+  x[1] = v.get_element(1);
+  x[2] = v.get_element(2);
+}
+
+/////////////////////////////////////////////////////////////////////////////
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/volume.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/volume.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/volume.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,91 @@
+//
+// volume.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifndef _math_isosurf_volume_h
+#define _math_isosurf_volume_h
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#include <math/optimize/function.h>
+#include <math/scmat/matrix.h>
+#include <math/scmat/vector3.h>
+
+namespace sc {
+
+/** A Volume is a Function of three variables. */
+class Volume: public Function {
+  private:
+    double _interp_acc;
+  protected:
+    double& interpolation_accuracy();
+
+    virtual void compute() = 0;
+
+    virtual void failure(const char*);
+  public:
+    Volume();
+    Volume(const Ref<KeyVal>&);
+    ~Volume();
+
+    void set_gradient(const SCVector3& g);
+    void set_gradient(RefSCVector& g);
+    void get_gradient(SCVector3& g);
+    void set_x(const SCVector3& x);
+    void set_x(const RefSCVector& x);
+    void get_x(SCVector3& x);
+
+    // find the corners of a bounding box which approximately
+    // contains all points with a value between valuemin and valuemax
+    // the result must satisfy p1[i] < p2[i]
+    virtual void boundingbox(double valuemin,
+                             double valuemax,
+                             SCVector3& p1, SCVector3& p2) = 0;
+
+    virtual void interpolate(const SCVector3& p1,
+                             const SCVector3& p2,
+                             double value,
+                             SCVector3& result);
+    virtual void solve(const SCVector3& p,
+                       const SCVector3& grad,
+                       double value,
+                       SCVector3& result);
+};
+
+}
+
+#ifdef INLINE_FUNCTIONS
+#include <math/isosurf/volume_i.h>
+#endif
+
+#endif
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.cc
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.cc	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,35 @@
+//
+// volume_i.cc
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma implementation
+#endif
+
+#undef INLINE_FUNCTIONS
+
+#include <math/isosurf/volume.h>
+#include <math/isosurf/volume_i.h>
Index: ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.h
===================================================================
--- ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
+++ ThirdParty/mpqc_open/src/lib/math/isosurf/volume_i.h	(revision 860145ede0311b3539e419277d8ca4a3091f40b2)
@@ -0,0 +1,51 @@
+//
+// volume_i.h
+//
+// Copyright (C) 1996 Limit Point Systems, Inc.
+//
+// Author: Curtis Janssen <cljanss@limitpt.com>
+// Maintainer: LPS
+//
+// This file is part of the SC Toolkit.
+//
+// The SC Toolkit is free software; you can redistribute it and/or modify
+// it under the terms of the GNU Library General Public License as published by
+// the Free Software Foundation; either version 2, or (at your option)
+// any later version.
+//
+// The SC Toolkit is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU Library General Public License for more details.
+//
+// You should have received a copy of the GNU Library General Public License
+// along with the SC Toolkit; see the file COPYING.LIB.  If not, write to
+// the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+//
+// The U.S. Government is granted a limited license as per AL 91-7.
+//
+
+#ifdef __GNUC__
+#pragma interface
+#endif
+
+#ifdef INLINE_FUNCTIONS
+#define INLINE inline
+#else
+#define INLINE
+#endif
+
+namespace sc {
+
+INLINE double&
+Volume::interpolation_accuracy()
+{
+  return _interp_acc;
+};
+
+}
+
+// Local Variables:
+// mode: c++
+// c-file-style: "CLJ"
+// End:
