ADDING NEW ELEMENTS IN CRISP
A GUIDE TO DEVELOPERS
Introduction:
Adding new elements in CRISP may seem difficult, but if you follow this guide you will find it a straightforward job. CRISP was developed in a way that allows easy addition of new elements.
You would need to have the source code of CRISPMP (the main program) as well as CRISPGP ( the geometry program). The Front Squasher (CRISPSQ) may not be required. To understand the function of these three modules of CRISP and how they interact with each other, please refer to SAGECRISP system description
If your new element resembles any of the existing elements in CRISP, ie having 4 vertex nodes as in a quadrilateral or 3 vertex nodes as in a traingle, then you may use SAGECRISP interface to generate the input data for the Front Squasher, run CRISPSQ then edit the resulting .gpd file to change the element type. Next you would need to run your newly modified geometry program to generate the remaining nodes (midside or internal) for your new element.
5 Steps for adding new element
Step 1 – Changes in BLOCK DATA
This fortran routine is present in all the three modules of CRISP. You would need to extend the following arrays to include your element details.
Array LIN(iparam,letp). This is also known as LINFO elsewhere in CRISP. LIN(iparam,letp) contains data for element type letp as follows
iparam 
Description 
Variable name 
1 
Total No. of NODES (displacements + pore pressure) 
NDPT 
2 
Total No. of VERTEX NODES 
NVN 
3 
Total No. of ELEMENT EDGES 
NEDG 
4 
Total No. of ELEMENT FACES (3D) 
NFAC 
5 
Total No. of DISPLACEMENT NODES 
NDN 
6 
Total No. of POREPRESSURE NODES 
NPN 
7 
No. of DISPLACEMENT NODES PER EDGE (execluding end nodes) 
NSDN 
8 
No. of POREPRESSURE NODES PER EDGE (execluding end nodes) 
NSDP 
9 
No. of INNER DISPLACEMENT NODES 
NIND 
10 
No. of INNER POREPRESSURE NODES 
NINP 
11 
No. of INTEGRATION POINTS 
NGP 
12 
Index to WEIGHTS AND INTEGRATION POINT coords 
INDX 
13 
NOT USED 

14 
Index to NODES ALONG EDGE (ARRAYS NP1, NP2) 
INDED 
15 
No. of AREA COORDINATES 
NL 
16 
Total No. of DEGREES OF FREEDOM (D.O.F.) IN ELEMENT 
MDFE 
17 
CENTROID INTEGRATION POINT NUMBER 
NCGP 
18 20 
NOT USED. 

