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:

Support for solid and interface elements in both 2D and 3D.
2D elements support plane strain, plane stress and axisymmetry.
Generation of stiffness, consistent or diagonal mass and geometric stiffness matrices for structural analysis.
Generation of conductance and consistent and diagonal capacitance matrices for thermal analysis.
Generation of consistent nodal force vectors for edge and/or face loading. Generation of body forces due to accelerations and rotations.
Robust element technologies that avoid locking mechanisms and/or rank deficiencies in both the infinitesimal and large deformation regimes.
Multipoint constraints for rigid links, distributed coupling forces, etc.
Methods for easy integration into new or existing finite element applications.
Element properties include variable thickness and laminated shells, tapered and general cross section beams.
Element independent corotational utilities for implementing large rotation capabilities. Geometric load stiffness may be computed. The corotation formulation is also useful for enhancing the response of linear element formulations which are sensitive to rigid body motions.
Abstract interface between elements and material models allowing users to add their own material models.
Linear material models include isotropic, orthotropic and anisotropic temperature dependent materials.
Nonlinear material models include large strain plasticity and hyperelasticity.
All material models are temperature dependent
Disciplines include structural and thermal analysis.
Object-based architecture, written in ANSI C with C++, FORTRAN, and C# language bindings.
Hardware and operating system independence.

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.

Elements

Solid2D     2D solid element.
Solid3D     3D solid element.
Shell3D     3D monocoque or composite shell element.
Mem3D       3D membrane element.
Beam3D      3D beam and shell stiffener element.
Truss3D     3D truss or cable element.
Inter2D     2D interface element.
Inter3D     3D interface element.
           
Bulk        Volumetric mass or heat capacitance
ConMass     Concentrated or degree of freedom mass
ConElas     Grounded or degree of freedom elastic spring
Gap         Point to point gap element
Spring      Spring element
           
MPC         Multi-point constraint

Material Models

LinMat      Linear isotropic, orthotropic or anisotropic materials
PlasMat     Isotropic, rate independent, large strain plasticity
HyperMat    Mooney-Rivlin and Ogden materials
           
ShellProp   Shell element wall properties
BeamProp    Beam element section properties
InterProp   Interface and gap element properties
           
MatlFun     Define function pointers to material model

Element and Node Utilities

ElemTran    Element transformations for rotated coordinate systems.
NodeGeom    Node accelerations, rigid body modes, inertia relief.
Corot       Corotational large rotation 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:

For each element type (solid, shell, beam, etc.) an element formulation object is instantiated. The current topology accepted by the Solid2D object may be changed at any time by calling vfe_Solid2DSetTopology. As a result, the same Solid2D object can be used for either, say, an isotropic 3 node triangular or an orthotropic 9 node quadrilateral element. The C source code to create a Solid2D object for a 4 node quadrilateral is given by
```
  vfe_Solid2D *solid2d;
       
  solid2d = vfe_Solid2DBegin ();
  vfe_Solid2DSetTopology (solid2d,SYS_SHAPEQUAD,2,0);
```

Prior to requesting computed quantities from an element object, a material model object must be created and the material properties defined. The LinMat module is designed to support linear isotopic, orthotropic or anisotropic materials. Material properties such as elastic moduli, coefficients of thermal expansion, density, etc. may be entered. The C source code to create a linear isotropic material model is given below
```
  vfe_LinMat *linmat;
  Vdouble young   = 10000000.0;
  Vdouble poisson = 0.3;
  Vdouble props[2];

  linmat = vfe_LinMatBegin ();
  vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
  props[0] = young;
  props[1] = poisson;
  vfe_LinMatSetElasProp (linmat,props);
  vfe_LinMatSetDensity (linmat,0.0000133);
```

Once the element and material model objects have been created, the material model is connected to the element formulation using a set of function pointers. The MatlFun module implements the abstract interface between the element formulation and the underlying material model. Each material model module, such as LinMat, has a function which fills a MatlFun object with the function pointers particular to that material model module. Typical functions which are registered in a MatlFun object allow the element formulation to query for quantities such as stress as a function of strain, material moduli, moduli sparsity, etc. The MatlFun object is then registered with the element formulation object completing the connection. The C source code illustrating this procedure for the element formulation and material model objects created above is given as follows
```
  vfe_MatlFun *matlfun;

  matlfun = vfe_MatlFunBegin ();
  vfe_LinMatMatlFun (linmat,matlfun);
  vfe_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun);
```

