Isogeometric Analysis 2011

Austin, Texas – January 13-15, 2011

Overview on Use of FEAP

www.ce.berkeley.edu/˜rlt/presentations/

Robert L. Taylor
Department of Civil & Environmental Engineering
University of California, Berkeley

12 January 2011
Solving Problems using FEAP

• This lecture presents:

– Summary of capabilities of FEAP.

– Input data file structure.

– Problem solution modes.

– Graphics capability.

– Use with T-splines meshing.

Background Reading

• Additional information available in:

– On line manuals for FEAP at

www.ce.berkeley.edu/feap
Manuals available: Installation, User, Contact, Parallel, Example,
Programmer).

– Books on FEM: O.C. Zienkiewicz and R.L. Taylor, The Finite

Element Method, 6th edition, Elsevier Butterworth-Heinemann,
Oxford, 2005.
See last chapter in Basis and Solids volumes.

• FEAP distributed under license by UC.

• Small version FEAPpv available free: www.berkeley.edu/feap/feappv.

FEAP Program Capabilities

• FEAP: Finite Element Analysis Program.

• Solves problems formulated by a finite element method.

• Problems can be:

– Linear or non-linear.
– Static, quasi-static or transient.
– Coupled (multi-physics) homogeneous or partitioned.
– Parallel (using PETSc)
– Multi-scale F E 2 form.

• Solution in batch or interactive mode.

• Output: Print, graphics (screen or PostScript), time history files.

FEAP Program Capabilities

• Finite element library includes:

– Small and finite deformation analysis of solids.

– Thermal analysis of solids.

– Small and large displacement frame (bending, shear & axial

deformation).

– Small and large displacement membrane.

– Small displacement shell.

• Elements can be used in 1-d, 2-d & 3-d analyses.

FEAP Program Capabilities

• Small deformation elements use material library containing:

– Elastic (Isotropic & orthotropic with thermal expansion)

– Visco-elastic (isotropic deviatoric only; complex moduli)

– Elasto-plastic – J2 with isotropic and kinematic hardening

– Generalized plasticity – J2 model.

– RVE - interface (multi-scale use)

– User model interface.

FEAP Program Capabilities

• Finite deformation elements use material library containing:

– Elastic (Neohookean, Mooney-Rivlin, St.Venant-Kirchhoff,

Ogden, Fung, Arruda-Boyce, Yeoh)

– Visco-elastic (isotropic deviatoric with damage)

– Elasto-plastic – J2 with isotropic and kinematic hardening

– RVE - interface (multi-scale use)

– User model interface.

FEAP Program Capabilities

• Startup (Windows):
FEAP Program Capabilities

• After plot initiation (Windows):

FEAP Program Capabilities

• Startup & Graphics (Unix/Mac):

FEAP Input Data Files

• Description of mesh by input data file(s).

• Basic structure of file:

– Control data

– Mesh inputs (nodes, elements, material data, loads, etc.)

∗ Nodes: 1:NUMNP; Elements 1:NUMEL; Materials 1:NUMMAT

– Mesh manipulation (e.g., merge parts, link nodes, etc.)

– Contact interface description (none with NURBS or T-splines).

– Solution using command language statements.

FEAP Input Data Files
• Control data:

FEAP * * title for output file

NUMNP NUMEL NUMMAT NDIM NDOF NEN

where
– FEAP - Indicates start of problem.
– NUMNP - Number of nodes
– NUMEL - Number of elements
– NUMMAT - Number of material sets
– NDIM - Spatial dimension of mesh
– NDOF - Maximum number d.o.f/node
– NEN - Maximum number nodes/element

• Often, NUMNP, NUMEL, NUMMAT can be input as zero (0).

For T-spline solution: all can be omitted.
FEAP Input Data Files

• Example: Mesh for a patch test.

FEAP * * Patch test
0 0 0 2 2 4
BOUNdary
COORdinates 1 0 1 1
1 0 0.0 0.0 4 0 1 0
2 0 10.0 0.0
3 0 10.0 10.0 FORCe
4 0 0.0 10.0 2 0 100.0 0.0
5 0 2.0 2.0 3 0 100.0 0.0
6 0 7.5 3.0
7 0 6.5 7.0 MATErial 1
8 0 2.5 6.5 solid
elastic isotropic 200e9 0.3
ELEMents
1 1 1 1 2 6 5 END
4 0 1 4 1 5 8 INTEractive
5 0 1 5 6 7 8 STOP

