Al, who holds a degree in systems design, is a developer for Rogue Wave Software. He can be contacted P.O. Box 2328, Corvallis, OR 97339 or at alv@roguewave.com.
Object-oriented programming is emerging as the dominant software methodology of the '90s. Even the traditional scientific-computing community is emerging from its Fortran shell to investigate the promise of C++ and other OO languages. Why? Because object-oriented languages allow you to code at a higher level of abstraction than traditional, procedural languages. Consequently, object-oriented code ends up being easier to read, write, maintain, and debug than traditional C or Fortran code.
Coding at a higher level of abstraction is, of course, relative. Many of us can remember when almost all scientific programs were coded in assembler (or even directly, via switches). In those days, numerical programmers had to translate mathematics down to the bit, byte, and register level. Then along came Fortran, and you could code formulas involving real and complex numbers directly, without having to translate to assembler. Today, object-oriented languages like C++ let you create even higher-level abstractions, and treat them as if they were built-in data types. You can, for example, model vectors and matrices, functions, and nonlinear systems of equations. Then you can use these abstractions (or classes) as building blocks for still-higher levels of abstraction.
In this article, I'll discuss code developed at the University of Waterloo to research spline-based, finite-element models. To keep the the discussion focused, I'll concentrate on the specific problem of determining the load-deflection behavior of a simple bicycle-wheel model like that shown in Figure 1.
The goal of the application is to simulate how a bicycle wheel behaves when force is applied to it--the stress in the spokes, the shape of the wheel, and the force it can withstand before buckling. All of these can be computed once you know how the applied force deforms the wheel. Mathematically, the deformation of a structure due to a load is modeled by partial differential equations; the equations used to model the spokes are shown in Figure 2. Unfortunately, for all but the simplest of cases, the exact solutions to the differential equations aren't known, and so we must be satisfied with an approximate solution.
One way of constructing an approximate solution is the finite-element method. I'll explain the key points of the method using a simple example: exerting downward force on a yardstick on end; see Figure 3. As you push on the yardstick, it compresses and the tip moves down slightly. As you continue to push harder and harder, the yardstick eventually bows outward, and the tip moves down much more easily. This phenomena is called buckling. Since most structures cannot withstand buckling, an engineer would often like to know the precise relationship between weight and maximum deflection; this relationship can be plotted on a force-deflection diagram like the one in Figure 3. (As you can see from this diagram, the yardstick buckles at about 16 pounds of force; feel free to test this out at home). Ideally, you'd create this plot by directly solving the differential equations for each different weight, except that this is usually too difficult. It turns out that for this simple example you can solve the equations exactly, but let's ignore that for the time being and focus on the finite-element method.
The key idea behind the finite-element method is to restrict the possible deformations to a simple form that approximates the true deflection. For example, you might restrict the deformation of the yardstick to be a cubic polynomial. In the variant of the finite-element method we're interested in, the deformation is restricted to be a piecewise polynomial spline function. A spline, as shown in Figure 5, is a smooth curve whose shape is determined by control points (the Xs in Figure 2). You can use a mathematical technique (such as the Rayleigh-Ritz or Bubnov-Galerkin method) to convert the governing equations to a new set of equations with the spline control points as unknowns. Looking at the yardstick example, you can represent the deflection of the yardstick as a spline with, say, five control points. The Rayleigh-Ritz method is used to convert the governing differential equation into a nonlinear system of equations with the control points as unknowns.
You've now approximated the governing system of partial differential equations with a system of nonlinear algebraic equations whose unknowns are control points in a spline representing the deflection. Unfortunately, these nonlinear equations are still too complicated to solve directly, so you need to use some iterative solution technique like Newton's method. To apply such a technique, you need to be able to compute two things:
If you're familiar with the finite-element method, you may be wondering why we're interested in using splines to approximate the deformation, rather than the conventional approach of breaking the structure down into small pieces (elements), assuming polynomial deflection on each element, and then stitching the elements together. One reason is efficiency: The use of splines can lead to considerably more efficient solutions, giving better accuracy for fewer degrees of freedom.; this can be shown both in practice and in theory. A second reason for preferring splines is simplicity: In computer-aided design systems, the geometry of structures is often already defined using splines. If so, why not use that geometry to guide the deformation analysis, rather than having to discretize the structure into elements?
Still, a drawback of the spline-based method as opposed to the conventional element-based method is that the spline formulation leads to more complex equations for the out-of-balance forces and tangent-stiffness matrix--and one of the reasons for using object-oriented programming is to help manage this complexity.
If you recall, the goal in this article is to construct the load-deflection curve for the two-dimensional bicycle-wheel model shown in Figure 1. The wheel consists of a hub, a rim, and 36 spokes. The hub is considered to be rigid--no amount of force will deform it--and fixed in space. The rim and spokes are modeled as nonlinear, thin, shallow beams. (The adjectives refer to assumptions made in deriving their governing equations.) One end of each spoke is fixed in space by the hub; the other is pin-jointed to the rim. The wheel is loaded by a point force pushing in on the rim. This force could represent the force of a wall smashing into the wheel, or the floor pushing up on the wheel. (If you're concerned about a world where moving walls collide with stationary wheels, note that by changing reference frames you get the physically identical, but more intuitive, situations of a moving wheel hitting a stationary wall, or a rider's weight pushing the wheel against the ground.)
Deflections of the components which deform--the rim and spokes--are determined using the spline finite-element method. To apply this method, you need to be able to construct an out-of-bounds force vector and tangent-stiffness matrix for each component. The complete formulas for the components of these vectors and matrices are long and complex. Still, the formula for one component in the spoke's tangent-stiffness matrix is shown in Figure 4. This formula is typical: An entry in the matrix is obtained as the integral of an expression involving the current deformation function and the spline-basis functions. I'll describe shortly how to write C++ code to compute this expression.
For a structure (like the bicycle wheel) with more than one part, the finite-element models for each individual part must somehow be tied together into a model for the whole. This is done by imposing constraints on the control points of the individual pieces so that the deformation is consistent between pieces according to how the pieces are joined. In the case of the wheel, the end-control point of each spoke indicates the deformation at that end--this control point is constrained to be equal to the deformation of the rim where the spoke is attached. The constraints between all the spokes and the rim are represented using a sparse matrix; this matrix is used to amalgamate the out-of-balance force vector and tangent-stiffness matrix of the component parts into a global out-of-bounds force vector and tangent-stiffness matrix.
Now comes the tough part: choosing which abstractions (classes) to use to model the problem. You do this in three stages: First try to identify components of other software packages which you can reuse in this design; then create components custom built for this type of application; and finally build any application-specific classes you need, and put the pieces together into an application. The key components of the design are shown in Figure 5. You may find it helpful to occasionally refer to this figure to see where the modules described fit in.
Perhaps the most fundamental module needed by any scientific code are classes for storing and manipulating arrays of numbers. Neither C nor Fortran provides good abstractions for arrays: In C, you can't treat an array as a first-class object, and Fortran lacks the capabilities for dynamic resizing of arrays. Neither language lets you express operations on arrays in a high-level fashion: You must use explicit loops to do even basic computations. These drawbacks are easily remedied using C++ classes. Operator overloading and a well-thought-out class interface provide excellent expressive capabilities. In fact, even if you use no classes beyond a set of vector, matrix, and array types, a good array class library can make the transition from C or Fortran worthwhile.
When work on the spline finite-element project began in 1988, there were only a handful of people doing numerics in C++; consequently almost no class libraries existed. Once we began writing and using our own basic vector classes, one of the truths of C++ became clear: While using C++ classes is easy, writing them is not. It demands an excellent knowledge of the language and more skill than we had at the time. One of the few other C++ numerical projects was the Data Analysis and Interactive Modeling System (DAIMS) at the University of Washington's Oceanography department, led by Thomas Keffer. They, too, had been struggling with building vector and matrix classes. After a lengthy exchange of ideas and experiences over the Internet, ideas from our library were incorporated into the DAIMS classes. The result was a high-performance, flexible library of vector and matrix classes. Tom went on to found Rogue Wave software, and the vector and matrix classes evolved into the core of Math.h++, the first commercially available C++ class library.
The concept of a function is fundamental to mathematics, and yet the only support available for this concept, in either C or Fortran, is the rather primitive pointer-to-function data type. Just to be sure we're on the same wavelength, I'm talking here about mathematical functions like f(t), not subroutines in a computer program. In our finite-element code, functions are used to describe just about all aspects of the problem: geometry, material properties, loading, and deformation. Since no function class library existed back in 1988, we developed classes to support the notion of univariate function, various special types of functions like polynomials and splines, and things related to functions, such as integrators.
The most important class in this module is Func, which embodies the idea of a univariate function. The key part of the header for Func is shown in Listing One (page 91). (In the interest of clarity, member functions unnecessary for this application have been omitted.) By overloading the function-call operator, (), to do evaluation, using instances of Func becomes very natural. Consider, for example, the simple plot subroutine Listing Two (page 91) that uses the Func class. Func is an example of an abstract (pure virtual, in C++ talk) class. You can't instantiate a Func itself, since Func represents an abstract idea, not a concrete one; instead you must instantiate a class derived from Func with an actual evaluation method. Two classes derived from Func are needed in the bicycle-wheel application: Poly to represent polynomials, and NUBFunc to represent a nonuniform B-spline function.
Besides evaluating, you need to be able to combine functions in expressions, like f*g+h. The problem is that in order to evaluate this expression later, you need to keep references to its constituent functions: f, g, and h. The solution is to use a smart pointer class, FuncHandle, for referring to functions on the heap. The FuncHandle takes care of reference-counting the functions and deleting them when they are no longer needed. Once you have handles referring to functions, the expression classes can keep handles to their constituent functions, and we can use expressions. Example 1 shows how you might use this to plot the function sin(x)cos(x): (Note that the CFunc class is used to encapsulate subprogram style functions.)
FuncHandle s(new CFunc (sin)); FuncHandle c(new CFunc (cos)); plot (s*c);
To form our finite-element, out-of-balance force vectors and tangent-stiffness matrices, and for many other applications throughout science and engineering, you need to be able to integrate functions. One way to think of integration is as a functional: a function of a function. The Integrator class embodies this idea: It is an abstract base class, like Func, but to evaluate it requires not a real number, but a function. The GQuad class is a subclass of integrator for doing integration using Gauss quadrature. Listing Three (page 91) shows how to use the integrator and the ability to form expressions of functions to compute an entry in a finite-element tangent-stiffness matrix using the formula in Figure 5. The code is beautiful: It mirrors almost exactly the mathematics: Think about how this would look in Fortran!
Often the structures you're interested in (including our bicycle-wheel model) are governed by a system of nonlinear equations. Unlike linear equations, nonlinear equations usually cannot be explicitly inverted, so you rely on some sort of approximation technique to solve them. The, nonlinear equations module provides classes for representing and solving nonlinear equations. Although the classes in this module were carefully designed to be useful in other contexts, I'll use structural analysis words, like "forces" and "displacements," to fit in with the rest of this article.
A general system of nonlinear equations can be written as f[i],(x[j]), where (in our context) x[j] are control points representing displacement, and f[i] are out-of-balance forces at that displacement. We are usually interested in the equilibrium position of a structural system, which is the displacement where there are no out-of-balance forces, f[i](x[j]) = 0. The nonlinear equation is represented by the abstract base class NLSys ( Listing Four, page 91). (In the interest of clarity, member functions unnecessary for the application have been omitted.) There are two things which we require of an NLSys object: Given a vector of variables x[j], we must be able to evaluate the out-of-balance forces, f[i], and solve the system of equations J[x]=f where J is the Jacobian matrix (tangent-stiffness matrix) of the system.
The point of representing a nonlinear system using the NLSys class is to be able to solve systems of equations and plot load-deflection curves. Since these sorts of algorithms have many parameters, you encapsulate them into objects. The NRSolve class uses Newton-Raphson iteration to find a solution to the equations given an initial guess. The YSContinue class uses Yang and Shieh's continuation method to compute a complete load-deflection curve.
At this point, the list of non-application-specific classes you could buy, beg, borrow, steal, or (as a last resort) write is just about exhausted. The next step is to use them to write classes for the domain of interest: finite-element structural analysis. These domain-specific classes are typically a little easier to write than the domain-independent classes, since you don't have to worry quite so much about other people using them for purposes you hadn't thought of.
Most of the domain-specific classes in the spline finite-element project represent types of structural models. Although we've written a large number of these classes, I'll consider only the two needed for the wheel: NLBeam and NLArch. NLArch represents a finite-element model of a shallow, thin, nonlinear curved beam. We'll use it to model the wheel's rim. The NLBeam class is a special case of NLArch where the curvature is allowed to be infinite, thus making it a straight beam. This class is used to represent the spokes. The job of these classes is to represent the nonlinear system of equations governing the control points in the spline approximation.
A domain-specific class needed (other than a structural component) is the StructureConstraint class. This is used to represent the constraints between substructures in a model. In this case, you need this class to act as the glue between the spokes and the rim.
With all the component parts in hand, all you need to do is snap them together to build the application, An instance of the bicycle-wheel class itself, Wheel, contains an array of 36 spoke objects (of class NLBeam), one NLArch object representing the rim, and a StructureConstraint object to tie the spokes to the rim.
The final step is to use the wheel model to generate a load-deflection curve. First you need to interface the Wheel object to the nonlinear equation solver module's NLSys object. This is done using multiple inheritance: A new object, WheelSys, is defined which inherits from both Wheel and NLSys. This is a common, simple way to tie pieces of an application together in C++. The WheelSys class uses a Double PDFact object internally to maintain a factorization of the current tangent-stiffness matrix; this is used to implement the NLSys solveTangentSystem() member function. The response of the wheel to a head-on collision force is shown in Figure 7, where you can see how strong a well-built bicycle wheel can be--buckling does not even begin to occur until past 25,000 pounds of force have been applied.
_OBJECT-ORIENTED FINITE-ELEMENT SOFTWARE_ by Al Vermeulen[LISTING ONE]
class Func {
public:
virtual ~Func() {};
virtual double evaluate(double x) const =0;
virtual double evaluate(double x, int d) const;
double operator()(double x) const {return evaluate(x);}
double operator()(double x,int d) const {return evaluate(x,d);}
virtual DoubleVec operator()(const DoubleVec& x) const;
virtual DoubleVec operator()(const DoubleVec& x, const IntVec& d) const;
virtual DoubleVec operator()(double x, const IntVec& d) const;
};
[LISTING TWO]
void plot(const Func& f, double a, double b, double step)
{
moveto(a,f(a));
for(double x=a+inc; x<=b; x+=step) {
drawto(x,f(x));
}
drawto(b,f(b));
}
[LISTING THREE]FuncHandle dudx = diff(u,x); FuncHandle dwdx = diff(w,x); FuncHandle dNidx = diff(Ni,x); FuncHandle integrand = EA*(dudx+dwdx*dwdx/2)*dNidx; double fi = I(integrand); // I is an instance of the GQuad class[LISTING FOUR]
class NLSys
{
public:
virtual DoubleVec value(const DoubleVec& x) =0;
// Calculate K(x).
virtual int setKt(const DoubleVec& x) =0;
// Evaluate the tangent stiffness matrix at the point "x".
// Subsequent calls to solveTangentSystem should use the
// tangent stiffness matrix as set here. If the tangent
// stiffness matrix is singular at this point then return a
// 0 and do nothing, otherwise return a 1.
virtual DoubleVec solveTangentSystem(const DoubleVec& b) =0;
// Solve the system Ku=b for u where K is the last calculated
// tangent stiffness matrix. u and b are increments.
};
Copyright © 1993, Dr. Dobb's Journal