Building Finite Element Models With Mathematica

C.J. Purcell
Defence Research Establishment Atlantic
PO Box 1012, 9 Grove St., Dartmouth, N.S., Canada B2Y3Z7

Abstract

The goal of this work was to assess the utility of Mathematica for finite element analysis, specifically in the areas of model generation and results display. This document will describe a prototype Mathematica package that can be used to create finite element (FE) models, both 2D and 3D. The application area is sonar transducer design, and the target finite element code is MAVART3D.  However, the package will readily generalize to other FE codes, for example, the node and element descriptions are not hardwired, but are user definable. Transducer models are usually simple enough that they can be built by command files, with 3D visualization most useful for debugging, and results display. Mathematica is well suited to this method of model building. A number of commercial mesh generators use a WYSIWYG approach where changes to a model often require starting from scratch. In the parametric approach used here, the model is  described by parameters contained in a notebook and revisions to a design are easily made.

1.  Introduction

Starting in 1976, Defence Research Establishment Atlantic (DREA) began to develop a finite element code called MAVART, which is an acronym for Model for the Analysis of the Vibrations And Radiation of Transducers. The main purpose of MAVART is to analyze the performance of piezoelectric underwater acoustic transducers, like hydrophones and sonar projectors. The first version of the code had an axisymmetric capability only but it played an important role in DREA's development of a number of sonar projectors over the last 20 years, including the free flooding ring, ringshell, resonant pipe and barrel stave projector (BSP). The BSP violates the symmetry assumption such that a fully 3D model is needed to predict its performance, and this led to development of MAVART3D, over the period 1991-96.

During the 1980's DREA developed model generation tools for MAVART, written in Fortran for mainframes and Tektronix graphics terminals.  These tools could not readily be extended to 3D. The solution adopted was to lease a commercial off the shelf modeller, PATRAN [1], which ran on Silicon Graphics workstations. Translators had to be written to interface PATRAN with MAVART. While PATRAN met our requirements, it tied us to a limited range of expensive computer platforms, just when it was becoming apparent that MAVART could run effectively on the increasingly powerful personal computers. There was no low cost (even if that meant reduced capability) way to use PATRAN, i.e., you couldn't use it on a PC or laptop,  it was high-end workstation or nothing. While this had little impact on DREA's in-house work, it inhibited transferring MAVART3D to industry, where modest computing resources are the norm. This provided the incentive to develop an affordable modeller for MAVART3D, because without it, the code would probably never leave the lab. After reviewing the features of a number of modellers a wish list for the bare essentials was drawn up.

Parametric modelling quickly rose to the top of the list. Parametric modelling means that the part dimensions, finite element meshing instructions, boundary conditions, etc., can be specified by parameters in a file of commands, the file can be processed by a command interpreter and new models built by the computer with minimal user intervention.  Modellers that are purely mouse-driven, or WYSIWYG (what you see is what you get) in style, are fine for making a single model of a given design. However they are next to useless when a new design is being optimized and numerous modifications are required. With a WYSIWYG modeller, to make a change of shape you often have to go back to the beginning and start clicking again. A hybrid approach is used by some modellers where the model can be built with mouse clicks, but in the background user input is converted into commands which are logged in a session file. This file is all that is needed to be stored, since the complete model can be generated from it, and it is usually tiny in size compared to the finite element model. With some computer operating systems, a shell script can instruct the modeller to process a session file into a finite element model, do the translation, finite element analysis, postprocessing, and return a perfomance measure that can guide the next edit of the session file. A high-end system like Pro/ENGINEER [2] can do all this from within the application.

After getting comfortable working with session files one realizes that interactivity is not essential during the process of creating simple models, such as are needed for sonar transducer design. Working from a hand-drawn sketch with occasional plotting of the preliminary stages of the model can be quite efficent. Truly interactive graphics, with live rotations and animations are most useful when debugging the model, and showing off results to the boss or clients.
 
With these considerations, the wish list for a 3D modeller became:
 1. parametric modelling
 2. command language driven - with GUI to be added later if affordable
 3. graphic capable, but real-time not essential
 4. complex variables built-in
 5. portable - should run on anything - MAC, PC, UNIX, ALPHA...
 6. all data files ASCII text only, e-mailable
 7. easy to maintain and modify
 8. uses standard file types (e.g., Postscript) for graphics.
 
With these goals established, the next step was to identify a suitable software development environment. After our first look at it in one of Nancy Blachman's courses,  Mathematica [3] was the obvious choice. The biggest unknown was whether the number crunching performance would suffice and the only way to check this was to prototype. Mathematica  is a great prototyper, but how far beyond prototyping can it be driven? Finite element modelling is definitely getting into the realm of serious number crunching, and therefore an interesting test of Mathematica's  capabilities. Our prototype finite element modelling package has the working title ModelMaker. 

The design of the package is described in Section 2.  In particular, Section 2.1 introduces points, and developes the database used to store them. This same database structure is used for all the other entities: curves, surfaces, solids, node types, element types, nodes, elements and materials, that are presented in subsequent parts of Section 2. Examples which illustrate possible applications are shown in Section 3. Comments on the way ahead for ModelMaker and conclusions are given in Sections 4 and 5.

2.  The ModelMaker Package

The ModelMaker package is based on standard finite element methodology -  first the geometrical description of a part is defined, then the geometry is meshed to create nodes and elements, and boundary conditions and material properties are applied to the mesh. The model is then written out for processing by the external analysis code. The results of analysis are read back for post processing and display of deformations, etc. The geometrical description of a part is defined by lists of entities: points, curves, surfaces and solids. The naming convention for entities uses arbitrary integer identifiers (IDs). Points are the only entities that need contain Cartesian coordinates. Curves are defined by giving a list of points and the name of an interpolation function, similarly surfaces are defined by a list of curves that define the cross sections of the surface and an interpolation function, and solids by a list of cross-sectional surfaces and an interpolation function. The default interpolation method for curves, surfaces and solids uses cubic parametric equations, but other functions can be substituted.

2.1  Points

The key software component of a modeller in terms of performance is the database. Mathematica  can be used to build an object oriented programming (OOP) environment [4] which was investigated for its applicability to managing the data. While this worked, it made debugging hard because the OOP approach hides the database from the user.  An easy way to build a random access database in Mathematica [5] is to assign a value to a function corresponding to each entry in a database.  For example, a function points can be assigned a value (here a list) corresponding to the ID  = 1, with 

points[1]  =  {x,y,z}; 

 and next 3 points with arbitrary IDs are defined.

Clear[points,x,y,z]; (* erase any old values *)
points[1] = {0,0,0}; (* definitions of points 1,4,7 *)
points[4] = {1,0,0};
points[7] = {x,y,z};

If we now ask for points[4] we recover the data corresponding to the ID of 4 with no surprises:

points[4]

Note what happens with the last point:

points[7]

where the coordinates remain as symbols whose values have not yet been defined. Later we can assign values to x,y,z.

{x,y,z} = {1,1,-1};

And if we now ask for points[7] we see an interesting possibility:

points[7]

Here the location of a point has been set after its logical definition. This could be useful, but even better, the definition for points[7] remembers its parentage. If we make yet another assignment:

{x,y,z} = {x-1,y-1,-z};

and ask once more for

points[7]

the coordinates have changed again. This means the database can contain symbols whose values we can change to modify, or morph a structure. You can find the complete definition for points by asking:

??points

where the symbols x,y,z are again revealed.

Next we need an easy way to locate points. One method suggested in [5] is to maintain a separate list of IDs, which in the example so far could be something like:

pointlist = {1, 4, 7};

To link these IDs to the point data, another field can be prepended to the point data record, containing a pointer back to pointlist. The next stage in the evolution of our database is a record structure illustrated by :

Clear[points];   (* erase any old values *)
points[1] = {1,{0,0,0}}; (* make definitions *)
points[4] = {2,{1,0,0}};
points[7] = {3,{x,y,z}};

Now to get all the points in the point database we can Map the function points over all the IDs in pointlist.

Map[points,pointlist]

So far we have just considered integer IDs but symbolic IDs would work fine.

points[center] = {4,{1,1,1}};
AppendTo[pointlist,center]; (* add the new id *)
(* now do another but find the next pointer automatically *)
points[top] = {1+Length[pointlist],{1,0,1}};
AppendTo[pointlist,top];    (* add another new id *)
Map[points,pointlist]

and we can extract the coordinates of the point named center just as for the integer IDs.

points[center]

because pointlist can contain symbols just as well as integers

pointlist

The process of creating new points thus is a definition, of the form 

points[id]  =  {1+Length[pointlist],data} 

followed by:

AppendTo[pointlist,id] 

to update the list of points.

Having created points we can plot them with something like:

Show[Graphics3D[{PointSize[.1],
Map[Point,Last[Transpose[Map[points,pointlist]]]]}]];

Now that we've created some points how do we get rid of them?  Let's say we want to delete the point named center.  We first must clean up pointlist using the built-in function Drop as shown:

pointlist = Drop[pointlist,{points[center][[1]]}]

Next we delete the definition points[center] :

points[center] = .

Where this starts to get messy is that all the pointers to the remaining points, in this example the pointer belonging to the point named top, must be updated.

points[top]

Note that points[top][[1]] is 5, but it must be updated to 4 to account for the deleted entry. When the database is large this becomes a time consuming operation, and not a nice database design.

To avoid this problem we consider another database design  - what for Mathematica  is obvious  - namely base everything on lists.  Mathematica manipulates lists using the functional programming paradigm much faster than if procedural constructs are used, so it is natural to consider a definition for points not as a function, but as a list. We no longer need to create and update pointlist, and can make a further simplification by discarding the pointers to pointlist. This is illustrated with the following examples:

points = {}; (* initialize as an empty list *)
AppendTo[points, {1,{0,0,0}}];(* make definitions as lists *)
AppendTo[points, {4,{1,0,0}}];
AppendTo[points, {7,{x,y,z}}];
AppendTo[points, {center,{1,1,1}}];

points

The concept of a pointlist as a list of all IDs is still useful but now it can be generated automatically from the list points using the following built-in routine:

pointlist := Transpose[points][[1]];

pointlist

Another tool that will be useful is a way to find the pointer that corresponds to a particular ID:

pointindex[id_] := Position[pointlist,id][[1,1]];

Using pointindex we discover where the data belonging to center are stored in the array points.

pointindex[center]

Thus to retrieve the data corresponding to an ID of center we evaluate

points[[ pointindex[center] ]]

which is equivalent to our previous method where we had defined points to be a function. What we have gained by defining points to be a list is automatic updating of pointlist and no requirement for pointers to pointlist in the points database. We can now apply Mathematica's built-in list manipulation tools directly to the database and for small models this improves speed in all areas, particularly deletes. When we get to larger databases and routines that need to randomly access the database and thus reference pointlist and pointindex a lot, the performance starts to suffer because of the time required to do the Transpose in pointlist and Position in pointindex. A hybrid approach is to use lists for the points database, and to treat pointlist and pointindex as assignments rather than delayed assignments,  only refreshing their definitions when the database is modified. This is the database design adopted throughout the latest iteration of the ModelMaker package, because it offered a reasonable compromise between time required for random access and time required for updates of the data.

Moving on to other points database manipulations, TRANSLATE, ROTATE and MIRROR  are easy to implement with functional programming, but EQUIVALENCE was more interesting. Press et al [5] presented a Fortran subroutine called ECLAZZ which was coded almost literally in Mathematica.  It runs slowly because it is written with Do loops. Lee [6] published a recursive equivalencer written using functional programming and it runs much faster than a procedural coding of ECLAZZ. The nonrecursive equivalencer in ModelMaker is faster yet, and is based on assigning points to voxels (volume equivalent of pixels) so the comparison on distance between points is done using integers.

Here are the point functions that have been coded thus far in ModelMaker.
The symbol id is a single identifier, ids means a list of identifiers.

To use these definitions first load the ModelMaker package with:

