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 SAGE-CRISP 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 SAGE-CRISP 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 (mid-side 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 21-50 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 |
3-NODED BAR |
(2-D Plane strain only) |
5 points |
2 |
6-NODED LST TRIANGLE |
(2-D) |
7 points |
3 |
6-NODED LST TRIANGLE |
(2-D Consolidation) |
7 points |
4 |
8-NODED QUADRILATERAL. |
(2-D) |
3x3 |
5 |
8-NODED QUADRILATERAL |
(2-D Consolidation) |
3x3 |
6 |
15-NODED CUST TRIANGLE |
(2-D) |
16 |
7 |
22-NODED CUST TRIANGLE |
(2-D Consolidation) |
16 |
8 |
20-NODED BRICK |
(3-D) |
3x3x3 |
9 |
20-NODED BRICK |
(3-D Consolidation) |
3x3x3 |
10 |
10-NODED TETRA-HEDRA |
(3-D) |
4 points |
11 |
10-NODED TETRA-HEDRA |
(3-D Consolidation) |
4 points |
12 |
3-NODED BEAM |
(2-D Plane strain only) |
5 points |
13 |
8-NODED SLIP ELEMENT |
(2-D Plane strain only) |
5 points |
14 |
2-NODED BAR |
(2-D Plane strain only) |
5 points |
15 |
2-NODED BEAM |
(2-D 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 X-Y plane |
RXY |
6 |
Not used |
As an example, for the 8-noded consolidation element, considering the first corner node only, we have:
MIN(1,1,5)=1 for X displacement
MIN(2,1,5)=2 for Y-displacement
MIN(3,1,5)=4 for pore pressure variable
Node 5, being a mid-side node with no pore pressure variable, would have:
MIN(1,5,5)=1 for X displacement
MIN(2,5,5)=2 for Y-displacement
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 8-noded 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 8-noded 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 8-noded 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 x-direction) 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 1-2,2-3,3-4 and 4-1. 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 mid-side 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.