Table of Contents

1. Introduction

The VfeTools library from Visual Kinematics, Inc. is an object-based software development toolkit designed for generating finite element component matrices and vectors that add functionality to new or existing general purpose finite element programs, flexible body simulations, or special purpose tools tailored for a particular application. VfeTools enables programmers to compute standard quantities arising from the finite element method such as stiffness and mass matrices, consistent load vectors, stress and strain, etc., that are generated at the element level. Designed to provide robust and efficient element formulations in a modular, object-based architecture, VfeTools imposes few restrictions on the specific data structures used within a host application to maintain the computational grid and/or solution results. The basic features of VfeTools are summarized below:

Table of Contents

1.1 Module Summary

VfeTools is a set of modules that allow for the generation of many element-related matrices and vectors for 2D and 3D solid, shell, and beam element technologies. It contains three basic sets of modules: element formulations, material models and element and node utilities.

All VfeTools modules can be used independently of each other, with the exception of the MatlFun module which must be used with each element formulation module. The MatlFun module forms the abstract interface between each element formulation and the underlying material model.

Each instance of a module is termed an object. Users can generate as many objects as needed. The internal state of an object may be changed at any time to reflect, for example, the element topology it is meant to accept. As an example, users may change the current topology accepted by a Solid2D object from a 4 node quadrilateral to a 6 node triangle.

VfeTools can be used in a straightforward manner by the host application. The general methodology for generating element stiffness and mass matrices can be described by the following steps:

Table of Contents

1.2 Element Library

VfeTools currently supports solid, shell, beam, truss, and contact/interface element types for the disciplines of structural and thermal analysis. In addition, specialized discrete elements such as springs, masses, bulk capacitance and multipoint constraints are supported. The basic solid, shell, beam and truss element types are illustrated in Figure 1-1.


Figure 1-1, Basic Element Types
Within the context of a particular discipline, element formulation modules are partitioned by spatial dimension (2D, 3D). Within the context of 2D space, options are provided in structural analysis for axisymmetric, plane strain and plane stress cases. In thermal analysis options are provided for axisymmetric and planar cases. As a result one would use, for example, the Solid2D module for axisymmetric solids and the Shell3D module for shells in 3D space. The element degrees of freedom typically occur at the element nodes. For structural analysis the basic degrees of freedom are translations and, for elements such as shells and beams, translations and rotations. For thermal analysis the basic degrees of freedom are temperatures.

A variety of element topologies are supported. Basic topologies include lines, triangles, quadrilaterals, tetrahedra, pyramids, wedges and hexahedra. VfeTools supports linear and quadratic element orders in general, however cubic 2D and 3D solid elements are supported. There are two distinct forms for each element topology:

Numerical integration is used to integrate various properties across the element geometry. In many situations Gaussian quadrature is used. The total number of integration points in an element is the product of the number of element integration points over the nodal geometry of the element times the number of integration points over the element "section" such as a shell thickness or beam cross section. As explained in the following section, options are provided to choose between full or reduced integration.

Table of Contents

1.3 Element Technologies

All element modules contain formulations engineered to avoid locking mechanisms and rank deficiencies. This ensures that a formulation is available for each element type so that each type can be used without qualification over a large range of modeling and simulation situations. The general types of element formulations used may be summarized as follows:

1.3.1 Isoparametric Elements

Isoparametric elements are considered to be the basic finite element technology without the addition of any enhancements. This technology assumes a displacement or temperature field which varies isoparametrically with the element geometry. The integration order is selected to ensure no rank deficiencies in the element stiffness matrix. As a result of this approach, the user must use special care when selecting this formulation. In particular, standard isoparametric elements should be avoided in two specific cases: a) when the material behaves in a quasi-incompressible model, such as when the Poisson's ratio approaches 0.5, and b) when a small number of elements are used through a section to simulate bending. This technology is said to lock in the incompressible limit, and to be too stiff in bending particularly when using the linear element orders. This technology is recommended for 6 node triangular 2D solid elements, 10 node tetrahedral and all pyramid 3D solid elements.

1.3.2 Uniform Reduced Integration Elements

Uniform reduced integration elements are a simple extension of the isoparametric element technology in which the integration order is usually reduced by one order in one or more element directions. The resulting integration order is not sufficient to eliminate zero energy deformation modes in the element stiffness matrix. As a result this technology should generally be avoided except in special situations. For certain element types the zero energy modes do not propagate in any practical assembly of elements. In these cases, the reduced integration results in an element which is computationally more efficient and often eliminates element locking mechanisms yielding better results with fewer elements. This technology can be used with relative safety with 8 node quadrilateral 2D solid elements, 20 node hexahedral and 15 node wedge 3D solid elements and 8 node quadrilateral 3D shell elements. The zero energy modes can be robustly eliminated if an element is connected to other elements across it's edges or faces. This is discussed in more detail in the the sections below.

1.3.3 Assumed Natural Strain Elements

This technology is available for 6 node triangular shell elements only. The sampled natural strains are optiomal with respect to the derivative of a prescribed cubic field in a subparametric element. It is recommended as the most robust and best performing technology for this element type.