• FEAP counts: 8 nodes, 5 elements & 1 material.

• Uses first 4-characters of text.

FEAP Input Data Files

• Above input statements give the mesh and after solution the plot
may be added.
4 3

3 DISPLACEMENT 1

7 0.00E+00
8 1.30E-10
2.60E-10
3.90E-10
5.20E-10
2 6.50E-10
4 5 7.80E-10
9.10E-10
Current View
6 Min = 0.00E+00
X = 0.00E+00
Y = 0.00E+00
5
Max = 9.10E-10
X = 1.00E+01
1 Y = 0.00E+00

Time = 0.00E+00
1 2

• Some basic solution commands described later.

FEAP Input Data Files

• Mesh data can be split into parts using INCLude option.

• Example: Patch test.

FEAP * * Patch test File Ipmesh contains

0 0 0 2 2 4
COORdinates
INCLude Ipmesh 1 0 0.0 0.0
2 0 10.0 0.0
MATErial 1 .... etc. ....
solid
elastic isotropic 200e9 0.3 FORCe
2 0 100.0 0.0
END 3 0 100.0 0.0
INTEractive
STOP etc.

• Similar form used later to interface T-spline file.

FEAP Input Data Files

• FEAP commands interpret 4 characters.

• FEAP numerical data given by:

– Numerical: 1, 1.0, -5.3e-3
– Parameters (limited to 2-characters):
aa, z1, e0
– Expressions: x0+r*sind(30), e1*b*h^3/12
Operations: +, -, *, /, ^
Functions: sin, abs, exp, atan, ...

• Parameters set by commands:

PARAmeter
x0 = 5.56
pi = acos(-1.0)
! End with blank line.

Note: Data following ! not used by FEAP.

FEAP Input Data Files

• All data read by parser.

• FEAP counts number of nodes, elements and material sets.

Can be zero on control statement.

• Standard data record structure:

– Unformatted input: Field widths limited to 15 characters
– No more than 16 data items per record.
– Data items separated by: comma (,); space ( ); equal (=).
– Blank characters ignored except in expressions.
Note: x + r*sind(3) would be read as 3 fields.
– Data after ! ignored (comment)
FEAP Input Data Files

• Nodal coordinate and element connection

– Coordinates by node number (COOR):
NUMBER N_GEN X_1 X_2 X_3
∗ NUMBER = node number
∗ N_GEN = increment to next node
∗ X_i = value (i = 1,NDIM)
– Elements by number (ELEM):
NUMBER N_GEN N_MAT IX_1 .. IX_NEL
∗ NUMBER = element number
∗ N_GEN = increment to nodes
∗ N_MAT = material set identifier
∗ IX_i = node number (i = 1,NEL; NEL ≤ NEN)
FEAP Input Data Files

• Nodal data specified by:

– Node number
FORCe
5 0 5.0 -2.3
sets node 5 force: F1 = 5, F2 = −2.3.
– Edge coordinate value
EBOUndary
2 5.0 1 0
sets BC code: ID1 = 1 and ID2 = 0 for nodes with x2 = 5.
Note: A zero has unknown solution value.
– Coordinate value
CFORce
node 5.0 0.0 5.0 -2.3
sets force at node closest to x1 = 5, x2 = 0 to F1 = 5 and
F2 = −2.3.
FEAP Input Data Files

• Material data sets

– Specifies element type: SOLId, PLATe, SHELl, FRAMe, TRUSs, GAP,

POINt, THERmal, USER.

– Defines associated element group.

– Defines degree-of-freedom assignments.

– Describes constitutive model property values.

– Defines finite element formulation (e.g., displacement, mixed,

small or finite deformation, etc.).

– Defines other element properties (e.g., quadrature order, body

loading, etc.).

– Specify NURBS or T-spline interpolation and quadratures.

FEAP Input Data Files

• Material set: Form for linear elastic T-spline solid element

MATErial ma
SOLId
ELAStic ISOTropic E nu
T-SPline interp q1 q2 q3 ! (or NURBS)

– Use standard displacement model to describe elements

where N MAT = MA.