<<mm`

point[] ;                  (* new point at the origin, next avbl id *)
point[NEW,{x,y,z}];        (* new point at x,y,z, next avbl id*)
point[NEW,{x,y,z},id];     (* new point at x,y,z, with id *)
point[PLOT];               (* plot all points with no options *)
point[PLOT,,options];      (* plot all points with Graphics3D options*)
point[PLOT,id,options] ;   (* plot one point *)
point[PLOT,ids,options];   (* plot a list of points *)
point[PLOT,{{id1,id2,id3:1}},options]; (* plot range of points step id3 *)
point[DELETE];             (* delete all points *)
point[DELETE,ids];         (* delete list of points *)
point[DELETE,id,option] ;  (* delete one point *)
point[DELETE,{{id1,id2,id3:1}}];     (* delete range of points *)
point[TRANSLATE,ids,{dx,dy,dz},n];   (* translate list of points n times *)
point[ROTATE,ids,{roll,pitch,yaw},{x0,y0,z0},n]; (* create n rotated points,
 Euler angles {roll, pitch, yaw} about center of rotation {x0,y0,z0} *)
point[MIRROR,plane,ids];   (* mirror points about plane="X", "Y" or "Z" *)
point[EQUIVALENCE];        (* equivalence all the points *)
point[EQUIVALENCE,ids];    (* equivalence a list of points *)
point[RENUMBER];           (* renumber all points *)
point[CURVELIST,ids];      (* list the points belonging to curves *)
point[INTERSECT,id1,id2];  (* new point at intersection of curves id1,id2*)
point[NEXT];               (* the next point id to be created *)
point[LIST,id1,id2,id3:1]; (* returns a list of point ids *)
point[UNUSED];             (* returns points not used by curves *)

2.2  Curves

Curves in ModelMaker are vector-valued functions of a parameter, say s, which is bounded in the interval 0<s<1. The curves are described entirely by a list of point IDs, and a user-specified function named by the curve option InterpolationFunction which describes how to interpolate between the points. Note if the default interpolation function (polynomial) is used, only 2, 3 or 4 points can be in the list of IDs (corresponding to linear, quadratic, or cubic space curves). For checking that the list of points is the correct length for a user defined interpolation function, user defined interpolation functions should set the option datasize.  For example,  if the function is defined by:

fun[x_,data_List]  :=  ...

then use:

Options[fun]  =  {datasize -> 3};

to require that the function fun needs 3 points to define it. When the curve[NEW, ...] function checks user input
for errors, it evaluates

datasize /. Options[fun] 

to determine if the required number of datapoints have been provided.

To introduce curves, start by creating the default curve, a unit length straight line segment along the X-axis, plot it,  and display the curves database.

geometry[DELETE]; 
curve[];
curve[PLOT];
curves

Just as for points, the ID is the first entry in each row of the array curves, followed by a field called owner whose purpose is to list any surfaces that use this curve (here owner is an empty list). The next entry is the head of the function that interpolates between the points listed in the last entry in the record.

We plot curves using the PLOT tasking. Note that labels and points are included, but we can switch these and other graphics features on and off with a combination of global variables, and Graphics3D options. The use of global variables is not good Mathematica  practise, and they will probably be phased out in future versions of the package.

If we now examine the points database:

points

we see 2 points both owned by curve 1 have been created.  When curves get deleted by the curve[DELETE, ...] function the corresponding point owner records will be updated, but by setting the curve option childdelete->True, the points will also be deleted if they are solely owned by the curve. The owner record is essential for keeping a consistent database since otherwise deleting points 1 or 2 would cause curve 1 to become undefined. Deleting points 1 or 2 will also cause curve 1 to be deleted if the point option parentdelete->True (the default).  Keeping track of owners slows down database operations noticeably, so there is a price to pay for protecting the innocent.

We create curves using the NEW task, as in the following example, where we have set a global variable plotpoints, and used the Graphics3D option PlotRange.

point[DELETE]; curve[DELETE]; plotpoints=15;
point[]; 
point[NEW,{1,0,0}];
point[NEW,{1,1,0}];
point[NEW,{1,1,1}];
curve[NEW,{1,2,3,4}]; 
curve[PLOT,1,PlotRange->All];

The default interpolation function is a parametric polynomial space curve, of order up to 3. The equation for a cubic parametric space curve was derived using Mathematica's routine LinearSolve to solve for the coefficents of a cubic through 4 arbitrary points in space, under the assumption that the points correspond to the parametric coordinates s=0,1/3,2/3,1. This choice is arbitrary and in fact has only a subtle effect on the shape of the curve between the points. To get the coefficents of the cubic the following procedure was used.

(* Solve for the equation of a cubic polynomial space curve x[s],
through any 4 points p[1]...p[4] in 3-space. *)

Clear[a,A,x,"p*",M,m,b,sol,s];
a:=Array[A,{4}];
(* This creates our description of
the cubic space curve x[s_], now must solve for
the coefficients A[i] *)

x[s_]=Sum[(A[i] * s^(i-1) ),{i,1,4}]; 

<<LinearAlgebra`Orthogonalization` (* get Identity matrix*)

Index[i_Integer]:=i ; (* trivial helper function 
to flatten array "a", in higher dimensions
this gets more interesting *)

(* Construct orthonormal basis *)
Do[(A[i]=IdentityMatrix[4][[Index[i]]]),{i,1,4}]; 

p[1]= x[0];  (* now define the points as symbols *)
p[2]= x[1/3];(* p[i_Integer] are lists of 3 coordinates*)
p[3]= x[2/3];
p[4]= x[1];

M[row_Integer]=Block[{i},
(i[n_Integer]:=n ;
Table[p[row].A[i[j]],{j,1,4}])]; 

m=Table[M[i],{i,1,4}]; 
b={p1,p2,p3,p4};
sol=LinearSolve[m,b]; (* Solve m.sol=b for the 
vector "sol" of 4 elements *)
Do[(A[i]=sol[[Index[i]]]),{i,1,4}]; 

Clear[x];
x[s_]=Sum[(A[i] * s^(i-1) ),{i,1,4}];
Simplify[x[s]] (* The result, right hand side of x[s] *)

The expressions for the quadratic and linear space curves can be derived by the above procedure by reducing the number of points, but their expressions are so simple they can be solved by inspection. The utility of the above method becomes clear when cubic surfaces and solids are derived, and the algebra would be very difficult to do by hand. These are discussed in Sections 2.3 and 2.4.

The cubic, quadratic and linear space curves are fully specified by the coordinates of the points through which they must pass. There are no free parameters available to control the slopes of the curves at the endpoints. Thus there may arise situations where 2 space curves of these types joined end to end will not result in an adequately fair curve through the join. Interpolation schemes that smoothly handle such cases use polynomials with additional control parameters, or overlapping polynomials. Mathematica is ideal for experimenting with different interpolation schemes and its built-in Interpolation function can be used for this purpose. User-defined interpolation functions to replace the default cubic parametric space curves can be integrated with ModelMaker as the following example shows.

Consider a function  fun[s_,{id1_,id2_,id3_}] yet to be defined, that needs 3 points to define it. Using the points from the previous example, set the option InterpolationFunction->fun when  invoking the NEW  curve task, and then check the database for the result. Of course this can't be plotted until the function fun has been defined, but you can TRANSLATE, ROTATE, etc., since those operations only affect the points referenced by the curves.

curve[NEW,{1,2,3}, 5,InterpolationFunction->fun]

curves[[cindex[5] ]]

For a more concrete illustation make some points by creating a box using:

Plot[0,{x,0,1}];

Then select the resulting graphic and Command-click in the box with the mouse several times, and Copy and Paste the two dimensional coordinates. These will be mapped to 3D coordinates by the point[NEW...] function automatically.

(* The right hand side of data below were pasted in from the Clipboard. *)

data={{0.0562377, 0.182123}, {0.0889814, 0.473515}, 
{0.21177, 0.665568}, {0.383674, 0.639078}, 
{0.485998, 0.347687}, {0.735668, 0.274839}, 
{0.829806, 0.506627}, {0.915758, 0.698681}};

geometry[DELETE];
point[NEW,data];

Next create a curve using the ModelMaker package function Interpolation3D (which applies the built-in function Interpolation to each of the 3 vector components):

curve[NEW,pointlist,InterpolationFunction->Interpolation3D];
curve[PLOT];

Since curves are specified by lists of points, in order to apply transformations like TRANSLATE, ROTATE, MIRROR, etc. to the curves these tasks apply the transformations to the points, so most of the work is done by the point transformation routines. An example of translation follows:

curve[TRANSLATE,{1},{0,0,1},2];
curve[PLOT];

Here are all the curve functions developed to date.  Most curve plotting is done by the built-in routines ParametricPlot3D and Graphics3D and options to them are passed into the plot routines with global variables and by using the package FilterOptions.

curve[];     (* new unit length straight line along x axis *)
curve[NEW,ids,option];(* if InterpolationFunction->Automatic the default interpolation function is polynomial with Length[ids]==2,3 or 4 *)
curve[NEW,ids,id,option];   (* new curve with specified ID *)
curve[DELETE];              (* delete all curves, surfaces, solids *)
curve[DELETE,id, option];   (* delete a curve *)
curve[DELETE,{{id1,id2,id3:1}}]; (* delete a range of curves *)
curve[DELETE,ids];          (* delete a list of curves *)
curve[PLOT];                (* plot all curves, no options*)
curve[PLOT,,options];       (* plot all curves + Graphics3D options *)
curve[PLOT,id,options] ;    (* plot a single curve *)
curve[PLOT,{{id1,id2,id3:1}},options];
curve[PLOT,ids,options];    (* plots random list of curves *)
curve[TRANSLATE,list,{dx,dy,dz},n];   (* create n translated curves *)
curve[ROTATE,list,{roll,pitch,yaw},{x0,y0,z0},n]; (* create n rotated curves, Euler angles {roll, pitch, yaw} about center {x0,y0,z0} *)
curve[MIRROR,plane,ids];    (* mirror curves, plane="X", "Y" or "Z" *)
curve[NEXT];                (* the next curve id to be created *)
curve[LENGTH,id,s];         (* arc length  of curve id from 0 to s *)
curve[REVERSE];             (* reverse all curves *)
curve[REVERSE,id];          (* reverse curve id *)
curve[REVERSE,{{id1,id2,id3:1}}]; (* reverse a range of curves *)
curve[REVERSE,ids];         (* reverse a list of curves *)
curve[RAISE,id,numpoints];  (* interpolate to higher order polynomial *)
curve[FIT,order,ids];       (* nonlinear least squares fit of a
 polynomial space curve of order 1-3 through the list of points ids *)
curve[AVERAGE,id1,id2,weight:1/2]; (* average 2 curves *)
curve[AVERAGE,id1,id2,weight,id];  (* average curves with specified ID *)
curve[SURFACELIST,ids];     (* returns all curves in surface list ids *)
curve[RENUMBER];            (* renumber all curves *)
curve[LIST,id1,id2,id3:1];  (* returns a list of curve ids *)
curve[{x},id];  (* point on curve id at parametric coordinate x *)

2.3  Surfaces

Next we consider surfaces, and just as curves only need to reference the IDs of points, surfaces only reference the IDs of curves.  The data structure of each row in the list surfaces looks exactly like the structure of the list curves, except that the IDs refer not to points, but to curves. The geometrical interpretation of these curves is that they are cross-sections of the surface. The interpolation function that is used to interpolate between curves will default to cubic parametrics, but user-supplied functions can be input as the option InterpolationFunction. There are situations where a surface, when created, will create curves and hence own them, so the owner field in the curves database entries refers to ownership by surfaces. Similarly the owner field in a surfaces database entry is a list of solids that own the surface.

Here are the surface tasks supported by ModelMaker. Most look very much like the corresponding curve tasks, and in their implementation they often just call the curve routines. For example, the routine to ROTATE a surface just calls the curve[ROTATE, ...] routine. This reusability of code means a time savings in writing, and in maintenance.

surface[];                   (* create a unit square surface *)
surface[NEW,ids,id, option]; (* new surface id made from list of curves ids *)
surface[DELETE];
surface[DELETE,id, options];
surface[DELETE,{{id1,id2,id3:1}}];
surface[PLOT];               
surface[PLOT,id,option]; 
surface[PLOT,ids,option];
surface[PLOT,{{id1,id2,id3}},option];
surface[AREA,id,s:1,t:1];  (* the area of a surface *)
surface[TRANSLATE,ids,{dx,dy,dz},n]; 
surface[ROTATE,ids,{roll,pitch,yaw},{x0,y0,z0},n];
surface[NEXT];           (* the next surface id to be created *)
surface[REVERSE];        (* reverse all surfaces *)
surface[REVERSE,id];     (* reverse surface id, replaces original *)
surface[REVERSE,{{id1,id2,id3:1}}];
surface[REVERSE,ids];
surface[MIRROR,plane,ids];
surface[RAISE,id,n];(* interpolate a polynomial surface to higher order *)
surface[NORMAL,id]; (* normal to a surface at parametric coordinates {0,0}*)
surface[AVERAGE,id1,id2,weight:1/2];  (* average 2 surfaces *)
surface[AVERAGE,id1,id2,weight:1/2,id];
surface[RENUMBER];            
surface[CENTROID,id];          (* finds the centroid of a surface *)
surface[LIST,id1,id2,id3:1]    (* returns list of surface ids *)
surface[PLOTCURVES,ids,option];(* plot curves that define surfaces ids *)
surface[{x,y},id];  (* coordinates of point on surface id 
at parametric coordinate {x,y} *)

The default interpolation function for surfaces is the cubic polynomial. Deriving the formulae for it by hand is sufficently difficult as to be impractical and Mathematica was used to generate the explicit form using the routine LinearSolve, so that it could be compiled for speed. The uncompiled version can be used when infinite precision evaluations are needed, and this is selected with the surface and solid option Compiled. Here we derive the expression for the cubic surface. Now there are 2 parameters, say { s, t } which we bound in the unit square. The unknown coefficents of the polynomials are found under the assumption that the points correspond to the parametric coordinates s = { 0, 1/3, 2/3, 1 } and t = { 0, 1/3, 2/3, 1 }. The quadratic and linear surfaces can be handled similarly,  by reducing the number of points.

(* Solve for the equation of a cubic polynomial space surface x[u,v], through 16 points p[1]...p[16] in 3-space assumed to lie on a set of 4 cross-sectional curves such that p[1]..p[4] defines the first curve, p[5]...p[8] the second and so on. *)

Clear[a,A,x,xs,"p*",M,m,b,sol,Index,i,j,u,v];
a:=Array[A,{4,4}];

x[u_,v_]=Sum[(A[i,j] * u^(i-1) * v^(j-1) ),
{i,1,4},{j,1,4}]; 
(* The cubic surface, now solve for coefficients A[i,j] *)

<<LinearAlgebra`Orthogonalization` (* get Identity matrix *)

Index[i_,j_]:=i +4(j-1) ; (* To flatten array "a" *)

Do[Do[(A[i,j]=
IdentityMatrix[16][[Index[i,j]]]),
{i,1,4}],{j,1,4}]; (* construct an orthonormal basis *)
(* now define the points as g[1]= x[0,    0]; 
                            g[2]= x[1/3,  0]; ... etc...
                            with the following *)
xs=Map[Apply[x,#]&,Flatten[Table[{u,v},
{v,0,1,1/3},{u,0,1,1/3}],1]];
Do[p[i]=xs[[i]],{i,1,16}];

M[row_Integer]=Block[{i,j},
(j[n_Integer]:=1+Quotient[n-1,4];
i[n_Integer]:=n - 4 (j[n]-1);
Table[p[row].A[i[n],j[n]],{n,1,16}])]; 

m=Table[M[i],{i,1,16}]; (* the complete m matrix *)
b={p1,p2,p3,p4,p5,p6,p7,p8,
p9,p10,p11,p12,p13,p14,p15,p16};
sol=LinearSolve[m,b]; (* Solve m.sol=b for the 
vector solution sol of 16 elements *)

Do[Do[(A[i,j]=sol[[Index[i,j]]]),
{i,1,4}],{j,1,4}]; 

Clear[x];
x[u_,v_]=Sum[(A[i,j] * u^(i-1) * v^(j-1)),{i,1,4},{j,1,4}];

(* Print results and simplify*)

j[n_Integer]:=1+Quotient[n-1,4];
i[n_Integer]:=n - 4 (j[n]-1);
Print["("];
Do[Print["(",InputForm[Simplify[A[i[n],j[n]]* u^(i[n]-1) *
v^(j[n]-1) ]],If[n==16,"))",")+"]],{n,1,16}]

(* The right hand side of cubic polynomial x[u,v] 
where the vectors p1, p2...p16 are the points 
on the surface. *)

2.4  Solids

Next we consider solids, and just as surfaces only needed to reference the IDs of curves, solids only reference the IDs of surfaces. The geometrical interpretation of these surfaces is that they are cross-sections of the solid. The function that is used to interpolate between surfaces will default to cubic parametrics, but user-supplied functions can be input as the option InterpolationFunction to the solid[NEW, ...] function.  The owner field in a solid database entry is presently not used, but is reserved for future addition of another class of entity,  to be called part or group, which may be useful to represent major subassemblies of a large model.

Here are the solid tasks that have thus far been coded. Most look very much like the corresponding surface tasks, and in their implementation they often just call the surface routines.

solid[];
solid[NEW,ids];
solid[NEW,ids,id, option];
solid[DELETE];
solid[DELETE,id, options];
solid[DELETE,{{id1,id2,id3:1}}];
solid[DELETE,ids];
solid[PLOT];
solid[PLOT,,option];
solid[PLOT,id,option];
solid[PLOT,{{id1,id2,id3}},option];
solid[PLOT,ids,option];
solid[TRANSLATE,ids,{dx,dy,dz},n]; 
solid[ROTATE,ids,{roll,pitch,yaw},{x0,y0,z0},n];
solid[MIRROR,plane,ids];
solid[RENUMBER];
solid[CENTROID,id]; (* location of centroid of a solid *)
solid[CHECK,id]; (* checks sign of volume *)
solid[NEXT];
solid[LIST,id1,id2,id3:1];
solid[{x,y,z},id];  (* coordinates of point in solid id 
at parametric coordinate {x,y,z} *)

The default interpolation function for solids is the cubic polynomial. Deriving the formula for it by hand is totally  impractical and Mathematica was used to generate the explicit form using the routine LinearSolve, so that it could be compiled for speed. The uncompiled version can be used when infinite precision evaluations are needed. Here we derive the expression for the cubic solid. Now there are 3 parameters, say { u,v,w } which are bounded in the unit cube. Our solution for the coefficents of the polynomials is found under the assumption that the 64 points required to specify the cubic solid correspond to all combinations of the parametric coordinates u = { 0, 1/3, 2/3, 1} and v = { 0, 1/3, 2/3, 1} and
 w = { 0, 1/3, 2/3, 1}. The quadratic and linear surfaces can be handled by reducing the number of points.

(* Solve for the equation of a cubic polynomial space solid x[u,v,w], through 64 points p[1]...p[64] in 3-space assumed to lie on a set of 4 cross-sectional surfaces such that p[1]..p[16] defines the first surface, p[17]...p[32] the second and so on. *)
Clear[a,A,x,xs,"p*",M,m,b,sol,Index,u,v,w];

a:=Array[A,{4,4,4}];

(* The cubic solid, now find unknown coefficients A[i,j,k]... *)
x[u_,v_,w_]=Sum[(A[i,j,k] * u^(i-1) * v^(j-1) * w^(k-1)),
{i,1,4},{j,1,4},{k,1,4}]; 

<<LinearAlgebra`Orthogonalization` (* Identity matrix *)

Index[i_,j_,k_]:=i +4(j-1) + 16(k-1); (* flattens array "a" *)

Do[Do[Do[(A[i,j,k]=
IdentityMatrix[64][[Index[i,j,k]]]),
{i,1,4}],{j,1,4}],{k,1,4}] (* construct orthonormal basis *)

(* now define the points as
follows: g[1]= x[0,    0,  0];
         g[2]= x[1/3,  0,  0];... etc. by the following: *)

Clear[g,xs];
xs=Map[Apply[x,#]&,Flatten[Table[{u,v,w},
{w,0,1,1/3},{v,0,1,1/3},{u,0,1,1/3}],2]];
Do[p[i]=xs[[i]],{i,1,64}];

M[row_Integer]=Block[{i,j,k},
(k[n_Integer]:=1+Quotient[n-1,16];
j[n_Integer]:=1+Quotient[(n-16 (k[n]-1))-1,4];
i[n_Integer]:=n-16 (k[n]-1) - 4 (j[n]-1);
Table[p[row].A[i[n],j[n],k[n]],{n,1,64}])]; 

m=Table[M[i],{i,1,64}]; 
b={p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12,p13,
  p14,p15,p16,p17,p18,p19,p20,p21,p22,p23,p24,p25,
  p26,p27,p28,p29,p30,p31,p32,p33,p34,p35,p36,p37,
  p38,p39,p40,p41,p42,p43,p44,p45,p46,p47,p48,p49,
  p50,p51,p52,p53,p54,p55,p56,p57,p58,p59,p60,p61,
  p62,p63,p64};
  
sol = LinearSolve[m,b]; 
 
Do[Do[Do[(A[i,j,k]=sol[[Index[i,j,k]]]), (* extract coefficents *)
{i,1,4}],{j,1,4}],{k,1,4}];  (* from the solution *)

Clear[X];
X[u_,v_,w_]=Sum[(A[i,j,k] * u^(i-1) * v^(j-1) * w^(k-1)),
{i,1,4},{j,1,4},{k,1,4}]; (* Assemble complete expression for
 cubic polynomial solid *)

Here is the expression for the cubic solid in short form (the full expression takes 20 pages to print out in full).

X[u,v,w]//Short

In the notebook version of this document you can look at the full expression for the cubic polynomial solid by expanding this closed group of cells.

X[u_,v_,w_]:=
(p1 + ((-11*p1)/2 + 9*p2 - (9*p3)/2 + p4)*u + 
 
   (9*p1 - (45*p2)/2 + 18*p3 - (9*p4)/2)*u^2 + 
 
   ((-9*p1)/2 + (27*p2)/2 - (27*p3)/2 + (9*p4)/2)*u^3 + 
 
   ((-11*p1)/2 + p13 + 9*p5 - (9*p9)/2)*v + 
 
   ((121*p1)/4 - (81*p10)/2 + (81*p11)/4 - (9*p12)/2 - 
 
      (11*p13)/2 + 9*p14 - (9*p15)/2 + p16 - (99*p2)/2 + 
 
      (99*p3)/4 - (11*p4)/2 - (99*p5)/2 + 81*p6 - 
 
      (81*p7)/2 + 9*p8 + (99*p9)/4)*u*v + 
 
   ((-99*p1)/2 + (405*p10)/4 - 81*p11 + (81*p12)/4 + 
 
      9*p13 - (45*p14)/2 + 18*p15 - (9*p16)/2 + 
 
      (495*p2)/4 - 99*p3 + (99*p4)/4 + 81*p5 - 
 
      (405*p6)/2 + 162*p7 - (81*p8)/2 - (81*p9)/2)*u^2*v + 
 
   ((99*p1)/4 - (243*p10)/4 + (243*p11)/4 - (81*p12)/4 - 
 
      (9*p13)/2 + (27*p14)/2 - (27*p15)/2 + (9*p16)/2 - 
 
      (297*p2)/4 + (297*p3)/4 - (99*p4)/4 - (81*p5)/2 + 
 
      (243*p6)/2 - (243*p7)/2 + (81*p8)/2 + (81*p9)/4)*u^3*
 
    v + (9*p1 - (9*p13)/2 - (45*p5)/2 + 18*p9)*v^2 + 
 
   ((-99*p1)/2 + 162*p10 - 81*p11 + 18*p12 + (99*p13)/4 - 
 
      (81*p14)/2 + (81*p15)/4 - (9*p16)/2 + 81*p2 - 
 
      (81*p3)/2 + 9*p4 + (495*p5)/4 - (405*p6)/2 + 
 
      (405*p7)/4 - (45*p8)/2 - 99*p9)*u*v^2 + 
 
   (81*p1 - 405*p10 + 324*p11 - 81*p12 - (81*p13)/2 + 
 
      (405*p14)/4 - 81*p15 + (81*p16)/4 - (405*p2)/2 + 
 
      162*p3 - (81*p4)/2 - (405*p5)/2 + (2025*p6)/4 - 
 
      405*p7 + (405*p8)/4 + 162*p9)*u^2*v^2 + 
 
   ((-81*p1)/2 + 243*p10 - 243*p11 + 81*p12 + (81*p13)/4 - 
 
      (243*p14)/4 + (243*p15)/4 - (81*p16)/4 + 
 
      (243*p2)/2 - (243*p3)/2 + (81*p4)/2 + (405*p5)/4 - 
 
      (1215*p6)/4 + (1215*p7)/4 - (405*p8)/4 - 81*p9)*u^3*
 
    v^2 + ((-9*p1)/2 + (9*p13)/2 + (27*p5)/2 - (27*p9)/2)*
 
    v^3 + ((99*p1)/4 - (243*p10)/2 + (243*p11)/4 - 
 
      (27*p12)/2 - (99*p13)/4 + (81*p14)/2 - (81*p15)/4 + 
 
      (9*p16)/2 - (81*p2)/2 + (81*p3)/4 - (9*p4)/2 - 
 
      (297*p5)/4 + (243*p6)/2 - (243*p7)/4 + (27*p8)/2 + 
 
      (297*p9)/4)*u*v^3 + 
 
   ((-81*p1)/2 + (1215*p10)/4 - 243*p11 + (243*p12)/4 + 
 
      (81*p13)/2 - (405*p14)/4 + 81*p15 - (81*p16)/4 + 
 
      (405*p2)/4 - 81*p3 + (81*p4)/4 + (243*p5)/2 - 
 
      (1215*p6)/4 + 243*p7 - (243*p8)/4 - (243*p9)/2)*u^2*
 
    v^3 + ((81*p1)/4 - (729*p10)/4 + (729*p11)/4 - 
 
      (243*p12)/4 - (81*p13)/4 + (243*p14)/4 - 
 
      (243*p15)/4 + (81*p16)/4 - (243*p2)/4 + (243*p3)/4 - 
 
      (81*p4)/4 - (243*p5)/4 + (729*p6)/4 - (729*p7)/4 + 
 
      (243*p8)/4 + (243*p9)/4)*u^3*v^3 + 
 
   ((-11*p1)/2 + 9*p17 - (9*p33)/2 + p49)*w + 
 
   ((121*p1)/4 - (99*p17)/2 + 81*p18 - (81*p19)/2 - 
 
      (99*p2)/2 + 9*p20 + (99*p3)/4 + (99*p33)/4 - 
 
      (81*p34)/2 + (81*p35)/4 - (9*p36)/2 - (11*p4)/2 - 
 
      (11*p49)/2 + 9*p50 - (9*p51)/2 + p52)*u*w + 
 
   ((-99*p1)/2 + 81*p17 - (405*p18)/2 + 162*p19 + 
 
      (495*p2)/4 - (81*p20)/2 - 99*p3 - (81*p33)/2 + 
 
      (405*p34)/4 - 81*p35 + (81*p36)/4 + (99*p4)/4 + 
 
      9*p49 - (45*p50)/2 + 18*p51 - (9*p52)/2)*u^2*w + 
 
   ((99*p1)/4 - (81*p17)/2 + (243*p18)/2 - (243*p19)/2 - 
 
      (297*p2)/4 + (81*p20)/2 + (297*p3)/4 + (81*p33)/4 - 
 
      (243*p34)/4 + (243*p35)/4 - (81*p36)/4 - (99*p4)/4 - 
 
      (9*p49)/2 + (27*p50)/2 - (27*p51)/2 + (9*p52)/2)*u^3*
 
    w + ((121*p1)/4 - (11*p13)/2 - (99*p17)/2 + 81*p21 - 
 
      (81*p25)/2 + 9*p29 + (99*p33)/4 - (81*p37)/2 + 
 
      (81*p41)/4 - (9*p45)/2 - (11*p49)/2 - (99*p5)/2 + 
 
      9*p53 - (9*p57)/2 + p61 + (99*p9)/4)*v*w + 
 
   ((-1331*p1)/8 + (891*p10)/4 - (891*p11)/8 + 
 
      (99*p12)/4 + (121*p13)/4 - (99*p14)/2 + (99*p15)/4 - 
 
      (11*p16)/2 + (1089*p17)/4 - (891*p18)/2 + 
 
      (891*p19)/4 + (1089*p2)/4 - (99*p20)/2 - 
 
      (891*p21)/2 + 729*p22 - (729*p23)/2 + 81*p24 + 
 
      (891*p25)/4 - (729*p26)/2 + (729*p27)/4 - 
 
      (81*p28)/2 - (99*p29)/2 - (1089*p3)/8 + 81*p30 - 
 
      (81*p31)/2 + 9*p32 - (1089*p33)/8 + (891*p34)/4 - 
 
      (891*p35)/8 + (99*p36)/4 + (891*p37)/4 - 
 
      (729*p38)/2 + (729*p39)/4 + (121*p4)/4 - 
 
      (81*p40)/2 - (891*p41)/8 + (729*p42)/4 - 
 
      (729*p43)/8 + (81*p44)/4 + (99*p45)/4 - (81*p46)/2 + 
 
      (81*p47)/4 - (9*p48)/2 + (121*p49)/4 + (1089*p5)/4 - 
 
      (99*p50)/2 + (99*p51)/4 - (11*p52)/2 - (99*p53)/2 + 
 
      81*p54 - (81*p55)/2 + 9*p56 + (99*p57)/4 - 
 
      (81*p58)/2 + (81*p59)/4 - (891*p6)/2 - (9*p60)/2 - 
 
      (11*p61)/2 + 9*p62 - (9*p63)/2 + p64 + (891*p7)/4 - 
 
      (99*p8)/2 - (1089*p9)/8)*u*v*w + 
 
   ((1089*p1)/4 - (4455*p10)/8 + (891*p11)/2 - 
 
      (891*p12)/8 - (99*p13)/2 + (495*p14)/4 - 99*p15 + 
 
      (99*p16)/4 - (891*p17)/2 + (4455*p18)/4 - 891*p19 - 
 
      (5445*p2)/8 + (891*p20)/4 + 729*p21 - (3645*p22)/2 + 
 
      1458*p23 - (729*p24)/2 - (729*p25)/2 + 
 
      (3645*p26)/4 - 729*p27 + (729*p28)/4 + 81*p29 + 
 
      (1089*p3)/2 - (405*p30)/2 + 162*p31 - (81*p32)/2 + 
 
      (891*p33)/4 - (4455*p34)/8 + (891*p35)/2 - 
 
      (891*p36)/8 - (729*p37)/2 + (3645*p38)/4 - 729*p39 - 
 
      (1089*p4)/8 + (729*p40)/4 + (729*p41)/4 - 
 
      (3645*p42)/8 + (729*p43)/2 - (729*p44)/8 - 
 
      (81*p45)/2 + (405*p46)/4 - 81*p47 + (81*p48)/4 - 
 
      (99*p49)/2 - (891*p5)/2 + (495*p50)/4 - 99*p51 + 
 
      (99*p52)/4 + 81*p53 - (405*p54)/2 + 162*p55 - 
 
      (81*p56)/2 - (81*p57)/2 + (405*p58)/4 - 81*p59 + 
 
      (4455*p6)/4 + (81*p60)/4 + 9*p61 - (45*p62)/2 + 
 
      18*p63 - (9*p64)/2 - 891*p7 + (891*p8)/4 + (891*p9)/4
 
      )*u^2*v*w + ((-1089*p1)/8 + (2673*p10)/8 - 
 
      (2673*p11)/8 + (891*p12)/8 + (99*p13)/4 - 
 
      (297*p14)/4 + (297*p15)/4 - (99*p16)/4 + 
 
      (891*p17)/4 - (2673*p18)/4 + (2673*p19)/4 + 
 
      (3267*p2)/8 - (891*p20)/4 - (729*p21)/2 + 
 
      (2187*p22)/2 - (2187*p23)/2 + (729*p24)/2 + 
 
      (729*p25)/4 - (2187*p26)/4 + (2187*p27)/4 - 
 
      (729*p28)/4 - (81*p29)/2 - (3267*p3)/8 + 
 
      (243*p30)/2 - (243*p31)/2 + (81*p32)/2 - 
 
      (891*p33)/8 + (2673*p34)/8 - (2673*p35)/8 + 
 
      (891*p36)/8 + (729*p37)/4 - (2187*p38)/4 + 
 
      (2187*p39)/4 + (1089*p4)/8 - (729*p40)/4 - 
 
      (729*p41)/8 + (2187*p42)/8 - (2187*p43)/8 + 
 
      (729*p44)/8 + (81*p45)/4 - (243*p46)/4 + 
 
      (243*p47)/4 - (81*p48)/4 + (99*p49)/4 + (891*p5)/4 - 
 
      (297*p50)/4 + (297*p51)/4 - (99*p52)/4 - 
 
      (81*p53)/2 + (243*p54)/2 - (243*p55)/2 + 
 
      (81*p56)/2 + (81*p57)/4 - (243*p58)/4 + 
 
      (243*p59)/4 - (2673*p6)/4 - (81*p60)/4 - (9*p61)/2 + 
 
      (27*p62)/2 - (27*p63)/2 + (9*p64)/2 + (2673*p7)/4 - 
 
      (891*p8)/4 - (891*p9)/8)*u^3*v*w + 
 
   ((-99*p1)/2 + (99*p13)/4 + 81*p17 - (405*p21)/2 + 
 
      162*p25 - (81*p29)/2 - (81*p33)/2 + (405*p37)/4 - 
 
      81*p41 + (81*p45)/4 + 9*p49 + (495*p5)/4 - 
 
      (45*p53)/2 + 18*p57 - (9*p61)/2 - 99*p9)*v^2*w + 
 
   ((1089*p1)/4 - 891*p10 + (891*p11)/2 - 99*p12 - 
 
      (1089*p13)/8 + (891*p14)/4 - (891*p15)/8 + 
 
      (99*p16)/4 - (891*p17)/2 + 729*p18 - (729*p19)/2 - 
 
      (891*p2)/2 + 81*p20 + (4455*p21)/4 - (3645*p22)/2 + 
 
      (3645*p23)/4 - (405*p24)/2 - 891*p25 + 1458*p26 - 
 
      729*p27 + 162*p28 + (891*p29)/4 + (891*p3)/4 - 
 
      (729*p30)/2 + (729*p31)/4 - (81*p32)/2 + 
 
      (891*p33)/4 - (729*p34)/2 + (729*p35)/4 - 
 
      (81*p36)/2 - (4455*p37)/8 + (3645*p38)/4 - 
 
      (3645*p39)/8 - (99*p4)/2 + (405*p40)/4 + 
 
      (891*p41)/2 - 729*p42 + (729*p43)/2 - 81*p44 - 
 
      (891*p45)/8 + (729*p46)/4 - (729*p47)/8 + 
 
      (81*p48)/4 - (99*p49)/2 - (5445*p5)/8 + 81*p50 - 
 
      (81*p51)/2 + 9*p52 + (495*p53)/4 - (405*p54)/2 + 
 
      (405*p55)/4 - (45*p56)/2 - 99*p57 + 162*p58 - 
 
      81*p59 + (4455*p6)/4 + 18*p60 + (99*p61)/4 - 
 
      (81*p62)/2 + (81*p63)/4 - (9*p64)/2 - (4455*p7)/8 + 
 
      (495*p8)/4 + (1089*p9)/2)*u*v^2*w + 
 
   ((-891*p1)/2 + (4455*p10)/2 - 1782*p11 + (891*p12)/2 + 
 
      (891*p13)/4 - (4455*p14)/8 + (891*p15)/2 - 
 
      (891*p16)/8 + 729*p17 - (3645*p18)/2 + 1458*p19 + 
 
      (4455*p2)/4 - (729*p20)/2 - (3645*p21)/2 + 
 
      (18225*p22)/4 - 3645*p23 + (3645*p24)/4 + 1458*p25 - 
 
      3645*p26 + 2916*p27 - 729*p28 - (729*p29)/2 - 
 
      891*p3 + (3645*p30)/4 - 729*p31 + (729*p32)/4 - 
 
      (729*p33)/2 + (3645*p34)/4 - 729*p35 + (729*p36)/4 + 
 
      (3645*p37)/4 - (18225*p38)/8 + (3645*p39)/2 + 
 
      (891*p4)/4 - (3645*p40)/8 - 729*p41 + (3645*p42)/2 - 
 
      1458*p43 + (729*p44)/2 + (729*p45)/4 - 
 
      (3645*p46)/8 + (729*p47)/2 - (729*p48)/8 + 81*p49 + 
 
      (4455*p5)/4 - (405*p50)/2 + 162*p51 - (81*p52)/2 - 
 
      (405*p53)/2 + (2025*p54)/4 - 405*p55 + (405*p56)/4 + 
 
      162*p57 - 405*p58 + 324*p59 - (22275*p6)/8 - 
 
      81*p60 - (81*p61)/2 + (405*p62)/4 - 81*p63 + 
 
      (81*p64)/4 + (4455*p7)/2 - (4455*p8)/8 - 891*p9)*u^2*
 
    v^2*w + ((891*p1)/4 - (2673*p10)/2 + (2673*p11)/2 - 
 
      (891*p12)/2 - (891*p13)/8 + (2673*p14)/8 - 
 
      (2673*p15)/8 + (891*p16)/8 - (729*p17)/2 + 
 
      (2187*p18)/2 - (2187*p19)/2 - (2673*p2)/4 + 
 
      (729*p20)/2 + (3645*p21)/4 - (10935*p22)/4 + 
 
      (10935*p23)/4 - (3645*p24)/4 - 729*p25 + 2187*p26 - 
 
      2187*p27 + 729*p28 + (729*p29)/4 + (2673*p3)/4 - 
 
      (2187*p30)/4 + (2187*p31)/4 - (729*p32)/4 + 
 
      (729*p33)/4 - (2187*p34)/4 + (2187*p35)/4 - 
 
      (729*p36)/4 - (3645*p37)/8 + (10935*p38)/8 - 
 
      (10935*p39)/8 - (891*p4)/4 + (3645*p40)/8 + 
 
      (729*p41)/2 - (2187*p42)/2 + (2187*p43)/2 - 
 
      (729*p44)/2 - (729*p45)/8 + (2187*p46)/8 - 
 
      (2187*p47)/8 + (729*p48)/8 - (81*p49)/2 - 
 
      (4455*p5)/8 + (243*p50)/2 - (243*p51)/2 + 
 
      (81*p52)/2 + (405*p53)/4 - (1215*p54)/4 + 
 
      (1215*p55)/4 - (405*p56)/4 - 81*p57 + 243*p58 - 
 
      243*p59 + (13365*p6)/8 + 81*p60 + (81*p61)/4 - 
 
      (243*p62)/4 + (243*p63)/4 - (81*p64)/4 - 
 
      (13365*p7)/8 + (4455*p8)/8 + (891*p9)/2)*u^3*v^2*w + 
 
   ((99*p1)/4 - (99*p13)/4 - (81*p17)/2 + (243*p21)/2 - 
 
      (243*p25)/2 + (81*p29)/2 + (81*p33)/4 - 
 
      (243*p37)/4 + (243*p41)/4 - (81*p45)/4 - (9*p49)/2 - 
 
      (297*p5)/4 + (27*p53)/2 - (27*p57)/2 + (9*p61)/2 + 
 
      (297*p9)/4)*v^3*w + 
 
   ((-1089*p1)/8 + (2673*p10)/4 - (2673*p11)/8 + 
 
      (297*p12)/4 + (1089*p13)/8 - (891*p14)/4 + 
 
      (891*p15)/8 - (99*p16)/4 + (891*p17)/4 - 
 
      (729*p18)/2 + (729*p19)/4 + (891*p2)/4 - 
 
      (81*p20)/2 - (2673*p21)/4 + (2187*p22)/2 - 
 
      (2187*p23)/4 + (243*p24)/2 + (2673*p25)/4 - 
 
      (2187*p26)/2 + (2187*p27)/4 - (243*p28)/2 - 
 
      (891*p29)/4 - (891*p3)/8 + (729*p30)/2 - 
 
      (729*p31)/4 + (81*p32)/2 - (891*p33)/8 + 
 
      (729*p34)/4 - (729*p35)/8 + (81*p36)/4 + 
 
      (2673*p37)/8 - (2187*p38)/4 + (2187*p39)/8 + 
 
      (99*p4)/4 - (243*p40)/4 - (2673*p41)/8 + 
 
      (2187*p42)/4 - (2187*p43)/8 + (243*p44)/4 + 
 
      (891*p45)/8 - (729*p46)/4 + (729*p47)/8 - 
 
      (81*p48)/4 + (99*p49)/4 + (3267*p5)/8 - (81*p50)/2 + 
 
      (81*p51)/4 - (9*p52)/2 - (297*p53)/4 + (243*p54)/2 - 
 
      (243*p55)/4 + (27*p56)/2 + (297*p57)/4 - 
 
      (243*p58)/2 + (243*p59)/4 - (2673*p6)/4 - 
 
      (27*p60)/2 - (99*p61)/4 + (81*p62)/2 - (81*p63)/4 + 
 
      (9*p64)/2 + (2673*p7)/8 - (297*p8)/4 - (3267*p9)/8)*
 
    u*v^3*w + ((891*p1)/4 - (13365*p10)/8 + (2673*p11)/2 - 
 
      (2673*p12)/8 - (891*p13)/4 + (4455*p14)/8 - 
 
      (891*p15)/2 + (891*p16)/8 - (729*p17)/2 + 
 
      (3645*p18)/4 - 729*p19 - (4455*p2)/8 + (729*p20)/4 + 
 
      (2187*p21)/2 - (10935*p22)/4 + 2187*p23 - 
 
      (2187*p24)/4 - (2187*p25)/2 + (10935*p26)/4 - 
 
      2187*p27 + (2187*p28)/4 + (729*p29)/2 + (891*p3)/2 - 
 
      (3645*p30)/4 + 729*p31 - (729*p32)/4 + (729*p33)/4 - 
 
      (3645*p34)/8 + (729*p35)/2 - (729*p36)/8 - 
 
      (2187*p37)/4 + (10935*p38)/8 - (2187*p39)/2 - 
 
      (891*p4)/8 + (2187*p40)/8 + (2187*p41)/4 - 
 
      (10935*p42)/8 + (2187*p43)/2 - (2187*p44)/8 - 
 
      (729*p45)/4 + (3645*p46)/8 - (729*p47)/2 + 
 
      (729*p48)/8 - (81*p49)/2 - (2673*p5)/4 + 
 
      (405*p50)/4 - 81*p51 + (81*p52)/4 + (243*p53)/2 - 
 
      (1215*p54)/4 + 243*p55 - (243*p56)/4 - (243*p57)/2 + 
 
      (1215*p58)/4 - 243*p59 + (13365*p6)/8 + 
 
      (243*p60)/4 + (81*p61)/2 - (405*p62)/4 + 81*p63 - 
 
      (81*p64)/4 - (2673*p7)/2 + (2673*p8)/8 + (2673*p9)/4)
 
     *u^2*v^3*w + ((-891*p1)/8 + (8019*p10)/8 - 
 
      (8019*p11)/8 + (2673*p12)/8 + (891*p13)/8 - 
 
      (2673*p14)/8 + (2673*p15)/8 - (891*p16)/8 + 
 
      (729*p17)/4 - (2187*p18)/4 + (2187*p19)/4 + 
 
      (2673*p2)/8 - (729*p20)/4 - (2187*p21)/4 + 
 
      (6561*p22)/4 - (6561*p23)/4 + (2187*p24)/4 + 
 
      (2187*p25)/4 - (6561*p26)/4 + (6561*p27)/4 - 
 
      (2187*p28)/4 - (729*p29)/4 - (2673*p3)/8 + 
 
      (2187*p30)/4 - (2187*p31)/4 + (729*p32)/4 - 
 
      (729*p33)/8 + (2187*p34)/8 - (2187*p35)/8 + 
 
      (729*p36)/8 + (2187*p37)/8 - (6561*p38)/8 + 
 
      (6561*p39)/8 + (891*p4)/8 - (2187*p40)/8 - 
 
      (2187*p41)/8 + (6561*p42)/8 - (6561*p43)/8 + 
 
      (2187*p44)/8 + (729*p45)/8 - (2187*p46)/8 + 
 
      (2187*p47)/8 - (729*p48)/8 + (81*p49)/4 + 
 
      (2673*p5)/8 - (243*p50)/4 + (243*p51)/4 - 
 
      (81*p52)/4 - (243*p53)/4 + (729*p54)/4 - 
 
      (729*p55)/4 + (243*p56)/4 + (243*p57)/4 - 
 
      (729*p58)/4 + (729*p59)/4 - (8019*p6)/8 - 
 
      (243*p60)/4 - (81*p61)/4 + (243*p62)/4 - 
 
      (243*p63)/4 + (81*p64)/4 + (8019*p7)/8 - 
 
      (2673*p8)/8 - (2673*p9)/8)*u^3*v^3*w + 
 
   (9*p1 - (45*p17)/2 + 18*p33 - (9*p49)/2)*w^2 + 
 
   ((-99*p1)/2 + (495*p17)/4 - (405*p18)/2 + (405*p19)/4 + 
 
      81*p2 - (45*p20)/2 - (81*p3)/2 - 99*p33 + 162*p34 - 
 
      81*p35 + 18*p36 + 9*p4 + (99*p49)/4 - (81*p50)/2 + 
 
      (81*p51)/4 - (9*p52)/2)*u*w^2 + 
 
   (81*p1 - (405*p17)/2 + (2025*p18)/4 - 405*p19 - 
 
      (405*p2)/2 + (405*p20)/4 + 162*p3 + 162*p33 - 
 
      405*p34 + 324*p35 - 81*p36 - (81*p4)/2 - 
 
      (81*p49)/2 + (405*p50)/4 - 81*p51 + (81*p52)/4)*u^2*
 
    w^2 + ((-81*p1)/2 + (405*p17)/4 - (1215*p18)/4 + 
 
      (1215*p19)/4 + (243*p2)/2 - (405*p20)/4 - 
 
      (243*p3)/2 - 81*p33 + 243*p34 - 243*p35 + 81*p36 + 
 
      (81*p4)/2 + (81*p49)/4 - (243*p50)/4 + (243*p51)/4 - 
 
      (81*p52)/4)*u^3*w^2 + 
 
   ((-99*p1)/2 + 9*p13 + (495*p17)/4 - (405*p21)/2 + 
 
      (405*p25)/4 - (45*p29)/2 - 99*p33 + 162*p37 - 
 
      81*p41 + 18*p45 + (99*p49)/4 + 81*p5 - (81*p53)/2 + 
 
      (81*p57)/4 - (9*p61)/2 - (81*p9)/2)*v*w^2 + 
 
   ((1089*p1)/4 - (729*p10)/2 + (729*p11)/4 - (81*p12)/2 - 
 
      (99*p13)/2 + 81*p14 - (81*p15)/2 + 9*p16 - 
 
      (5445*p17)/8 + (4455*p18)/4 - (4455*p19)/8 - 
 
      (891*p2)/2 + (495*p20)/4 + (4455*p21)/4 - 
 
      (3645*p22)/2 + (3645*p23)/4 - (405*p24)/2 - 
 
      (4455*p25)/8 + (3645*p26)/4 - (3645*p27)/8 + 
 
      (405*p28)/4 + (495*p29)/4 + (891*p3)/4 - 
 
      (405*p30)/2 + (405*p31)/4 - (45*p32)/2 + 
 
      (1089*p33)/2 - 891*p34 + (891*p35)/2 - 99*p36 - 
 
      891*p37 + 1458*p38 - 729*p39 - (99*p4)/2 + 162*p40 + 
 
      (891*p41)/2 - 729*p42 + (729*p43)/2 - 81*p44 - 
 
      99*p45 + 162*p46 - 81*p47 + 18*p48 - (1089*p49)/8 - 
 
      (891*p5)/2 + (891*p50)/4 - (891*p51)/8 + 
 
      (99*p52)/4 + (891*p53)/4 - (729*p54)/2 + 
 
      (729*p55)/4 - (81*p56)/2 - (891*p57)/8 + 
 
      (729*p58)/4 - (729*p59)/8 + 729*p6 + (81*p60)/4 + 
 
      (99*p61)/4 - (81*p62)/2 + (81*p63)/4 - (9*p64)/2 - 
 
      (729*p7)/2 + 81*p8 + (891*p9)/4)*u*v*w^2 + 
 
   ((-891*p1)/2 + (3645*p10)/4 - 729*p11 + (729*p12)/4 + 
 
      81*p13 - (405*p14)/2 + 162*p15 - (81*p16)/2 + 
 
      (4455*p17)/4 - (22275*p18)/8 + (4455*p19)/2 + 
 
      (4455*p2)/4 - (4455*p20)/8 - (3645*p21)/2 + 
 
      (18225*p22)/4 - 3645*p23 + (3645*p24)/4 + 
 
      (3645*p25)/4 - (18225*p26)/8 + (3645*p27)/2 - 
 
      (3645*p28)/8 - (405*p29)/2 - 891*p3 + (2025*p30)/4 - 
 
      405*p31 + (405*p32)/4 - 891*p33 + (4455*p34)/2 - 
 
      1782*p35 + (891*p36)/2 + 1458*p37 - 3645*p38 + 
 
      2916*p39 + (891*p4)/4 - 729*p40 - 729*p41 + 
 
      (3645*p42)/2 - 1458*p43 + (729*p44)/2 + 162*p45 - 
 
      405*p46 + 324*p47 - 81*p48 + (891*p49)/4 + 729*p5 - 
 
      (4455*p50)/8 + (891*p51)/2 - (891*p52)/8 - 
 
      (729*p53)/2 + (3645*p54)/4 - 729*p55 + (729*p56)/4 + 
 
      (729*p57)/4 - (3645*p58)/8 + (729*p59)/2 - 
 
      (3645*p6)/2 - (729*p60)/8 - (81*p61)/2 + 
 
      (405*p62)/4 - 81*p63 + (81*p64)/4 + 1458*p7 - 
 
      (729*p8)/2 - (729*p9)/2)*u^2*v*w^2 + 
 
   ((891*p1)/4 - (2187*p10)/4 + (2187*p11)/4 - 
 
      (729*p12)/4 - (81*p13)/2 + (243*p14)/2 - 
 
      (243*p15)/2 + (81*p16)/2 - (4455*p17)/8 + 
 
      (13365*p18)/8 - (13365*p19)/8 - (2673*p2)/4 + 
 
      (4455*p20)/8 + (3645*p21)/4 - (10935*p22)/4 + 
 
      (10935*p23)/4 - (3645*p24)/4 - (3645*p25)/8 + 
 
      (10935*p26)/8 - (10935*p27)/8 + (3645*p28)/8 + 
 
      (405*p29)/4 + (2673*p3)/4 - (1215*p30)/4 + 
 
      (1215*p31)/4 - (405*p32)/4 + (891*p33)/2 - 
 
      (2673*p34)/2 + (2673*p35)/2 - (891*p36)/2 - 
 
      729*p37 + 2187*p38 - 2187*p39 - (891*p4)/4 + 
 
      729*p40 + (729*p41)/2 - (2187*p42)/2 + 
 
      (2187*p43)/2 - (729*p44)/2 - 81*p45 + 243*p46 - 
 
      243*p47 + 81*p48 - (891*p49)/8 - (729*p5)/2 + 
 
      (2673*p50)/8 - (2673*p51)/8 + (891*p52)/8 + 
 
      (729*p53)/4 - (2187*p54)/4 + (2187*p55)/4 - 
 
      (729*p56)/4 - (729*p57)/8 + (2187*p58)/8 - 
 
      (2187*p59)/8 + (2187*p6)/2 + (729*p60)/8 + 
 
      (81*p61)/4 - (243*p62)/4 + (243*p63)/4 - 
 
      (81*p64)/4 - (2187*p7)/2 + (729*p8)/2 + (729*p9)/4)*
 
    u^3*v*w^2 + (81*p1 - (81*p13)/2 - (405*p17)/2 + 
 
      (2025*p21)/4 - 405*p25 + (405*p29)/4 + 162*p33 - 
 
      405*p37 + 324*p41 - 81*p45 - (81*p49)/2 - 
 
      (405*p5)/2 + (405*p53)/4 - 81*p57 + (81*p61)/4 + 
 
      162*p9)*v^2*w^2 + 
 
   ((-891*p1)/2 + 1458*p10 - 729*p11 + 162*p12 + 
 
      (891*p13)/4 - (729*p14)/2 + (729*p15)/4 - 
 
      (81*p16)/2 + (4455*p17)/4 - (3645*p18)/2 + 
 
      (3645*p19)/4 + 729*p2 - (405*p20)/2 - 
 
      (22275*p21)/8 + (18225*p22)/4 - (18225*p23)/8 + 
 
      (2025*p24)/4 + (4455*p25)/2 - 3645*p26 + 
 
      (3645*p27)/2 - 405*p28 - (4455*p29)/8 - (729*p3)/2 + 
 
      (3645*p30)/4 - (3645*p31)/8 + (405*p32)/4 - 
 
      891*p33 + 1458*p34 - 729*p35 + 162*p36 + 
 
      (4455*p37)/2 - 3645*p38 + (3645*p39)/2 + 81*p4 - 
 
      405*p40 - 1782*p41 + 2916*p42 - 1458*p43 + 324*p44 + 
 
      (891*p45)/2 - 729*p46 + (729*p47)/2 - 81*p48 + 
 
      (891*p49)/4 + (4455*p5)/4 - (729*p50)/2 + 
 
      (729*p51)/4 - (81*p52)/2 - (4455*p53)/8 + 
 
      (3645*p54)/4 - (3645*p55)/8 + (405*p56)/4 + 
 
      (891*p57)/2 - 729*p58 + (729*p59)/2 - (3645*p6)/2 - 
 
      81*p60 - (891*p61)/8 + (729*p62)/4 - (729*p63)/8 + 
 
      (81*p64)/4 + (3645*p7)/4 - (405*p8)/2 - 891*p9)*u*
 
    v^2*w^2 + (729*p1 - 3645*p10 + 2916*p11 - 729*p12 - 
 
      (729*p13)/2 + (3645*p14)/4 - 729*p15 + (729*p16)/4 - 
 
      (3645*p17)/2 + (18225*p18)/4 - 3645*p19 - 
 
      (3645*p2)/2 + (3645*p20)/4 + (18225*p21)/4 - 
 
      (91125*p22)/8 + (18225*p23)/2 - (18225*p24)/8 - 
 
      3645*p25 + (18225*p26)/2 - 7290*p27 + (3645*p28)/2 + 
 
      (3645*p29)/4 + 1458*p3 - (18225*p30)/8 + 
 
      (3645*p31)/2 - (3645*p32)/8 + 1458*p33 - 3645*p34 + 
 
      2916*p35 - 729*p36 - 3645*p37 + (18225*p38)/2 - 
 
      7290*p39 - (729*p4)/2 + (3645*p40)/2 + 2916*p41 - 
 
      7290*p42 + 5832*p43 - 1458*p44 - 729*p45 + 
 
      (3645*p46)/2 - 1458*p47 + (729*p48)/2 - 
 
      (729*p49)/2 - (3645*p5)/2 + (3645*p50)/4 - 729*p51 + 
 
      (729*p52)/4 + (3645*p53)/4 - (18225*p54)/8 + 
 
      (3645*p55)/2 - (3645*p56)/8 - 729*p57 + 
 
      (3645*p58)/2 - 1458*p59 + (18225*p6)/4 + 
 
      (729*p60)/2 + (729*p61)/4 - (3645*p62)/8 + 
 
      (729*p63)/2 - (729*p64)/8 - 3645*p7 + (3645*p8)/4 + 
 
      1458*p9)*u^2*v^2*w^2 + 
 
   ((-729*p1)/2 + 2187*p10 - 2187*p11 + 729*p12 + 
 
      (729*p13)/4 - (2187*p14)/4 + (2187*p15)/4 - 
 
      (729*p16)/4 + (3645*p17)/4 - (10935*p18)/4 + 
 
      (10935*p19)/4 + (2187*p2)/2 - (3645*p20)/4 - 
 
      (18225*p21)/8 + (54675*p22)/8 - (54675*p23)/8 + 
 
      (18225*p24)/8 + (3645*p25)/2 - (10935*p26)/2 + 
 
      (10935*p27)/2 - (3645*p28)/2 - (3645*p29)/8 - 
 
      (2187*p3)/2 + (10935*p30)/8 - (10935*p31)/8 + 
 
      (3645*p32)/8 - 729*p33 + 2187*p34 - 2187*p35 + 
 
      729*p36 + (3645*p37)/2 - (10935*p38)/2 + 
 
      (10935*p39)/2 + (729*p4)/2 - (3645*p40)/2 - 
 
      1458*p41 + 4374*p42 - 4374*p43 + 1458*p44 + 
 
      (729*p45)/2 - (2187*p46)/2 + (2187*p47)/2 - 
 
      (729*p48)/2 + (729*p49)/4 + (3645*p5)/4 - 
 
      (2187*p50)/4 + (2187*p51)/4 - (729*p52)/4 - 
 
      (3645*p53)/8 + (10935*p54)/8 - (10935*p55)/8 + 
 
      (3645*p56)/8 + (729*p57)/2 - (2187*p58)/2 + 
 
      (2187*p59)/2 - (10935*p6)/4 - (729*p60)/2 - 
 
      (729*p61)/8 + (2187*p62)/8 - (2187*p63)/8 + 
 
      (729*p64)/8 + (10935*p7)/4 - (3645*p8)/4 - 729*p9)*
 
    u^3*v^2*w^2 + ((-81*p1)/2 + (81*p13)/2 + (405*p17)/4 - 
 
      (1215*p21)/4 + (1215*p25)/4 - (405*p29)/4 - 81*p33 + 
 
      243*p37 - 243*p41 + 81*p45 + (81*p49)/4 + 
 
      (243*p5)/2 - (243*p53)/4 + (243*p57)/4 - 
 
      (81*p61)/4 - (243*p9)/2)*v^3*w^2 + 
 
   ((891*p1)/4 - (2187*p10)/2 + (2187*p11)/4 - 
 
      (243*p12)/2 - (891*p13)/4 + (729*p14)/2 - 
 
      (729*p15)/4 + (81*p16)/2 - (4455*p17)/8 + 
 
      (3645*p18)/4 - (3645*p19)/8 - (729*p2)/2 + 
 
      (405*p20)/4 + (13365*p21)/8 - (10935*p22)/4 + 
 
      (10935*p23)/8 - (1215*p24)/4 - (13365*p25)/8 + 
 
      (10935*p26)/4 - (10935*p27)/8 + (1215*p28)/4 + 
 
      (4455*p29)/8 + (729*p3)/4 - (3645*p30)/4 + 
 
      (3645*p31)/8 - (405*p32)/4 + (891*p33)/2 - 729*p34 + 
 
      (729*p35)/2 - 81*p36 - (2673*p37)/2 + 2187*p38 - 
 
      (2187*p39)/2 - (81*p4)/2 + 243*p40 + (2673*p41)/2 - 
 
      2187*p42 + (2187*p43)/2 - 243*p44 - (891*p45)/2 + 
 
      729*p46 - (729*p47)/2 + 81*p48 - (891*p49)/8 - 
 
      (2673*p5)/4 + (729*p50)/4 - (729*p51)/8 + 
 
      (81*p52)/4 + (2673*p53)/8 - (2187*p54)/4 + 
 
      (2187*p55)/8 - (243*p56)/4 - (2673*p57)/8 + 
 
      (2187*p58)/4 - (2187*p59)/8 + (2187*p6)/2 + 
 
      (243*p60)/4 + (891*p61)/8 - (729*p62)/4 + 
 
      (729*p63)/8 - (81*p64)/4 - (2187*p7)/4 + 
 
      (243*p8)/2 + (2673*p9)/4)*u*v^3*w^2 + 
 
   ((-729*p1)/2 + (10935*p10)/4 - 2187*p11 + 
 
      (2187*p12)/4 + (729*p13)/2 - (3645*p14)/4 + 
 
      729*p15 - (729*p16)/4 + (3645*p17)/4 - 
 
      (18225*p18)/8 + (3645*p19)/2 + (3645*p2)/4 - 
 
      (3645*p20)/8 - (10935*p21)/4 + (54675*p22)/8 - 
 
      (10935*p23)/2 + (10935*p24)/8 + (10935*p25)/4 - 
 
      (54675*p26)/8 + (10935*p27)/2 - (10935*p28)/8 - 
 
      (3645*p29)/4 - 729*p3 + (18225*p30)/8 - 
 
      (3645*p31)/2 + (3645*p32)/8 - 729*p33 + 
 
      (3645*p34)/2 - 1458*p35 + (729*p36)/2 + 2187*p37 - 
 
      (10935*p38)/2 + 4374*p39 + (729*p4)/4 - 
 
      (2187*p40)/2 - 2187*p41 + (10935*p42)/2 - 4374*p43 + 
 
      (2187*p44)/2 + 729*p45 - (3645*p46)/2 + 1458*p47 - 
 
      (729*p48)/2 + (729*p49)/4 + (2187*p5)/2 - 
 
      (3645*p50)/8 + (729*p51)/2 - (729*p52)/8 - 
 
      (2187*p53)/4 + (10935*p54)/8 - (2187*p55)/2 + 
 
      (2187*p56)/8 + (2187*p57)/4 - (10935*p58)/8 + 
 
      (2187*p59)/2 - (10935*p6)/4 - (2187*p60)/8 - 
 
      (729*p61)/4 + (3645*p62)/8 - (729*p63)/2 + 
 
      (729*p64)/8 + 2187*p7 - (2187*p8)/4 - (2187*p9)/2)*
 
    u^2*v^3*w^2 + ((729*p1)/4 - (6561*p10)/4 + 
 
      (6561*p11)/4 - (2187*p12)/4 - (729*p13)/4 + 
 
      (2187*p14)/4 - (2187*p15)/4 + (729*p16)/4 - 
 
      (3645*p17)/8 + (10935*p18)/8 - (10935*p19)/8 - 
 
      (2187*p2)/4 + (3645*p20)/8 + (10935*p21)/8 - 
 
      (32805*p22)/8 + (32805*p23)/8 - (10935*p24)/8 - 
 
      (10935*p25)/8 + (32805*p26)/8 - (32805*p27)/8 + 
 
      (10935*p28)/8 + (3645*p29)/8 + (2187*p3)/4 - 
 
      (10935*p30)/8 + (10935*p31)/8 - (3645*p32)/8 + 
 
      (729*p33)/2 - (2187*p34)/2 + (2187*p35)/2 - 
 
      (729*p36)/2 - (2187*p37)/2 + (6561*p38)/2 - 
 
      (6561*p39)/2 - (729*p4)/4 + (2187*p40)/2 + 
 
      (2187*p41)/2 - (6561*p42)/2 + (6561*p43)/2 - 
 
      (2187*p44)/2 - (729*p45)/2 + (2187*p46)/2 - 
 
      (2187*p47)/2 + (729*p48)/2 - (729*p49)/8 - 
 
      (2187*p5)/4 + (2187*p50)/8 - (2187*p51)/8 + 
 
      (729*p52)/8 + (2187*p53)/8 - (6561*p54)/8 + 
 
      (6561*p55)/8 - (2187*p56)/8 - (2187*p57)/8 + 
 
      (6561*p58)/8 - (6561*p59)/8 + (6561*p6)/4 + 
 
      (2187*p60)/8 + (729*p61)/8 - (2187*p62)/8 + 
 
      (2187*p63)/8 - (729*p64)/8 - (6561*p7)/4 + 
 
      (2187*p8)/4 + (2187*p9)/4)*u^3*v^3*w^2 + 
 
   ((-9*p1)/2 + (27*p17)/2 - (27*p33)/2 + (9*p49)/2)*w^3 + 
 
   ((99*p1)/4 - (297*p17)/4 + (243*p18)/2 - (243*p19)/4 - 
 
      (81*p2)/2 + (27*p20)/2 + (81*p3)/4 + (297*p33)/4 - 
 
      (243*p34)/2 + (243*p35)/4 - (27*p36)/2 - (9*p4)/2 - 
 
      (99*p49)/4 + (81*p50)/2 - (81*p51)/4 + (9*p52)/2)*u*
 
    w^3 + ((-81*p1)/2 + (243*p17)/2 - (1215*p18)/4 + 
 
      243*p19 + (405*p2)/4 - (243*p20)/4 - 81*p3 - 
 
      (243*p33)/2 + (1215*p34)/4 - 243*p35 + (243*p36)/4 + 
 
      (81*p4)/4 + (81*p49)/2 - (405*p50)/4 + 81*p51 - 
 
      (81*p52)/4)*u^2*w^3 + 
 
   ((81*p1)/4 - (243*p17)/4 + (729*p18)/4 - (729*p19)/4 - 
 
      (243*p2)/4 + (243*p20)/4 + (243*p3)/4 + 
 
      (243*p33)/4 - (729*p34)/4 + (729*p35)/4 - 
 
      (243*p36)/4 - (81*p4)/4 - (81*p49)/4 + (243*p50)/4 - 
 
      (243*p51)/4 + (81*p52)/4)*u^3*w^3 + 
 
   ((99*p1)/4 - (9*p13)/2 - (297*p17)/4 + (243*p21)/2 - 
 
      (243*p25)/4 + (27*p29)/2 + (297*p33)/4 - 
 
      (243*p37)/2 + (243*p41)/4 - (27*p45)/2 - 
 
      (99*p49)/4 - (81*p5)/2 + (81*p53)/2 - (81*p57)/4 + 
 
      (9*p61)/2 + (81*p9)/4)*v*w^3 + 
 
   ((-1089*p1)/8 + (729*p10)/4 - (729*p11)/8 + 
 
      (81*p12)/4 + (99*p13)/4 - (81*p14)/2 + (81*p15)/4 - 
 
      (9*p16)/2 + (3267*p17)/8 - (2673*p18)/4 + 
 
      (2673*p19)/8 + (891*p2)/4 - (297*p20)/4 - 
 
      (2673*p21)/4 + (2187*p22)/2 - (2187*p23)/4 + 
 
      (243*p24)/2 + (2673*p25)/8 - (2187*p26)/4 + 
 
      (2187*p27)/8 - (243*p28)/4 - (297*p29)/4 - 
 
      (891*p3)/8 + (243*p30)/2 - (243*p31)/4 + 
 
      (27*p32)/2 - (3267*p33)/8 + (2673*p34)/4 - 
 
      (2673*p35)/8 + (297*p36)/4 + (2673*p37)/4 - 
 
      (2187*p38)/2 + (2187*p39)/4 + (99*p4)/4 - 
 
      (243*p40)/2 - (2673*p41)/8 + (2187*p42)/4 - 
 
      (2187*p43)/8 + (243*p44)/4 + (297*p45)/4 - 
 
      (243*p46)/2 + (243*p47)/4 - (27*p48)/2 + 
 
      (1089*p49)/8 + (891*p5)/4 - (891*p50)/4 + 
 
      (891*p51)/8 - (99*p52)/4 - (891*p53)/4 + 
 
      (729*p54)/2 - (729*p55)/4 + (81*p56)/2 + 
 
      (891*p57)/8 - (729*p58)/4 + (729*p59)/8 - 
 
      (729*p6)/2 - (81*p60)/4 - (99*p61)/4 + (81*p62)/2 - 
 
      (81*p63)/4 + (9*p64)/2 + (729*p7)/4 - (81*p8)/2 - 
 
      (891*p9)/8)*u*v*w^3 + 
 
   ((891*p1)/4 - (3645*p10)/8 + (729*p11)/2 - 
 
      (729*p12)/8 - (81*p13)/2 + (405*p14)/4 - 81*p15 + 
 
      (81*p16)/4 - (2673*p17)/4 + (13365*p18)/8 - 
 
      (2673*p19)/2 - (4455*p2)/8 + (2673*p20)/8 + 
 
      (2187*p21)/2 - (10935*p22)/4 + 2187*p23 - 
 
      (2187*p24)/4 - (2187*p25)/4 + (10935*p26)/8 - 
 
      (2187*p27)/2 + (2187*p28)/8 + (243*p29)/2 + 
 
      (891*p3)/2 - (1215*p30)/4 + 243*p31 - (243*p32)/4 + 
 
      (2673*p33)/4 - (13365*p34)/8 + (2673*p35)/2 - 
 
      (2673*p36)/8 - (2187*p37)/2 + (10935*p38)/4 - 
 
      2187*p39 - (891*p4)/8 + (2187*p40)/4 + 
 
      (2187*p41)/4 - (10935*p42)/8 + (2187*p43)/2 - 
 
      (2187*p44)/8 - (243*p45)/2 + (1215*p46)/4 - 
 
      243*p47 + (243*p48)/4 - (891*p49)/4 - (729*p5)/2 + 
 
      (4455*p50)/8 - (891*p51)/2 + (891*p52)/8 + 
 
      (729*p53)/2 - (3645*p54)/4 + 729*p55 - (729*p56)/4 - 
 
      (729*p57)/4 + (3645*p58)/8 - (729*p59)/2 + 
 
      (3645*p6)/4 + (729*p60)/8 + (81*p61)/2 - 
 
      (405*p62)/4 + 81*p63 - (81*p64)/4 - 729*p7 + 
 
      (729*p8)/4 + (729*p9)/4)*u^2*v*w^3 + 
 
   ((-891*p1)/8 + (2187*p10)/8 - (2187*p11)/8 + 
 
      (729*p12)/8 + (81*p13)/4 - (243*p14)/4 + 
 
      (243*p15)/4 - (81*p16)/4 + (2673*p17)/8 - 
 
      (8019*p18)/8 + (8019*p19)/8 + (2673*p2)/8 - 
 
      (2673*p20)/8 - (2187*p21)/4 + (6561*p22)/4 - 
 
      (6561*p23)/4 + (2187*p24)/4 + (2187*p25)/8 - 
 
      (6561*p26)/8 + (6561*p27)/8 - (2187*p28)/8 - 
 
      (243*p29)/4 - (2673*p3)/8 + (729*p30)/4 - 
 
      (729*p31)/4 + (243*p32)/4 - (2673*p33)/8 + 
 
      (8019*p34)/8 - (8019*p35)/8 + (2673*p36)/8 + 
 
      (2187*p37)/4 - (6561*p38)/4 + (6561*p39)/4 + 
 
      (891*p4)/8 - (2187*p40)/4 - (2187*p41)/8 + 
 
      (6561*p42)/8 - (6561*p43)/8 + (2187*p44)/8 + 
 
      (243*p45)/4 - (729*p46)/4 + (729*p47)/4 - 
 
      (243*p48)/4 + (891*p49)/8 + (729*p5)/4 - 
 
      (2673*p50)/8 + (2673*p51)/8 - (891*p52)/8 - 
 
      (729*p53)/4 + (2187*p54)/4 - (2187*p55)/4 + 
 
      (729*p56)/4 + (729*p57)/8 - (2187*p58)/8 + 
 
      (2187*p59)/8 - (2187*p6)/4 - (729*p60)/8 - 
 
      (81*p61)/4 + (243*p62)/4 - (243*p63)/4 + 
 
      (81*p64)/4 + (2187*p7)/4 - (729*p8)/4 - (729*p9)/8)*
 
    u^3*v*w^3 + ((-81*p1)/2 + (81*p13)/4 + (243*p17)/2 - 
 
      (1215*p21)/4 + 243*p25 - (243*p29)/4 - (243*p33)/2 + 
 
      (1215*p37)/4 - 243*p41 + (243*p45)/4 + (81*p49)/2 + 
 
      (405*p5)/4 - (405*p53)/4 + 81*p57 - (81*p61)/4 - 
 
      81*p9)*v^2*w^3 + ((891*p1)/4 - 729*p10 + 
 
      (729*p11)/2 - 81*p12 - (891*p13)/8 + (729*p14)/4 - 
 
      (729*p15)/8 + (81*p16)/4 - (2673*p17)/4 + 
 
      (2187*p18)/2 - (2187*p19)/4 - (729*p2)/2 + 
 
      (243*p20)/2 + (13365*p21)/8 - (10935*p22)/4 + 
 
      (10935*p23)/8 - (1215*p24)/4 - (2673*p25)/2 + 
 
      2187*p26 - (2187*p27)/2 + 243*p28 + (2673*p29)/8 + 
 
      (729*p3)/4 - (2187*p30)/4 + (2187*p31)/8 - 
 
      (243*p32)/4 + (2673*p33)/4 - (2187*p34)/2 + 
 
      (2187*p35)/4 - (243*p36)/2 - (13365*p37)/8 + 
 
      (10935*p38)/4 - (10935*p39)/8 - (81*p4)/2 + 
 
      (1215*p40)/4 + (2673*p41)/2 - 2187*p42 + 
 
      (2187*p43)/2 - 243*p44 - (2673*p45)/8 + 
 
      (2187*p46)/4 - (2187*p47)/8 + (243*p48)/4 - 
 
      (891*p49)/4 - (4455*p5)/8 + (729*p50)/2 - 
 
      (729*p51)/4 + (81*p52)/2 + (4455*p53)/8 - 
 
      (3645*p54)/4 + (3645*p55)/8 - (405*p56)/4 - 
 
      (891*p57)/2 + 729*p58 - (729*p59)/2 + (3645*p6)/4 + 
 
      81*p60 + (891*p61)/8 - (729*p62)/4 + (729*p63)/8 - 
 
      (81*p64)/4 - (3645*p7)/8 + (405*p8)/4 + (891*p9)/2)*
 
    u*v^2*w^3 + ((-729*p1)/2 + (3645*p10)/2 - 1458*p11 + 
 
      (729*p12)/2 + (729*p13)/4 - (3645*p14)/8 + 
 
      (729*p15)/2 - (729*p16)/8 + (2187*p17)/2 - 
 
      (10935*p18)/4 + 2187*p19 + (3645*p2)/4 - 
 
      (2187*p20)/4 - (10935*p21)/4 + (54675*p22)/8 - 
 
      (10935*p23)/2 + (10935*p24)/8 + 2187*p25 - 
 
      (10935*p26)/2 + 4374*p27 - (2187*p28)/2 - 
 
      (2187*p29)/4 - 729*p3 + (10935*p30)/8 - 
 
      (2187*p31)/2 + (2187*p32)/8 - (2187*p33)/2 + 
 
      (10935*p34)/4 - 2187*p35 + (2187*p36)/4 + 
 
      (10935*p37)/4 - (54675*p38)/8 + (10935*p39)/2 + 
 
      (729*p4)/4 - (10935*p40)/8 - 2187*p41 + 
 
      (10935*p42)/2 - 4374*p43 + (2187*p44)/2 + 
 
      (2187*p45)/4 - (10935*p46)/8 + (2187*p47)/2 - 
 
      (2187*p48)/8 + (729*p49)/2 + (3645*p5)/4 - 
 
      (3645*p50)/4 + 729*p51 - (729*p52)/4 - 
 
      (3645*p53)/4 + (18225*p54)/8 - (3645*p55)/2 + 
 
      (3645*p56)/8 + 729*p57 - (3645*p58)/2 + 1458*p59 - 
 
      (18225*p6)/8 - (729*p60)/2 - (729*p61)/4 + 
 
      (3645*p62)/8 - (729*p63)/2 + (729*p64)/8 + 
 
      (3645*p7)/2 - (3645*p8)/8 - 729*p9)*u^2*v^2*w^3 + 
 
   ((729*p1)/4 - (2187*p10)/2 + (2187*p11)/2 - 
 
      (729*p12)/2 - (729*p13)/8 + (2187*p14)/8 - 
 
      (2187*p15)/8 + (729*p16)/8 - (2187*p17)/4 + 
 
      (6561*p18)/4 - (6561*p19)/4 - (2187*p2)/4 + 
 
      (2187*p20)/4 + (10935*p21)/8 - (32805*p22)/8 + 
 
      (32805*p23)/8 - (10935*p24)/8 - (2187*p25)/2 + 
 
      (6561*p26)/2 - (6561*p27)/2 + (2187*p28)/2 + 
 
      (2187*p29)/8 + (2187*p3)/4 - (6561*p30)/8 + 
 
      (6561*p31)/8 - (2187*p32)/8 + (2187*p33)/4 - 
 
      (6561*p34)/4 + (6561*p35)/4 - (2187*p36)/4 - 
 
      (10935*p37)/8 + (32805*p38)/8 - (32805*p39)/8 - 
 
      (729*p4)/4 + (10935*p40)/8 + (2187*p41)/2 - 
 
      (6561*p42)/2 + (6561*p43)/2 - (2187*p44)/2 - 
 
      (2187*p45)/8 + (6561*p46)/8 - (6561*p47)/8 + 
 
      (2187*p48)/8 - (729*p49)/4 - (3645*p5)/8 + 
 
      (2187*p50)/4 - (2187*p51)/4 + (729*p52)/4 + 
 
      (3645*p53)/8 - (10935*p54)/8 + (10935*p55)/8 - 
 
      (3645*p56)/8 - (729*p57)/2 + (2187*p58)/2 - 
 
      (2187*p59)/2 + (10935*p6)/8 + (729*p60)/2 + 
 
      (729*p61)/8 - (2187*p62)/8 + (2187*p63)/8 - 
 
      (729*p64)/8 - (10935*p7)/8 + (3645*p8)/8 + (729*p9)/2
 
      )*u^3*v^2*w^3 + ((81*p1)/4 - (81*p13)/4 - 
 
      (243*p17)/4 + (729*p21)/4 - (729*p25)/4 + 
 
      (243*p29)/4 + (243*p33)/4 - (729*p37)/4 + 
 
      (729*p41)/4 - (243*p45)/4 - (81*p49)/4 - 
 
      (243*p5)/4 + (243*p53)/4 - (243*p57)/4 + 
 
      (81*p61)/4 + (243*p9)/4)*v^3*w^3 + 
 
   ((-891*p1)/8 + (2187*p10)/4 - (2187*p11)/8 + 
 
      (243*p12)/4 + (891*p13)/8 - (729*p14)/4 + 
 
      (729*p15)/8 - (81*p16)/4 + (2673*p17)/8 - 
 
      (2187*p18)/4 + (2187*p19)/8 + (729*p2)/4 - 
 
      (243*p20)/4 - (8019*p21)/8 + (6561*p22)/4 - 
 
      (6561*p23)/8 + (729*p24)/4 + (8019*p25)/8 - 
 
      (6561*p26)/4 + (6561*p27)/8 - (729*p28)/4 - 
 
      (2673*p29)/8 - (729*p3)/8 + (2187*p30)/4 - 
 
      (2187*p31)/8 + (243*p32)/4 - (2673*p33)/8 + 
 
      (2187*p34)/4 - (2187*p35)/8 + (243*p36)/4 + 
 
      (8019*p37)/8 - (6561*p38)/4 + (6561*p39)/8 + 
 
      (81*p4)/4 - (729*p40)/4 - (8019*p41)/8 + 
 
      (6561*p42)/4 - (6561*p43)/8 + (729*p44)/4 + 
 
      (2673*p45)/8 - (2187*p46)/4 + (2187*p47)/8 - 
 
      (243*p48)/4 + (891*p49)/8 + (2673*p5)/8 - 
 
      (729*p50)/4 + (729*p51)/8 - (81*p52)/4 - 
 
      (2673*p53)/8 + (2187*p54)/4 - (2187*p55)/8 + 
 
      (243*p56)/4 + (2673*p57)/8 - (2187*p58)/4 + 
 
      (2187*p59)/8 - (2187*p6)/4 - (243*p60)/4 - 
 
      (891*p61)/8 + (729*p62)/4 - (729*p63)/8 + 
 
      (81*p64)/4 + (2187*p7)/8 - (243*p8)/4 - (2673*p9)/8)*
 
    u*v^3*w^3 + ((729*p1)/4 - (10935*p10)/8 + 
 
      (2187*p11)/2 - (2187*p12)/8 - (729*p13)/4 + 
 
      (3645*p14)/8 - (729*p15)/2 + (729*p16)/8 - 
 
      (2187*p17)/4 + (10935*p18)/8 - (2187*p19)/2 - 
 
      (3645*p2)/8 + (2187*p20)/8 + (6561*p21)/4 - 
 
      (32805*p22)/8 + (6561*p23)/2 - (6561*p24)/8 - 
 
      (6561*p25)/4 + (32805*p26)/8 - (6561*p27)/2 + 
 
      (6561*p28)/8 + (2187*p29)/4 + (729*p3)/2 - 
 
      (10935*p30)/8 + (2187*p31)/2 - (2187*p32)/8 + 
 
      (2187*p33)/4 - (10935*p34)/8 + (2187*p35)/2 - 
 
      (2187*p36)/8 - (6561*p37)/4 + (32805*p38)/8 - 
 
      (6561*p39)/2 - (729*p4)/8 + (6561*p40)/8 + 
 
      (6561*p41)/4 - (32805*p42)/8 + (6561*p43)/2 - 
 
      (6561*p44)/8 - (2187*p45)/4 + (10935*p46)/8 - 
 
      (2187*p47)/2 + (2187*p48)/8 - (729*p49)/4 - 
 
      (2187*p5)/4 + (3645*p50)/8 - (729*p51)/2 + 
 
      (729*p52)/8 + (2187*p53)/4 - (10935*p54)/8 + 
 
      (2187*p55)/2 - (2187*p56)/8 - (2187*p57)/4 + 
 
      (10935*p58)/8 - (2187*p59)/2 + (10935*p6)/8 + 
 
      (2187*p60)/8 + (729*p61)/4 - (3645*p62)/8 + 
 
      (729*p63)/2 - (729*p64)/8 - (2187*p7)/2 + 
 
      (2187*p8)/8 + (2187*p9)/4)*u^2*v^3*w^3 + 
 
   ((-729*p1)/8 + (6561*p10)/8 - (6561*p11)/8 + 
 
      (2187*p12)/8 + (729*p13)/8 - (2187*p14)/8 + 
 
      (2187*p15)/8 - (729*p16)/8 + (2187*p17)/8 - 
 
      (6561*p18)/8 + (6561*p19)/8 + (2187*p2)/8 - 
 
      (2187*p20)/8 - (6561*p21)/8 + (19683*p22)/8 - 
 
      (19683*p23)/8 + (6561*p24)/8 + (6561*p25)/8 - 
 
      (19683*p26)/8 + (19683*p27)/8 - (6561*p28)/8 - 
 
      (2187*p29)/8 - (2187*p3)/8 + (6561*p30)/8 - 
 
      (6561*p31)/8 + (2187*p32)/8 - (2187*p33)/8 + 
 
      (6561*p34)/8 - (6561*p35)/8 + (2187*p36)/8 + 
 
      (6561*p37)/8 - (19683*p38)/8 + (19683*p39)/8 + 
 
      (729*p4)/8 - (6561*p40)/8 - (6561*p41)/8 + 
 
      (19683*p42)/8 - (19683*p43)/8 + (6561*p44)/8 + 
 
      (2187*p45)/8 - (6561*p46)/8 + (6561*p47)/8 - 
 
      (2187*p48)/8 + (729*p49)/8 + (2187*p5)/8 - 
 
      (2187*p50)/8 + (2187*p51)/8 - (729*p52)/8 - 
 
      (2187*p53)/8 + (6561*p54)/8 - (6561*p55)/8 + 
 
      (2187*p56)/8 + (2187*p57)/8 - (6561*p58)/8 + 
 
      (6561*p59)/8 - (6561*p6)/8 - (2187*p60)/8 - 
 
      (729*p61)/8 + (2187*p62)/8 - (2187*p63)/8 + 
 
      (729*p64)/8 + (6561*p7)/8 - (2187*p8)/8 - (2187*p9)/8
 
      )*u^3*v^3*w^3);

One might ask why bother with the cubic solid when the built-in Interpolation function is available. In spite of the huge length of the preceeding expression for the cubic solid, it evaluates 6 times faster when compiled than the corresponding Interpolate3D function. Even when the cubic solid is evaluated uncompiled it is still twice as fast, so the advantage from using Interpolate3D comes into play when you need to build curves with more than 4 defining points, or surfaces from more than 16 points, or solids from more than 64 points.

2.5  Nodes and Elements

Next we present how nodes and elements are stored, created, and manipulated. Each finite element solver has its own unique node and element types, and to allow for easy introduction of different solvers, or adding new elements, the node and element descriptions in ModelMaker are not hard-wired, but are loaded in node and element type databases (called nodetypes and elementtypes) with a NEW command, just like the geometric entities. The list nodetypes stores an integer ID for the type, the name of the analysis code, short and long names for the node type, and the number of degrees of freedom. The list elementtypes stores an integer ID for the type, the name of the analysis code, short and long names for the element type, the geometry that the element can mesh, nodal coordinates, instructions on how to fill the unit cell (the unit cube in 3D, the unit square in 2D) with the minimum number of elements that fill the cell, plotting instructions and a list of required node types. Here are the element type commands (the nodetype commands are similar):

elementtype[NEXT];
elementtype[NEW,{Analysis,shortname,fullname,
 geometry,data,edges,fills,nodetype},id];
elementtype[NEW,{Analysis,shortname,fullname,
 geometry,data,edges,fills,nodetype}];
elementtype[DELETE];
elementtype[DELETE,id];
elementtype[DELETE,{{id1,id2,id3}}];
elementtype[DELETE,ids];

ModelMaker's mesh command uses uniform parametric meshing - nodes are created in the parameter space of each geometric entity in the model and then these are mapped into 3-space using the interpolation function of the entities. The mesh command also creates the list elements by calling the element function. This is followed by nodal equivalencing to tie all the elements together. The node types and element types required for a solver named e.g., Analysis, get loaded by the command 

load[Analysis]

which must be defined for each solver. Currently values for Analysis of only "mavart" and "mavart3d" are defined. The nodes are stored in a data structure similar to points, with additional fields for fixities, and loads. To apply boundary conditions to the nodes, suitably named wrappers, like node[VOLTAGE,v,ids] invoke the more generic node command node[FIXITY,fixities,data,ids] to apply voltage fixities to e.g., MAVART elements. These wrappers will have to be customized for various analysis codes and included in the definition of the load command.

Plotting of nodes and elements relies on the same techniques previously developed for geometry entities, with the addition of symbology to show boundary conditions and apply deformations read back from the analysis results.

Here are the node, element and mesh commands:

node[];
node[NEW];
node[NEW,{x,y,z}];
node[NEW,{x,y,z},nodetype];
node[NEW,{x,y,z},nodetype,fixity,data];
node[NEW,{x,y,z},nodetype,fixity,data,id];
node[DELETE,id];
node[DELETE,{{id1,id2,id3}}];
node[DELETE];
node[PLOT];
node[PLOT,,option];
node[PLOT,ids,option];
node[PLOT,{{id1,id2,id3}},option];
node[ROTATE,ids,{r,p,y},{x0,y0,z0},n];
node[TRANSLATE,ids,{dx,dy,dz},n];
node[MIRROR,plane,ids];
node[TYPE,nodetype,ids]; (* assign node type to ids *)
node[FIXITY,jfix,valfaz,ids]; (* assign fixity flags and complex values of fixities to the list of nodes ids *)
node[FORCE,{fx,fy,fz},ids]; (* apply a complex force to ids *)
node[DISPLACEMENT,{dx,dy,dz},ids]; (* set complex displacement fixities to ids *)
node[VOLTAGE,v,ids]; (* apply a voltage fixity v to ids *)
node[PRESSURE,p,ids]; (* apply a pressure fixity p to ids *)
node[EQUIVALENCE];
node[EQUIVALENCE,ids];
node[RENUMBER];
node[ELEMLIST,ids]; (* return list of nodes belonging to element list ids *)
node[NEXT];
node[LIST,id1,id2,id3];
node[UNUSED];  (* return list of nodes not used by any element *)
node[SYMMETRIZE,axis]; (* list nodes on n-fold symmetry boundaries *)

element[NEW,elementtypeid,materialid,nodeids,data];
element[NEW,elementtypeid,materialid,nodeids,data,id];
element[DELETE];
element[DELETE,id];
element[DELETE,{{id1,id2,id3}}];
element[DELETE,ids];
element[PLOT];
element[PLOT,,option];
element[PLOT,id,option];
element[PLOT,{{id1,id2,id3}},option];
element[PLOT,ids,option];
element[TRANSLATE,ids,{dx,dy,dz},n];
element[ROTATE,ids,{r,p,y},{x0,y0,z0},n];
element[MATERIAL,id,ids]; (* assign material id to elements ids *)
element[PROPERTY,data,ids]; (* assign element properties to ids *)
element[RENUMBER];
element[NEXT];

mesh[id,ids,{nx},option];      (* mesh with elementtype id, the curves ids *)
mesh[id,ids,{nx,ny},option];   (* mesh with elementtype id, the surfaces ids *)
mesh[id,ids,{nx,ny,nz},option];(* mesh with elementtype id, the solids ids *)

load[Analysis]; (* loads element types and mesh definitions 
for an analysis type. *)
fem[DELETE];  (* delete all finite element data *)

2.6  Materials

Materials are created and stored with a similar database design as all the other components. Presently only a minimal set of functions have been developed, and these will need customization to support a variety of solvers.

material[NEXT];
material[NEW,type,density,thermal,stiffness,electroelastic,dielectric];
material[NEW,type,density,thermal,stiffness,electroelastic,dielectric,id];
material[ISOTROPIC,type,density,YoungsModulus,PoissonRatio];
material[DELETE];
material[DELETE,id];
material[DELETE,{{id1,id2,id3}}];
material[DELETE,ids];

2.7  Input and Output

Lastly, the file input and output routines are described. Finite element files are usually fixed Fortran format, and require translation into a form Mathematica can handle. The translation tools all exist within Mathematica, but ModelMaker's translation routine is slow and this section of the Package needs a cleanup. These file I/O routines will have to be customized for different FE codes since there are no universal standards for file structures for finite elements. The readresult function reads the PATRAN neutral nodal displacement results file format. To display results they are loaded into the nodes database as fixities. The element plotting routines can then be used to plot displacement results. Mathematica has a built-in animation capability so dynamic visualization of displacement results is possible.  Here are the file I/O routines developed to date:

read[];              (* read a MAVART3D data file *)
readresult[];        (* read a PATRAN neutral nodeal displacement file *)
write[title];        (* write a MAVART3D data file *)
printmavart[];       (* print a MAVART3D data file on the screen *)

To bring all this together we now use ModelMaker to build some finite element models.

3.  Examples

3.1  Morphing Simple Beam

This example shows the construction of a simple beam fixed at one end, with a force of 4 N in the -Z direction applied to the free end. The beam's shape can be changed or morphed,  even after the model is completely assembled. The parameters of the beam namely height, width, and length are inserted into the model with Hold wrapped around them. Functions like EQUIVALENCE work with the Holds left in place. Using Hold makes it particularly easy to do parametric modelling. You can change the values assigned to height, width, and length at any time and not permanently affect the database. Without the Holds you would have to be careful not to assign values to the symbols before the model was complete.

Clear[height,width,length];
geometry[DELETE];
point[NEW,{0,0,0},1];
point[TRANSLATE,{1},{0,0,Hold[height]},1];
curve[NEW,{1,2}];
curve[TRANSLATE,{1},{Hold[length],0,0},1];
surface[NEW,{1,2}];
surface[TRANSLATE,{1},{0,Hold[width],0},1];
solid[NEW,{1,2}];

If we look at the points database, we see the held symbols for the beam dimensions are embedded, and we can apply ReleaseHold at any time.

points

To plot the solid we must set some values for height,width and length. The PLOT tasking applies ReleaseHold to the database before plotting. Here we set the parameters to turn the beam into a cube.

{height,width,length}={1,1,1};
solid[PLOT];

Next we load the definitions needed for our analysis type (MAVART3D), mesh the beam with 5 hexahedral (MAVART3D type 1) elements, and examine the nodes. The mesh task does not apply ReleaseHold.

fem[DELETE];
load["mavart3d"];
mesh[1,{1},{1,5,1},Compiled->False];

nodes//Short

We see the held values have propagated through to the nodes and so if the values assigned to the symbols are changed, the output forms of the elements will immediately change, as the following shows:

{height,width,length}={.5,.5,2};
element[PLOT];

(* Next apply a unit force to one end of the beam and 
fix the other, setting material properties of aluminum *)
node[FORCE,{0,0,-1},{21,22,23,24}];
node[DISPLACEMENT,{0,0,0},{1,2,5,6}];
material[DELETE];
material[ISOTROPIC,1,2700,100 10^9,.31];
element[MATERIAL,1,elementlist];

(* Set some fancy global plot symbols, so the fixities
and loads are now shown.*)
nodesize=.02;
nonzerofixity=True; 
zerodispfixity=True;
vectorsize=.5;conesize=.05;
elementshrink=10;elementform="FILL";

(* Change the shape, and the fixities and loads track the logical definitions (node IDs), and are not affected by the change in node coordinates *)
{height,width,length}={.2,1,5};
element[PLOT];

3.2  Morphing Airplane

Here is a more interesting morphing example than the simple beam. This is a model of an airplane complete with parameterized wing span, sweep, dihedral, fusilage length, width, tail geometry, etc. In the notebook version of this document you can animate this image to show a series of designs as all these parameters are varied. Here the full benefit of the parametric approach to modelling is apparent.

3.3  Boat Hull

This example shows a shell model of a boat hull built using curves and surfaces defined with Interpolate3D. Using the parametric cubic equation as the interpolation function did not produce fair curves through the defining points - because more than 4 points were needed to describe the hull cross sections accurately. In the notebook version of this document this image can be rotated.

3.4  Modal Analysis Sonar Projector Stave

This example shows that results from an external solver can be imported back into Mathematica using the readresult function and then animated. In the notebook version of this document you can double click the image below and see the lowest torsional mode of a single stave of the barrel stave projector, as calculated by MAVART3D.

4.  Future Directions for the Code

A number of necessary improvements to ModelMaker have been identified throughout the preceding sections.  Future enhancements needed include reducing dependency on global variables, making the mesh routines more capable, improving the file I/O performance,  and adding a graphical user interface. Mathematica 3.0 [8] offers tools for building buttons and palettes, and a palettized version of ModelMaker has already been attempted. Designing an effective GUI will take considerable effort with the payback largely for users unfamiliar with Mathematica. Experienced users may prefer the command language approach in the Mathematica Front End environment.

5.  Conclusions

This project started from a requirement for a low cost, portable 3D finite element modeller compatible with MAVART3D, and Mathematica was ideal for turning these requirements into a prototype. The prototype was developed in far less time than would have been needed had conventional tools (C, C++, X) been used. The ModelMaker package runs fast enough on personal computers that it is clear  Mathematica can host a finite element modeller suitable for sonar transducer design. It should be noted that ModelMaker can produce 2D models as easily as 3D, and thus can also support the axisymmetric MAVART.  How far beyond simple models ModelMaker  can be applied is hard to predict from the limited testing done to date. Considering that the present code has not been optimized for speed, and that the speed of personal computers is increasing rapidly, the ModelMaker package has considerable growth potential, probably enough to justify commercial development.

6.  References

1. MacNeil Schwendler Corp., PATRAN, Los Angeles, CA.
2. Parametric Technology Inc., Pro/ENGINEER, Waltham, MA.
3. Wolfram Research Inc., Mathematica, Version 2.2, Champaign, IL(1993).
4. Roman Maeder, The Mathematica Programmer, AP Professional, Cambridge, MA, 1994.
5. Roman Maeder, Programming in Mathematica, Addison-Wesley, NY, 1990.
6. W. Press et al., Numerical Recipes, The Art of Scientific Computing, Cambridge Univ. Press, NY, 1989.
7. Xah Lee, Distinct Elements, http://www.best.com/~xah/
8. Wolfram Research Inc., Mathematica, Version 3.0, Champaign, Illinois (1996).