The host application must then prepare itself to receive the output quantities, such as stiffness, computed by the element formulation. The first task is to determine the degree of freedom (dof) structure, termed the degree of freedom "map", or dof map of the element. Elements can be queried for the number of degrees of freedom they contain and the so-called dof map. Degrees of freedom may be either nodal and/or elemental, depending upon the element technology. The dof map consists of a location and tag for each degree of freedom in the element. The location is a one-based nodal index into the element connectivity (zero implies an elemental rather than a nodal freedom). The tag is an associated enumerated type such as translation in the x direction, temperature, etc. Note that the defined constant VFE_MAXDOF is equal to the maximum number of degrees of freedom used by any element formulation in VfeTools (about 80). Source code is as follows:
```
  Vint nedofs;
  Vint loc[VFE_MAXDOF];
  Vint tag[VFE_MAXDOF];

  vfe_Solid2DNumDof (solid2d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
  vfe_Solid2DDofMap (solid2d,SYS_ANALYSIS_STRUCTURAL,loc,tag);
```

The dof map produced by the call to vfe_Solid2DDofMap above is equivalent to the following:

  Vint loc[8] = { 1,1, 2,2, 3,3, 4,4 };
  Vint tag[8] = { SYS_DOF_TX,SYS_DOF_TY, SYS_DOF_TX,SYS_DOF_TY,
                  SYS_DOF_TX,SYS_DOF_TY, SYS_DOF_TX,SYS_DOF_TY };

Once the degree of freedom map has been determined for a given element topology, the host application may loop over all the elements of the given topology and form element stiffness, mass, force, etc. These functions typically require the host application to gather the global coordinates for each element as input. Symmetric element stiffness and consistent mass matrices are output as lower triangular matrices by rows. For example, a lower triangular stiffness matrix for an element with 4 degrees of freedom would be output in the order k11, k12, k22, k13, k23, k33, k14 ,k24, k34, k44.
```
  --               --
  | k11             |
  | k12 k22         |
  | k13 k23 k33     |
  | k14 k24 k34 k44 |
  --               --
```
Element force and diagonal mass are output as vectors. The host application is responsible for the memory required to hold the output stiffness, mass, force, etc. Note that the defined constant VFE_MAXNODE is equal to the maximum number of nodes used by any element formulation in VfeTools (about 30).
```
  Vdouble x[VFE_MAXNODE][3];
  Vdouble stiff[VFE_MAXDOF*(VFE_MAXDOF)/2], mdiag[VFE_MAXDOF];

              /* localize coordinates, x, from host application */
  .
  .
              /* generate element stiffness and diagonal mass */ 
  vfe_Solid2DStiff (solid2d,x,stiff);
  vfe_Solid2DMassDiag (solid2d,x mdiag);
```

Once computed, these quantities can be used by the host application in a variety of ways: the stiffness matrix can be assembled into a global stiffness matrix; loads can be assembled into a global load vector, etc. For those users that do not have a finite element host application, the VfsTools development library complements the functionality available in VfeTools to provide global system assembly, solution and degree of freedom management operations.

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.

In 3D analysis elements are defined in the x,y,z global Cartesian coordinate system.
In 2D plane stress and plane strain analysis the elements are defined in the x,y global Cartesian plane and the z axis must be a principal material direction.
In 2D axisymmetric analysis the elements are assumed to be bodies of revolution described in the cylindrical polar coordinates r,z,theta. The radial and axial coordinates are denoted by r and z respectively. At the theta = 0. plane, the r and z coordinates coincide with the x and y global Cartesian coordinates and the theta axis must be a principal material direction. All element properties are integrated through the entire body of revolution.

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:

Serendipity finite elements which are characterized by having corner nodes and nodes along element edges.
Lagrange finite elements which are complete products of Lagrange polynomials and are characterized by having corner nodes and nodes at uniform intervals in each element natural coordinate direction. Therefore quadratic Lagrange elements may have face centered and element centered nodes.

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.

Reaction Forces - The finite element equations are a discrete representation of a set of balance laws. In the case of structural mechanics, these laws translate into a set of element reaction forces that are conjugate to each nodal and/or elemental degree of freedom. As the displacement and rotation fields can be described by the nodal and elemental degrees of freedom, so can the strains, stresses, and the resulting reaction forces. As equilibrium sets in, the sum of all these reaction forces originating from adjacent elements at a node must equal the sum of the externally applied loads and the inertia loads. This is represented by the equation
```
  F_reaction = F_applied + F_inertia
```