– Elastic properties are isotropic with parameters set by E (Young’s

modulus) and nu (Poisson’s ratio).

– Specifies T-spline interpolation with q1,q2,q3 quadrature.

– Small deformation element.

FEAP Input Data Files

• Finite deformation solids

– Set explicitly:
MATErial ma
SOLId
ELAStic ISOTropic E nu
T-Spline quadratrue q q q
FINIte

– Set implicitly by model

MATErial ma
SOLId
ELAStic NEOHook E nu
T-Spline quadratrue q q q

• Small deformation set by SMALl.

NURBS Isogeometric Modeling
NURBS Analysis procedure in FEAP

• Define coarse set of control points, knots, 1-d knot-point list,

side-patch description:

• Example: Curved beam – specification of NURBS and knots

y
NURBs − Control points
1 0 5.0 0.0 1.00 6 5
2 0 5.0 5.0 sqrt(2)/2
3 0 0.0 5.0 1.00
4 0 10.0 0.0 1.00 b
3 2
5 0 10.0 10.0 sqrt(2)/2
6 0 0.0 10.0 1.00
a r

KNOTs θ x
1 4
knot 1 0.00 0.00 1.00 1.00 P
knot 2 0.00 0.00 0.00 1.00 1.00 1.00
NURBS Isogeometric Modeling

• Example: Continued – Specification of sides and patch

y
− Control points
NSIDes
6 5
side 1 0 2 1 2 3
side 2 0 2 4 5 6
side 3 0 1 1 4
b
3 2
NBLOck
block 2 1 3 a r
θ x
1 4

P
• Need to add material properties, loading and boundary conditions.
Use standard FEAP commands for most.
NURBS Isogeometric Modeling
k-refinement in circumferential direction

(a) Elevate radial knot (b) Insert knot 1

(c) Insert knot 2 (d) Insert knot 3

NURBS Isogeometric Modeling
• Elevate NURBS block 1, direction 3, 2 orders.
include Ispini
batch
elevate init
elevate block 1 3 2
elevate end
end
• Knot insertions for NURBS block 1 in direction 3.
parameter
kk = 0

loop,4
parameter
kk = kk + 0.25

include Ispini
batch
insert init
insert block 1 3 kk 1 ! Last entry is number of times
insert end
end
next
NURBS Isogeometric Modeling
k-refinement in FEAP (Cont.)
• Each elevation or insertion creates a flat mesh NURB mesh
• Can be used to describe problems recursively in input file using
INCLude NURB_mesh
• Use loops (in input file) to perform repeated insertions
LOOP,9
PARAmeter
d = d + 0.1

INCLude NURB_mesh
BATCh
INSErt INITialize
INSErt KNOT 1 d 1
INSErt END
END
NEXT

inserts knot 1 nine times at intervals of 0.1 units.

NURBS Isogeometric Modeling
• Curved beam under end shear results, exact energy
1
Eex = [log 2 − 06] = 0.02964966844238
2
π
10

0
10

−2
Energy Error

10

−4
10

Q4
Q9
−6
10 Q16
H−Nurb: p=1,2
H−Nurb: p=2,2
P−Nurb: h=1,2
−8
10 1 2 3 4
10 10 10 10
Number Equations

Curved beam subjected to end shear. Energy error.

NURBS Isogeometric Modeling
Confined tension strip

• Consider rectangular 8 × 8 with unit circular hole.

• Lateral boundaries restrained and strip stretched by

uniform displacement

• Material given by modified Neo-Hookean model

J −2/3 b :
   
W =1
4K J 2 − 1 − 2 log J +1
2G 1−3
where b = F FT and properties set to:

K = 1.6 × 108 and G = 3.2 × 104

Results in a nearly incompressible behavior.

• Model by Q9/P1 standard mixed elements and F̄ NURBS

NURBS Isogeometric Modeling
Confined tension strip: Q2/P1 Elements

_________________
DISPLACEMENT 1 _________________
DISPLACEMENT 2

0.00E+00 0.00E+00
3.15E-02 3.85E-02
6.30E-02 7.69E-02
9.44E-02 1.15E-01
1.26E-01 1.54E-01
1.57E-01 1.92E-01
1.89E-01 2.31E-01
2.20E-01 2.69E-01
2.52E-01 3.08E-01
2.83E-01 3.46E-01
3.15E-01 3.85E-01
3.46E-01 4.23E-01
3.78E-01 4.61E-01

