# Using the pre-processor

Preliminaries

## Preliminaries

Before starting to use PREFEL, you must define the problem, and draw up a coarse finite element mesh. The first thing to do, is to decide which type of problem you are interested in using FELIPE to solve. Is your problem one in which the primary unknown quantity is a vector with x and y components (for example a field of displacements or flow velocities), or one in which the primary unknown is a scalar variable (such as temperature or pressure)? FELIPE uses the elasticity main engine' ELAST.FOR as an example of the vector' class of problems, and the main engine' POISS.FOR (solving Poisson's equation over a 2D region) to illustrate the scalar' class of problems. In the analysis of beams, there are three primary unknowns: the x and y displacements and the rotation at each node.

New users of FELIPE should first practise on solving Poisson, elasticity and beam problems at Basic level, before moving to develop their own main engine' and produce data files for it at Advanced level. A brief description of each standard problem now follows:

• Elasticity. In elasticity problems, the material properties which must be defined for each element, are the Young's modulus E, the Poisson's ratio \nu and the material thickness t. Loading consists of point loads applied at nodes, and surface tractions applied on element boundary edges. The mesh represents a linear elastic material in either plane stress (zero out-of-plane stress) or plane strain ( i.e. no movement in the out-of-plane direction).

• Poisson. Here, we seek to solve numerically the Poisson equation (a.k.a. the quasi-harmonic equation):

$$- a_x \frac{\partial^2 u}{\partial x^2} - a_y \frac{\partial^2 u}{\partial y^2} = f(x,y)$$

over a 2D region, with appropriate boundary conditions; these may be:

• Dirichlet boundary: fixed values of U on the boundary;
• reflective boundary: zero gradient or flow across the boundary;
• radiating boundary: flow fixed or proportional to the temperature.
The parameters a_x, a_y must be defined for each element, in the same way as the material parameters E, \nu in elasticity problems. The right-hand-side function f(x,y) will be defined in a FUNCTION routine in the program POISS.FOR.

• Beam. The axial stiffness of the beam is defined by its Young's modulus E, and its cross-sectional area A. Its flexural stiffness is defined by the Moment of Inertia I. Apart from the material properties, the preparation of a mesh for a beam or frame problem is a simplified version of the procedure for 2D elasticity problems, and will not be dealt with separately.

The next step is to draw up a coarse mesh. If you will be using the programs ELAST.FOR or POISS.FOR , then the mesh must consist of 8-noded quadrilateral elements (for elasticity problems) or of 3-noded linear triangles (for Poisson problems) respectively. In the latter case, however, it is possible to create a mesh of 4-noded linear quadrilaterals, which can be refined (see \S2.4.3), and finally triangulated (cutting each quadrilateral into two triangles). For BEAM.FOR, each member of the frame is represented by a two-noded beam element. The range of elements available is illustrated in the Appendix.

Determine the coordinates of all nodes in your mesh. Note that if using 8-noded quadrilaterals, you can create elements with curved sides, when the three nodes making up the side are not colinear. The curve is determined by the coordinates of the three nodes, and the quadratic basis functions in one dimension.

Number all the elements, and all the nodes; note that the numbering will be changed later, in the pre-processor.

Values of the parameters E, \nu, t (for Elasticity problems), a_x and a_y (for Poisson problems) must be associated with each element. These values will henceforth be referred to as the material properties of the element.

The mesh boundary conditions must also be specified --- that is, the parts of the boundary which are fixed in one or both directions (for elasticity) or which have Dirichlet boundary conditions (for Poisson problems) must be specified --- in the main mesh input section.

The above data defining the mesh will be referred to as Part A of the full data required for input to a main engine’. Part B data consists of the applied loads (for elasticity -- see \S2.6), any boundary inputs (e.g. specified values of U on Dirichlet boundary planes) and any incremental loading or timestepping regime data.

The recommended steps for preparing a full data file (which have been followed in the examples provided in the package), are:

1. produce a file containing the Part A data describing the coarse mesh, at Basic level, and save it with a filename such as test1.dat;
2. modify the Part A data by mesh refinement, element and node renumbering, etc. (if necessary changing to Advanced level), to produce the Part A data describing the final mesh, in a file named test2.dat;
3. add Part B data for the desired applied loading, saving the full data file as test3.dat, and run the appropriate main engine’ with this as the input file;
4. Use the pre-processor to add new Part B loads to test2.dat --- or to edit the Part B loads in test3.dat --- producing new datafiles test4.dat and so on, to see the effects of different load regimes.
It is however possible to proceed to add Part B data directly to the mesh input in test1.dat, if desired.

Note that whenever using the pre-processor to read in and edit an existing data file, the new file being produced must be given a different name. This is because PREFEL reads the input file and writes to the output file simultaneously.}

As mentioned above, once the coarse mesh has been defined, it will be displayed graphically and the user can refine it in various ways. The element and node numbering of the resulting mesh can then be re-numbered to reduce the bandwidth of the associated global stiffness matrix.

## Creating a data file

Run the pre-processor PREFEL by double-clicking on its icon in the C:\FELIPE directory.

