svn merge ^/branches/released . --username marha

18 files changed, 6949 insertions, 0 deletions

diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.cc b/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.cc
new file mode 100644
index 000000000..b414f535f
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.cc
@@ -0,0 +1,260 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+**
+** http://oss.sgi.com/projects/FreeB
+**
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+**
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+**
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <math.h>
+#include "bezierEval.h"
+
+#ifdef __WATCOMC__
+#pragma warning 14  10
+#endif
+
+#define TOLERANCE 0.0001
+
+#ifndef MAX_ORDER
+#define MAX_ORDER 16
+#endif
+
+#ifndef MAX_DIMENSION
+#define MAX_DIMENSION 4
+#endif
+
+static void normalize(float vec[3]);
+static void crossProduct(float x[3], float y[3], float ret[3]);
+#if 0 // UNUSED
+static void bezierCurveEvalfast(float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retpoint[]);
+#endif
+
+static float binomialCoefficients[8][8] = {
+  {1,0,0,0,0,0,0,0},
+  {1,1,0,0,0,0,0,0},
+  {1,2,1,0,0,0,0,0},
+  {1,3,3,1,0,0,0,0},
+  {1,4,6,4,1,0,0,0},
+  {1,5,10,10,5,1,0,0},
+  {1,6,15,20,15,6,1,0},
+  {1,7,21,35,35,21,7,1}
+};
+
+void bezierCurveEval(float u0, float u1, int order, float *ctlpoints, int stride, int dimension, float u, float retpoint[])
+{
+  float uprime = (u-u0)/(u1-u0);
+  float *ctlptr = ctlpoints;
+  float oneMinusX = 1.0f-uprime;
+  float XPower = 1.0f;
+
+  int i,k;
+  for(k=0; k<dimension; k++)
+    retpoint[k] = (*(ctlptr + k));
+
+  for(i=1; i<order; i++){
+    ctlptr += stride;
+    XPower *= uprime;
+    for(k=0; k<dimension; k++) {
+      retpoint[k] = retpoint[k]*oneMinusX + ctlptr[k]* binomialCoefficients[order-1][i] * XPower;
+    }
+  }
+}
+
+
+#if 0 // UNUSED
+/*order = degree +1 >=1.
+ */
+void bezierCurveEvalfast(float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retpoint[])
+{
+  float uprime = (u-u0)/(u1-u0);
+  float buf[MAX_ORDER][MAX_ORDER][MAX_DIMENSION];
+  float* ctlptr = ctlpoints;
+  int r, i,j;
+  for(i=0; i<order; i++) {
+    for(j=0; j<dimension; j++)
+      buf[0][i][j] = ctlptr[j];
+    ctlptr += stride;
+  }
+  for(r=1; r<order; r++){
+    for(i=0; i<order-r; i++) {
+      for(j=0; j<dimension; j++)
+	buf[r][i][j] = (1-uprime)*buf[r-1][i][j] + uprime*buf[r-1][i+1][j];
+    }
+  }
+
+  for(j=0; j<dimension; j++)
+    retpoint[j] = buf[order-1][0][j];
+}
+#endif
+
+
+/*order = degree +1 >=1.
+ */
+void bezierCurveEvalDer(float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retDer[])
+{
+  int i,k;
+  float width = u1-u0;
+  float *ctlptr = ctlpoints;
+
+  float buf[MAX_ORDER][MAX_DIMENSION];
+  if(order == 1){
+    for(k=0; k<dimension; k++)
+      retDer[k]=0;
+  }
+  for(i=0; i<order-1; i++){
+    for(k=0; k<dimension; k++) {
+      buf[i][k] = (ctlptr[stride+k] - ctlptr[k])*(order-1)/width;
+    }
+    ctlptr += stride;
+  }
+
+  bezierCurveEval(u0, u1, order-1, (float*) buf, MAX_DIMENSION,  dimension, u, retDer);
+}
+
+void bezierCurveEvalDerGen(int der, float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retDer[])
+{
+  int i,k,r;
+  float *ctlptr = ctlpoints;
+  float width=u1-u0;
+  float buf[MAX_ORDER][MAX_ORDER][MAX_DIMENSION];
+  if(der<0) der=0;
+  for(i=0; i<order; i++){
+    for(k=0; k<dimension; k++){
+      buf[0][i][k] = ctlptr[k];
+    }
+    ctlptr += stride;
+  }
+
+
+  for(r=1; r<=der; r++){
+    for(i=0; i<order-r; i++){
+      for(k=0; k<dimension; k++){
+	buf[r][i][k] = (buf[r-1][i+1][k] - buf[r-1][i][k])*(order-r)/width;
+      }
+    }
+  }
+
+  bezierCurveEval(u0, u1, order-der, (float *) (buf[der]), MAX_DIMENSION, dimension, u, retDer);
+}
+
+/*the Bezier bivarite polynomial is:
+ * sum[i:0,uorder-1][j:0,vorder-1] { ctlpoints[i*ustride+j*vstride] * B(i)*B(j)
+ * where B(i) and B(j) are basis functions
+ */
+void bezierSurfEvalDerGen(int uder, int vder, float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float ret[])
+{
+  int i;
+  float newPoints[MAX_ORDER][MAX_DIMENSION];
+
+  for(i=0; i<uorder; i++){
+
+    bezierCurveEvalDerGen(vder, v0, v1, vorder, ctlpoints+ustride*i, vstride, dimension, v, newPoints[i]);
+
+  }
+
+  bezierCurveEvalDerGen(uder, u0, u1, uorder, (float *) newPoints, MAX_DIMENSION, dimension, u, ret);
+}
+
+
+/*division by w is performed*/
+void bezierSurfEval(float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float ret[])
+{
+  bezierSurfEvalDerGen(0, 0, u0, u1, uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u, v, ret);
+  if(dimension == 4) /*homogeneous*/{
+    ret[0] /= ret[3];
+    ret[1] /= ret[3];
+    ret[2] /= ret[3];
+  }
+}
+
+void bezierSurfEvalNormal(float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float retNormal[])
+{
+  float partialU[4];
+  float partialV[4];
+  assert(dimension>=3 && dimension <=4);
+  bezierSurfEvalDerGen(1,0, u0, u1, uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u, v, partialU);
+  bezierSurfEvalDerGen(0,1, u0, u1, uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u, v, partialV);
+
+  if(dimension == 3){/*inhomogeneous*/
+    crossProduct(partialU, partialV, retNormal);
+
+    normalize(retNormal);
+
+    return;
+  }
+  else { /*homogeneous*/
+    float val[4]; /*the point coordinates (without derivative)*/
+    float newPartialU[MAX_DIMENSION];
+    float newPartialV[MAX_DIMENSION];
+    int i;
+    bezierSurfEvalDerGen(0,0, u0, u1, uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u, v, val);
+
+    for(i=0; i<=2; i++){
+      newPartialU[i] = partialU[i] * val[3] - val[i] * partialU[3];
+      newPartialV[i] = partialV[i] * val[3] - val[i] * partialV[3];
+    }
+    crossProduct(newPartialU, newPartialV, retNormal);
+    normalize(retNormal);
+  }
+}
+
+/*if size is 0, then nothing is done*/
+static void normalize(float vec[3])
+{
+  float size = (float)sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
+
+  if(size < TOLERANCE)
+    {
+#ifdef DEBUG
+      fprintf(stderr, "Warning: in oglBSpline.c normal is 0\n");
+#endif
+      return;
+    }
+  else {
+    vec[0] = vec[0]/size;
+    vec[1] = vec[1]/size;
+    vec[2] = vec[2]/size;
+  }
+}
+
+
+static void crossProduct(float x[3], float y[3], float ret[3])
+{
+  ret[0] = x[1]*y[2] - y[1]*x[2];
+  ret[1] = x[2]*y[0] - y[2]*x[0];
+  ret[2] = x[0]*y[1] - y[0]*x[1];
+
+}
+
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.h b/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.h
new file mode 100644
index 000000000..93e8c69d1
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierEval.h
@@ -0,0 +1,48 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+/*
+*/
+
+#ifndef _BEZIEREVAL_H
+#define _BEZIEREVAL_H
+
+void bezierCurveEval(float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retpoint[]);
+void bezierCurveEvalDer(float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retDer[]);
+void bezierCurveEvalDerGen(int der, float u0, float u1, int order, float *ctlpoints, int stride,  int dimension, float u, float retDer[]);
+
+
+void bezierSurfEvalDerGen(int uder, int vder, float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float ret[]);
+
+void bezierSurfEval(float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float ret[]);
+
+void bezierSurfEvalNormal(float u0, float u1, int uorder, float v0, float v1, int vorder, int dimension, float *ctlpoints, int ustride, int vstride, float u, float v, float retNormal[]);
+
+
+#endif
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.cc b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.cc
new file mode 100644
index 000000000..dbab3a1a2
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.cc
@@ -0,0 +1,206 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+** 
+** http://oss.sgi.com/projects/FreeB
+** 
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+** 
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+** 
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include "gluos.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <GL/glu.h> /*for drawing bzier patch*/
+#include "bezierPatch.h"
+#include "bezierEval.h"
+
+/*
+ *allocate an instance of bezierPatch. The control points are unknown. But
+ *the space of this array is allocated with size of 
+ *   uorder*vorder*dimension
+ *
+ */
+bezierPatch* bezierPatchMake(float umin, float vmin, float umax, float vmax, int uorder, int vorder, int dimension)
+{
+  bezierPatch* ret = (bezierPatch*) malloc(sizeof(bezierPatch));
+  assert(ret);
+  ret->umin = umin;
+  ret->vmin = vmin;
+  ret->umax = umax;
+  ret->vmax = vmax;
+  ret->uorder = uorder;
+  ret->vorder = vorder;
+  ret->dimension = dimension;
+  ret->ctlpoints = (float*) malloc(sizeof(float) * dimension * uorder * vorder);
+  assert(ret->ctlpoints);
+
+  ret->next = NULL;
+
+  return ret;
+}
+
+bezierPatch* bezierPatchMake2(float umin, float vmin, float umax, float vmax, int uorder, int vorder, int dimension, int ustride, int vstride,  float* ctlpoints)
+{
+  bezierPatch* ret = (bezierPatch*) malloc(sizeof(bezierPatch));
+  assert(ret);
+  ret->umin = umin;
+  ret->vmin = vmin;
+  ret->umax = umax;
+  ret->vmax = vmax;
+  ret->uorder = uorder;
+  ret->vorder = vorder;
+  ret->dimension = dimension;
+  ret->ctlpoints = (float*) malloc(sizeof(float) * dimension * uorder * vorder);
+  assert(ret->ctlpoints);
+
+  /*copy the control points there*/
+  int the_ustride = vorder * dimension;
+  int the_vstride = dimension;
+  for(int i=0; i<uorder; i++)
+    for(int j=0; j<vorder; j++)
+      for(int k=0; k<dimension; k++)
+	ret->ctlpoints[i * the_ustride + j*the_vstride+k] = ctlpoints[i*ustride+j*vstride+k];
+  
+  ret->next = NULL;
+
+  return ret;
+}
+
+/*
+ *deallocate the space as allocated by Make
+ */
+void bezierPatchDelete(bezierPatch *b)
+{
+  free(b->ctlpoints);
+  free(b);
+}
+
+/*delete the whole linked list
+ */
+void bezierPatchDeleteList(bezierPatch *b)
+{
+  bezierPatch *temp;
+  while (b != NULL) {
+    temp = b;
+    b = b->next;
+    bezierPatchDelete(temp);    
+  }
+}
+
+bezierPatch* bezierPatchInsert(bezierPatch *list, bezierPatch *b)
+{
+  b->next = list;
+  return b;
+}
+
+/*print the data stored in this patch*/
+void bezierPatchPrint(bezierPatch *b)
+{
+  printf("bezierPatch:\n");
+  printf("umin,umax=(%f,%f), (vmin, vmax)=(%f,%f)\n", b->umin, b->umax, b->vmin, b->vmax);
+  printf("uorder=%i, vorder=%i\n", b->uorder, b->vorder);
+  printf("idmension = %i\n", b->dimension);
+}
+       
+/*print the whole list*/
+void bezierPatchPrintList(bezierPatch *list)
+{
+  bezierPatch* temp;
+  for(temp=list; temp != NULL; temp = temp->next)
+    bezierPatchPrint(temp);
+}
+
+void bezierPatchEval(bezierPatch *b, float u, float v, float ret[])
+{
+  if(   u >= b->umin && u<= b->umax
+     && v >= b->vmin && v<= b->vmax)
+    {
+
+      bezierSurfEval(b->umin, b->umax, b->uorder, b->vmin, b->vmax, b->vorder, b->dimension, b->ctlpoints, b->dimension * b->vorder, b->dimension, u, v, ret);
+
+     }
+  else if(b->next != NULL)
+    bezierPatchEval(b->next, u,v, ret);
+  else 
+    bezierSurfEval(b->umin, b->umax, b->uorder, b->vmin, b->vmax, b->vorder, b->dimension, b->ctlpoints, b->dimension * b->vorder, b->dimension, u, v, ret);    
+}
+
+/*the returned normal is normlized
+ */
+void bezierPatchEvalNormal(bezierPatch *b, float u, float v, float ret[])
+{
+  bezierSurfEvalNormal(b->umin, b->umax, b->uorder, b->vmin, b->vmax, b->vorder, b->dimension, b->ctlpoints, b->dimension * b->vorder, b->dimension, u, v, ret);  
+
+  if(   u >= b->umin && u<= b->umax
+     && v >= b->vmin && v<= b->vmax)
+    {
+      bezierSurfEvalNormal(b->umin, b->umax, b->uorder, b->vmin, b->vmax, b->vorder, b->dimension, b->ctlpoints, b->dimension * b->vorder, b->dimension, u, v, ret);
+     }
+  else if(b->next != NULL)
+    bezierPatchEvalNormal(b->next, u,v, ret);
+  else 
+    bezierSurfEvalNormal(b->umin, b->umax, b->uorder, b->vmin, b->vmax, b->vorder, b->dimension, b->ctlpoints, b->dimension * b->vorder, b->dimension, u, v, ret);    
+
+}
+
+void bezierPatchDraw(bezierPatch *bpatch, int u_reso, int v_reso)
+{
+  if(bpatch->dimension == 3)
+    glMap2f(GL_MAP2_VERTEX_3, bpatch->umin, bpatch->umax, 3*bpatch->vorder, bpatch->uorder, bpatch->vmin, bpatch->vmax,3, bpatch->vorder, (GLfloat*) bpatch->ctlpoints);
+  else
+    glMap2f(GL_MAP2_VERTEX_4, bpatch->umin, bpatch->umax, 4*bpatch->vorder, bpatch->uorder, bpatch->vmin, bpatch->vmax,3, bpatch->vorder, (GLfloat*) bpatch->ctlpoints);    
+
+  glMapGrid2f(u_reso, bpatch->umin, bpatch->umax, 
+	      v_reso, bpatch->vmin, bpatch->vmax);
+  glEvalMesh2(GL_LINE, 0, u_reso, 0, v_reso);
+}
+
+void bezierPatchListDraw(bezierPatch *list, int u_reso, int v_reso)
+{
+  bezierPatch *temp;
+glEnable(GL_LIGHTING);
+glEnable(GL_LIGHT0);
+glEnable(GL_MAP2_VERTEX_3);
+glEnable(GL_AUTO_NORMAL);
+glEnable(GL_NORMALIZE);
+glColor3f(1,0,0);
+#ifdef DEBUG
+printf("mapmap\n");
+#endif
+
+
+  for(temp = list; temp != NULL; temp = temp->next)
+    bezierPatchDraw(temp, u_reso, v_reso);
+}
+	    
+
+
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.h b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.h
new file mode 100644
index 000000000..981183afa
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatch.h
@@ -0,0 +1,104 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+/*
+*/
+
+#ifndef _BEZIERPATCH_H
+#define _BEZIERPATCH_H
+
+typedef struct bezierPatch{
+  float umin, vmin, umax, vmax;
+  int uorder; /*order= degree + 1*/
+  int vorder; 
+
+  /*
+   *the control points are stored in a one dimensional  array.
+   *the surface is defined as:
+   *      s(u,v) = sum_{i,j} P(i,j) * B_i(u) * B_j(v).
+   *where P(i,j) are the control points, B_i(.) are Bezier
+   *basis functions.
+   *Each control point can have dimension 3 or 4: (x,y,z,w).
+   *The components of P(i,j) are stored in a one dimensional 
+   *array: 
+   *       ctlpoints[]
+   *in the order of:
+   * P[0,0], P[0,1], ..., P[0,vorder-1],
+   * P[1,0], P[1,1], ..., P[1,vorder-1],
+   *  ...             
+   * P[uorder-1,0], P[uorder-1,1], ..., P[uorder-1,vorder-1].
+   */
+  int dimension;
+  float* ctlpoints;
+
+  /*
+   *in case we have to manage multiple bezierPatches.
+   */
+  struct bezierPatch  *next; 
+
+} bezierPatch;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+bezierPatch* bezierPatchMake(float umin, float vmin, float umax, float vmax, int urder, int vorder, int dimension);
+
+bezierPatch* bezierPatchMake2(float umin, float vmin, float umax, float vmax, int urder, int vorder, int dimension, int ustride, int vstride, float *ctlpoints);
+
+
+bezierPatch* bezierPatchInsert(bezierPatch *list, bezierPatch *b);
+
+void bezierPatchDelete(bezierPatch *b);
+
+void bezierPatchDeleteList(bezierPatch *b);
+
+void bezierPatchPrint(bezierPatch *b);
+
+void bezierPatchPrintList(bezierPatch *list);
+
+void bezierPatchEval(bezierPatch *b, float u, float v, float ret[]);
+
+void bezierPatchEvalNormal(bezierPatch *b, float u, float v, float retNormal[]);
+
+void bezierPatchEval(bezierPatch *b, float u, float v, float ret[]);
+
+void bezierPatchEvalNormal(bezierPatch *b, float u, float v, float ret[]);
+
+
+void bezierPatchDraw(bezierPatch *bpatch, int u_reso, int v_reso);
+
+void bezierPatchListDraw(bezierPatch *list, int u_reso, int v_reso);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.cc b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.cc
new file mode 100644
index 000000000..ac7ff84fc
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.cc
@@ -0,0 +1,610 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+** 
+** http://oss.sgi.com/projects/FreeB
+** 
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+** 
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+** 
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include "gluos.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+#include <GL/gl.h>
+#include "bezierEval.h"
+#include "bezierPatchMesh.h"
+
+static int isDegenerate(float A[2], float B[2], float C[2]);
+
+void drawStrips(float *vertex_array, float *normal_array, int *length_array, GLenum *type_array, int num_strips)
+{
+  int i,j,k;
+  k=0;
+  /*k is the index of the first component of the current vertex*/
+  for(i=0; i<num_strips; i++)
+    {
+      glBegin(type_array[i]);
+      for(j=0; j<length_array[i]; j++)
+	{
+	  glNormal3fv(normal_array+k);
+	  glVertex3fv(vertex_array+k);
+	  k += 3;
+	}
+      glEnd();
+    }
+}
+
+void bezierPatchMeshListDelDeg(bezierPatchMesh* list)
+{
+  bezierPatchMesh* temp;
+  for(temp=list; temp != NULL; temp = temp->next)
+    {
+      bezierPatchMeshDelDeg(temp);
+    }
+}
+
+void bezierPatchMeshListDelete(bezierPatchMesh *list)
+{
+  if(list == NULL) return;
+  bezierPatchMeshListDelete(list->next);
+  bezierPatchMeshDelete(list);  
+}
+
+
+
+
+bezierPatchMesh* bezierPatchMeshListReverse(bezierPatchMesh* list)
+{
+ bezierPatchMesh* ret=NULL;
+ bezierPatchMesh* temp;
+ bezierPatchMesh* nextone;
+  for(temp = list; temp != NULL; temp = nextone)
+    {
+      nextone = temp->next;
+      ret=bezierPatchMeshListInsert(ret, temp);
+    }
+ return ret;
+}
+
+/*maptype is either GL_MAP2_VERTEX_3 or GL_MAP2_VERTEX_4
+ */
+bezierPatchMesh *bezierPatchMeshMake(int maptype, float umin, float umax, int ustride, int uorder, float vmin, float vmax, int vstride, int vorder, float *ctlpoints,  int size_UVarray, int size_length_array)
+{
+  int i,j,k;
+  int dimension;
+  int the_ustride;
+  int the_vstride;
+
+  if(maptype == GL_MAP2_VERTEX_3) dimension = 3;
+  else if (maptype==GL_MAP2_VERTEX_4) dimension = 4;
+  else {
+    fprintf(stderr, "error in inMap2f, maptype=%i is wrong, maptype,map is invalid\n", maptype);
+    return NULL;
+  }
+
+  bezierPatchMesh *ret = (bezierPatchMesh*) malloc(sizeof(bezierPatchMesh));
+  assert(ret);
+
+  ret->bpatch_normal = NULL;
+  ret->bpatch_color  = NULL;
+  ret->bpatch_texcoord = NULL;
+  ret->bpatch = bezierPatchMake(umin, vmin, umax, vmax, uorder, vorder, dimension);
+
+  /*copy the control points there*/
+  the_ustride = vorder * dimension;
+  the_vstride = dimension;
+  for(i=0; i<uorder; i++)
+    for(j=0; j<vorder; j++)
+      for(k=0; k<dimension; k++)
+	ret->bpatch->ctlpoints[i * the_ustride + j*the_vstride+k] = ctlpoints[i*ustride+j*vstride+k];
+  
+
+  ret->size_UVarray = size_UVarray;
+  ret->size_length_array = size_length_array;
+  ret->UVarray = (float*) malloc(sizeof(float) * size_UVarray);
+  assert(ret->UVarray);
+  ret->length_array = (int *)malloc(sizeof(int) * size_length_array);
+  assert(ret->length_array);
+  ret->type_array = (GLenum *)malloc(sizeof(GLenum) * size_length_array);
+  assert(ret->type_array);
+
+  ret->index_UVarray = 0;
+  ret->index_length_array = 0;
+
+  ret->vertex_array = NULL;
+  ret->normal_array = NULL;
+  ret->color_array  = NULL;
+  ret->texcoord_array = NULL;
+
+  ret->next = NULL;
+  return ret;
+}
+
+bezierPatchMesh *bezierPatchMeshMake2(int size_UVarray, int size_length_array)
+{
+  bezierPatchMesh *ret = (bezierPatchMesh*) malloc(sizeof(bezierPatchMesh));
+  assert(ret);
+
+  ret->bpatch = NULL;
+  ret->bpatch_normal = NULL;
+  ret->bpatch_color  = NULL;
+  ret->bpatch_texcoord = NULL;
+
+  ret->size_UVarray = size_UVarray;
+  ret->size_length_array = size_length_array;
+  ret->UVarray = (float*) malloc(sizeof(float) * size_UVarray);
+  assert(ret->UVarray);
+  ret->length_array = (int *)malloc(sizeof(int) * size_length_array);
+  assert(ret->length_array);
+  ret->type_array = (GLenum *)malloc(sizeof(GLenum) * size_length_array);
+  assert(ret->type_array);
+
+  ret->index_UVarray = 0;
+  ret->index_length_array = 0;
+
+  ret->vertex_array = NULL;
+  ret->normal_array = NULL;
+  ret->color_array  = NULL;
+  ret->texcoord_array = NULL;
+
+  ret->next = NULL;
+  return ret;
+}
+
+void bezierPatchMeshPutPatch(bezierPatchMesh *bpm, int maptype, float umin, float umax, int ustride, int uorder, float vmin, float vmax, int vstride, int vorder, float *ctlpoints)
+{
+  switch(maptype){
+  case GL_MAP2_VERTEX_3:
+    bpm->bpatch = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 3, ustride, vstride, ctlpoints);
+    break;
+  case GL_MAP2_VERTEX_4:
+    bpm->bpatch = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 4,ustride, vstride, ctlpoints );
+    break;
+  case GL_MAP2_NORMAL:
+    bpm->bpatch_normal = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 3, ustride, vstride, ctlpoints);
+    break;
+  case GL_MAP2_INDEX:
+    bpm->bpatch_color = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 1, ustride, vstride, ctlpoints);
+    break;
+  case GL_MAP2_COLOR_4:
+    bpm->bpatch_color = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 4, ustride, vstride, ctlpoints);
+    break;
+  case GL_MAP2_TEXTURE_COORD_1:
+    bpm->bpatch_texcoord = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 1, ustride, vstride, ctlpoints);
+    break;
+  case GL_MAP2_TEXTURE_COORD_2:
+    bpm->bpatch_texcoord = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 2, ustride, vstride, ctlpoints);
+    break;    
+  case GL_MAP2_TEXTURE_COORD_3:
+    bpm->bpatch_texcoord = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 3, ustride, vstride, ctlpoints);
+    break;    
+  case GL_MAP2_TEXTURE_COORD_4:
+    bpm->bpatch_texcoord = bezierPatchMake2(umin, vmin, umax, vmax, uorder, vorder, 4, ustride, vstride, ctlpoints);
+    break;    
+  default:
+    fprintf(stderr, "error in bezierPatchMeshPutPatch, maptype=%i is wrong, maptype,map is invalid\n", maptype);
+  }
+}
+  
+
+/*delete everything including the arrays. So if you want to output the
+ *pointers of the arrays, you should not use this function to deallocate space.
+ *you should dealocate manually
+ */
+void bezierPatchMeshDelete(bezierPatchMesh *bpm)
+{
+  if(bpm->bpatch != NULL)
+    bezierPatchDelete(bpm->bpatch);
+  if(bpm->bpatch_normal != NULL)
+    bezierPatchDelete(bpm->bpatch_normal);
+  if(bpm->bpatch_color != NULL)
+    bezierPatchDelete(bpm->bpatch_color);
+  if(bpm->bpatch_texcoord != NULL)
+    bezierPatchDelete(bpm->bpatch_texcoord);
+  
+  free(bpm->UVarray);
+  free(bpm->length_array);
+  free(bpm->vertex_array);
+  free(bpm->normal_array);
+  free(bpm->type_array);
+  free(bpm);
+}
+ 
+/*begin a strip
+ *type is the primitive type:
+ */
+void bezierPatchMeshBeginStrip(bezierPatchMesh *bpm, GLenum type)
+{
+  bpm->counter = 0;
+  bpm->type = type;
+}
+
+/*signal the end of the current strip*/
+void bezierPatchMeshEndStrip(bezierPatchMesh *bpm)
+{
+  int i;
+  
+  /*if there are no vertices in this strip, then nothing needs to be done*/
+  if(bpm->counter == 0) return;
+  
+  /*if the length_array is full, it should be expanded*/
+  if(bpm->index_length_array >= bpm->size_length_array)
+    {
+      int *temp = (int*) malloc(sizeof(int) * (bpm->size_length_array*2 + 1));
+      assert(temp);
+      GLenum *temp_type = (GLenum*) malloc(sizeof(GLenum) * (bpm->size_length_array*2 + 1));
+      assert(temp_type);
+      /*update the size*/
+      bpm->size_length_array = bpm->size_length_array*2 + 1;
+      
+      /*copy*/
+      for(i=0; i<bpm->index_length_array; i++)
+	{
+	  temp[i] = bpm->length_array[i];
+	  temp_type[i] = bpm->type_array[i];
+	}
+      
+      /*deallocate old array*/
+      free(bpm->length_array);
+      free(bpm->type_array);
+      
+      /*point to the new array which is twice as bigger*/
+      bpm->length_array = temp;
+      bpm->type_array = temp_type;
+    }
+  bpm->type_array[bpm->index_length_array] = bpm->type;
+  bpm->length_array[bpm->index_length_array++] = bpm->counter;
+
+}
+
+/*insert (u,v) */
+void bezierPatchMeshInsertUV(bezierPatchMesh *bpm, float u, float v)
+{
+  int i;
+  /*if the UVarray is full, it should be expanded*/
+  if(bpm->index_UVarray+1 >= bpm->size_UVarray)
+    {
+      float *temp = (float*) malloc(sizeof(float) * (bpm->size_UVarray * 2 + 2));
+      assert(temp);
+      
+      /*update the size*/
+      bpm->size_UVarray = bpm->size_UVarray*2 + 2;
+      
+      /*copy*/
+      for(i=0; i<bpm->index_UVarray; i++)
+	{
+	  temp[i] = bpm->UVarray[i];
+	}
+      
+      /*deallocate old array*/
+      free(bpm->UVarray);
+      
+      /*pointing to the new arrays*/
+      bpm->UVarray = temp;
+    }
+  /*insert the new UV*/
+  bpm->UVarray[bpm->index_UVarray] = u;
+  bpm->index_UVarray++;
+  bpm->UVarray[bpm->index_UVarray] = v;
+  bpm->index_UVarray++;
+
+  /*update counter: one more vertex*/
+  bpm->counter++;
+
+
+}
+
+void bezierPatchMeshPrint(bezierPatchMesh *bpm)
+{
+  int i;
+  printf("the bezier patch is\n");
+  bezierPatchPrint(bpm->bpatch);
+  printf("index_length_array= %i\n", bpm->index_length_array);
+  printf("size_length_array =%i\n", bpm->size_length_array);
+  printf("index_UVarray =%i\n", bpm->index_UVarray);
+  printf("size_UVarray =%i\n", bpm->size_UVarray);
+  printf("UVarray is\n");
+  for(i=0; i<bpm->index_UVarray; i++)
+    printf("%f ", bpm->UVarray[i]);
+
+  printf("length_array is\n");
+  for(i=0; i<bpm->index_length_array; i++)
+    printf("%i ", bpm->length_array[i]);
+  printf("\n");
+
+}
+
+/*insert a new patch in front of the current linked list and return the new list*/
+bezierPatchMesh* bezierPatchMeshListInsert(bezierPatchMesh* list, bezierPatchMesh* bpm)
+{
+  bpm->next=list;
+  return bpm;
+}
+
+/*print all the patches*/
+void bezierPatchMeshListPrint(bezierPatchMesh* list)
+{
+  bezierPatchMesh *temp;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      bezierPatchMeshPrint(temp);
+    }
+}
+
+int bezierPatchMeshListTotalStrips(bezierPatchMesh* list)
+{
+  int sum=0;
+  bezierPatchMesh *temp;
+  for(temp=list; temp != NULL; temp = temp->next)
+    {
+      sum += temp->index_length_array;
+    }
+  return sum;
+}
+
+int bezierPatchMeshListTotalVert(bezierPatchMesh* list)
+{
+  int sum=0;
+  bezierPatchMesh *temp;
+  for(temp=list; temp != NULL; temp = temp->next)
+    {
+      sum += temp->index_UVarray;
+    }
+  return sum/2;
+}
+
+int bezierPatchMeshListNumTriangles(bezierPatchMesh* list)
+{
+  int sum=0;
+  bezierPatchMesh* temp;
+  for(temp=list; temp != NULL; temp = temp->next)
+    {
+      sum +=  bezierPatchMeshNumTriangles(temp);
+    }
+  return sum;
+}
+
+int bezierPatchMeshNumTriangles(bezierPatchMesh* bpm)
+{
+  int i;
+  int sum=0;
+  for(i=0; i<bpm->index_length_array; i++)
+    {
+      switch(bpm->type_array[i])
+	{
+	case GL_TRIANGLES:
+	  sum += bpm->length_array[i]/3;
+	  break;
+	case GL_TRIANGLE_FAN:
+	  if(bpm->length_array[i] > 2)
+	    sum += bpm->length_array[i]-2;
+	  break;
+	case GL_TRIANGLE_STRIP:
+	  if(bpm->length_array[i] > 2)
+	    sum += bpm->length_array[i]-2;
+	  break;
+	case GL_QUAD_STRIP:
+	  if(bpm->length_array[i]>2)
+	    sum += (bpm->length_array[i]-2);
+	  break;
+	default:
+	  fprintf(stderr,"error in bezierPatchMeshListNumTriangles, type invalid\n");
+	}
+    }
+  return sum;
+}
+
+/*delete degenerate triangles*/
+void bezierPatchMeshDelDeg(bezierPatchMesh* bpm)
+{
+  if(bpm == NULL) return;
+  int i,j,k;
+  int *new_length_array;
+  GLenum *new_type_array;
+  int index_new_length_array;
+  float *new_UVarray;
+  int index_new_UVarray;
+
+  new_length_array = (int*)malloc(sizeof(int) * bpm->index_length_array);
+  assert(new_length_array);
+  new_type_array = (GLenum*)malloc(sizeof(GLenum) * bpm->index_length_array);
+  assert(new_length_array);
+  new_UVarray = (float*) malloc(sizeof(float) * bpm->index_UVarray);
+  assert(new_UVarray);
+
+  index_new_length_array = 0;
+  index_new_UVarray=0;
+  k=0;
+  for(i=0; i<bpm->index_length_array; i++){
+    
+    /*(if not degenerate, we have to copy*/
+    if( (bpm->length_array[i] != 3) || (!isDegenerate(bpm->UVarray+k, bpm->UVarray+k+2, bpm->UVarray+k+4)))
+	  {
+	    for(j=0; j<2* bpm->length_array[i]; j++)
+	      new_UVarray[index_new_UVarray++] = bpm->UVarray[k++];
+
+	    new_length_array[index_new_length_array] = bpm->length_array[i];
+	    new_type_array[index_new_length_array] = bpm->type_array[i];
+	    index_new_length_array++;
+	  }
+    else
+      {
+	k += 6;
+      }
+  }  
+  free(bpm->UVarray);
+  free(bpm->length_array);
+  free(bpm->type_array);
+  bpm->UVarray=new_UVarray;
+  bpm->length_array=new_length_array;
+  bpm->type_array=new_type_array;
+  bpm->index_UVarray = index_new_UVarray;
+  bpm->index_length_array = index_new_length_array;
+  
+}
+
+/*(u,v) to XYZ
+ *the xyz and normals are stored in vertex_array, 
+ *and normal_array. the spaces of both are allocated here
+ */
+void bezierPatchMeshEval(bezierPatchMesh* bpm)
+{
+  int i,j,k,l;
+  float u,v;
+  float u0 = bpm->bpatch->umin;
+  float u1 = bpm->bpatch->umax;
+  int uorder = bpm->bpatch->uorder;
+  float v0 = bpm->bpatch->vmin;
+  float v1 = bpm->bpatch->vmax;
+  int vorder = bpm->bpatch->vorder;
+  int dimension = bpm->bpatch->dimension;
+  int ustride = dimension * vorder;
+  int vstride = dimension;
+  float *ctlpoints = bpm->bpatch->ctlpoints;
+  
+  bpm->vertex_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3);
+  assert(bpm->vertex_array);
+  bpm->normal_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3);
+  assert(bpm->normal_array);
+
+  k=0;
+  l=0;
+  for(i=0; i<bpm->index_length_array; i++)
+    {
+      for(j=0; j<bpm->length_array[i]; j++)
+	{
+	  u = bpm->UVarray[k];
+	  v = bpm->UVarray[k+1];
+	  bezierSurfEval(u0,u1,uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u,v, bpm->vertex_array+l);
+	  bezierSurfEvalNormal(u0,u1,uorder, v0, v1, vorder, dimension, ctlpoints, ustride, vstride, u,v, bpm->normal_array+l);
+	  k += 2;
+	  l += 3;
+	}
+    }
+}
+    
+void bezierPatchMeshListEval(bezierPatchMesh* list)
+{
+  bezierPatchMesh* temp;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      bezierPatchMeshEval(temp);
+    }
+}
+
+void bezierPatchMeshDraw(bezierPatchMesh* bpm)
+{
+  int i,j,k;
+  k=0;
+  /*k is the index of the first component of the current vertex*/
+  for(i=0; i<bpm->index_length_array; i++)
+    {
+      glBegin(bpm->type_array[i]);
+      for(j=0; j<bpm->length_array[i]; j++)
+	{
+	  glNormal3fv(bpm->normal_array+k);
+	  glVertex3fv(bpm->vertex_array+k);
+	  k+= 3;
+	}
+      glEnd();
+    }
+}
+
+void bezierPatchMeshListDraw(bezierPatchMesh* list)
+{
+  bezierPatchMesh* temp;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      bezierPatchMeshDraw(temp);
+    }
+}
+
+void bezierPatchMeshListCollect(bezierPatchMesh* list, float **vertex_array, float **normal_array, int **length_array, GLenum **type_array, int *num_strips)
+{
+  int i,j,k,l;
+  bezierPatchMesh *temp;
+  int total_num_vertices = bezierPatchMeshListTotalVert(list);
+  (*vertex_array) = (float *) malloc(sizeof(float) * total_num_vertices*3);
+  assert(*vertex_array);
+  (*normal_array) = (float *) malloc(sizeof(float) * total_num_vertices*3);
+  assert(*normal_array);
+
+  *num_strips = bezierPatchMeshListTotalStrips(list);
+   
+  *length_array = (int*) malloc(sizeof(int) * (*num_strips));
+  assert(*length_array);
+
+  *type_array = (GLenum*) malloc(sizeof(GLenum) * (*num_strips));
+  assert(*type_array);
+  
+  k=0;
+  l=0;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      int x=0;
+      for(i=0; i<temp->index_length_array; i++)
+	{
+	  for(j=0; j<temp->length_array[i]; j++)
+	    {
+	      (*vertex_array)[k] = temp->vertex_array[x];
+	      (*vertex_array)[k+1] = temp->vertex_array[x+1];
+	      (*vertex_array)[k+2] = temp->vertex_array[x+2];
+
+	      (*normal_array)[k] = temp->normal_array[x];
+	      (*normal_array)[k+1] = temp->normal_array[x+1];
+	      (*normal_array)[k+2] = temp->normal_array[x+2];
+
+	      x += 3;
+	      k += 3;
+	    }
+	  (*type_array)[l]  = temp->type_array[i];
+	  (*length_array)[l++] = temp->length_array[i];
+	}
+    }
+}
+
+
+
+static int isDegenerate(float A[2], float B[2], float C[2])
+{
+  if( (A[0] == B[0] && A[1]==B[1]) ||
+      (A[0] == C[0] && A[1]==C[1]) ||
+      (B[0] == C[0] && B[1]==C[1])
+     )
+    return 1;
+  else
+    return 0;
+}
+
+
+
+
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.h b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.h
new file mode 100644
index 000000000..449329665
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/bezierPatchMesh.h
@@ -0,0 +1,120 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+/*
+*/
+
+#ifndef _BEZIERPATCHMESH_H
+#define _BEZIERPATCHMESH_H
+
+#include "bezierPatch.h"
+
+typedef struct bezierPatchMesh{
+  bezierPatch *bpatch; /*vertex*/
+  bezierPatch *bpatch_normal; 
+  bezierPatch *bpatch_texcoord; /*s,t,r,q*/
+  bezierPatch *bpatch_color; /*RGBA*/
+
+  float *UVarray; /*all UV components of all vertices of all strips*/
+  int   *length_array; /*[i] is the number of vertices in the ith strip*/
+  GLenum *type_array;  /*[i] is the type of the ith primitive*/
+
+  /*to support dynamic insertion*/
+  int size_UVarray;
+  int index_UVarray;
+  int size_length_array;
+  int index_length_array;
+  
+  int counter; /*track the current strip size*/
+  GLenum type; /*track the current type: 0: GL_TRIANGLES, 1: GL_TRIANGLE_STRIP*/
+  
+  /*we eventually want to evaluate from (u,v) to (x,y,z) and draw them*/
+  float *vertex_array; /*each vertex contains three components*/
+  float *normal_array; /*each normal contains three components*/
+  float *color_array;
+  float *texcoord_array;
+
+  /*in case we need a linked list*/
+  struct bezierPatchMesh *next;
+} bezierPatchMesh;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+
+bezierPatchMesh *bezierPatchMeshMake(int maptype, float umin, float umax, int ustride, int uorder, float vmin, float vmax, int vstride, int vorder, float *ctlpoints,  int size_UVarray, int size_length_array);
+
+/*initilize patches to be null*/
+bezierPatchMesh *bezierPatchMeshMake2(int size_UVarray, int size_length_array);
+
+void bezierPatchMeshPutPatch(bezierPatchMesh *bpm, int maptype, float umin, float umax, int ustride, int uorder, float vmin, float vmax, int vstride, int vorder, float *ctlpoints);
+
+void bezierPatchMeshDelete(bezierPatchMesh *bpm);
+
+void bezierPatchMeshBeginStrip(bezierPatchMesh *bpm, GLenum type);
+
+void bezierPatchMeshEndStrip(bezierPatchMesh *bpm);
+
+void bezierPatchMeshInsertUV(bezierPatchMesh *bpm, float u, float v);
+
+void bezierPatchMeshPrint(bezierPatchMesh *bpm);
+
+bezierPatchMesh* bezierPatchMeshListInsert(bezierPatchMesh* list, bezierPatchMesh* bpm);
+
+void bezierPatchMeshListPrint(bezierPatchMesh* list);
+
+int bezierPatchMeshListTotalStrips(bezierPatchMesh* list);
+
+int bezierPatchMeshListTotalVert(bezierPatchMesh* list);
+int bezierPatchMeshNumTriangles(bezierPatchMesh* bpm);
+int bezierPatchMeshListNumTriangles(bezierPatchMesh* list);
+
+void bezierPatchMeshDelDeg(bezierPatchMesh* bpm);
+
+
+void bezierPatchMeshEval(bezierPatchMesh* bpm);
+  
+void bezierPatchMeshDraw(bezierPatchMesh* bpm);
+
+void bezierPatchMeshListDraw(bezierPatchMesh* list);
+void bezierPatchMeshListEval(bezierPatchMesh* list);
+void bezierPatchMeshListCollect(bezierPatchMesh* list, float **vertex_array, float **normal_array, int **length_array, GLenum **type_array, int *num_strips);
+
+void bezierPatchMeshListDelDeg(bezierPatchMesh* list);
+void bezierPatchMeshListDelete(bezierPatchMesh *list);
+bezierPatchMesh* bezierPatchMeshListReverse(bezierPatchMesh* list);
+void drawStrips(float *vertex_array, float *normal_array, int *length_array, GLenum *type_array, int num_strips);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.cc b/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.cc
new file mode 100644
index 000000000..b6591dba0
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.cc
@@ -0,0 +1,402 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+**
+** http://oss.sgi.com/projects/FreeB
+**
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+**
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+**
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+*/
+
+/*
+ * glcurveval.c++
+ *
+ */
+
+/* Polynomial Evaluator Interface */
+
+#include "gluos.h"
+#include "glimports.h"
+#include "glrenderer.h"
+#include "glcurveval.h"
+#include "nurbsconsts.h"
+ 
+OpenGLCurveEvaluator::OpenGLCurveEvaluator(void) 
+{
+  //no default callback functions
+  beginCallBackN = NULL;
+  endCallBackN = NULL;
+  vertexCallBackN = NULL;
+  normalCallBackN = NULL;
+  colorCallBackN = NULL;
+  texcoordCallBackN = NULL;
+  beginCallBackData = NULL;
+  endCallBackData = NULL;
+  vertexCallBackData = NULL;
+  normalCallBackData = NULL;
+  colorCallBackData = NULL;
+  texcoordCallBackData = NULL;
+
+  userData = NULL;
+  
+  vertex_flag = 0;
+  normal_flag = 0;
+  color_flag = 0;
+  texcoord_flag = 0;
+
+  em_vertex.uprime = -1.0;
+  em_normal.uprime = -1.0;
+  em_color.uprime = -1.0;
+  em_texcoord.uprime = -1.0;
+  output_triangles = 0; // don't output triangles by default
+}
+
+OpenGLCurveEvaluator::~OpenGLCurveEvaluator(void) 
+{ 
+}
+
+/* added nonsense to avoid the warning messages at compile time */
+void
+OpenGLCurveEvaluator::addMap(CurveMap *m)
+{
+	m = m;
+}
+
+void
+OpenGLCurveEvaluator::range1f(long type, REAL *from, REAL *to)
+{
+	type = type;
+	from = from;
+	to = to;
+}
+
+void
+OpenGLCurveEvaluator::domain1f(REAL ulo, REAL uhi)
+{
+	ulo = ulo;
+	uhi = uhi;
+}
+
+void
+OpenGLCurveEvaluator::bgnline(void)
+{
+  if(output_triangles)
+    beginCallBack(GL_LINE_STRIP, userData);
+  else
+    glBegin((GLenum) GL_LINE_STRIP);
+}
+
+void
+OpenGLCurveEvaluator::endline(void)
+{
+  if(output_triangles)
+    endCallBack(userData);
+  else
+    glEnd();
+}
+
+/*---------------------------------------------------------------------------
+ * disable - turn off a curve map
+ *---------------------------------------------------------------------------
+ */
+void
+OpenGLCurveEvaluator::disable(long type)
+{
+    glDisable((GLenum) type);
+}
+
+/*---------------------------------------------------------------------------
+ * enable - turn on a curve map
+ *---------------------------------------------------------------------------
+ */
+void
+OpenGLCurveEvaluator::enable(long type)
+{
+    glEnable((GLenum) type);
+}
+
+/*-------------------------------------------------------------------------
+ * mapgrid1f - define a lattice of points with origin and offset
+ *-------------------------------------------------------------------------
+ */
+void 
+OpenGLCurveEvaluator::mapgrid1f(long nu, REAL u0, REAL u1)
+{
+  if(output_triangles)
+    {
+      global_grid_u0 = u0;
+      global_grid_u1 = u1;
+      global_grid_nu = (int) nu;
+    }
+  else
+    glMapGrid1f((GLint) nu, (GLfloat) u0, (GLfloat) u1);
+}
+
+/*-------------------------------------------------------------------------
+ * bgnmap1 - preamble to curve definition and evaluations
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLCurveEvaluator::bgnmap1f(long)
+{
+  if(output_triangles)
+    {
+      //initialized so that no maps are set initially
+      vertex_flag = 0;
+      normal_flag = 0;
+      color_flag = 0;
+      texcoord_flag = 0;
+      //no need to worry about gl states when doing callback
+    }
+  else
+    glPushAttrib((GLbitfield) GL_EVAL_BIT);
+}
+
+/*-------------------------------------------------------------------------
+ * endmap1 - postamble to a curve map
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLCurveEvaluator::endmap1f(void)
+{
+  if(output_triangles)
+    {
+      
+    }
+  else
+    glPopAttrib();
+}
+
+/*-------------------------------------------------------------------------
+ * map1f - pass a desription of a curve map
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLCurveEvaluator::map1f(
+    long type,		 	/* map type */
+    REAL ulo,			/* lower parametric bound */
+    REAL uhi,			/* upper parametric bound */
+    long stride, 		/* distance to next point in REALS */
+    long order,			/* parametric order */
+    REAL *pts 			/* control points */
+)
+{
+  if(output_triangles)
+    {
+      int dimension = 0;
+      int which = 0;
+      switch(type){
+      case GL_MAP1_VERTEX_3:
+	which = 0;
+	dimension = 3;
+	break;
+      case GL_MAP1_VERTEX_4:
+	which=0;
+	dimension = 4;
+	break;
+      case GL_MAP1_INDEX:
+	which=2;
+	dimension = 1;
+	break;
+      case GL_MAP1_COLOR_4:
+	which=2;
+	dimension = 4;
+	break;
+      case GL_MAP1_NORMAL:
+	which=1;
+	dimension = 3;
+	break;
+      case GL_MAP1_TEXTURE_COORD_1:
+	which=3;
+	dimension = 1;
+	break;
+      case GL_MAP1_TEXTURE_COORD_2:
+	which=3;
+	dimension = 2;
+	break;
+	
+      case GL_MAP1_TEXTURE_COORD_3:
+	which=3;
+	dimension = 3;
+	break;
+      case GL_MAP1_TEXTURE_COORD_4:
+	which=3;
+	dimension = 4;
+	break;	
+      }
+      inMap1f(which, dimension, ulo, uhi, stride, order, pts); 	      
+    }       
+  else
+    glMap1f((GLenum) type, (GLfloat) ulo, (GLfloat) uhi, (GLint) stride, 
+	    (GLint) order, (const GLfloat *) pts);
+}
+
+/*-------------------------------------------------------------------------
+ * mapmesh1f - evaluate a mesh of points on lattice
+ *-------------------------------------------------------------------------
+ */
+void OpenGLCurveEvaluator::mapmesh1f(long style, long from, long to)
+{
+  if(output_triangles)
+    {
+      inMapMesh1f((int) from, (int) to);      
+    }
+  else
+    {
+      switch(style) {
+      default:
+      case N_MESHFILL:
+      case N_MESHLINE:
+	glEvalMesh1((GLenum) GL_LINE, (GLint) from, (GLint) to);
+	break;
+      case N_MESHPOINT:
+	glEvalMesh1((GLenum) GL_POINT, (GLint) from, (GLint) to);
+	break;
+      }
+    }
+}
+
+/*-------------------------------------------------------------------------
+ * evalpoint1i - evaluate a point on a curve
+ *-------------------------------------------------------------------------
+ */
+void OpenGLCurveEvaluator::evalpoint1i(long i)
+{
+    glEvalPoint1((GLint) i);
+}
+
+/*-------------------------------------------------------------------------
+ * evalcoord1f - evaluate a point on a curve
+ *-------------------------------------------------------------------------
+ */
+void OpenGLCurveEvaluator::evalcoord1f(long, REAL u)
+{
+    glEvalCoord1f((GLfloat) u);
+}
+
+void
+#ifdef _WIN32
+OpenGLCurveEvaluator::putCallBack(GLenum which, void (GLAPIENTRY *fn)())
+#else
+OpenGLCurveEvaluator::putCallBack(GLenum which, _GLUfuncptr fn)
+#endif
+{
+  switch(which)
+  {
+    case GLU_NURBS_BEGIN:
+      beginCallBackN = (void (GLAPIENTRY *) (GLenum)) fn;
+      break;
+    case GLU_NURBS_END:
+      endCallBackN = (void (GLAPIENTRY *) (void)) fn;
+      break;
+    case GLU_NURBS_VERTEX:
+      vertexCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_NORMAL:
+      normalCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_COLOR:
+      colorCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_TEXTURE_COORD:
+      texcoordCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_BEGIN_DATA:
+      beginCallBackData = (void (GLAPIENTRY *) (GLenum, void*)) fn;
+      break;
+    case GLU_NURBS_END_DATA:
+      endCallBackData = (void (GLAPIENTRY *) (void*)) fn;
+      break;
+    case GLU_NURBS_VERTEX_DATA:
+      vertexCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_NORMAL_DATA:
+      normalCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_COLOR_DATA:
+      colorCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_TEXTURE_COORD_DATA:
+      texcoordCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+  }
+}
+
+void
+OpenGLCurveEvaluator::beginCallBack(GLenum which, void *data)
+{
+  if(beginCallBackData)
+    beginCallBackData(which, data);
+  else if(beginCallBackN)
+    beginCallBackN(which);
+}
+
+void
+OpenGLCurveEvaluator::endCallBack(void *data)
+{
+  if(endCallBackData)
+    endCallBackData(data);
+  else if(endCallBackN)
+    endCallBackN();
+}
+
+void
+OpenGLCurveEvaluator::vertexCallBack(const GLfloat *vert, void* data)
+{
+  if(vertexCallBackData)
+    vertexCallBackData(vert, data);
+  else if(vertexCallBackN)
+    vertexCallBackN(vert);
+}
+
+
+void
+OpenGLCurveEvaluator::normalCallBack(const GLfloat *normal, void* data)
+{
+  if(normalCallBackData)
+    normalCallBackData(normal, data);
+  else if(normalCallBackN)
+    normalCallBackN(normal);
+}
+
+void
+OpenGLCurveEvaluator::colorCallBack(const GLfloat *color, void* data)
+{
+  if(colorCallBackData)
+    colorCallBackData(color, data);
+  else if(colorCallBackN)
+    colorCallBackN(color);
+}
+
+void
+OpenGLCurveEvaluator::texcoordCallBack(const GLfloat *texcoord, void* data)
+{
+  if(texcoordCallBackData)
+    texcoordCallBackData(texcoord, data);
+  else if(texcoordCallBackN)
+    texcoordCallBackN(texcoord);
+}
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.h b/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.h
new file mode 100644
index 000000000..8fa493034
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glcurveval.h
@@ -0,0 +1,157 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * glcurveval.h
+ *
+ */
+
+#ifndef __gluglcurveval_h_
+#define __gluglcurveval_h_
+
+#include "gluos.h"
+#include <GL/gl.h>
+#include <GL/glu.h>
+#include "basiccrveval.h"
+
+class CurveMap;
+
+/*for internal evaluator callback stuff*/
+#ifndef IN_MAX_BEZIER_ORDER
+#define IN_MAX_BEZIER_ORDER 40 /*XXX should be bigger than machine order*/
+#endif
+			
+#ifndef IN_MAX_DIMENSION
+#define IN_MAX_DIMENSION 4 
+#endif
+
+typedef struct curveEvalMachine{
+  REAL uprime; //cached previously evaluated uprime
+  int k; //the dimension
+  REAL u1;
+  REAL u2;
+  int ustride;
+  int uorder;
+  REAL ctlpoints[IN_MAX_BEZIER_ORDER*IN_MAX_DIMENSION];
+  REAL ucoeff[IN_MAX_BEZIER_ORDER];//cache the polynomial values
+} curveEvalMachine;
+
+class OpenGLCurveEvaluator : public BasicCurveEvaluator  {  
+public:
+			OpenGLCurveEvaluator(void);
+			virtual ~OpenGLCurveEvaluator(void);
+    void		range1f(long, REAL *, REAL *);
+    void		domain1f(REAL, REAL);
+    void		addMap(CurveMap *);
+
+    void		enable(long);
+    void		disable(long);
+    void		bgnmap1f(long);
+    void		map1f(long, REAL, REAL, long, long, REAL *);
+    void		mapgrid1f(long, REAL, REAL);
+    void		mapmesh1f(long, long, long);
+    void		evalpoint1i(long);
+    void		evalcoord1f(long, REAL);
+    void		endmap1f(void);
+
+    void		bgnline(void);
+    void		endline(void);
+
+    void                put_vertices_call_back(int flag)
+      {
+	output_triangles = flag;
+      }
+#ifdef _WIN32
+    void               putCallBack(GLenum which, void (GLAPIENTRY *fn)() );
+#else
+    void               putCallBack(GLenum which, _GLUfuncptr fn );
+#endif
+    void               set_callback_userData(void *data)
+      {
+	userData = data;
+      }
+
+/*------------------begin for curveEvalMachine------------*/
+curveEvalMachine em_vertex;
+curveEvalMachine em_normal;
+curveEvalMachine em_color;
+curveEvalMachine em_texcoord;
+int vertex_flag; //whether there is a vertex map or not
+int normal_flag; //whether there is a normal map or not
+int color_flag; //whether there is a color map or not
+int texcoord_flag; //whether there is a texture map or not
+
+REAL global_grid_u0;
+REAL global_grid_u1;
+int global_grid_nu;
+
+void inMap1f(int which, //0: vert, 1: norm, 2: color, 3: tex
+	     int dimension,
+	     REAL ulower,
+	     REAL uupper,
+	     int ustride,
+	     int uorder,
+	     REAL *ctlpoints);
+
+void inPreEvaluate(int order, REAL vprime, REAL *coeff);
+void inDoDomain1(curveEvalMachine *em, REAL u, REAL *retPoint);
+void inDoEvalCoord1(REAL u);
+void inMapMesh1f(int umin, int umax);
+
+void     (GLAPIENTRY *beginCallBackN) (GLenum type);
+void     (GLAPIENTRY *endCallBackN)   (void);
+void     (GLAPIENTRY *vertexCallBackN) (const GLfloat *vert);
+void     (GLAPIENTRY *normalCallBackN) (const GLfloat *normal);
+void     (GLAPIENTRY *colorCallBackN) (const GLfloat *color);
+void     (GLAPIENTRY *texcoordCallBackN) (const GLfloat *texcoord);
+
+void     (GLAPIENTRY *beginCallBackData) (GLenum type, void* data);
+void     (GLAPIENTRY *endCallBackData)   (void* data);
+void     (GLAPIENTRY *vertexCallBackData) (const GLfloat *vert, void* data);
+void     (GLAPIENTRY *normalCallBackData) (const GLfloat *normal, void* data);
+void     (GLAPIENTRY *colorCallBackData) (const GLfloat *color, void* data);
+void     (GLAPIENTRY *texcoordCallBackData) (const GLfloat *texcoord, void* data);
+
+void* userData; //the opaque pointer for Data callback functions
+void  beginCallBack(GLenum type, void* data);
+void endCallBack(void* data);
+void vertexCallBack(const GLfloat *vert, void *data);
+void normalCallBack(const GLfloat *normal, void* data);
+void colorCallBack(const  GLfloat *color, void* data);
+void texcoordCallBack(const GLfloat *texcoord, void* data);
+
+
+/*------------------end   for curveEvalMachine------------*/
+
+private:
+    int output_triangles; //true 1; false 0
+};
+
+#endif /* __gluglcurveval_h_ */
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glimports.h b/mesalib/src/glu/sgi/libnurbs/interface/glimports.h
new file mode 100644
index 000000000..6e69feb40
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glimports.h
@@ -0,0 +1,42 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * glimports.h
+ *
+ */
+
+#ifndef __gluimports_h_
+#define __gluimports_h_
+
+#include "mystdlib.h"
+#include "mystdio.h"
+
+#endif /* __gluimports_h_ */
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glinterface.cc b/mesalib/src/glu/sgi/libnurbs/interface/glinterface.cc
new file mode 100644
index 000000000..ba64bcd2d
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glinterface.cc
@@ -0,0 +1,469 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+** 
+** http://oss.sgi.com/projects/FreeB
+** 
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+** 
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+** 
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include "gluos.h"
+#include <GL/gl.h>
+#include <GL/glu.h>
+#include <stdio.h>
+#include "glimports.h"
+#include "glrenderer.h"
+#include "nurbsconsts.h"
+
+//#define DOWN_LOAD_NURBS
+#ifdef DOWN_LOAD_NURBS
+
+#include "oglTrimNurbs.h"
+static int surfcount = 0;
+static oglTrimNurbs* otn = NULL;
+nurbSurf* tempNurb = NULL;
+oglTrimLoops* tempTrim = NULL;
+#endif
+
+
+//for LOD
+extern "C" {void glu_LOD_eval_list(GLUnurbs *nurb, int level);}
+
+void glu_LOD_eval_list(GLUnurbs *nurb, int level)
+{
+	nurb->LOD_eval_list(level);
+}
+
+GLUnurbs * GLAPIENTRY
+gluNewNurbsRenderer(void)
+{
+  GLUnurbs *t;
+  
+  t = new GLUnurbs();
+  return t;
+}
+
+void GLAPIENTRY
+gluDeleteNurbsRenderer(GLUnurbs *r)
+{
+    delete r;
+}
+
+extern "C"
+void GLAPIENTRY
+
+gluDeleteNurbsTessellatorEXT(GLUnurbsObj *r)
+{
+  delete r;
+}
+
+void GLAPIENTRY
+gluBeginSurface(GLUnurbs *r)
+{
+#ifdef DOWN_LOAD_NURBS
+surfcount++;
+tempTrim = OTL_make(10,10);
+#endif
+    r->bgnsurface(0); 
+}
+
+void GLAPIENTRY
+gluBeginCurve(GLUnurbs *r)
+{
+    r->bgncurve(0); 
+}
+
+void GLAPIENTRY
+gluEndCurve(GLUnurbs *r)
+{
+    r->endcurve(); 
+}
+
+void GLAPIENTRY
+gluEndSurface(GLUnurbs *r)
+{
+#ifdef DOWN_LOAD_NURBS
+if(surfcount == 1)
+  otn = OTN_make(1);
+OTN_insert(otn, tempNurb, tempTrim);
+if(surfcount  >= 1)
+{
+#ifdef DEBUG
+printf("write file\n");
+#endif
+OTN_write(otn, "out.otn");
+
+}
+#endif
+
+    r->endsurface(); 
+}
+
+void GLAPIENTRY
+gluBeginTrim(GLUnurbs *r)
+{
+#ifdef DOWN_LOAD_NURBS
+OTL_bgnTrim(tempTrim);
+#endif
+
+    r->bgntrim(); 
+}
+
+void GLAPIENTRY
+gluEndTrim(GLUnurbs *r)
+{
+#ifdef DOWN_LOAD_NURBS
+OTL_endTrim(tempTrim);
+#endif
+    r->endtrim(); 
+}
+
+void GLAPIENTRY
+gluPwlCurve(GLUnurbs *r, GLint count, INREAL array[], 
+		GLint stride, GLenum type)
+{
+#ifdef DOWN_LOAD_NURBS
+OTL_pwlCurve(tempTrim, count, array, stride, type);
+#endif
+
+    int realType;
+    switch(type) {
+      case GLU_MAP1_TRIM_2:
+	realType = N_P2D;
+	break;
+      case GLU_MAP1_TRIM_3:
+	realType = N_P2DR;
+	break;
+      default:
+	realType = type;
+	break;
+    }
+    r->pwlcurve(count, array, sizeof(INREAL) * stride, realType);
+}
+
+void GLAPIENTRY
+gluNurbsCurve(GLUnurbs *r, GLint nknots, INREAL knot[], GLint stride, 
+		  INREAL ctlarray[], GLint order, GLenum type)
+{
+#ifdef DOWN_LOAD_NURBS
+OTL_nurbsCurve(tempTrim, nknots, knot, stride, ctlarray, order, type);
+#endif
+
+    int realType;
+
+    switch(type) {
+      case GLU_MAP1_TRIM_2:
+	realType = N_P2D;
+	break;
+      case GLU_MAP1_TRIM_3:
+	realType = N_P2DR;
+	break;
+      default:
+	realType = type;
+	break;
+    }
+
+    r->nurbscurve(nknots, knot, sizeof(INREAL) * stride, ctlarray, order, 
+	    realType);
+}
+
+void GLAPIENTRY
+gluNurbsSurface(GLUnurbs *r, GLint sknot_count, GLfloat *sknot, 
+			    GLint tknot_count, GLfloat *tknot, 
+			    GLint s_stride, GLint t_stride, 
+			    GLfloat *ctlarray, GLint sorder, GLint torder, 
+			    GLenum type)
+{
+#ifdef DOWN_LOAD_NURBS
+  {
+    int dimension;
+    switch(type){
+    case GL_MAP2_VERTEX_3:
+      dimension = 3;
+      break;
+    case GL_MAP2_VERTEX_4:
+      dimension = 4;
+      break;
+    default:
+      fprintf(stderr, "error in glinterface.c++, type no implemented\n");
+      exit(1);
+    }
+tempNurb = nurbSurfMake(sknot_count, sknot,
+			tknot_count, tknot,
+			sorder, torder,
+			dimension,
+			ctlarray,
+			s_stride, t_stride);
+			
+  }
+#endif
+
+    r->nurbssurface(sknot_count, sknot, tknot_count, tknot, 
+	    sizeof(INREAL) * s_stride, sizeof(INREAL) * t_stride, 
+	    ctlarray, sorder, torder, type);
+}
+
+void GLAPIENTRY
+gluLoadSamplingMatrices(GLUnurbs *r, const GLfloat modelMatrix[16],
+			    const GLfloat projMatrix[16], 
+			    const GLint viewport[4])
+{
+    r->useGLMatrices(modelMatrix, projMatrix, viewport);
+}
+
+void GLAPIENTRY
+gluNurbsProperty(GLUnurbs *r, GLenum property, GLfloat value)
+{
+    GLfloat nurbsValue;
+    
+    switch (property) {
+      case GLU_AUTO_LOAD_MATRIX:      
+        r->setautoloadmode(value);
+	return;
+
+      case GLU_CULLING:
+	if (value != 0.0) {
+	    nurbsValue = N_CULLINGON;
+	} else {
+	    nurbsValue = N_NOCULLING;
+	}
+	r->setnurbsproperty(GL_MAP2_VERTEX_3, N_CULLING, nurbsValue);
+	r->setnurbsproperty(GL_MAP2_VERTEX_4, N_CULLING, nurbsValue);
+	r->setnurbsproperty(GL_MAP1_VERTEX_3, N_CULLING, nurbsValue);
+	r->setnurbsproperty(GL_MAP1_VERTEX_4, N_CULLING, nurbsValue);
+        return;
+
+      case GLU_SAMPLING_METHOD:
+	if (value == GLU_PATH_LENGTH) {
+	    nurbsValue = N_PATHLENGTH;
+	} else if (value == GLU_PARAMETRIC_ERROR) {
+	    nurbsValue = N_PARAMETRICDISTANCE;
+	} else if (value == GLU_DOMAIN_DISTANCE) {
+	    nurbsValue = N_DOMAINDISTANCE;
+            r->set_is_domain_distance_sampling(1); //optimzing untrimmed case
+
+	} else if (value == GLU_OBJECT_PARAMETRIC_ERROR) {
+	    nurbsValue = N_OBJECTSPACE_PARA;
+	    r->setautoloadmode( 0.0 ); 
+	    r->setSamplingMatrixIdentity();
+	} else if (value == GLU_OBJECT_PATH_LENGTH) {
+	    nurbsValue = N_OBJECTSPACE_PATH;
+	    r->setautoloadmode( 0.0 ); 
+	    r->setSamplingMatrixIdentity();
+	} else {
+            r->postError(GLU_INVALID_VALUE);
+            return;
+        }
+
+	r->setnurbsproperty(GL_MAP2_VERTEX_3, N_SAMPLINGMETHOD, nurbsValue);
+	r->setnurbsproperty(GL_MAP2_VERTEX_4, N_SAMPLINGMETHOD, nurbsValue);
+	r->setnurbsproperty(GL_MAP1_VERTEX_3, N_SAMPLINGMETHOD, nurbsValue);
+	r->setnurbsproperty(GL_MAP1_VERTEX_4, N_SAMPLINGMETHOD, nurbsValue);
+	return;
+
+      case GLU_SAMPLING_TOLERANCE:
+	r->setnurbsproperty(GL_MAP2_VERTEX_3, N_PIXEL_TOLERANCE, value);
+	r->setnurbsproperty(GL_MAP2_VERTEX_4, N_PIXEL_TOLERANCE, value);
+	r->setnurbsproperty(GL_MAP1_VERTEX_3, N_PIXEL_TOLERANCE, value);
+	r->setnurbsproperty(GL_MAP1_VERTEX_4, N_PIXEL_TOLERANCE, value);
+	return;
+
+      case GLU_PARAMETRIC_TOLERANCE:
+	r->setnurbsproperty(GL_MAP2_VERTEX_3, N_ERROR_TOLERANCE, value);
+        r->setnurbsproperty(GL_MAP2_VERTEX_4, N_ERROR_TOLERANCE, value);
+        r->setnurbsproperty(GL_MAP1_VERTEX_3, N_ERROR_TOLERANCE, value);
+        r->setnurbsproperty(GL_MAP1_VERTEX_4, N_ERROR_TOLERANCE, value);
+        return;
+	
+
+      case GLU_DISPLAY_MODE:
+	
+	if (value == GLU_FILL) {
+	  nurbsValue = N_FILL;
+	} else if (value == GLU_OUTLINE_POLYGON) {
+	  nurbsValue = N_OUTLINE_POLY;
+	} else if (value == GLU_OUTLINE_PATCH) {
+	  nurbsValue = N_OUTLINE_PATCH;
+	} else {
+	  r->postError(GLU_INVALID_VALUE);
+	  return;
+	}
+	r->setnurbsproperty(N_DISPLAY, nurbsValue);
+	
+	break;
+
+      case GLU_U_STEP:
+    	r->setnurbsproperty(GL_MAP1_VERTEX_3, N_S_STEPS, value);
+    	r->setnurbsproperty(GL_MAP1_VERTEX_4, N_S_STEPS, value);
+    	r->setnurbsproperty(GL_MAP2_VERTEX_3, N_S_STEPS, value);
+    	r->setnurbsproperty(GL_MAP2_VERTEX_4, N_S_STEPS, value);
+	
+	//added for optimizing untrimmed case
+        r->set_domain_distance_u_rate(value);
+	break;
+
+      case GLU_V_STEP:
+        r->setnurbsproperty(GL_MAP1_VERTEX_3, N_T_STEPS, value);
+        r->setnurbsproperty(GL_MAP1_VERTEX_4, N_T_STEPS, value);
+        r->setnurbsproperty(GL_MAP2_VERTEX_3, N_T_STEPS, value);
+        r->setnurbsproperty(GL_MAP2_VERTEX_4, N_T_STEPS, value);
+
+	//added for optimizing untrimmed case
+        r->set_domain_distance_v_rate(value);
+	break;
+	
+      case GLU_NURBS_MODE:
+	if(value == GLU_NURBS_RENDERER)
+	  r->put_callbackFlag(0);
+	else if(value == GLU_NURBS_TESSELLATOR)
+	  r->put_callbackFlag(1);
+	else
+	  r->postError(GLU_INVALID_ENUM);
+	break;
+
+      default:
+	r->postError(GLU_INVALID_ENUM);
+	return;	
+    }
+}
+
+void GLAPIENTRY
+gluGetNurbsProperty(GLUnurbs *r, GLenum property, GLfloat *value)
+{
+    GLfloat nurbsValue;
+
+    switch(property) {
+      case GLU_AUTO_LOAD_MATRIX:
+	if (r->getautoloadmode()) {
+	    *value = GL_TRUE;
+	} else {
+	    *value = GL_FALSE;
+	}
+	break;
+      case GLU_CULLING:
+	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_CULLING, &nurbsValue);
+	if (nurbsValue == N_CULLINGON) {
+	    *value = GL_TRUE;
+	} else {
+	    *value = GL_FALSE;
+	}
+	break;
+      case GLU_SAMPLING_METHOD:
+	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_SAMPLINGMETHOD, value);
+	if(*value == N_PATHLENGTH)
+	  *value = GLU_PATH_LENGTH;
+	else if(*value == N_PARAMETRICDISTANCE)
+	  *value = GLU_PARAMETRIC_ERROR;
+	else if(*value == N_DOMAINDISTANCE)
+	  *value = GLU_DOMAIN_DISTANCE;
+	else if(*value == N_OBJECTSPACE_PATH)
+	  *value = GLU_OBJECT_PATH_LENGTH;
+	else if(*value == N_OBJECTSPACE_PARA)
+	  *value = GLU_OBJECT_PARAMETRIC_ERROR;	
+	break;
+      case GLU_SAMPLING_TOLERANCE:
+	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_PIXEL_TOLERANCE, value);
+	break;
+      case GLU_PARAMETRIC_TOLERANCE:
+	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_ERROR_TOLERANCE, value);
+        break;
+
+      case GLU_U_STEP:
+    	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_S_STEPS, value);
+	break;
+      case GLU_V_STEP:
+    	r->getnurbsproperty(GL_MAP2_VERTEX_3, N_T_STEPS, value);
+	break;
+      case GLU_DISPLAY_MODE:
+	r->getnurbsproperty(N_DISPLAY, &nurbsValue);
+	if (nurbsValue == N_FILL) {
+	    *value = GLU_FILL;
+	} else if (nurbsValue == N_OUTLINE_POLY) {
+	    *value = GLU_OUTLINE_POLYGON;
+	} else {
+	    *value = GLU_OUTLINE_PATCH;
+	}
+	break;
+
+      case GLU_NURBS_MODE:
+	if(r->is_callback())
+	  *value = GLU_NURBS_TESSELLATOR;
+	else
+	  *value = GLU_NURBS_RENDERER;
+	break;
+	
+      default:
+	r->postError(GLU_INVALID_ENUM);
+	return;
+    }
+}
+
+extern "C" void GLAPIENTRY
+gluNurbsCallback(GLUnurbs *r, GLenum which, _GLUfuncptr fn )
+{
+    switch (which) {
+    case GLU_NURBS_BEGIN:
+    case GLU_NURBS_END:
+    case GLU_NURBS_VERTEX:
+    case GLU_NURBS_NORMAL:
+    case GLU_NURBS_TEXTURE_COORD:
+    case GLU_NURBS_COLOR:
+    case GLU_NURBS_BEGIN_DATA:
+    case GLU_NURBS_END_DATA:
+    case GLU_NURBS_VERTEX_DATA:
+    case GLU_NURBS_NORMAL_DATA:
+    case GLU_NURBS_TEXTURE_COORD_DATA:
+    case GLU_NURBS_COLOR_DATA: 
+	r->putSurfCallBack(which, fn);
+	break;
+
+    case GLU_NURBS_ERROR:
+	r->errorCallback = (void (APIENTRY *)( GLenum e )) fn;
+	break;
+    default:
+	r->postError(GLU_INVALID_ENUM);
+	return;
+    }
+}
+
+extern "C"
+void GLAPIENTRY
+gluNurbsCallbackDataEXT(GLUnurbs* r, void* userData)
+{
+  r->setNurbsCallbackData(userData);
+}
+
+extern "C"
+void GLAPIENTRY
+gluNurbsCallbackData(GLUnurbs* r, void* userData)
+{
+  gluNurbsCallbackDataEXT(r,userData);
+}
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.cc b/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.cc
new file mode 100644
index 000000000..17123fbc6
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.cc
@@ -0,0 +1,301 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+**
+** http://oss.sgi.com/projects/FreeB
+**
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+**
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+**
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include "gluos.h"
+#include "glimports.h"
+#include "glrenderer.h"
+
+GLUnurbs::GLUnurbs()
+	: NurbsTessellator(curveEvaluator, surfaceEvaluator)
+{
+    redefineMaps();
+    defineMap(GL_MAP2_NORMAL, 0, 3);
+    defineMap(GL_MAP1_NORMAL, 0, 3);
+    defineMap(GL_MAP2_TEXTURE_COORD_1, 0, 1);
+    defineMap(GL_MAP1_TEXTURE_COORD_1, 0, 1);
+    defineMap(GL_MAP2_TEXTURE_COORD_2, 0, 2);
+    defineMap(GL_MAP1_TEXTURE_COORD_2, 0, 2);
+    defineMap(GL_MAP2_TEXTURE_COORD_3, 0, 3);
+    defineMap(GL_MAP1_TEXTURE_COORD_3, 0, 3);
+    defineMap(GL_MAP2_TEXTURE_COORD_4, 1, 4);
+    defineMap(GL_MAP1_TEXTURE_COORD_4, 1, 4);
+    defineMap(GL_MAP2_VERTEX_4, 1, 4);
+    defineMap(GL_MAP1_VERTEX_4, 1, 4);
+    defineMap(GL_MAP2_VERTEX_3, 0, 3);
+    defineMap(GL_MAP1_VERTEX_3, 0, 3);
+    defineMap(GL_MAP2_COLOR_4, 0, 4);
+    defineMap(GL_MAP1_COLOR_4, 0, 4);
+    defineMap(GL_MAP2_INDEX, 0, 1);
+    defineMap(GL_MAP1_INDEX, 0, 1);
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_SAMPLINGMETHOD, (float) N_PATHLENGTH);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_SAMPLINGMETHOD, (float) N_PATHLENGTH);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_SAMPLINGMETHOD, (float) N_PATHLENGTH);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_SAMPLINGMETHOD, (float) N_PATHLENGTH);
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_PIXEL_TOLERANCE, (float) 50.0);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_PIXEL_TOLERANCE, (float) 50.0);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_PIXEL_TOLERANCE, (float) 50.0);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_PIXEL_TOLERANCE, (float) 50.0);
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_ERROR_TOLERANCE, (float) 0.50);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_ERROR_TOLERANCE, (float) 0.50);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_ERROR_TOLERANCE, (float) 0.50);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_ERROR_TOLERANCE, (float) 0.50);
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_S_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_S_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_S_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_S_STEPS, (float) 100.0);
+
+    //added for optimizing untrimmed case
+    set_domain_distance_u_rate(100.0);
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_T_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_T_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_T_STEPS, (float) 100.0);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_T_STEPS, (float) 100.0);
+
+    //added for optimizing untrimmed case
+    set_domain_distance_v_rate(100.0);
+    set_is_domain_distance_sampling(0); //since the default is path_length
+
+    //default autoloadmode is true
+    autoloadmode = 1;
+
+    //default callbackFlag is 0
+    callbackFlag = 0;
+
+    errorCallback = NULL;
+}
+
+void
+GLUnurbs::bgnrender(void)
+{
+    if (autoloadmode) {
+	loadGLMatrices();
+    }
+}
+
+void
+GLUnurbs::endrender(void)
+{
+}
+
+void
+GLUnurbs::errorHandler(int i)
+{
+    int gluError;
+
+    gluError = i + (GLU_NURBS_ERROR1 - 1);
+    postError( gluError );
+}
+
+void
+GLUnurbs::loadGLMatrices(void)
+{
+    GLfloat vmat[4][4];
+    GLint viewport[4];
+
+    grabGLMatrix((GLfloat (*)[4]) vmat);
+    loadCullingMatrix((GLfloat (*)[4]) vmat);
+    ::glGetIntegerv((GLenum) GL_VIEWPORT, (GLint *) viewport);
+    loadSamplingMatrix((const GLfloat (*)[4]) vmat, (const GLint *) viewport);
+}
+
+void
+GLUnurbs::useGLMatrices(const GLfloat modelMatrix[16],
+			  const GLfloat projMatrix[16],
+			  const GLint viewport[4])
+{
+    GLfloat vmat[4][4];
+
+    multmatrix4d(vmat, (const GLfloat (*)[4]) modelMatrix,
+	    (const GLfloat (*)[4]) projMatrix);
+    loadCullingMatrix((GLfloat (*)[4]) vmat);
+    loadSamplingMatrix((const GLfloat (*)[4]) vmat, (const GLint *) viewport);
+}
+
+/*--------------------------------------------------------------------------
+ * grabGLMatrix
+ *--------------------------------------------------------------------------
+ */
+
+void
+GLUnurbs::grabGLMatrix(GLfloat vmat[4][4])
+{
+    GLfloat m1[4][4], m2[4][4];
+
+    ::glGetFloatv((GLenum) GL_MODELVIEW_MATRIX, (GLfloat *) &(m1[0][0]));
+    ::glGetFloatv((GLenum) GL_PROJECTION_MATRIX, (GLfloat *) &(m2[0][0]));
+    multmatrix4d((GLfloat (*)[4]) vmat,
+	    (const GLfloat (*)[4]) m1, (const GLfloat (*)[4]) m2);
+}
+
+//for object space tesselation: view independent
+void
+GLUnurbs::setSamplingMatrixIdentity( void )
+{
+  INREAL smat[4][4] = {
+    {1,0,0,0},
+    {0,1,0,0},
+    {0,0,1,0},
+    {0,0,0,1}
+  };
+  const long rstride = sizeof(smat[0]) / sizeof(smat[0][0]);
+  const long cstride = 1;
+
+  setnurbsproperty(GL_MAP1_VERTEX_3, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+		   cstride);
+  setnurbsproperty(GL_MAP1_VERTEX_4, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+		   cstride);
+  setnurbsproperty(GL_MAP2_VERTEX_3, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+		   cstride);
+  setnurbsproperty(GL_MAP2_VERTEX_4, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+		   cstride);
+}
+
+
+void
+GLUnurbs::loadSamplingMatrix(const GLfloat vmat[4][4],
+			       const GLint viewport[4])
+{
+
+    /* rescale the mapping to correspond to pixels in x/y */
+    REAL xsize = 0.5 * (REAL) (viewport[2]);
+    REAL ysize = 0.5 * (REAL) (viewport[3]);
+
+    INREAL smat[4][4];
+    smat[0][0] = vmat[0][0] * xsize;
+    smat[1][0] = vmat[1][0] * xsize;
+    smat[2][0] = vmat[2][0] * xsize;
+    smat[3][0] = vmat[3][0] * xsize;
+
+    smat[0][1] = vmat[0][1] * ysize;
+    smat[1][1] = vmat[1][1] * ysize;
+    smat[2][1] = vmat[2][1] * ysize;
+    smat[3][1] = vmat[3][1] * ysize;
+
+    smat[0][2] = 0.0;
+    smat[1][2] = 0.0;
+    smat[2][2] = 0.0;
+    smat[3][2] = 0.0;
+
+    smat[0][3] = vmat[0][3];
+    smat[1][3] = vmat[1][3];
+    smat[2][3] = vmat[2][3];
+    smat[3][3] = vmat[3][3];
+
+    const long rstride = sizeof(smat[0]) / sizeof(smat[0][0]);
+    const long cstride = 1;
+
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+	    cstride);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+	    cstride);
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+	    cstride);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_SAMPLINGMATRIX, &smat[0][0], rstride,
+	    cstride);
+}
+
+void
+GLUnurbs::loadCullingMatrix(GLfloat vmat[4][4])
+{
+    INREAL cmat[4][4];
+
+    cmat[0][0] = vmat[0][0];
+    cmat[0][1] = vmat[0][1];
+    cmat[0][2] = vmat[0][2];
+    cmat[0][3] = vmat[0][3];
+
+    cmat[1][0] = vmat[1][0];
+    cmat[1][1] = vmat[1][1];
+    cmat[1][2] = vmat[1][2];
+    cmat[1][3] = vmat[1][3];
+
+    cmat[2][0] = vmat[2][0];
+    cmat[2][1] = vmat[2][1];
+    cmat[2][2] = vmat[2][2];
+    cmat[2][3] = vmat[2][3];
+
+    cmat[3][0] = vmat[3][0];
+    cmat[3][1] = vmat[3][1];
+    cmat[3][2] = vmat[3][2];
+    cmat[3][3] = vmat[3][3];
+
+    const long rstride = sizeof(cmat[0]) / sizeof(cmat[0][0]);
+    const long cstride = 1;
+
+    setnurbsproperty(GL_MAP2_VERTEX_3, N_CULLINGMATRIX, &cmat[0][0], rstride,
+	    cstride);
+    setnurbsproperty(GL_MAP2_VERTEX_4, N_CULLINGMATRIX, &cmat[0][0], rstride,
+	    cstride);
+	//added for curves by zl
+    setnurbsproperty(GL_MAP1_VERTEX_3, N_CULLINGMATRIX, &cmat[0][0], rstride,
+	    cstride);
+    setnurbsproperty(GL_MAP1_VERTEX_4, N_CULLINGMATRIX, &cmat[0][0], rstride,
+	    cstride);
+}
+
+/*---------------------------------------------------------------------
+ * A = B * MAT ; transform a 4d vector through a 4x4 matrix
+ *---------------------------------------------------------------------
+ */
+void
+GLUnurbs::transform4d(GLfloat A[4], GLfloat B[4], GLfloat mat[4][4])
+{
+
+    A[0] = B[0]*mat[0][0] + B[1]*mat[1][0] + B[2]*mat[2][0] + B[3]*mat[3][0];
+    A[1] = B[0]*mat[0][1] + B[1]*mat[1][1] + B[2]*mat[2][1] + B[3]*mat[3][1];
+    A[2] = B[0]*mat[0][2] + B[1]*mat[1][2] + B[2]*mat[2][2] + B[3]*mat[3][2];
+    A[3] = B[0]*mat[0][3] + B[1]*mat[1][3] + B[2]*mat[2][3] + B[3]*mat[3][3];
+}
+
+/*---------------------------------------------------------------------
+ * new = [left][right] ; multiply two matrices together
+ *---------------------------------------------------------------------
+ */
+void
+GLUnurbs::multmatrix4d (GLfloat n[4][4], const GLfloat left[4][4],
+		 const GLfloat right[4][4])
+{
+    transform4d ((GLfloat *) n[0],(GLfloat *) left[0],(GLfloat (*)[4]) right);
+    transform4d ((GLfloat *) n[1],(GLfloat *) left[1],(GLfloat (*)[4]) right);
+    transform4d ((GLfloat *) n[2],(GLfloat *) left[2],(GLfloat (*)[4]) right);
+    transform4d ((GLfloat *) n[3],(GLfloat *) left[3],(GLfloat (*)[4]) right);
+}
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.h b/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.h
new file mode 100644
index 000000000..3b72d44d1
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glrenderer.h
@@ -0,0 +1,146 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * glrenderer.h
+ *
+ */
+
+#ifndef __gluglrenderer_h_
+#define __gluglrenderer_h_
+
+#include <GL/gl.h>
+#include <GL/glu.h>
+#include "nurbstess.h"
+#include "glsurfeval.h"
+#include "glcurveval.h"
+
+extern "C" {
+      typedef void (APIENTRY *errorCallbackType)( GLenum );
+}
+
+class GLUnurbs : public NurbsTessellator {
+
+public:
+		GLUnurbs( void );
+    void 	loadGLMatrices( void );
+    void	useGLMatrices( const GLfloat modelMatrix[16], 
+			       const GLfloat projMatrix[16],
+			       const GLint viewport[4] );
+    void		setSamplingMatrixIdentity( void );
+
+    void 	errorHandler( int );
+    void	bgnrender( void );
+    void	endrender( void );
+    void	setautoloadmode( INREAL value )
+      		    { 
+		      
+		      if (value) autoloadmode = GL_TRUE; 
+		      else autoloadmode = GL_FALSE; 
+			
+		    }
+    GLboolean	getautoloadmode( void ) { return autoloadmode; }
+
+    errorCallbackType errorCallback;
+    void	postError( int which ) 
+		    { if (errorCallback) (errorCallback)( (GLenum)which ); }
+#ifdef _WIN32
+    void        putSurfCallBack(GLenum which, void (GLAPIENTRY *fn)() )
+#else
+    void        putSurfCallBack(GLenum which, _GLUfuncptr fn )
+#endif
+      {
+	curveEvaluator.putCallBack(which, fn);
+	surfaceEvaluator.putCallBack(which, fn);
+      }
+
+    int         get_vertices_call_back()
+      {
+	return surfaceEvaluator.get_vertices_call_back();
+      }
+    
+    void        put_vertices_call_back(int flag)
+      {
+	surfaceEvaluator.put_vertices_call_back(flag);
+      }
+
+    int         get_callback_auto_normal()
+      {
+        return surfaceEvaluator.get_callback_auto_normal();
+      }
+ 
+    void         put_callback_auto_normal(int flag)
+      {
+        surfaceEvaluator.put_callback_auto_normal(flag);
+      }
+
+    void       setNurbsCallbackData(void* userData)
+      {
+       curveEvaluator.set_callback_userData(userData);
+       surfaceEvaluator.set_callback_userData(userData);
+     }
+
+
+    //for LOD
+    void LOD_eval_list(int level)
+      {
+	surfaceEvaluator.LOD_eval_list(level);
+      }
+
+    //NEWCALLBACK
+    int        is_callback()
+      {
+	return callbackFlag;
+      }
+    void       put_callbackFlag(int flag)
+      {
+	callbackFlag = flag;
+	surfaceEvaluator.put_vertices_call_back(flag);
+	curveEvaluator.put_vertices_call_back(flag);
+      }
+
+private:
+    GLboolean			autoloadmode;
+    OpenGLSurfaceEvaluator	surfaceEvaluator;
+    OpenGLCurveEvaluator	curveEvaluator;
+
+    void		loadSamplingMatrix( const GLfloat vmat[4][4], 
+			        const GLint viewport[4] );
+    void		loadCullingMatrix( GLfloat vmat[4][4] );
+    static void		grabGLMatrix( GLfloat vmat[4][4] );
+    static void		transform4d( GLfloat A[4], GLfloat B[4], 
+				GLfloat mat[4][4] );
+    static void		multmatrix4d( GLfloat n[4][4], const GLfloat left[4][4],
+				const GLfloat right[4][4] );
+
+   int                  callbackFlag;
+};
+
+#endif /* __gluglrenderer_h_ */
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.cc b/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.cc
new file mode 100644
index 000000000..b5bfab1e2
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.cc
@@ -0,0 +1,1293 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+**
+** http://oss.sgi.com/projects/FreeB
+**
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+**
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+**
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+*/
+
+/*
+ * glsurfeval.c++
+ *
+ */
+
+/* Polynomial Evaluator Interface */
+#include "gluos.h"
+#include <stdio.h>
+#include "glimports.h"
+#include "glrenderer.h"
+#include "glsurfeval.h"
+#include "nurbsconsts.h"
+#include "bezierPatchMesh.h"
+
+
+//extern int surfcount;
+//int surfcount=0;
+
+/*#define USE_INTERNAL_EVAL*/ //use internal evaluator
+
+/*whether do evaluation or not*/
+/*#define NO_EVALUATION*/
+
+//#define USE_LOD //for LOD test, have to turn on USE_LOD in insurfeval.c++ too
+
+/*for statistics*/
+//#define STATISTICS
+#ifdef STATISTICS
+static int STAT_num_of_triangles=0;
+static int STAT_num_of_eval_vertices=0;
+static int STAT_num_of_quad_strips=0;
+#endif
+
+/*for output triangles*/
+/*#define OUTPUT_TRIANGLES*/
+
+
+/*#define FOR_CHRIS*/
+#ifdef FOR_CHRIS
+extern "C"  {  void                evalUStripExt(int n_upper, REAL v_upper, REAL* upper_val,
+				   int n_lower, REAL v_lower, REAL* lower_val);}
+
+extern "C" {   void                evalVStripExt(int n_left, REAL u_left, REAL* left_val,
+				   int n_right, REAL u_right, REAL* right_val);
+	     }
+#endif
+
+
+/**************begin for LOD_eval_list***********/
+void OpenGLSurfaceEvaluator::LOD_eval_list(int level)
+{
+  if(level == 0)
+    LOD_eval_level = 1;
+  else if(level == 1)
+    LOD_eval_level = 2;
+  else if(level == 2)
+    LOD_eval_level = 4;
+  else
+    LOD_eval_level = 8;
+
+  inBPMListEvalEM(global_bpm);
+}
+
+
+OpenGLSurfaceEvaluator::OpenGLSurfaceEvaluator()
+{
+    int i;
+
+    for (i=0; i<VERTEX_CACHE_SIZE; i++) {
+	vertexCache[i] = new StoredVertex;
+    }
+    tmeshing = 0;
+    which = 0;
+    vcount = 0;
+
+    global_uorder = 0;
+    global_vorder = 0;
+    global_uprime = -1.0;
+    global_vprime = -1.0;
+    global_vprime_BV = -1.0;
+    global_uprime_BU = -1.0;
+    global_uorder_BU = 0;
+    global_vorder_BU = 0;
+    global_uorder_BV = 0;
+    global_vorder_BV = 0;
+    global_baseData = NULL;
+        
+    global_bpm = NULL;
+    output_triangles = 0; //don't output triangles by default
+
+    //no default callback functions
+    beginCallBackN = NULL;
+    endCallBackN = NULL;
+    vertexCallBackN = NULL;
+    normalCallBackN = NULL;
+    colorCallBackN = NULL;
+    texcoordCallBackN = NULL;
+    beginCallBackData = NULL;
+    endCallBackData = NULL;
+    vertexCallBackData = NULL;
+    normalCallBackData = NULL;
+    colorCallBackData = NULL;
+    texcoordCallBackData = NULL;
+
+    userData = NULL;
+
+    auto_normal_flag = 0;
+    callback_auto_normal = 0; //default of GLU_CALLBACK_AUTO_NORMAL is 0
+    vertex_flag = 0;
+    normal_flag = 0;
+    color_flag = 0;
+    texcoord_flag = 0;
+
+    em_vertex.uprime = -1.0;
+    em_vertex.vprime = -1.0;
+    em_normal.uprime = -1.0;
+    em_normal.vprime = -1.0;
+    em_color.uprime = -1.0;
+    em_color.vprime = -1.0;
+    em_texcoord.uprime = -1.0;
+    em_texcoord.vprime = -1.0;
+
+#ifdef USE_LOD
+    LOD_eval_level = 1;
+#endif
+}
+
+OpenGLSurfaceEvaluator::~OpenGLSurfaceEvaluator()
+{
+   for (int ii= 0; ii< VERTEX_CACHE_SIZE; ii++) {
+      delete vertexCache[ii];
+      vertexCache[ii]= 0;
+   }
+}
+
+/*---------------------------------------------------------------------------
+ * disable - turn off a map
+ *---------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::disable(long type)
+{
+    glDisable((GLenum) type);
+}
+
+/*---------------------------------------------------------------------------
+ * enable - turn on a map
+ *---------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::enable(long type)
+{
+    glEnable((GLenum) type);
+}
+
+/*-------------------------------------------------------------------------
+ * mapgrid2f - define a lattice of points with origin and offset
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::mapgrid2f(long nu, REAL u0, REAL u1, long nv, REAL v0, REAL v1)
+{
+#ifdef USE_INTERNAL_EVAL
+  inMapGrid2f((int) nu, (REAL) u0, (REAL) u1, (int) nv,
+	      (REAL) v0, (REAL) v1);
+#else
+
+  if(output_triangles)  
+    {
+      global_grid_u0 = u0;
+      global_grid_u1 = u1;
+      global_grid_nu = nu;
+      global_grid_v0 = v0;
+      global_grid_v1 = v1;
+      global_grid_nv = nv;
+    }
+  else
+    glMapGrid2d((GLint) nu, (GLdouble) u0, (GLdouble) u1, (GLint) nv,
+	    (GLdouble) v0, (GLdouble) v1);
+
+#endif
+}
+
+void
+OpenGLSurfaceEvaluator::polymode(long style)
+{
+  if(! output_triangles)
+    {
+      switch(style) {
+      default:
+      case N_MESHFILL:
+        
+	glPolygonMode((GLenum) GL_FRONT_AND_BACK, (GLenum) GL_FILL);
+	break;
+      case N_MESHLINE:
+	glPolygonMode((GLenum) GL_FRONT_AND_BACK, (GLenum) GL_LINE);
+	break;
+      case N_MESHPOINT:
+	glPolygonMode((GLenum) GL_FRONT_AND_BACK, (GLenum) GL_POINT);
+	break;
+      }
+    }
+}
+
+void
+OpenGLSurfaceEvaluator::bgnline(void)
+{
+  if(output_triangles)
+    bezierPatchMeshBeginStrip(global_bpm, GL_LINE_STRIP);
+  else
+    glBegin((GLenum) GL_LINE_STRIP);
+}
+
+void
+OpenGLSurfaceEvaluator::endline(void)
+{
+  if(output_triangles)
+    bezierPatchMeshEndStrip(global_bpm);
+  else
+    glEnd();
+}
+
+void
+OpenGLSurfaceEvaluator::range2f(long type, REAL *from, REAL *to)
+{
+}
+
+void
+OpenGLSurfaceEvaluator::domain2f(REAL ulo, REAL uhi, REAL vlo, REAL vhi)
+{
+}
+
+void
+OpenGLSurfaceEvaluator::bgnclosedline(void)
+{
+  if(output_triangles)
+    bezierPatchMeshBeginStrip(global_bpm, GL_LINE_LOOP);
+  else
+    glBegin((GLenum) GL_LINE_LOOP);
+}
+
+void
+OpenGLSurfaceEvaluator::endclosedline(void)
+{
+  if(output_triangles)
+    bezierPatchMeshEndStrip(global_bpm);
+  else
+    glEnd();
+}
+
+
+
+
+
+void
+OpenGLSurfaceEvaluator::bgntmesh(void)
+{
+
+    tmeshing = 1;
+    which = 0;
+    vcount = 0;
+
+    if(output_triangles)
+      bezierPatchMeshBeginStrip(global_bpm, GL_TRIANGLES);
+    else
+      glBegin((GLenum) GL_TRIANGLES);
+
+}
+
+void
+OpenGLSurfaceEvaluator::swaptmesh(void)
+{
+    which = 1 - which;
+
+}
+
+void
+OpenGLSurfaceEvaluator::endtmesh(void)
+{
+    tmeshing = 0;
+
+
+    if(output_triangles)
+      bezierPatchMeshEndStrip(global_bpm);
+    else
+      glEnd();
+}
+
+void
+OpenGLSurfaceEvaluator::bgntfan(void)
+{
+
+  if(output_triangles)
+    bezierPatchMeshBeginStrip(global_bpm, GL_TRIANGLE_FAN);
+  else
+    glBegin((GLenum) GL_TRIANGLE_FAN);
+
+}
+void
+OpenGLSurfaceEvaluator::endtfan(void)
+{
+  if(output_triangles)
+	bezierPatchMeshEndStrip(global_bpm);
+  else
+    glEnd();
+}
+
+void
+OpenGLSurfaceEvaluator::evalUStrip(int n_upper, REAL v_upper, REAL* upper_val, int n_lower, REAL v_lower, REAL* lower_val)
+{
+#ifdef USE_INTERNAL_EVAL
+  inEvalUStrip(n_upper, v_upper, upper_val,
+	n_lower, v_lower, lower_val);
+#else
+
+#ifdef FOR_CHRIS
+  evalUStripExt(n_upper, v_upper, upper_val,
+		 n_lower, v_lower, lower_val);
+  return;
+
+#endif
+  int i,j,k,l;
+  REAL leftMostV[2];
+
+  /*
+   *the algorithm works by scanning from left to right.
+   *leftMostV: the left most of the remaining verteces (on both upper and lower).
+   *	       it could an element of upperVerts or lowerVerts.
+   *i: upperVerts[i] is the first vertex to the right of leftMostV on upper line
+   *j: lowerVerts[j] is the first vertex to the right of leftMostV on lower line
+   */
+
+  /*initialize i,j,and leftMostV
+   */
+  if(upper_val[0] <= lower_val[0])
+    {
+      i=1;
+      j=0;
+
+      leftMostV[0] = upper_val[0];
+      leftMostV[1] = v_upper;
+    }
+  else
+    {
+      i=0;
+      j=1;
+
+      leftMostV[0] = lower_val[0];
+      leftMostV[1] = v_lower;
+
+    }
+
+  /*the main loop.
+   *the invariance is that:
+   *at the beginning of each loop, the meaning of i,j,and leftMostV are
+   *maintained
+   */
+  while(1)
+    {
+      if(i >= n_upper) /*case1: no more in upper*/
+	{
+	  if(j<n_lower-1) /*at least two vertices in lower*/
+	    {
+	      bgntfan();
+	      coord2f(leftMostV[0], leftMostV[1]);
+//	      glNormal3fv(leftMostNormal);
+//		glVertex3fv(leftMostXYZ);
+
+	      while(j<n_lower){
+		coord2f(lower_val[j], v_lower);
+//		glNormal3fv(lowerNormal[j]);
+//		glVertex3fv(lowerXYZ[j]);
+		j++;
+
+	      }
+	      endtfan();
+	    }
+	  break; /*exit the main loop*/
+	}
+      else if(j>= n_lower) /*case2: no more in lower*/
+	{
+	  if(i<n_upper-1) /*at least two vertices in upper*/
+	    {
+	      bgntfan();
+	      coord2f(leftMostV[0], leftMostV[1]);
+//	      glNormal3fv(leftMostNormal);
+//	      glVertex3fv(leftMostXYZ);
+        
+	      for(k=n_upper-1; k>=i; k--) /*reverse order for two-side lighting*/
+		{
+		  coord2f(upper_val[k], v_upper);
+//		  glNormal3fv(upperNormal[k]);
+//		  glVertex3fv(upperXYZ[k]);
+		}
+
+	      endtfan();
+	    }
+	  break; /*exit the main loop*/
+	}
+      else /* case3: neither is empty, plus the leftMostV, there is at least one triangle to output*/
+	{
+	  if(upper_val[i] <= lower_val[j])
+	    {
+	      bgntfan();
+	      coord2f(lower_val[j], v_lower);
+//	      glNormal3fv(lowerNormal[j]);
+//	      glVertex3fv(lowerXYZ[j]);
+
+	      /*find the last k>=i such that
+	       *upperverts[k][0] <= lowerverts[j][0]
+	       */
+	      k=i;
+
+	      while(k<n_upper)
+		{
+		  if(upper_val[k] > lower_val[j])
+		    break;
+		  k++;
+
+		}
+	      k--;
+
+
+	      for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
+		{
+		  coord2f(upper_val[l], v_upper);
+//		  glNormal3fv(upperNormal[l]);
+//		  glVertex3fv(upperXYZ[l]);
+
+		}
+	      coord2f(leftMostV[0], leftMostV[1]);
+//	      glNormal3fv(leftMostNormal);
+//	      glVertex3fv(leftMostXYZ);
+
+	      endtfan();
+
+	      /*update i and leftMostV for next loop
+	       */
+	      i = k+1;
+
+	      leftMostV[0] = upper_val[k];
+	      leftMostV[1] = v_upper;
+//	      leftMostNormal = upperNormal[k];
+//	      leftMostXYZ = upperXYZ[k];
+	    }
+	  else /*upperVerts[i][0] > lowerVerts[j][0]*/
+	    {
+	      bgntfan();
+	      coord2f(upper_val[i], v_upper);
+//	      glNormal3fv(upperNormal[i]);
+//	      glVertex3fv(upperXYZ[i]);
+        
+	      coord2f(leftMostV[0], leftMostV[1]);
+//		glNormal3fv(leftMostNormal);
+//	      glVertex3fv(leftMostXYZ);
+        
+
+	      /*find the last k>=j such that
+	       *lowerverts[k][0] < upperverts[i][0]
+	       */
+	      k=j;
+	      while(k< n_lower)
+		{
+		  if(lower_val[k] >= upper_val[i])
+		    break;
+		  coord2f(lower_val[k], v_lower);
+//		  glNormal3fv(lowerNormal[k]);
+//		  glVertex3fv(lowerXYZ[k]);
+
+		  k++;
+		}
+	      endtfan();
+
+	      /*update j and leftMostV for next loop
+	       */
+	      j=k;
+	      leftMostV[0] = lower_val[j-1];
+	      leftMostV[1] = v_lower;
+
+//	      leftMostNormal = lowerNormal[j-1];
+//	      leftMostXYZ = lowerXYZ[j-1];
+	    }
+	}
+    }
+  //clean up
+//  free(upperXYZ);
+//  free(lowerXYZ);
+//  free(upperNormal);
+//  free(lowerNormal);
+#endif
+
+}
+
+
+void
+OpenGLSurfaceEvaluator::evalVStrip(int n_left, REAL u_left, REAL* left_val, int n_right, REAL u_right, REAL* right_val)
+{
+#ifdef USE_INTERNAL_EVAL
+	inEvalVStrip(n_left, u_left, left_val,
+	n_right, u_right, right_val);
+#else
+
+#ifdef FOR_CHRIS
+	evalVStripExt(n_left, u_left, left_val,
+		      n_right, u_right, right_val);
+	return;
+
+#endif
+
+  int i,j,k,l;
+  REAL botMostV[2];
+  /*
+   *the algorithm works by scanning from bot to top.
+   *botMostV: the bot most of the remaining verteces (on both left and right).
+   *	       it could an element of leftVerts or rightVerts.
+   *i: leftVerts[i] is the first vertex to the top of botMostV on left line
+   *j: rightVerts[j] is the first vertex to the top of botMostV on rightline
+   */
+
+  /*initialize i,j,and botMostV
+   */
+  if(left_val[0] <= right_val[0])
+    {
+      i=1;
+      j=0;
+
+      botMostV[0] = u_left;
+      botMostV[1] = left_val[0];
+    }
+  else
+    {
+      i=0;
+      j=1;
+
+      botMostV[0] = u_right;
+      botMostV[1] = right_val[0];
+    }
+
+  /*the main loop.
+   *the invariance is that:
+   *at the beginning of each loop, the meaning of i,j,and botMostV are
+   *maintained
+   */
+  while(1)
+    {
+      if(i >= n_left) /*case1: no more in left*/
+	{
+	  if(j<n_right-1) /*at least two vertices in right*/
+	    {
+	      bgntfan();
+	      coord2f(botMostV[0], botMostV[1]);
+	      while(j<n_right){
+		coord2f(u_right, right_val[j]);
+//		glNormal3fv(rightNormal[j]);
+//		glVertex3fv(rightXYZ[j]);
+		j++;
+
+	      }
+	      endtfan();
+	    }
+	  break; /*exit the main loop*/
+	}
+      else if(j>= n_right) /*case2: no more in right*/
+	{
+	  if(i<n_left-1) /*at least two vertices in left*/
+	    {
+	      bgntfan();
+	      coord2f(botMostV[0], botMostV[1]);
+//	      glNormal3fv(botMostNormal);
+//	      glVertex3fv(botMostXYZ);
+        
+	      for(k=n_left-1; k>=i; k--) /*reverse order for two-side lighting*/
+		{
+		  coord2f(u_left, left_val[k]);
+//		  glNormal3fv(leftNormal[k]);
+//		  glVertex3fv(leftXYZ[k]);
+		}
+
+	      endtfan();
+	    }
+	  break; /*exit the main loop*/
+	}
+      else /* case3: neither is empty, plus the botMostV, there is at least one triangle to output*/
+	{
+	  if(left_val[i] <= right_val[j])
+	    {
+	      bgntfan();
+	      coord2f(u_right, right_val[j]);
+//	      glNormal3fv(rightNormal[j]);
+//	      glVertex3fv(rightXYZ[j]);
+
+	      /*find the last k>=i such that
+	       *leftverts[k][0] <= rightverts[j][0]
+	       */
+	      k=i;
+
+	      while(k<n_left)
+		{
+		  if(left_val[k] > right_val[j])
+		    break;
+		  k++;
+
+		}
+	      k--;
+
+
+	      for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
+		{
+		  coord2f(u_left, left_val[l]);
+//		  glNormal3fv(leftNormal[l]);
+//		  glVertex3fv(leftXYZ[l]);
+
+		}
+	      coord2f(botMostV[0], botMostV[1]);
+//	      glNormal3fv(botMostNormal);
+//	      glVertex3fv(botMostXYZ);
+
+	      endtfan();
+
+	      /*update i and botMostV for next loop
+	       */
+	      i = k+1;
+
+	      botMostV[0] = u_left;
+	      botMostV[1] = left_val[k];
+//	      botMostNormal = leftNormal[k];
+//	      botMostXYZ = leftXYZ[k];
+	    }
+	  else /*left_val[i] > right_val[j])*/
+	    {
+	      bgntfan();
+	      coord2f(u_left, left_val[i]);
+//	      glNormal3fv(leftNormal[i]);
+//	      glVertex3fv(leftXYZ[i]);
+        
+	      coord2f(botMostV[0], botMostV[1]);
+//	      glNormal3fv(botMostNormal);
+//	      glVertex3fv(botMostXYZ);
+        
+
+	      /*find the last k>=j such that
+	       *rightverts[k][0] < leftverts[i][0]
+	       */
+	      k=j;
+	      while(k< n_right)
+		{
+		  if(right_val[k] >= left_val[i])
+		    break;
+		  coord2f(u_right, right_val[k]);
+//		  glNormal3fv(rightNormal[k]);
+//		  glVertex3fv(rightXYZ[k]);
+
+		  k++;
+		}
+	      endtfan();
+
+	      /*update j and botMostV for next loop
+	       */
+	      j=k;
+	      botMostV[0] = u_right;
+	      botMostV[1] = right_val[j-1];
+
+//	      botMostNormal = rightNormal[j-1];
+//	      botMostXYZ = rightXYZ[j-1];
+	    }
+	}
+    }
+  //clean up
+//  free(leftXYZ);
+//  free(leftNormal);
+//  free(rightXYZ);
+//  free(rightNormal);
+#endif
+}
+
+
+void
+OpenGLSurfaceEvaluator::bgnqstrip(void)
+{
+  if(output_triangles)
+    bezierPatchMeshBeginStrip(global_bpm, GL_QUAD_STRIP);
+  else
+    glBegin((GLenum) GL_QUAD_STRIP);
+
+#ifdef STATISTICS
+	STAT_num_of_quad_strips++;
+#endif
+}
+
+void
+OpenGLSurfaceEvaluator::endqstrip(void)
+{
+  if(output_triangles)
+    bezierPatchMeshEndStrip(global_bpm);
+  else
+    glEnd();
+
+}
+
+/*-------------------------------------------------------------------------
+ * bgnmap2f - preamble to surface definition and evaluations
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::bgnmap2f(long)
+{
+  if(output_triangles)
+    {
+      /*deallocate the space which may has been
+       *allocated by global_bpm previously
+       */
+      if(global_bpm != NULL) {
+	bezierPatchMeshListDelete(global_bpm);
+	global_bpm = NULL;
+      }
+
+
+      /*
+	auto_normal_flag = 1; //always output normal in callback mode.
+			      //we could have used the following code,
+			      //but Inspector doesn't have gl context
+			      //before it calls tessellator.
+			      //this way is temporary.
+	*/
+      //NEWCALLBACK
+      //if one of the two normal callback functions are set,
+      //then set
+      if(normalCallBackN != NULL ||
+	 normalCallBackData != NULL)
+	auto_normal_flag = 1;
+      else
+	auto_normal_flag = 0;
+
+      //initialize so that no maps initially
+      vertex_flag = 0;
+      normal_flag = 0;
+      color_flag = 0;
+      texcoord_flag = 0;
+
+      /*
+      if(glIsEnabled(GL_AUTO_NORMAL) == GL_TRUE)
+	auto_normal_flag = 1;
+      else if (callback_auto_normal == 1)
+	auto_normal_flag = 1;
+      else
+	auto_normal_flag = 0;
+	*/
+	  glPushAttrib((GLbitfield) GL_EVAL_BIT);
+
+    }
+  else
+    {
+      glPushAttrib((GLbitfield) GL_EVAL_BIT);
+
+      /*to avoid side effect, we restor the opengl state for GL_POLYGON_MODE
+       */       
+      glGetIntegerv(GL_POLYGON_MODE, gl_polygon_mode);
+    }
+
+}
+
+/*-------------------------------------------------------------------------
+ * endmap2f - postamble to a map
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::endmap2f(void)
+{
+
+  if(output_triangles)
+    {
+      //bezierPatchMeshListDelDeg(global_bpm);
+
+      //    bezierPatchMeshListEval(global_bpm);
+
+      //surfcount++;
+      //printf("surfcount=%i\n", surfcount);
+      //if(surfcount == 8) exit(0);
+
+      inBPMListEvalEM(global_bpm);
+
+
+
+/*
+    global_bpm = bezierPatchMeshListReverse(global_bpm);
+    {
+      float *vertex_array;
+      float *normal_array;
+      int *length_array;
+      int *type_array;
+      int num_strips;
+      bezierPatchMeshListCollect(global_bpm, &vertex_array, &normal_array, &length_array, &type_array, &num_strips);
+      drawStrips(vertex_array, normal_array, length_array, type_array, num_strips);
+      free(vertex_array);
+      free(normal_array);
+      free(length_array);
+      free(type_array);
+    }
+*/
+
+    //bezierPatchMeshListPrint(global_bpm);
+    //bezierPatchMeshListDraw(global_bpm);
+
+//	  printf("num triangles=%i\n", bezierPatchMeshListNumTriangles(global_bpm));
+
+#ifdef USE_LOD
+#else
+    bezierPatchMeshListDelete(global_bpm);
+    global_bpm = NULL;
+#endif
+	glPopAttrib();
+  }
+else
+  {
+#ifndef USE_LOD
+    glPopAttrib();
+#endif
+
+#ifdef STATISTICS
+    fprintf(stderr, "num_vertices=%i,num_triangles=%i,num_quads_strips=%i\n", STAT_num_of_eval_vertices,STAT_num_of_triangles,STAT_num_of_quad_strips);
+#endif
+
+    /*to restore the gl_polygon_mode
+     */
+#ifndef USE_LOD
+    glPolygonMode( GL_FRONT, (GLenum) gl_polygon_mode[0]);
+    glPolygonMode( GL_BACK,  (GLenum) gl_polygon_mode[1]);
+#endif
+}
+
+}
+
+/*-------------------------------------------------------------------------
+ * map2f - pass a desription of a surface map
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::map2f(
+    long _type,
+    REAL _ulower,	/* u lower domain coord 	*/
+    REAL _uupper,	/* u upper domain coord 	*/
+    long _ustride,	/* interpoint distance		*/
+    long _uorder,	/* parametric order		*/
+    REAL _vlower,	/* v lower domain coord 	*/
+    REAL _vupper,	/* v upper domain coord 	*/
+    long _vstride,	/* interpoint distance		*/
+    long _vorder,	/* parametric order		*/
+    REAL *pts)	/* control points		*/
+{
+#ifdef USE_INTERNAL_EVAL
+   inMap2f((int) _type, (REAL) _ulower, (REAL) _uupper,
+	    (int) _ustride, (int) _uorder, (REAL) _vlower,
+	    (REAL) _vupper, (int) _vstride, (int) _vorder,
+	    (REAL *) pts);
+#else
+
+
+
+   if(output_triangles)
+     {
+       if(global_bpm == NULL)
+	 global_bpm = bezierPatchMeshMake2(10,10);
+       if(
+	  (global_bpm->bpatch == NULL &&
+	  (_type == GL_MAP2_VERTEX_3 || _type == GL_MAP2_VERTEX_4))
+	  ||
+	  (global_bpm->bpatch_normal == NULL &&
+	   (_type == GL_MAP2_NORMAL))
+	  ||
+	  (global_bpm->bpatch_color == NULL &&
+	   (_type == GL_MAP2_INDEX || _type == GL_MAP2_COLOR_4))
+	  ||
+	  (global_bpm->bpatch_texcoord == NULL &&
+	   (_type == GL_MAP2_TEXTURE_COORD_1 ||
+	    _type == GL_MAP2_TEXTURE_COORD_2 ||
+	    _type == GL_MAP2_TEXTURE_COORD_3 ||
+	    _type == GL_MAP2_TEXTURE_COORD_4 )
+	   ))
+	 {
+	   bezierPatchMeshPutPatch(global_bpm, (int) _type, _ulower, _uupper,(int)  _ustride,(int) _uorder,_vlower, _vupper, (int) _vstride, (int) _vorder, pts);
+	 }
+       else /*new surface patch (with multiple maps) starts*/
+	 {
+	   bezierPatchMesh *temp = bezierPatchMeshMake2(10,10);
+	   bezierPatchMeshPutPatch(temp, (int) _type, _ulower, _uupper,(int)  _ustride,(int) _uorder,_vlower, _vupper, (int) _vstride, (int) _vorder, pts);
+	   global_bpm = bezierPatchMeshListInsert(global_bpm, temp);
+
+	   /*
+	   global_bpm = bezierPatchMeshListInsert(global_bpm,
+						  bezierPatchMeshMake(
+								      (int) _type, _ulower, _uupper,(int)  _ustride, (int) _uorder, _vlower, _vupper, (int) _vstride, (int) _vorder, pts, 10, 10));
+	   */
+	 }
+     }
+   else /*not output triangles*/
+     {
+       glMap2f((GLenum) _type, (GLfloat) _ulower, (GLfloat) _uupper,
+	       (GLint) _ustride, (GLint) _uorder, (GLfloat) _vlower,
+	       (GLfloat) _vupper, (GLint) _vstride, (GLint) _vorder,
+	       (const GLfloat *) pts);
+     }
+
+#endif
+}
+
+
+/*-------------------------------------------------------------------------
+ * mapmesh2f - evaluate a mesh of points on lattice
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::mapmesh2f(long style, long umin, long umax, long vmin, long vmax)
+{
+#ifdef NO_EVALUATION
+return;
+#endif
+
+#ifdef USE_INTERNAL_EVAL
+    inEvalMesh2((int)umin, (int)vmin, (int)umax, (int)vmax);
+#else
+
+
+
+if(output_triangles)
+{
+#ifdef USE_LOD
+  bezierPatchMeshBeginStrip(global_bpm, GL_POLYGON);
+  bezierPatchMeshInsertUV(global_bpm, global_grid_u0, global_grid_v0);
+  bezierPatchMeshInsertUV(global_bpm, global_grid_u1, global_grid_v1);
+  bezierPatchMeshInsertUV(global_bpm, (REAL)global_grid_nu, (REAL)global_grid_nv);
+  bezierPatchMeshInsertUV(global_bpm, (REAL)umin, (REAL)vmin);
+  bezierPatchMeshInsertUV(global_bpm, (REAL)umax, (REAL)vmax);
+  bezierPatchMeshEndStrip(global_bpm);
+
+#else
+
+  REAL du, dv;
+  long i,j;
+  if(global_grid_nu == 0 || global_grid_nv == 0)
+    return; /*no points need to be output*/
+  du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
+  dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;
+
+  if(global_grid_nu >= global_grid_nv){
+
+    for(i=umin; i<umax; i++){
+      REAL u1 = (i==global_grid_nu)? global_grid_u1:(global_grid_u0 + i*du);
+      REAL u2 = ((i+1) == global_grid_nu)? global_grid_u1: (global_grid_u0+(i+1)*du);
+
+      bgnqstrip();
+      for(j=vmax; j>=vmin; j--){
+	REAL v1 = (j == global_grid_nv)? global_grid_v1: (global_grid_v0 +j*dv);
+        
+	coord2f(u1, v1);
+	coord2f(u2, v1);
+      }
+      endqstrip();
+    }
+  }
+  else{
+
+    for(i=vmin; i<vmax; i++){
+      REAL v1 = (i==global_grid_nv)? global_grid_v1:(global_grid_v0 + i*dv);
+      REAL v2 = ((i+1) == global_grid_nv)? global_grid_v1: (global_grid_v0+(i+1)*dv);
+
+      bgnqstrip();
+      for(j=umax; j>=umin; j--){
+	REAL u1 = (j == global_grid_nu)? global_grid_u1: (global_grid_u0 +j*du);        
+	coord2f(u1, v2);
+	coord2f(u1, v1);
+      }
+      endqstrip();
+    }
+  }
+#endif
+}
+else
+{
+    switch(style) {
+    default:
+    case N_MESHFILL:
+	glEvalMesh2((GLenum) GL_FILL, (GLint) umin, (GLint) umax,
+		(GLint) vmin, (GLint) vmax);
+	break;
+    case N_MESHLINE:
+	glEvalMesh2((GLenum) GL_LINE, (GLint) umin, (GLint) umax,
+		(GLint) vmin, (GLint) vmax);
+	break;
+    case N_MESHPOINT:
+	glEvalMesh2((GLenum) GL_POINT, (GLint) umin, (GLint) umax,
+		(GLint) vmin, (GLint) vmax);
+	break;
+    }
+  }
+
+#endif
+
+#ifdef STATISTICS
+	STAT_num_of_quad_strips += (umax-umin)*(vmax-vmin);
+#endif
+}
+
+/*-------------------------------------------------------------------------
+ * evalcoord2f - evaluate a point on a surface
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::evalcoord2f(long, REAL u, REAL v)
+{
+
+
+#ifdef NO_EVALUATION
+return;
+#endif
+
+
+    newtmeshvert(u, v);
+}
+
+/*-------------------------------------------------------------------------
+ * evalpoint2i - evaluate a grid point
+ *-------------------------------------------------------------------------
+ */
+void
+OpenGLSurfaceEvaluator::evalpoint2i(long u, long v)
+{
+#ifdef NO_EVALUATION
+return;
+#endif
+
+    newtmeshvert(u, v);
+}
+
+void
+OpenGLSurfaceEvaluator::point2i( long u, long v )
+{
+#ifdef NO_EVALUATION
+return;
+#else
+
+#ifdef USE_INTERNAL_EVAL
+    inEvalPoint2( (int)u,  (int)v);
+#else
+
+
+if(output_triangles)
+{
+
+  REAL du, dv;
+  REAL fu,fv;
+  du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
+  dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;
+  fu = (u==global_grid_nu)? global_grid_u1:(global_grid_u0 + u*du);
+  fv = (v == global_grid_nv)? global_grid_v1: (global_grid_v0 +v*dv);
+  coord2f(fu,fv);
+}
+else
+    glEvalPoint2((GLint) u, (GLint) v);
+
+
+#endif
+
+#ifdef STATISTICS
+  STAT_num_of_eval_vertices++;
+#endif
+
+#endif
+
+}
+
+void
+OpenGLSurfaceEvaluator::coord2f( REAL u, REAL v )
+{
+#ifdef NO_EVALUATION
+return;
+#else
+
+#ifdef USE_INTERNAL_EVAL
+    inEvalCoord2f( u, v);
+#else
+
+
+if(output_triangles)
+    bezierPatchMeshInsertUV(global_bpm, u,v);
+else
+    glEvalCoord2f((GLfloat) u, (GLfloat) v);
+
+
+#endif
+
+
+#ifdef STATISTICS
+  STAT_num_of_eval_vertices++;
+#endif
+
+#endif
+}
+
+void
+OpenGLSurfaceEvaluator::newtmeshvert( long u, long v )
+{
+#ifdef NO_EVALUATION
+return;
+#else
+
+    if (tmeshing) {
+
+	if (vcount == 2) {
+	    vertexCache[0]->invoke(this);
+	    vertexCache[1]->invoke(this);
+	    point2i( u,  v);
+
+	} else {
+	    vcount++;
+	}
+
+	vertexCache[which]->saveEvalPoint(u, v);
+	which = 1 - which;
+    } else {
+	point2i( u,  v);
+    }
+#endif
+}
+
+void
+OpenGLSurfaceEvaluator::newtmeshvert( REAL u, REAL v )
+{
+#ifdef NO_EVALUATION
+return;
+#else
+    if (tmeshing) {
+
+
+	if (vcount == 2) {
+	    vertexCache[0]->invoke(this);
+	    vertexCache[1]->invoke(this);
+	    coord2f(u,v);
+
+	} else {
+	    vcount++;
+	}
+
+	vertexCache[which]->saveEvalCoord(u, v);
+	which = 1 - which;
+    } else {
+
+	coord2f( u,  v);
+    }
+#endif
+
+}
+
+#ifdef _WIN32
+void OpenGLSurfaceEvaluator::putCallBack(GLenum which, void (GLAPIENTRY *fn)() )
+#else
+void OpenGLSurfaceEvaluator::putCallBack(GLenum which, _GLUfuncptr fn )
+#endif
+{
+  switch(which)
+    {
+    case GLU_NURBS_BEGIN:
+      beginCallBackN = (void (GLAPIENTRY *) (GLenum)) fn;
+      break;
+    case GLU_NURBS_END:
+      endCallBackN = (void (GLAPIENTRY *) (void)) fn;
+      break;
+    case GLU_NURBS_VERTEX:
+      vertexCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_NORMAL:
+      normalCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_COLOR:
+      colorCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_TEXTURE_COORD:
+      texcoordCallBackN = (void (GLAPIENTRY *) (const GLfloat*)) fn;
+      break;
+    case GLU_NURBS_BEGIN_DATA:
+      beginCallBackData = (void (GLAPIENTRY *) (GLenum, void*)) fn;
+      break;
+    case GLU_NURBS_END_DATA:
+      endCallBackData = (void (GLAPIENTRY *) (void*)) fn;
+      break;
+    case GLU_NURBS_VERTEX_DATA:
+      vertexCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_NORMAL_DATA:
+      normalCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_COLOR_DATA:
+      colorCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+    case GLU_NURBS_TEXTURE_COORD_DATA:
+      texcoordCallBackData = (void (GLAPIENTRY *) (const GLfloat*, void*)) fn;
+      break;
+
+    }
+}
+
+
+void
+OpenGLSurfaceEvaluator::beginCallBack(GLenum which, void *data)
+{
+  if(beginCallBackData)
+    beginCallBackData(which, data);
+  else if(beginCallBackN)
+    beginCallBackN(which);
+}
+
+void
+OpenGLSurfaceEvaluator::endCallBack(void *data)
+{
+  if(endCallBackData)
+    endCallBackData(data);
+  else if(endCallBackN)
+    endCallBackN();
+}
+
+void
+OpenGLSurfaceEvaluator::vertexCallBack(const GLfloat *vert, void* data)
+{
+  if(vertexCallBackData)
+    vertexCallBackData(vert, data);
+  else if(vertexCallBackN)
+    vertexCallBackN(vert);
+}
+
+
+void
+OpenGLSurfaceEvaluator::normalCallBack(const GLfloat *normal, void* data)
+{
+  if(normalCallBackData)
+    normalCallBackData(normal, data);
+  else if(normalCallBackN)
+    normalCallBackN(normal);
+}
+
+void
+OpenGLSurfaceEvaluator::colorCallBack(const GLfloat *color, void* data)
+{
+  if(colorCallBackData)
+    colorCallBackData(color, data);
+  else if(colorCallBackN)
+    colorCallBackN(color);
+}
+
+void
+OpenGLSurfaceEvaluator::texcoordCallBack(const GLfloat *texcoord, void* data)
+{
+  if(texcoordCallBackData)
+    texcoordCallBackData(texcoord, data);
+  else if(texcoordCallBackN)
+    texcoordCallBackN(texcoord);
+}
+
+
+
+
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.h b/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.h
new file mode 100644
index 000000000..1567c6b09
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/glsurfeval.h
@@ -0,0 +1,404 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * glsurfeval.h
+ *
+ */
+
+#ifndef __gluglsurfeval_h_
+#define __gluglsurfeval_h_
+
+#include "basicsurfeval.h"
+#include "bezierPatchMesh.h" //in case output triangles
+#include <GL/gl.h>
+#include <GL/glu.h>
+
+class SurfaceMap;
+class OpenGLSurfaceEvaluator;
+class StoredVertex;
+
+#define TYPECOORD	1
+#define TYPEPOINT	2
+
+/* Cache up to 3 vertices from tmeshes */
+#define VERTEX_CACHE_SIZE	3
+
+/*for internal evaluator callback stuff*/
+#ifndef IN_MAX_BEZIER_ORDER
+#define IN_MAX_BEZIER_ORDER 40 /*XXX should be bigger than machine order*/
+#endif
+			
+#ifndef IN_MAX_DIMENSION
+#define IN_MAX_DIMENSION 4 
+#endif
+
+typedef struct surfEvalMachine{
+  REAL uprime;//cached previusly evaluated uprime.
+  REAL vprime;
+  int k; /*the dimension*/
+  REAL u1;
+  REAL u2;
+  int ustride;
+  int uorder;
+  REAL v1;
+  REAL v2;
+  int vstride;
+  int vorder;
+  REAL ctlPoints[IN_MAX_BEZIER_ORDER*IN_MAX_BEZIER_ORDER*IN_MAX_DIMENSION];
+  REAL ucoeff[IN_MAX_BEZIER_ORDER]; /*cache the polynomial values*/
+  REAL vcoeff[IN_MAX_BEZIER_ORDER];
+  REAL ucoeffDeriv[IN_MAX_BEZIER_ORDER]; /*cache the polynomial derivatives*/
+  REAL vcoeffDeriv[IN_MAX_BEZIER_ORDER];
+} surfEvalMachine;
+  
+  
+
+class StoredVertex {
+public:
+    		StoredVertex() { type = 0; }
+		~StoredVertex(void) {}
+    void	saveEvalCoord(REAL x, REAL y) 
+		    {coord[0] = x; coord[1] = y; type = TYPECOORD; }
+    void	saveEvalPoint(long x, long y)
+		    {point[0] = x; point[1] = y; type = TYPEPOINT; }
+    void	invoke(OpenGLSurfaceEvaluator *eval);
+
+private:
+    int		type;
+    REAL	coord[2];
+    long	point[2];
+};
+
+class OpenGLSurfaceEvaluator : public BasicSurfaceEvaluator {
+public:
+			OpenGLSurfaceEvaluator();
+    			virtual ~OpenGLSurfaceEvaluator( void );
+    void		polymode( long style );
+    void		range2f( long, REAL *, REAL * );
+    void		domain2f( REAL, REAL, REAL, REAL );
+    void		addMap( SurfaceMap * ) { }
+
+    void		enable( long );
+    void		disable( long );
+    void		bgnmap2f( long );
+    void		map2f( long, REAL, REAL, long, long, 
+				     REAL, REAL, long, long, REAL * );
+    void		mapgrid2f( long, REAL, REAL, long, REAL, REAL );
+    void		mapmesh2f( long, long, long, long, long );
+    void		evalcoord2f( long, REAL, REAL );
+    void		evalpoint2i( long, long );
+    void		endmap2f( void );
+
+    void	 	bgnline( void );
+    void	 	endline( void );
+    void	 	bgnclosedline( void );
+    void	 	endclosedline( void );
+    void	 	bgntmesh( void );
+    void	 	swaptmesh( void );
+    void	 	endtmesh( void );
+    void	 	bgnqstrip( void );
+    void	 	endqstrip( void );
+
+    void                bgntfan( void );
+    void                endtfan( void );
+    void                evalUStrip(int n_upper, REAL v_upper, REAL* upper_val,
+                                   int n_lower, REAL v_lower, REAL* lower_val);
+    void                evalVStrip(int n_left, REAL u_left, REAL* left_val,
+                                   int n_right, REAL u_right, REAL* right_val);
+
+    void		coord2f( REAL, REAL );
+    void		point2i( long, long );
+
+    void		newtmeshvert( REAL, REAL );
+    void		newtmeshvert( long, long );
+
+#ifdef _WIN32
+    void 	        putCallBack(GLenum which, void (GLAPIENTRY *fn)() );
+#else
+    void 	        putCallBack(GLenum which, _GLUfuncptr fn );
+#endif
+
+    int                 get_vertices_call_back()
+      {
+	return output_triangles;
+      }
+    void                put_vertices_call_back(int flag)
+      {
+	output_triangles = flag;
+      }
+
+    void                 put_callback_auto_normal(int flag)
+      {
+        callback_auto_normal = flag;
+      } 
+
+   int                   get_callback_auto_normal()
+     {
+        return callback_auto_normal;
+      }
+
+   void                  set_callback_userData(void* data)
+     {
+       userData = data;
+     }
+
+    /**************begin for LOD_eval_list***********/
+    void LOD_eval_list(int level);
+
+
+
+   
+private:
+    StoredVertex	*vertexCache[VERTEX_CACHE_SIZE];
+    int			tmeshing;
+    int			which;
+    int			vcount;
+
+    GLint              gl_polygon_mode[2];/*to save and restore so that
+					 *no side effect
+					 */
+    bezierPatchMesh        *global_bpm; //for output triangles
+    int                output_triangles; //true 1 or false 0
+
+
+
+    void (GLAPIENTRY *beginCallBackN) (GLenum type);
+    void (GLAPIENTRY *endCallBackN)   (void);
+    void (GLAPIENTRY *vertexCallBackN) (const GLfloat *vert);
+    void (GLAPIENTRY *normalCallBackN) (const GLfloat *normal);
+    void (GLAPIENTRY *colorCallBackN) (const GLfloat *color);
+    void (GLAPIENTRY *texcoordCallBackN) (const GLfloat *texcoord);
+
+    void (GLAPIENTRY *beginCallBackData) (GLenum type, void* data);
+    void (GLAPIENTRY *endCallBackData)   (void* data);
+    void (GLAPIENTRY *vertexCallBackData) (const GLfloat *vert, void* data);
+    void (GLAPIENTRY *normalCallBackData) (const GLfloat *normal, void* data);
+    void (GLAPIENTRY *colorCallBackData) (const GLfloat *color, void* data);
+    void (GLAPIENTRY *texcoordCallBackData) (const GLfloat *texcoord, void* data);
+
+    void               beginCallBack (GLenum type, void* data);
+    void               endCallBack   (void* data);
+    void               vertexCallBack (const GLfloat *vert, void* data);
+    void               normalCallBack (const GLfloat *normal, void* data);
+    void               colorCallBack (const GLfloat *color, void* data);
+    void               texcoordCallBack (const GLfloat *texcoord, void* data);
+
+
+    void* userData; //the opaque pointer for Data callback functions.
+
+   /*LOD evaluation*/
+   void LOD_triangle(REAL A[2], REAL B[2], REAL C[2],
+		     int level);
+   void LOD_eval(int num_vert, REAL* verts, int type, int level);
+		     
+  int LOD_eval_level; //set by LOD_eval_list()
+
+   /*************begin for internal evaluators*****************/
+			
+ /*the following global variables are only defined in this file. 
+ *They are used to cache the precomputed Bezier polynomial values.
+ *These calues may be used consecutively in which case we don't have 
+ *recompute these values again.
+ */
+ int global_uorder; /*store the uorder in the previous evaluation*/
+ int global_vorder; /*store the vorder in the previous evaluation*/
+ REAL global_uprime;
+ REAL global_vprime;
+ REAL global_vprime_BV;
+ REAL global_uprime_BU;
+ int global_uorder_BV; /*store the uorder in the previous evaluation*/
+ int global_vorder_BV; /*store the vorder in the previous evaluation*/
+ int global_uorder_BU; /*store the uorder in the previous evaluation*/
+ int global_vorder_BU; /*store the vorder in the previous evaluation*/
+
+ REAL global_ucoeff[IN_MAX_BEZIER_ORDER]; /*cache the polynomial values*/
+ REAL global_vcoeff[IN_MAX_BEZIER_ORDER];
+ REAL global_ucoeffDeriv[IN_MAX_BEZIER_ORDER]; /*cache the polynomial derivatives*/
+ REAL global_vcoeffDeriv[IN_MAX_BEZIER_ORDER];
+
+ REAL global_BV[IN_MAX_BEZIER_ORDER][IN_MAX_DIMENSION];
+ REAL global_PBV[IN_MAX_BEZIER_ORDER][IN_MAX_DIMENSION];
+ REAL global_BU[IN_MAX_BEZIER_ORDER][IN_MAX_DIMENSION];
+ REAL global_PBU[IN_MAX_BEZIER_ORDER][IN_MAX_DIMENSION];
+ REAL* global_baseData;
+
+ int    global_ev_k; /*the dimension*/
+ REAL global_ev_u1;
+ REAL global_ev_u2;
+ int    global_ev_ustride;
+ int    global_ev_uorder;
+ REAL global_ev_v1;
+ REAL global_ev_v2;
+ int    global_ev_vstride;
+ int    global_ev_vorder;
+ REAL global_ev_ctlPoints[IN_MAX_BEZIER_ORDER*IN_MAX_BEZIER_ORDER*IN_MAX_DIMENSION];
+
+ REAL  global_grid_u0;
+ REAL  global_grid_u1;
+ int     global_grid_nu;
+ REAL  global_grid_v0;
+ REAL  global_grid_v1;
+ int     global_grid_nv;
+
+/*functions*/
+ void inDoDomain2WithDerivs(int k, REAL u, REAL v, 
+				REAL u1, REAL u2, int uorder, 
+				REAL v1,  REAL v2, int vorder, 
+				REAL *baseData,
+				REAL *retPoint, REAL *retdu, REAL *retdv);
+ void inPreEvaluate(int order, REAL vprime, REAL *coeff);
+ void inPreEvaluateWithDeriv(int order, REAL vprime, REAL *coeff, REAL *coeffDeriv);
+ void inComputeFirstPartials(REAL *p, REAL *pu, REAL *pv);
+ void inComputeNormal2(REAL *pu, REAL *pv, REAL *n);
+ void inDoEvalCoord2(REAL u, REAL v,
+		     REAL *retPoint, REAL *retNormal);
+ void inDoEvalCoord2NOGE(REAL u, REAL v,
+		     REAL *retPoint, REAL *retNormal);
+ void inMap2f(int k,
+	      REAL ulower,
+	      REAL uupper,
+	      int ustride,
+	      int uorder,
+	      REAL vlower,
+	      REAL vupper,
+	      int vstride,
+	      int vorder,
+	      REAL *ctlPoints);
+
+ void inMapGrid2f(int nu, REAL u0, REAL u1, 
+		  int nv, REAL v0, REAL v1);
+
+ void inEvalMesh2(int lowU, int lowV, int highU, int highV);
+ void inEvalPoint2(int i, int j);
+ void inEvalCoord2f(REAL u, REAL v);
+
+void inEvalULine(int n_points, REAL v, REAL* u_vals, 
+	int stride, REAL ret_points[][3], REAL ret_normals[][3]);
+
+void inEvalVLine(int n_points, REAL u, REAL* v_vals, 
+	int stride, REAL ret_points[][3], REAL ret_normals[][3]);
+
+void inEvalUStrip(int n_upper, REAL v_upper, REAL* upper_val, 
+                       int n_lower, REAL v_lower, REAL* lower_val
+                       );
+void inEvalVStrip(int n_left, REAL u_left, REAL* left_val, int n_right, REAL u_right, REAL* right_val);
+
+void inPreEvaluateBV(int k, int uorder, int vorder, REAL vprime, REAL *baseData);
+void inPreEvaluateBU(int k, int uorder, int vorder, REAL uprime, REAL *baseData);
+void inPreEvaluateBV_intfac(REAL v )
+  {
+   inPreEvaluateBV(global_ev_k, global_ev_uorder, global_ev_vorder, (v-global_ev_v1)/(global_ev_v2-global_ev_v1), global_ev_ctlPoints);
+  }
+
+void inPreEvaluateBU_intfac(REAL u)
+  {
+    inPreEvaluateBU(global_ev_k, global_ev_uorder, global_ev_vorder, (u-global_ev_u1)/(global_ev_u2-global_ev_u1), global_ev_ctlPoints); 
+  }
+
+void inDoDomain2WithDerivsBV(int k, REAL u, REAL v,
+			     REAL u1, REAL u2, int uorder,
+			     REAL v1, REAL v2, int vorder,
+			     REAL *baseData,
+			     REAL *retPoint, REAL* retdu, REAL *retdv);
+
+void inDoDomain2WithDerivsBU(int k, REAL u, REAL v,
+			     REAL u1, REAL u2, int uorder,
+			     REAL v1, REAL v2, int vorder,
+			     REAL *baseData,
+			     REAL *retPoint, REAL* retdu, REAL *retdv);
+
+
+void inDoEvalCoord2NOGE_BV(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal);
+
+void inDoEvalCoord2NOGE_BU(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal);
+
+void inBPMEval(bezierPatchMesh* bpm);
+void inBPMListEval(bezierPatchMesh* list);
+
+/*-------------begin for surfEvalMachine -------------*/
+surfEvalMachine em_vertex;
+surfEvalMachine em_normal;
+surfEvalMachine em_color;
+surfEvalMachine em_texcoord;
+
+int auto_normal_flag; //whether to output normla or not in callback
+                      //determined by GL_AUTO_NORMAL and callback_auto_normal
+int callback_auto_normal; //GLU_CALLBACK_AUTO_NORMAL_EXT
+int vertex_flag;
+int normal_flag;
+int color_flag;
+int texcoord_flag;
+
+void inMap2fEM(int which, //0:vert,1:norm,2:color,3:tex
+	       int dimension,
+	      REAL ulower,
+	      REAL uupper,
+	      int ustride,
+	      int uorder,
+	      REAL vlower,
+	      REAL vupper,
+	      int vstride,
+	      int vorder,
+	      REAL *ctlPoints);
+
+void inDoDomain2WithDerivsEM(surfEvalMachine *em, REAL u, REAL v, 
+				REAL *retPoint, REAL *retdu, REAL *retdv);
+void inDoDomain2EM(surfEvalMachine *em, REAL u, REAL v, 
+				REAL *retPoint);
+ void inDoEvalCoord2EM(REAL u, REAL v);
+
+void inBPMEvalEM(bezierPatchMesh* bpm);
+void inBPMListEvalEM(bezierPatchMesh* list);
+
+/*-------------end for surfEvalMachine -------------*/
+
+
+   /*************end for internal evaluators*****************/
+		       
+};
+
+inline void StoredVertex::invoke(OpenGLSurfaceEvaluator *eval)
+{
+    switch(type) {
+      case TYPECOORD:
+	eval->coord2f(coord[0], coord[1]);
+	break;
+      case TYPEPOINT:
+	eval->point2i(point[0], point[1]);
+	break;
+      default:
+	break;
+    }
+}
+
+#endif /* __gluglsurfeval_h_ */
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/incurveeval.cc b/mesalib/src/glu/sgi/libnurbs/interface/incurveeval.cc
new file mode 100644
index 000000000..96ea8896a
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/incurveeval.cc
@@ -0,0 +1,206 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+** 
+** http://oss.sgi.com/projects/FreeB
+** 
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+** 
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+** 
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include "glcurveval.h"
+
+
+/*
+ *compute the Bezier polynomials C[n,j](v) for all j at v with 
+ *return values stored in coeff[], where 
+ *  C[n,j](v) = (n,j) * v^j * (1-v)^(n-j),
+ *  j=0,1,2,...,n.
+ *order : n+1
+ *vprime: v
+ *coeff : coeff[j]=C[n,j](v), this array store the returned values.
+ *The algorithm is a recursive scheme:
+ *   C[0,0]=1;
+ *   C[n,j](v) = (1-v)*C[n-1,j](v) + v*C[n-1,j-1](v), n>=1
+ *This code is copied from opengl/soft/so_eval.c:PreEvaluate
+ */
+void OpenGLCurveEvaluator::inPreEvaluate(int order, REAL vprime, REAL *coeff)
+{
+  int i, j;
+  REAL oldval, temp;
+  REAL oneMinusvprime;
+  
+  /*
+   * Minor optimization
+   * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to
+     * their i==1 loop values to avoid the initialization and the i==1 loop.
+     */
+  if (order == 1) {
+    coeff[0] = 1.0;
+    return;
+  }
+  
+  oneMinusvprime = 1-vprime;
+  coeff[0] = oneMinusvprime;
+  coeff[1] = vprime;
+  if (order == 2) return;
+  
+  for (i = 2; i < order; i++) {
+    oldval = coeff[0] * vprime;
+    coeff[0] = oneMinusvprime * coeff[0];
+    for (j = 1; j < i; j++) {
+      temp = oldval;
+      oldval = coeff[j] * vprime;
+	    coeff[j] = temp + oneMinusvprime * coeff[j];
+    }
+    coeff[j] = oldval;
+  }
+}
+
+void OpenGLCurveEvaluator::inMap1f(int which, //0: vert, 1: norm, 2: color, 3: tex
+				   int k, //dimension
+				   REAL ulower,
+				   REAL uupper,
+				   int ustride,
+				   int uorder,
+				   REAL *ctlpoints)
+{
+  int i,x;
+  curveEvalMachine *temp_em;
+  switch(which){
+  case 0: //vertex
+    vertex_flag = 1;
+    temp_em = &em_vertex;
+    break;
+  case 1: //normal
+    normal_flag = 1;
+    temp_em = &em_normal;
+    break;
+  case 2: //color
+    color_flag = 1;
+    temp_em = &em_color;
+    break;
+  default:
+    texcoord_flag = 1;
+    temp_em = &em_texcoord;
+    break;
+  }
+  
+  REAL *data = temp_em->ctlpoints;
+  temp_em->uprime = -1; //initialized
+  temp_em->k = k;
+  temp_em->u1 = ulower;
+  temp_em->u2 = uupper;
+  temp_em->ustride = ustride;
+  temp_em->uorder = uorder;
+  /*copy the control points*/
+  for(i=0; i<uorder; i++){
+    for(x=0; x<k; x++){
+      data[x] = ctlpoints[x];
+    }
+    ctlpoints += ustride;
+    data += k;
+  }     
+}
+
+void OpenGLCurveEvaluator::inDoDomain1(curveEvalMachine *em, REAL u, REAL *retPoint)
+{
+  int j, row;
+  REAL the_uprime;
+  REAL *data;
+  
+  if(em->u2 == em->u1)
+    return;
+  the_uprime = (u-em->u1) / (em->u2-em->u1);
+  /*use already cached values if possible*/
+  if(em->uprime != the_uprime){
+    inPreEvaluate(em->uorder, the_uprime, em->ucoeff);
+    em->uprime = the_uprime;
+  }
+  
+  for(j=0; j<em->k; j++){
+    data = em->ctlpoints+j;
+    retPoint[j] = 0.0;
+    for(row=0; row<em->uorder; row++)
+      {
+	retPoint[j] += em->ucoeff[row] * (*data);
+	data += em->k;
+      }
+  } 
+}
+
+void  OpenGLCurveEvaluator::inDoEvalCoord1(REAL u)
+{
+  REAL temp_vertex[4];
+  REAL temp_normal[3];
+  REAL temp_color[4];
+  REAL temp_texcoord[4];
+  if(texcoord_flag) //there is a texture map
+    {
+      inDoDomain1(&em_texcoord, u, temp_texcoord);
+      texcoordCallBack(temp_texcoord, userData);
+    }
+#ifdef DEBUG
+printf("color_flag = %i\n", color_flag);
+#endif
+  if(color_flag) //there is a color map
+    {
+      inDoDomain1(&em_color, u, temp_color);
+      colorCallBack(temp_color, userData);
+    }
+  if(normal_flag) //there is a normal map
+    {
+      inDoDomain1(&em_normal, u, temp_normal);
+      normalCallBack(temp_normal, userData);
+    }
+  if(vertex_flag)
+    {
+      inDoDomain1(&em_vertex, u, temp_vertex);
+      vertexCallBack(temp_vertex, userData);
+    }
+}
+
+void OpenGLCurveEvaluator::inMapMesh1f(int umin, int umax)
+{
+  REAL du, u;
+  int i;
+  if(global_grid_nu == 0)
+    return; //no points to output
+  du = (global_grid_u1 - global_grid_u0) / (REAL) global_grid_nu;
+  bgnline();
+  for(i=umin; i<= umax; i++){
+    u = (i==global_grid_nu)? global_grid_u1: global_grid_u0 + i*du;
+    inDoEvalCoord1(u);
+  }
+  endline();
+}
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/insurfeval.cc b/mesalib/src/glu/sgi/libnurbs/interface/insurfeval.cc
new file mode 100644
index 000000000..9d0c82a91
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/insurfeval.cc
@@ -0,0 +1,2064 @@
+/*
+** License Applicability. Except to the extent portions of this file are
+** made subject to an alternative license as permitted in the SGI Free
+** Software License B, Version 1.1 (the "License"), the contents of this
+** file are subject only to the provisions of the License. You may not use
+** this file except in compliance with the License. You may obtain a copy
+** of the License at Silicon Graphics, Inc., attn: Legal Services, 1600
+** Amphitheatre Parkway, Mountain View, CA 94043-1351, or at:
+** 
+** http://oss.sgi.com/projects/FreeB
+** 
+** Note that, as provided in the License, the Software is distributed on an
+** "AS IS" basis, with ALL EXPRESS AND IMPLIED WARRANTIES AND CONDITIONS
+** DISCLAIMED, INCLUDING, WITHOUT LIMITATION, ANY IMPLIED WARRANTIES AND
+** CONDITIONS OF MERCHANTABILITY, SATISFACTORY QUALITY, FITNESS FOR A
+** PARTICULAR PURPOSE, AND NON-INFRINGEMENT.
+** 
+** Original Code. The Original Code is: OpenGL Sample Implementation,
+** Version 1.2.1, released January 26, 2000, developed by Silicon Graphics,
+** Inc. The Original Code is Copyright (c) 1991-2000 Silicon Graphics, Inc.
+** Copyright in any portions created by third parties is as indicated
+** elsewhere herein. All Rights Reserved.
+** 
+** Additional Notice Provisions: The application programming interfaces
+** established by SGI in conjunction with the Original Code are The
+** OpenGL(R) Graphics System: A Specification (Version 1.2.1), released
+** April 1, 1999; The OpenGL(R) Graphics System Utility Library (Version
+** 1.3), released November 4, 1998; and OpenGL(R) Graphics with the X
+** Window System(R) (Version 1.3), released October 19, 1998. This software
+** was created using the OpenGL(R) version 1.2.1 Sample Implementation
+** published by SGI, but has not been independently verified as being
+** compliant with the OpenGL(R) version 1.2.1 Specification.
+**
+*/
+/*
+*/
+
+#include "gluos.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <GL/gl.h>
+#include <math.h>
+#include <assert.h>
+
+#include "glsurfeval.h"
+
+//extern int surfcount;
+
+//#define CRACK_TEST
+
+#define AVOID_ZERO_NORMAL
+
+#ifdef AVOID_ZERO_NORMAL
+#define myabs(x)  ((x>0)? x: (-x))
+#define MYZERO 0.000001
+#define MYDELTA 0.001
+#endif
+
+//#define USE_LOD
+#ifdef USE_LOD
+//#define LOD_EVAL_COORD(u,v) inDoEvalCoord2EM(u,v)
+#define LOD_EVAL_COORD(u,v) glEvalCoord2f(u,v)
+
+static void LOD_interpolate(REAL A[2], REAL B[2], REAL C[2], int j, int k, int pow2_level,
+			    REAL& u, REAL& v)
+{
+  REAL a,a1,b,b1;
+
+  a = ((REAL) j) / ((REAL) pow2_level);
+  a1 = 1-a;
+
+  if(j != 0)
+    {
+      b = ((REAL) k) / ((REAL)j);
+      b1 = 1-b;
+    }
+  REAL x,y,z;
+  x = a1;
+  if(j==0)
+    {
+      y=0; z=0;
+    }
+  else{
+    y = b1*a;
+    z = b *a;
+  }
+
+  u = x*A[0] + y*B[0] + z*C[0];
+  v = x*A[1] + y*B[1] + z*C[1];
+}
+
+void OpenGLSurfaceEvaluator::LOD_triangle(REAL A[2], REAL B[2], REAL C[2], 
+			 int level)     
+{
+  int k,j;
+  int pow2_level;
+  /*compute 2^level*/
+  pow2_level = 1;
+
+  for(j=0; j<level; j++)
+    pow2_level *= 2;
+  for(j=0; j<=pow2_level-1; j++)
+    {
+      REAL u,v;
+
+/*      beginCallBack(GL_TRIANGLE_STRIP);*/
+glBegin(GL_TRIANGLE_STRIP);
+      LOD_interpolate(A,B,C, j+1, j+1, pow2_level, u,v);
+#ifdef USE_LOD
+      LOD_EVAL_COORD(u,v);
+//      glEvalCoord2f(u,v);
+#else
+      inDoEvalCoord2EM(u,v);
+#endif
+
+      for(k=0; k<=j; k++)
+	{
+	  LOD_interpolate(A,B,C,j,j-k,pow2_level, u,v);
+#ifdef USE_LOD
+          LOD_EVAL_COORD(u,v);
+//	  glEvalCoord2f(u,v);
+#else
+	  inDoEvalCoord2EM(u,v);
+#endif
+
+	  LOD_interpolate(A,B,C,j+1,j-k,pow2_level, u,v);
+
+#ifdef USE_LOD
+	  LOD_EVAL_COORD(u,v);
+//	  glEvalCoord2f(u,v);
+#else
+	  inDoEvalCoord2EM(u,v);
+#endif
+	}
+//      endCallBack();	
+glEnd();
+    }
+}
+
+void OpenGLSurfaceEvaluator::LOD_eval(int num_vert, REAL* verts, int type,
+		     int level
+		     )
+{
+  int i,k;
+  switch(type){
+  case GL_TRIANGLE_STRIP:
+  case GL_QUAD_STRIP:
+    for(i=2, k=4; i<=num_vert-2; i+=2, k+=4)
+      {
+	LOD_triangle(verts+k-4, verts+k-2, verts+k,
+		     level
+		     );
+	LOD_triangle(verts+k-2, verts+k+2, verts+k,
+		     level
+		     );
+      }
+    if(num_vert % 2 ==1) 
+      {
+	LOD_triangle(verts+2*(num_vert-3), verts+2*(num_vert-2), verts+2*(num_vert-1),
+		     level
+		     );
+      }
+    break;      
+  case GL_TRIANGLE_FAN:
+    for(i=1, k=2; i<=num_vert-2; i++, k+=2)
+      {
+	LOD_triangle(verts,verts+k, verts+k+2,
+		     level
+		     );
+      }
+    break;
+  
+  default:
+    fprintf(stderr, "typy not supported in LOD_\n");
+  }
+}
+	
+
+#endif //USE_LOD
+
+//#define  GENERIC_TEST
+#ifdef GENERIC_TEST
+extern float xmin, xmax, ymin, ymax, zmin, zmax; /*bounding box*/
+extern int temp_signal;
+
+static void gTessVertexSphere(float u, float v, float temp_normal[3], float temp_vertex[3])
+{
+  float r=2.0;
+  float Ox = 0.5*(xmin+xmax);
+  float Oy = 0.5*(ymin+ymax);
+  float Oz = 0.5*(zmin+zmax);
+  float nx = cos(v) * sin(u);
+  float ny = sin(v) * sin(u);
+  float nz = cos(u);
+  float x= Ox+r * nx;
+  float y= Oy+r * ny;
+  float z= Oz+r * nz;
+
+  temp_normal[0] = nx;
+  temp_normal[1] = ny;
+  temp_normal[2] =  nz;
+  temp_vertex[0] = x;
+  temp_vertex[1] = y;
+  temp_vertex[2] = z;
+
+//  glNormal3f(nx,ny,nz);
+//  glVertex3f(x,y,z);
+}
+
+static void gTessVertexCyl(float u, float v, float temp_normal[3], float temp_vertex[3])
+{
+   float r=2.0;
+  float Ox = 0.5*(xmin+xmax);
+  float Oy = 0.5*(ymin+ymax);
+  float Oz = 0.5*(zmin+zmax);
+  float nx = cos(v);
+  float ny = sin(v);
+  float nz = 0;
+  float x= Ox+r * nx;
+  float y= Oy+r * ny;
+  float z= Oz - 2*u;
+
+  temp_normal[0] = nx;
+  temp_normal[1] = ny;
+  temp_normal[2] =  nz;
+  temp_vertex[0] = x;
+  temp_vertex[1] = y;
+  temp_vertex[2] = z;
+
+/*  
+  glNormal3f(nx,ny,nz);
+  glVertex3f(x,y,z);
+*/
+}
+
+#endif //GENERIC_TEST
+
+void OpenGLSurfaceEvaluator::inBPMListEval(bezierPatchMesh* list)
+{
+  bezierPatchMesh* temp;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      inBPMEval(temp);
+    }
+}
+
+void OpenGLSurfaceEvaluator::inBPMEval(bezierPatchMesh* bpm)
+{
+  int i,j,k,l;
+  float u,v;
+
+  int ustride = bpm->bpatch->dimension * bpm->bpatch->vorder;
+  int vstride = bpm->bpatch->dimension;
+  inMap2f( 
+	  (bpm->bpatch->dimension == 3)? GL_MAP2_VERTEX_3 : GL_MAP2_VERTEX_4,
+	  bpm->bpatch->umin,
+	  bpm->bpatch->umax,
+	  ustride,
+	  bpm->bpatch->uorder,
+	  bpm->bpatch->vmin,
+	  bpm->bpatch->vmax,
+	  vstride,
+	  bpm->bpatch->vorder,
+	  bpm->bpatch->ctlpoints);
+  
+  bpm->vertex_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3+1); /*in case the origional dimenion is 4, then we need 4 space to pass to evaluator.*/
+  assert(bpm->vertex_array);
+  bpm->normal_array = (float*) malloc(sizeof(float)* (bpm->index_UVarray/2) * 3);
+  assert(bpm->normal_array);
+#ifdef CRACK_TEST
+if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
+  && global_ev_v1 ==2 &&   global_ev_v2 == 3)
+{
+REAL vertex[4];
+REAL normal[4];
+#ifdef DEBUG
+printf("***number 1\n");
+#endif
+
+beginCallBack(GL_QUAD_STRIP, NULL);
+inEvalCoord2f(3.0, 3.0);
+inEvalCoord2f(2.0, 3.0);
+inEvalCoord2f(3.0, 2.7);
+inEvalCoord2f(2.0, 2.7);
+inEvalCoord2f(3.0, 2.0);
+inEvalCoord2f(2.0, 2.0);
+endCallBack(NULL);
+
+
+beginCallBack(GL_TRIANGLE_STRIP, NULL);
+inEvalCoord2f(2.0, 3.0);
+inEvalCoord2f(2.0, 2.0);
+inEvalCoord2f(2.0, 2.7);
+endCallBack(NULL);
+
+}
+
+/*
+if(  global_ev_u1 ==2 &&   global_ev_u2 == 3
+  && global_ev_v1 ==1 &&   global_ev_v2 == 2)
+{
+#ifdef DEBUG
+printf("***number 2\n");
+#endif
+beginCallBack(GL_QUAD_STRIP);
+inEvalCoord2f(2.0, 2.0);
+inEvalCoord2f(2.0, 1.0);
+inEvalCoord2f(3.0, 2.0);
+inEvalCoord2f(3.0, 1.0);
+endCallBack();
+}
+*/
+if(  global_ev_u1 ==1 &&   global_ev_u2 == 2
+  && global_ev_v1 ==2 &&   global_ev_v2 == 3)
+{
+#ifdef DEBUG
+printf("***number 3\n");
+#endif
+beginCallBack(GL_QUAD_STRIP, NULL);
+inEvalCoord2f(2.0, 3.0);
+inEvalCoord2f(1.0, 3.0);
+inEvalCoord2f(2.0, 2.3);
+inEvalCoord2f(1.0, 2.3);
+inEvalCoord2f(2.0, 2.0);
+inEvalCoord2f(1.0, 2.0);
+endCallBack(NULL);
+
+beginCallBack(GL_TRIANGLE_STRIP, NULL);
+inEvalCoord2f(2.0, 2.3);
+inEvalCoord2f(2.0, 2.0);
+inEvalCoord2f(2.0, 3.0);
+endCallBack(NULL);
+
+}
+return;
+#endif
+
+  k=0;
+  l=0;
+
+  for(i=0; i<bpm->index_length_array; i++)
+    {
+      beginCallBack(bpm->type_array[i], userData);
+      for(j=0; j<bpm->length_array[i]; j++)
+	{
+	  u = bpm->UVarray[k];
+	  v = bpm->UVarray[k+1];
+	  inDoEvalCoord2NOGE(u,v,
+			     bpm->vertex_array+l,
+			     bpm->normal_array+l);
+
+	  normalCallBack(bpm->normal_array+l, userData);
+	  vertexCallBack(bpm->vertex_array+l, userData);
+
+	  k += 2;
+	  l += 3;
+	}
+      endCallBack(userData);
+    }
+}
+
+void OpenGLSurfaceEvaluator::inEvalPoint2(int i, int j)
+{
+  REAL du, dv;
+  REAL point[4];
+  REAL normal[3];
+  REAL u,v;
+  du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
+  dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;
+  u = (i==global_grid_nu)? global_grid_u1:(global_grid_u0 + i*du);
+  v = (j == global_grid_nv)? global_grid_v1: (global_grid_v0 +j*dv);
+  inDoEvalCoord2(u,v,point,normal);
+}
+
+void OpenGLSurfaceEvaluator::inEvalCoord2f(REAL u, REAL v)
+{
+
+  REAL point[4];
+  REAL normal[3];
+  inDoEvalCoord2(u,v,point, normal);
+}
+
+
+
+/*define a grid. store the values into the global variabls:
+ * global_grid_*
+ *These values will be used later by evaluating functions
+ */
+void OpenGLSurfaceEvaluator::inMapGrid2f(int nu, REAL u0, REAL u1,
+		 int nv, REAL v0, REAL v1)
+{
+ global_grid_u0 = u0;
+ global_grid_u1 = u1;
+ global_grid_nu = nu;
+ global_grid_v0 = v0;
+ global_grid_v1 = v1;
+ global_grid_nv = nv;
+}
+
+void OpenGLSurfaceEvaluator::inEvalMesh2(int lowU, int lowV, int highU, int highV)
+{
+  REAL du, dv;
+  int i,j;
+  REAL point[4];
+  REAL normal[3];
+  if(global_grid_nu == 0 || global_grid_nv == 0)
+    return; /*no points need to be output*/
+  du = (global_grid_u1 - global_grid_u0) / (REAL)global_grid_nu;
+  dv = (global_grid_v1 - global_grid_v0) / (REAL)global_grid_nv;  
+  
+  if(global_grid_nu >= global_grid_nv){
+    for(i=lowU; i<highU; i++){
+      REAL u1 = (i==global_grid_nu)? global_grid_u1:(global_grid_u0 + i*du);
+      REAL u2 = ((i+1) == global_grid_nu)? global_grid_u1: (global_grid_u0+(i+1)*du);
+      
+      bgnqstrip();
+      for(j=highV; j>=lowV; j--){
+	REAL v1 = (j == global_grid_nv)? global_grid_v1: (global_grid_v0 +j*dv);
+	
+	inDoEvalCoord2(u1, v1, point, normal);
+	inDoEvalCoord2(u2, v1, point, normal);
+      }
+      endqstrip();
+    }
+  }
+  
+  else{
+    for(i=lowV; i<highV; i++){
+      REAL v1 = (i==global_grid_nv)? global_grid_v1:(global_grid_v0 + i*dv);
+      REAL v2 = ((i+1) == global_grid_nv)? global_grid_v1: (global_grid_v0+(i+1)*dv);
+      
+      bgnqstrip();
+      for(j=highU; j>=lowU; j--){
+	REAL u1 = (j == global_grid_nu)? global_grid_u1: (global_grid_u0 +j*du);	
+	inDoEvalCoord2(u1, v2, point, normal);
+	inDoEvalCoord2(u1, v1, point, normal);
+      }
+      endqstrip();
+    }
+  }
+    
+}
+
+void OpenGLSurfaceEvaluator::inMap2f(int k,
+	     REAL ulower,
+	     REAL uupper,
+	     int ustride,
+	     int uorder,
+	     REAL vlower,
+	     REAL vupper,
+	     int vstride,
+	     int vorder,
+	     REAL *ctlPoints)
+{
+  int i,j,x;
+  REAL *data = global_ev_ctlPoints;
+  
+
+
+  if(k == GL_MAP2_VERTEX_3) k=3;
+  else if (k==GL_MAP2_VERTEX_4) k =4;
+  else {
+    printf("error in inMap2f, maptype=%i is wrong, k,map is not updated\n", k);
+    return;
+  }
+  
+  global_ev_k = k;
+  global_ev_u1 = ulower;
+  global_ev_u2 = uupper;
+  global_ev_ustride = ustride;
+  global_ev_uorder = uorder;
+  global_ev_v1 = vlower;
+  global_ev_v2 = vupper;
+  global_ev_vstride = vstride;
+  global_ev_vorder = vorder;
+
+  /*copy the contrl points from ctlPoints to global_ev_ctlPoints*/
+  for (i=0; i<uorder; i++) {
+    for (j=0; j<vorder; j++) {
+      for (x=0; x<k; x++) {
+	data[x] = ctlPoints[x];
+      }
+      ctlPoints += vstride;
+      data += k;
+    }
+    ctlPoints += ustride - vstride * vorder;
+  }
+
+}
+
+
+/*
+ *given a point p with homegeneous coordiante (x,y,z,w), 
+ *let pu(x,y,z,w) be its partial derivative vector with
+ *respect to u
+ *and pv(x,y,z,w) be its partial derivative vector with repect to v.
+ *This function returns the partial derivative vectors of the
+ *inhomegensous coordinates, i.e., 
+ * (x/w, y/w, z/w) with respect to u and v.
+ */
+void OpenGLSurfaceEvaluator::inComputeFirstPartials(REAL *p, REAL *pu, REAL *pv)
+{
+    pu[0] = pu[0]*p[3] - pu[3]*p[0];
+    pu[1] = pu[1]*p[3] - pu[3]*p[1];
+    pu[2] = pu[2]*p[3] - pu[3]*p[2];
+
+    pv[0] = pv[0]*p[3] - pv[3]*p[0];
+    pv[1] = pv[1]*p[3] - pv[3]*p[1];
+    pv[2] = pv[2]*p[3] - pv[3]*p[2];
+}
+
+/*compute the cross product of pu and pv and normalize.
+ *the normal is returned in retNormal
+ * pu: dimension 3
+ * pv: dimension 3
+ * n: return normal, of dimension 3
+ */
+void OpenGLSurfaceEvaluator::inComputeNormal2(REAL *pu, REAL *pv, REAL *n)
+{
+  REAL mag; 
+
+  n[0] = pu[1]*pv[2] - pu[2]*pv[1];
+  n[1] = pu[2]*pv[0] - pu[0]*pv[2];
+  n[2] = pu[0]*pv[1] - pu[1]*pv[0];  
+
+  mag = sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]);
+
+  if (mag > 0.0) {
+     n[0] /= mag; 
+     n[1] /= mag;
+     n[2] /= mag;
+  }
+}
+ 
+
+
+/*Compute point and normal
+ *see the head of inDoDomain2WithDerivs
+ *for the meaning of the arguments
+ */
+void OpenGLSurfaceEvaluator::inDoEvalCoord2(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal)
+{
+
+  REAL du[4];
+  REAL dv[4];
+
+ 
+  assert(global_ev_k>=3 && global_ev_k <= 4);
+  /*compute homegeneous point and partial derivatives*/
+  inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
+
+#ifdef AVOID_ZERO_NORMAL
+
+  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
+    {
+
+      REAL tempdu[4];
+      REAL tempdata[4];
+      REAL u1 = global_ev_u1;
+      REAL u2 = global_ev_u2;
+      if(u-MYDELTA*(u2-u1) < u1)
+	u = u+ MYDELTA*(u2-u1);
+      else
+	u = u-MYDELTA*(u2-u1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
+    }
+  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
+    {
+      REAL tempdv[4];
+      REAL tempdata[4];
+      REAL v1 = global_ev_v1;
+      REAL v2 = global_ev_v2;
+      if(v-MYDELTA*(v2-v1) < v1)
+	v = v+ MYDELTA*(v2-v1);
+      else
+	v = v-MYDELTA*(v2-v1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
+    }
+#endif
+
+
+  /*compute normal*/
+  switch(global_ev_k){
+  case 3:
+    inComputeNormal2(du, dv, retNormal);
+
+    break;
+  case 4:
+    inComputeFirstPartials(retPoint, du, dv);
+    inComputeNormal2(du, dv, retNormal);
+    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
+    retPoint[0] /= retPoint[3];
+    retPoint[1] /= retPoint[3];
+    retPoint[2] /= retPoint[3];
+    break;
+  }
+  /*output this vertex*/
+/*  inMeshStreamInsert(global_ms, retPoint, retNormal);*/
+
+
+
+  glNormal3fv(retNormal);
+  glVertex3fv(retPoint);
+
+
+
+
+  #ifdef DEBUG
+  printf("vertex(%f,%f,%f)\n", retPoint[0],retPoint[1],retPoint[2]);
+  #endif
+  
+
+
+}
+
+/*Compute point and normal
+ *see the head of inDoDomain2WithDerivs
+ *for the meaning of the arguments
+ */
+void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BU(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal)
+{
+
+  REAL du[4];
+  REAL dv[4];
+
+ 
+  assert(global_ev_k>=3 && global_ev_k <= 4);
+  /*compute homegeneous point and partial derivatives*/
+//   inPreEvaluateBU(global_ev_k, global_ev_uorder, global_ev_vorder, (u-global_ev_u1)/(global_ev_u2-global_ev_u1), global_ev_ctlPoints);
+  inDoDomain2WithDerivsBU(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
+
+
+#ifdef AVOID_ZERO_NORMAL
+
+  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
+    {
+
+      REAL tempdu[4];
+      REAL tempdata[4];
+      REAL u1 = global_ev_u1;
+      REAL u2 = global_ev_u2;
+      if(u-MYDELTA*(u2-u1) < u1)
+	u = u+ MYDELTA*(u2-u1);
+      else
+	u = u-MYDELTA*(u2-u1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
+    }
+  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
+    {
+      REAL tempdv[4];
+      REAL tempdata[4];
+      REAL v1 = global_ev_v1;
+      REAL v2 = global_ev_v2;
+      if(v-MYDELTA*(v2-v1) < v1)
+	v = v+ MYDELTA*(v2-v1);
+      else
+	v = v-MYDELTA*(v2-v1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
+    }
+#endif
+
+  /*compute normal*/
+  switch(global_ev_k){
+  case 3:
+    inComputeNormal2(du, dv, retNormal);
+    break;
+  case 4:
+    inComputeFirstPartials(retPoint, du, dv);
+    inComputeNormal2(du, dv, retNormal);
+    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
+    retPoint[0] /= retPoint[3];
+    retPoint[1] /= retPoint[3];
+    retPoint[2] /= retPoint[3];
+    break;
+  }
+}
+
+/*Compute point and normal
+ *see the head of inDoDomain2WithDerivs
+ *for the meaning of the arguments
+ */
+void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE_BV(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal)
+{
+
+  REAL du[4];
+  REAL dv[4];
+
+ 
+  assert(global_ev_k>=3 && global_ev_k <= 4);
+  /*compute homegeneous point and partial derivatives*/
+//   inPreEvaluateBV(global_ev_k, global_ev_uorder, global_ev_vorder, (v-global_ev_v1)/(global_ev_v2-global_ev_v1), global_ev_ctlPoints);
+
+  inDoDomain2WithDerivsBV(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
+
+
+#ifdef AVOID_ZERO_NORMAL
+
+  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
+    {
+
+      REAL tempdu[4];
+      REAL tempdata[4];
+      REAL u1 = global_ev_u1;
+      REAL u2 = global_ev_u2;
+      if(u-MYDELTA*(u2-u1) < u1)
+	u = u+ MYDELTA*(u2-u1);
+      else
+	u = u-MYDELTA*(u2-u1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
+    }
+  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
+    {
+      REAL tempdv[4];
+      REAL tempdata[4];
+      REAL v1 = global_ev_v1;
+      REAL v2 = global_ev_v2;
+      if(v-MYDELTA*(v2-v1) < v1)
+	v = v+ MYDELTA*(v2-v1);
+      else
+	v = v-MYDELTA*(v2-v1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
+    }
+#endif
+
+  /*compute normal*/
+  switch(global_ev_k){
+  case 3:
+    inComputeNormal2(du, dv, retNormal);
+    break;
+  case 4:
+    inComputeFirstPartials(retPoint, du, dv);
+    inComputeNormal2(du, dv, retNormal);
+    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
+    retPoint[0] /= retPoint[3];
+    retPoint[1] /= retPoint[3];
+    retPoint[2] /= retPoint[3];
+    break;
+  }
+}
+ 
+
+/*Compute point and normal
+ *see the head of inDoDomain2WithDerivs
+ *for the meaning of the arguments
+ */
+void OpenGLSurfaceEvaluator::inDoEvalCoord2NOGE(REAL u, REAL v,
+			   REAL *retPoint, REAL *retNormal)
+{
+
+  REAL du[4];
+  REAL dv[4];
+
+ 
+  assert(global_ev_k>=3 && global_ev_k <= 4);
+  /*compute homegeneous point and partial derivatives*/
+  inDoDomain2WithDerivs(global_ev_k, u, v, global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, retPoint, du, dv);
+
+
+#ifdef AVOID_ZERO_NORMAL
+
+  if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
+    {
+
+      REAL tempdu[4];
+      REAL tempdata[4];
+      REAL u1 = global_ev_u1;
+      REAL u2 = global_ev_u2;
+      if(u-MYDELTA*(u2-u1) < u1)
+	u = u+ MYDELTA*(u2-u1);
+      else
+	u = u-MYDELTA*(u2-u1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, tempdu, dv);
+    }
+  if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
+    {
+      REAL tempdv[4];
+      REAL tempdata[4];
+      REAL v1 = global_ev_v1;
+      REAL v2 = global_ev_v2;
+      if(v-MYDELTA*(v2-v1) < v1)
+	v = v+ MYDELTA*(v2-v1);
+      else
+	v = v-MYDELTA*(v2-v1);
+      inDoDomain2WithDerivs(global_ev_k, u,v,global_ev_u1, global_ev_u2, global_ev_uorder, global_ev_v1, global_ev_v2, global_ev_vorder, global_ev_ctlPoints, tempdata, du, tempdv);
+    }
+#endif
+
+  /*compute normal*/
+  switch(global_ev_k){
+  case 3:
+    inComputeNormal2(du, dv, retNormal);
+    break;
+  case 4:
+    inComputeFirstPartials(retPoint, du, dv);
+    inComputeNormal2(du, dv, retNormal);
+    /*transform the homegeneous coordinate of retPoint into inhomogenous one*/
+    retPoint[0] /= retPoint[3];
+    retPoint[1] /= retPoint[3];
+    retPoint[2] /= retPoint[3];
+    break;
+  }
+//  glNormal3fv(retNormal);
+//  glVertex3fv(retPoint);
+}
+ 
+void OpenGLSurfaceEvaluator::inPreEvaluateBV(int k, int uorder, int vorder, REAL vprime, REAL *baseData)
+{
+  int j,row,col;
+  REAL p, pdv;
+  REAL *data;
+
+  if(global_vprime != vprime || global_vorder != vorder) {      
+    inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
+    global_vprime = vprime;
+    global_vorder = vorder;
+  }
+
+  for(j=0; j<k; j++){
+    data = baseData+j;
+    for(row=0; row<uorder; row++){
+      p = global_vcoeff[0] * (*data);
+      pdv = global_vcoeffDeriv[0] * (*data);
+      data += k;
+      for(col = 1; col < vorder; col++){
+	p += global_vcoeff[col] *  (*data);
+	pdv += global_vcoeffDeriv[col] * (*data);
+	data += k;
+      }
+      global_BV[row][j]  = p;
+      global_PBV[row][j]  = pdv;
+    }
+  }
+}
+
+void OpenGLSurfaceEvaluator::inPreEvaluateBU(int k, int uorder, int vorder, REAL uprime, REAL *baseData)
+{
+  int j,row,col;
+  REAL p, pdu;
+  REAL *data;
+
+  if(global_uprime != uprime || global_uorder != uorder) {      
+    inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
+    global_uprime = uprime;
+    global_uorder = uorder;
+  }
+
+  for(j=0; j<k; j++){
+    data = baseData+j;
+    for(col=0; col<vorder; col++){
+      data = baseData+j + k*col;
+      p = global_ucoeff[0] * (*data);
+      pdu = global_ucoeffDeriv[0] * (*data);
+      data += k*uorder;
+      for(row = 1; row < uorder; row++){
+	p += global_ucoeff[row] *  (*data);
+	pdu += global_ucoeffDeriv[row] * (*data);
+	data += k * uorder;
+      }
+      global_BU[col][j]  = p;
+      global_PBU[col][j]  = pdu;
+    }
+  }
+}
+ 
+void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBU(int k, REAL u, REAL v,
+						      REAL u1, REAL u2, int uorder,
+						      REAL v1, REAL v2, int vorder,
+						      REAL *baseData,
+						      REAL *retPoint, REAL* retdu, REAL *retdv)
+{
+  int j, col;
+
+  REAL vprime;
+
+
+  if((u2 == u1) || (v2 == v1))
+    return;
+
+  vprime = (v - v1) / (v2 - v1);
+
+
+  if(global_vprime != vprime || global_vorder != vorder) {
+    inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
+    global_vprime = vprime;
+    global_vorder = vorder;
+  }
+
+
+  for(j=0; j<k; j++)
+    {
+      retPoint[j] = retdu[j] = retdv[j] = 0.0;
+      for (col = 0; col < vorder; col++)  {
+	retPoint[j] += global_BU[col][j] * global_vcoeff[col];
+	retdu[j] += global_PBU[col][j] * global_vcoeff[col];
+	retdv[j] += global_BU[col][j] * global_vcoeffDeriv[col];
+      }
+    }
+}    
+   
+void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsBV(int k, REAL u, REAL v,
+						      REAL u1, REAL u2, int uorder,
+						      REAL v1, REAL v2, int vorder,
+						      REAL *baseData,
+						      REAL *retPoint, REAL* retdu, REAL *retdv)
+{
+  int j, row;
+  REAL uprime;
+
+
+  if((u2 == u1) || (v2 == v1))
+    return;
+  uprime = (u - u1) / (u2 - u1);
+
+
+  if(global_uprime != uprime || global_uorder != uorder) {
+    inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
+    global_uprime = uprime;
+    global_uorder = uorder;
+  }
+
+
+  for(j=0; j<k; j++)
+    {
+      retPoint[j] = retdu[j] = retdv[j] = 0.0;
+      for (row = 0; row < uorder; row++)  {
+	retPoint[j] += global_BV[row][j] * global_ucoeff[row];
+	retdu[j] += global_BV[row][j] * global_ucoeffDeriv[row];
+	retdv[j] += global_PBV[row][j] * global_ucoeff[row];
+      }
+    }
+}
+  
+
+/*
+ *given a Bezier surface, and parameter (u,v), compute the point in the object space,
+ *and the normal
+ *k: the dimension of the object space: usually 2,3,or 4.
+ *u,v: the paramter pair.
+ *u1,u2,uorder: the Bezier polynomial of u coord is defined on [u1,u2] with order uorder.
+ *v1,v2,vorder: the Bezier polynomial of v coord is defined on [v1,v2] with order vorder.
+ *baseData: contrl points. arranged as: (u,v,k).
+ *retPoint:  the computed point (one point) with dimension k.
+ *retdu: the computed partial derivative with respect to u.
+ *retdv: the computed partial derivative with respect to v.
+ */
+void OpenGLSurfaceEvaluator::inDoDomain2WithDerivs(int k, REAL u, REAL v, 
+				REAL u1, REAL u2, int uorder, 
+				REAL v1,  REAL v2, int vorder, 
+				REAL *baseData,
+				REAL *retPoint, REAL *retdu, REAL *retdv)
+{
+    int j, row, col;
+    REAL uprime;
+    REAL vprime;
+    REAL p;
+    REAL pdv;
+    REAL *data;
+
+    if((u2 == u1) || (v2 == v1))
+	return;
+    uprime = (u - u1) / (u2 - u1);
+    vprime = (v - v1) / (v2 - v1);
+    
+    /* Compute coefficients for values and derivs */
+
+    /* Use already cached values if possible */
+    if(global_uprime != uprime || global_uorder != uorder) {
+        inPreEvaluateWithDeriv(uorder, uprime, global_ucoeff, global_ucoeffDeriv);
+	global_uorder = uorder;
+	global_uprime = uprime;
+    }
+    if (global_vprime != vprime || 
+	  global_vorder != vorder) {
+	inPreEvaluateWithDeriv(vorder, vprime, global_vcoeff, global_vcoeffDeriv);
+	global_vorder = vorder;
+	global_vprime = vprime;
+    }
+
+    for (j = 0; j < k; j++) {
+	data=baseData+j;
+	retPoint[j] = retdu[j] = retdv[j] = 0.0;
+	for (row = 0; row < uorder; row++)  {
+	    /* 
+	    ** Minor optimization.
+	    ** The col == 0 part of the loop is extracted so we don't
+	    ** have to initialize p and pdv to 0.
+	    */
+	    p = global_vcoeff[0] * (*data);
+	    pdv = global_vcoeffDeriv[0] * (*data);
+	    data += k;
+	    for (col = 1; col < vorder; col++) {
+		/* Incrementally build up p, pdv value */
+		p += global_vcoeff[col] * (*data);
+		pdv += global_vcoeffDeriv[col] * (*data);
+		data += k;
+	    }
+	    /* Use p, pdv value to incrementally add up r, du, dv */
+	    retPoint[j] += global_ucoeff[row] * p;
+	    retdu[j] += global_ucoeffDeriv[row] * p;
+	    retdv[j] += global_ucoeff[row] * pdv;
+	}
+    }  
+}
+
+
+/*
+ *compute the Bezier polynomials C[n,j](v) for all j at v with 
+ *return values stored in coeff[], where 
+ *  C[n,j](v) = (n,j) * v^j * (1-v)^(n-j),
+ *  j=0,1,2,...,n.
+ *order : n+1
+ *vprime: v
+ *coeff : coeff[j]=C[n,j](v), this array store the returned values.
+ *The algorithm is a recursive scheme:
+ *   C[0,0]=1;
+ *   C[n,j](v) = (1-v)*C[n-1,j](v) + v*C[n-1,j-1](v), n>=1
+ *This code is copied from opengl/soft/so_eval.c:PreEvaluate
+ */
+void OpenGLSurfaceEvaluator::inPreEvaluate(int order, REAL vprime, REAL *coeff)
+{
+  int i, j;
+  REAL oldval, temp;
+  REAL oneMinusvprime;
+  
+  /*
+   * Minor optimization
+   * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to
+     * their i==1 loop values to avoid the initialization and the i==1 loop.
+     */
+  if (order == 1) {
+    coeff[0] = 1.0;
+    return;
+  }
+  
+  oneMinusvprime = 1-vprime;
+  coeff[0] = oneMinusvprime;
+  coeff[1] = vprime;
+  if (order == 2) return;
+  
+  for (i = 2; i < order; i++) {
+    oldval = coeff[0] * vprime;
+    coeff[0] = oneMinusvprime * coeff[0];
+    for (j = 1; j < i; j++) {
+      temp = oldval;
+      oldval = coeff[j] * vprime;
+	    coeff[j] = temp + oneMinusvprime * coeff[j];
+    }
+    coeff[j] = oldval;
+  }
+}
+
+/*
+ *compute the Bezier polynomials C[n,j](v) and derivatives for all j at v with 
+ *return values stored in coeff[] and coeffDeriv[].
+ *see the head of function inPreEvaluate for the definition of C[n,j](v)
+ *and how to compute the values. 
+ *The algorithm to compute the derivative is:
+ *   dC[0,0](v) = 0.
+ *   dC[n,j](v) = n*(dC[n-1,j-1](v) - dC[n-1,j](v)).
+ *
+ *This code is copied from opengl/soft/so_eval.c:PreEvaluateWidthDeriv
+ */
+void OpenGLSurfaceEvaluator::inPreEvaluateWithDeriv(int order, REAL vprime, 
+    REAL *coeff, REAL *coeffDeriv)
+{
+  int i, j;
+  REAL oldval, temp;
+  REAL oneMinusvprime;
+  
+  oneMinusvprime = 1-vprime;
+  /*
+   * Minor optimization
+   * Compute orders 1 and 2 outright, and set coeff[0], coeff[1] to 
+   * their i==1 loop values to avoid the initialization and the i==1 loop.
+   */
+  if (order == 1) {
+    coeff[0] = 1.0;
+    coeffDeriv[0] = 0.0;
+    return;
+  } else if (order == 2) {
+    coeffDeriv[0] = -1.0;
+    coeffDeriv[1] = 1.0;
+    coeff[0] = oneMinusvprime;
+    coeff[1] = vprime;
+    return;
+  }
+  coeff[0] = oneMinusvprime;
+  coeff[1] = vprime;
+  for (i = 2; i < order - 1; i++) {
+    oldval = coeff[0] * vprime;
+    coeff[0] = oneMinusvprime * coeff[0];
+    for (j = 1; j < i; j++) {
+      temp = oldval;
+      oldval = coeff[j] * vprime;
+      coeff[j] = temp + oneMinusvprime * coeff[j];
+    }
+    coeff[j] = oldval;
+  }
+  coeffDeriv[0] = -coeff[0];
+  /*
+   ** Minor optimization:
+   ** Would make this a "for (j=1; j<order-1; j++)" loop, but it is always
+   ** executed at least once, so this is more efficient.
+   */
+  j=1;
+  do {
+    coeffDeriv[j] = coeff[j-1] - coeff[j];
+    j++;
+  } while (j < order - 1);
+  coeffDeriv[j] = coeff[j-1];
+  
+  oldval = coeff[0] * vprime;
+  coeff[0] = oneMinusvprime * coeff[0];
+  for (j = 1; j < i; j++) {
+    temp = oldval;
+    oldval = coeff[j] * vprime;
+    coeff[j] = temp + oneMinusvprime * coeff[j];
+  }
+  coeff[j] = oldval;
+}
+
+void OpenGLSurfaceEvaluator::inEvalULine(int n_points, REAL v, REAL* u_vals, 
+	int stride, REAL ret_points[][3], REAL ret_normals[][3])
+{
+  int i,k;
+  REAL temp[4];
+inPreEvaluateBV_intfac(v);
+
+  for(i=0,k=0; i<n_points; i++, k += stride)
+    {
+      inDoEvalCoord2NOGE_BV(u_vals[k],v,temp, ret_normals[i]);
+
+      ret_points[i][0] = temp[0];
+      ret_points[i][1] = temp[1];
+      ret_points[i][2] = temp[2];
+
+    }
+
+}
+
+void OpenGLSurfaceEvaluator::inEvalVLine(int n_points, REAL u, REAL* v_vals, 
+	int stride, REAL ret_points[][3], REAL ret_normals[][3])
+{
+  int i,k;
+  REAL temp[4];
+inPreEvaluateBU_intfac(u);
+  for(i=0,k=0; i<n_points; i++, k += stride)
+    {
+      inDoEvalCoord2NOGE_BU(u, v_vals[k], temp, ret_normals[i]);
+      ret_points[i][0] = temp[0];
+      ret_points[i][1] = temp[1];
+      ret_points[i][2] = temp[2];
+    }
+}
+      
+
+/*triangulate a strip bounded by two lines which are parallel  to U-axis
+ *upperVerts: the verteces on the upper line
+ *lowerVertx: the verteces on the lower line
+ *n_upper >=1
+ *n_lower >=1
+ */
+void OpenGLSurfaceEvaluator::inEvalUStrip(int n_upper, REAL v_upper, REAL* upper_val, int n_lower, REAL v_lower, REAL* lower_val)
+{
+  int i,j,k,l;
+  REAL leftMostV[2];
+ typedef REAL REAL3[3];
+
+  REAL3* upperXYZ = (REAL3*) malloc(sizeof(REAL3)*n_upper);
+  assert(upperXYZ);
+  REAL3* upperNormal = (REAL3*) malloc(sizeof(REAL3) * n_upper);
+  assert(upperNormal);
+  REAL3* lowerXYZ = (REAL3*) malloc(sizeof(REAL3)*n_lower);
+  assert(lowerXYZ);
+  REAL3* lowerNormal = (REAL3*) malloc(sizeof(REAL3) * n_lower);
+  assert(lowerNormal);
+  
+  inEvalULine(n_upper, v_upper, upper_val,  1, upperXYZ, upperNormal);
+  inEvalULine(n_lower, v_lower, lower_val,  1, lowerXYZ, lowerNormal);
+
+
+
+  REAL* leftMostXYZ;
+  REAL* leftMostNormal;
+
+  /*
+   *the algorithm works by scanning from left to right.
+   *leftMostV: the left most of the remaining verteces (on both upper and lower).
+   *           it could an element of upperVerts or lowerVerts.
+   *i: upperVerts[i] is the first vertex to the right of leftMostV on upper line   *j: lowerVerts[j] is the first vertex to the right of leftMostV on lower line   */
+
+  /*initialize i,j,and leftMostV
+   */
+  if(upper_val[0] <= lower_val[0])
+    {
+      i=1;
+      j=0;
+
+      leftMostV[0] = upper_val[0];
+      leftMostV[1] = v_upper;
+      leftMostXYZ = upperXYZ[0];
+      leftMostNormal = upperNormal[0];
+    }
+  else
+    {
+      i=0;
+      j=1;
+
+      leftMostV[0] = lower_val[0];
+      leftMostV[1] = v_lower;
+
+      leftMostXYZ = lowerXYZ[0];
+      leftMostNormal = lowerNormal[0];
+    }
+  
+  /*the main loop.
+   *the invariance is that: 
+   *at the beginning of each loop, the meaning of i,j,and leftMostV are 
+   *maintained
+   */
+  while(1)
+    {
+      if(i >= n_upper) /*case1: no more in upper*/
+        {
+          if(j<n_lower-1) /*at least two vertices in lower*/
+            {
+              bgntfan();
+	      glNormal3fv(leftMostNormal);
+              glVertex3fv(leftMostXYZ);
+
+              while(j<n_lower){
+		glNormal3fv(lowerNormal[j]);
+		glVertex3fv(lowerXYZ[j]);
+		j++;
+
+              }
+              endtfan();
+            }
+          break; /*exit the main loop*/
+        }
+      else if(j>= n_lower) /*case2: no more in lower*/
+        {
+          if(i<n_upper-1) /*at least two vertices in upper*/
+            {
+              bgntfan();
+	      glNormal3fv(leftMostNormal);
+	      glVertex3fv(leftMostXYZ);
+	      
+              for(k=n_upper-1; k>=i; k--) /*reverse order for two-side lighting*/
+		{
+		  glNormal3fv(upperNormal[k]);
+		  glVertex3fv(upperXYZ[k]);
+		}
+
+              endtfan();
+            }
+          break; /*exit the main loop*/
+        }
+      else /* case3: neither is empty, plus the leftMostV, there is at least one triangle to output*/
+        {
+          if(upper_val[i] <= lower_val[j])
+            {
+	      bgntfan();
+
+	      glNormal3fv(lowerNormal[j]);
+	      glVertex3fv(lowerXYZ[j]);
+
+              /*find the last k>=i such that 
+               *upperverts[k][0] <= lowerverts[j][0]
+               */
+              k=i;
+
+              while(k<n_upper)
+                {
+                  if(upper_val[k] > lower_val[j])
+                    break;
+                  k++;
+
+                }
+              k--;
+
+
+              for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
+                {
+		  glNormal3fv(upperNormal[l]);
+		  glVertex3fv(upperXYZ[l]);
+
+                }
+	      glNormal3fv(leftMostNormal);
+	      glVertex3fv(leftMostXYZ);
+
+              endtfan();
+
+              /*update i and leftMostV for next loop
+               */
+              i = k+1;
+
+	      leftMostV[0] = upper_val[k];
+	      leftMostV[1] = v_upper;
+	      leftMostNormal = upperNormal[k];
+	      leftMostXYZ = upperXYZ[k];
+            }
+          else /*upperVerts[i][0] > lowerVerts[j][0]*/
+            {
+	      bgntfan();
+	      glNormal3fv(upperNormal[i]);
+	      glVertex3fv(upperXYZ[i]);
+	      
+              glNormal3fv(leftMostNormal);
+	      glVertex3fv(leftMostXYZ);
+	      
+
+              /*find the last k>=j such that
+               *lowerverts[k][0] < upperverts[i][0]
+               */
+              k=j;
+              while(k< n_lower)
+                {
+                  if(lower_val[k] >= upper_val[i])
+                    break;
+		  glNormal3fv(lowerNormal[k]);
+		  glVertex3fv(lowerXYZ[k]);
+
+                  k++;
+                }
+              endtfan();
+
+              /*update j and leftMostV for next loop
+               */
+              j=k;
+	      leftMostV[0] = lower_val[j-1];
+	      leftMostV[1] = v_lower;
+
+	      leftMostNormal = lowerNormal[j-1];
+	      leftMostXYZ = lowerXYZ[j-1];
+            }     
+        }
+    }
+  //clean up 
+  free(upperXYZ);
+  free(lowerXYZ);
+  free(upperNormal);
+  free(lowerNormal);
+}
+
+/*triangulate a strip bounded by two lines which are parallel  to V-axis
+ *leftVerts: the verteces on the left line
+ *rightVertx: the verteces on the right line
+ *n_left >=1
+ *n_right >=1
+ */
+void OpenGLSurfaceEvaluator::inEvalVStrip(int n_left, REAL u_left, REAL* left_val, int n_right, REAL u_right, REAL* right_val)
+{
+  int i,j,k,l;
+  REAL botMostV[2];
+  typedef REAL REAL3[3];
+
+  REAL3* leftXYZ = (REAL3*) malloc(sizeof(REAL3)*n_left);
+  assert(leftXYZ);
+  REAL3* leftNormal = (REAL3*) malloc(sizeof(REAL3) * n_left);
+  assert(leftNormal);
+  REAL3* rightXYZ = (REAL3*) malloc(sizeof(REAL3)*n_right);
+  assert(rightXYZ);
+  REAL3* rightNormal = (REAL3*) malloc(sizeof(REAL3) * n_right);
+  assert(rightNormal);
+  
+  inEvalVLine(n_left, u_left, left_val,  1, leftXYZ, leftNormal);
+  inEvalVLine(n_right, u_right, right_val,  1, rightXYZ, rightNormal);
+
+
+
+  REAL* botMostXYZ;
+  REAL* botMostNormal;
+
+  /*
+   *the algorithm works by scanning from bot to top.
+   *botMostV: the bot most of the remaining verteces (on both left and right).
+   *           it could an element of leftVerts or rightVerts.
+   *i: leftVerts[i] is the first vertex to the top of botMostV on left line   
+   *j: rightVerts[j] is the first vertex to the top of botMostV on rightline   */
+
+  /*initialize i,j,and botMostV
+   */
+  if(left_val[0] <= right_val[0])
+    {
+      i=1;
+      j=0;
+
+      botMostV[0] = u_left;
+      botMostV[1] = left_val[0];
+      botMostXYZ = leftXYZ[0];
+      botMostNormal = leftNormal[0];
+    }
+  else
+    {
+      i=0;
+      j=1;
+
+      botMostV[0] = u_right;
+      botMostV[1] = right_val[0];
+
+      botMostXYZ = rightXYZ[0];
+      botMostNormal = rightNormal[0];
+    }
+  
+  /*the main loop.
+   *the invariance is that: 
+   *at the beginning of each loop, the meaning of i,j,and botMostV are 
+   *maintained
+   */
+  while(1)
+    {
+      if(i >= n_left) /*case1: no more in left*/
+        {
+          if(j<n_right-1) /*at least two vertices in right*/
+            {
+              bgntfan();
+	      glNormal3fv(botMostNormal);
+              glVertex3fv(botMostXYZ);
+
+              while(j<n_right){
+		glNormal3fv(rightNormal[j]);
+		glVertex3fv(rightXYZ[j]);
+		j++;
+
+              }
+              endtfan();
+            }
+          break; /*exit the main loop*/
+        }
+      else if(j>= n_right) /*case2: no more in right*/
+        {
+          if(i<n_left-1) /*at least two vertices in left*/
+            {
+              bgntfan();
+	      glNormal3fv(botMostNormal);
+	      glVertex3fv(botMostXYZ);
+	      
+              for(k=n_left-1; k>=i; k--) /*reverse order for two-side lighting*/
+		{
+		  glNormal3fv(leftNormal[k]);
+		  glVertex3fv(leftXYZ[k]);
+		}
+
+              endtfan();
+            }
+          break; /*exit the main loop*/
+        }
+      else /* case3: neither is empty, plus the botMostV, there is at least one triangle to output*/
+        {
+          if(left_val[i] <= right_val[j])
+            {
+	      bgntfan();
+
+	      glNormal3fv(rightNormal[j]);
+	      glVertex3fv(rightXYZ[j]);
+
+              /*find the last k>=i such that 
+               *leftverts[k][0] <= rightverts[j][0]
+               */
+              k=i;
+
+              while(k<n_left)
+                {
+                  if(left_val[k] > right_val[j])
+                    break;
+                  k++;
+
+                }
+              k--;
+
+
+              for(l=k; l>=i; l--)/*the reverse is for two-side lighting*/
+                {
+		  glNormal3fv(leftNormal[l]);
+		  glVertex3fv(leftXYZ[l]);
+
+                }
+	      glNormal3fv(botMostNormal);
+	      glVertex3fv(botMostXYZ);
+
+              endtfan();
+
+              /*update i and botMostV for next loop
+               */
+              i = k+1;
+
+	      botMostV[0] = u_left;
+	      botMostV[1] = left_val[k];
+	      botMostNormal = leftNormal[k];
+	      botMostXYZ = leftXYZ[k];
+            }
+          else /*left_val[i] > right_val[j])*/
+            {
+	      bgntfan();
+	      glNormal3fv(leftNormal[i]);
+	      glVertex3fv(leftXYZ[i]);
+	      
+              glNormal3fv(botMostNormal);
+	      glVertex3fv(botMostXYZ);
+	      
+
+              /*find the last k>=j such that
+               *rightverts[k][0] < leftverts[i][0]
+               */
+              k=j;
+              while(k< n_right)
+                {
+                  if(right_val[k] >= left_val[i])
+                    break;
+		  glNormal3fv(rightNormal[k]);
+		  glVertex3fv(rightXYZ[k]);
+
+                  k++;
+                }
+              endtfan();
+
+              /*update j and botMostV for next loop
+               */
+              j=k;
+	      botMostV[0] = u_right;
+	      botMostV[1] = right_val[j-1];
+
+	      botMostNormal = rightNormal[j-1];
+	      botMostXYZ = rightXYZ[j-1];
+            }     
+        }
+    }
+  //clean up 
+  free(leftXYZ);
+  free(rightXYZ);
+  free(leftNormal);
+  free(rightNormal);
+}
+
+/*-----------------------begin evalMachine-------------------*/
+void OpenGLSurfaceEvaluator::inMap2fEM(int which, int k,
+	     REAL ulower,
+	     REAL uupper,
+	     int ustride,
+	     int uorder,
+	     REAL vlower,
+	     REAL vupper,
+	     int vstride,
+	     int vorder,
+	     REAL *ctlPoints)
+{
+  int i,j,x;
+  surfEvalMachine *temp_em;
+  switch(which){
+  case 0: //vertex
+    vertex_flag = 1;
+    temp_em = &em_vertex;
+    break;
+  case 1: //normal
+    normal_flag = 1;
+    temp_em = &em_normal;
+    break;
+  case 2: //color
+    color_flag = 1;
+    temp_em = &em_color;
+    break;
+  default:
+    texcoord_flag = 1;
+    temp_em = &em_texcoord;
+    break;
+  }
+
+  REAL *data = temp_em->ctlPoints;
+  
+  temp_em->uprime = -1;//initilized
+  temp_em->vprime = -1;
+
+  temp_em->k = k;
+  temp_em->u1 = ulower;
+  temp_em->u2 = uupper;
+  temp_em->ustride = ustride;
+  temp_em->uorder = uorder;
+  temp_em->v1 = vlower;
+  temp_em->v2 = vupper;
+  temp_em->vstride = vstride;
+  temp_em->vorder = vorder;
+
+  /*copy the contrl points from ctlPoints to global_ev_ctlPoints*/
+  for (i=0; i<uorder; i++) {
+    for (j=0; j<vorder; j++) {
+      for (x=0; x<k; x++) {
+	data[x] = ctlPoints[x];
+      }
+      ctlPoints += vstride;
+      data += k;
+    }
+    ctlPoints += ustride - vstride * vorder;
+  }
+}
+
+void OpenGLSurfaceEvaluator::inDoDomain2WithDerivsEM(surfEvalMachine *em, REAL u, REAL v, 
+				REAL *retPoint, REAL *retdu, REAL *retdv)
+{
+    int j, row, col;
+    REAL the_uprime;
+    REAL the_vprime;
+    REAL p;
+    REAL pdv;
+    REAL *data;
+
+    if((em->u2 == em->u1) || (em->v2 == em->v1))
+	return;
+    the_uprime = (u - em->u1) / (em->u2 - em->u1);
+    the_vprime = (v - em->v1) / (em->v2 - em->v1);
+    
+    /* Compute coefficients for values and derivs */
+
+    /* Use already cached values if possible */
+    if(em->uprime != the_uprime) {
+        inPreEvaluateWithDeriv(em->uorder, the_uprime, em->ucoeff, em->ucoeffDeriv);
+	em->uprime = the_uprime;
+    }
+    if (em->vprime != the_vprime) {
+	inPreEvaluateWithDeriv(em->vorder, the_vprime, em->vcoeff, em->vcoeffDeriv);
+	em->vprime = the_vprime;
+    }
+
+    for (j = 0; j < em->k; j++) {
+	data=em->ctlPoints+j;
+	retPoint[j] = retdu[j] = retdv[j] = 0.0;
+	for (row = 0; row < em->uorder; row++)  {
+	    /* 
+	    ** Minor optimization.
+	    ** The col == 0 part of the loop is extracted so we don't
+	    ** have to initialize p and pdv to 0.
+	    */
+	    p = em->vcoeff[0] * (*data);
+	    pdv = em->vcoeffDeriv[0] * (*data);
+	    data += em->k;
+	    for (col = 1; col < em->vorder; col++) {
+		/* Incrementally build up p, pdv value */
+		p += em->vcoeff[col] * (*data);
+		pdv += em->vcoeffDeriv[col] * (*data);
+		data += em->k;
+	    }
+	    /* Use p, pdv value to incrementally add up r, du, dv */
+	    retPoint[j] += em->ucoeff[row] * p;
+	    retdu[j] += em->ucoeffDeriv[row] * p;
+	    retdv[j] += em->ucoeff[row] * pdv;
+	}
+    }  
+}  
+
+void OpenGLSurfaceEvaluator::inDoDomain2EM(surfEvalMachine *em, REAL u, REAL v, 
+				REAL *retPoint)
+{
+    int j, row, col;
+    REAL the_uprime;
+    REAL the_vprime;
+    REAL p;
+    REAL *data;
+
+    if((em->u2 == em->u1) || (em->v2 == em->v1))
+	return;
+    the_uprime = (u - em->u1) / (em->u2 - em->u1);
+    the_vprime = (v - em->v1) / (em->v2 - em->v1);
+    
+    /* Compute coefficients for values and derivs */
+
+    /* Use already cached values if possible */
+    if(em->uprime != the_uprime) {
+        inPreEvaluate(em->uorder, the_uprime, em->ucoeff);
+	em->uprime = the_uprime;
+    }
+    if (em->vprime != the_vprime) {
+	inPreEvaluate(em->vorder, the_vprime, em->vcoeff);
+	em->vprime = the_vprime;
+    }
+
+    for (j = 0; j < em->k; j++) {
+	data=em->ctlPoints+j;
+	retPoint[j] = 0.0;
+	for (row = 0; row < em->uorder; row++)  {
+	    /* 
+	    ** Minor optimization.
+	    ** The col == 0 part of the loop is extracted so we don't
+	    ** have to initialize p and pdv to 0.
+	    */
+	    p = em->vcoeff[0] * (*data);
+	    data += em->k;
+	    for (col = 1; col < em->vorder; col++) {
+		/* Incrementally build up p, pdv value */
+		p += em->vcoeff[col] * (*data);
+		data += em->k;
+	    }
+	    /* Use p, pdv value to incrementally add up r, du, dv */
+	    retPoint[j] += em->ucoeff[row] * p;
+	}
+    }  
+}  
+
+
+void OpenGLSurfaceEvaluator::inDoEvalCoord2EM(REAL u, REAL v)
+{
+  REAL temp_vertex[5];
+  REAL temp_normal[3];
+  REAL temp_color[4];
+  REAL temp_texcoord[4];
+
+  if(texcoord_flag)
+    {
+      inDoDomain2EM(&em_texcoord, u,v, temp_texcoord);
+      texcoordCallBack(temp_texcoord, userData);
+    }
+  if(color_flag)
+    {
+      inDoDomain2EM(&em_color, u,v, temp_color);
+      colorCallBack(temp_color, userData);
+    }
+
+  if(normal_flag) //there is a normla map
+    {
+      inDoDomain2EM(&em_normal, u,v, temp_normal);
+      normalCallBack(temp_normal, userData);
+    
+      if(vertex_flag)
+	{
+	  inDoDomain2EM(&em_vertex, u,v,temp_vertex);
+	  if(em_vertex.k == 4)
+	    {
+	      temp_vertex[0] /= temp_vertex[3];
+	      temp_vertex[1] /= temp_vertex[3];
+	      temp_vertex[2] /= temp_vertex[3];	      
+	    }
+          temp_vertex[3]=u;
+          temp_vertex[4]=v;	  
+	  vertexCallBack(temp_vertex, userData);
+	}
+    }
+  else if(auto_normal_flag) //no normal map but there is a normal callbackfunctin
+    {
+      REAL du[4];
+      REAL dv[4];
+      
+      /*compute homegeneous point and partial derivatives*/
+      inDoDomain2WithDerivsEM(&em_vertex, u,v,temp_vertex,du,dv);
+
+      if(em_vertex.k ==4)
+	inComputeFirstPartials(temp_vertex, du, dv);
+
+#ifdef AVOID_ZERO_NORMAL
+      if(myabs(dv[0]) <= MYZERO && myabs(dv[1]) <= MYZERO && myabs(dv[2]) <= MYZERO)
+	{
+	  
+	  REAL tempdu[4];
+	  REAL tempdata[4];
+	  REAL u1 = em_vertex.u1;
+	  REAL u2 = em_vertex.u2;
+	  if(u-MYDELTA*(u2-u1) < u1)
+	    u = u+ MYDELTA*(u2-u1);
+	  else
+	    u = u-MYDELTA*(u2-u1);
+	  inDoDomain2WithDerivsEM(&em_vertex,u,v, tempdata, tempdu, dv);
+
+	  if(em_vertex.k ==4)
+	    inComputeFirstPartials(temp_vertex, du, dv);	  
+	}
+      else if(myabs(du[0]) <= MYZERO && myabs(du[1]) <= MYZERO && myabs(du[2]) <= MYZERO)
+	{
+	  REAL tempdv[4];
+	  REAL tempdata[4];
+	  REAL v1 = em_vertex.v1;
+	  REAL v2 = em_vertex.v2;
+	  if(v-MYDELTA*(v2-v1) < v1)
+	    v = v+ MYDELTA*(v2-v1);
+	  else
+	    v = v-MYDELTA*(v2-v1);
+	  inDoDomain2WithDerivsEM(&em_vertex,u,v, tempdata, du, tempdv);
+
+	  if(em_vertex.k ==4)
+	    inComputeFirstPartials(temp_vertex, du, dv);
+	}
+#endif
+
+      /*compute normal*/
+      switch(em_vertex.k){
+      case 3:
+
+	inComputeNormal2(du, dv, temp_normal);
+	break;
+      case 4:
+
+//	inComputeFirstPartials(temp_vertex, du, dv);
+	inComputeNormal2(du, dv, temp_normal);
+
+	/*transform the homegeneous coordinate of retPoint into inhomogenous one*/
+	temp_vertex[0] /= temp_vertex[3];
+	temp_vertex[1] /= temp_vertex[3];
+	temp_vertex[2] /= temp_vertex[3];
+	break;
+      }
+      normalCallBack(temp_normal, userData);
+      temp_vertex[3] = u;
+      temp_vertex[4] = v;
+      vertexCallBack(temp_vertex, userData);
+      
+    }/*end if auto_normal*/
+  else //no normal map, and no normal callback function
+    {
+      if(vertex_flag)
+	{
+	  inDoDomain2EM(&em_vertex, u,v,temp_vertex);
+	  if(em_vertex.k == 4)
+	    {
+	      temp_vertex[0] /= temp_vertex[3];
+	      temp_vertex[1] /= temp_vertex[3];
+	      temp_vertex[2] /= temp_vertex[3];	      
+	    }
+          temp_vertex[3] = u;
+          temp_vertex[4] = v;
+	  vertexCallBack(temp_vertex, userData);
+	}
+    }
+}
+
+
+void OpenGLSurfaceEvaluator::inBPMEvalEM(bezierPatchMesh* bpm)
+{
+  int i,j,k;
+  float u,v;
+
+  int ustride;
+  int vstride;
+
+#ifdef USE_LOD
+  if(bpm->bpatch != NULL)
+    {
+      bezierPatch* p=bpm->bpatch;
+      ustride = p->dimension * p->vorder;
+      vstride = p->dimension;
+
+      glMap2f( (p->dimension == 3)? GL_MAP2_VERTEX_3 : GL_MAP2_VERTEX_4,
+	      p->umin,
+	      p->umax,
+	      ustride,
+	      p->uorder,
+	      p->vmin,
+	      p->vmax,
+	      vstride,
+	      p->vorder,
+	      p->ctlpoints);
+
+
+/*
+    inMap2fEM(0, p->dimension,
+	  p->umin,
+	  p->umax,
+	  ustride,
+	  p->uorder,
+	  p->vmin,
+	  p->vmax,
+	  vstride,
+	  p->vorder,
+	  p->ctlpoints);
+*/
+    }
+#else
+
+  if(bpm->bpatch != NULL){
+    bezierPatch* p = bpm->bpatch;
+    ustride = p->dimension * p->vorder;
+    vstride = p->dimension;
+    inMap2fEM(0, p->dimension,
+	  p->umin,
+	  p->umax,
+	  ustride,
+	  p->uorder,
+	  p->vmin,
+	  p->vmax,
+	  vstride,
+	  p->vorder,
+	  p->ctlpoints);
+  }
+  if(bpm->bpatch_normal != NULL){
+    bezierPatch* p = bpm->bpatch_normal;
+    ustride = p->dimension * p->vorder;
+    vstride = p->dimension;
+    inMap2fEM(1, p->dimension,
+	  p->umin,
+	  p->umax,
+	  ustride,
+	  p->uorder,
+	  p->vmin,
+	  p->vmax,
+	  vstride,
+	  p->vorder,
+	  p->ctlpoints);
+  }
+  if(bpm->bpatch_color != NULL){
+    bezierPatch* p = bpm->bpatch_color;
+    ustride = p->dimension * p->vorder;
+    vstride = p->dimension;
+    inMap2fEM(2, p->dimension,
+	  p->umin,
+	  p->umax,
+	  ustride,
+	  p->uorder,
+	  p->vmin,
+	  p->vmax,
+	  vstride,
+	  p->vorder,
+	  p->ctlpoints);
+  }
+  if(bpm->bpatch_texcoord != NULL){
+    bezierPatch* p = bpm->bpatch_texcoord;
+    ustride = p->dimension * p->vorder;
+    vstride = p->dimension;
+    inMap2fEM(3, p->dimension,
+	  p->umin,
+	  p->umax,
+	  ustride,
+	  p->uorder,
+	  p->vmin,
+	  p->vmax,
+	  vstride,
+	  p->vorder,
+	  p->ctlpoints);
+  }
+#endif
+
+
+  k=0;
+  for(i=0; i<bpm->index_length_array; i++)
+    {
+#ifdef USE_LOD
+      if(bpm->type_array[i] == GL_POLYGON) //a mesh
+	{
+	  GLfloat *temp = bpm->UVarray+k;
+	  GLfloat u0 = temp[0];
+	  GLfloat v0 = temp[1];
+	  GLfloat u1 = temp[2];
+	  GLfloat v1 = temp[3];
+	  GLint nu = (GLint) ( temp[4]);
+	  GLint nv = (GLint) ( temp[5]);
+	  GLint umin = (GLint) ( temp[6]);
+	  GLint vmin = (GLint) ( temp[7]);
+	  GLint umax = (GLint) ( temp[8]);
+	  GLint vmax = (GLint) ( temp[9]);
+
+	  glMapGrid2f(LOD_eval_level*nu, u0, u1, LOD_eval_level*nv, v0, v1);
+	  glEvalMesh2(GL_FILL, LOD_eval_level*umin, LOD_eval_level*umax, LOD_eval_level*vmin, LOD_eval_level*vmax);
+	}
+      else
+	{
+	  LOD_eval(bpm->length_array[i], bpm->UVarray+k, bpm->type_array[i],
+		   0
+		   );
+	}
+	  k+= 2*bpm->length_array[i];       
+    
+#else //undef  USE_LOD
+
+#ifdef CRACK_TEST
+if(  bpm->bpatch->umin == 2 &&   bpm->bpatch->umax == 3
+  && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
+{
+REAL vertex[4];
+REAL normal[4];
+#ifdef DEBUG
+printf("***number ****1\n");
+#endif
+
+beginCallBack(GL_QUAD_STRIP, NULL);
+inDoEvalCoord2EM(3.0, 3.0);
+inDoEvalCoord2EM(2.0, 3.0);
+inDoEvalCoord2EM(3.0, 2.7);
+inDoEvalCoord2EM(2.0, 2.7);
+inDoEvalCoord2EM(3.0, 2.0);
+inDoEvalCoord2EM(2.0, 2.0);
+endCallBack(NULL);
+
+beginCallBack(GL_TRIANGLE_STRIP, NULL);
+inDoEvalCoord2EM(2.0, 3.0);
+inDoEvalCoord2EM(2.0, 2.0);
+inDoEvalCoord2EM(2.0, 2.7);
+endCallBack(NULL);
+
+}
+if(  bpm->bpatch->umin == 1 &&   bpm->bpatch->umax == 2
+  && bpm->bpatch->vmin ==2 &&    bpm->bpatch->vmax == 3)
+{
+#ifdef DEBUG
+printf("***number 3\n");
+#endif
+beginCallBack(GL_QUAD_STRIP, NULL);
+inDoEvalCoord2EM(2.0, 3.0);
+inDoEvalCoord2EM(1.0, 3.0);
+inDoEvalCoord2EM(2.0, 2.3);
+inDoEvalCoord2EM(1.0, 2.3);
+inDoEvalCoord2EM(2.0, 2.0);
+inDoEvalCoord2EM(1.0, 2.0);
+endCallBack(NULL);
+
+beginCallBack(GL_TRIANGLE_STRIP, NULL);
+inDoEvalCoord2EM(2.0, 2.3);
+inDoEvalCoord2EM(2.0, 2.0);
+inDoEvalCoord2EM(2.0, 3.0);
+endCallBack(NULL);
+
+}
+return;
+#endif //CRACK_TEST
+
+      beginCallBack(bpm->type_array[i], userData);
+
+      for(j=0; j<bpm->length_array[i]; j++)
+	{
+	  u = bpm->UVarray[k];
+	  v = bpm->UVarray[k+1];
+#ifdef USE_LOD
+          LOD_EVAL_COORD(u,v);
+//	  glEvalCoord2f(u,v);
+#else
+
+#ifdef  GENERIC_TEST
+          float temp_normal[3];
+          float temp_vertex[3];
+          if(temp_signal == 0)
+	    {
+	      gTessVertexSphere(u,v, temp_normal, temp_vertex);
+//printf("normal=(%f,%f,%f)\n", temp_normal[0], temp_normal[1], temp_normal[2])//printf("veretx=(%f,%f,%f)\n", temp_vertex[0], temp_vertex[1], temp_vertex[2]);
+              normalCallBack(temp_normal, userData);
+	      vertexCallBack(temp_vertex, userData);
+	    }
+          else if(temp_signal == 1)
+	    {
+	      gTessVertexCyl(u,v, temp_normal, temp_vertex);
+//printf("normal=(%f,%f,%f)\n", temp_normal[0], temp_normal[1], temp_normal[2])//printf("veretx=(%f,%f,%f)\n", temp_vertex[0], temp_vertex[1], temp_vertex[2]);
+              normalCallBack(temp_normal, userData);
+	      vertexCallBack(temp_vertex, userData);
+	    }
+	  else
+#endif //GENERIC_TEST
+
+	    inDoEvalCoord2EM(u,v);
+     
+#endif //USE_LOD
+
+	  k += 2;
+	}
+      endCallBack(userData);
+
+#endif //USE_LOD
+    }
+}
+
+void OpenGLSurfaceEvaluator::inBPMListEvalEM(bezierPatchMesh* list)
+{
+  bezierPatchMesh* temp;
+  for(temp = list; temp != NULL; temp = temp->next)
+    {
+      inBPMEvalEM(temp);
+    }
+}
+
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/mystdio.h b/mesalib/src/glu/sgi/libnurbs/interface/mystdio.h
new file mode 100644
index 000000000..5ab49effd
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/mystdio.h
@@ -0,0 +1,60 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * mystdio.h
+ *
+ */
+
+#ifndef __glumystdio_h_
+#define __glumystdio_h_
+
+#ifdef STANDALONE
+inline void _glu_dprintf( const char *, ... ) { }
+#endif
+
+#ifdef LIBRARYBUILD
+#ifndef NDEBUG
+#include <stdio.h>
+#define _glu_dprintf printf
+#else
+inline void _glu_dprintf( const char *, ... ) { }
+#endif
+#endif
+
+#ifdef GLBUILD
+inline void _glu_dprintf( const char *, ... ) { }
+#endif
+
+#ifndef NULL
+#define NULL		0
+#endif
+
+#endif /* __glumystdio_h_ */
diff --git a/mesalib/src/glu/sgi/libnurbs/interface/mystdlib.h b/mesalib/src/glu/sgi/libnurbs/interface/mystdlib.h
new file mode 100644
index 000000000..ab7a3b273
--- /dev/null
+++ b/mesalib/src/glu/sgi/libnurbs/interface/mystdlib.h
@@ -0,0 +1,57 @@
+/*
+ * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008)
+ * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved.
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice including the dates of first publication and
+ * either this permission notice or a reference to
+ * http://oss.sgi.com/projects/FreeB/
+ * shall be included in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+ * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
+ * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ *
+ * Except as contained in this notice, the name of Silicon Graphics, Inc.
+ * shall not be used in advertising or otherwise to promote the sale, use or
+ * other dealings in this Software without prior written authorization from
+ * Silicon Graphics, Inc.
+ */
+
+/*
+ * mystdlib.h
+ *
+ */
+
+#ifndef __glumystdlib_h_
+#define __glumystdlib_h_
+
+#ifdef STANDALONE
+typedef unsigned int size_t;
+extern "C" void 	abort( void );
+extern "C" void *	malloc( size_t );
+extern "C" void 	free( void * );
+#endif
+
+#ifdef LIBRARYBUILD
+#include <stdlib.h>
+#endif
+
+#ifdef GLBUILD
+typedef unsigned int size_t;
+extern "C" void 	abort( void );
+extern "C" void *	malloc( size_t );
+extern "C" void 	free( void * );
+#endif
+
+#endif /* __glumystdlib_h_ */