Time = 2.00E+01 Time = 2.00E+01

(a) u-displacement (b) v-displacement

_________________
Mises Stress _________________
DISPLACEMENT 2

6.99E+01 0.00E+00
4.08E+03 3.85E-02
8.09E+03 7.69E-02
1.21E+04 1.15E-01
1.61E+04 1.54E-01
2.01E+04 1.92E-01
2.41E+04 2.31E-01
2.82E+04 2.69E-01
3.22E+04 3.08E-01
3.62E+04 3.46E-01
4.02E+04 3.85E-01
4.42E+04 4.23E-01
4.82E+04 4.61E-01

Time = 2.00E+01 Time = 2.00E+01

(c) Mises stress (d) Deformed shape

NURBS Isogeometric Modeling
Confined tension strip: Q2/Q1/Q1 u − p − θ NURBS

_________________
DISPLACEMENT 1 _________________
DISPLACEMENT 2

4.22E-15 0.00E+00
3.15E-02 3.84E-02
6.29E-02 7.69E-02
9.44E-02 1.15E-01
1.26E-01 1.54E-01
1.57E-01 1.92E-01
1.89E-01 2.31E-01
2.20E-01 2.69E-01
2.52E-01 3.08E-01
2.83E-01 3.46E-01
3.15E-01 3.84E-01
3.46E-01 4.23E-01
3.77E-01 4.61E-01

Time = 2.00E+01 Time = 2.00E+01

(a) u-displacement (b) v-displacement

_________________
Mises Stress _________________
STRESS 2

2.36E+02 1.79E+03
3.71E+03 5.73E+03
7.18E+03 9.68E+03
1.07E+04 1.36E+04
1.41E+04 1.76E+04
1.76E+04 2.15E+04
2.11E+04 2.54E+04
2.46E+04 2.94E+04
2.80E+04 3.33E+04
3.15E+04 3.73E+04
3.50E+04 4.12E+04
3.84E+04 4.52E+04
4.19E+04 4.91E+04

Time = 2.00E+01 Time = 2.00E+01

(c) Mises stress (d) Deformed shape

NURBS Isogeometric Modeling
Confined tension strip: Q2/Q1/P1 u − p − θ NURBS

_________________
DISPLACEMENT 1 _________________
DISPLACEMENT 2

4.12E-15 0.00E+00
3.15E-02 3.85E-02
6.29E-02 7.69E-02
9.44E-02 1.15E-01
1.26E-01 1.54E-01
1.57E-01 1.92E-01
1.89E-01 2.31E-01
2.20E-01 2.69E-01
2.52E-01 3.08E-01
2.83E-01 3.46E-01
3.15E-01 3.85E-01
3.46E-01 4.23E-01
3.78E-01 4.61E-01

Time = 2.00E+01 Time = 2.00E+01

(a) u-displacement (b) v-displacement

_________________
Mises Stress _________________
DISPLACEMENT 2

2.44E+02 0.00E+00
3.72E+03 3.85E-02
7.19E+03 7.69E-02
1.07E+04 1.15E-01
1.41E+04 1.54E-01
1.76E+04 1.92E-01
2.11E+04 2.31E-01
2.46E+04 2.69E-01
2.80E+04 3.08E-01
3.15E+04 3.46E-01
3.50E+04 3.85E-01
3.84E+04 4.23E-01
4.19E+04 4.61E-01

Time = 2.00E+01 Time = 2.00E+01

(c) Mises stress (d) Deformed shape

NURBS Isogeometric Modeling
Elastic-Plastic Necking

• Plane strain strip under uniform extension

• Finite strain elastic-plastic model with J2 yield function.

• Material properties are:

K = 164.206 ; G = 80.1938 ; σ0 = 0.45 ; σ∞ = 0.715 ; h = 0.12924 ;
Uniaxial yield stress given by s
2
σy = σ∞ + (σ0 − σ∞) exp(−βep) + h ep
3
where ep effective plastic strain and h isotropic hardening.

• Geometry: Length is 53.334 and width is 12.826. Symmetry used

for one quadrant model (reduce center to 0.982 of width).

• Use quadratic NURBS with 10 knots in width and 20 in length.

