OpenGL Optimizer allows you to interact with an abstract object (a representation or rep) and treat rendering as a separate operation. A simple example of a rep is a sphere, defined by a radius and a center. After defining a sphere, you can implement how it is rendered in several ways: by tessellating, by a sphere-specific draw routine, or conceivably by hardware. Member functions of geometric-primitive classes allow you to implement the most appropriate way of rendering. The fundamental rendering step of tessellating a representation is discussed in Chapter 11, “Rendering Higher-Order Primitives: Tessellators.”

NURBS curves and surfaces are included in the set of OpenGL Optimizer reps. OpenGL also has these, but OpenGL Optimizer NURBS have two advantages:

See Chapter 10, “Creating and Maintaining Surface Topology,” and “opTessNurbSurfaceAction”.

OpenGL Optimizer uses the value for π held in the variable M_PI, declared in csBasic.h: 3.14159265358979323846.

The classes opVec2, opVec3, and opVec4 define two-, three-, and four-dimensional vectors, and include common operations of vector algebra such as addition, scalar multiplication, cross products, and so on. See the header file opVec.h for a list of operations defined for each vector.

The important distinction between these vector classes and csVec of Cosmo3D is that OpenGL Optimizer vectors are declared opReal and so can be single or double precision, depending on the version of the OpenGL Optimizer library you use.

The opScalar class is the base class for defining scalar functions; it allows you to conveniently read and write functions. The class provides a virtual evaluation method.

class opScalar : public csObject
{ public:
// Creating and destroying
opScalar();
virtual ~opScalar();
virtual opReal eval(opReal u); 
virtual csObject* clone(csNode::CloneEnum);
}; 

class opCompositeScalar : public opScalar
{ public:
// Creating and destroying
opCompositeScalar( );
opCompositeScalar(opScalar *outFun, opScalar *inFun);
virtual ~opCompositeScalar();

// Accessor functions
opScalar *getOutF() 
opScalar *getInF() 
void      setOutF(opScalar *outF);
void      setInF (opScalar *inF);

opReal eval(opReal t);

//Copy
virtual csObject clone(csNode::CloneEnum what);
};

class opPolyScalar : public opScalar
{
public:
    // Creating and destroying
opPolyScalar( void );
opPolyScalar( int degree, opReal* coef); 
virtual ~opPolyScalar();
    // Accessor functions
void set( int degree, opReal* coef);
int getDegree() 
opReal getCoef( int i) 
    // Evaluators
opReal eval(opReal u);
    //Copy
virtual csObject* clone ( csNode::CloneEnum what);
}:

Each geometric primitive is defined with respect to its own coordinate system. The elementary definition of an object gives a particular orientation and location with respect to the origin. This reference frame can, in turn, be manipulated by a csTransform to place it in a scene or manipulate it (see Chapter 7, “Interactive Highlighting and Manipulating”).

The base class for higher-order primitives has methods that allow you to locate and orient a primitive with respect to its own reference frame. These methods make insertion of csTransform nodes whenever you want to define the location or orientation of an object or to change the shape of an object unnecessary.

The location is defined by an opVec2 or opVec3, and the orientation is controlled by a 3 x 3 matrix, held in the class opFrame. If the matrix is not a rotation matrix, you can change the shape of an object, for example, you can distort a sphere into an ellipsoid.

class opFrame
{
public:
opReal m[3][3];
bool   identity;

void setIdentity()
opFrame();
};

The class has the following main methods:

class opRep : public csShape
{
public:
opRep( );
virtual ~opRep( );

// Accessor functions
void setOrigin( const opVec3&  org );
void setOrient( const opFrame& mat );

opVec3  getOrigin();
opFrame getOrient();

// Utility methods
virtual int getMemSize();

public: 
// Comso3d virtual functions
virtual void draw(csDrawAction* action);
virtual void isect(csIsectAction* ia);
};

opRep's subclasses typically include evaluator methods to determine coordinates of points, tangents, and normals. If you do not want the values corresponding to the default position, do not call these methods before you use setOrient() and setOrigin() to locate an opRep. Thus, for example, when defining points on a circle, first set the center and the radius, then call setOrient() to set the orientation, and then evaluate points.

Planar curves consist of sets of points, described by two-dimensional vectors, opVec2s. They are parameterized by the opReal variable t; as t varies, a point “moves” along the curve. t can be thought of as the amount of time that has passed as a point moves along the curve. Or, t can measure the distance traveled.

More precisely, each component of a point on the curve is a function of t, which lies in the parameter interval (t₀, t₁) on the real line. Points on the curve are described by a pair of functions of t: (x(t), y(t)).

Figure 9-2. Parametric Curve: Parameter Interval (0,1).

Classes derived from opCurve2d inherit a set of evaluator functions which, for a given value of t, evaluate a point on the curve, the tangent and normal vectors at the point, and the curvature. Naturally, the base-class evaluator that locates points on the curve is a pure virtual function.

To evaluate tangent and normal vectors at a point, opCurve2d provides virtual functions that, by default, use finite-central-difference calculations. To compute the tangent to the curve at p[t], a point on the curve, the tangent evaluator function takes the vector connecting two “nearby” points on the curve, p[t+Dt] - p[t-Dt] where Δτ is “small,” and divides by 2Δτ. Similarly, a finite-central-difference calculation of the normal vector uses the difference between two nearby tangent vectors: n[t] = (t[t+Dt] -t[t-Dt])/2Δτ.

class opCurve2d : public opRep
{
public:
// Creating and destroying
opCurve2d( );
opCurve2d( opReal beginT, opReal endT  );

virtual ~opCurve2d();

// Accessor functions
void setBeginT( opReal beginT );
void setEndT(   opReal endT );

opReal getBeginT();
opReal getEndT();

opVec2 getBeginPt();
opVec2 getEndPt();

opVec2 getBeginTan();
opVec2 getEndTan();

void   setClosed( opLoop loopVal );
opLoop getClosed();

void  setClosedTol( opReal tol );
opReal getClosedTol();

// Evaluators
virtual void evalPt(    opReal t, opVec2 &pnt ) = 0;
virtual void evalTan(   opReal t, opVec2 &tan );
virtual void evalNorm(  opReal t, opVec2 &norm );
virtual void evalCurv(  opReal t, opReal &curv );
virtual void eval(      opReal t, opVec2 &pnt, 
                                  opVec2 &tan, 
                                  opReal &curv, 
                                  opVec2 &norm );
};