1.3.4 Mixed Elements

Mixed element technology is a Hu-Washizu based formulation that projects both strains and stresses into an optimal subspace with the goal of improving the overall performance of the element. This technology is available for all 2D solid element quadrilateral topologies and 3D solid hexahedral topologies and all 3D shell element quadrilateral topologies. This technology is recommended as the most robust and best performing technology for general use.

1.3.5 Enhanced Elements

Enhanced element technology is an "enhancement" of one of the above element formulations in which internal element freedoms are added to primarily improve the response of low order solid elements. The enhancement may improve the performance of the element in the incompressible limit, improve bending behavior or support accurate representation of transverse shear. Enhanced technology is currently available for 4 node quadrilateral 2D solid and 8 node hexahedral 3D solid, and 6 node pentahedral (wedge) 3D solid elements. For 4 node quadrilaterals, 2 internal freedoms are added and for 8 node hexahedra, 9 internal freedoms are added, and for 6 node pentahedra, 5 internal freedoms are added. Enhanced technology is used to implement transverse shear deformations in 2 node Kirchhoff 3D beam elements.

1.3.6 Element Naming Convention

The following naming convention is used to refer to element type, dimensionality, number of nodes and technology:
     Module_Name(number_of_nodes,technology)
The number_of_nodes and technology fields are optional.

The module name defines the element type and dimensionality. For example 2D solid elements are implemented in the Solid2D module. The number of nodes serves as a convenient and unique notation for the element shape and polynomial order within any module. The technologies, such as "isoparametric", are abbreviated as follows: isoparametric, ISOP; uniform reduced, URED. For example, the notation for the 2D solid Serendipity parabolic quadrilateral using uniform reduced technology would be as follows:

     Solid2D(8,URED)
If an asterisk, *, appears in place of the number of nodes or technology, it is meant to apply for all numbers of nodes or all technologies.

1.3.7 Solid Elements

Solid elements are volumetric elements in VfeTools. The modules Solid2D and Solid3D provide general purpose solid element formulations in 2D and 3D respectively.

Coordinates

     Solid2D(*)  x, y
     Solid3D(*)  x, y, z

Degrees of freedom

     Solid2D(*)  tx, ty
     Solid3D(*)  tx, ty, tz

Stresses and strains

     Solid2D(*)  sxx,syy,szz,sxy,syz,szx
                 exx,eyy,ezz,gxy,gyz,gzx
     Solid3D(*)  sxx,syy,szz,sxy,syz,szx
                 exx,eyy,ezz,gxy,gyz,gzx
In all 2D analysis syz, szx, gyz and gzx are zero. In 2D plane stress szz is zero, in 2D plane strain ezz is zero.

Integration and polynomial order

                  Normal  Reduced
     Solid2D(3)   1       1       linear triangle
     Solid2D(4)   2x2     1       linear quadrilateral
     Solid2D(6)   3       1       parabolic triangle
     Solid2D(8)   3x3     2x2     parabolic Serendipity quadrilateral
     Solid2D(9)   3x3     2x2     parabolic Lagrange quadrilateral

                  Normal  Reduced
     Solid3D(4)   1       1       linear tetrahedron
     Solid3D(5)   5       1       linear pyramid
     Solid3D(6)   1x2     1       linear wedge
     Solid3D(8)   2x2x2   1       linear hexahedron
     Solid3D(10)  4       1       parabolic tetrahedron
     Solid3D(13)  3x3x3   5       parabolic Serendipity pyramid
     Solid3D(14)  3x3x3   5       parabolic Lagrange pyramid
     Solid3D(15)  3x3     3x2     parabolic Serendipity wedge
     Solid3D(18)  3x3     3x2     parabolic Lagrange wedge
     Solid3D(20)  3x3x3   2x2x2   parabolic Serendipity hexahedron
     Solid3D(27)  3x3x3   2x2x2   parabolic Lagrange hexahedron
Gaussian quadrature is used for all quadrilateral and hexahedral elements and in the prismatic direction of wedge elements.

1.3.8 Shell Elements

Shell elements are degenerate solid elements in VfeTools in which a given direction, the shell thickness, may be assumed to be small. The Shell3D module provides general purpose shell element formulations in 3D space.

Coordinates

     Shell3D(*)  x, y, z

Degrees of freedom

     Shell3D(*)  tx, ty, tz, rx, ry, rz

Stress resultants and midsurface strains and curvatures in shell local coordinate system

     Shell3D(*)  Nxx,Nyy,Nxy,Mxx,Myy,Mxy,Qxz,Qyz
                 exx,eyy,gxy,kxx,kyy,kxy,gxz,gyz
The number of integration points used by each shell element type appears below. The integration orders are for the shell surface integration. Gaussian quadrature is used for all quadrilateral elements.
                  Normal  Reduced
     Shell3D(3)   1       1       linear triangle
     Shell3D(4)   2x2     1       linear quadrilateral
     Shell3D(6)   3       1       parabolic triangle
     Shell3D(8)   3x3     2x2     parabolic Serendipity quadrilateral
     Shell3D(9)   3x3     2x2     parabolic Lagrange quadrilateral