Analysis also performed with Q1/P0 & Q2/P1 elements.
NURBS Isogeometric Modeling
Elastic-Plastic Necking Example

_________________
Mises Stress

5.71E-01
5.95E-01
6.18E-01
6.42E-01
6.66E-01
6.89E-01
7.13E-01
7.37E-01
7.60E-01
7.84E-01
8.07E-01
8.31E-01
8.55E-01

Time = 5.00E+00

(a) Mesh (b) Mises stress at full extension

 NURBS Isogeometric Modeling
Elastic-Plastic Necking Example

3 5
FE − Q1
FE − Q1/P0 4.5
FE − Q2/P1
2.5
NURB − p = 2,2 4

3.5
Radial Displacement

1/2 − Reaction Sum

2
3

1.5 2.5

2
1
1.5

1 FE − Q1
0.5
FE − Q1/P0
0.5 FE − Q2/P1
NURB − p = 2,2
0 0
0 1 2 3 4 5 0 1 2 3 4 5
t − Time Y−Displacement

(c) Necking displacement (d) Load-displacement

FEAP Command Language Solution

• After FE mesh described, users responsible for solution steps.

• May be given in BATCh or INTEractive mode.

• Example of batch solution:
BATCh
TANGent ! Form "K" Tangent
FORM ! Form "R = F - P(u)" Residual
SOLVe ! Solve "K du = R"
END
or
BATCh
TANGent,,1 ! Means same as above
END
• Interactive solution with:
INTEractive
• Can give more than one batch and/or interactive set.
FEAP Command Language Solution
• LOOP-NEXT can be used in solution mode also.

• Useful for iterative or time stepping solutions.

• Time stepping for linear transient solution:
TRANSient ! Activate Newmark
DT,,dt ! Set time step
LOOP,,nt ! DO ’nt’ steps
TIME ! Advance time
TANGent,,1 ! Solve step
.... ! Outputs/plots
NEXT ! End of loop
...
• Iterative solution for non-linear problem
LOOP,,ni ! DO ’ni’ iterations
TANGent,,1 ! Solve iterate
NEXT ! End of loop
...
FEAP Command Language Solution
• Transient non-linear solution combines both:
TRANSient
DT,,dt
LOOP,,nt
TIME
LOOP,,ni
TANGent,,1
NEXT
....
NEXT
...
• Note: All command language statements have form:
TEXT TEXT Value Value Value
• Example: Solution step with line search:
TANGent,LINE_search,1
• FEAP can have unsymmetric tangent by:
UTANgent,,1
FEAP Command Language Solution
• Using FEAP in interactive mode.

• Permits easy view of solution at each command step.

• Graphical form using plot options.

Give only PLOT then individual commands.

• Look at properties of an element matrix.

– Use: EPRInt to print last element array.
– Use: EIGEn,<VECTor>,kk to print eigen-pairs for element ’kk’.

• Use: HELP to see commands (MANUal,,3 to get all).

• Use: SHOW,<option> to see size, dictionary, individual arrays.

• Use: HIST to see previous commands used; also to re-execute
previous command sets.

HIST,,3,12 ! Redo commands 3 to 12

FEAP Command Language Solution

• Example of EIGE from 8-node brick:

List 16 Command 1> eige,,5

*Command 1 * eige v: 5.00 0.00 0.00
t= 8.44 0.27
Eigenvalues for STIFFNESS of element 5

Matrix: Eigenvalue
row/col 1 2 3 4 5 6
1 6.174E+01 2.472E+01 2.465E+01 2.389E+01 2.323E+01 2.290E+01
2 2.179E+01 2.051E+01 1.891E+01 1.550E+01 1.218E+01 1.195E+01
3 1.165E+01 7.279E+00 6.691E+00 6.134E+00 4.060E+00 3.971E+00

• Note: No zero eigenvalues. Due to geometric stiffness!

FEAP Graphical Features

• Use: PLOT to enter graphics mode.

(Exit back to command language using ’e’).

• Use: HELP to see list of commands

• Plot contours:
CONTour 1 ! Display contour of DOF 1
STREss 1 ! Display stress component 1
• USE: POSTscript to output file for later use
POST ! Start PostScript file
MESH
OUTLine
POST ! End PostScript file
Produces set of files: FeapAAAA.eps, etc.
FEAP Graphical Features