21  50 are the No. of D.O.F. OF EACH NODE OF ELEMENT 
NDFN 
Subscripts 2150 assume that no element type has more than 30 nodes. Otherwise increase the limiting size for iparam to accommodate extra nodes.
All of these parameters can be defined easily. The exception is INDX and INDED. For information on how to define INDX see section on arrays W and L below. Defining INDED is described in step 4 below.
In CRISP release 97.1, the following elements are available
letp 
Element description 
Dimensions & d. of f 
Integration 
1 
3NODED BAR 
(2D Plane strain only) 
5 points 
2 
6NODED LST TRIANGLE 
(2D) 
7 points 
3 
6NODED LST TRIANGLE 
(2D Consolidation) 
7 points 
4 
8NODED QUADRILATERAL. 
(2D) 
3x3 
5 
8NODED QUADRILATERAL 
(2D Consolidation) 
3x3 
6 
15NODED CUST TRIANGLE 
(2D) 
16 
7 
22NODED CUST TRIANGLE 
(2D Consolidation) 
16 
8 
20NODED BRICK 
(3D) 
3x3x3 
9 
20NODED BRICK 
(3D Consolidation) 
3x3x3 
10 
10NODED TETRAHEDRA 
(3D) 
4 points 
11 
10NODED TETRAHEDRA 
(3D Consolidation) 
4 points 
12 
3NODED BEAM 
(2D Plane strain only) 
5 points 
13 
8NODED SLIP ELEMENT 
(2D Plane strain only) 
5 points 
14 
2NODED BAR 
(2D Plane strain only) 
5 points 
15 
2NODED BEAM 
(2D Plane strain only) 
5 points 
The following elements were also added for a future release (early 1998)
16 
8 NODED QUADRILATERAL 
(2 D Seepage) 
3x3 
17 
4 NODED QUAD 
(2 D Seepage) 
2x2 
18 
4 NODED QUAD 
(2 D) 
2x2 
In addition to array LIN, you would also need to extend array MIN(ivarb,inode,letp) to add your element details. This contains the Unique Program Variable Number (UPVN) for each of the variable of a node as follows:
UPVN 
Description 
Variable name 
1 
Displacement in X direction 
U 
2 
Displacement in Y direction. 
V 
3 
Displacement in Z direction 
W 
4 
Excess pore pressure 
P 
5 
Rotation in XY plane 
RXY 
6 
Not used 
As an example, for the 8noded consolidation element, considering the first corner node only, we have:
MIN(1,1,5)=1 for X displacement
MIN(2,1,5)=2 for Ydisplacement
MIN(3,1,5)=4 for pore pressure variable
Node 5, being a midside node with no pore pressure variable, would have:
MIN(1,5,5)=1 for X displacement
MIN(2,5,5)=2 for Ydisplacement
Finally, it may be necessary to extend arrays W(indx) and L(idim,indx) to include the integration point details for your new element if these are not available in the current arrays W and L.
Array W(indx) contains the weights of the integration points. The position of the new weights is associated with the index indx defined in array LIN above. So for example, the 8noded consolidation element has indx=12 as defined in LIN. Therefore, its integration points must start from 13. These points would fill up the positions W(13..to..21)
Similarly, array L(idim,indx) may have to be extended accordingly. This contains the integration point coordinates, bearing in mind that idim for triangles refers to area coordinates
A final note concerns the sizes of LIN, and MIN, which would have to be adjusted throughout CRISPMP and CRISPGP. The maximum number of element types is defined in the geometry program CRISPGP by the variable LTZ (set as =16 in release 97.1 of CRISP). This would have to be adjusted too in order to accommodate the new element.
Step 2 – Changes to conditional statements on element type
You would now need to check the program for IF conditions related to element types. The element type is referred to by the variable LT (would be called LETP in new release). Search for all of these IF conditions containing LT and adjust accordingly.
Step 3 – Changes in LSTIFF
This concerns element details specific for elements with pore pressure nodes. The data for these elements are in subroutine LSTIFF in the main program. The aim here is to define two arrays KP and KD which contain indexes which would help in the assembly of the element stiffness. If your element does not contain pore pressure variables, then you needn't worry about this step.
Array KP is associated with array NXP(letp) which is also defined in LSTIFF. Similarly, array KD is associated with NXD(letp).
NXP(letp) contains indexes for each element type which are associated with the starting point in the corresponding array KP where the variable locations are defined.
For example, with the 8noded consolidation element (type 5) NXP(5)=3. Therefore the starting point for array KP is 3+1 giving:
KP(4)=3, KP(5)=6, KP(6)=9, KP(7)=12
where 3,6,9 and 12 are the locations of pore pressure variables for this consolidation element
Similarly, for the 8noded consolidation element, NXD(5)=41. The starting point for array KD is 41+1 giving:
KD(42)=1, KD(43)=4, KD(44)=7, KD(45)=10,
KD(46)=13, KD(47)=15, KD(48)=17, KD(49)=19
where 1,4,7,10,13,15,17 and 19 are the locations of the first displacement variable (ie xdirection) for this consolidation element. The locations for Y and Z displacement may follow sequentially by making use of NDIM.
If your new element contains pore pressure nodes, then you may need to extend arrays KP,KD,NXP, and NXD to include your variable locations.
Step 4 – Defining INDED
This concerns the definition of INDED in array LIN in BLOCK DATA. INDED is an index to nodes along the edge of an element. It is used by routine SETNP in the main program in order to help define arrays NP1 and NP2. Arrays NP1(iedge) and NP2(iedge) ,in turn, give the index to the element nodal connectivity array NCORR for nodes at either end of each element edge. Arrays NP1 and NP2 are defined in routine SETNP as follows:
SUBSCRIPT 
1 
2 
3 
4 
5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
NP1 
1 
2 
3 
1 
2 
3 
4 
5 
6 
7 
8 
1 
2 
3 
4 
1 
2 
3 
1 
2 
3 
NP2 
2 
3 
1 
2 
3 
4 
1 
6 
7 
8 
5 
5 
6 
7 
8 
2 
3 
1 
4 
4 
4 
For any quadrilateral element, the vertex nodes are numbered as 1,2 3 and 4. This leads to each of the four edges being numbered as 12,23,34 and 41. These can be extracted from the above table giving NP1(4 to 7)=1,2,3,4 and NP2(4 to 7)=2,3,4,1. The following loop shows how these numbers can be obtained from arrays NP1 and NP2:
DO 20 IEDGE = 1, NEDG ! NEDG = 4 For a quad element
JNODE1 = NP1(IEDGE+INDED)
JNODE2 = NP2(IEDGE+INDED)
By setting INDED=3 for this quadrilateral element we obtain the vertex node numbers as follows
Loop on iedge 1 gives JNODE1=1, JNODE2=2
Loop on iedge 2 gives JNODE1=2, JNODE2=3
Loop on iedge 3 gives JNODE1=3, JNODE2=4
Loop on iedge 4 gives JNODE1=4, JNODE2=1
JNODE1 and JNODE2 may be used in the above loop to obtain the actual node numbers in the mesh using NCORR.
In CRISP, all triangular elements have the vertex nodes numbered as 1,2 and 3 followed by midside and internal nodes. This results in INDED=0. Similarly, for bar and beam elements, with the end nodes as 1 and 2, INDED=0. Quadrilateral, Brick and slip elements have INDED=3, while tetrahedra elements have INDED=15.
Step 5 – Changes in UPOUT
Finally, you would need to decide what changes you need to make in UPOUT to accommodate the new element. It is worth bearing in mind that array NQ(nnodes) stores the number of variables for each node in your mesh. NQ is set in routine MAKENZ and it extracts these degrees of freedom from array MIN which was set in BLOCK DATA. So it may be convenient to use NQ to define your output statements according to the number of variables you have for each node.