1.3.9 Beam Elements

The Beam3D module provides general purpose beam element formulations in 3D space.

Coordinates

     Beam3D(*)  x, y, z

Degrees of freedom

     Beam3D(*)  tx, ty, tz, rx, ry, rz

Stress resultants and reference axis strains, curvatures and twist in beam local coordinate system

     Beam3D(*)  Nxx,Myy,Mzz,Txx,Qxy,Qzx
                exx,kyy,kzz,twist,gxy,gzx

The number of integration points used by each beam element type appears below. The integration orders are for the beam axis integration. Gaussian quadrature is used for all elements.

                  Normal  Reduced
     Beam3D(2)    1       1       linear line
     Beam3D(3)    2       1       parabolic line

Table of Contents

1.4 Accuracy of Finite Element Formulations

The purpose of this section is to present the results of the solid, shell and beam element formulations on a set of standard problems which illustrate element accuracy. The results which appear below are a subset of the so-called MacNeal-Harder tests which are described in "A Proposed Standard Set of Problems to Test Finite Element Accuracy", by R. H. MacNeal and R. L. Harder, published in Vol. 1, No. 1 of Finite Elements in Analysis and Design, 1985, North-Holland Publishing Co., Amsterdam, The Netherlands. In each case the element formulation used is the default for each element type. For comparison, the published results of other well known element technologies are included if available. Unless specified, all finite element results are normalized by the reference solution of MacNeal and Harder.

1.4.1 Straight Cantilever Beam

This problem tests the accuracy of rectangular, trapezoidal and parallelogram shaped elements for a cantilever beam under extension, in-plane and out-of-plane bending and twisting. The problem is illustrated in Figures 1-1a, 1-1b and 1-1c. The length of the beam is 6.0, the width is 0.2 and the depth is 0.1 . Young's modulus is 1.0E+07, and Poisson's ratio is 0.30 . The element edges of the trapezoidal and parallelogram meshes are inclined at 45 degrees. A constant mesh of 6 by 1 elements is used. The loading is unit forces at the free end.

The value for the normalized tip displacement in the direction of the load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 3.0E-05 for extension, 0.1081 for in-plane shear, 0.4321 for out-of-plane shear and 0.03208 for twist.


Figure 1-1a, Straight Cantilever Beam, Rectangular Elements

Figure 1-1a, Straight Cantilever Beam, Trapezoidal Elements

Figure 1-1a, Straight Cantilever Beam, Parallelogram Elements
Tip loading
direction
                   (a) Rectangular elements
                   Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
Extension          0.995       0.997       0.999       0.999
In-plane shear     0.903       0.983       0.987       0.994
Out-of-plane shear 0.980       0.992       0.994       0.994
Twist              0.946       0.954       0.949       0.949

NASTRAN            QUAD4                   QUAD8
Extension          0.995                   0.999
In-plane shear     0.904                   0.987
Out-of-plane shear 0.986                   0.991
Twist              0.941                   0.950

                   (b) Trapezoidal elements
                   Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
Extension          0.988       0.994       0.999       0.999
In-plane shear     0.043       0.891       0.965       0.989
Out-of-plane shear 0.971       0.979       0.995       0.994
Twist              0.927       0.947       0.950       0.952

NASTRAN            QUAD4                   QUAD8
Extension          0.996                   0.999
In-plane shear     0.071                   0.946
Out-of-plane shear 0.968                   0.998
Twist              0.951                   0.943

                   (c) Parallelogram elements
                   Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
Extension          1.011       0.993       0.999       0.999
In-plane shear     0.017       0.888       0.987       0.987
Out-of-plane shear 0.980       0.986       0.994       0.995
Twist              0.855       0.956       0.949       0.949

NASTRAN            QUAD4                   QUAD8
Extension          0.996                   0.999
In-plane shear     0.080                   0.995
Out-of-plane shear 0.977                   0.985
Twist              0.945                   0.965

1.4.2 Curved Beam

This problem tests the accuracy of in-plane and out-of-plane loading for a curved beam. The problem is illustrated in Figure 1-2. The inner radius of the beam is 4.12, the outer radius os 4.32 making an arc of 90 degrees. The thickness is 0.1 . Young's modulus is 1.0E+07, and Poisson's ratio is 0.25 . A constant mesh of 6 by 1 elements is used. The loading is unit forces at the tip in the in-plane and out-of-plane directions.

The value for the normalized tip displacement in direction of load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.08734 for the in-plane load and 0.5022 for the out-of-plane load.


Figure 1-2, Curved Beam
Tip loading
direction       Normalized tip displacement in direction of load
                Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
In-plane        0.833       0.998       1.021       1.011
Out-of-plane    0.936       0.938       0.957       0.957

NASTRAN         QUAD4                   QUAD8
In-plane        0.833                   1.007
Out-of-plane    0.951                   0.971

                Solid3D(8)  Solid3D(10) Solid3D(15) Solid3D(20) Solid3D(27)
