Table of Contents

1. Introduction

The VfsTools library from Visual Kinematics, Inc. is an object-based software development toolkit designed to enable the solution of the systems of equations which typically occur in finite element simulations. VfsTools enables programmers to manage global system degrees of freedom, assemble and solve linear systems and compute eigenvalues and eigenvectors for large linear systems. The basic features of VfsTools are summarized below: Table of Contents

1.1 Module Summary

VfsTools is a set of modules that allow for the solution of large linear systems which result from finite element simulations. It contains two basic sets of modules: degree of freedom management and system vectors and matrices.

Each instance of a module is termed an object. Users can generate as many objects as needed. For example two SysMatrix objects would be required to hold a stiffness and a mass matrix. Several SysVector objects would be required to hold several computed normal mode of vibration eigenvectors.

The general methodology for solving a system of equations, for example in structural analysis, involves managing the global degrees of freedom, assembling the individual element matrices, such as element stiffness matrices, into the global system stiffness matrix, and performing a solution given a right hand side load vector. The procedure is as follows:

Table of Contents

1.2 Degree of Freedom Management

It is the case in finite element analysis that the large system matrix is assembled from a set of small element matrix contributions. A significant portion of the logic of performing this assembly is in the management of the assignment of global degree of freedom (dof) numbers. The DofTab module is designed to perform this function. The unique global freedoms at a node point, such as displacements, used by each element connected to the node is determined. To this set of nodal freedoms any internal elemental freedom must be added such as internal pressure freedoms, etc. The DofTab module processes the node connectivity and degree of freedom requirements of each element to assign global dof numbers to the unique freedoms which represent the assembled system.

Once the global dof have been assigned, the global dof numbers may then be queried for each element so that the element matrix contributions may be correctly assembled into the system matrix.

Table of Contents

1.3 System Vectors and Matrices

The two modules SysVector and SysMatrix are designed to manage and operate upon global degree of freedom (dof) system vectors and matrices. System vectors result, for example, from the assembly of individual finite element vectors, such as consistent applied load vectors. The SysVector module provides functions for assembly, query and a variety of arithmetic operations. System matrices result from the assembly of finite element matrices, such as stiffness matrices. System matrices may be symmetric, unsymmetric or diagonal and held in sparse or full data structures. A diagonal matrix is provided for diagonal matrices which result from assembling finite element diagonal mass matrices. Functions exist in the SysMatrix module for assembly, direct and iterative solution, eigenvalue and eigenvector extraction, etc.

The memory requirements of system vectors and system matrices is usually significant and is dependent upon the number of equations maintained and, in the case of the SysMatrix module, on the solution procedures performed.

The procedure for using the iterative solvers is much like that for the direct solver with the addition of the node coordinate and element topology and connectivity information. This information is formally entered into the DofTab object.

Table of Contents

1.4 Direct and Iterative Solution of a Linear System

The SysMatrix module solves a linear system of equations of the form:

Ax = b
where A is a symmetric or unsymmetric matrix and b is an input right hand side vector and x is an unknown vector. The vectors b and x are represented by SysVector objects. Two solution procedures are available, 1) direct solution and 2) iterative solution using a preconditioned conjugate gradient or preconditioned generalized minimum residual method.

The direct solution method is the most general solution procedure. It is applicable to indefinite matrices and the solution time is independent of the conditioning of the matrix A. The matrix is factorized into the form

A = L D L(T)
where L is a lower triangular matrix with unit diagonal, D is a diagonal matrix and (T) stands for the transpose. Direct solution can require a considerable amount of memory to hold the factorized matrix L. It is the method of choice when more than one right hand side vector is to be solved for or the matrix is indefinite.

Two iterative solution methods are available: 1) preconditioned conjugate gradient procedure which is applicable if the matrix is symmetric positive definite and 2) preconditioned generalized minimum residual procedure if the matrix is unsymmetric or indefinite. In general, iterative solution is effective if a small number of right hand sides is to be solved for and the problem size is relatively large. An advantage of the method is that the memory requirements are considerably less than the direct solution method. Iterative solution is generally applicable to static analysis of large 3D solid element models in which a small number of load cases are applied. A solution may be obtained to a large problem which can not be solved by a direct method because of available memory constraints.