A Text dialogue window will appear. This is used for text input. Note that in answering yes/no questions, it suffices for you to press the y or n keys; you do not need to press Return or Enter. The same applies in instances where only a single digit is expected. Also note that if you press the Esc key in answer to a question, you will be able to exit the program without saving. It is worth pointing out that when running the compiled main engines’, you do need to press Return in answer to questions. This is because the main engines’ use only standard Fortran77 which can be compiled by any compiler, and this does not include a getkey’’ routine for intercepting keyboard input.

Once some basic data has been entered in the Text window, a Graphics window will appear. Dialogue will then switch between the Text and Graphics windows; you will be prompted when to bring the appropriate window to the top of the desktop (by clicking on its caption bar). In the Graphics window the available options will appear as a menu on the right-hand side of the window. You choose a menu item by left- clicking on it; if you right-click on a menu item a Help screen will appear. Next, you will be asked if you wish to create a new .dat file; answering y will invoke the CREATE subroutine, which will now be described. (By answering n you can modify an existing file, as described in the next section)

First, you give the name for the data file you will create (Do not type the .dat filename extension). The name should have a maximum of 8 characters (excluding the .dat extension); if longer, it will be truncated. A useful convention is to end it with a digit indicating if it is the 1st, 2nd, etc. mesh in the problem definition process, e.g. myprob1.dat. It must be a valid MS-DOS name - no spaces or non-alphanumeric characters.

By default, the datafile will be created in the same folder as that in which the pre-processor is stored. You will be asked if this is OK; if you answer n' you will be able to select a new folder from a Windows selection box.

Next, you can enter an identifying title for the mesh (max. 50 characters, including spaces).

Now you specify the type of problem to be solved, as described above:

• type 1 for a scalar-variable problem: at Basic level this will be a Poisson analysis;
• or type 2 for a vector-variable problem: an elasticity analysis at Basic level.
• or type 3 for a beam or frame problem.
(The simplest analysis to start with, is a Poisson problem). Note that the number you type (1, 2 or 3) indicates the number of basic unknowns at each node, used in creating the main graphical display in the post-processor.

You then specify whether you will be performing a Basic or Advanced level analysis. The main engines' ELAST.FOR, POISS.FOR and FRAME.FOR are written for Basic-level analyses only, so if you plan to use one of these main engines' you should reply b. The main additional features of Advanced-level analyses are that all the 1D and 2D element types can be used in the mesh, and there is a wider range of material properties and loading r\'egimes in elasticity analyses. A student exercise might be to modify the ELAST.FOR program to use the four-noded linear quadrilateral element; meshes with this element would be prepared using the Advanced option. Type a for Advanced or b for Basic -- new users should use Basic level.

The type of analysis (scalar or vector variable), and its level (Basic or Advanced), is stored in an integer variable named NDOFN which denotes the number of degrees of freedom per node. In a vector-variable mesh there are two degrees of freedom per node (in elasticity the x- and y- displacements) while in a scalar-variable mesh there is only one degree of freedom per node: the potential U. In a frame, there are three degrees of freedom: the x- and y- displacements, plus the rotation of the beam at the node.

An Advanced-level mesh is denoted by setting NDOFN to -1 or -2 or -3 rather than 1 or 2 or 3. The value of NDOFN can be seen in the first line of data in the resulting .dat file. (At Advanced level, additional degrees of freedom can be defined at a later stage) At advanced level, you can also define additional parameters to use in your own customized 'main engines'.

You now have the option of defining the basic mesh data (nodal coordinates and element-node connections) either graphically using the mouse, or by keyboard entry.

If you have chosen the keyboard option for describing the mesh, you will be prompted to enter the total number of nodes and the total number of elements in the finite element mesh; these integers are stored in NPOIN and NELEM. These are also written in the first data line of the .dat file --- see \S2.8.

Then you enter the x and y coordinates of each node in turn (Type the two coordinates, separated by a space or comma, and press Return). Once all the nodal coordinates are entered, the screen will list them, and give you the chance to correct any errors. Then it will open a Graphics window, and plot the nodes in the x,y-plane so that you can check them. A pop-up box asks if you want to make any changes, i.e. to alter any of the coordinates.

The focus reverts to the Text window, either for you to correct a nodal coordinate or to proceed to list, for each element in turn, the node numbers of the nodes around the element. To see the Text window you will need to click on its caption bar to bring it to the foreground. During the mesh-creation process, you will be prompted when you need to click on the Text or Graphics windows to make them visible.

It is important that the nodes around each element should be listed anticlockwise, starting at a corner of the element, as illustrated in Appendix A. A full list of the element types available is given in the Table below, but if you are performing a Basic-level analysis, the pre-processor will only permit you to use the corresponding element types: linear triangles or quadrilaterals (to be triangulated later) for Poisson problems (NDOFN=1), 8-noded quadrilaterals for elasticity (NDOFN=2), 2-noded line elements for beams (NDOFN=3).

 NSIDE NNODE Element type 1 2 2-noded 1D line element 1 3 3-noded 1D line element 2 6 Quadratic joint element 2 4 infinite joint element 3 3 Linear triangle 4 4 Linear quadrilateral 3 6 Quadratic (6-noded) triangle 4 6 Thin (6-noded) quadrilateral 4 8 Serendipity (8-noded) quadrilateral 3 4 Linear infinite element 3 5 Quadratic infinite element 2 3 Corner infinite element
