A First Vfs Program

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.

CC++FORTRANC#

A First Program - C Version

#include <stdio.h>
#include <math.h>
#include "sam/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