The procedure for using the iterative solver is very similar to that for the direct solver with the addition of the specification of the node coordinates and element topology. This information is formally entered into the DofTab object.

Table of Contents

1.5 Compiling and Linking a VfsTools Application

To use VfsTools on a particular computer platform, the VfsTools source must be compiled and linked to an application. VfsTools is written in ANSI C and C++.

VfsTools is basically platform independent and as a result no user defined C preprocessor directives are required to compile VfsTools on any supported platform. For performance reasons VfsTools utilizes the BLAS for certain operations. In order to benefit from the speedup provided by the BLAS routines users are encouraged to utilize an external BLAS library. The currently supported external BLAS libraries are MKL from Intel and OpenBLAS. Note that if an external BLAS library is not used then certain VfsTools functionality will not be available, in particular the Block Lanczos and AMLS eigensolvers. Other solver technologies will have significantly degraded performance.

If MKL BLAS are enabled the modules in the base, vfs, and vfx directories must be compiled by defining either VKI_LIBAPI_BLASMKL_SEQUENTIAL or VKI_LIBAPI_BLASMKL_THREAD. If the 64-bit integer version of the MKL BLAS is used, the compiler flag VKI_LIBAPI_BLASMKL_ILP64 must be defined. Otherwise the 32-bit integer interface to the MKL BLAS will be assumed. If OpenBLAS are enabled then the base, vfs and vfx directories must be compiled with VKI_LIBAPI_OPENBLAS defined. The 32-bit or 64-bit version will be automatically detected during compilation from the OpenBLAS header file.

It is suggested to use the Intel BLAS on both Windows and LINUX platforms. See the following URL.

http://www.intel.com
In addition to the native direct sparse solver implementation in VfsTools, two external direct sparse solver packages are supported, Intel Pardiso and MUMPS. The use of these external solvers is not recommended. They have been included primarily for benchmarking purposes.

In order to enable the function calls to either of these packages VfsTools source must be compiled by defining VKI_LIBAPI_PARDISO and/or VKI_LIBAPI_MUMPS respectively. The Intel Pardiso package is included in the Intel MKL library as discussed above. MUMPS is written using MPI. In order to avoid the complications of full distributed memory parallelization, compile and link MUMPS using the sequential MPI library provided in the MUMPS distribution. The MUMPS package is available from the following URL.

http://mumps.enseeiht.fr
Currently, use of MUMPS and Pardiso is limited to symmetric matrices. In core and out of core versions are supported for both solvers. If out of core operation is to be used, or if the solution requires a large number of right hand side backsubstitutions, such as in eigensolution or the PCG iterative solver, it is suggested to use MUMPS. If in core operation is used with very few right hand side backsubstitutions, Pardiso is recommended.

Table of Contents

1.6 Parallelizing a VfsTools Application

A number of VfsTools modules support parallel processing. The current technology is designed to use threads in a shared memory environment. All the multi-threading technology is built upon C++ threading, OpenMP and threaded MKL BLAS. All parallelization done with OpenMP requires that VKI_LIBAPI_OPENMP is defined. All parallelization done with MKL BLAS requires that VKI_LIBAPI_BLASMKL_THREAD is defined. The SysMatrix and DofTab modules support internal parallelization using multi-threading. This form of parallelization is coarse grained.

It is possible to parallelize the element matrix assembly. Unlike other forms of parallelization that occur internally within a SysMatrix or DofTab function, assembly parallelization requires the host application be threaded and assemble groups of elements which have been computed by the DofTab object.

Table of Contents

1.7 A First Program - C Version

As an example of a simple VfsTools application the following program computes the deflection of two truss elements under a load. The resulting deflections are printed to standard output.
#include <stdio.h>
#include <math.h>
#include "vfs/vfs.h"

                   /* Stiffnesses (Lower Triangle) for Two Trusses */ 