In-plane        0.822       0.887       0.947       1.016       1.021
Out-of-plane    0.790       0.868       0.908       0.929       0.932

NASTRAN         HEXA(8)                             HEX20(R)
In-plane        0.880                               1.006
Out-of-plane    0.849                               0.959

                Beam3D(2)   Beam3D(3)
In-plane        1.022       1.042
Out-of-plane    0.991       1.003

1.4.3 Twisted Beam

This problem tests the accuracy of in-plane and out-of-plane loading for a twisted beam. It tests the effect of warp on shell elements. The problem is illustrated in Figure 1-3. The length of the beam is 12.0, the width is 1.1, the thickness (depth) is 0.32. The twist of the beam from the fixed end to the free end is 90 degrees. Young's modulus is 29.0E+06, and Poisson's ratio is 0.22 . A constant mesh of 12 by 2 nodal spaces is used. The loading is unit forces at the tip in the in-plane and out-of-plane directions.

The value for the normalized tip displacement in direction of load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.005424 for the in-plane load and 0.001754 for the out-of-plane load.


Figure 1-3, Twisted Beam
Tip loading
direction       Normalized tip displacement in direction of load
                Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
In-plane        0.987       0.984       0.997       1.001
Out-of-plane    0.981       0.987       1.003       1.003

NASTRAN         QUAD4                   QUAD8
In-plane        0.985                   0.998
Out-of-plane    0.993                   0.998

1.4.4 Rectangular Plate

This problem tests plate bending for a rectangular plate. The test is run using with two different aspect ratios of the plate 1 by 1 and 1 by 5. The problem is illustrated in Figure 1-4. The dimensions of the plate are a = 1.0, b = 1.0 or b = 5.0 . The thickness of the plate is 0.0001 for shell element testing and 0.01 for solid element testing. Young's modulus is 1.7472E+07, and Poisson's ratio is 0.3 . The boundary conditions on the plate edges are simply supported or clamped. The loading is a uniform pressure of 1.0E-04 or a central load at the intersection of the two symmetry boundaries of 4.0E-04.

The value for the normalized displacement at the center (intersection of two symmetry boundaries) of the plate is used as a measure of solution accuracy. We have assumed the solutions for the normalization of results to be the following for shell elements, multiply by 1.0e-06 for solid elements since the thickness for solid elements is 100 times that used for shell elements.

Boundary        Aspect      Displacement at center of plate
supports        ratio b/a   Uniform pressure      Concentrated force

Simple          1.0         4.062                 11.60
Simple          5.0        12.97                  16.96
Clamped         1.0         1.26                   5.60
Clamped         5.0         2.56                   7.23


Figure 1-4, Rectangular Plate
Simple support uniform load, aspect ratio = 1.0
Normalized lateral deflection at center
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       0.987       1.093       0.985       1.003
 4       0.997       1.029       1.001       1.003
 6       0.999       1.014       1.001       1.001
 8       0.999       1.008       1.000       1.000

NASTRAN  QUAD4                   QUAD8
 2       0.981                   0.927
 4       1.004                   0.996
 6       1.003                   0.999
 8       1.002                   1.000

Simple support uniform load, aspect ratio = 5.0
Normalized lateral deflection at center
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       0.890       0.595       1.279       1.245
 4       0.941       0.779       0.986       0.995
 6       0.970       0.844       0.990       0.996
 8       0.982       0.890       0.994       0.998

NASTRAN  QUAD4                   QUAD8
 2       1.052                   1.223
 4       0.991                   1.003
 6       0.997                   1.000
 8       0.998                   1.000

X

Clamped support concentrated load, aspect ratio = 1.0
Normalized lateral deflection at center
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       0.881       0.001       1.116       1.116
 4       0.971       0.879       0.902       1.011
 6       0.988       0.966       0.990       1.004
 8       0.994       0.987       0.984       1.003

NASTRAN  QUAD4                   QUAD8
 2       0.934                   1.076
 4       1.010                   0.969
 6       1.012                   0.992
 8       1.010                   0.997

Clamped support concentrated load, aspect ratio = 5.0
Normalized lateral deflection at center
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       0.319       0.001       0.423       0.423
 4       0.841       0.401       0.942       0.950
 6       0.922       0.556       0.957       0.988
 8       0.954       0.672       0.969       0.994

NASTRAN  QUAD4                   QUAD8
 2       0.519                   0.542
 4       0.863                   0.754
 6       0.940                   0.932
 8       0.972                   0.975

1.4.5 Scordelis-Lo Roof

This problem tests the accuracy of shell elements for a cylindrical shell. The problem is illustrated in Figure 1-5. The radius of the shell is 25.0, the length is 25.0, the thickness is 0.25, Young's modulus is 4.32E+08, and Poisson's ratio is 0.0 . The loading = 90.0 per unit are in the -Z direction. The constraint on the unlabelled curved edge is tx = tz = 0.

The value for the midside vertical displacement is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.3024 .


