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.
DofTab Manage and number global degrees of freedom
SysVector System vector operations SysMatrix System matrix operations
IPardiso PARDISO-like interface
SolProc Time integration
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:
vfs_DofTab *doftab;
doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);
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);
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); }
vfs_SysMatrixProcess (sysmatrix); vfs_SysMatrixFactor (sysmatrix);
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);
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]);
}
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.
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.
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.Ax = b
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
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.A = L D L(T)
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.
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.comIn 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.frCurrently, 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.
Direct and iterative linear equation solver. Lanczos and AMLS eigensolution Matrix reduction
Equation reordering for sparse factorization. Creation of matrix structure. Coloring node/elements for parallel assembly
#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
#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
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
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);
}
}