static Vdouble kt[2][10] = {
   { 100.,
     0.,   0.,
    -100., 0., 100.,
     0.,   0.,   0.,   0. },
   { 0.,
     0., 100.,
     0.,   0.,   0.,
     0.,-100.,   0., 100. },
   };
                   /* Connectivity for Two Trusses */ 
static Vint ix[2][2] = {
   {1,2}, {2,3} };

                   /* Generic Truss Dof Map */
static Vint loc[4] = {
   1, 1, 2, 2 };
static Vint tag[4] = {
   SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TX, SYS_DOF_TY };

/*
                _ 10. load at 45 degrees
                /|
               /
     1 -/\/\- 2
    pin       |
              /
              \
              |
              3
             pin
*/

/*----------------------------------------------------------------------
                      Solve for Displacement of Two Trusses.
----------------------------------------------------------------------*/
int
main()
{
   vfs_SysMatrix *sysmatrix;
   vfs_SysVector *load, *soln;
   vfs_DofTab    *doftab;
   Vint numnp, numel;
   Vint nix, nndofs, nedofs;
   Vint eid, nid;
   Vint restag[2];
   Vint neq, nre;
   Vint lodtag[2];
   Vdouble lodval[2];
   Vint lm[4];
   Vint soltag[2];
   Vdouble solval[2];

                   /* Set number of nodes, elements */
   numnp = 3;
   numel = 2;
                   /* Set number of nodes per element */
   nix = 2;
                   /* Set dof per node, dof per element */
   nndofs = 2;
   nedofs = nix * nndofs;

                   /* Create degree of freedom table */ 
   doftab = vfs_DofTabBegin();
   vfs_DofTabDef (doftab,numnp,numel);

                   /* Loop over elements to set degree of freedom use */ 
   vfs_DofTabSetMap (doftab,nedofs,loc,tag);
   for(eid = 1; eid <= numel; eid++) {
      vfs_DofTabElemUse (doftab,eid,nix,ix[eid-1]);
   }

                   /* Suppress all translations at nodes 1 and 3 */
   restag[0] = SYS_DOF_TX;
   restag[1] = SYS_DOF_TY;
   nid = 1;
   vfs_DofTabNodePerm (doftab,nid,nndofs,restag);
   nid = 3;
   vfs_DofTabNodePerm (doftab,nid,nndofs,restag);

                   /* Process dof table, count equations and reactions */
   vfs_DofTabProcess (doftab, &neq, &nre);

                   /* Create system stiffness matrix */ 
   sysmatrix = vfs_SysMatrixBegin();
   vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
   vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab);

                   /* Initialize */ 
   vfs_SysMatrixPreProcess (sysmatrix);
   vfs_SysMatrixInit (sysmatrix,0.);
                   /* Assemble actual stiffness */ 
   for(eid = 1; eid <= numel; eid++) {
      vfs_DofTabElemDof (doftab,eid,nix,ix[eid-1],&nedofs,lm);
      vfs_SysMatrixAssem (sysmatrix,nedofs,lm,kt[eid-1]);
   }

                   /* Factorize matrix */ 
   vfs_SysMatrixProcess (sysmatrix);
   vfs_SysMatrixFactor (sysmatrix);

                   /* Create system vector for assembled load */ 
   load = vfs_SysVectorBegin ();
   vfs_SysVectorDef (load,neq);

                   /* Define load at node 2 */ 
   nid = 2;
   lodtag[0] = SYS_DOF_TX;
   lodtag[1] = SYS_DOF_TY;
   vfs_DofTabNodeDof (doftab,nid,nndofs,lodtag,lm);
   lodval[0] = .5*sqrt(2.)*10.;
   lodval[1] = .5*sqrt(2.)*10.;
   vfs_SysVectorAssem (load,nndofs,lm,lodval);

                   /* Create system vector for solution */ 
   soln = vfs_SysVectorBegin ();
   vfs_SysVectorDef (soln,neq);

                   /* Backsubstitution for solution */ 
   vfs_SysMatrixSolve (sysmatrix,load,soln);

                   /* Gather and print solution */ 
   printf("Displacements\n");
   soltag[0] = SYS_DOF_TX;
   soltag[1] = SYS_DOF_TY;
   for(nid = 1; nid <= numnp; nid++) {
      vfs_DofTabNodeDof (doftab,nid,nndofs,soltag,lm);
      solval[0] = 0.;
      solval[1] = 0.;
      vfs_SysVectorGather (soln,nndofs,lm,solval);
      printf(" nid= %d  ux= %10f  uy= %10f\n",nid,solval[0],solval[1]);
   }
   
                   /* Delete objects */ 
   vfs_DofTabEnd (doftab);
   vfs_SysVectorEnd (load);
   vfs_SysVectorEnd (soln);
   vfs_SysMatrixEnd (sysmatrix);
   return 0;
}