Figure 1-5, Scordelis-Lo Roof
Normalized vertical deflection at midpoint of free edge
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       1.387       0.872       0.832       0.863
 4       1.048       0.863       0.955       0.961
 6       1.017       0.970       0.987       0.988
 8       1.010       0.992       0.993       0.993
10       1.009       0.997       0.994       0.994

NASTRAN  QUAD4                   QUAD8
 2       1.376                   1.021
 4       1.050                   0.984
 6       1.018                   1.002
 8       1.008                   0.997
10       1.004                   0.996

1.4.6 Spherical Shell

This problem tests the accuracy of shell elements for a spherical shell. The problem is illustrated in Figure 1-6. The radius of the shell is 10.0, the thickness is 0.04, Young's modulus is 6.825E+07 and Poisson's ratio is 0.3 . The loading consists of concentrated forces on the shell quarter model as shown. For triangular elements each quadrilateral is subdivided into two triangles.

The radial displacement under the load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.094 .


Figure 1-6, Spherical Shell
Normalized radial deflection at load point
         Shell3D(4)  Shell3D(6)  Shell3D(8)  Shell3D(9)
 2       0.909       0.910       1.333       1.774
 4       0.999       0.860       0.646       0.889
 8       0.990       1.012       0.961       0.989
12       0.990       1.019       0.987       0.991

NASTRAN  QUAD4                   QUAD8
 2       0.972                   0.025
 4       1.024                   0.121
 8       1.005                   0.823
12       0.998                   0.992

ABAQUS   S4          STRI65      S8R5        S9R5
 2       0.444       0.229       0.796       0.669
 4       0.925       0.916       0.951       0.944
 8       0.974       0.967       0.982       0.982
12       0.980       0.974       0.982       0.982

Normalized radial deflection at load point
         Solid3D(8)  Solid3D(20) Solid3D(27)
 4       0.034       0.022       0.367
 8       0.702       0.571       0.980
12       0.951       0.919       1.007

NASTRAN  HEXA(8)     HEX20(R)
 4       0.039       0.162
 8       0.730       0.776
12       0.955       0.972
Table of Contents

1.5 Computed Element Quantities

All element modules are designed to compute a series of quantities that are of interest to the user. These quantities are outlined below along with a brief explanation of their origin and use.

Table of Contents

1.6 Degrees of Freedom

Finite Element simulations consist of decomposing the computational domain into a collection of discrete entities denoted "finite elements". Inside these elements, some field is approximated in a known way. For example, in a 3 node triangular element, most fields are linearly approximated in the element. This approximation can be uniquely described if one knows the values attained by the field at the nodes of the triangle. At times, the element approximation technology demands more than a simple linear or quadratic interpolation, and a complete field representation may require additional internal (also called elemental) unknowns. As these elemental degrees of freedom are connected to a single element only, they can be statically condensed out at the element matrix level.

The collection of all nodal and elemental unknowns is denoted as the degree of freedom map or "dof map" of the element. Typical degrees of freedom in structural mechanics are the Cartesian displacements or translations. Shell and beam elements require nodal rotational degrees of freedom in addition to the translations.

In heat transfer simulations, typical degrees of freedom can be the temperature and/or its gradient. As the analysis type and the element technology are selected along with the element topology, the number and type of degrees of freedom associated with the element are determined. Note that the number of degrees of freedom and dof map depends on the discipline (e.g., stress analysis, thermal analysis transfer, etc.), the technology utilized (e.g., isoparametric, enhanced strains, etc.), and the element topology and order (e.g. quadratic triangle, linear quadrilateral).

Every element in VfeTools provides a method to query the element for the number of degrees of freedom, as well as the dof map. The number of degrees of freedom, nedofs, is queried as follows:

       vfe_Solid2DNumDof (solid2d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
while the dof map, vectors loc and tag (described below), are queried as follows:

       vfe_Solid2DDofMap (solid2d,SYS_ANALYSIS_STRUCTURAL,loc,tag);
The vector loc is dimensioned nedofs long, and contains for every degree of freedom in the element, the node index into the element connectivity associated with the degree of freedom. A node index of 1 indicates that the degree of freedom is associated with the first node in the element connection list. If the degree of freedom happens to be internal to the element, the node index is zero.

The array tag is also dimensioned nedofs, and describes the type of each of the degrees of freedom. The degree of freedom types are specified by a defined constant as follows:

In 2D axisymmetric analysis SYS_DOF_TX coincides with the radial coordinate r, SYS_DOF_TY coincides with the axial coordinate z and SYS_DOF_TZ coincides with the circumferential coordinate theta, Note that SYS_DOF_Si always equals SYS_DOF_S0 + i.

Table of Contents

1.7 Material Models

Material models in VfeTools are of two basic types: 1) primitive materials and 2) element materials which contain one or more primitive materials and to which element properties such as shell element thickness have been added.

Primitive materials are characterized by a material response which excludes any element geometry. Examples include stress - strain relations, heat flux - temperature gradient relations, etc. Within VfeTools all modules which implement primitive materials are given names with the string "Mat" appended. As an example, the LinMat module implements linear isotropic, orthotropic and anisotropic materials. Primitive materials include:

Certain primitive materials, such as plasticity, require additional history information to be supplied at each step in the analysis. The type and quantity of history information is dependent upon the character of the primitive material. The user is required to manage the memory required for the material history.

Element materials are material responses which include the effects of element properties. Examples include layered composites which include element thickness and offset information as well as primitive material definitions for each layer. An element material need not necessarily include a primitive material. For example an element material for an elastic spring may simply contain a single scalar spring constant. Within VfeTools all modules which implement element materials are given names with the string "Prop" appended. As an example, the ShellProp module implements layered composites for shell element formulations. Element materials include:

Table of Contents

1.8 Error Estimation and Mesh Adaptation

In the finite element method, the displacement field is usually explicitly assumed to be continuous between elements. However, under such an assumption, the computed stress field will generally be discontinuous between elements to some degree depending upon the accuracy of the finite element discretization. This fact can be exploited to compute an estimate of the error of a given discretization and suggest a new discretization size to obtain an improved solution.

An effective method for estimating error is to compute the strain energy error on an element by element basis and compare it to the total strain energy. This computation requires an integration over the element of the difference between the exact and computed nodal stresses. The computed nodal stresses are determined by the element given the displacement field. An approximation to the exact nodal stresses may be obtained by averaging the computed stresses between adjacent elements where the stresses are expected to be continuous. This approximation process is termed "stress recovery" and the method suggested above for recovering stress is one of many possible methods.

A number of element modules in VfeTools provide methods to compute the element energy and energy error. Supported element types include solid, shell, beam and truss element types. For structural analysis the StrsAdapt function is provided. This function computes the element strain energy, setot and strain energy error, seerr. In addition, useful quantities such as the characteristic element length, h, effective polynomial order, p, and dimension of the element, d, are returned.

       vfe_Solid2DStrsAdapt (solid2d,x,u,strss,setot,seerr,&h,&p,&d);
The element characteristic length is the edge length of an element with the same area or volume assuming equal length, straight edges and maximum minimum angles between edges. The element polynomial order is approximately the polynomial order of displacement interpolation along an element edge. The element dimension is the number of local element dimensions spanned by the element connectivity. For example, 3D shells and 2D solids have an element dimensionality of 2, 3D solids have an element dimensionality of 3.

A similar function is available in heat transfer analysis for computing total heat energy and heat energy error given the nodal temperatures and recovered nodal flux.

       vfe_Solid2DHFlxAdapt (solid2d,x,u,hflxs,hetot,heerr,&h,&p,&d);
To obtain an estimate of the global error the elemental energies and energy errors are computed using the functions above and summed for all applicable elements. Let E represent the total strain energy and e the energy error obtained in this manner, then the total solution error err may be obtained as:

   err = sqrt(e/E)
A solution is considered converged when the solution error is less than some predefined tolerance, tol, usually 5%. It is interesting to note that is has been shown that if seerr is equal for all elements, then the finite element model using a given number of elements is the optimal one. This fact suggests a adaptation process in which a new characteristic element length is computed for each element to equilibrate the solution error and reduce it to the convergence tolerance.

First an estimate of the number of elements in the adapted mesh is determined which yields a solution error equal to the convergence tolerance. This is done by assuming that the mesh size will be uniformly scaled by a scaling factor, f. Since the energy error, seerr, in any element is of the order h**p where h and p are the element characteristic length and polynomial order, for a mesh characterized by p, the energy error in the new mesh scaled by f is related to the current energy error by

   e_new = (f**p)*e
Dividing by E and solving for f yields

   f = (tol/err)**(1/p)
The number of elements in the new mesh is given by

   N_new = N/(f**d)
where N is the current number of elements and d is the element dimensionality in the mesh. Assuming that the element energy error is equilibrated in the new mesh, the target element energy error is given as

   seerr_tar = tol*E/sqrt(N_new)
Although the target element energy error is computed assuming a uniform scaling of the current mesh, it is more realistic and efficient to adjust the mesh size locally according to the current contribution of each element to the solution error. This suggests that the new characteristic element size for each element be computed as

   h_new = (seerr_tar/seerr)**(1/p)*h
New element sizes computed in this fashion may be given to an automatic mesh generator to generate an entirely new mesh which will be refined or derefined in a very localized way.

Table of Contents

1.9 Performance Issues

VfeTools pays special attention to performance, as finite element calculations are often compute intensive and can take from a few seconds to several hours, depending on problem size.

To address performance issues from a strict operation count, cache size and processor architecture point of view, VfeTools modules avoid zero arithmetic.

Table of Contents

1.10 Compiling and Linking a VfeTools Application

To use VfeTools on a particular computer platform, the VfeTools source must be compiled and linked to an application. Either the object files may be used directly or they may be installed in an object library archive so that the loader may selectively relocate only the objects which are required. VfeTools is written in ANSI C.