For each element in turn, enter NSIDE and NNODE --- as each of these will be a single digit, you do not need to press Return --- followed by the node numbers around the element (separated by commas or spaces). You will be asked to enter a digit specifying the material type LMTYPE of the element. This is independent of the material properties, and can be used to indicate the type of behaviour (e.g. 1=elastic, 2=elasto-plastic). LMTYPE is used in the consolidation main engine' CONSL.FOR, where LMTYPE=1 for structural elements and LMTYPE=2 for soil elements (with pore pressure degrees of freedom). For Basic level analyses, just enter 1 for LMTYPE.

When all elements have been described, the element-node connections will be listed. You can make changes to correct errors, and the data is also checked to make sure that all the nodes which were given coordinates, have been used in the mesh. If an element is to be deleted, the pre-processor will automatically delete any nodes which thus become unused, and re-number the remaining nodes.

As an alternative to the above, a simple mesh can be defined graphically. If you choose this option, you first enter the extent of the mesh in the $x$ and $y$ directions. A graphics screen with this region is then displayed, and you click on it to define the node positions. When the nodes have been defined, you define the elements. For each element you define in a popup window its type (by {\tt NSIDE} and {\tt NNODE}) and material type; then click on the nodes around the element anticlockwise, in the graphics screen. Note that you must start at a corner node and proceed anticlockwise, as shown in the chart at the end of this chapter. This method can be used for beam, finite and infinite elements, but not for joint elements. Curved sides can be defined. If you subsequently wish to make changes, this can be done as for the keyboard-entry method.

Once a consistent mesh has been accepted, it will be displayed in the Graphics window, and you can now examine the mesh by clicking the mouse on the menu bar at the right-hand-side of the display. Here, there are options to display the element numbers and node numbers on the mesh, or to pick a particular node or element by the mouse or by typing in its number. There is also a zoom facility; click on Zoom' and then specify a rectangular window by clicking on two opposite corners of the window on the display. Note that if the zoom area is too small, and does not have any elements fully within it, a blank graph will result! An error message in this case will advise you to restore the full mesh and try a larger zoom area. When you are satisfied with the mesh, click on Continue' in the bottom right-hand-corner of the screen.

The next job is to specify the nodes at which fixed or Dirichlet boundary conditions apply. This is done graphically, by clicking the mouse on the appropriate nodes, for the conditions specified in the yellow pop-up box at the top of the screen. Note that for elasticity problems, there are three types of fixity which can be applied at a node:

1. fixed in the x-direction, but free to move in the y-direction;
2. free to move in the y-direction, but fixed in the x-direction;
3. fixed in both directions.
The three types of fixity are recorded by a two-digit fixity code: 10, 01 and 11 respectively. (A 1 in the i'th digit indicates that the displacement in the i'th coordinate direction is fixed). In a beam analysis (NDOFN=3) you will also be asked whether to fix the rotation; the displacement and rotation fixities are combined in a three-digit fixity code where the third digit refers to the rotation.

For Poisson problems, the Neumann or reflective boundary condition (no outward normal flow on the boundary) is the natural boundary condition for the problem, and parts of the boundary with this condition need no special treatment. Nodes with Dirichlet boundary conditions (specified values of the variable U) should be clicked in this section; they will be given a fixity code of 1. The actual value of U which is taken at the node, will be determined by defining values on Dirichlet planes in the "loading" section at the end; alternatively, a function g(x,y) in the main engine' POISS.FOR may be used to fix these at run-time.

The whole mesh is displayed again, but this time the parts of the boundary between fixed nodes are displayed in grey or black --- grey if fixed in one direction only, black if fixed in both. (Note that this does not strictly apply when there are more than two d.o.f.s per node!)

Next, the material property sets are entered, in the Text window. You may enter sets of properties, each corresponding to a different material; in the next data block you will associate these material property sets with the elements in the mesh. Note that set 1 will be the default material property set; any elements not associated with a specific set will be given the properties of set 1. If you have defined more than one property set, you will, for each set in turn, be prompted to click on the elements associated with that set. (To choose an element, click the mouse on the centre-of-gravity of the element, which is highlighted by a white dot).

In Basic-level analyses, there are two, three or four components in each property set;

• Young's modulus E, Moment of Inertia I and cross-sectional area A for beams;
• the elastic parameters E and \nu, and material thickness T for elasticity; you can also define a tensile strength \sigma_{ten} as component 4;
• and the coefficients a_x and a_y for Poisson's equation.
These will be stored as components 1,2 and 3 respectively. In Advanced-level analyses you are prompted for values of nine material property components; the additional components can be used for nonlinear analyses, for example in elasto-plasticity analyses the strength and coefficient of friction can be input (Note that the plasticity programs are for plane strain, and do not require input of material thickness --- the plasticity parameters are input to property components 3 onwards.) Details of the material constants used in the Advanced-level programs supplied with FELIPE, are given in Chapter 7.