The output of this example program appears below.

Displacements
 nid= 1  ux=   0.000000  uy=   0.000000
 nid= 2  ux=   0.070711  uy=   0.070711
 nid= 3  ux=   0.000000  uy=   0.000000

Table of Contents

1.8 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 <math.h>
#include "base/base.h"
#include "vfs/vfs.h"

                   /* Stiffnesses (Lower Triangle) for Two Trusses */ 
static Vdouble kt[2][10] = {
   { 100.,
     0.,   0.,
    -100., 0., 100.,
     0.,   0.,   0.,   0. },
   { 0.,
     0., 100.,
     0.,   0.,   0.,
     0.,-100.,   0., 100. },
   };
                   /* Connectivity for Two Trusses */ 
static Vint ix[2][2] = {
   {1,2}, {2,3} };

                   /* Generic Truss Dof Map */
static Vint loc[4] = {
   1, 1, 2, 2 };
static Vint tag[4] = {
   SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TX, SYS_DOF_TY };

/*
                _ 10. load at 45 degrees
                /|
               /
     1 -/\/\- 2
    pin       |
              /
              \
              |
              3
             pin
*/

/*----------------------------------------------------------------------
                      Solve for Displacement of Two Trusses.
----------------------------------------------------------------------*/
int
main()
{
   vfs_SysMatrix *sysmatrix;
   vfs_SysVector *load, *soln;
   vfs_DofTab    *doftab;
   Vint numnp, numel;
   Vint nix, nndofs, nedofs;
   Vint eid, nid;
   Vint restag[2];
   Vint neq, nre;
   Vint lodtag[2];
   Vdouble lodval[2];
   Vint lm[4];
   Vint soltag[2];
   Vdouble solval[2];

                   /* Set number of nodes, elements */
   numnp = 3;
   numel = 2;
                   /* Set number of nodes per element */
   nix = 2;
                   /* Set dof per node, dof per element */
   nndofs = 2;
   nedofs = nix * nndofs;

                   /* Create degree of freedom table */ 
   doftab = new vfs_DofTab;
   doftab->Def (numnp,numel);

                   /* Loop over elements to set degree of freedom use */ 
   doftab->SetMap (nedofs,loc,tag);
   for(eid = 1; eid <= numel; eid++) {
      doftab->ElemUse (eid,nix,ix[eid-1]);
   }

                   /* Suppress all translations at nodes 1 and 3 */
   restag[0] = SYS_DOF_TX;
   restag[1] = SYS_DOF_TY;
   nid = 1;
   doftab->NodePerm (nid,nndofs,restag);
   nid = 3;
   doftab->NodePerm (nid,nndofs,restag);

                   /* Process dof table, count equations and reactions */
   doftab->Process (&neq, &nre);

                   /* Create system stiffness matrix */ 
   sysmatrix = new vfs_SysMatrix;
   sysmatrix->Def (SYSMATRIX_SYMM_SPARSE,neq);
   sysmatrix->SetObject (VFS_DOFTAB,(Vobject*)doftab);

                   /* Initialize */ 
   sysmatrix->PreProcess ();
   sysmatrix->Init (0.);

                   /* Assemble actual stiffness */ 
   for(eid = 1; eid <= numel; eid++) {
      doftab->ElemDof (eid,nix,ix[eid-1],&nedofs,lm);
      sysmatrix->Assem (nedofs,lm,kt[eid-1]);
   }

                   /* Factorize matrix */ 
   sysmatrix->Process ();
   sysmatrix->Factor ();

                   /* Create system vector for assembled load */ 
   load = new vfs_SysVector;
   load->Def (neq);

                   /* Define load at node 2 */ 
   nid = 2;
   lodtag[0] = SYS_DOF_TX;
   lodtag[1] = SYS_DOF_TY;
   doftab->NodeDof (nid,nndofs,lodtag,lm);
   lodval[0] = .5*sqrt(2.)*10.;
   lodval[1] = .5*sqrt(2.)*10.;
   load->Assem (nndofs,lm,lodval);

                   /* Create system vector for solution */ 
   soln = new vfs_SysVector;
   soln->Def (neq);

                   /* Backsubstitution for solution */ 
   sysmatrix->Solve (load,soln);

                   /* Gather and print solution */ 
   printf("Displacements\n");
   soltag[0] = SYS_DOF_TX;
   soltag[1] = SYS_DOF_TY;
   for(nid = 1; nid <= numnp; nid++) {
      doftab->NodeDof (nid,nndofs,soltag,lm);
      solval[0] = 0.;
      solval[1] = 0.;
      soln->Gather (nndofs,lm,solval);
      printf(" nid= %d  ux= %10f  uy= %10f\n",nid,solval[0],solval[1]);
   }
   
                   /* Delete objects */ 
   delete doftab;
   delete load;
   delete soln;
   delete sysmatrix;
   return 0;
}
Table of Contents