• In command language mode can also produce files of time history

data
BATCh
TPLOt,,int ! Output each ’int’ steps
END
DISPl node dof
STREss element component
SHOW ! List to output file
Problem Solving with FEAP

• Running multiple problems:

• Prepare normal FEAP input files: Iprobx

(Here x denotes a problem number)

• Prepare master input file with data:

INCLude Iprob1
INCLude Iprob2
...
INCLude Iprobn
STOP

• Any STOP commands in Iprobx will be ignored.

• File names must be given in case sensitive mode.

Problem Solving with FEAP

• Example: Problem with mesh refinement

Master file Problem file Isquare

PARAmeter FEAP * * Mesh for square
n = 2 0 0 0 2 2 4
BLOCk
INCLude Isquare CART n n
PARAmeter 1 0 0
n = 4 2 5 0
3 5 5
INCLude Isquare 4 0 5
... ...
PARAmeter END
n = 64 BATCh
INCLude Isquare TANG,,1
STOP DISP,,1
END
Problem Solving with FEAP
• Produces sequence of meshes:
Problem Solving with FEAP

• To use this feature mesh must

– Either not have explicit definition of a node number or master
file must provide all node numbers as parameters.
– Let FEAP do counting for number of elements and nodes.

• Parameters may describe any data (including node and element

numbers). For example, range of Poisson ratios may be parame-
ters.

• Run one coarse mesh in INTEractive mode to establish the

command language statements needed.

• Test the file on two coarse meshes (or two parameters) to ensure
all works as planned.

• A NOPRint in mesh & batch execution file reduces output file size.
Programming Elements for FEAP in Fortran

• User elements all have single subprogram to interface

SUBROUTINE ELMTnn(D,UL,XL,IX,TL,S,R,NDF,NDM,NST,ISW)
where nn ranges from 01 to 50.

• In addition information is passed using COMMON statements.

INTEGER NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR
COMMON /CDATA/ NUMNP,NUMEL,NUMMAT,NEN,NEQ,IPR
Example: NEN is used to dimension some arrays.

• Best to use an include statement, e.g.:

INCLUDE ’cdata.h’
N.B. Windows case insensitive for file name, UNIX is not!
Programming Elements for FEAP in Fortran

• Main control switch (ISW) values:

ISW Description Command

0 Describe element SHOW,ELEM,nn
1 Input D(*) MATE ma
3 Compute stiffness TANG or
and residual UTAN
4 Output values STRE
5 Compute mass MASS
6 Compute residual FORM
8 Project values PLOT,STRE

• Currently options go to 26 (see Programmer Manual).

Programming Elements for FEAP in Fortran

• Typical displacement formulation element module

Call Quadrature()

Loop: L = 1,LINT IF: Isoparametric

Call Interp() IF: Isogeometric
Call Constitutive()
Form Resid
Form Tangent
End Loop

• Use quadrature and interpolation modules.

Programming for FEAP in Fortran

• Other user module options exist to modify FEAP. Examples:

– UMESH subprograms for user input data. Permits interfacing to

mesh generation programs (e.g., T-Splines & GiD from CIMNE
in Barcelona, Spain).

– UMACR subprograms for user command language statements.

– UMODEL subprograms for user stress-strain equations.

• UMESH2 used to interface T-spline refinement file to FEAP.

T-spline Solutions with FEAP

• Control records:

FEAP * * Problem Title

0 0 0 0 0 0

No number of nodes, etc. needed (sometime NDF needed).

• Material model records:

MATErial number
SOLId or MEMBrane only
ELASTIC .....
T-SPline ,, q1 q2 q3
....
! End of material

where q1, q2, q3 are number of quadrature points.

Input of T-spline File

• To input the T-spline mesh the statements are:

T-SPline
PLOT INTErvals = nint
MATErial number = ma
FILE = xxxxxxx.ext <or xxxxxxx.txt>
! End of T-spline inputs

• The PLOT and MATE records are optional.

If given they must precede the FILE record.
Input of T-spline File

• The plot interval value divides each T-Element into nint

subintervals in each direction for graphics display.

• Default value for nint is 1. Permitted range is 1 to 7.

• The material statement assigns all elements in the T-spline file

the material number ma (Default is 1).

• Multiple sets of T-SPline files may be included.