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.