1.9 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 vfs_ is replaced by vfsf_.

C 
C               _ 10. load at 45 degrees
C               /|
C              /
C    1 -/\/\- 2
C   pin       |
C             /
C             \
C             |
C             3
C            pin
C 

C-----------------------------------------------------------------------
C                   Solve for Displacement of Two Trusses.
C-----------------------------------------------------------------------
      program main
      include 'base/fortran/base.inc'
      include 'vfs/fortran/vfs.inc'
C                  Stiffnesses (Lower Triangle) for Two Trusses
      double precision kt(10,2)
      data kt /
     &  100.,
     &    0.,   0.,
     & -100., 0., 100.,
     &    0.,   0.,   0.,   0.,
     &    0.,
     &    0., 100.,
     &    0.,   0.,   0.,
     &    0.,-100.,   0., 100. /
C                  Connectivity for Two Trusses
      integer ix(2,2)
      data ix / 1,2, 2,3 /
C                  Generic Truss Dof Map
      integer loc(4)
      data loc / 1,1, 2,2 /
      integer tag(4)
      data tag / SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TX, SYS_DOF_TY /

      double precision sysmatrix
      double precision load, soln
      double precision doftab
      integer numnp, numel
      integer nix, nndofs, nedofs
      integer eid, nid
      integer restag(2)
      integer neq, nre
      integer lodtag(2)
      double precision lodval(2)
      integer lm(4)
      integer soltag(2)
      double precision solval(2)

C                  Set number of nodes, elements
      numnp = 3
      numel = 2
C                  Set number of nodes per element
      nix = 2
C                  Set dof per node, dof per element
      nndofs = 2
      nedofs = nix * nndofs

C                  Create degree of freedom table
      call vfsf_DofTabBegin (doftab)
      call vfsf_DofTabDef (doftab,numnp,numel)

C                  Loop over elements to set degree of freedom use
      call vfsf_DofTabSetMap (doftab,nedofs,loc,tag)
      do eid = 1,numel
         call vfsf_DofTabElemUse (doftab,eid,nix,ix(1,eid))
      enddo

C                  Suppress all translations at nodes 1 and 3
      restag(1) = SYS_DOF_TX
      restag(2) = SYS_DOF_TY
      nid = 1
      call vfsf_DofTabNodePerm (doftab,nid,nndofs,restag)
      nid = 3
      call vfsf_DofTabNodePerm (doftab,nid,nndofs,restag)

C                  Process dof table, count equations and reactions
      call vfsf_DofTabProcess (doftab,neq,nre)