Stiffness Matrix - In a linear simulation, the reaction forces F_reaction are computed given the stiffness matrix, K, for the element, and the vector of displacements, u, for the element. In a nonlinear analysis, the stiffness matrix is the negative of the derivative of the forces F_reaction with respect to the displacements. In the linear regime, the forces F_reaction, displacement u, and stiffness matrix K are related by the following expression:
```
  F_reaction = K u
```

Mass Matrix - The mass matrix M is used to compute the inertia loads from the element accelerations, a. VfeTools provides two distinct mass matrices: a lumped, or diagonal mass matrix, where a lumped mass is attached to each of the element degrees of freedom; and a consistent mass matrix, where the approximating field is used to compute a mass matrix that couples the various degrees of freedom in the element.
```
  F_inertia = - M a
```

Stress and Strain - the stresses and strains may be computed at the element nodes. In the case of shell and beam elements, depending upon the type of element formulation selected, the stress resultants and strains and curvatures are computed at element nodes. These quantities are provided along with the reaction forces at the degrees of freedom.

Geometric Stiffness - The geometric stiffness is a linear approximation to the changing stiffness if large deformations are used in the analysis. It can be shown that this linear approximation is proportional to the existing stress field, and in case of a compressive stress field, this additional stiffness can be used to estimate the linear buckling load of a structure.

Loads - Loads may be applied as concentrated loads or distributed loads. Concentrated loads are applied at specified nodes. Distributed loads may be applied throughout the volume or on the surface or edge of an element. The resulting loads at the element degrees of freedom can be calculated in a consistent manner according to the topology and element technology utilized.

Thermal Analysis - All of the above quantities have their counterpart in heat transfer simulations. Reaction forces become Power; the stiffness matrix becomes the Conductance Matrix; the mass becomes the Capacitance Matrix; the degrees of freedom u are now associated with the temperature; and stresses and strains are replaced by heat flux and temperature gradient, respectively. There is no counterpart to the geometric stiffness in a thermal analysis.

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:

SYS_DOF_TX - translation, tx, in the global x direction.
SYS_DOF_TY - translation, ty, in the global y direction.
SYS_DOF_TZ - translation, tz, in the global z direction
SYS_DOF_RX - rotation, rx, about the global x direction.
SYS_DOF_RY - rotation, ry, about the global y direction.
SYS_DOF_RZ - rotation, rz, about the global z direction.
SYS_DOF_TEMP - temperature, temp
SYS_DOF_VOLT - voltage, volt
SYS_DOF_WARP - warp, warp
SYS_DOF_PRES - pressure, pres
SYS_DOF_ELEM - internal elemental freedom
SYS_DOF_LAGM - Lagrange multiplier freedom
SYS_DOF_VX - vector component, vx, in the global x direction.
SYS_DOF_VY - vector component, vy, in the global y direction.
SYS_DOF_VX - vector component, vz, in the global z direction.
SYS_DOF_S0 - scalar 0, s0
SYS_DOF_S1 - scalar 1, s1
SYS_DOF_S2 - scalar 2, s2
SYS_DOF_S3 - scalar 3, s3
SYS_DOF_S4 - scalar 4, s4
SYS_DOF_S5 - scalar 5, s5
SYS_DOF_S6 - scalar 6, s6
SYS_DOF_S7 - scalar 7, s7

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:

Elastic properties, these properties are required for structural analysis and relate stress to strain.
Thermal properties, these properties are required for thermal analysis and relate heat flux to temperature gradient.
Density, the mass density is used to compute mass in structural analysis and, together with the specific heat, is used to compute heat capacitance in thermal analysis.
Specific heat, this property is used along with the density to compute heat capacitance in thermal analysis.
Coefficients of thermal expansion, these coefficients are used with the reference temperature to compute thermal strain in structural analysis.
Reference temperature, this property is used with the coefficients of thermal expansion to compute thermal strain in structural analysis.

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:

Shell properties, for shell element formulations. These properties include preintegrated material stiffness matrices for a given shell wall construction, general layered composite geometry and nodal thicknesses and offsets.
Beam properties, for beam and shell stiffener element formulations. These properties include preintegrated material stiffness matrices for a given beam section construction, general beam section geometry and nodal thicknesses and offsets.
Interface properties, for interface and contact element formulations. These properties include initial gaps, contact stiffnesses, film coefficients, etc.

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:

SGI:

cc -c -O2 -OPT:ieee_arith=3 -OPT:roundoff=3 -OPT:alias=restrict *.c

SUN:
```
cc -c -fast -xO3 -xrestrict *.c
```

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);
   }
}