VfeTools is basically platform independent and as a result no user defined C preprocessor directives are required to compile VfeTools on any supported platform. However for performance reasons certain platform specific compiler options may be necessary to notify the compiler that pointers are not "aliased", etc.

It is suggested that during the development cycle that the source be conditionally compiled with error checking by defining VKI_CHECK_ERRS as described in base library.

For example, on SGI systems, create a directory SGI under lib to hold the final devtools.a archive file. Then from the devtools/src/vfe directory compile using

cc -c -O2 -I.. *.c
creating .o files in the vfe directory. As mentioned above it is suggested to use the appropriate compiler option to specify that pointers are not aliased plus additional optimizations as follows:

To place the object files in an archive file issue

ar rs ../../lib/SGI/devtools.a *.o
The object files may be deleted after the devtools.a archive is successfully created.

To compile the base source, change directory to devtools/src/base and compile using

cc -c -O2 *.c
creating .o files in the base directory. To add these objects to the previously created devtools.a archive issue

ar rs ../../lib/SGI/devtools.a *.o
To compile the vis source, change directory to devtools/src/vis and compile using

cc -c -O2 *.c
creating .o files in the vis directory. To add these objects to the previously created devtools.a archive issue

ar rs ../../lib/SGI/devtools.a *.o
Again object files may deleted at this time. At this point the devtools.a archive contains all vfe, vis, and base objects. Place the devtools.a archive immediately before the graphics subsystem libraries in the load line. A devtools.a archive must be built for each computer platform using the methodology outlined above.

Table of Contents

1.11 A First Program - C Version

As an example of a simple VfeTools application the following program computes the stiffness matrix for a 2D plane stress 3 node triangular element. The material model is assumed isotropic, linear elastic, with Young's modulus and Poisson's ratio of 1.0 and 0.3, respectively. The resulting stiffness matrix is printed to standard output.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"

static Vdouble x[3][3] = {
  {0.,0.,0.}, {1.,0.,0.}, {0.,1.,0.} };

/*----------------------------------------------------------------------
                      Generate Stiffness Matrix for 3 Node Triangle
----------------------------------------------------------------------*/
int
main()
{
   Vint i,j;
   Vdouble prop[2];
   Vdouble kl[21];
   vfe_Solid2D *solid2d;
   vfe_MatlFun *matlfun;
   vfe_LinMat  *linmat;

                   /* Create material object and set properties */ 
   linmat  = vfe_LinMatBegin ();
   vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
   prop[0] = 1.0;
   prop[1] =  .3;
   vfe_LinMatSetElasProp (linmat,prop);

                   /* Create MatlFun object and load material functions */ 
   matlfun = vfe_MatlFunBegin ();
   vfe_LinMatMatlFun (linmat, matlfun);

                   /* Create 2D solid element formulation */ 
   solid2d = vfe_Solid2DBegin ();
   vfe_Solid2DSetTopology (solid2d,SYS_SHAPETRI,2,0);

                   /* Register material functions */ 
   vfe_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun);

                   /* Form and print stiffness matrix */ 
   vfe_Solid2DStiff (solid2d,x,kl);
   printf("Lower triangle of stiffness matrix\n");
   for(i = 0; i < 6; i++) {
      for(j = 0; j <= i; j++) {
         printf("%12.4e ",kl[i*(i+1)/2+j]);
      }
      printf("\n");
   }
                   /* Delete objects */ 
   vfe_Solid2DEnd (solid2d);
   vfe_MatlFunEnd (matlfun);
   vfe_LinMatEnd  (linmat);
   return 0;
}

The output of this example program appears below. Only the lower half of the matrix is printed.

Lower triangle of stiffness matrix
  7.4176e-01 
  3.5714e-01   7.4176e-01 
 -5.4945e-01  -1.6484e-01   5.4945e-01 
 -1.9231e-01  -1.9231e-01   0.0000e+00   1.9231e-01 
 -1.9231e-01  -1.9231e-01   0.0000e+00   1.9231e-01   1.9231e-01 
 -1.6484e-01  -5.4945e-01   1.6484e-01   0.0000e+00   0.0000e+00   5.4945e-01 

Table of Contents

1.12 A First Program - C++ Version

The following program is a listing of the C++ version of the same "A First Program" listed above which used C language bindings.

#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"

static Vdouble x[3][3] = {
  {0.,0.,0.}, {1.,0.,0.}, {0.,1.,0.} };