C                  Create system stiffness matrix
      call vfsf_SysMatrixBegin (sysmatrix)
      call vfsf_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq)

C                  Preprocess stiffness matrix
      call vfsf_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab)
      call vfsf_SysMatrixPreProcess (sysmatrix)

C                  Initialize
      call vfsf_SysMatrixInit (sysmatrix,0.)

C                  Assemble actual stiffness
      do eid = 1,numel
         call vfsf_DofTabElemDof (doftab,eid,nix,ix(1,eid),nedofs,lm)
         call vfsf_SysMatrixAssem (sysmatrix,nedofs,lm,kt(1,eid))
      enddo

C                  Factorize matrix
      call vfsf_SysMatrixProcess (sysmatrix)
      call vfsf_SysMatrixFactor (sysmatrix)

C                  Create system vector for assembled load
      call vfsf_SysVectorBegin (load)
      call vfsf_SysVectorDef (load,neq)

C                  Define load at node 2
      nid = 2
      lodtag(1) = SYS_DOF_TX
      lodtag(2) = SYS_DOF_TY
      call vfsf_DofTabNodeDof (doftab,nid,nndofs,lodtag,lm)
      lodval(1) = .5*sqrt(2.)*10.
      lodval(2) = .5*sqrt(2.)*10.
      call vfsf_SysVectorAssem (load,nndofs,lm,lodval)

C                  Create system vector for solution
      call vfsf_SysVectorBegin (soln)
      call vfsf_SysVectorDef (soln,neq)

C                  Backsubstitution for solution
      call vfsf_SysMatrixSolve (sysmatrix,load,soln)

C                  Gather and print solution
      write(6,*) "Displacements"
      soltag(1) = SYS_DOF_TX
      soltag(2) = SYS_DOF_TY
      do nid = 1,numnp
         call vfsf_DofTabNodeDof (doftab,nid,nndofs,soltag,lm)
         solval(1) = 0.
         solval(2) = 0.
         call vfsf_SysVectorGather (soln,nndofs,lm,solval)
         write(6,1000) nid,solval(1),solval(2)
      enddo
   
C                  Delete objects
      call vfsf_DofTabEnd (doftab)
      call vfsf_SysVectorEnd (load)
      call vfsf_SysVectorEnd (soln)
      call vfsf_SysMatrixEnd (sysmatrix)
1000  format(' nid= ',i2,'  ux= ',f11.6,'  uy= ',f11.6)
      stop
      end
Table of Contents

1.10 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 vfs_ is replaced by vfs..

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

public class intro1 {

                   /* Stiffnesses (Lower Triangle) for Two Trusses */ 
   public static double [][] kt = {
      new double[] { 100.0,
                     0.0,   0.0,
                    -100.0, 0.0, 100.0,
                     0.0,   0.0,   0.0,   0.0 },
      new double[] { 0.0,
                     0.0, 100.0,
                     0.0,   0.0,   0.0,
                     0.0,-100.0,   0.0, 100.0 }
   };
                   /* Connectivity for Two Trusses */ 
   public static int [][] ix = {
      new int[] {1,2},
      new int[] {2,3}
   };

                   /* Generic Truss Dof Map */
   public static int [] loc = { 1, 1, 2, 2 };

