1. Introduction
The CEETRON Solve Vfs library - also known as VfsTools - 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:
Global degree of freedom management.
Matrix assembly in sparse, full or diagonal formats.
Direct and iterative solution of large linear systems of equations. Solution techniques include direct sparse solvers and preconditioned conjugate gradient (PCG) and preconditioned generalized minimum residual (GMRES) iterative solvers.
Direct sparse solver supports real symmetric, real symmetric structure with non-symmetric coefficients, Hermitian, complex symmetric, and complex symmetric structure with non-symmetric coefficients.
Compute the generalized eigenvalue problem for large linear systems. Block Lanczos and AMLS (Automatic Multi-Level Substructuring) technigues are available. This enables the user to compute buckling loads and normal modes of vibration.
Perform static matrix condensation for large systems.
Perform general Ritz reductions of system matrices.
Compute rigid body modes
Object-based architecture, written in ANSI C and C++ with C++, FORTRAN, and C# language bindings.
Hardware and O/S independence.
All compute intensive operations are extensively parallelized using C++ threading and OpenMP.
Out-of-core options exist for direct, iterative and eigen solvers.
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.
Degree of Freedom Management
DofTab - Manage and number global degrees of freedom
System Vectors and Matrices
SysVector - System vector operations
SysMatrix - System matrix operations
Emulator
IPardiso - PARDISO-like interface
Analysis Procedures
SolProc - Time integration
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:
The element degrees of freedom must be registered to determine the degrees of freedom which are active in the global system. This information is used to count the number of active degrees of freedom and assign a number to each so that the individual element matrices and vectors can be assembled into the global system. The degree of freedom management module, DofTab is designed for this purpose. The C source code to create a DofTab object and declare it to maintain the degrees of freedom for a given number of nodes and elements, numnp and numel is as follows:
vfs_DofTab *doftab; doftab = vfs_DofTabBegin(); vfs_DofTabDef (doftab,numnp,numel);
The degree of freedom requirements of each element must now be registered with the DofTab object. This is done by first informing the DofTab using
vfs_DofTabSetMap()of the degree of freedom pattern or “map” of the type of element to be registered. Here nedofs is the number of freedoms used by the element and loc and tag are vectors containing the element node index and degree of freedom type for each degree of freedom. The VfeTools library contains functions to compute nedofs, loc and tag for all element types. Then the element connectivity of all elements which conform to the given map are registered in the DofTab usingvfs_DofTabElemUse(). The arguments tovfs_DofTabElemUse()are the element index, eid, the number of connected nodes, nix, and the list of connected nodes, ix. Note that the index space for all nodes and elements starts from one. Any permanently restrained degrees of freedom at nodes are specified usingvfs_DofTabNodePerm(). In the example below a single x translation at node 16 is restrained. If multipoint constraints (MPC) exist, they must be defined usingvfs_DofTabSetMPC(). Finally the global degrees of freedom are numbered usingvfs_DofTabProcess(). The total number of equations and reactions are returned.vfs_DofTabSetMap (doftab,nedofs,loc,tag); for(eid = 1; eid <= numel; eid++ { /* localize element connectivity ix[] from host application */ . . vfs_DofTabElemUse (doftab,eid,nix,ix); } tag[0] = SYS_DOF_TX; vfs_DofTabNodePerm (doftab,16,1,tag); vfs_DofTabProcess (doftab, &neq, &nre);
Now the host application must assemble and factorize the system matrix. A SysMatrix, system matrix is instanced and declared to be a symmetric sparse matrix of order neq. The DofTab object must be set into the SysMatrix object using
vfs_SysMatrixSetObject()so that the system matrix object can access the degree of freedom information and the structure of the sparse matrix. This is done with a call to thevfs_SysMatrixPreProcess()function. This function must be called before the matrix values are initialized usingvfs_SysMatrixInit()and the subsequent assembly of the matrix values can be performed. Assembly is performed using the DofTab object functionvfs_DofTabElemDof()and a subsequent call tovfs_SysMatrixAssem()in a loop over all elements to be assembled. Usevfs_DofTabElemDof()to gather the global dofs for each element and return them in integer vector, lm, for this purpose. Note that this process can be parallelized by the user.sysmatrix = vfs_SysMatrixBegin(); vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq); vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab); /* establish matrix structure */ vfs_SysMatrixPreProcess (sysmatrix); /* initialize values to 0. */ vfs_SysMatrixInit (sysmatrix,0.); /* assemble all elements */ for(eid = 1; eid <= numel; eid++) { /* localize element connectivity ix[] */ /* generate element stiffness matrix kt[] */ . . vfs_DofTabElemDof (doftab,1,&eid,nix,ix,nedofs,lm); vfs_SysMatrixAssem (sysmatrix,1,nedofs,lm,kt); }
If the matrix is to be factored a symbolic factorization must be performed using
vfs_SysMatrixProcess()before factorization withvfs_SysMatrixFactor().vfs_SysMatrixProcess (sysmatrix); vfs_SysMatrixFactor (sysmatrix);
The right hand side load vector is instanced and a load is applied at node 5 to the x and y translations. A solution vector is instanced to contain the displacement solution. The system is solved by calling
vfs_SysMatrixSolve()with the load and solution system vectors.load = vfs_SysVectorBegin (); vfs_SysVectorDef (load,neq); tag[0] = SYS_DOF_TX; tag[1] = SYS_DOF_TY; vfs_DofTabNodeDof (doftab,5,2,tag,lm); lval[0] = .707; lval[1] = .707; vfs_SysVectorAssem (load,1,2,lm,lval); soln = vfs_SysVectorBegin (); vfs_SysVectorDef (soln,neq); vfs_SysMatrixSolve (sysmatrix,load,soln);
Finally the x and y translational displacement values are gathered from the solution system vector for printing. The values are gathered into vector sval. The DofTab function
vfs_DofTabNodeDof()returns a degree of freedom number of zero in the lm vector if the requested node degree of freedom was restrained. The functionvfs_SysVectorGather()does not set a value into sval if the requested degree of freedom number is zero.printf("Displacements\n"); tag[0] = SYS_DOF_TX; tag[1] = SYS_DOF_TY; for(nid = 1; nid <= numnp; nid++) { vfs_DofTabNodeDof (doftab,nid,2,tag,lm); sval[0] = 0.; sval[1] = 0.; vfs_SysVectorGather (soln,2,lm,sval); printf(" nid= %d ux= %10f uy= %10f\n",nid,sval[0],sval[1]); }
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.
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.
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.
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.
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.
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.
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.
SysMatrix using
vfs_SysMatrixSetNumThreads()Direct and iterative linear equation solver.
Lanczos and AMLS eigensolution
Matrix reductionDofTab using
vfs_DofTabSetNumThreads()Equation reordering for sparse factorization.
Creation of matrix structure.
Coloring node/elements for parallel assembly
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.
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
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;
}
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
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);
}
}