The number of property sets you describe, will be written as an integer parameter NPROS, along with NELEM (no. of elements), NPOIN (no. of nodes) and NDOFN (no. of degrees of freedom per node), as the basic mesh details in the first data line after the title.

Once this has been done, the basic data set is complete. A print file (with filename test1.prt if the data file is test1.dat) will also be written, echoing the mesh details in a readable format (read using Notepad or a word processor). You can now exit the pre-processor, saving this data file, or you can modify it by for example refining the mesh and adding loading, to produce a new data file. You can also choose to add loads directly to the mesh you have created. When you modify an existing data file, you must give the new file a different name from that of the existing file.

## Modifying a data file

If you have chosen to modify the file you have created (or if you have answered n to the Create a new file?' prompt and selected a .dat file from the directory), the existing mesh will be displayed, and by clicking on Continue' you will see the Main Menu of modification options. These are:
1. Modify mesh
2. Double mesh
3. Refine mesh
4. Triangulate
5. Save mesh
6. Re-start from saved mesh
The Text window will also disappear at this point, and all subsequent interaction with the user will be through dialogue windows.

You choose an option by left-clicking on it. If you right-click on a menu option, a Help screen appears describing the option. The facilities in each of these options will now be described.

## Modify mesh

If you are modifying a datafile written at Basic level, you will be asked if you want to write the new file at Advanced level. If you do want to add Advanced features to the data, you should make this change first. In this case, there are two new facilities offered:
• you will be invited to specify extra mesh parameters (up to four reals and four integers) which will be written with the basic mesh data (NELEM, NPOIN, NDOFN, NPROS) at the start of the datafile. Your main engine' may for example require the value of the gravitational constant or the pore water bulk modulus, which can be written here.
• you will be asked if you want to add additional degrees of freedom. There may be a maximum of four degrees of freedom (d.o.f.s) per node. You can also have d.o.f.s defined only at corner nodes of quadratic elements. For example, the CONSL.FOR program models the saturated soil using 8-noded quadrilaterals with (u,v) displacement d.of.s at the four midside nodes and (u,v,p) d.o.f.s at the corners. We would want to increase to 3 the number of d.o.f.s at corner nodes of soil elements only; to avoid applying this to nodes of structural elements, you are prompted to say which element material type LMTYPE the new d.o.f.s should be created for.
You will then be offered the following menu:
1. Delete existing elements
3. Change nodal coordinates
4. Alter constraints ( i.e. nodal boundary conditions.) You can also change the number of degrees of freedom at an individual node, here.
5. Continue
These options work graphically: you choose an option, then click on the node or element to be processed. In most of the options, you are able to zoom in on the part of the mesh you wish to modify.

To add a new element, you specify its type by NSIDE and NNODE, then list the nodes around it anticlockwise. You will be prompted for the coordinates and fixity code of any new nodes. If they are midside nodes of the element their coordinates will be calculated by linear interpolation from the corner nodes of the side if you give their coordinates as 0,0, thus making the side straight.

When deleting elements, note that a maximum of 20 elements can be deleted in one batch, but there is a further limitation that not more than 50 nodes can thus become redundant and be deleted also. It is therefore best to delete only a few elements at a time. When nodes are no longer used in the mesh, they are deleted and the remaining nodes automatically renumbered. You can choose an individual element to delete, or you can delete a group of elements by clicking the right mouse-button first, then by clicking on the vertices of a polygon enclosing the group, as described in the next paragraph.

If you have created new degrees of freedom in the mesh, they will not have any boundary conditions associated with them. You can define Dirichlet nodes using alter constraints'. For example, if the d.o.f.s have been increased from (u,v) to u,v,p), clicking on the appropriate node will show a three-digit fixity code; the first two digits give the fixities of the original u,v d.o.f.s, and you can change the third digit from 0 to 1 to make the node lie on a Dirichlet boundary for pore-pressure or temperature.

When moving a node (changing its coordinates), you can click on the new position on the graph (zoomed if necessary). The coordinates of the point you have clicked are displayed, and you can edit these to fix the final position of the node. Note that if you wish to move a node beyond the window limits in a zoomed display, you will need to move it to the edge, then restore the full mesh, and zoom again if needed for the further move.

Clicking on Continue', you are given the chance to alter the values in the material property sets, including adding new sets (the Text window will reappear for this purpose), and then to change the assignment of elements to property sets, graphically. In this last task, you can click on the individual elements to be reassigned, or select a group of elements by enclosing their centres of gravity in a polygon; this is done by first right-clicking anywhere in the mesh, then left-clicking to define each of the vertices around the polygon. Right-click anywhere on the mesh again to join the last vertex to the first, completing the polygon. The elements whose centres lie within the polygon will be reassigned to the specified property set. This facility for defining groups of elements is used again later, in body force loading, for example.

The following Main Menu items allow the user, with a few mouse-clicks and keystrokes, to generate a complex mesh with hundreds of elements and nodes from his/her coarse mesh. Note that when the mesh is refined or doubled, the element material types and fixities get propagated to the newly-created nodes and elements; extra degrees of freedom also get assigned according to the rules that have been specified.

## Double mesh

Doubling a mesh produces a new mesh with twice the number of elements. The new half of the new mesh is created either by reflecting the existing mesh in a horizontal or vertical line (doubling by reflection) or by joining a copy of the mesh onto one side (doubling by translation). You need to select the type of doubling, and then to select the direction: for example, choosing the positive x-direction will produce a new half of the mesh joined to the existing half on its right-hand side. Node and element numbering in the pre-existing half of the mesh are unchanged by this process.

## Refine mesh

You can choose to refine the mesh Automatically or Manually. Under Automatic refinement, each quadrilateral element will be divided into 4 quadrilaterals, and each triangular element will be divided into 3 triangles. The new elements and nodes are numbered, but will need to be re-numbered before use in order to reduce the bandwidth of the global stiffness matrix -- see below.

Manual refinement only works for quadrilateral elements. Thus, to prepare a complex mesh for Poisson problems, it is most convenient to create the mesh from linear quadrilaterals, refine it using this option, and finally to subdivide each quadrilateral into two linear triangles using the Triangulate option. You have the opportunity to zoom to the area of interest before applying refinement; note that the refinement will still be propagated throughout the mesh as appropriate, not just in the zoom region.

Manual division of a quadrilateral element is performed by first specifying an edge by clicking on its corner nodes. Then dialogue windows will appear asking you the number of divisions, whether you want equal divisions, and -- if you want unequal divisions -- asking you the relative weights (proportional lengths) for each division. Be careful that you are entering the weights in the correct order, from the first end-node to the second. The weightings are relative only; typing 1,2,4,8 or 0.5,1.0,2.0,4.0 will have the same effect.

These divisions will then be used to subdivide the element, and this subdivision will be continue into neighbouring elements, and then throughout the mesh to ensure its consistency. Allocation of fixity codes for the new nodes, and property sets for the new subelements, is performed automatically.

Now a new edge can be chosen and subdivided. Note that there is a restriction on the number of individual element edges in a mesh in this option, so don't try anything too fancy! Also, it is better to subdivide a coarse mesh and then double the refined mesh, rather than doubling first and then refining, to avoid the risk of hitting this maximum on the number of edges.

Although this manual subdivision routine is restricted to quadrilateral elements, it is able to cope with linear and quadratic sides appearing in the same mesh, and will also subdivide one-dimensional and infinite elements (although the infinite' sides of infinite elements cannot be subdivided).

## Triangulate mesh

With this option, each quadrilateral element will be cut into two triangles --- a linear quadrilateral will produce linear triangles, and an 8-noded quadrilateral will produce quadratic triangles (with a new node created at the centrepoint). At the Basic level, the whole mesh will be triangulated, with no input from the user. At Advanced level, you can triangulate the whole mesh, a single element, or a group for elements (by defining a region around them using the right mouse-button, as described above).

## Save mesh, and Restart from saved mesh

The user should periodically save the current mesh to the output file, by clicking on the Save Mesh option. This will overwrite any previously-saved mesh. Should you then make a mistake in, say, refining the mesh, you can return to the main menu and click the Restart from Save option. This reads back the last-saved mesh, and you can re-start from this point.

## Reordering the elements and nodes

When the modifications of the previous section have been completed, the new mesh will look to be ready for use in the finite element analysis. However, it will be very inefficient, since the numbering of the nodes and/or elements will create a global stiffness matrix with a very large bandwidth or frontwidth. (For direct solvers such as Choleski decomposition, the bandwidth depends on the node numbering; for a frontal solver the frontwidth depends on the element numbering.)

You should therefore answer y to the prompt to re-order the elements. The reordering algorithm (a simple Cuthill-McKee algorithm) works by asking you to pick (with the mouse) an element at which to start the renumbering; this becomes element 1, then the neighbouring elements are numbered 2,3,... , and then their neighbours are renumbered, and so on through the mesh. The fronwidth of the resulting renumbered mesh is calculated, and if this is smaller than the frontwidth of the existing mesh, this renumbering is retained. You then have the chance to pick another element at which to start the renumbering, and after each such trial it is the renumbering producing the smallest frontwidth so far which is retained. When you are satisfied that you have achieved a well-numbered mesh, click on End'.

This has renumbered the elements, but not the nodes. There is now the option of renumbering the nodes as well; first the nodes around the new element 1 are numbered 1,2,3,... , then the as-yet-unrenumbered nodes around element 2, and so on. The consequent reduction in bandwidth is reported, and there is an option of displaying graphically the structure of the global stiffness matrix before and after renumbering. This is an extremely helpful feature for students needing to understand the relationship between node-numbering and stiffness matrix structure.

For elasticity problems, the loading r\'egime has to be added, and you will be prompted whether you wish to proceed to this, or to save the new file without adding loading data. For Laplace and Poisson problems, you will be invited to add data defining the values of u on the Dirichlet part of the mesh boundary.

This and the next section describe the Part B data which is now required to produce a complete data file for input to a main engine’.

In Basic-level elasticity and frame analyses, there are two types of loads which can be applied: point loads at nodes, and surface tractions on element edges. In Advanced-level elasticity analyses there are two further types of loads available: specified displacements at nodes, and body forces over elements. As in the Modify mesh' option, the user is able to zoom on the part of the mesh to be loaded, before applying each type of load.

To apply a point load at a node, click on the node and then type in the components of the force in the x- and y-directions. (A negative value of the component means that it acts in the negative x or y direction) The same process is used to apply specified displacements at nodes, at Advanced level.

At Advanced level, body forces can now be input. The components of the force are typed in, and then the user clicks on the elements to which this force applies, or defines a group of elements by clicking on the vertices of a bounding polygon as described earlier.

Surface tractions along element edges are entered by first creating surface traction sets defining the normal and tangential the components of the pressure acting at end nodes of the edge; each such pressure is allocated a traction set number. The first component of the pressure is the component acting normal to the edge; a positive value means that it acts onto the element, compressing it. The second component is the tangential component of the pressure, i.e. acting along the edge; a positive sign means that it acts anticlockwise around the element. Note that these are spot values of the pressure field (typically in Kilonewtons per square metre) distributed over the edge. In the case of line elements (representing trusses or beams), the body' of the element is imagined as if the node numbering continued, anti-clockwise, around a two-dimensional element with the line as the first edge.

Once the surface traction sets have been entered, they are assigned to nodes defining the edges on which those sets act. First click on an element with loaded side. Then click on the two corner nodes defining the edge to be loaded. Then type the surface traction sets corresponding to each corner node. If the same set is assigned to each corner node, it means that a uniform surface traction with this value is applied along the edge. If two different sets are assigned, the surface tractions change linearly from one value to the other along the side. If you have only defined one traction set, this will automatically be assigned as a uniform load on each loaded edge.

The next block of data, defining various types of boundary loading', is used for Poisson analysis and advanced elasticity problems. For Basic-level elasticity and beam analyses, a header for the variables nbdry (no. of boundary data sets) and nbdir (no. of Dircihlet boundary planes), followed by a blank line (i.e. nbdry = nbdir = 0) is written in the data file.

The second set of boundary loading' data defines planes on which the unknown potential field takes a specified value, i.e. Dirichlet boundaries. This is important in Poisson-type problems (NDOFN=1), and is also used for the scalar temperature field in the coupled elasticity program THERM.FOR (NDOFN=-2).

We first deal with radiating boundaries, for Poisson problems (NDOFN=1). Here, the boundary condition is

$$q_n = \bar{q} + \alpha U$$

where q_n is the normal flow (gradient of U normal to the boundary), \bar{q} is a specified flow, and \alpha is a radiation constant. Up to nine sets of radiation boundaries can be entered. The values of \bar{q} and \alpha are entered, and then you click on the nodes in the mesh associated with this boundary condition set.

For Advanced elasticity problems you will be asked if you wish to define excavation boundaries; these are only used in conjunction with program ELADV.FOR (or your own main engine' if you have included this feature in it). First, an in situ stress field can be prescribed. The type of stress field is decided by ISTYP, read from the first dataline in this section. The two parameters defining the stress field are also read from this line, and are:

 ISTYP = Stress field type (1 or 2) ISET = Set no. REAL1 = horizontal stress \sigma_x (ISTYP=1) or soil density \gamma (ISTYP=2) REAL2 = vertical stress \sigma_y (ISTYP=1) or lateral stress ratio K_0 (ISTYP=2)

Thus, for ISTYP=1 the stress at all Gauss points of soil elements (with material type LMTYP(L)>1) is (\sigma_x , \sigma_y). For ISTYP=2, the stress state at a soil element Gauss point with coordinates (x,y), y<0 is (-K_0 \gamma y , -\gamma y). Note that in ELADV, material type LMTYPE=1 is assumed to mean the element is a structural element with no initial stresses; soil/rock elements have LMTYPE=2,3,.... Also, with an overburden stress field, the plane y=0 is taken as the ground level. See the description of the program ELADV.FOR, for further information.

Finally, for Poisson problems or thermoelasticity analyses, you will be invited to input data defining the values to be assigned to the variable U on the Dirichlet part of the mesh boundary, which was specified by the fixity codes. (Parts of the boundary which have not been given such fixed values, will have a natural Neumann boundary condition -- e.g. insulated, in a heatflow application). This data section is also offered when preparing vector-variable meshes at Advanced level, since you may have created extra degrees of freedom for which scalar-variable-type boundary conditions are needed (the temperature field in thermoelasticity, for example). At Advanced level, you are therefore asked which degree of freedom each boundary condition set applies to. Note that this data is only handled in POISS.FOR and THERM.FOR.

You specify Dirichlet boundary values by defining a series of planes, horizontal or vertical, and for each plane assigning a value of u which all Dirichlet nodes which lie on that plane, will take. A vertical plane is defined by specifying the x-coordinate; a horizontal plane is defined by specifying the y-coordinate. In the text screen, you are prompted to type x or y, then the coordinate value, and finally the value of u associated with that plane. The program POISS.FOR is written to accept values for up to 8 such planes. The conditions will be applied in order, so that if a point lies on the intersection of two planes, it is the later condition which will apply.

For any Dirichlet nodes whose value has not been fixed by this method, the program POISS.FOR will use its in-built boundary value function bdryfn(x,y) to assign values at the Dirichlet nodes. The fixed-planes and the boundary-function techniques can thus be used in combination, to produce a complex boundary r\'{e}gime.

## Editing existing Part B data

If the input datafile already contains Part B data, the user may opt to edit this rather than create completely new Part B data. If the editing option is chosen, a warning is given that this will only make sense if you have not changed the mesh topology!’’. That is, in processing the Part A data earlier, any deletion of elements/nodes, or mesh doubling or refinement, will have resulted in renumbering of the elements/nodes, so that the load data in the input file will refer to elements/nodes whose numbers have changed.

Each applied load will be presented separately to the user, who can choose whether to include it in the new data file or to ignore it. There is also scope to edit it, by changing the load magnitude, for example, before including; additional new loads can also be specified. The locations of the elements or nodes referred to in each data item will be displayed on the mesh graph, so the user can see if the item still applies to the correct location. In the display of elements associated with a particular body force set, displayed elements can be deleted from the set by clicking on them, and new ones added. Part B data items which no longer make sense --- element/node numbers which now exceed NELEM/NPOIN, or surface tractions on elements where the node numbers no longer define an element edge --- will automatically be discarded. This load modification process also works with any boundary set or Dirichlet plane data from the original datafile.

## Exiting the pre-processor

At Advanced level, you are also invited to provide data specifying incremental loading and timestepping regimes - see below.

Once all the Part B data have been added, the pre-processor reports that the new data file has been written. The Part B data in a readable format will be appended to the print file (with a .prt filename extension) which contains the input data in readable form. You can also generate a PostScript graphics file of the mesh, which will have a .ps filename extension. Then you can exit the pre-processor, or return to modifying this or another data file. If you select a data-file to modify (usually the one which has just been written), the Text window will reappear.

## Format of the final data file

In the .dat data file produced by the pre-processor, there are one-line headers before each block of data --- these are ignored by the main program. The data file has the following block structure:

1: Title
2: Basic mesh constants
3: Element data: element-node connections, and property set
4: Nodal data: fixity code and coordinates
5: Material properties

The format of each block is as follows:
 Block Format Variables 1 A50 TITLE (max. 50 characters) 2 4I5 NPOIN NELEM NDOFN NPROS 3 I1, I4, I2, I3, 9I5 LCARD NSIDE LMTYPE LSET L (LNODS(L,N),N=1,NNODE) 4 I1,I4,I5,2G10.3 LCARD KODE NP COORD(NP,1) COORD(NP,2) NDOF(NP) 5 I1,I4,I5,10X,G10.3 LCARD NSET ICOMP VALUE
where

NPOIN = No. of nodes in mesh
NELEM = No. of elements in mesh
NDOFN = No. of degrees of freedom per node, see \S2.3
NPROS = No. of material property sets
NSIDE = No. of sides in element L
LMTYPE = Material type of element L
NNODE = No. of nodes in element L
LSET = Material property set number of element L
LNODS(L,N) = Global node number of the N'th node of element L
KODE = Fixity code of node NP; see above
COORD(NP,I) = I'th coordinate of node NP
NDOF(NP) = No. of degrees of freedom of node NP
VALUE = value of component ICOMP of material property set NSET

LCARD is a character which indicates the last line of a data block by holding the value 1; otherwise it is zero or blank.

At Advanced level, extra mesh constants (two reals, four integers) will be appended to block 2, so this line will have format 4I5,2G10.3,4I5.

Note that all variable names follow the Fortran naming convention: integers begin with a character in the range I--N.

For vector-variable problems, the mesh data is followed by data lines describing the applied loading. This data contains the following blocks (separated by header lines):

 Block Data type L1 Number of body force and surface traction sets L2 Nodal forces and specified displacements L3 Body force sets L4 Elements in each body force set L5 Surface traction sets L6 Edges with applied surface tractions

The format of each block is as follows:
 Block Format Variables L1 2I5 NBSET  NPSET L2 I1,I4,I5,2E10.3 LCARD LD NP FX FY L3 I1,I4,3X,I2,10X,E10.3 LCARD ISET I12 BXY L4 I1,I4,14I5 LCARD ISET elements with body force set ISET L5 I1,I4,3X,2I2,8X,2E10.3 LCARD NSET I1 I2 PN PT L6 I1,I4,I5,2(I5,I2,1X) LCARD NNAS L N1 NSET1 N2 NSET2

where:

NBSET = No. of body force sets
NPSET = No. of surface traction sets
LD = 1 for nodal force; 2 for specified displacement
FX,FY = Forces in x and y directions at node NP
BXY = Body force in x (if I12=1) or y (if I12=2) direction in body force set ISET
I1,I2 not used
PN,PT = Normal & tangential components of surface tractions in set NSET
N1,N2 = corner nodes of element L defining a loaded edge
NNAS = Number of nodes (2 or 3) along edge of element L
NSET1,NSET2 = surface traction sets associated with nodes N1 and N2

Blocks L3 and L4 do not appear in Basic-level data sets, and otherwise only appear if NBSET>0. Similarly, blocks L5 and L6 only appear if NPSET>0.

In block L5, the variables I1 and I2 take the values 1 and 2 respectively, to indicate that the following two values are the normal and tangential components of the traction, respectively.

Note: Earlier versions of FELIPE, up to and including version 3, used the format I1,I4,I5,3(I2,I3) for block L6. This has been changed in version 3.1 to allow for 4- and 5-digit node numbers in large meshes. This unfortunately means that data sets prepared under version 3, which involve surface tractions, must be modified manually to work with version 3.1.

For Poisson problems (and for vector-variable files written at Advanced level, i.e. for files with NDOFN less than 2), this data is followed by data lines describing the boundary conditions for scalar d.o.f.s. A header is followed by a line (format 2I5) stating the number of radiating (for NDOFN=1) or excavation (NDOFN=-2) boundary sets -- parameter NBDRY and the number of Dirichlet boundary planes -- parameter NBDIR. As with NBSET and NPSET above, if either of these parameters is zero then the corresponding block described below does not appear.

After a header line, radiating boundary sets (if NDOFN=1 and NBDRY is not 0) are defined in lines of the form:

LCARD,INDOF,ISET,ALPHA,QBAR, nodes in this set

in format I1,2I2,2E10.3,10I5, where

 INDOF = Degree of freedom to which this condition applies (default 1) ISET = Set no. ALPHA = radiation constant QBAR = specified normal flow

In the case of excavation boundaries (NDOFN=-2 and NBDRY not 0), the post-header lines are of the form:

LCARD,ISTYP,ISET,SX,SY, nodes in this set (if ISTYP=1: uniform stress field)

or

LCARD, ISTYP,ISET,GAMMA,K0, nodes in this set (if ISTYP=2: overburden stress field)

in format I1,2I2,2E10.3,10I5. See \S 2.6 for the definition of these variables.

The second block in this section, if NBDIR is not 0, defines Dirichlet boundary values (see \S 2.6 above). The data for the value of u on each plane defined by an x (for vertical plane) or y (for horizontal plane) coordinate, is given in a line:

LCARD,component, =',COORD,VALUE,INDOF

where component is either x' or y', in format: I1,A2,A2,E10.3,5X,E10.3,I5. INDOF is the degree of freedom to which the condition applies. For example, the line:

0 x =      4.000         60.000    1

means that U = 60.0 for Dirichlet nodes on the plane x = 4.0, where U is the 1st degree of freedom at the node.

Note that in all FELIPE programs, the standard Fortran naming convention (integers start with characters between I and N) is followed.

Nonlinear analyses normally involve adding the loading in increments. Thus once the total loads have been specified, the user is asked at Advanced level if these should be added in increments. The user specifies the number of increments (maximum 9), and then the size of each increment; this data is appended to the data file in a line of format I5,9F5.2 .

Parabolic p.d.e.s, and time-dependent material models such as elasto-viscoplasticity and poroelasticity, require a timestepping algorithm in their main engines', and the data for this can be input as the final section of load data at Advanced level. The parameters which can be specified are:

1. {Initial timestep size DTINT - default 0.001}
2. {Timestep increase factor FTIME - default 1.2}
3. {Max. allowable timestep increase ALLOW - default 10.0}
4. {Time discretization parameter THETA - default 0.0}
5. {Final time TTEND - default 100.0}
6. {Timestep for initial solution DTONE - default 0.0}
The purpose of these parameters is to define the following timestep control algorithms:

## Algorithm used in VPLAS.FOR:

The timestepping will start with a timestep of \Delta t_0=DTINT, the next timestep will be \Delta t_1=DTINT*FTIME, and thereafter the timestep will be increased by a factor of FTIME at each iteration (unless this is overridden by a timestep control algorithm within the main engine' itself). When the timestep would be greater than DTINT*ALLOW, it is capped at \Delta t=ALLOW. Whenever a new load increment is added, the timestep reverts to \Delta t=DTINT.

The parameter THETA, 0 <= \theta <= 1, defining the time discretization algorithm is a common feature of many schemes; \theta=0.0 defines an explicit discretization, \theta=1.0 defines a fully-implicit scheme.

TTEND and DTONE are not used in the viscoplasticity timestepping r\'egime.

## Algorithm used in CONSL.FOR

For soil consolidation, an initial solution is obtained using the timestep DTONE - for the two-level element in the CONSL.FOR program, a value of zero is acceptable for this. Then, the timestep DTINT is used for five timesteps before it is increased by a factor of FTIME, for a further five steps, and so on. Timestepping halts when the time TTEND is reached, or earlier if requested by the user at run-time. THETA is also used in this program; a fully-implicit scheme where THETA=1, is recommended.

The timestepping data DTINT,FTIME,ALLOW,THETA,TTEND,DTONE are written in a line with format 6E10.3.

Of course, you are free to modify the timestepping r\'egimes in the main engines' to your particular needs, and recompile the programs.