Parametric lines in the plane are defined by beginning and ending points. The parameterization is such that as t varies from t1 to t2, a point on the line “moves,” at a uniform rate, from the beginning to the ending point.

Figure 9-3. Line in the Plane Parameterization

class opLine2d : public opCurve2d
{
public:
// Creating and destroying
opLine2d();
opLine2d( opReal x1, opReal y1, opReal t1,
          opReal x2, opReal y2, opReal t2 );
virtual ~opLine2d();

// Accessor functions
void setPoint1( opReal x1, opReal y1, opReal t1 );
void setPoint2( opReal x2, opReal y2, opReal t2 );
void getPoint1( opReal *x1, opReal *y1, opReal *t1 );
void getPoint2( opReal *x2, opReal *y2, opReal *t2 );
// Evaluators
void evalPt(   opReal t, opVec2 &pnt );
//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

Use the class opCircle2d to define a parametric circle in the plane. The parameterization is such that t is the angular displacement, in radians, in a counterclockwise direction from the x axis. Figure 9-4 illustrates the parameterization of the circle.

Figure 9-4. Circle in the Plane Parameterization

Class Declaration for opCircle2d

The class has the following main methods:

class opCircle2d : public opCurve2d
{
public:
    // Creating amd destroying
opCircle2d();
opCircle2d( opReal rad, opVec2 *org );
virtual ~opCircle2d();
    // Accessor functions
void   setRadius( opReal rad ) ;
opReal getRadius()             ;
    // Evaluator
void evalPt(   opReal t, opVec2 &pnt );
void evalTan(  opReal t, opVec2 &tan );
void evalCurv( opReal t, opReal &curv );
void evalNorm( opReal t, opVec2 &norm );
void eval(     opReal t, 
               opVec2 &pnt, 
               opVec2 &tan, 
               opReal &curv,
               opVec2& norm );
    //Copy
virtual csNode* clone(csNode::CloneEnum what)
};

The class opSuperQuadCurve2d provides methods to define a generalization of a circle that, when used for constructing a swept surface, is convenient for generating rounded, nearly square surfaces, or surfaces with sharp cusps (see “Swept Surfaces”). Two examples of superquadrics appear in repTest.

The position along the curve is specified by an angle from the x axis, in the same as for an opCircle2d. The shape of the curve is controlled by a second parameter.

A superquadric is the set of points (x,y) given by the following equation that clearly expresses the relationship to the equation of a circle:

The above equation can be written in a parametric form:

The family of curves generated by these equations as the quantity α varies can be described as follows (see Figure 9-5).

Four points are always on the curve for any value of a: (\ξβ1 ρ, 0) and (0, \ξβ1 ρ).

class opSuperQuadCurve2d : public opCurve2d
{
public:
// Creating and destroying
opSuperQuadCurve2d();
opSuperQuadCurve2d( opReal radius, 
                    opVec2 *origin,
                    opReal exponent );
virtual ~opSuperQuadCurve2d();

// Accessor functions
void   setRadius( opReal _radius );
opReal getRadius();

void   setExponent( opReal _exponent );
opReal getExponent();

// Evaluator
void evalPt(   opReal t, opVec2 &pnt );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

A spline is a mathematical technique for generating a single geometric object from pieces. An advantage of breaking a curve into pieces is greater flexibility when you have many points controlling the shape: changes to one piece of the curve do not have significant effects on remote pieces. To define a spline curve for a range of values for the parameter t, say from 0 to 3, you “tie” together pieces of curves defined over intervals of values for t. For example, you might assign curve pieces to the three intervals 0 to 1, 1 to 2, and 2 to 3. The four points in the set of parameters, 0, 1, 2, and 3, define the endpoints of the intervals and are called knots.

A Hermite-spline curve is a curve whose segments are cubic polynomials of the parameter t, where the coefficients of the polynomials are determined by the position and tangent to the curve at each knot point. Thus the curve passes through each of a set of specified points with a specified tangent vector. The set of knot points must be increasing values of the parameter t.

Figure 9-6. Hermite Spline Curve Parameterization

class opHsplineCurve2d : public opCurve2d
{
public:
// Creating and destroying
opHsplineCurve2d( opReal tBegin = 0.0, opReal tEnd = 1.0 );
virtual ~opHsplineCurve2d( );

// Accessor functions  
void setPoint(   int i, opVec2 &p );
void setTangent( int i, opVec2 &tng );
void setKnot(    int i, opReal t );

int getKnotCount();
opVec2* getPoint( int i );
opVec2* getTangent( int i );
opReal  getKnot( int i );

virtual int getMemSize( );

// Evaluator
virtual void evalPt( opReal t, opVec2 &pnt );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

The acronym NURBS stands for “nonuniform rational B-splines.” NURBS define a set of curves and surfaces that generalizes Bezier curves. Both NURBS curves and Bezier curves are “smooth” curves that are well suited for CAD design work. They are essentially determined by a set of points that controls the shape of the curves, although the points do not lie on the curves.

Because NURBS properties are not widely known, a discussion of their features precedes details of how to create instances of them. The discussion is necessarily brief and is intended to provide the minimum amount of information needed to start using OpenGL Optimizer NURBS classes.

This general discussion of NURBS is presented in the following sections:

For more information, consult the following sources, which are listed in “Recommended Background Reading”:

Equation Used to Calculate a NURBS Curve

The equation that defines the NURBS curve is

Alternative Equation for a NURBS Curve

If you have a surface developed from the alternative expression for a NURBS surface:

you must change the coordinates of the control points to get the same surface from OpenGL Optimizer; you convert the coordinates of the control points from (x,y,w) to (wx,wy,w).

The class opNurbCurve2d defines a nonuniform rational B-spline curve in the plane, the simplest NURBS object provided by OpenGL Optimizer.

class opNurbCurve2d : public opCurve2d
{
public:
// Creating and destroying
opNurbCurve2d( opReal tBegin = 0.0, opReal tEnd = 1.0 );
virtual ~opNurbCurve2d( );

// Accessor functions
void setControlHull(     int i, opVec2 &p );
void setControlHull(     int i, opVec3 &p );
void setWeight(          int i, opReal w  );
void setKnot(            int i, opReal t  );
void setControlHullSize( int s            );

opVec2* getControlHull( int i );
opReal  getWeight( int i );
int     getControlHullSize( );
int     getKnotCount( );
opReal  getKnot( int i );
int     getOrder( );

void removeControlHullPnt(int i);
void removeKnot(int i);

// Evaluator
virtual void evalPt( opReal t, opVec2 &pnt );

// Memory foot print
virtual int getMemSize( );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

A piecewise polynomial curve consists of an array of polynomial curves. Each polynomial curve is a polynomial mapping from t to UV plane, where the domain is a subinterval of [0,1]. The polynomial coefficients are set by setControlHull().

Notice that an opPieceWisePolyCurve2d is a subclass of opCurve2d. The domain of a opPieceWisePolyCurve2d is defined to be [0, n] where n is the number of pieces.

If reverse is 0, then for any given t in [0, n], its corresponding uv is evaluated in the following way: The index of the piece that corresponds to t is floor(t), and the polynomial of that piece is evaluated at w1 + (t-floor(t)) * (w2-w1) to get the (u,v), where [w1, w2] is the domain interval (set by setLimitParas()) of this piece.

If reverse is 1, then for any given t in [0,n], we first transform t into n-t, then perform the normal evaluation (at n-t) as described in the preceding paragraph.

class opPieceWisePolyCurve2d : public opCurve3d
{
public:
// Creating and destroying
opPieceWisePolyCurve2d ( );
~opPieceWisePolyCurve2d ( );

//Accessor functions
void setControlHull ( int piece, int i, const opVec2& p);
opVec2& getControlHull ( int piece, int i);
void setLimitParas ( int piece, opReal w1, opReal w2);
void setReverse ( int _reverse);
opVec2& getLimitParas ( int piece);
int getReverse ( );
int getPatchCount ( );
int getOrder ( int piece);
virtual void evalPt ( opReal t, opVec2& pnt);
virtual void evalBreakPoints ( opParaSurface* sur);

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

The class opDisCurve2d is the base class for making a discrete curve from line segments connecting a sequence of points in the plane. Because opDisCurve2d is not derived from opCurve2d, it does not inherit that class's finite difference functions for calculating derivatives, therefore, opDisCurve2d includes member functions that calculate arc length, tangents, principal normals, and curvatures using finite central differences. Figure 9-7 illustrates the definition of the curve by a set of points.

Figure 9-7. Discrete Curve Definition

class opDisCurve2d : public opRep
{
public:
// Creating and destroying
opDisCurve2d( void );
opDisCurve2d( int nPoints, opReal *points );

virtual ~opDisCurve2d( void );

// Accessor functions
void set (int nPoints, opReal* points);
opVec2 getBeginPt();
opVec2 getEndPt();

opLoop getClosed();
void   setClosed( opLoop c );

void setPoint( int i, const opVec2& pnt ); 
opVec2 getPoint( int i);
int getPointCount(); 
opVec2 getTangent(int i); 
opVec2 getNormal(int i); 
opReal getCurvature(int i);

// Evaluators
void computeTangents( );
void computeNormals( );
void computeCurvatures( );
void computeDerivatives( );

// Copy
virtual csNode* clone(csNode::CloneEnum::what);
};

The base class for lines in space, opLine3d, is essentially the same as opLine2d, discussed in “Lines in the Plane”. The main differences are due to the need to manage three-dimensional vectors. Thus all vector variables are opVec3 and the constructor takes six variables to define the endpoints of the line.

The default orientation of the curve is identical to that for the planar curve opLine2d; you can translate and rotate the line in three-dimensional space with the methods setOrigin() and setOrient() inherited from opRep.

The class opCircle3d defines a parametric circle with an arbitrary location and orientation in space. The parameterization of the circle, before you change its location or orientation, is such that t is the angular displacement, in radians, in a counterclockwise direction from the x axis.

The class declaration for opCircle3d is identical to that for opCircle2d, discussed in “Circles in the Plane”, except for the changes from opVec2 to opVec3. The member functions perform the same operations. For more information, see the discussion in the section “Circles in the Plane”.

If the matrix you use to orient an opCircle3d does not correspond to a rotation about an axis—that is, the matrix is not orthonormal— you not only change the tilt of the plane in which the circle lies but also change the radius, and may distort the circle into an ellipse. For a discussion of useful matrices, see the book by J. D. Foley, et al., in “Recommended Background Reading”.

The class opSuperQuadCurve3d provides methods to define a superquadric in space (see “Superquadric Curves: opSuperQuadCurve2d”). The class declaration is identical to that for opSuperQuad2d except that position on the curve is defined by an opVec3.

The default orientation of the curve is identical to that for the planar curve opSuperQuad2d; you can translate and rotate the curve in three-dimensional space with the methods setOrigin() and setOrient() inherited from opRep.

The class opHsplineCurve3d provides methods to define a Hermite spline curve in space. The definition of the curve is the same as that for a Hermite spline curve in the plane, discussed in “Hermite-Spline Curves in the Plane”. The class declaration is the same as that for opHsplineCurve2d, but the position and tangent vectors are opVec3s.

The basic properties of NURBS are discussed in the section “NURBS Overview”. In an effort to keep things as simple as possible, the discussion in that section has a bias toward curves in the plane. But the principles and control parameters are, with one difference, the same for NURBS curves in space.

The difference is that control points for NURBS curves in space can be anywhere in space instead of being restricted to a plane. The section “Examples of NURBS Curves” in Chapter 8 of The Inventor Mentor presents illustrations of NURBS curves in space, along with their control parameters.

The class opNurbCurve3d is the base class for NURBS curves in space. Its class declaration is practically identical to that for opNurbCurve2d but all occurrences of opVec2 are changed to opVec3. In addition, the vector argument of setControlHull() can be an opVec3, if you just want to specify control point locations, or an opVec4, if you want to append weighting information as a fourth component. See the discussion in the section “NURBS Curves in the Plane”.

A planar curve in the u-v plane describes a curve on the surface, given a parameterized surface (see the section “Parametric Surfaces”). Each point on the curve in the parameter plane is “lifted up” to the surface. Such curves are known as composite curves because they are described mathematically as the composition of the function describing the curve and the function describing the surface. The edge of a surface defined by a trim curve is a composite curve.

opCompositeCurve3d is the base class for composite curves. This class is useful for defining trim curves and surface silhouettes in the parametric surface's coordinate system.

class opCompositeCurve3d : public opCurve3d
{
public:
// Creating and destroying
opCompositeCurve3d( );
opCompositeCurve3d( opParaSurface *sur, opCurve2d *cur );
~opCompositeCurve3d( );

// Accessor functions
void set( opParaSurface *sur, opCurve2d *cur );
opParaSurface* getParaSurface() {return s;}
opCurve2d* getCurve2d() {return c;}

// Evaluator
void evalPt( opReal u, opVec3 &pnt );
};

The class opDisCurve3d is the base class for making a discrete curve of line segments connecting a sequence of points in space. The class declaration for opDisCurve3d is identical to that for opDisCurve2d, discussed in “Discrete Curves in the Plane”, but opVec2 changes to opVec3. The member functions perform the same operations.

One application of an opDisCurve3d and opHsplineCurve3d is to use them to interactively specify routing for tubing. These are the operations to perform:

To locate a point on a parametric surface, you need two parameters, referred to as u and v in OpenGL Optimizer. The set of u and v values that describe the surface are known as the parameter space, or coordinate system, of the surface (see Figure 9-8).

More precisely, the coordinates of the points in space that define a parametric surface are described by a set of three functions of two parameters: (x(u,v), y(u,v), z(u,v)).

Well-known examples of a parametric surface are a sphere or a globe. On a globe you can locate points with two parameters: latitude and longitude. The rectangular grid of latitude and longitudes is the coordinate system that describes points on the globe.

Figure 9-8. Parametric Surface: Unit-Square Coordinate System

To define the extent of a parametric curve, pick an interval. For accurate trimming of a parametric surface, you need more complex tools. You are likely to need:

OpenGL Optimizer keeps the trim loop side on the left as you look down on the u-v plane while a point moves along the curve in the direction of increasing t; you can hold on to the surface with your left hand as you go along the trim loop. Thus a clockwise loop removes a hole; a counterclockwise loop keeps the enclosed region and eliminates everything outside. Do not create a trim loop that crosses itself like a figure eight.

OpenGL Optimizer allows you to maintain curves to define the edges of a surface. These curves are opCurve2d objects defined in the u-v plane that are “lifted” to the surface by the parameterization. The main use of these curves is to eliminate a portion of the surface on one side of the curve. The name of a curve in the coordinate system that is used to define (possibly a piece of) such a surface edge is a trim curve. One or more joined trim curves form a sequence called a trim loop. To be of use, trim curves should form a closed loop or reach the edges of the coordinate system for the surface. Figure 9-9 illustrates trim loops and their effect on a surface.

Figure 9-9. Trim Loops and Trimmed Surface: Both Trim Loops Made of Four Trim Curves

An opEdge defines a trim curve in u, v space. opEdge holds information about a surface's adjacency. Each opEdge identifies an opBoundary, which the class opTopo uses to keep track of surface connectivity, and continuous and discrete versions of the trim curve associated with the boundary. The members of an opEdge are set by the toplogy building tools; the methods of opEdge access the members. Topology building and the classes opTopo and opBoundary are discussed further in Chapter 10, “Creating and Maintaining Surface Topology.”

The information held in opEdge allows tessellators to determine whether a set of vertices has already been developed for points shared with other surfaces. You can also find other surfaces that have the same edge or trim-curve endpoint as that defined by a given trim curve.

The set*() methods are mainly used when reading surface data from a file and creating OpenGL Optimizer data structures.

class opEdge
{
public:
// Creating and destroying
opEdge();
~opEdge();

opCurve2d *getContCurve();
void setContCurve(opCurve2d *c);

opDisCurve2d *getDisCurve();
void setDisCurve( opDisCurve2d *d);

int getBoundary();

void setBoundaryDir( int dir );
int getBoundaryDir();

opEdge* clone( csNode::CloneEnum what);
};

opParaSurface is the base class for parametric surfaces in OpenGL Optimizer. As for the base classes opCurve2d and opCurve3d, opParaSurface includes a pure virtual function to evaluate points on the surface and default evaluator functions that calculate derivatives using finite central differences. The surface normal at a point is the cross product of the partial derivatives.

For parametric curves whose extent is defined by the interval of values for t, the extent of an opParaSurface is, initially, defined by all the points in its parameter space.

class opParaSurface : public opRep
{
public:
// Creating and destroying
opParaSurface();
opParaSurface( opReal _beginU = 0, opReal _endU = 1, 
               opReal _beginV = 0, opReal _endV = 1,
               int    _topoId = 0, int    _solid_id = -1 );

~opParaSurface();

// Accessor functions
void setBeginU( opReal u );
void setEndU(   opReal u );
void setBeginV( opReal v );
void setEndV(   opReal v );
void setSolidId( int solidId);
void SetTopoId( int topoId);
void setSurfaceId (int surfaceId);
opReal getBeginU()
opReal getEndU() 
opReal getBeginV()
opReal getEndV()

int     getTrimLoopCount();
opLoop  getTrimLoopClosed( int loopNum );
int     getTrimCurveCount( int loopNum );
opEdge* getTrimCurve( int loopNum, int curveNum );

int     getTopoId();
int     getSolidId();
int     getSurfaceId();
  
void setHandednessHint( bool _clockWise ) 
bool getHandednessHint()
  
void insertTrimCurve( int loopNum, opCurve2d *c, opDisCurve2d *d );

// Explicit add a trim curve to a trim loop
void addTrimCurve(int loopNum, opCurve2d *c, opDisCurve2d *d );
void setTrimLoopClosed(  int loopNum, opLoop closed );

// Surface evaluators
virtual void evalPt(   opReal u, opReal v, opVec3 &pnt ) = 0;
virtual void evalDu(   opReal u, opReal v, opVec3 &Du );
virtual void evalDv(   opReal u, opReal v, opVec3 &Dv );
virtual void evalDuu(  opReal u, opReal v, opVec3 &Duu );
virtual void evalDvv(  opReal u, opReal v, opVec3 &Dvv );
virtual void evalDuv(  opReal u, opReal v, opVec3 &Duv );
virtual void evalNorm( opReal u, opReal v, opVec3 &norm );

// Directional derivative evaluators
virtual void evalD( opReal u, opReal v, opReal theta, opVec3 &D );
virtual void evalDD( opReal u, opReal v, opReal theta, opVec3 &DD );

virtual void eval( opReal u, opReal v,
opVec3 &p,            // The point
opVec3 &Du,           // The derivative in the u direction
opVec3 &Dv,           // The derivative in the v direction
opVec3 &Duu,          // The 2nd derivative in the u direction
opVec3 &Dvv,          // The 2nd derivative in the v direction
opVec3 &Duv,          // The cross derivative
opReal &s,            // Texture coordinates
opReal &t  );

void clearTessallation();
};

The simplest parametric surface is a plane. The class opPlane defines a plane by two parameter intervals and three points that define the two coordinate directions. Figure 9-10 illustrates the parameterization of an opPlane.

Figure 9-10. Plane Parameterization

class opPlane : public opParaSurface
{
public:
// Creating and destroying
opPlane();
opPlane( opReal x1, opReal y1, opReal z1, opReal u1, opReal v1,
         opReal x2, opReal y2, opReal z2, opReal u2,
         opReal x3, opReal y3, opReal z3, opReal v3 );
virtual ~opPlane();

// Accessor functions
void setPoint1( opReal x1, opReal y1, opReal z1, opReal u1, opReal v1);
void setPoint2( opReal x2, opReal y2, opReal z2, opReal u2 );
void setPoint3( opReal x3, opReal y3, opReal z3, opReal v3 );
 
void getPoint1( opReal *x1, opReal *y1, opReal *z1, 
                opReal *u1, opReal *v1 );
void getPoint2( opReal *x2, opReal *y2, opReal *z2, opReal *u2 );
void getPoint3( opReal *x3, opReal *y3, opReal *z3, opReal *v3 );

// Evaluators
void evalPt( opReal u, opReal v, opVec3 &pnt );
void evalDu( opReal u, opReal v, opVec3 &Du );
void evalDv( opReal u, opReal v, opVec3 &Dv );
void evalNorm( opReal u, opReal v, opVec3 &norm );

virtual csNode* clone( csNode::CloneEnum what);

The surface of the sphere is parameterized by angles, in radians, for latitude and longitude; v corresponds to longitude, u to latitude. Figure 9-11 illustrates the parameterization of an opSphere.

Figure 9-11. Sphere Parameterization

class opSphere : public opParaSurface
{
public:
// Creating and destroying
opSphere( );
opSphere( opReal radius );
~opSphere( );
      
// Accessor functions
void setRadius( opReal radiusVal ) 
opReal getRadius( )              

// Evaluators
void evalPt(   opReal  u, opReal v, opVec3 &pnt );
void evalNorm( opReal  u, opReal v, opVec3 &norm );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
}

x = radius * cos(u) * sin(v)
y = radius * sin(u) * sin(v)
z = radius * cos(v)

opSphere *sphere = new opSphere( 3 );

// under certain conditions, a trim curve is added that keeps only the // portion of the surface above a circle 
if ( nVersions <= 0 )
{
opCircle2d  *trimCircle2d = 
                   new opCircle2d( 1.0, new opVec2(M_PI/2.0,M_PI) );
sphere->addTrimCurve( 0, trimCircle2d );
}
setUpShape( sphere, OP_XDIST*numObject++, Y, OP_VIEWDIST );

The opCylinder class provides methods for describing a cylinder. 

A cylinder can be defined geometrically as the surface in space that is swept by moving a circle along an axis that is perpendicular to the plane of the circle and passes through the center of the circle.

The parameterization of an opCylinder is as follows: u represents the position on the circle and that v represents the position along the axis.

Figure 9-12. Cylinder Parameterization

class opCylinder : public opParaSurface
{
public:
// Creating and destroying
opCylinder( void );
opCylinder( opReal radius, opReal height );
~opCylinder();

// Accessor functions
void setRadius( opReal radiusVal ) ;
void setHeight( opReal heightVal );

opReal getRadius( ) 
opReal getHeight( ) 

// Evaluators
void evalPt(   opReal  u, opReal v, opVec3 &pnt );
void evalNorm( opReal  u, opReal v, opVec3 &norm );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

The opTorus class provides methods to describe a torus. Figure 9-13 illustrates a torus, and how it is parameterized in opTorus.

A torus can be defined geometrically as the surface in space that is swept by moving a circle, the minor circle, through space such that its center lies on a second circle, the major circle, and the planes of the two circles are always perpendicular to each other, with the plane of the minor circle aligned along radii of the major circle. The parametrization of the surface is that u represents a position on the major circle and v represents a position on the minor circle.

Figure 9-13. Torus Parameterization

class opTorus : public opParaSurface
{
public:
// Creating and destroying
opTorus( );
opTorus( opReal majorRadius, opReal minorRadius );
~opTorus();
      
// Accessor functions
void setMajorRadius(  opReal majorRadiusVal )
void setMinorRadius( opReal minorRadiusVal ) 
opReal getMajorRadius( )
opReal getMinorRadius( )

// Evaluators
virtual void evalPt(   opReal  u, opReal v, opVec3 &pnt );
virtual void evalNorm( opReal  u, opReal v, opVec3 &norm );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
}

You can define a cone geometrically by sweeping a circle along an axis in a way similar to the way a cylinder is defined; however, as the circle is swept along the axis, the radius changes linearly with distance.

The parameterization of a point on an opCone is that u measures the angle, in radians, of the point on the circle, and that v measures the distance along the axis from the origin. To truncate a cone, yielding a frustum, adjust the value for v.

Figure 9-14. Cone Parameterization

class opCone : public opParaSurface
{
public:
// Creating and destroying
opCone( void );
opCone( opReal radius, opReal height );
~opCone();
      
// Accessor functions
void setRadius( opReal radius ) ;
void setHeight( opReal height );
opReal getRadius( )              
opReal getHeight( )              

// Evaluators
void evalPt(   opReal  u, opReal v, opVec3 &pnt );
void evalNorm( opReal  u, opReal v, opVec3 &norm ); 

// Copy
virtual csNode* clone(csNode::CloneEnum what);
}

The class opSweptSurface provides methods to describe a general swept surface. Three examples of swept surfaces have been presented: a cylinder, a torus, and a cone. In the first two cases a simple cross-section, a circle of constant radius, was swept along a path. For a cone, the radius of the circle varied according to a simple profile.

To describe a swept surface, you specify a path, a cross section, and a coordinate frame in which the graph of the cross section is drawn at each point on the path. The parameterization of the surface is that u denotes the position along the path and v denotes the position on the cross-section curve. You can also specify a profile, which adjusts the size of the cross-section curve. Thus, for example, with a simple profile method you could generate a sphere from a straight-line path and a circular cross section. Figure 9-15 illustrates the feature of a swept surface.

Figure 9-15. Swept Surface: Moving Reference Frame and Effect of Profile Function

class opSweptSurface : public opParaSurface
{
public:
// Creating and destroying
opSweptSurface( void );
opSweptSurface( opCurve3d *crossSection,
                opCurve3d *_path,
                opCurve3d *_t,
                opCurve3d *_b,
                opScalar  *_profile  );
~opSweptSurface( );

// Accessor functions
void setCrossSection( opCurve3d *_crossSection );
void setPath( opCurve3d *_path );
void setT( opCurve3d *_tng );
void setB( opCurve3d *_b );
void setProf( opScalar  *_profile );
opCurve3d *getCrossSection() ;
opCurve3d *getPath() ;
opCurve3d *getT() ;
opCurve3d *getB() ;
opScalar  *getProf();

virtual void evalPt( opReal u, opReal v, opVec3 &pnt );

virtual csNode clone(csNode::CloneEnum what);
};

As a convenience, the class opFrenetSweptSurface allows you to use the Frenet frame of the path to define the orientation vectors in a swept surface. The Frenet frame is defined by the three unit vectors derived from the tangent, the principal normal, and their cross product. This set of vectors facilitates orienting the cross section perpendicularly to the path at every point.

---

Note: The path for an opFrenetSweptSurface must be at least a cubic to allow for the principal normal calculation, which requires a second derivative.

---

class opFrenetSweptSurface : public opSweptSurface
{
public:
// Accessor functions
opFrenetSweptSurface( void );
opFrenetSweptSurface( opCurve3d *crossSection,
                      opCurve3d *path,
                      opScalar  *profile  );
~opFrenetSweptSurface( );

// Accessor functions
void set( opCurve3d *crossSection,
          opCurve3d *path,
          opScalar  *profile  );

// Copy
virtual csNode clone(csNode::CloneEnum what);
};

// Scalar curve used by the swept surface primitive
static opReal profile( opReal t )
{
return 0.5*cos(t*5.0) + 1.25;
};
 
opCircle3d  *cross = 
                    new opCircle3d( 0.75, new opVec3( 0.0, 0.0, 0.0) );
opCircle3d  *path  = 
                    new opCircle3d( 1.75, new opVec3( 0.0, 0.0, 0.0) );
opFrenetSweptSurface *fswept = 
                      new opFrenetSweptSurface( cross, path, profile );
fswept->setHandednessHint( true );

A ruled surface is generated from two curves in space, both parameterized by the same variable, u. A particular value of u specifies a point on both curves. A ruled surface is defined by connecting the two points with a straight line parameterized by v. The parameterization of the resulting surface is always the unit square in the u-v plane, regardless of the parameterizations of the original curves.

Figure 9-16. Ruled Surface Parameterization

A bilinear interpolation of four points is perhaps the simplest example of a ruled surface, one for which the “curves” that define the surface are in fact straight lines. Thus, you connect two pairs of points in space with lines and then develop the ruled surface. For a bilinear interpolation, the parameterization by u and v is such that, if one of them is held constant, a point “moves” along the connecting straight line at a uniform speed as the other parameter is varied.

class opRuled : public opParaSurface
{
public:
// Creating and destroying
opRuled();
opRuled( opCurve3d *c1, opCurve3d *c2 );
~opRuled();
      
// Accessor functions
void setCurve1( opCurve3d *_c1 );
void setCurve2( opCurve3d *_c2 );
opCurve3d *getCurve1( )   
opCurve3d *getCurve2( )   

// Evaluators
void evalPt(   opReal  u, opReal v, opVec3 &pnt );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

A Coons patch is arguably the simplest surface you can define from four curves whose endpoints match and form a closed loop. Think of the four curves as defining the four sides of the patch, with one pair on opposite sides of the patch defining the top and bottom curves and the other pair defining the left and right curves (see Figure 9-17). The top and bottom curves are parameterized by u, and the left and right curves by v. Thus, u is the “horizontal” coordinate and v the “vertical” coordinate.

The patch is made by

Figure 9-17 illustrates the construction. To understand the result, notice that, after you add the two ruled surfaces, each side of the boundary of the resulting surface is the sum of the original bounding curve and the straight line connecting the bounding curve's endpoints. The straight line was introduced by the construction of the ruled surface that did not include the boundary curve. Subtracting the bilinear interpolation eliminates the straight-line components of the sum, leaving just the original four curves as the boundary of the resulting surface.

Figure 9-17. Coons Patch Construction

class opCoons : public opParaSurface
{
public:
opCoons( );
opCoons( opCurve3d *right,   opCurve3d *left,
         opCurve3d *bottom,  opCurve3d *top );
~opCoons( );

// Accessor functions
void setRight(  opCurve3d *right );
void setLeft( opCurve3d *left );
void setBottom( opCurve3d *bottom );
void setTop( opCurve3d *top );

opCurve3d* getTop();
opCurve3d* getBottom();
opCurve3d* getLeft();
opCurve3d* getRight();

// Surface point evaluator
void evalPt( opReal u, opReal v, opVec3 &pnt );

//Copy
virtual csNode* clone(csNode::CloneEnum what);
};

Just as a NURBS curve consists of Bezier curves, a NURBS surface consists of Bezier surfaces. The set of control parameters is essentially the same for the curves and surfaces: a set of knots, a control hull, and a set of weights. However, for a NURBS surface, the knots form a grid in the coordinate system of the surface; that is, in the u-v plane, and the control hull is a grid of points in space that loosely defines the surface.

Understanding a Bezier surface helps you understand and use a NURBS surface. A Bezier surface is defined essentially as the surface formed by sweeping a Bezier cross section curve through space, along a path defined by a Bezier curve. But, unlike an opSweptSurface, the shape of the cross-section can be changed.

You define a Bezier surface as follows:

A more rigorous discussion appears in the book Curves and Surface for Computer Aided Geometric Design, listed in the section “Recommended Background Reading”.

A NURBS surface joins Bezier surfaces in a smooth way, similar to NURBS curves joining Bezier curves. The class opNurbSurface provides methods to describe a NURBS surface.

class opNurbSurface : public opParaSurface
{
public:
// Creating and destroying
opNurbSurface( void );
~opNurbSurface( void );

// Accessor functions
void setControlHull( int iu, int iv, opVec3 &p );
void setControlHull( int iu, int iv, opVec4 &p );
void setWeight( int iu, int iv, opReal w ); 
void setUknot( int iu, opReal u );
void setVknot( int iv, opReal v );
void setControlHullUSize( int s );
void setControlHullVSize( int s );

// Get the same parameters
opVec3& getControlHull( int iu, int iv) ;
int     getControlHullUSize( void );
int     getControlHullVSize( void );
opReal  getWeight( int iu, int iv)
opReal& getUknot( int iu);
opReal& getVknot( int iv);
int     getUknotCount( void );
int     getVknotCount( void );
int     getUorder( void ) ;
int     getVorder( void ) ;
void    removeControlHullElm(int ui, int iv);
void    removeUknot(int iu);
void    removeVknow(int iv);
void    flipUV();

// Evaluator
virtual void evalPt(   opReal u, opReal v, opVec3 &pnt );
virtual void evalDu(   opReal u, opReal v, opVec3 &Du );
virtual void evalDv(   opReal u, opReal v, opVec3 &Du );
virtual void evalNorm( opReal u, opReal v, opVec3 &norm );

int getMemSize();
virtual csNode* clone(csNode::CloneEnum what);
};

Methods in opNurbSurface

The member functions are essentially the same as those for opNurbCurve3d (see “NURBS Curves in Space”), however:

• The hull is a grid of opVec3s indexed by i and j.

• The set of knots is defined by points on the u and v axes.

• There are B-spline basis functions (of possibly differing orders) associated with each coordinate direction.

---

Note: opNurbSurface redefines the virtual evaluators inherited from opParaSurface for tangent and normal vectors; the methods use the NURBS equation rather than finite, central differences.

---

Equation Used to Calculate a NURBS Surface

Indexing is determined by the following equation that OpenGL Optimizer uses to calculate a NURBS surface (the index i corresponds to iu in the API, and j corresponds to iv):

Alternative Equation for a NURBS Surface

A NURBS surface can also be developed from the following alternative expression:

For this case, you must change the coordinates of the control points to get the same surface from OpenGL Optimizer. You convert the coordinates of the control points from (x,y,z,w) to (wx,wy,wz,w).

int i, j;

opNurbSurface *nurb = new opNurbSurface;

// Control hull dimensions
#define USIZE 4
#define VSIZE 5

// Set up the control hull size because we know a priori how big
// the nurb is.  The next two lines are used for space
// efficiency but are functionally unnecessary.
nurb->setControlHullUSize(USIZE);
nurb->setControlHullVSize(VSIZE);

// Make the control hull be an oscillating grid
for ( i = 0; i < VSIZE; i++ )
{
opReal y = i/(float)(VSIZE - 1) * 2*M_PI - M_PI;
      
for ( j = 0; j < USIZE; j++ )
{
opReal x = j/(float)(USIZE - 1) * 2*M_PI - M_PI;
opReal val = 6*pow( cos(sqrt(x*x + y*y)), 2.0); 

// Make the control hull a box, j maps to u and i maps to v
nurb->setControlHull( i, j, opVec3( x, y, val));

// Add the weights
nurb->setWeight( i, j, 1.0 );
}
}

// Add the knot points
nurb->setUknot( 0,  0.0 );
nurb->setUknot( 1,  0.0 );
nurb->setUknot( 2,  0.0 );
nurb->setUknot( 3,  0.0 );
nurb->setUknot( 4,  1.0 );
nurb->setUknot( 5,  1.0 );
nurb->setUknot( 6,  1.0 );
nurb->setUknot( 7,  1.0 );

nurb->setVknot( 0,  0.0 );
nurb->setVknot( 1,  0.0 );
nurb->setVknot( 2,  0.0 );
nurb->setVknot( 3,  0.0 );
nurb->setVknot( 4,  1.0 );
nurb->setVknot( 5,  1.0 );
nurb->setVknot( 6,  1.0 );
nurb->setVknot( 7,  1.0 );

// Only trim reps in the first row
if ( nVersions <= 0 )
{
// Add a super quadric trim curve 
opSuperQuadCurve2d  *trimCircle0 = new opSuperQuadCurve2d( 0.25, new opVec2(0.25, 0.50), 2.0 );
nurb->addTrimCurve( 0, trimCircle0, NULL );

// make a 4-th order nurb trim curve
opNurbCurve2d *l = new opNurbCurve2d;

l->setKnot(0,0.0);
l->setKnot(1,0.0);
l->setKnot(2,0.0);
l->setKnot(3,0.0);
l->setKnot(4,1.0);
l->setKnot(5,1.0);
l->setKnot(6,1.0);
l->setKnot(7,1.0);
l->setControlHull(0,opVec2(0.50,0.50));
l->setControlHull(1,opVec2(0.90,0.10));
l->setControlHull(2,opVec2(0.90,0.90));
l->setControlHull(3,opVec2(0.50,0.50));

nurb->addTrimCurve( 1, l, NULL  );
}

Hermite-spline surfaces interpolate a grid of points; that is, they pass through the set of specified points under the constraint that you supply the tangents at each point in the u and v directions and the mixed partial derivative at each point. This surface definition is the natural generalization of Hermite-spline curves, discussed in “Hermite-Spline Curves in the Plane”.

Figure 9-19. Hermite Spline Surface With Derivatives Specified at Knot Points

Hermite-spline surfaces are made of Hermite patches (see Figure 9-19). A bicubic Hermite patch expands the definition of a bilinear interpolation to include specification of first derivatives and mixed partial derivatives of the surface at each of the four corners. The adjective “bicubic” in the name of the patches refers to the mathematical definition, which includes products of the cubic Hermite polynomials that define a Hermite-spline curve.

An advantage of including the derivatives to constrain the surface is that it is simple to combine the patches into a smooth composite surface, that is, into a Hermite-spline surface. A more formal discussion of these objects appears in the book Curves and Surface for Computer Aided Geometric Design listed in the section “Recommended Background Reading”.

class opHsplineSurface : public opParaSurface
{
public:
// Creating and destroying
opHsplineSurface();
opHsplineSurface( opReal *_p, 
                  opReal *_tu,  opReal *_tv, opReal *_tuv, 
                  opReal *_uu, opReal *_vv, 
                  int uKnotCount, int vKnotCount );
~opHsplineSurface();

// Accessor functions
opVec3& getP( int i, int j );
opVec3& getTu( int i, int j );
opVec3& getTv( int i, int j );
opVec3& getTuv( int i, int j );
opReal  getUknot( int i );
opReal  getVknot( int j );
int     getUknotCount();
int     getVknotCount();
opReal  getCylindrical();

void setAll( opReal *p, 
             opReal *tu,
             opReal *tv, 
             opReal *tuv,
             opReal *uu,
             opReal *vv,
             int uKnotCount,
             int vKnotCount );
void setCylindrical(opReal cylinderical);

// Surface point evaluator
void evalPt( opReal u, opReal v, opVec3 &pnt );
// Copy
virtual csNode* clone ( csNode::CloneEnum what);
};

The class has the following main methods:

class opCuboid : public opRep
{
public:
// Creating and destroying
opCuboid( );
opCuboid( opReal width, opReal height, opReal depth );
~opCuboid();

// Accessor functions
void setWidth( opReal widthVal);
opReal getWidth( )

void setHeight( opReal heightVal );
opReal getHeight( )

void setDepth( opReal depthVal );
opReal getDepth( );

// Copy
virtual csNode* clone(csNode::CloneEnum what);
};

opDisSurface is the base for the all discrete surfaces and, more generally, higher-dimensional meshes. A discrete surface is described as a set of discrete points interconnected by a specific topology. An example of such a topology is a planar grid structure. The base class provides methods only for discrete trim curves.

The opRegMesh template class describes a vector-valued function over a rectangular mesh. Thus, an opRegMesh is the natural object for visualizing many data sets or scientific modeling calculations.

The type of the template is determined by the return value of the mesh function you define. For example, you can describe a discrete surface with a two-dimensional grid and a mesh function that returns csVec3f positions of points on the surface. Thus the mesh would be of type csVec3f. A surface tiling is developed by the member function evalPt(), which interpolates values of the mesh function.

A mesh can have an arbitrary number of dimensions, although opRegMesh provides special operations for two-, three-, or four-dimensional meshes. A mesh can have regular or variable spacing in all dimensions. In general, if you specify a mesh by an array of grid points, then the argument of the mesh function must be the same data type as the grid points.

template <class T>
class opRegMesh : public opDisSurface
{
public:
opRegMesh( );
opRegMesh( int Xres, int Yres );
opRegMesh( int Xres, int Yres, int Zres );
opRegMesh( int Xres, int Yres, int Zres, int Tres );
opRegMesh( int d,    int *res );
   
~opRegMesh( );

// Set and get the dimensionality of the mesh
void setDim (int _dim)  
int  getDim (void)      

// Set and get the dimension of the mesh
void setRes( int Xres, int Yres );
void setRes( int Xres, int Yres, int Zres );
void setRes( int Xres, int Yres, int Zres, int Tres );
void setRes( int d,    int *res );

int *getRes( ) ;

// Set and get the type
void setType( opRegMeshType meshType ) 
opRegMeshType getType( )

// Set and get the origin
void setOrigin( opReal *Origin )
opReal& getOrigin( ) 
 
// Set and get the delta spacing
void setSpacing( opReal *Delta );

opReal& getSpacing( ) ;
opReal& getSpacing( int i )
opReal& getSpacing( int i, int j );
opReal& getSpacing( int i, int j, int k );
opReal& getSpacing( int i, int j, int k, int l );
   
// Arbitary indexing via an index vector
opReal& getSpacing( int *index );

// Set and get the mesh function
void  setFunction( T *function ) 
T    *getFunction( )        

// Set and get variable spacing grid 
// (memory maintained by calling program
// assumes sizeof(grid) = 
//                  ndim*sizeof(opReal) * res[0]*res[1]*...*res[ndim-1]
void     setGrid (opReal *_grid) 
opReal  *getGrid ()              

// Single index subscripting operator
T& operator[]( int i );

// One, two, three and four dimensional indexing operators
T& operator()( int i );
T& operator()( int i, int j );
T& operator()( int i, int j, int k );
T& operator()( int i, int j, int k, int l );
   
// Arbitary indexing via an index vector
T& operator()( int *index );
   
// Point interpolated evaluators
void evalPt( T& pt, opReal x, opReal y );
void evalPt( T& pt, opReal x, opReal y, opReal z );
void evalPt( T& pt, opReal x, opReal y, opReal z, opReal t );   
// Extract positional information out of grid
opReal  gridVal ( int i );
csVec2f gridVal ( int i, int j );
csVec3f gridVal ( int i, int j, int k );
csVec4f gridVal ( int i, int j, int k, int l );

// Can set extents if you know them, or compute them
void setExtents (T _min,  T _max);
void getExtents (T *_min, T *_max);
// compute min/max over all data points
bool computeExtents (bool force);
};

make_data_cube (&data, dims);
ndim = 3;
      
...

//     Set origin and mesh spacing
opReal orig[3]  = ;
opReal delta[3] = ;

// --- Allocate opRegMesh to contain raw data
opRegMesh<opReal>  *rm = new opRegMesh<opReal>

//     Load parameters of the opRegMesh rm:
rm->setType (opConstant);
rm->setRes (ndim, dims);
rm->setDim (ndim);  
// do this after setRes because setRes(d,res) will reset dim=4
rm->setOrigin (orig);
rm->setSpacing (delta);

//     Load function values:
rm->setFunction (data);

densityMesh = new opRegMesh<opReal>;

densityMesh->setType (opVariable);
densityMesh->setDim (3);
densityMesh->setRes (dims[0], dims[1], dims[2]);
densityMesh->setOrigin (orig);
densityMesh->setGrid (grid);
densityMesh->setFunction (real_rho);

momentumMesh = new opRegMesh<csVec3f>;

momentumMesh->setType (opVariable);
momentumMesh->setDim (3);
momentumMesh->setRes (dims[0], dims[1], dims[2]);
momentumMesh->setOrigin (orig);
momentumMesh->setGrid (grid);
momentumMesh->setFunction (momentum);