/*----------------------------------------------------------------------
                      Generate Stiffness Matrix for 3 Node Triangle
----------------------------------------------------------------------*/
int
main()
{
   Vint i,j;
   Vdouble prop[2];
   Vdouble kl[21];
   vfe_Solid2D *solid2d;
   vfe_MatlFun *matlfun;
   vfe_LinMat  *linmat;

                   /* Create material object and set properties */ 
   linmat  = new vfe_LinMat ();
   linmat->Def (SYS_MAT_ISOTROPIC);
   prop[0] = 1.0;
   prop[1] =  .3;
   linmat->SetElasProp (prop);

                   /* Create MatlFun object and load material functions */ 
   matlfun = new vfe_MatlFun ();
   linmat->MatlFun (matlfun);

                   /* Create 2D solid element formulation */ 
   solid2d = new vfe_Solid2D ();
   solid2d->SetTopology (SYS_SHAPETRI,2,0);

                   /* Register material functions */ 
   solid2d->SetObject (VFE_MATLFUN,matlfun);

                   /* Form and print stiffness matrix */ 
   solid2d->Stiff (x,kl);
   printf("Lower triangle of stiffness matrix\n");
   for(i = 0; i < 6; ++i ) {
      for(j = 0; j <= i; ++j) {
         printf("%12.4e ",kl[i*(i+1)/2+j]);
      }
      printf("\n");
   }
                   /* Delete objects */ 
   delete solid2d;
   delete matlfun;
   delete linmat;
   return 0;
}
Table of Contents

1.13 A First Program - FORTRAN Version

The following program is a listing of the FORTRAN version of the same "A First Program" listed above which used C language bindings. Although the FORTRAN language does not distinguish between upper and lower case characters, the example retains the conventions employed in the C version for clarity. Note that the prefix vfe_ is replaced by vfef_.

C-----------------------------------------------------------------------
C                   Generate Stiffness Matrix for 3 Node Triangle
C-----------------------------------------------------------------------
      program main
      include 'base/fortran/base.inc'
      include 'vfe/fortran/vfe.inc'
      integer i,j
      double precision prop(2)
      double precision kl(21)
      double precision solid2d
      double precision matlfun
      double precision linmat
      double precision x(3,3)
      data x / 0.,0.,0., 1.,0.,0., 0.,1.,0. /

C                  Create material object and set properties
      call vfef_LinMatBegin (linmat)
      call vfef_LinMatDef (linmat,SYS_MAT_ISOTROPIC)
      prop(1) = 1.0
      prop(2) =  .3
      call vfef_LinMatSetElasProp (linmat,prop)

C                  Create MatlFun object and load material functions
      call vfef_MatlFunBegin (matlfun)
      call vfef_LinMatMatlFun (linmat, matlfun)

C                  Create 2D solid element formulation
      call vfef_Solid2DBegin (solid2d)
      call vfef_Solid2DSetTopology (solid2d,SYS_SHAPETRI,2,0)

C                  Register material functions
      call vfef_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun)

C                  Form and print stiffness matrix
      call vfef_Solid2DStiff (solid2d,x,kl)
      write(6,*) "Lower triangle of stiffness matrix"
      do i = 0,5
         write(6,1000) (kl(i*(i+1)/2+j+1),j=0,i)
      enddo
C                  Delete objects
      call vfef_Solid2DEnd (solid2d)
      call vfef_MatlFunEnd (matlfun)
      call vfef_LinMatEnd  (linmat)
C
1000  format(1x,6e12.4)
      stop
      end
Table of Contents

1.14 A First Program - C# Version

The following program is a listing of the C# version of the same "A First Program" listed above which used C language bindings. Note that the prefix vfe_ is replaced by vfe..

using System;
using System.Runtime.InteropServices;
using System.Reflection;
using System.Text;
using DevTools;

public class intro1 {

   public static double [] x = { 0.0,0.0,0.0, 1.0,0.0,0.0, 0.0,1.0,0.0 };

/*----------------------------------------------------------------------
                      Generate Stiffness Matrix for 3 Node Triangle
----------------------------------------------------------------------*/
   public static void Main() {

      int i = 0,j = 0;
      double [] prop = new double[2];
      double [] kl   = new double[21];
      IntPtr solid2d;
      IntPtr matlfun;
      IntPtr linmat;

                   /* Create material object and set properties */ 
      linmat  = vfe.LinMatBegin ();
      vfe.LinMatDef (linmat,vsy.SYS_MAT_ISOTROPIC);
      prop[0] = 1.0;
      prop[1] = 0.3;
      vfe.LinMatSetElasProp (linmat,prop);

                   /* Create MatlFun object and load material functions */ 
      matlfun = vfe.MatlFunBegin ();
      vfe.LinMatMatlFun (linmat, matlfun);

                   /* Create 2D solid element formulation */ 
      solid2d = vfe.Solid2DBegin ();
      vfe.Solid2DSetTopology (solid2d,vsy.SYS_SHAPETRI,2,0);

                   /* Register material functions */ 
      vfe.Solid2DSetObject (solid2d,vfe.VFE_MATLFUN,matlfun);

                   /* Form and print stiffness matrix */ 
      vfe.Solid2DStiff (solid2d,x,kl);
      Console.Write("Lower triangle of stiffness matrix\n");
      for(i = 0; i < 6; i++) {
         for(j = 0; j <= i; j++) {
            Console.Write("{0} ",kl[i*(i+1)/2+j]);
         }
         Console.Write("\n");
      }
                   /* Delete objects */ 
      vfe.Solid2DEnd (solid2d);
      vfe.MatlFunEnd (matlfun);
      vfe.LinMatEnd  (linmat);
   }
}