   public static int [] tag = {
      vsy.SYS_DOF_TX, vsy.SYS_DOF_TY, vsy.SYS_DOF_TX, vsy.SYS_DOF_TY };

/*
                _ 10. load at 45 degrees
                /|
               /
     1 -/\/\- 2
    pin       |
              /
              \
              |
              3
             pin
*/

/*----------------------------------------------------------------------
                      Solve for Displacement of Two Trusses.
----------------------------------------------------------------------*/
   public static void Main() {

      IntPtr sysmatrix;
      IntPtr load, soln;
      IntPtr doftab;
      int numnp = 0, numel = 0;
      int nix = 0, nndofs = 0, nedofs = 0;
      int eid = 0, nid = 0;
      int [] restag = new int[2];
      int neq = 0, nre = 0;
      int [] lodtag = new int[2];
      double [] lodval = new double[2];
      int [] lm = new int[4];
      int [] soltag = new int[2];
      double [] solval = new double[2];

                   /* Set number of nodes, elements */
      numnp = 3;
      numel = 2;
                   /* Set number of nodes per element */
      nix = 2;
                   /* Set dof per node, dof per element */
      nndofs = 2;
      nedofs = nix * nndofs;

                   /* Create degree of freedom table */ 
      doftab = vfs.DofTabBegin();
      vfs.DofTabDef (doftab,numnp,numel);

                   /* Loop over elements to set degree of freedom use */ 
      vfs.DofTabSetMap (doftab,nedofs,loc,tag);
      for(eid = 1; eid <= numel; eid++) {
         vfs.DofTabElemUse (doftab,eid,nix,ix[eid-1]);
      }

                   /* Suppress all translations at nodes 1 and 3 */
      restag[0] = vsy.SYS_DOF_TX;
      restag[1] = vsy.SYS_DOF_TY;
      nid = 1;
      vfs.DofTabNodePerm (doftab,nid,nndofs,restag);
      nid = 3;
      vfs.DofTabNodePerm (doftab,nid,nndofs,restag);

                   /* Process dof table, count equations and reactions */
      vfs.DofTabProcess (doftab,ref neq,ref nre);

                   /* Create system stiffness matrix */ 
      sysmatrix = vfs.SysMatrixBegin();
      vfs.SysMatrixDef (sysmatrix,vfs.SYSMATRIX_SYMM_SPARSE,neq);
      vfs.SysMatrixSetObject (sysmatrix,vfs.VFS_DOFTAB,doftab);

                   /* PreProcess, to form matrix structure */ 
      vfs.SysMatrixPreProcess (sysmatrix);

                   /* Initialize */ 
      vfs.SysMatrixInit (sysmatrix,0.0);

                   /* Assemble actual stiffness */ 
      for(eid = 1; eid <= numel; eid++) {
         vfs.DofTabElemDof (doftab,eid,nix,ix[eid-1],ref nedofs,lm);
         vfs.SysMatrixAssem (sysmatrix,nedofs,lm,kt[eid-1]);
      }

                   /* Factorize matrix */ 
      vfs.SysMatrixProcess (sysmatrix);
      vfs.SysMatrixFactor (sysmatrix);

                   /* Create system vector for assembled load */ 
      load = vfs.SysVectorBegin ();
      vfs.SysVectorDef (load,neq);

                   /* Define load at node 2 */ 
      nid = 2;
      lodtag[0] = vsy.SYS_DOF_TX;
      lodtag[1] = vsy.SYS_DOF_TY;
      vfs.DofTabNodeDof (doftab,nid,nndofs,lodtag,lm);
      lodval[0] = 0.5*Math.Sqrt(2.0)*10.0;
      lodval[1] = 0.5*Math.Sqrt(2.0)*10.0;
      vfs.SysVectorAssem (load,nndofs,lm,lodval);

                   /* Create system vector for solution */ 
      soln = vfs.SysVectorBegin ();
      vfs.SysVectorDef (soln,neq);

                   /* Backsubstitution for solution */ 
      vfs.SysMatrixSolve (sysmatrix,load,soln);

                   /* Gather and print solution */ 
      Console.Write("Displacements\n");
      soltag[0] = vsy.SYS_DOF_TX;
      soltag[1] = vsy.SYS_DOF_TY;
      for(nid = 1; nid <= numnp; nid++) {
         vfs.DofTabNodeDof (doftab,nid,nndofs,soltag,lm);
         solval[0] = 0.0;
         solval[1] = 0.0;
         vfs.SysVectorGather (soln,nndofs,lm,solval);
         Console.Write(" nid= {0}  ux= {1}  uy= {2}\n",nid,solval[0],solval[1]);
      }
   
                   /* Delete objects */ 
      vfs.DofTabEnd (doftab);
      vfs.SysVectorEnd (load);
      vfs.SysVectorEnd (soln);
      vfs.SysMatrixEnd (sysmatrix);
   }
}