6. Examples
This chapter contains example VfsTools programs. The source code for all examples resides in devtools/src/vfs/exam.
6.1. Example 1, Solve a System of Linear Equations
This example illustrates solving a small symmetric system of linear equations. The order of the matrix is 6, ie the number of global degrees of freedom. Four objects are required; a SysMatrix object to hold the symmetric matrix, two SysVector objects to hold the right hand side vector and solution vector and a DofTab object to manage the assignment and processing of the global degrees of freedom.
The DofTab object is documented with reference to “nodes” and “elements”.
In this case of assembling a given global matrix structure, each entry
in the matrix is an element and each row is a node.
The matrix structure is establish in DofTab by entering each off diagonal
entry as a 2x2 element and each diagonal entry as a 1x1 element.
The function vfs_DofTabProcess() will determine the matrix structure
and global degree of freedom ordering. It returns the number of global
degrees of freedom and number of reactions.
The symmetric matrix is denoted by A, and is represented in a row, column
index format. A SysMatrix object is instanced and defined to hold
a symmetric sparse matrix. The order, ie. number of global degrees of freedom,
is returned as neq by vfs_DofTabProcess() - in this case 6.
The DofTab object must be set in the SysMatrix object
using vfs_SysMatrixSetObject(), and vfs_SysMatrixPreProcess() called
to access the matrix structure from DofTab.
The matrix values are initialized to zero
using vfs_SysMatrixInit() and the actual matrix values are assembled
using vfs_SysMatrixAssem().
The matrix is then factored by calling vfs_SysMatrixFactor(). The default
solution procedure is a direct solver so a full factorization is performed
by this function.
The right hand side vector, denoted by B, is entered into a SysVector
object using vfs_SysVectorScatter(). Use vfs_DofTabNodeDof() to find the
global degree of freedom number assigned to each row.
A solution system vector is instanced and the system is solved by calling
vfs_SysMatrixSolve() with the right hand side and solution system vectors
as arguments. The entries in the solution vector are queried using
vfs_SysVectorGather() and printed.
#include <stdio.h>
#include "vfs/vfs.h"
/*
A Matrix B Vector
4 0 1 1 0 0 11
0 4 1 0 1 0 16
1 1 4 0 1 0 20
1 0 0 4 1 0 22
0 1 1 1 4 0 29
0 0 0 0 0 4 24
*/
static Vdouble aval[12] = {
4.0, 1.0, 1.0, 4.0, 1.0, 1.0,
4.0, 1.0, 4.0, 1.0, 4.0, 4.0 };
static Vint icol[12] = {
1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6 };
static Vint irow[12] = {
1, 3, 4, 2, 3, 5, 3, 5, 4, 5, 5, 6 };
static Vint loc[2] = {1,2};
static Vint tag[2] = {SYS_DOF_TEMP,SYS_DOF_TEMP};
/* B Vector */
static Vdouble bval[6] = {
11.0, 16.0, 20.0, 22.0, 29.0, 24.0 };
/*----------------------------------------------------------------------
Solve a System of Linear Equations
----------------------------------------------------------------------*/
int
main()
{
Vint i;
vfs_SysMatrix *amat;
vfs_SysVector *bvec, *xvec;
vfs_DofTab *doftab;
Vint nedofs, nix;
Vint neq, nre;
Vint nnequn, nnzero;
Vint lm[2], ix[2];
Vdouble vm[3];
Vdouble xval;
/* Set degree of freedom table */
nnequn = 6;
nnzero = 12;
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,nnequn,nnzero);
/* Specify active temperature degrees of freedom */
for(i = 0; i < nnzero; ++i) {
if(icol[i] == irow[i]) {
nedofs = 1;
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
ix[0] = icol[i];
nix = 1;
vfs_DofTabElemUse (doftab,i+1,nix,ix);
} else {
nedofs = 2;
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
ix[0] = icol[i];
ix[1] = irow[i];
nix = 2;
vfs_DofTabElemUse (doftab,i+1,nix,ix);
}
}
/* Process degrees of freedom */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create system stiffness matrix */
amat = vfs_SysMatrixBegin();
vfs_SysMatrixDef (amat,SYSMATRIX_SYMM_SPARSE,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixSetObject (amat,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (amat);
/* Initialize */
vfs_SysMatrixInit (amat,0.);
/* Assemble actual values */
for(i = 0; i < nnzero; ++i) {
if(irow[i] == icol[i]) {
ix[0] = icol[i];
nix = 1;
vfs_DofTabElemDof (doftab,i+1,nix,ix,&nedofs,lm);
vm[0] = aval[i];
} else {
ix[0] = icol[i];
ix[1] = irow[i];
nix = 2;
vfs_DofTabElemDof (doftab,i+1,nix,ix,&nedofs,lm);
vm[0] = 0.;
vm[1] = aval[i];
vm[2] = 0.;
}
vfs_SysMatrixAssem (amat,nedofs,lm,vm);
}
/* Symbolic factorization */
vfs_SysMatrixProcess (amat);
/* Factorize matrix */
vfs_SysMatrixFactor (amat);
/* Create system vector for assembled load */
bvec = vfs_SysVectorBegin ();
vfs_SysVectorDef (bvec,neq);
/* Set load vector */
for(i = 1; i <= neq; ++i) {
vfs_DofTabNodeDof (doftab,i,1,tag,lm);
vfs_SysVectorScatter (bvec,1,lm,&bval[i-1]);
}
/* Create system vector for solution */
xvec = vfs_SysVectorBegin ();
vfs_SysVectorDef (xvec,neq);
/* Backsubstitution for solution */
vfs_SysMatrixSolve (amat,bvec,xvec);
/* Gather and print solution */
printf("Solution vector\n");
for(i = 1; i <= neq; i++) {
vfs_DofTabNodeDof (doftab,i,1,tag,lm);
vfs_SysVectorGather (xvec,1,lm,&xval);
printf("x[%d]= %e\n",i,xval);
}
/* Delete objects */
vfs_SysVectorEnd (bvec);
vfs_SysVectorEnd (xvec);
vfs_SysMatrixEnd (amat);
vfs_DofTabEnd (doftab);
return 0;
}
The output of this example program appears below.
Solution vector
x[1]= 1.000000e+00
x[2]= 2.000000e+00
x[3]= 3.000000e+00
x[4]= 4.000000e+00
x[5]= 5.000000e+00
x[6]= 6.000000e+00
6.2. Example 1a, Solve for Vibration Eigenvalues and Eigenvectors
This example illustrates solving a “vibration” eigenproblem.
The stiffness matrix is sparse symmetric and the mass matrix is diagonal.
The stiffness matrix is assembled as in Example 1 and processed
using vfs_SysMatrixProcess() since it is the matrix to be factored
during the solution of the eigenproblem. The diagonal mass matrix
is assembled using vfs_SysMatrixAssemDiag(). The DofTab object
does not need to be set in a diagonal SysMatrix object.
The 4 lowest vibration modes are requested
using the Lanczos eigensolver. The lower bound of the interval is 0.
and the upper bound is unset and assumed to be infinity.
The actual eigensolution occurs with the call to vfs_SysMatrixEigen().
The sysmatrix amat is factored and an iteration is performed to solve
the requested eigenproblem. Upon completion the number of eigenvalues
computed and the eigenvalues themselves are returned. The eigenvectors
themselves are retrieved in SysVector objects associated with
each individual eigenvalue using the function vfs_SysMatrixEigenMode().
After all desired eigenvectors are retrieved, the internal storage
for the eigenvectors should be released using vfs_SysMatrixEigenFree().
The eigenvalues and eigenvectors are then printed.
#include <stdio.h>
#include "vfs/vfs.h"
/*
A Matrix B Matrix
4 0 1 1 0 0 10
0 4 1 0 1 0 12
1 1 4 0 1 0 14
1 0 0 4 1 0 16
0 1 1 1 4 0 18
0 0 0 0 0 4 20
*/
/* A matrix */
static Vdouble aval[12] = {
4.0, 1.0, 1.0, 4.0, 1.0, 1.0,
4.0, 1.0, 4.0, 1.0, 4.0, 4.0 };
static Vint icol[12] = {
1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6 };
static Vint irow[12] = {
1, 3, 4, 2, 3, 5, 3, 5, 4, 5, 5, 6 };
static Vint loc[2] = {1,2};
static Vint tag[2] = {SYS_DOF_TEMP,SYS_DOF_TEMP};
/* B matrix */
static Vdouble bval[6] = {
10.0, 12.0, 14.0, 16.0, 18.0, 20.0 };
/*----------------------------------------------------------------------
Solve for Vibration Eigenvalues and Eigenvectors
----------------------------------------------------------------------*/
int
main()
{
Vint i, j;
vfs_SysMatrix *amat, *bmat;
vfs_DofTab *doftab;
Vint nedofs, nix;
Vint neq, nre, neigen;
Vint nnzero;
Vint lm[2], ix[2];
Vdouble vm[3];
Vdouble ev[4];
vfs_SysVector *evec[4];
Vdouble xval;
/* Set degree of freedom table */
nnzero = 12;
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,nnzero,nnzero);
/* Specify active temperature degrees of freedom */
for(i = 0; i < nnzero; ++i) {
if(icol[i] == irow[i]) {
nedofs = 1;
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
ix[0] = icol[i];
nix = 1;
vfs_DofTabElemUse (doftab,i+1,nix,ix);
} else {
nedofs = 2;
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
ix[0] = icol[i];
ix[1] = irow[i];
nix = 2;
vfs_DofTabElemUse (doftab,i+1,nix,ix);
}
}
/* Process degrees of freedom */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create system stiffness matrix */
amat = vfs_SysMatrixBegin();
vfs_SysMatrixDef (amat,SYSMATRIX_SYMM_SPARSE,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixSetObject (amat,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (amat);
/* Initialize */
vfs_SysMatrixInit (amat,0.);
/* Assemble actual values */
for(i = 0; i < nnzero; ++i) {
if(irow[i] == icol[i]) {
ix[0] = icol[i];
nix = 1;
vfs_DofTabElemDof (doftab,i+1,nix,ix,&nedofs,lm);
vm[0] = aval[i];
} else {
ix[0] = icol[i];
ix[1] = irow[i];
nix = 2;
vfs_DofTabElemDof (doftab,i+1,nix,ix,&nedofs,lm);
vm[0] = 0.;
vm[1] = aval[i];
vm[2] = 0.;
}
vfs_SysMatrixAssem (amat,nedofs,lm,vm);
}
/* Symbolic factorization */
vfs_SysMatrixProcess (amat);
/* Create system diagonal mass matrix */
bmat = vfs_SysMatrixBegin();
vfs_SysMatrixDef (bmat,SYSMATRIX_DIAG,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixPreProcess (amat);
/* Assemble actual values */
for(i = 0; i < neq; ++i) {
vfs_DofTabNodeDof (doftab,i+1,1,tag,lm);
nedofs = 1;
vm[0] = bval[i];
vfs_SysMatrixAssemDiag (bmat,nedofs,lm,vm);
}
/* set eigen problem to 4 lowest eigenvalues */
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN,SYSMATRIX_EIGEN_VIBE);
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN_METHOD,SYS_EIGEN_LANCZOS);
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN_NUM,4);
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN_TYPE,SYS_EIGEN_LOWEST);
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN_LOWERTEST,SYS_ON);
vfs_SysMatrixSetParamd (amat,SYSMATRIX_EIGEN_LOWER,0.);
vfs_SysMatrixSetParami (amat,SYSMATRIX_EIGEN_UPPERTEST,SYS_OFF);
/* Solve eigenproblem */
vfs_SysMatrixEigen (amat,bmat,&neigen,ev);
/* Create system vectors eigenvectors */
for(i = 0; i < neigen; i++) {
evec[i] = vfs_SysVectorBegin ();
vfs_SysVectorDef (evec[i],neq);
vfs_SysMatrixEigenMode (amat,i+1,evec[i]);
}
/* free internal eigenmode storage */
vfs_SysMatrixEigenFree (amat);
/* Gather and print eigenvectors */
for(i = 0; i < neigen; i++) {
printf("Eigenvector %d, eigenvalue= %e\n",i+1,ev[i]);
for(j = 1; j <= neq; j++) {
vfs_SysVectorGather (evec[i],1,&j,&xval);
printf("x[%d]= %e\n",j,xval);
}
}
/* Delete objects */
for(i = 0; i < neigen; i++) {
vfs_SysVectorEnd (evec[i]);
}
vfs_SysMatrixEnd (amat);
vfs_SysMatrixEnd (bmat);
vfs_DofTabEnd (doftab);
return 0;
}
6.3. Example 2, Solve for Unknowns of Four Elements
This example illustrates assembling four identical finite element matrices together and solving the system for a given load vector. The element matrix is held in array ke in lower triangular form. The global degrees of freedom for each element are held in array lm. The elements are meant to be interpreted as 4 node quadrilateral elements with a single degree of freedom at each node.
#include <stdio.h>
#include "base/base.h"
#include "vfs/vfs.h"
/* Element Stiffness (Lower Triangle) */
static Vdouble ke[10] = {
4.0,
-1.0, 4.0,
-1.0, -1.0, 4.0,
-1.0, -1.0, -1.0, 4.0 };
/* Global Degrees of Freedom for Four Elements
assumed to be connected as follows:
7 --- 8 --- 9
| 3 | 4 |
4 --- 5 --- 6
| 1 | 2 |
1 --- 2 --- 3
*/
static Vint ix[4][4] = {
{ 1, 2, 5, 4 },
{ 2, 3, 6, 5 },
{ 4, 5, 8, 7 },
{ 5, 6, 9, 8 } };
/* Load vector for Each Global Degree of Freedom */
static Vdouble lval[9] = {
-7.0, -8.0, -1.0, 4.0, 20.0, 16.0, 11.0, 28.0, 17.0 };
static Vint loc[4] = {1,2,3,4};
static Vint tag[4] = {SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP};
/*----------------------------------------------------------------------
Solve for Unknowns of Four Elements
----------------------------------------------------------------------*/
int
main()
{
Vint i;
vfs_DofTab *doftab;
vfs_SysMatrix *sysmatrix;
vfs_SysVector *load, *soln;
Vint numel, numnp;
Vint nedofs;
Vint neq, nre, nix, idof, lm[4];
Vdouble xval;
/* Set number of nodes and elements */
numnp = 9;
numel = 4;
/* Set degree of freedom table */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
/* Set number of degrees of freedom/nodes per element */
nedofs = 4;
nix = 4;
/* Specify constant element map */
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
/* Specify active temperature degrees of freedom */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemUse (doftab,i+1,nix,ix[i]);
}
/* Process degrees of freedom */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create system stiffness matrix */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize */
vfs_SysMatrixInit (sysmatrix,0.);
/* Assemble actual stiffness */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemDof (doftab,i+1,nix,ix[i],&nedofs,lm);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}
/* Factorize matrix */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixFactor (sysmatrix);
/* Create system vector for load */
load = vfs_SysVectorBegin ();
vfs_SysVectorDef (load,neq);
vfs_SysVectorInit (load,0.);
/* Set load vector */
for(i = 1; i <= numnp; ++i) {
vfs_DofTabNodeDof (doftab,i,1,tag,&idof);
vfs_SysVectorAssem (load,1,&idof,&lval[i-1]);
}
/* 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("Solution vector\n");
for(i = 1; i <= numnp; i++) {
vfs_DofTabNodeDof (doftab,i,1,tag,&idof);
vfs_SysVectorGather (soln,1,&idof,&xval);
printf("x[%d]= %e\n",i,xval);
}
/* Delete objects */
vfs_SysVectorEnd (load);
vfs_SysVectorEnd (soln);
vfs_SysMatrixEnd (sysmatrix);
vfs_DofTabEnd (doftab);
return 0;
}
The output of this example program appears below.
Solution vector
x[1]= 1.000000e+00
x[2]= 2.000000e+00
x[3]= 3.000000e+00
x[4]= 4.000000e+00
x[5]= 5.000000e+00
x[6]= 6.000000e+00
x[7]= 7.000000e+00
x[8]= 8.000000e+00
x[9]= 9.000000e+00
6.4. Example 2a, Solve for Unknowns of Four Elements with MPC
This example illustrates applying a multipoint constraint (MPC)
to the finite element model of Example 2.
The MPC equation is .5*x8 + .5*x5 - x6 = 0., the solution at node 6
is constrained to be the average of the solution at node 5 and 8.
The MPC equation must be specified using vfs_DofTabSetMPC() before
the function vfs_DofTabProcess() is called.
#include <stdio.h>
#include "vfs/vfs.h"
/* Element Stiffness (Lower Triangle) */
static Vdouble ke[10] = {
4.0,
-1.0, 4.0,
-1.0, -1.0, 4.0,
-1.0, -1.0, -1.0, 4.0 };
/* Global Degrees of Freedom for Four Elements
assumed to be connected as follows:
7 --- 8 --- 9
| 3 | 4 |
4 --- 5 --- 6
| 1 | 2 |
1 --- 2 --- 3
*/
static Vint ix[4][4] = {
{ 1, 2, 5, 4 },
{ 2, 3, 6, 5 },
{ 4, 5, 8, 7 },
{ 5, 6, 9, 8 } };
/* Load vector for Each Global Degree of Freedom */
static Vdouble lval[9] = {
-7.0, -8.0, -1.0, 4.0, 20.0, 16.0, 11.0, 28.0, 17.0 };
static Vint loc[4] = {1,2,3,4};
static Vint tag[4] = {SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP};
/* MPC dofs and coefficients */
/* .5*x8 + .5*x5 - x6 = 0 */
static Vint mpcix[3] = { 8, 5, 6 };
static Vdouble mpccoef[3] = {.5,.5,-1.};
/*----------------------------------------------------------------------
Solve for Unknowns of Four Elements with MPC
----------------------------------------------------------------------*/
int
main()
{
Vint i;
vfs_DofTab *doftab;
vfs_SysMatrix *sysmatrix;
vfs_SysVector *load, *soln;
Vint numel, numnp;
Vint nedofs;
Vint neq, nre, nix, idof, lm[4];
Vdouble xval;
/* Set number of nodes and elements */
numnp = 9;
numel = 4;
/* Set degree of freedom table */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
/* Set number of degrees of freedom/nodes per element */
nedofs = 4;
nix = 4;
/* Specify constant element map */
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
/* Specify active temperature degrees of freedom */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemUse (doftab,i+1,nix,ix[i]);
}
/* Set MPC */
vfs_DofTabSetMPC (doftab,3,mpcix,tag,mpccoef,0.);
/* Process degrees of freedom */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create system stiffness matrix */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize */
vfs_SysMatrixInit (sysmatrix,0.);
/* Assemble actual stiffness */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemDof (doftab,i+1,nix,ix[i],&nedofs,lm);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}
/* Factorize matrix */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixFactor (sysmatrix);
/* Create system vector for load */
load = vfs_SysVectorBegin ();
vfs_SysVectorDef (load,neq);
vfs_SysVectorInit (load,0.);
/* Set load vector */
for(i = 1; i <= numnp; ++i) {
vfs_DofTabNodeDof (doftab,i,1,tag,&idof);
vfs_SysVectorAssem (load,1,&idof,&lval[i-1]);
}
/* 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("Solution vector\n");
for(i = 1; i <= numnp; i++) {
vfs_DofTabNodeDof (doftab,i,1,tag,&idof);
vfs_SysVectorGather (soln,1,&idof,&xval);
printf("x[%d]= %e\n",i,xval);
}
/* Delete objects */
vfs_SysVectorEnd (load);
vfs_SysVectorEnd (soln);
vfs_SysMatrixEnd (sysmatrix);
vfs_DofTabEnd (doftab);
return 0;
}
6.5. Example 2b, Multiple Solutions of Four Elements with Restraints
This example illustrates explicitly restraining a degree of freedom to be zero and solving for multiple load vectors. The element stiffness matrices are characteristic of a thermal conduction matrix, with a single degree of freedom (the temperature) at each node. The temperature at node 3 is restrained to zero and the global degree of freedom numbering reflects this restraint. Note that restrained freedoms are flagged with a 0 in the global degree of freedom vector for an element.
Two load cases are solved, first a concentrated flux at node 7 and
second a distributed flux on elements 1 and 3.
The first case illustrates setting the loads at each degree of freedom
into the system vector using vfs_SysVectorScatter(). The second
case illustrates assembling loads using vfs_SysVectorAssem().
The difference being that vfs_SysVectorScatter() will set the values
while vfs_SysVectorAssem() will sum the values to existing values
in SysVector.
The two load cases are solved using vfs_SysMatrixSolveMult(). They
could also be solved individually using vfs_SysMatrixSolve().
The solution values are queried and printed using vfs_SysVectorGather().
Note that the outofcoreflag variable can be set to
perform an out-of-core direct factorization and solution with SysMatrix.
The default high performance multi-frontal solver is used.
#include <stdio.h>
#include "vfs/vfs.h"
/* 4-node quad thermal stiffness matrix (Lower Triangle) */
static Vdouble ke[10] = {
0.666667,
-0.166667, 0.666667,
-0.333333, -0.166667, 0.666667,
-0.166667, -0.333333, -0.166667, 0.666667 };
/* Node connectivity for four elements,
connectivity is counterclockwise from lower left.
7 --- 8 --- 9
| 3 | 4 |
4 --- 5 --- 6
| 1 | 2 |
1 --- 2 --- 3
Assume temperature at node 3 is restrained to zero,
*/
/* Global degree of freedom vectors for each element */
static Vint ix[4][4] = {
{ 1, 2, 5, 4 },
{ 2, 3, 6, 5 },
{ 4, 5, 8, 7 },
{ 5, 6, 9, 8 } };
static Vint loc[4] = {1,2,3,4};
static Vint tag[4] = {SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP};
/* Load vectors for each global degree of freedom */
/* Case 1, concentrated load at node 7 (degree of freedom 6) */
static Vdouble cval1 = 100.;
/* Case 2, distributed loads for left face elements 1 and 3 */
static Vdouble dval2[4] = { 50., 0., 0., 50. };
/*----------------------------------------------------------------------
Multiple Solutions of Four Elements with Restraints
----------------------------------------------------------------------*/
int
main()
{
Vint i;
vfs_DofTab *doftab;
vfs_SysMatrix *sysmatrix;
vfs_SysVector *load[2], *soln[2];
Vint numnp, numel, nix;
Vint nedofs, idof;
Vint neq, nre, lm[4];
Vdouble xval1, xval2;
Vint factorized;
Vlong factornonz;
Vint outofcoreflag;
/* Set number of nodes and elements */
numnp = 9;
numel = 4;
/* Set degree of freedom table */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
/* Set number of degrees of freedom/nodes per element */
nedofs = 4;
nix = 4;
/* Specify constant element map */
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
/* Specify active temperature degrees of freedom */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemUse (doftab,i+1,nix,ix[i]);
}
/* Create permanent constraint */
vfs_DofTabNodePerm (doftab,3,1,tag);
/* Process degrees of freedom */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create system stiffness matrix */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
/* PreProcess, to determine matrix structure */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize */
vfs_SysMatrixInit (sysmatrix,0.);
/* Assemble actual stiffness */
for(i = 0; i < numel; ++i) {
vfs_DofTabElemDof (doftab,i+1,nix,ix[i],&nedofs,lm);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}
/* set out-of-core options */
outofcoreflag = 1;
if(outofcoreflag) {
vfs_SysMatrixSetParami (sysmatrix,SYSMATRIX_FACTOR_OOC,SYS_ON);
/* set file to be used for temporary storage of factor */
vfs_SysMatrixSetOocFile (sysmatrix,SYSMATRIX_KFACTOR_FILE,"exam2b.ooc");
}
/* Process matrix, symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);
/* Factorize matrix */
vfs_SysMatrixFactor (sysmatrix);
/* check for factorization failure */
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_FACTORIZED,&factorized);
if(factorized == 0) {
printf("Factorization failure\n");
exit(0);
}
/* Query and print size of factor */
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_FACTOR_NONZ,&factornonz);
printf("Factorization complete, number of non zeros= %lld\n",factornonz);
/* Create system vectors for loads */
load[0] = vfs_SysVectorBegin ();
load[1] = vfs_SysVectorBegin ();
vfs_SysVectorDef (load[0],neq);
vfs_SysVectorDef (load[1],neq);
/* Set case 1 load vector */
vfs_SysVectorInit (load[0],0.);
vfs_DofTabNodeDof (doftab,7,1,tag,&idof);
vfs_SysVectorScatter (load[0],1,&idof,&cval1);
/* Set case 2 load vector */
vfs_SysVectorInit (load[1],0.);
vfs_DofTabElemDof (doftab,1,nix,ix[0],&nedofs,lm);
vfs_SysVectorAssem (load[1],nedofs,lm,dval2);
vfs_DofTabElemDof (doftab,3,nix,ix[2],&nedofs,lm);
vfs_SysVectorAssem (load[1],nedofs,lm,dval2);
/* Create system vectors for solutions */
soln[0] = vfs_SysVectorBegin ();
soln[1] = vfs_SysVectorBegin ();
vfs_SysVectorDef (soln[0],neq);
vfs_SysVectorDef (soln[1],neq);
/* Backsubstitution for solutions */
vfs_SysMatrixSolveMult (sysmatrix,2,load,soln);
/* Gather and print solutions */
printf("Solution vector\n");
for(i = 1; i <= numnp; i++) {
vfs_DofTabNodeDof (doftab,i,1,tag,lm);
if(lm[0]) {
vfs_SysVectorGather (soln[0],1,lm,&xval1);
vfs_SysVectorGather (soln[1],1,lm,&xval2);
printf("x1[%d]= %e, x2[%d]= %e\n",i,xval1,i,xval2);
} else {
printf("x1[%d]= %e, x2[%d]= %e\n",i,0.,i,0.);
}
}
/* Delete objects */
vfs_SysVectorEnd (load[0]);
vfs_SysVectorEnd (load[1]);
vfs_SysVectorEnd (soln[0]);
vfs_SysVectorEnd (soln[1]);
vfs_SysMatrixEnd (sysmatrix);
vfs_DofTabEnd (doftab);
return 0;
}
6.6. Example 3, Direct and Iterative Solver
This example illustrates all steps required to assemble a set of finite element matrices, and solve the resulting linear system with both the direct or iterative methods. For simplicity, the element matrix originates from a 4-noded quadrilateral element whose values are pre-defined as ke in the example file. The total number of element degrees of freedom for this problem is 4, with one degree of freedom per node, SYS_DOF_TEMP.
The toggle variable iterative_solver can be used to specify either the iterative solver (1), or the direct solver (0).
The mesh size is controlled by the parameter N. There are N x N elements and (N+1) x (N+1) nodes. A single boundary condition on the first node is specified as a temporary restraint, and its value is set to 2.0. A unit concentrated heat flow is applied on the last node.
The example also uses an instance of the DofTab object. This is required both to number the active equations, establish the form of the sparse matrix, as well as a helper object to form the preconditioner for the iterative solver.
If the iterative solver is to be used, the DofTab object must be
informed of that fact using
the function vfs_DofTabSetParami() with type DOFTAB_REDUCE
set to DOFTAB_REDUCE_ITER.
In addition, it needs to be augmented with the nodal coordinates
and element topologies for all elements. Coordinates are added with
vfs_DofTabSetCoords(), while element topologies are set with
vfs_DofTabSetTopology().
The iterative solver also requires that DofTab have advance knowledge of all
nodes where a temporary restraint is being defined in SysMatrix.
This is done with vfs_DofTabNodeDofSpec() or vfs_DofTabNodeSpec().
Once DofTab has been processed, the number of equations is used to define the SysMatrix object.
This example also illustrates an easy methodology to enforce non-zero boundary conditions. The parameter SYSMATRIX_SOLVETEMP, when enabled, indicates that initial values in the solution vector where temporary restraints are applied indicate will be used to compute an additional internal right-hand-side vector that is computed from these values. If not enabled, it is the user’s responsibility to add these loads to the SysVector object defined the right-hand-side vector.
Once again, additional information can be specified for the iterative solver in SysMatrix: the residual convergence tolerances - SYSMATRIX_ITER_UTOL and SYSMATRIX_ITER_FTOL. The last two default to 5.e-03 and 5.e-07, respectively.
Once the matrix has been factorized, the number of non-zero entries in the
matrix and in the factor are retrieved with vfs_SysMatrixGetLong(), using
the types SYSMATRIX_ASSEM_NONZ and SYSMATRIX_FACTOR_NONZ, respectively.
The total memory usage for the solution of the linear system is retrieved
with vfs_SysMatrixGetInteger() and the type SYSMATRIX_MEMORY. Memory
returned is in megabytes.
After the system has been solved, if the iterative method has been used,
the number of iterations is retrieved with vfs_SysMatrixGetInteger(), using
the type SYSMATRIX_NUMITER.
Finally, the solution on the first and last nodes are retrieved and printed.
#include <stdio.h>
#include "vfs/vfs.h"
/* 4-node quad thermal stiffness matrix */
static Vdouble ke[10] = {
0.666667,
-0.166667, 0.666667,
-0.333333, -0.166667, 0.666667,
-0.166667, -0.333333, -0.166667, 0.666667 };
/* element dof map arrays for 4-node quad thermal element */
static Vint nedofs = 4;
static Vint loc[4] = { 1, 2, 3, 4 };
static Vint tag[4] = { SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP };
/* topology information */
static Vint shape = SYS_SHAPEQUAD;
static Vint maxi = 2;
static Vint maxj = 0;
static Vint maxk = 0;
static Vint nix = 4;
#define MAX_DOF_ELEMENT 4
#define MAX_NODE_ELEMENT 4
/* flag for iterative solver */
static Vint iterative_solver = 1;
/* mesh size parameter */
#define N 50
/* load value */
static Vdouble one = 1.;
/* boundary condition value */
static Vdouble two = 2.;
/* dof for solution */
static Vint dof = SYS_DOF_TEMP;
/*----------------------------------------------------------------------
compute element connectivity in simple grid
----------------------------------------------------------------------*/
static void
get_elementconnectivity(Vint n, Vint ix[])
{
Vint i, j;
j = (n-1)/N+1;
i = n - (j-1)*N;
ix[0] = (N+1)*(j-1) + i;
ix[1] = ix[0]+1;
ix[2] = ix[1]+N+1;
ix[3] = ix[2]-1;
}
/*----------------------------------------------------------------------
Direct and Iterative Solver
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n, numiter, memory;
Vlong anonz, fnonz;
/* declare objects */
vfs_SysMatrix *sysmatrix;
vfs_SysVector *sysvectorb, *sysvectorx;
vfs_DofTab *doftab;
/* number of nodes, elements */
Vint numnp, numel;
/* number of equations and reactions */
Vint neq, nre;
/* number of degrees of freedom in element,
element degree of freedom numbers */
Vint lm[MAX_DOF_ELEMENT];
/* number of nodes in element and connectivity */
Vint ix[MAX_NODE_ELEMENT];
/* node coordinates */
Vdouble x[3];
/* solution value */
Vdouble value;
/* compute number of nodes and elements */
numnp = (N+1)*(N+1);
numel = N*N;
doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);
/* compute and set node coordinates */
for(n = 1; n <= numnp; n++) {
j = (n-1)/(N+1)+1;
i = n - (j-1)*(N+1);
x[0] = (Vdouble)(i-1);
x[1] = (Vdouble)(j-1);
x[2] = 0.;
vfs_DofTabSetCoords (doftab,n,x);
}
}
/* gather and set element topology,
degree of freedom type and connectivity */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,n,nix,ix);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetTopology (doftab,n,shape,maxi,maxj,maxk);
}
}
/* register boundary condition node */
if(iterative_solver) {
ix[0] = 1;
vfs_DofTabNodeSpec (doftab,1,1,ix);
}
/* process DofTab to determine number of equations */
vfs_DofTabProcess (doftab,&neq,&nre);
/* create sysmatrix object */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixSetParami (sysmatrix,SYSMATRIX_SOLVETEMP,SYS_ON);
if(iterative_solver) {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_ITER,neq);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_UTOL,5.e-7);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_FTOL,5.e-3);
} else {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
}
/* register DofTab with SysMatrix object */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab);
/* create load and displacement sysvector objects */
sysvectorb = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorb,neq);
sysvectorx = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorx,neq);
/* preprocess sysmatrix */
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize matrix values to zero */
vfs_SysMatrixInit (sysmatrix,0.);
/* assembly for each element */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}
/* set boundary condition value and temporary restraints */
vfs_SysVectorInit (sysvectorx,0.);
vfs_DofTabNodeDof (doftab,1,1,&dof,lm);
vfs_SysVectorScatter (sysvectorx,1,lm,&two);
vfs_SysMatrixSetTemp (sysmatrix,1,lm);
/* add true point load on last node */
vfs_SysVectorInit (sysvectorb,0.);
vfs_DofTabNodeDof (doftab,numnp,1,&dof,lm);
vfs_SysVectorAssem (sysvectorb,1,lm,&one);
/* equation ordering and/or symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_ASSEM_NONZ,&anonz);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_FACTOR_NONZ,&fnonz);
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_MEMORY,&memory);
printf("Assembled non-zeros= %lld\n",anonz);
printf("Factor non-zeros= %lld\n",fnonz);
printf("Total memory= %dMB\n",memory);
/* factorize for direct solution */
/* form preconditioner for iterative solution */
vfs_SysMatrixFactor (sysmatrix);
/* solve */
vfs_SysMatrixSolve (sysmatrix,sysvectorb,sysvectorx);
if(iterative_solver) {
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_NUMITER,&numiter);
printf("Number of iterations= %d\n",numiter);
}
/* access displacement vector using vfs_SysVectorGather */
printf ("Solution at first and last node:\n");
n = 1;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
n = numnp;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
/* delete objects */
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (sysmatrix);
vfs_SysVectorEnd (sysvectorb);
vfs_SysVectorEnd (sysvectorx);
return 0;
}
The output of this example program if the direct method is used appears below.
Assembled non-zeros= 12701
Factor non-zeros= 65032
Total memory= 0MB
Solution at first and last node:
Node= 1: 2.0000e+000
Node= 2601: 9.3942e+000
If the iterative method is used, the output is given by
Assembled non-zeros= 22700
Factor non-zeros= 4750
Total memory= 0MB
Number of iterations= 9
Solution at first and last node:
Node= 1: 2.0000e+000
Node= 2601: 9.3942e+000
6.7. Example 3a, Direct and Iterative Solver with Monitor Function
This example illustrates the use of a monitor function to capture data
generated by the SysMatrix as it executes either a vfs_SysMatrixFactor()
or a vfs_SysMatrixSolve() operation. The example is an extension of the
previous example, and only the additions to the previous example will be
highlighted.
SysMatrix calls a user-defined monitor function periodically. During a direct factorization, the monitor function is called after about every 1000 equations have been processed. This gives the user an opportunity to abort the factorization in case the estimated completion time exceeds the user’s expectations.
A monitor function can be set with vfs_SysMatrixSetFunction(). In addition
to the SysMatrix object, this function takes an enumerated function type
- in this case, SYSMATRIX_FUN_MONITOR - the monitor function pointer,
and a pointer to a user object. This user object is then made available to the
user every time the monitor function is called, and can be used to pass
information back and forth between SysMatrix and the calling application.
The monitor function itself, written by the user, has two arguments: SysMatrix and the user-defined pointer, cast to (Vobject*). The user should cast this pointer inside the monitor function to reveal the structure of the user data.
If the iterative solver is used, the example demonstrates how a monitor
function can be used to retrieve an iteration number, and the residual and
solution increment convergence values at that iteration. Note that these
values are only available while vfs_SysMatrixSolve() is executing.
If the direct solver is used, the example retrieves the equation number currently being processed by the factorization and prints this value to the screen.
The application is responsible for the user data, including memory allocation and deallocation.
#include <stdio.h>
#include "vfs/vfs.h"
/* 4-node quad thermal stiffness matrix */
static Vdouble ke[10] = {
0.666667,
-0.166667, 0.666667,
-0.333333, -0.166667, 0.666667,
-0.166667, -0.333333, -0.166667, 0.666667 };
/* element dof map arrays for 4-node quad thermal element */
static Vint nedofs = 4;
static Vint loc[4] = { 1, 2, 3, 4 };
static Vint tag[4] = { SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP };
/* topology information */
static Vint shape = SYS_SHAPEQUAD;
static Vint maxi = 2;
static Vint maxj = 0;
static Vint maxk = 0;
static Vint nix = 4;
#define MAX_DOF_ELEMENT 4
#define MAX_NODE_ELEMENT 4
/* flag for iterative solver */
static Vint iterative_solver = 1;
/* mesh size parameter */
#define N 50
/* load value */
static Vdouble one = 1.;
/* boundary condition value */
static Vdouble two = 2.;
/* dof for solution */
static Vint dof = SYS_DOF_TEMP;
/*
* user struct for monitor function
* whereami == 0 -> not in SysMatrixFactor
* whereami == 1 -> in SysMatrixFactor
*/
typedef struct {
Vint whereami;
} UserObj;
/*----------------------------------------------------------------------
monitoring function
----------------------------------------------------------------------*/
static void
monitor(vfs_SysMatrix *p, Vobject *obj)
{
UserObj *user = (UserObj*)obj;
Vint eqn, numiter, phase;
Vdouble uconv, fconv, econv;
if(iterative_solver == 0 && user->whereami == 1) {
/* retrieve and print factorization equation number */
vfs_SysMatrixGetInteger (p,SYSMATRIX_FACTOR_EQN,&eqn);
printf("...Factorization at equation= %10d\n",eqn);
} else if(iterative_solver == 1 && user->whereami == 2) {
vfs_SysMatrixGetInteger (p,SYSMATRIX_PHASE,&phase);
if(phase == SYSMATRIX_PHASE_SOLVE_ITER) {
vfs_SysMatrixGetInteger (p,SYSMATRIX_NUMITER,&numiter);
vfs_SysMatrixGetDouble (p,SYSMATRIX_ITER_FCONV,&fconv);
vfs_SysMatrixGetDouble (p,SYSMATRIX_ITER_UCONV,&uconv);
vfs_SysMatrixGetDouble (p,SYSMATRIX_ITER_ECONV,&econv);
printf("Iter %3d: Res conv= %11.4e, Sol conv= %11.4e, Ene conv= %11.4e\n",
numiter,fconv,uconv,econv);
}
}
}
/*----------------------------------------------------------------------
compute element connectivity in simple grid
----------------------------------------------------------------------*/
static void
get_elementconnectivity(Vint n, Vint ix[])
{
Vint i, j;
j = (n-1)/N+1;
i = n - (j-1)*N;
ix[0] = (N+1)*(j-1) + i;
ix[1] = ix[0]+1;
ix[2] = ix[1]+N+1;
ix[3] = ix[2]-1;
}
/*----------------------------------------------------------------------
Direct and Iterative Solver with Monitor Function
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n, numiter, memory;
Vlong anonz, fnonz;
Vdouble fconv, uconv, econv;
UserObj user;
/* declare objects */
vfs_SysMatrix *sysmatrix;
vfs_SysVector *sysvectorb, *sysvectorx;
vfs_DofTab *doftab;
/* number of nodes, elements */
Vint numnp, numel;
/* number of equations and reactions */
Vint neq, nre;
/* number of degrees of freedom in element,
element degree of freedom numbers */
Vint lm[MAX_DOF_ELEMENT];
/* number of nodes in element and connectivity */
Vint ix[MAX_NODE_ELEMENT];
/* node coordinates */
Vdouble x[3];
/* solution value */
Vdouble value;
/* compute number of nodes and elements */
numnp = (N+1)*(N+1);
numel = N*N;
doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);
/* gather and set node coordinates */
for(n = 1; n <= numnp; n++) {
j = (n-1)/(N+1)+1;
i = n - (j-1)*(N+1);
x[0] = (Vdouble)(i-1);
x[1] = (Vdouble)(j-1);
x[2] = 0.;
vfs_DofTabSetCoords (doftab,n,x);
}
}
/* gather and set element topology,
degree of freedom type and connectivity */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,n,nix,ix);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetTopology (doftab,n,shape,maxi,maxj,maxk);
}
}
/* register boundary condition node */
if(iterative_solver) {
ix[0] = 1;
vfs_DofTabNodeSpec (doftab,1,1,ix);
}
/* process DofTab to determine number of equations */
vfs_DofTabProcess (doftab,&neq,&nre);
/* create sysmatrix object */
sysmatrix = vfs_SysMatrixBegin();
/* set monitor function */
vfs_SysMatrixSetFunction (sysmatrix,SYSMATRIX_FUN_MONITOR,monitor,
(Vobject**)&user);
vfs_SysMatrixSetParami (sysmatrix,SYSMATRIX_SOLVETEMP,SYS_ON);
user.whereami = 0;
if(iterative_solver) {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_ITER,neq);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_UTOL,5.e-7);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_FTOL,5.e-3);
} else {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
}
/* register DofTab with SysMatrix object */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab);
/* create load and displacement sysvector objects */
sysvectorb = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorb,neq);
sysvectorx = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorx,neq);
/* preprocess sysmatrix */
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize matrix values to zero */
vfs_SysMatrixInit (sysmatrix,0.);
/* assembly for each element */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ke);
}
/* set boundary condition value */
vfs_SysVectorInit (sysvectorx,0.);
vfs_DofTabNodeDof (doftab,1,1,&dof,lm);
vfs_SysVectorScatter (sysvectorx,1,lm,&two);
/* set temporary restraints */
vfs_SysMatrixSetTemp (sysmatrix,1,lm);
/* initialize vector to zero */
vfs_SysVectorInit (sysvectorb,0.);
/* add point load on last node */
vfs_DofTabNodeDof (doftab,numnp,1,&dof,lm);
vfs_SysVectorAssem (sysvectorb,1,lm,&one);
/* equation ordering and/or symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_ASSEM_NONZ,&anonz);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_FACTOR_NONZ,&fnonz);
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_MEMORY,&memory);
printf("Assembled non-zeros= %lld\n",anonz);
printf("Factor non-zeros= %lld\n",fnonz);
printf("Total memory= %dMB\n",memory);
/* factorize for direct solution */
/* form preconditioner for iterative solution */
user.whereami = 1;
vfs_SysMatrixFactor (sysmatrix);
/* solve */
user.whereami = 2;
vfs_SysMatrixSolve (sysmatrix,sysvectorb,sysvectorx);
if(iterative_solver) {
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_NUMITER,&numiter);
vfs_SysMatrixGetDouble (sysmatrix,SYSMATRIX_ITER_FCONV,&fconv);
vfs_SysMatrixGetDouble (sysmatrix,SYSMATRIX_ITER_UCONV,&uconv);
vfs_SysMatrixGetDouble (sysmatrix,SYSMATRIX_ITER_ECONV,&econv);
printf("Number of iterations= %d\n",numiter);
printf("Final residual convergence= %11.4e\n",fconv);
printf("Final solution convergence= %11.4e\n",uconv);
printf("Final energy convergence= %11.4e\n",econv);
}
/* access displacement vector using vfs_SysVectorGather */
printf ("Solution at first and last node:\n");
n = 1;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
n = numnp;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
/* delete objects */
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (sysmatrix);
vfs_SysVectorEnd (sysvectorb);
vfs_SysVectorEnd (sysvectorx);
return 0;
}
The output of this example program if the direct method is used appears below.
Assembled non-zeros= 12701
Factor non-zeros= 71201
Total memory= 0MB
...Factorization at equation= 1005
...Factorization at equation= 2006
Solution at first and last node:
Node= 1: 2.0000e+000
Node= 2601: 9.3942e+000
If the iterative method is used, the output is given by
Assembled non-zeros= 22700
Factor non-zeros= 4750
Total memory= 0MB
Iter 0: Res conv= 1.0000e+00, Sol conv= 1.0000e+00, Ene conv= 1.0000e+00
Iter 1: Res conv= 4.1746e-01, Sol conv= 3.7973e-01, Ene conv= 1.0000e+00
Iter 2: Res conv= 1.2697e-01, Sol conv= 3.6108e-02, Ene conv= 1.0944e-01
Iter 3: Res conv= 6.5288e-02, Sol conv= 5.7108e-04, Ene conv= 1.6905e-02
Iter 4: Res conv= 2.0303e-02, Sol conv= 4.0489e-03, Ene conv= 3.0851e-03
Iter 5: Res conv= 8.9219e-03, Sol conv= 2.0185e-04, Ene conv= 4.3417e-04
Iter 6: Res conv= 3.0466e-03, Sol conv= 5.0965e-05, Ene conv= 6.0284e-05
Iter 7: Res conv= 1.3660e-03, Sol conv= 5.6594e-05, Ene conv= 1.0571e-05
Iter 8: Res conv= 4.3549e-04, Sol conv= 1.9956e-05, Ene conv= 1.5217e-06
Iter 9: Res conv= 2.0234e-04, Sol conv= 1.7523e-05, Ene conv= 2.0992e-07
Number of iterations= 9
Final residual convergence= 2.0234e-04
Final solution convergence= 1.7523e-05
Final energy convergence= 2.0992e-07
Solution at first and last node:
Node= 1: 2.0000e+00
Node= 2601: 9.3942e+00
6.8. Example 3b, Direct and Iterative Solver with Parallel Assembly
This example illustrates the assembling the model in Example 3a
in parallel using OpenMP. The DofTab object here is
used to compute the groups of elements to be assembled in parallel.
Use the function vfs_DofTabSetParami() with type DOFTAB_THREADASSEM
set to SYS_ON to enable the generation of the
assembly groups when vfs_DofTabProcess() is called.
This example hardwires the maximum number of OpenMP threads to 8.
The SysMatrix object must be setup for parallel assembly
by calling vfs_SysMatrixSetNumThreadsAssem() with
the maximum number of threads.
Use vfs_DofTabThreadAssemNum() and vfs_DofTabThreadAssemGroups()
to get the assembly group information from DofTab.
To assemble, enter a double loop, the outer loop over the element groups
and the inner loop over the element in each group. Note
the use of omp_get_thread_num() to get the thread number to
be passed to vfs_SysMatrixThreadAssem().
#include <stdio.h>
#include <omp.h>
#include "vfs/vfs.h"
/* 4-node quad thermal stiffness matrix */
static Vdouble ke[10] = {
0.666667,
-0.166667, 0.666667,
-0.333333, -0.166667, 0.666667,
-0.166667, -0.333333, -0.166667, 0.666667 };
/* element dof map arrays for 4-node quad thermal element */
static Vint nedofs = 4;
static Vint loc[4] = { 1, 2, 3, 4 };
static Vint tag[4] = { SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP, SYS_DOF_TEMP };
/* topology information */
static Vint shape = SYS_SHAPEQUAD;
static Vint maxi = 2;
static Vint maxj = 0;
static Vint maxk = 0;
static Vint nix = 4;
#define MAX_DOF_ELEMENT 4
#define MAX_NODE_ELEMENT 4
/* flag for iterative solver */
static Vint iterative_solver = 0;
/* mesh size parameter */
#define N 50
/* load value */
static Vdouble one = 1.;
/* boundary condition value */
static Vdouble two = 2.;
/* dof for solution */
static Vint dof = SYS_DOF_TEMP;
/* number of OpenMP threads to use during assembly */
static Vint nomp = 8;
/*----------------------------------------------------------------------
compute element connectivity in simple grid
----------------------------------------------------------------------*/
static void
get_elementconnectivity(Vint n, Vint ix[])
{
Vint i, j;
j = (n-1)/N+1;
i = n - (j-1)*N;
ix[0] = (N+1)*(j-1) + i;
ix[1] = ix[0]+1;
ix[2] = ix[1]+N+1;
ix[3] = ix[2]-1;
}
/*----------------------------------------------------------------------
Direct and Iterative Solver with Parallel Assembly
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n, numiter, memory, num;
Vint numgroups, *numelem, *blockid, *color, in, ig, istart, bid;
Vlong anonz, fnonz;
/* declare objects */
vfs_SysMatrix *sysmatrix;
vfs_SysVector *sysvectorb, *sysvectorx;
vfs_DofTab *doftab;
/* number of nodes, elements */
Vint numnp, numel;
/* number of equations and reactions */
Vint neq, nre;
/* number of degrees of freedom in element,
element degree of freedom numbers */
Vint lm[MAX_DOF_ELEMENT];
/* number of nodes in element and connectivity */
Vint ix[MAX_NODE_ELEMENT];
/* node coordinates */
Vdouble x[3];
/* solution value */
Vdouble value;
/* compute number of nodes and elements */
numnp = (N+1)*(N+1);
numel = N*N;
doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);
/* flag to compute element colors */
vfs_DofTabSetParami (doftab,DOFTAB_THREADASSEM,SYS_ON);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);
/* gather and set node coordinates */
for(n = 1; n <= numnp; n++) {
j = (n-1)/(N+1)+1;
i = n - (j-1)*(N+1);
x[0] = (Vdouble)(i-1);
x[1] = (Vdouble)(j-1);
x[2] = 0.;
vfs_DofTabSetCoords (doftab,n,x);
}
}
/* assume all elements belong to blockid 1 */
vfs_DofTabSetAssemBlock (doftab,1);
/* gather and set element topology,
degree of freedom type and connectivity */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,n,nix,ix);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetTopology (doftab,n,shape,maxi,maxj,maxk);
}
}
/* register boundary condition node */
if(iterative_solver) {
ix[0] = 1;
vfs_DofTabNodeSpec (doftab,1,1,ix);
}
/* process DofTab to determine number of equations */
vfs_DofTabProcess (doftab,&neq,&nre);
/* create sysmatrix object */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixSetParami (sysmatrix,SYSMATRIX_SOLVETEMP,SYS_ON);
if(iterative_solver) {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_ITER,neq);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_UTOL,5.e-7);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_FTOL,5.e-3);
} else {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
}
vfs_SysMatrixSetNumThreadsAssem (sysmatrix,nomp);
/* register DofTab with SysMatrix object */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab);
/* create load and displacement sysvector objects */
sysvectorb = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorb,neq);
sysvectorx = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorx,neq);
/* preprocess sysmatrix */
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize matrix values to zero */
vfs_SysMatrixInit (sysmatrix,0.);
/* flag that element assembly will be done in parallel */
vfs_DofTabThreadAssemNum (doftab,&numgroups);
/* retrieve coloring information */
numelem = (Vint*)malloc(numgroups*sizeof(Vint));
blockid = (Vint*)malloc(numgroups*sizeof(Vint));
color = (Vint*)malloc(numel*sizeof(Vint));
vfs_DofTabThreadAssemGroups (doftab,&numgroups,numelem,blockid,color);
/* assembly for each element */
istart = 0;
for(ig = 0; ig < numgroups; ++ig) {
Vint numthreads;
num = numelem[ig];
bid = blockid[ig];
numthreads = nomp;
if(bid < 0) {
numthreads = 1;
}
#pragma omp parallel for num_threads(numthreads)
for(in = 0; in < num; in++) {
Vint _n, _ix[4], _lm[4], _nedofs;
_n = color[istart+in];
get_elementconnectivity (_n,_ix);
vfs_DofTabElemDof (doftab,_n,nix,_ix,&_nedofs,_lm);
vfs_SysMatrixThreadAssem (sysmatrix,omp_get_thread_num()+1,
_nedofs,_lm,ke);
}
/* move to next group of elements */
istart += num;
}
free (numelem);
free (blockid);
free (color);
/* set boundary condition value and temporary restraints */
vfs_SysVectorInit (sysvectorx,0.);
vfs_DofTabNodeDof (doftab,1,1,&dof,lm);
vfs_SysVectorScatter (sysvectorx,1,lm,&two);
/* set temporary restraints */
vfs_SysMatrixSetTemp (sysmatrix,1,lm);
/* initialize vector to zero */
vfs_SysVectorInit (sysvectorb,0.);
/* add point load on last node */
vfs_DofTabNodeDof (doftab,numnp,1,&dof,lm);
vfs_SysVectorAssem (sysvectorb,1,lm,&one);
/* equation ordering and/or symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_ASSEM_NONZ,&anonz);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_FACTOR_NONZ,&fnonz);
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_MEMORY,&memory);
printf("Assembled non-zeros= %lld\n",anonz);
printf("Factor non-zeros= %lld\n",fnonz);
printf("Total memory= %dMB\n",memory);
/* factorize for direct solution */
/* form preconditioner for iterative solution */
vfs_SysMatrixFactor (sysmatrix);
/* solve */
vfs_SysMatrixSolve (sysmatrix,sysvectorb,sysvectorx);
if(iterative_solver) {
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_NUMITER,&numiter);
printf("Number of iterations= %d\n",numiter);
}
/* access displacement vector using vfs_SysVectorGather */
printf ("Solution at first and last node:\n");
n = 1;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
n = numnp;
vfs_DofTabNodeDof (doftab,n,1,&dof,lm);
vfs_SysVectorGather (sysvectorx,1,lm,&value);
printf(" Node= %4d: %11.4e\n",n,value);
/* delete objects */
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (sysmatrix);
vfs_SysVectorEnd (sysvectorb);
vfs_SysVectorEnd (sysvectorx);
return 0;
}
6.9. Example 3c, Direct and Iterative Solver with Rotated Systems
This example illustrates a 2D plane strain structural model with
rotated coordinate systems used for the x and y displacement components.
The model consists of a simple N by N grid of 2D plain strain
finite elements in the xy plane. Each elemeent is unit square.
Just for purposes of example, a rotated coordinate system is applied
to the in-plane displacements at each node not on a restrained
boundary or at the point of load application.
The rotated system is 90 degrees in the xy plane.
The rotated system at each affected point must be input
to the DofTab object using the function vfs_DofTabNodeDirCos().
#include <stdio.h>
#include <math.h>
#include "vfs/vfs.h"
/*
* 4-node quad solid plane strain stiffness matrix
* E = 1.e+7, nu = 0.3
* unit square element
* Isoparametric technology
*/
static Vdouble ke[36] = {
5.76923e+06,
2.40385e+06, 5.76923e+06,
-3.84615e+06,-4.80769e+05, 5.76923e+06,
4.80769e+05, 9.61538e+05,-2.40385e+06, 5.76923e+06,
-2.88462e+06,-2.40385e+06, 9.61538e+05,-4.80769e+05, 5.76923e+06,
-2.40385e+06,-2.88462e+06, 4.80769e+05,-3.84615e+06, 2.40385e+06, 5.76923e+06,
9.61538e+05, 4.80769e+05,-2.88462e+06, 2.40385e+06,-3.84615e+06,-4.80769e+05, 5.76923e+06,
-4.80769e+05,-3.84615e+06, 2.40385e+06,-2.88462e+06, 4.80769e+05, 9.61538e+05,-2.40385e+06, 5.76923e+06 };
/* element dof map arrays for 4-node quad thermal element */
static Vint nedofs = 8;
static Vint loc[8] = { 1, 1, 2, 2, 3, 3, 4, 4 };
static 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 };
/* topology information */
static Vint shape = SYS_SHAPEQUAD;
static Vint maxi = 2;
static Vint maxj = 0;
static Vint nix = 4;
static Vint xdof = SYS_DOF_TX;
static Vint ydof = SYS_DOF_TY;
static Vint soldof[2] = {SYS_DOF_TX, SYS_DOF_TY};
#define MAX_DOF_ELEMENT 8
#define MAX_NODE_ELEMENT 4
/* flag for iterative solver */
static Vint iterative_solver = 1;
/* mesh size parameter, elements in each direction */
#define N 50
/* load value */
static Vdouble onethousand = 1000.;
/*----------------------------------------------------------------------
Compute Element Connectivity in Simple Grid
/ F = point load
>------------------------------
>| |
>| |
>| |
>| |
>| |
>| |
>| |
>| |
>| |
>| |
>| |
>------------------------------
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
^ = y symmetry constraint
> = x symmetry constraint
----------------------------------------------------------------------*/
static void
get_elementconnectivity(Vint n, Vint ix[])
{
Vint i, j;
j = (n-1)/N+1;
i = n - (j-1)*N;
ix[0] = (N+1)*(j-1) + i;
ix[1] = ix[0]+1;
ix[2] = ix[1]+N+1;
ix[3] = ix[2]-1;
}
static void
transform_stiff(Vint ix[4], Vdouble ket[36])
{
Vdouble R[8][8], tm[2][2], Ke[8][8], T[8][8], x;
Vint i, j, k, n, numnp;
numnp = (N+1)*(N+1);
/* initialize stiffness transformation matrix */
for(i = 0; i < 8; ++i) {
for(j = 0; j < 8; ++j) {
R[i][j] = 0.;
}
R[i][i] = 1.;
}
/* update transformation matrix with nodal rotations */
for(k = 0; k < 4; ++k) {
n = ix[k];
j = (n-1)/(N+1)+1;
i = n - (j-1)*(N+1);
x = (Vdouble)(i-1);
if(i > 1 && j > 1 && n != numnp) {
tm[0][0] = 0.;
tm[0][1] = 1.;
tm[1][0] = -1.;
R[2*k ][2*k ] = R[2*k+1][2*k+1] = tm[0][0];
R[2*k ][2*k+1] = tm[0][1];
R[2*k+1][2*k ] = tm[1][0];
}
}
/* form full stiffness */
for(n = 0, j = 0; j < 8; ++j) {
for(i = 0; i <= j; ++i, ++n) {
Ke[i][j] = ke[n];
if(i != j) {
Ke[j][i] = ke[n];
}
}
}
/* form T = Ke*R^T */
for(i = 0; i < 8; ++i) {
for(j = 0; j < 8; ++j) {
T[i][j] = 0.;
for(k = 0; k < 8; ++k) {
T[i][j] += Ke[i][k]*R[j][k];
}
}
}
/* ket = R*Ke*R^T */
for(n = 0, j = 0; j < 8; ++j) {
for(i = 0; i <= j; ++i, ++n) {
ket[n] = 0.;
for(k = 0; k < 8; ++k) {
ket[n] += R[i][k]*T[k][j];
}
}
}
}
/*----------------------------------------------------------------------
Direct and Iterative Solver with Rotated Systems
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n, numiter, memory;
Vlong anonz, fnonz;
/* declare objects */
vfs_SysMatrix *sysmatrix;
vfs_SysVector *sysvectorb, *sysvectorx;
vfs_DofTab *doftab;
/* number of nodes, elements */
Vint numnp, numel;
/* number of equations and reactions */
Vint neq, nre;
/* number of degrees of freedom in element,
element degree of freedom numbers */
Vint lm[MAX_DOF_ELEMENT];
/* number of nodes in element and connectivity */
Vint ix[MAX_NODE_ELEMENT];
/* node coordinates */
Vdouble x[3];
/* solution value */
Vdouble value[2];
/* nodal transformation */
Vdouble tm[3][3];
/* transformed stiffness */
Vdouble ket[36];
/* compute number of nodes and elements */
numnp = (N+1)*(N+1);
numel = N*N;
doftab = vfs_DofTabBegin();
vfs_DofTabDef (doftab,numnp,numel);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);
/* compute and set node coordinates */
for(n = 1; n <= numnp; n++) {
j = (n-1)/(N+1)+1;
i = n - (j-1)*(N+1);
x[0] = (Vdouble)(i-1);
x[1] = (Vdouble)(j-1);
x[2] = 0.;
vfs_DofTabSetCoords (doftab,n,x);
/* set nodal coordinate transformation on interior nodes */
if(i > 1 && j > 1 && n != numnp) {
tm[0][2] = tm[1][2] = tm[2][0] = tm[2][1] = 0.;
tm[2][2] = 1.;
tm[0][0] = tm[1][1] = 0.;
tm[0][1] = 1.;
tm[1][0] = -1.;
vfs_DofTabNodeDirCos (doftab,n,tm);
}
}
}
/* gather and set element topology,
degree of freedom type and connectivity */
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabElemUse (doftab,n,nix,ix);
/* modification for iterative solver */
if(iterative_solver) {
vfs_DofTabSetTopology (doftab,n,shape,maxi,maxj,0);
}
}
/* set boundary conditions */
for(n = 1; n <= N+1; ++n) {
vfs_DofTabNodePerm (doftab,n,1,&ydof);
}
for(n = 1; n <= numnp; n += N+1) {
vfs_DofTabNodePerm (doftab,n,1,&xdof);
}
/* process DofTab to determine number of equations */
vfs_DofTabProcess (doftab,&neq,&nre);
/* create sysmatrix object */
sysmatrix = vfs_SysMatrixBegin();
if(iterative_solver) {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_ITER,neq);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_UTOL,5.e-7);
vfs_SysMatrixSetParamd (sysmatrix,SYSMATRIX_ITER_FTOL,5.e-3);
} else {
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_SYMM_SPARSE,neq);
}
/* register DofTab with SysMatrix object */
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,doftab);
/* create load and displacement sysvector objects */
sysvectorb = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorb,neq);
sysvectorx = vfs_SysVectorBegin();
vfs_SysVectorDef (sysvectorx,neq);
/* preprocess sysmatrix */
vfs_SysMatrixPreProcess (sysmatrix);
/* Initialize matrix values to zero */
vfs_SysMatrixInit (sysmatrix,0.);
/* assembly for each element */
for(n = 1; n <= numel; n++) {
get_elementconnectivity (n,ix);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
transform_stiff (ix,ket);
vfs_SysMatrixAssem (sysmatrix,nedofs,lm,ket);
}
/* add point load on last node */
vfs_SysVectorInit (sysvectorb,0.);
vfs_DofTabNodeDof (doftab,numnp,1,&xdof,lm);
vfs_SysVectorAssem (sysvectorb,1,lm,&onethousand);
vfs_DofTabNodeDof (doftab,numnp,1,&ydof,lm);
vfs_SysVectorAssem (sysvectorb,1,lm,&onethousand);
/* equation ordering and/or symbolic factorization */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_ASSEM_NONZ,&anonz);
vfs_SysMatrixGetLong (sysmatrix,SYSMATRIX_FACTOR_NONZ,&fnonz);
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_MEMORY,&memory);
printf("Assembled non-zeros= %lld\n",anonz);
printf("Factor non-zeros= %lld\n",fnonz);
printf("Total memory= %dMB\n",memory);
/* factorize for direct solution */
/* form preconditioner for iterative solution */
vfs_SysMatrixFactor (sysmatrix);
/* solve */
vfs_SysMatrixSolve (sysmatrix,sysvectorb,sysvectorx);
if(iterative_solver) {
vfs_SysMatrixGetInteger (sysmatrix,SYSMATRIX_NUMITER,&numiter);
printf("Number of iterations= %d\n",numiter);
}
/* access displacement vector using vfs_SysVectorGather */
printf ("Solution at last node:\n");
n = numnp;
vfs_DofTabNodeDof (doftab,n,2,soldof,lm);
vfs_SysVectorGather (sysvectorx,2,lm,value);
printf(" Node= %4d: (%11.4e, %11.4e)\n",n,value[0],value[1]);
/* delete objects */
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (sysmatrix);
vfs_SysVectorEnd (sysvectorb);
vfs_SysVectorEnd (sysvectorx);
return 0;
}
6.10. Example 4, Direct factorization solution of a non-symmetric system
This example illustrates the assembly, factorization, and solution of a
non-symmetric matrix. The matrix is defined as SYSMATRIX_NSYMM_SPARSE, and
rather than using vfs_SysMatrixAssem() for the
element assembly, the method vfs_SysMatrixAssemFull() is used instead; it
contains the full non-symmetric element matrix to be assembled, defined
row-wise.
#include <stdio.h>
#include "vfs/vfs.h"
/* Full Element Stiffness */
static Vdouble ke[16] = {
4.0, -2.0, -2.0, -2.0,
-1.0, 4.0, -2.0, -2.0,
-1.0, -1.0, 4.0, -2.0,
-1.0, -1.0, -1.0, 4.0 };
/* Global Degrees of Freedom for Four Elements
assumed to be connected as follows:
7 --- 8 --- 9
| 3 | 4 |
4 --- 5 --- 6
| 1 | 2 |
1 --- 2 --- 3
*/
static Vint ix[4][4] = {
{ 1, 2, 4, 5 },
{ 2, 3, 5, 6 },
{ 4, 5, 7, 8 },
{ 5, 6, 8, 9 } };
/* Load vector for Each node */
static Vdouble lval[9] = {
-7.0, -8.0, -1.0, 4.0, 20.0, 16.0, 11.0, 28.0, 17.0 };
/* temperature dof */
static Vint idof = SYS_DOF_TEMP;
/* locs and tags */
static Vint loc[4] = {1,2,3,4};
static Vint tag[4] = {SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP,SYS_DOF_TEMP};
/*----------------------------------------------------------------------
Solve for Unknowns of Four Elements
----------------------------------------------------------------------*/
int
main()
{
Vint i,j;
vfs_SysMatrix *sysmatrix;
vfs_SysVector *load, *soln, *col;
vfs_DofTab *doftab;
Vint numnp, numel;
Vint nedofs, nen;
Vint neq, nre, lm[4];
Vint inod, itag, jnod, jtag;
Vdouble xval,fval, *matrix;
/* Set number of nodes, elements */
numel = 4;
numnp = 9;
/* Set number of freedom and nodes per element */
nedofs = 4;
nen = 4;
/* Create doftab object */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
for(i = 1; i <= numel; ++i) {
vfs_DofTabElemUse (doftab,i,4,ix[i-1]);
}
/* process DofTab to determine number of equations */
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create non-symmetric system stiffness matrix */
sysmatrix = vfs_SysMatrixBegin();
vfs_SysMatrixDef (sysmatrix,SYSMATRIX_NSYMM_SPARSE,neq);
vfs_SysMatrixSetObject (sysmatrix,VFS_DOFTAB,(Vobject*)doftab);
/* preprocess stiffness */
vfs_SysMatrixPreProcess (sysmatrix);
/* Assemble actual stiffness */
vfs_SysMatrixInit (sysmatrix,0.);
for(i = 0; i < numel; ++i) {
vfs_DofTabElemDof (doftab,i+1,nen,ix[i],&nedofs,lm);
vfs_SysMatrixAssemFull (sysmatrix,nedofs,lm,ke);
}
/* Factorize matrix */
vfs_SysMatrixProcess (sysmatrix);
vfs_SysMatrixFactor (sysmatrix);
/* Create system vector for load */
load = vfs_SysVectorBegin ();
vfs_SysVectorDef (load,neq);
/* Set load vector */
for(i = 1; i <= neq; ++i) {
vfs_DofTabNodeDof (doftab,i,1,&idof,lm);
vfs_SysVectorAssem (load,1,lm,&lval[i-1]);
}
/* 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("Solution vector\n");
for(i = 1; i <= numnp; i++) {
vfs_DofTabNodeDof (doftab,i,1,&idof,lm);
vfs_SysVectorGather (soln,1,lm,&xval);
printf("x[%d]= %e\n",i,xval);
}
/* Verification */
vfs_SysMatrixProd (sysmatrix,soln,load);
printf("Original load vector\n");
for(i = 1; i <= numnp; i++) {
vfs_DofTabNodeDof (doftab,i,1,&idof,lm);
vfs_SysVectorGather (load,1,lm,&fval);
printf("f[%d]= %e\n",i,fval);
}
/* Print Matrix */
printf("Assembled Matrix\n");
col = vfs_SysVectorBegin ();
vfs_SysVectorDef (col,neq);
matrix = (Vdouble*)malloc(neq*neq*sizeof(Vdouble));
for(i = 1; i <= neq; i++) {
vfs_SysMatrixColumn (sysmatrix,i,col);
vfs_DofTabInvDof (doftab,1,&i,&inod,&itag);
for(j = 1; j <= neq; j++) {
vfs_DofTabInvDof (doftab,1,&j,&jnod,&jtag);
vfs_SysVectorGather (col,1,&j,&fval);
matrix[neq*(jnod-1)+inod-1] = fval;
}
}
for(i = 0; i < neq; ++i) {
for(j = 0; j < neq; ++j) {
printf("%5.2f ",matrix[neq*i+j]);
}
printf("\n");
}
free (matrix);
/* Delete objects */
vfs_SysVectorEnd (col);
vfs_DofTabEnd (doftab);
vfs_SysVectorEnd (load);
vfs_SysVectorEnd (soln);
vfs_SysMatrixEnd (sysmatrix);
return 0;
}
The output of this example program appears below.
Solution vector
x[1]= -5.066421e+01
x[2]= -4.098059e+01
x[3]= -3.291091e+01
x[4]= -3.170109e+01
x[5]= -2.514674e+01
x[6]= -1.968478e+01
x[7]= -1.853746e+01
x[8]= -1.415101e+01
x[9]= -1.049563e+01
Original load vector
f[1]= -7.000000e+00
f[2]= -8.000000e+00
f[3]= -1.000000e+00
f[4]= 4.000000e+00
f[5]= 2.000000e+01
f[6]= 1.600000e+01
f[7]= 1.100000e+01
f[8]= 2.800000e+01
f[9]= 1.700000e+01
Assembled Matrix
4.00 -2.00 0.00 -2.00 -2.00 0.00 0.00 0.00 0.00
-1.00 8.00 -2.00 -2.00 -4.00 -2.00 0.00 0.00 0.00
0.00 -1.00 4.00 0.00 -2.00 -2.00 0.00 0.00 0.00
-1.00 -1.00 0.00 8.00 -4.00 0.00 -2.00 -2.00 0.00
-1.00 -2.00 -1.00 -2.00 16.00 -4.00 -2.00 -4.00 -2.00
0.00 -1.00 -1.00 0.00 -2.00 8.00 0.00 -2.00 -2.00
0.00 0.00 0.00 -1.00 -1.00 0.00 4.00 -2.00 0.00
0.00 0.00 0.00 -1.00 -2.00 -1.00 -1.00 8.00 -2.00
0.00 0.00 0.00 0.00 -1.00 -1.00 0.00 -1.00 4.00
6.11. Example 5, Linear Static Solution Procedure using SolProc
This example illustrates the use of SolProc to solve a linear static structural problem, a one-dimensional truss structure consisting of 10 elements aligned in the x-direction, fixed at one end, and subjected to gravity load in the x-direction.
SolProc uses a reverse communication methodology: once the problem
type is defined, calls to vfs_SolProcExecute() return the next action to be
taken.
The actions described below are used in this example. VFS_ACTION_INITIAL: the
user is able to set an initial value for the solution vector. The example
initializes the solution vector to zero. VFS_ACTION_PREDICTOR: a
predictor phase is used to specify boundary conditions. Since only
homogeneous boundary conditions are used in this example, this action is a
no-op. VFS_ACTION_STIFF: Form the stiffness matrix. A simple element
stiffness matrix is assembled into the vfs_SysMatrix() object.
VFS_ACTION_FACTOR: The stiffness matrix is factored. VFS_ACTION_REACT:
the reaction forces for this problem are assembled. The SysVector
object to be used is retrieved from SolProc; both internal and
external loads are assembled. VFS_ACTION_SOLVE: the linear system of
equations is solved. Both the right-hand-side and the solution vectors are
retrieved from SolProc. VFS_ACTION_UPDATE: the total solution
vector is updated in the application. In this example, this is also a
no-op. Finally, VFS_ACTION_ENDSTEP, an indication that the solution is
complete and the final solution can be retrieved from SolProc. Since
a single step/increment is requested, the next action is
VFS_ACTION_COMPLETE, indicating that the process ha terminated.
#include <stdio.h>
#include "vfs/vfs.h"
/*
Use vfs_SolProc for a unit-length truss subjected to gravity load in
the x-direction
N=1 N=2 N=3...
(fixed) |> x----------x----------x
E=[1] E=[2]...
Example: Linear static analysis, Young=1., Area=1., Length=1.
*/
#define NUMEL 10
static Vint idof = SYS_DOF_TX;
static Vint loc[2] = {1,2};
static Vint tag[2] = {SYS_DOF_TX,SYS_DOF_TX};
/*----------------------------------------------------------------------
Linear Static Solution Procedure using SolProc
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff;
vfs_SysVector *u,*react;
vfs_DofTab *doftab;
Vint action,n,lm[2], numnp, ix[2], nedofs;
Vint neq, nre;
Vdouble sol,ul[2],fl[2],ke[3];
/* Construct DofTab object */
doftab = vfs_DofTabBegin ();
numnp = NUMEL + 1;
vfs_DofTabDef (doftab,numnp,NUMEL);
vfs_DofTabSetMap (doftab,2,loc,tag);
for(n = 1; n <= NUMEL; ++n) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,1,&idof);
vfs_DofTabProcess (doftab,&neq,&nre);
/* Construct element stiffness matrix */
ke[0] = (Vdouble)NUMEL;
ke[1] = -(Vdouble)NUMEL;
ke[2] = (Vdouble)NUMEL;
/* Instance and prepare SysMatrix object */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,0);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_STATIC,SOLPROC_LINEAR,1.);
/* Infinite loop for reverse communication */
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
printf("Form initial conditions\n");
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SysVectorInit (u,0.);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: no action in this application */
printf("Form predictor\n");
} else if(action == VFS_ACTION_STIFF) {
/* Form stiffness matrix */
printf("Form stiffness matrix\n");
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
printf("Factor stiffness matrix\n");
vfs_SysMatrixFactor (stiff);
} else if(action == VFS_ACTION_REACT) {
/* Form reaction force and loads */
printf("Form total reaction force\n");
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
/* applied gravity force */
fl[0] = fl[1] = 1./NUMEL;
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
/* element reaction force */
ul[0] = ul[1] = 0.;
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorGather (u,nedofs,lm,ul);
fl[0] = -(ke[0]*ul[0] + ke[1]*ul[1]);
fl[1] = -(ke[1]*ul[0] + ke[2]*ul[1]);
vfs_SysVectorAssem (react,2,lm,fl);
}
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
printf("Solve linear system\n");
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* Update solution */
printf("Update solution in application\n");
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
printf("End time step\n");
printf(" Displacements:\n");
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
printf(" %11.4e\n",0.);
for(n = 1; n <= numnp; n++) {
sol = 0.;
vfs_DofTabNodeDof (doftab,n,1,&idof,lm);
vfs_SysVectorGather (u,1,lm,&sol);
printf(" %11.4e\n",sol);
}
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (stiff);
vfs_DofTabEnd (doftab);
return 0;
}
The output of this example program appears below.
Form initial conditions
Form predictor
Form stiffness matrix
Factor stiffness matrix
Form total reaction force
Solve linear system
Update solution in application
End time step
Displacements:
0.0000e+00
1.9000e-01
3.6000e-01
5.1000e-01
6.4000e-01
7.5000e-01
8.4000e-01
9.1000e-01
9.6000e-01
9.9000e-01
1.0000e+00
Normal termination
6.12. Example 5a, Linear Dynamic Solution Procedure
This example illustrates the use of SolProc to solve a linear dynamic
structural problem. It complements example 5 in the sense that additional
actions are required from the user as vfs_SolProcExecute() is called. The
same structure is initially stretched and is then released, allowing it to
vibrate like a linear spring-mass system.
Dynamic structural problems are characterized by the solution of three state variables: displacements, velocities, and accelerations. Initial conditions are given for displacements and velocities only. As a result, the initial acceleration must be computed from the balance of linear momentum. This calculation requires that the mass matrix only be used in the linear equation solver. For this reason, it is customary to employ a diagonalized mass matrix, since its factorization is readily available. This linear solution against the mass matrix is only required once, during initialization. It is accomplished by two required actions: VFS_ACTION_FACTORMASS and VFS_ACTION_SOLVEMASS.
In transient problems, both structural and thermal, a modified stiffness
matrix is required: to the stiffness matrix one must add the mass matrix
which is multiplied by a scalar parameter dependent on time integration
parameters, a mass factor. This parameter can be extracted with
vfs_SolProcGetDouble(), the parameter being SOLPROC_MASSFACTOR.
SolProc also requires that the inertia forces be computed, in VFS_ACTION_INERTIA. These are given by the product of the mass matrix with the acceleration vector.
#include <stdio.h>
#include "vfs/vfs.h"
/*
Use vfs_SolProc for the transient analysis of a unit-length truss
subjected to an initial strain in the x-direction
Example: Linear static analysis, Young=1., Area=1., Length=1.
*/
#define NUMEL 10
static Vint idof = SYS_DOF_TX;
static Vint loc[2] = {1,2};
static Vint tag[2] = {SYS_DOF_TX,SYS_DOF_TX};
/*----------------------------------------------------------------------
Linear Dynamic Solution Procedure
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff,*mass;
vfs_SysVector *u,*v,*a,*react,*inert;
vfs_DofTab *doftab;
Vint action,n,lm[2], numnp, ix[2], istep, nedofs;
Vint neq, nre;
Vdouble sol,vel,accel,ul[2],fl[2],time,massfactor,ke[3],km[2];
/* Construct DofTab object */
doftab = vfs_DofTabBegin ();
numnp = NUMEL + 1;
vfs_DofTabDef (doftab,numnp,NUMEL);
vfs_DofTabSetMap (doftab,2,loc,tag);
for(n = 1; n <= NUMEL; ++n) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,1,&idof);
vfs_DofTabProcess (doftab,&neq,&nre);
/* Construct element stiffness matrix */
ke[0] = (Vdouble)NUMEL;
ke[1] = -(Vdouble)NUMEL;
ke[2] = (Vdouble)NUMEL;
/* Construct element diagonal mass matrix */
km[0] = 0.5/(Vdouble)NUMEL;
km[1] = 0.5/(Vdouble)NUMEL;
/* Instance and prepare stiffness matrix object */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance and prepare mass matrix object */
mass = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (mass,SYSMATRIX_DIAG,neq);
vfs_SysMatrixPreProcess (mass);
vfs_SysMatrixProcess (mass);
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysMatrixAssemDiag (mass,nedofs,lm,km);
}
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,0);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy(solproc,SYS_SOL_TRANSIENT,SOLPROC_LINEAR,16.*3.14159);
vfs_SolProcInitTimeStep (solproc,3.14159/12.);
/* Infinite loop for reverse communication */
istep = 0;
printf(" Time Disp Vel Accel:\n");
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_V0,&v);
vfs_SysVectorInit (u,0.);
for(n = 1; n <= numnp; n++) {
sol = 0.1*(Vdouble)(n-1)/(Vdouble)neq;
vfs_DofTabNodeDof (doftab,n,1,&idof,lm);
vfs_SysVectorScatter (v,1,lm,&sol);
}
} else if(action == VFS_ACTION_FACTORMASS) {
/* Factor mass matrix */
vfs_SysMatrixFactor (mass);
} else if(action == VFS_ACTION_SOLVEMASS) {
/* Solve linear system for mass matrix */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (mass,react,u);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: no action */
} else if(action == VFS_ACTION_STIFF) {
/* Form stiffness matrix */
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
vfs_SolProcGetDouble (solproc,SOLPROC_MASSFACTOR,&massfactor);
vfs_SysMatrixSum (stiff,massfactor,mass);
vfs_SysMatrixFactor (stiff);
} else if(action == VFS_ACTION_REACT) {
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
/* reaction force */
for(n = 1; n <= NUMEL; n++) {
ul[0] = ul[1] = 0.;
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorGather (u,nedofs,lm,ul);
fl[0] = -(ke[0]*ul[0] + ke[1]*ul[1]);
fl[1] = -(ke[1]*ul[0] + ke[2]*ul[1]);
vfs_SysVectorAssem (react,2,lm,fl);
}
} else if(action == VFS_ACTION_INERTIA) {
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_INERTIA,&inert);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_A1,&u);
vfs_SysMatrixProd (mass,u,inert);
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* Update solution */
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
++istep;
vfs_SolProcGetDouble (solproc,SOLPROC_TIME,&time);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_V1,&v);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_A1,&a);
vfs_DofTabNodeDof (doftab,numnp,1,&idof,lm);
vfs_SysVectorGather (u,1,lm,&sol);
vfs_SysVectorGather (v,1,lm,&vel);
vfs_SysVectorGather (a,1,lm,&accel);
printf(" %11.4e %11.4e %11.4e %11.4e\n",time,sol,vel,accel);
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (mass);
vfs_SysMatrixEnd (stiff);
vfs_DofTabEnd (doftab);
return 0;
}
The output of this example program is not reproduced here. It simply displays the displacement, velocity, and acceleration of the end node.
6.13. Example 5bvfe, Linear Dynamic Solution Procedure using VfeTools
This example is almost identical to example 5a, except that it uses Truss3D, MatlFun, and LinMat to generate the element stiffness, mass, reaction, and inertia. It illustrates how to combine DofTab, SolProc, SysMatrix, and SysVector along with element generating modules in VfeTools.
This example also illustrates stresses and strains being retrieved for output at the end of a time step.
#include <stdio.h>
#include "vfe/vfe.h"
#include "vfs/vfs.h"
/*
Use solproc for a NUMEL-elements linear transient structural problem
N=1 N=2 N=3
(fixed) |> x----------x----------x ...
E=[1] E=[2]
*/
#define NUMEL 10
static Vdouble elasprop[2] = {1.,0.};
static Vdouble density = 1.;
static Vdouble length = 1.;
static Vint itx = SYS_DOF_TX;
static Vint ity = SYS_DOF_TY;
static Vint itz = SYS_DOF_TZ;
/*----------------------------------------------------------------------
Linear Dynamic Solution Procedure using VfeTools
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff,*mass;
vfs_SysVector *u,*v,*a,*react,*inert;
vfs_DofTab *doftab;
vfe_Truss3D *truss3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
Vint action,n,lm[6],istep,nedofs, numnp, neq, nre, ix[2];
Vint loc[12], tag[12];
Vdouble sol,vel,accel,ul[12],fl[12],time,mf,ke[78],km[12],x[2][3],h;
Vdouble strs[12],strn[12];
/* Construct element stiffness matrix */
truss3d = vfe_Truss3DBegin ();
vfe_Truss3DSetTopology (truss3d,2);
vfe_Truss3DNumDof (truss3d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Truss3DDofMap (truss3d,SYS_ANALYSIS_STRUCTURAL,loc,tag);
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vfe_LinMatSetElasProp (linmat,elasprop);
vfe_LinMatSetDensity (linmat,density);
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat,matlfun);
vfe_Truss3DSetObject (truss3d,VFE_MATLFUN,matlfun);
/* Mesh size */
h = length/NUMEL;
numnp = NUMEL+1;
/* Construct DofTab object */
doftab = vfs_DofTabBegin ();
numnp = NUMEL + 1;
vfs_DofTabDef (doftab,numnp,NUMEL);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
for(n = 1; n <= NUMEL; ++n) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,1,&itx);
for(n = 1; n <= numnp; ++n) {
vfs_DofTabNodePerm (doftab,n,1,&ity);
vfs_DofTabNodePerm (doftab,n,1,&itz);
}
vfs_DofTabProcess (doftab,&neq,&nre);
/* Instance and prepare stiffness matrix object */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance and form mass matrix object */
mass = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (mass,SYSMATRIX_DIAG,neq);
vfs_SysMatrixPreProcess (mass);
vfs_SysMatrixProcess (mass);
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DMassDiag (truss3d,x,km);
vfs_SysMatrixAssemDiag (mass,nedofs,lm,km);
}
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,0);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_TRANSIENT,SOLPROC_LINEAR,16.*3.1416);
vfs_SolProcInitTimeStep (solproc,3.14159/12.);
/* Infinite loop for reverse communication */
istep = 0;
printf(
"Time Disp Vel Accel Stress Strain\n");
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_V0,&v);
vfs_SysVectorInit (v,0.);
for(n = 1; n <= numnp; n++) {
sol = (Vdouble)(n-1)/(Vdouble)neq;
vfs_DofTabNodeDof (doftab,n,1,&itx,lm);
vfs_SysVectorScatter (u,1,lm,&sol);
}
} else if(action == VFS_ACTION_FACTORMASS) {
/* Factor mass matrix */
vfs_SysMatrixFactor (mass);
} else if(action == VFS_ACTION_SOLVEMASS) {
/* Solve linear system for mass matrix */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (mass,react,u);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: no action */
} else if(action == VFS_ACTION_STIFF) {
/* Form stiffness matrix */
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= NUMEL; n++) {
ix[0] = n;
ix[1] = n+1;
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DStiff (truss3d,x,ke);
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
vfs_SolProcGetDouble (solproc,SOLPROC_MASSFACTOR,&mf);
vfs_SysMatrixSum (stiff,mf,mass);
vfs_SysMatrixFactor (stiff);
} else if(action == VFS_ACTION_REACT) {
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
/* reaction force */
for(n = 1; n <= NUMEL; n++) {
ul[0] = ul[1] = ul[2] = ul[3] = ul[4] = ul[5] = 0.;
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DReact (truss3d,x,ul,fl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
} else if(action == VFS_ACTION_INERTIA) {
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_INERTIA,&inert);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_A1,&u);
vfs_SysMatrixProd (mass,u,inert);
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* Update solution */
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
++istep;
vfs_SolProcGetDouble (solproc,SOLPROC_TIME,&time);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_V1,&v);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_A1,&a);
vfs_DofTabNodeDof (doftab,numnp,1,&itx,lm);
vfs_SysVectorGather (u,1,lm,&sol);
vfs_SysVectorGather (v,1,lm,&vel);
vfs_SysVectorGather (a,1,lm,&accel);
n = NUMEL;
ul[0] = ul[1] = ul[2] = ul[3] = ul[4] = ul[5] = 0.;
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DStrsStrn (truss3d,x,ul,strs,strn);
printf("%11.4e %11.4e %11.4e %11.4e %11.4e %11.4e\n",
time,sol,vel,accel,strs[0],strn[0]);
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (mass);
vfs_SysMatrixEnd (stiff);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
vfe_Truss3DEnd (truss3d);
vfs_DofTabEnd (doftab);
return 0;
}
6.14. Example 5c, Nonlinear Static Solution Procedure
This example augments example 5 by adding a nonlinear material model. This requires that the nonlinear driver in SolProc be used. A nonlinear problem can be driven with or without line search. The linesearch computes a scalar that multiplies the solution increment at each iteration and is mostly useful discontinuous nonlinearities such as gap/contact problems.
A nonlinear problem distinguishes itself from a linear problem in many ways. One of them is the fact that a nonlinear stiffness matrix is needed. This indicates the fact that the tangent matrix depends on the solution. Unless specifically requested, SolProc will try to reduce the number of times a new stiffness matrix is needed in order to reduce the total CPU time - usually associated with the factorization process. SolProc requests a nonlinear stiffness matrix with the action VFS_ACTION_REACTSTIFF, indicating that both the reaction and the nonlinear stiffness matrix are being requested.
#include <stdio.h>
#include <math.h>
#include "vfs/vfs.h"
/*
Use solproc for a NUMEL-element structural problem of a unit-length truss
subjected to gravity load in the x-direction
N=1 N=2 N=3...
(fixed) |> x----------x----------x
E=[1] E=[2]...
Example: Linear static analysis, Young=1.+e_xx, Area=1., Length=1.
(stiffening example, as Young's modulus varies linearly with strain)
*/
/*
model: 0 => Linear
1 => Quadratic
*/
static Vint model = 1;
static Vint idof = SYS_DOF_TX;
static Vint loc[2] = {1,2};
static Vint tag[2] = {SYS_DOF_TX,SYS_DOF_TX};
#define NUMEL 10
static void
vki_Reaction(vfs_SolProc *solproc, vfs_DofTab *doftab);
static void
vki_ElemReact(Vint n, vfs_SysVector *u, Vdouble fl[2], vfs_DofTab *doftab);
static void
vki_ElemStiff(Vint n, vfs_SysVector *u, Vdouble ke[3], vfs_DofTab *doftab);
/*----------------------------------------------------------------------
Nonlinear Static Solution Procedure
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff;
vfs_SysVector *u,*react;
vfs_DofTab *doftab;
Vint action,n,lm[2],factorized, neq, nre, ix[2], numnp, nedofs;
Vdouble sol,ke[3];
numnp = NUMEL+1;
/* Construct DofTab object */
doftab = vfs_DofTabBegin ();
numnp = NUMEL + 1;
vfs_DofTabDef (doftab,numnp,NUMEL);
vfs_DofTabSetMap (doftab,2,loc,tag);
for(n = 1; n <= NUMEL; ++n) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,1,&idof);
vfs_DofTabProcess (doftab,&neq,&nre);
/* Construct element stiffness matrix */
ke[0] = (Vdouble)NUMEL;
ke[1] = -(Vdouble)NUMEL;
ke[2] = (Vdouble)NUMEL;
/* Instance and prepare SysMatrix object */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
vfs_SolProcSetParamd (solproc,SOLPROC_FTOL,1.e-10);
vfs_SolProcSetParami (solproc,SOLPROC_LINESEARCH_FLAG,SYS_ON);
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,0);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_STATIC,SOLPROC_NONLINEAR,1.);
/* Infinite loop for reverse communication */
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SysVectorInit (u,0.);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: no action */
} else if(action == VFS_ACTION_REACTSTIFF) {
/* Form stiffness matrix */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= NUMEL; n++) {
vki_ElemStiff(n,u,ke,doftab);
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
vki_Reaction(solproc,doftab);
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
vfs_SysMatrixFactor (stiff);
vfs_SysMatrixGetInteger (stiff,SYSMATRIX_FACTORIZED,&factorized);
if(!factorized) {
printf("Matrix is singular\n");
return 0;
}
} else if(action == VFS_ACTION_REACT) {
/* Form reaction force and loads */
vki_Reaction(solproc,doftab);
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* Update solution */
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
printf("End time step\n");
printf(" Displacements:\n");
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
printf(" %11.4e %11.4e\n",0.,0.);
for(n = 1; n <= numnp; n++) {
sol = 0.;
vfs_DofTabNodeDof (doftab,n,1,&idof,lm);
vfs_SysVectorGather (u,1,lm,&sol);
printf(" %11.4e %11.4e\n",(Vdouble)n/(Vdouble)NUMEL,sol);
}
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (stiff);
vfs_DofTabEnd (doftab);
return 0;
}
static void
vki_Reaction(vfs_SolProc *solproc, vfs_DofTab *doftab)
{
vfs_SysVector *react,*u;
Vint n,lm[2],ix[2], nedofs;
Vdouble fl[2],norm;
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
/* applied gravity force */
for(n = 1; n <= NUMEL; n++) {
fl[0] = fl[1] = 1./NUMEL;
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
/* reaction force */
for(n = 1; n <= NUMEL; n++) {
vki_ElemReact(n,u,fl,doftab);
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
vfs_SysVectorNorm (react,SYSVECTOR_NORM_L2,&norm);
printf(" Reaction norm= %11.4e\n",norm);
}
static void
vki_ElemReact(Vint n, vfs_SysVector *u, Vdouble fl[2], vfs_DofTab *doftab)
{
Vint lm[2], ix[2], nedofs;
Vdouble ul[2],exx,E,sxx,h;
/* gather local displacements */
ix[0] = n; ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
ul[0] = 0.;ul[1] = 0.;
vfs_SysVectorGather (u,nedofs,lm,ul);
/* element size */
h = 1./NUMEL;
/* compute strain */
exx = (ul[1]-ul[0])/h;
/* form Young's modulus */
if(model == 0) {
E = 1.;
} else if(model == 1) {
E = 1.+exx;
} else if(model == 2) {
E = exp(exx);
}
/* form stress */
sxx = E*exx;
/* form element reaction force */
sxx = E*exx;
fl[0] = sxx;
fl[1] = -sxx;
}
static void
vki_ElemStiff(Vint n, vfs_SysVector *u, Vdouble ke[3], vfs_DofTab *doftab)
{
Vint lm[2], ix[2], nedofs;
Vdouble ul[2],exx,E,h;
/* gather local displacements */
ix[0] = n; ix[1] = n+1;
ul[0] = 0.;ul[1] = 0.;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
vfs_SysVectorGather (u,nedofs,lm,ul);
/* element size */
h = 1./NUMEL;
/* compute strain */
exx = (ul[1]-ul[0])/h;
/* form derivative of stress w.r.t. strain (NOT Young's modulus) */
if(model == 0) {
E = 1.;
} else if(model == 1) {
E = 1.+2.*exx;
}
/* form element stiffness matrix */
ke[0] = E/h;
ke[1] = -E/h;
ke[2] = E/h;
}
The output of this example program appears below.
Reaction norm= 6.0828e-01
Reaction norm= 1.3506e+00
Reaction norm= 1.0848e+00
Reaction norm= 2.0317e+00
Reaction norm= 2.0317e+00
Reaction norm= 2.5797e+00
Reaction norm= 1.1780e+00
Reaction norm= 1.1780e+00
Reaction norm= 1.1115e+00
Reaction norm= 1.0730e+00
Reaction norm= 1.0730e+00
Reaction norm= 2.7614e-01
Reaction norm= 1.2271e-01
Reaction norm= 1.2271e-01
Reaction norm= 1.8569e-03
Reaction norm= 5.5774e-05
Reaction norm= 1.6628e-06
Reaction norm= 4.9587e-08
Reaction norm= 1.4787e-09
Reaction norm= 4.4094e-11
End time step
Displacements:
0.0000e+00 0.0000e+00
1.0000e-01 9.6629e-02
2.0000e-01 1.8627e-01
3.0000e-01 2.6856e-01
4.0000e-01 3.4306e-01
5.0000e-01 4.0925e-01
6.0000e-01 4.6649e-01
7.0000e-01 5.1395e-01
8.0000e-01 5.5056e-01
9.0000e-01 5.7472e-01
1.0000e+00 5.8388e-01
Normal termination
6.15. Example 5dvfe, Pre-Stressed Nonlinear Static Solution Procedure
This example illustrates the introduction of a pre-stressed state in the model. These are particularly useful in truss structures which are unable to carry moment loads across. Hence, unless a pre-stress is used to generate an added geometric stiffness, the resulting stiffness matrix could be singular.
This example takes advantage of the fact that for all element technologies that support large strains in VfeTools the geometric stiffness is automatically added to the stiffness matrix in a large strain analysis.
An artificially large gravity load is used to demonstrate SolProc’s ability to recover from a lack of convergence by reducing the time step in the problem. Since this is a static analysis, the time parameter is simply used to scale the total load being applied.
The residual norm is printed when available - as indicated by the integer parameter SOLPROC_NEWREACT, and shows how the residual is driven to zero during each increment.
#include <stdio.h>
#include <math.h>
#include "vfe/vfe.h"
#include "vfs/vfs.h"
/*
Use solproc for a NUMEL-elements nonlinear structural problem
N=1 N=2 N=3
(fixed) |> x----------x----------x ...
E=[1] E=[2]
*/
static Vdouble elasprop[2] = {2.e+11,0.};
static Vdouble density = 6000.;
static Vdouble length = 1.;
static Vint bctags[3] = {SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ};
static Vdouble prestress[12] = {2.e+6,0.,0.,0.,0.,0.,2.e+6,0.,0.,0.,0.,0.};
static Vdouble area[2] = {1.e-4,1.e-4};
static Vdouble gravity[2][3] = { {0.,-1.e+4,0.},{0.,-1.e+4,0.} };
#define NUMEL 10
/*----------------------------------------------------------------------
Pre-Stressed Nonlinear Static Solution Procedure
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff;
vfs_SysVector *u,*react;
vfs_DofTab *doftab;
vfe_Truss3D *truss3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
Vint action,n,lm[6],istep,nedofs;
Vdouble sol[3],ul[6],fl[6],time,ke[21],x[2][3],h,norm;
Vint numnp,numel,locs[6],tags[6],ix[2];
Vint neq,nre,singular,newreact;
Vdouble g[2][3];
/* Construct VfeTools objects */
truss3d = vfe_Truss3DBegin ();
vfe_Truss3DSetTopology (truss3d,2);
vfe_Truss3DSetPropPtr (truss3d,VFE_PROP_AREA,area);
vfe_Truss3DSetPropPtr (truss3d,VFE_PROP_PRESTRESS,prestress);
vfe_Truss3DNumDof (truss3d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Truss3DDofMap (truss3d,SYS_ANALYSIS_STRUCTURAL,locs,tags);
vfe_Truss3DSetParami (truss3d,VFE_STRAINTYPE,VFE_LARGESTRAIN);
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vfe_LinMatSetElasProp (linmat,elasprop);
vfe_LinMatSetDensity (linmat,density);
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat,matlfun);
vfe_Truss3DSetObject (truss3d,VFE_MATLFUN,matlfun);
/* Mesh size */
h = length/NUMEL;
/* current gravity value */
g[0][0] = gravity[0][0]; g[0][1] = gravity[0][1]; g[0][2] = gravity[0][2];
g[1][0] = gravity[1][0]; g[1][1] = gravity[1][1]; g[1][2] = gravity[0][2];
/* Instance and prepare DofTab object */
numnp = NUMEL+1;
numel = NUMEL;
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
vfs_DofTabSetMap (doftab,nedofs,locs,tags);
for(n = 1; n <= numel; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,3,bctags);
vfs_DofTabNodePerm (doftab,numnp,3,bctags);
vfs_DofTabProcess (doftab,&neq,&nre);
/* Instance and prepare sysmatrix object */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,nre);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_STATIC,SOLPROC_NONLINEAR,1.);
vfs_SolProcInitTimeStep (solproc,1.);
/* Infinite loop for reverse communication */
istep = 0;
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
vfs_SolProcGetInteger (solproc,SOLPROC_NEWREACT,&newreact);
if(newreact) {
vfs_SolProcGetDouble (solproc,SOLPROC_REACTNORM,&norm);
printf(" Reaction norm= %11.4e\n",norm);
}
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SysVectorInit (u,0.);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: extract current time */
vfs_SolProcGetDouble (solproc,SOLPROC_TIME,&time);
printf("time= %11.4e\n",time);
g[0][1] = time*gravity[0][1];
g[1][1] = time*gravity[1][1];
} else if(action == VFS_ACTION_REACTSTIFF) {
/* Form stiffness matrix */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysMatrixInit (stiff,0.);
vfs_SysVectorInit (react,0.);
for(n = 1; n <= numel; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
ul[0] = ul[1] = ul[2] = ul[3] = ul[4] = ul[5] = 0.;
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DReactStiff (truss3d,x,ul,SYS_ON,fl,ke);
vfs_SysVectorAssem (react,nedofs,lm,fl);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
/* add gravity load */
vfe_Truss3DElemLoad (truss3d,x,g,fl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
vfs_SysMatrixFactor (stiff);
vfs_SysMatrixGetInteger (stiff,SYSMATRIX_SINGULAR,&singular);
if(singular) {
printf("Singular matrix detected... abort\n");
return 0;
}
} else if(action == VFS_ACTION_REACT) {
/* Form reaction force and loads */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
/* reaction force */
for(n = 1; n <= numel; n++) {
ix[0] = n;
ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
ul[0] = ul[1] = ul[2] = ul[3] = ul[4] = ul[5] = 0.;
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*h; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*h; x[1][1] = 0.; x[1][2] = 0.;
vfe_Truss3DReact (truss3d,x,ul,fl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
/* add gravity load */
vfe_Truss3DElemLoad (truss3d,x,g,fl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* Update solution */
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
++istep;
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,3,bctags,lm);
sol[0] = sol[1] = sol[2] = 0.;
vfs_SysVectorGather (u,3,lm,sol);
printf(" %3d %11.4e\n",n,sol[1]);
}
} else {
printf("missing action= %d\n",action);
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (stiff);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
vfe_Truss3DEnd (truss3d);
vfs_DofTabEnd (doftab);
return 0;
}
The output of this example program appears below.
time= 1.0000e+00
Reaction norm= 1.8000e+03
Reaction norm= 2.1750e+10
Reaction norm= 2.3731e+31
time= 1.0000e+00
Reaction norm= 1.8000e+03
Reaction norm= 2.1750e+10
Reaction norm= 6.4380e+09
Reaction norm= 1.9020e+09
Reaction norm= 5.5926e+08
Reaction norm= 1.6298e+08
Reaction norm= 4.6359e+07
Reaction norm= 1.2776e+07
Reaction norm= 3.6805e+06
Reaction norm= 1.3161e+06
Reaction norm= 5.6559e+05
Reaction norm= 2.5451e+05
Reaction norm= 1.2795e+05
Reaction norm= 6.4395e+04
Reaction norm= 3.4672e+04
Reaction norm= 1.7577e+04
Reaction norm= 5.9829e+03
time= 5.0000e-01
Reaction norm= 9.0000e+02
Reaction norm= 2.7355e+09
Reaction norm= 8.0941e+08
Reaction norm= 2.3598e+08
Reaction norm= 6.7668e+07
Reaction norm= 1.8670e+07
Reaction norm= 5.0604e+06
Reaction norm= 1.6888e+06
Reaction norm= 6.8326e+05
Reaction norm= 3.0545e+05
Reaction norm= 1.4965e+05
Reaction norm= 7.5108e+04
Reaction norm= 3.9982e+04
Reaction norm= 2.1372e+04
Reaction norm= 1.0731e+04
Reaction norm= 3.4205e+03
Reaction norm= 4.1956e+02
time= 2.5000e-01
Reaction norm= 4.5000e+02
Reaction norm= 3.5016e+08
Reaction norm= 1.1995e+08
Reaction norm= 3.1952e+07
Reaction norm= 8.4185e+06
Reaction norm= 2.4041e+06
Reaction norm= 8.5680e+05
Reaction norm= 3.6899e+05
Reaction norm= 1.6932e+05
Reaction norm= 8.5767e+04
Reaction norm= 4.4041e+04
Reaction norm= 2.3816e+04
Reaction norm= 1.2824e+04
Reaction norm= 6.2380e+03
Reaction norm= 1.8576e+03
Reaction norm= 1.9975e+02
Reaction norm= 2.2284e+00
time= 1.2500e-01
Reaction norm= 2.2500e+02
Reaction norm= 4.7663e+07
Reaction norm= 4.3337e+07
Reaction norm= 1.0834e+07
Reaction norm= 2.6846e+06
Reaction norm= 8.2433e+05
Reaction norm= 3.1147e+05
Reaction norm= 1.2832e+05
Reaction norm= 5.7081e+04
Reaction norm= 2.7632e+04
Reaction norm= 1.3709e+04
Reaction norm= 7.1491e+03
Reaction norm= 3.3596e+03
Reaction norm= 9.1149e+02
Reaction norm= 5.5672e+01
Reaction norm= 2.1125e-01
1 0.0000e+00
2 -4.3200e-03
3 -7.6812e-03
4 -1.0083e-02
5 -1.1524e-02
6 -1.2004e-02
7 -1.1524e-02
8 -1.0083e-02
9 -7.6812e-03
10 -4.3200e-03
11 0.0000e+00
time= 2.5000e-01
Reaction norm= 2.2500e+02
Reaction norm= 1.6410e+03
Reaction norm= 6.1724e+01
Reaction norm= 1.4552e+00
Reaction norm= 1.1374e-04
1 0.0000e+00
2 -5.4590e-03
3 -9.7074e-03
4 -1.2743e-02
5 -1.4565e-02
6 -1.5173e-02
7 -1.4565e-02
8 -1.2743e-02
9 -9.7074e-03
10 -5.4590e-03
11 0.0000e+00
time= 3.7500e-01
Reaction norm= 2.2500e+02
Reaction norm= 6.4620e+02
Reaction norm= 1.0439e+01
Reaction norm= 2.5229e-02
1 0.0000e+00
2 -6.2558e-03
3 -1.1125e-02
4 -1.4605e-02
5 -1.6694e-02
6 -1.7391e-02
7 -1.6694e-02
8 -1.4605e-02
9 -1.1125e-02
10 -6.2558e-03
11 0.0000e+00
time= 5.0000e-01
Reaction norm= 2.2500e+02
Reaction norm= 3.7493e+02
Reaction norm= 3.2531e+00
Reaction norm= 1.8706e-03
1 0.0000e+00
2 -6.8892e-03
3 -1.2252e-02
4 -1.6086e-02
5 -1.8387e-02
6 -1.9155e-02
7 -1.8387e-02
8 -1.6086e-02
9 -1.2252e-02
10 -6.8892e-03
11 0.0000e+00
time= 6.2500e-01
Reaction norm= 2.2500e+02
Reaction norm= 2.5487e+02
Reaction norm= 1.3695e+00
Reaction norm= 2.7486e-04
1 0.0000e+00
2 -7.4234e-03
3 -1.3203e-02
4 -1.7335e-02
5 -1.9816e-02
6 -2.0643e-02
7 -1.9816e-02
8 -1.7335e-02
9 -1.3203e-02
10 -7.4234e-03
11 0.0000e+00
time= 7.5000e-01
Reaction norm= 2.2500e+02
Reaction norm= 1.8893e+02
Reaction norm= 6.8884e-01
Reaction norm= 6.0140e-05
1 0.0000e+00
2 -7.8900e-03
3 -1.4034e-02
4 -1.8427e-02
5 -2.1064e-02
6 -2.1944e-02
7 -2.1064e-02
8 -1.8427e-02
9 -1.4034e-02
10 -7.8900e-03
11 0.0000e+00
time= 8.7500e-01
Reaction norm= 2.2500e+02
Reaction norm= 1.4793e+02
Reaction norm= 1.7690e+01
Reaction norm= 3.4706e+00
Reaction norm= 5.8881e-01
Reaction norm= 1.0887e-01
1 0.0000e+00
2 -8.3069e-03
3 -1.4777e-02
4 -1.9403e-02
5 -2.2180e-02
6 -2.3106e-02
7 -2.2180e-02
8 -1.9403e-02
9 -1.4777e-02
10 -8.3069e-03
11 0.0000e+00
time= 1.0000e+00
Reaction norm= 2.2499e+02
Reaction norm= 4.2877e+02
Reaction norm= 1.0694e+02
Reaction norm= 4.1244e+01
Reaction norm= 1.4282e+01
Reaction norm= 5.2934e+00
Reaction norm= 1.9356e+00
Reaction norm= 7.1269e-01
Reaction norm= 2.6163e-01
Reaction norm= 9.5951e-02
1 0.0000e+00
2 -8.6856e-03
3 -1.5451e-02
4 -2.0289e-02
5 -2.3194e-02
6 -2.4162e-02
7 -2.3194e-02
8 -2.0289e-02
9 -1.5451e-02
10 -8.6856e-03
11 0.0000e+00
Normal termination
6.16. Example 5evfevis, Large Rotation Cantilever Beam with Corot
This example illustrates the use of the Corot object to solve a cantilever beam problem subjected to an end-moment load. A critical step in this example is the need to modify the solution update done inside SolProc. The co-rotational formulation employs an exponential map - a map from incremental to large rotations. Since SolProc is unable to determine which degrees of freedom form a nodal rotation, it is up to the user to perform this map during the VFS_ACTION_UPDATE step of the analysis.
The solution to this problem is to turn the straight beam into a perfect circle. In order to help the user to visualize this result, the mesh and the nodal output at each step is placed in an Ideas Universal file.
#include <stdio.h>
#include <math.h>
#include "vis/vis.h"
#include "vfe/vfe.h"
#include "vfs/vfs.h"
/*---------------------------------------------------------------------------
Use solproc for a 10-element cantilever beam bending under
large rotations
N=1 N=2 N=3
(fixed) |> x----------x----------x ...
E=[1] E=[2]
---------------------------------------------------------------------------*/
static Vdouble elasprop[2] = {1.e+7,0.};
static Vint bctags[6] = { SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ,
SYS_DOF_RX,SYS_DOF_RY,SYS_DOF_RZ };
static Vint rots[3] = {SYS_DOF_RX,SYS_DOF_RY,SYS_DOF_RZ};
static Vint numnp = 101;
static Vint numel = 100;
/*----------------------------------------------------------------------
Large Rotation Cantilever Beam with Corot
----------------------------------------------------------------------*/
int
main()
{
vfs_SolProc *solproc;
vfs_SysMatrix *stiff;
vfs_SysVector *u,*react,*un1,*un,*inc;
vfs_DofTab *doftab;
vfe_Beam3D *beam3d;
vfe_MatlFun *matlfun,*matlfun2;
vfe_BeamProp *beamprop;
vfe_LinMat *linmat;
vfe_Corot *corot;
vis_RProp *rprop;
vis_Connect *connect;
vis_State *state;
Vint action,n,i,lm[18],istep,nedofs;
Vdouble sol[6],ul[18],fl[18],ke[171],x[2][3],bprop[10],time,moment;
Vdouble tm[3][3],unv[3],un1v[3],uinc[3],norm;
Vint locs[18],tags[18],ix[2];
Vint neq,nre,singular,newreact;
/* Construct VisTools objects for output purposes */
rprop = vis_RPropBegin ();
vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
vis_RPropSetType(rprop,SYS_RES_D);
connect = vis_ConnectBegin ();
vis_ConnectPre (connect,SYS_DOUBLE);
for(n = 1; n <= numel; n++) {
ix[0] = n; ix[1] = n+1;
vis_ConnectSetTopology (connect,n,SYS_SHAPELINE,2,0,0);
vis_ConnectSetElemNode (connect,n,ix);
}
for(n = 1; n <= numnp; n++) {
x[0][0] = (n-1)*0.1;
x[0][1] = 0.;
x[0][2] = 0.;
vis_ConnectSetCoordsdv (connect,n,x[0]);
}
vis_ConnectWrite (connect,SYS_SDRC_UNIVERSAL,"exam5evfevis.unv");
state = vis_StateBegin ();
vis_StatePre (state,SYS_DOUBLE);
vis_StateDef (state,numnp,SYS_NODE,SYS_NONE,SYS_VECTOR);
/* Construct VfeTools objects */
beam3d = vfe_Beam3DBegin ();
vfe_Beam3DSetTopology (beam3d,2);
vfe_Beam3DSetParami (beam3d,VFE_TECH,VFE_TECH_ENHANCED);
vfe_Beam3DNumDof (beam3d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Beam3DDofMap (beam3d,SYS_ANALYSIS_STRUCTURAL,locs,tags);
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vfe_LinMatSetElasProp (linmat,elasprop);
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat,matlfun);
beamprop = vfe_BeamPropBegin ();
vfe_BeamPropDef (beamprop,BEAMPROP_SECTION);
bprop[0] = 1.1310e-02;
bprop[1] = 1.0179e-05;
bprop[2] = 1.0179e-05;
bprop[3] = 0.;
bprop[4] = 2.0358e-05;
bprop[5] = 5./6.;
bprop[6] = 5./6.;
bprop[7] = 0.;
bprop[8] = 0.;
bprop[9] = 0.;
vfe_BeamPropSetSection (beamprop,bprop,1);
vfe_BeamPropSetMatlFun (beamprop,1,matlfun);
matlfun2 = vfe_MatlFunBegin ();
vfe_BeamPropMatlFun (beamprop,matlfun2);
vfe_Beam3DSetObject (beam3d,VFE_MATLFUN,matlfun2);
corot = vfe_CorotBegin();
vfe_CorotSetTopology (corot,SYS_SHAPELINE,2,0,0);
vfe_CorotSetMap (corot,nedofs,locs,tags);
/* Instance and prepare DofTab object */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
vfs_DofTabSetMap (doftab,nedofs,locs,tags);
for(n = 1; n <= numel; n++) {
ix[0] = n; ix[1] = n+1;
vfs_DofTabElemUse (doftab,n,2,ix);
}
vfs_DofTabNodePerm (doftab,1,6,bctags);
vfs_DofTabProcess (doftab,&neq,&nre);
/* Instance and prepare sysmatrix objects */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,nre);
/* Specify a static, linear analysis for 1 unit of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_STATIC,SOLPROC_NONLINEAR,1.);
vfs_SolProcInitTimeStep (solproc,0.1);
/* Infinite loop for reverse communication */
istep = 0;
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
vfs_SolProcGetInteger (solproc,SOLPROC_NEWREACT,&newreact);
if(newreact) {
vfs_SolProcGetDouble (solproc,SOLPROC_REACTNORM,&norm);
printf(" Reaction norm= %11.4e\n",norm);
}
if(action == VFS_ACTION_INITIAL) {
/* Specify initial conditions */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&u);
vfs_SysVectorInit (u,0.);
} else if(action == VFS_ACTION_PREDICTOR) {
/* Predictor: extract current time and load */
vfs_SolProcGetDouble (solproc,SOLPROC_TIME,&time);
printf("time= %11.4e\n",time);
moment = 0.66e+2*time;
} else if(action == VFS_ACTION_REACTSTIFF) {
/* Form element stiffness and reaction */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= numel; n++) {
ix[0] = n; ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
for(i = 0; i < nedofs; i++) {
ul[i] = 0.;
}
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*0.1; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*0.1; x[1][1] = 0.; x[1][2] = 0.;
vfe_CorotPull (corot,x,ul,tm);
vfe_Beam3DReactStiff (beam3d,x,ul,SYS_ON,fl,ke);
vfe_CorotPush (corot,SYS_ON,fl,SYS_ON,ke);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
/* Add point cantilever moment */
vfs_DofTabNodeDof (doftab,numnp,1,&bctags[5],lm);
vfs_SysVectorAssem (react,1,lm,&moment);
} else if(action == VFS_ACTION_FACTOR) {
/* Factor stiffness matrix */
vfs_SysMatrixFactor (stiff);
vfs_SysMatrixGetInteger (stiff,SYSMATRIX_SINGULAR,&singular);
if(singular) {
printf("Singular matrix detected... abort\n");
return 0;
}
} else if(action == VFS_ACTION_REACT) {
/* Form element stiffness and reaction */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
vfs_SysVectorInit (react,0.);
for(n = 1; n <= numel; n++) {
ix[0] = n; ix[1] = n+1;
vfs_DofTabElemDof (doftab,n,2,ix,&nedofs,lm);
for(i = 0; i < nedofs; i++) {
ul[i] = 0.;
}
vfs_SysVectorGather (u,nedofs,lm,ul);
x[0][0] = (n-1)*0.1; x[0][1] = 0.; x[0][2] = 0.;
x[1][0] = n*0.1; x[1][1] = 0.; x[1][2] = 0.;
vfe_CorotPull (corot,x,ul,tm);
vfe_Beam3DReact (beam3d,x,ul,fl);
vfe_CorotPush (corot,SYS_ON,fl,SYS_OFF,NULL);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
/* Add point cantilever moment */
vfs_DofTabNodeDof (doftab,11,1,&bctags[5],lm);
vfs_SysVectorAssem (react,1,lm,&moment);
} else if(action == VFS_ACTION_SOLVE) {
/* Solve linear system */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&u);
vfs_SysMatrixSolve (stiff,react,u);
} else if(action == VFS_ACTION_UPDATE) {
/* modify rotations */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&un1);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&un);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_INC,&inc);
for(n=1; n<=numnp; ++n) {
unv[0] = unv[1] = unv[2] = 0.;
un1v[0] = un1v[1] = un1v[2] = 0.;
uinc[0] = uinc[1] = uinc[2] = 0.;
vfs_DofTabNodeDof (doftab,n,3,rots,lm);
vfs_SysVectorGather (un,3,lm,unv);
vfs_SysVectorGather (un1,3,lm,un1v);
vfs_SysVectorGather (inc,3,lm,uinc);
vfe_CorotUpdateRot (corot,unv,uinc,un1v);
vfs_SysVectorScatter (un1,3,lm,un1v);
}
/* Update solution */
} else if(action == VFS_ACTION_ENDSTEP) {
/* Terminate step */
++istep;
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&u);
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,3,bctags,lm);
sol[0] = sol[1] = sol[2] = 0.;
vfs_SysVectorGather (u,3,lm,sol);
vis_StateSetDatadv (state,n,sol);
if(n == numnp) {
printf(" %3d %11.4e %11.4e\n",n,sol[0],sol[1]);
}
}
vis_RPropSetIds(rprop,istep,0,0);
vis_StateWrite (state,rprop,SYS_SDRC_UNIVERSAL,"exam5evfevis.unv");
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (stiff);
vfe_MatlFunEnd (matlfun);
vfe_MatlFunEnd (matlfun2);
vfe_LinMatEnd (linmat);
vfe_Beam3DEnd (beam3d);
vfe_BeamPropEnd (beamprop);
vfe_CorotEnd (corot);
vfs_DofTabEnd (doftab);
vis_ConnectEnd (connect);
vis_StateEnd (state);
vis_RPropEnd (rprop);
return 0;
}
The output appears below:
time= 1.0000e-01
Reaction norm= 6.6000e+00
Reaction norm= 2.1871e+04
Reaction norm= 3.0139e+08
time= 1.0000e-01
Reaction norm= 6.6000e+00
Reaction norm= 2.1871e+04
Reaction norm= 7.9631e+02
Reaction norm= 1.0914e+01
Reaction norm= 2.6630e-04
101 -6.8609e-01 3.1300e+00
time= 2.0000e-01
Reaction norm= 1.4758e+01
Reaction norm= 1.4371e+04
Reaction norm= 6.9928e+04
Reaction norm= 9.0774e+03
Reaction norm= 9.9998e+02
Reaction norm= 3.8369e+00
Reaction norm= 1.9385e-04
101 -2.5763e+00 5.6248e+00
time= 3.0000e-01
Reaction norm= 2.3797e+01
Reaction norm= 5.5966e+04
Reaction norm= 1.6090e+05
Reaction norm= 3.5402e+04
Reaction norm= 1.1186e+04
Reaction norm= 4.8195e+02
Reaction norm= 2.7453e+00
Reaction norm= 3.3419e-05
101 -5.2151e+00 7.0211e+00
time= 4.0000e-01
Reaction norm= 3.3000e+01
Reaction norm= 1.1040e+05
Reaction norm= 2.7732e+05
Reaction norm= 8.3204e+04
Reaction norm= 4.2300e+04
Reaction norm= 5.4689e+03
Reaction norm= 3.1182e+02
Reaction norm= 4.7957e-01
Reaction norm= 2.4769e-06
101 -7.9911e+00 7.1469e+00
time= 5.0000e-01
Reaction norm= 4.2261e+01
Reaction norm= 1.7066e+05
Reaction norm= 4.0709e+05
Reaction norm= 1.4839e+05
Reaction norm= 9.4735e+04
Reaction norm= 2.0500e+04
Reaction norm= 3.6310e+03
Reaction norm= 6.6263e+01
Reaction norm= 4.5334e-02
Reaction norm= 1.6593e-08
101 -1.0309e+01 6.1538e+00
time= 6.0000e-01
Reaction norm= 5.1548e+01
Reaction norm= 2.3421e+05
Reaction norm= 5.4460e+05
Reaction norm= 2.2439e+05
Reaction norm= 1.6133e+05
Reaction norm= 4.6502e+04
Reaction norm= 1.4738e+04
Reaction norm= 9.8285e+02
Reaction norm= 9.6385e+00
Reaction norm= 4.9677e-04
101 -1.1750e+01 4.4537e+00
time= 7.0000e-01
Reaction norm= 6.0849e+01
Reaction norm= 3.0011e+05
Reaction norm= 6.8748e+05
Reaction norm= 3.0702e+05
Reaction norm= 2.3619e+05
Reaction norm= 8.1740e+04
Reaction norm= 3.5741e+04
Reaction norm= 4.8045e+03
Reaction norm= 2.1502e+02
Reaction norm= 2.5895e-01
Reaction norm= 6.3083e-07
101 -1.2170e+01 2.5841e+00
time= 8.0000e-01
Reaction norm= 7.0159e+01
Reaction norm= 3.6804e+05
Reaction norm= 8.3477e+05
Reaction norm= 3.9412e+05
Reaction norm= 3.1617e+05
Reaction norm= 1.2366e+05
Reaction norm= 6.5677e+04
Reaction norm= 1.3163e+04
Reaction norm= 1.4542e+03
Reaction norm= 1.2007e+01
Reaction norm= 1.3462e-03
101 -1.1715e+01 1.0467e+00
time= 9.0000e-01
Reaction norm= 7.9475e+01
Reaction norm= 4.3929e+05
Reaction norm= 5.4227e+06
Reaction norm= 2.3663e+06
Reaction norm= 2.5740e+06
Reaction norm= 1.1134e+06
Reaction norm= 1.1510e+06
Reaction norm= 4.7216e+05
Reaction norm= 4.1722e+05
Reaction norm= 1.3649e+05
Reaction norm= 7.3173e+04
Reaction norm= 1.1383e+04
Reaction norm= 9.5925e+02
Reaction norm= 3.9270e+00
Reaction norm= 1.4163e-04
101 -1.0742e+01 1.6887e-01
time= 1.0000e+00
Reaction norm= 8.8794e+01
Reaction norm= 5.0984e+05
Reaction norm= 5.9201e+06
Reaction norm= 2.5770e+06
Reaction norm= 2.7857e+06
Reaction norm= 1.2135e+06
Reaction norm= 1.2417e+06
Reaction norm= 5.1882e+05
Reaction norm= 4.6129e+05
Reaction norm= 1.6788e+05
Reaction norm= 1.0742e+05
Reaction norm= 2.4856e+04
Reaction norm= 5.0424e+03
Reaction norm= 1.2549e+02
Reaction norm= 1.5393e-01
Reaction norm= 1.0804e-07
101 -9.6924e+00 3.0979e-02
Normal termination
6.17. Example 5fvfevis, Nonlinear Elasto-plastic Transient Solution Procedure
This example illustrates the use of the SolProc object to solve a cantilever beam problem modeled with three-dimensional elements subjected to an initial velocity that varies linearly from the fixed end. An elasto-plastic material model is employed. The nonlinear behavior is obvious from the non-harmonic solution obtained as the stress moves beyond the material’s yield point.
A MapMesh object is first employed to mesh a cantilever beam with 8-noded hexahedral elements. Nodes at the fixed face of the beam are marked with a VIS_MISCID1 association value of 5. In addition, a corner node at the free end is marked with a VIS_MISCID2 association of 2. Because of the bending characteristics of this problem the element technology is set to VFE_TECH_MIXED.
An elasto-plastic material model is chosen with a piecewise-linear-varying hardening curve. Properties approximating steel are chosen. The plasticity model requires history variables. The number of history variables per element is obtained by multiplying the number of history variables at each integration point, nhist, by the number of integration points in the element, nintpnt. Two sets of history variables are needed: values at the beginning of a time step, N, and values at the end of the time step, N+1. Rather than copying the values at N+1 to N at the end of a time step, these pointers are simply swapped with each other. During element loops - e.g., forming the stiffness matrix - the pointer to the element’s history variable must be set for each element during the loop.
DofTab is used to number the equations. DofTabNodePerm is called for each node on the fixed end, thereby eliminating all the displacement components.
A diagonal mass matrix is employed for the transient analysis. This is common practice in structural analysis.
The elasto-plastic analysis demands that SolProc use a nonlinear strategy. If no information is available, a short-run analysis should be performed to determine an initial and a minimum time step. In this problem, a small initial time step was chosen with the understanding that it would increase over time. A minimum time step 10 times smaller was set to prevent the problem from taking too long to terminate. No maximum time step was set, as in transient analysis the half-step tolerance takes care of controlling it.
#include <stdio.h>
#include <math.h>
#include "vis/vismesh.h"
#include "vdm/vdm.h"
#include "vfe/vfe.h"
#include "vfs/vfs.h"
/*
Use solproc for a transient nonlinear structural problem
*/
static Vdouble hard[44][2] = {
{ 0.00000e+00, 0.000000e+00 },
{ 0.21052e-03, 0.009996e+09 },
{ 0.52686e-03, 0.019996e+09 },
{ 0.97685e-03, 0.029996e+09 },
{ 0.15923e-02, 0.039996e+09 },
{ 0.24090e-02, 0.049996e+09 },
{ 0.34674e-02, 0.059996e+09 },
{ 0.48120e-02, 0.069996e+09 },
{ 0.64921e-02, 0.079996e+09 },
{ 0.85618e-02, 0.089996e+09 },
{ 0.11080e-01, 0.099996e+09 },
{ 0.14110e-01, 0.109996e+09 },
{ 0.17723e-01, 0.119996e+09 },
{ 0.21991e-01, 0.129996e+09 },
{ 0.26994e-01, 0.139996e+09 },
{ 0.32819e-01, 0.149996e+09 },
{ 0.39556e-01, 0.159996e+09 },
{ 0.47301e-01, 0.169996e+09 },
{ 0.56159e-01, 0.179996e+09 },
{ 0.66236e-01, 0.189996e+09 },
{ 0.77648e-01, 0.199996e+09 },
{ 0.90516e-01, 0.209996e+09 },
{ 0.10497e+00, 0.219996e+09 },
{ 0.12114e+00, 0.229996e+09 },
{ 0.13916e+00, 0.239996e+09 },
{ 0.15919e+00, 0.249996e+09 },
{ 0.18138e+00, 0.259996e+09 },
{ 0.20588e+00, 0.269996e+09 },
{ 0.23287e+00, 0.279996e+09 },
{ 0.26252e+00, 0.289996e+09 },
{ 0.29502e+00, 0.299996e+09 },
{ 0.33054e+00, 0.309996e+09 },
{ 0.36929e+00, 0.319996e+09 },
{ 0.41147e+00, 0.329996e+09 },
{ 0.45729e+00, 0.339996e+09 },
{ 0.50696e+00, 0.349996e+09 },
{ 0.56073e+00, 0.359996e+09 },
{ 0.61881e+00, 0.369996e+09 },
{ 0.68145e+00, 0.379996e+09 },
{ 0.74890e+00, 0.389996e+09 },
{ 0.82142e+00, 0.399996e+09 },
{ 0.89928e+00, 0.409996e+09 },
{ 0.98274e+00, 0.419996e+09 },
{ 0.10721e+01, 0.429996e+09 } };
#define L 3.
#define B 0.1
#define H 0.1
static Vdouble corners[8][3] = {
{0.,0.,0.},
{L, 0.,0.},
{L, H ,0.},
{0.,H ,0.},
{0.,0.,B },
{L, 0.,B },
{L, H ,B },
{0.,H ,B } };
static Vint patch[8] = {1,2,3,4,5,6,7,8};
static Vint bctags[3] = {SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ};
static Vdouble v_max = 6.;
/*----------------------------------------------------------------------
Nonlinear Elasto-plastic Transient Solution Procedure
----------------------------------------------------------------------*/
int
main()
{
vis_MapMesh *mapmesh;
vis_Connect *connect;
vfe_PlasMat *plasmat;
vfe_MatlFun *matlfun;
vfe_Solid3D *solid3d;
vfs_DofTab *doftab;
vfs_SolProc *solproc;
vfs_SysMatrix *stiff, *mass;
vfs_SysVector *sysv, *react;
Vint n, numnp, numel, nhist, nintpnt, outputnode;
Vint nedofs, loc[24], tag[24];
Vint nix, ix[8], lm[24], aid;
Vint neq, nre;
Vint action, singular, stepid;
Vdouble elasprop[2], yield;
Vdouble *hist1, *hist2, *tmphist;
Vdouble coords[8][3];
Vdouble md[24];
Vdouble value, ul[24], fl[24], sl[300];
Vdouble massfactor, time;
Vdouble strs[24*6], strn[24*6];
/* Generate simple mesh */
mapmesh = vis_MapMeshBegin ();
vis_MapMeshDef (mapmesh,8,1);
for(n = 1; n <= 8; n++) {
vis_MapMeshSetPoint (mapmesh,n,corners[n-1]);
}
vis_MapMeshSetPatch (mapmesh,1,SYS_SHAPEHEX,2,0,0,patch);
vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID1,1,SYS_FACE,5,5);
vis_MapMeshSetPointAssoc (mapmesh,VIS_MISCID2,2,2);
vis_MapMeshSetPatchParami (mapmesh,1,VIS_MESH_SHAPE,SYS_SHAPEHEX);
vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMI,20);
vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMJ,5);
vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMK,1);
connect = vis_ConnectBegin ();
vis_ConnectPre (connect,SYS_DOUBLE);
vis_MapMeshGenerate (mapmesh,connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
/* Locate output node */
for(n = 1; n <= numnp; n++) {
vis_ConnectNodeAssoc (connect,VIS_MISCID2,1,&n,&aid);
if(aid == 2) {
outputnode = n;
break;
}
}
/* Construct PlasMat object using exam6.c data */
plasmat = vfe_PlasMatBegin ();
vfe_PlasMatSetParami (plasmat,VFE_STRAINTYPE,VFE_SMALLSTRAIN);
elasprop[0] = 2.1e+11;
elasprop[1] = 0.3;
vfe_PlasMatSetElasProp (plasmat,elasprop);
vfe_PlasMatDef (plasmat,PLASMAT_PIECEWISE);
yield = 0.91294e+08;
vfe_PlasMatSetPiecewise (plasmat,yield,44,hard);
vfe_PlasMatSetDensity (plasmat,7850.);
matlfun = vfe_MatlFunBegin ();
vfe_PlasMatMatlFun (plasmat,matlfun);
vfe_MatlFunNumHist(matlfun,&nhist);
/* Construct Solid3D object */
solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetTopology (solid3d,SYS_SHAPEHEX,2,0,0);
vfe_Solid3DSetParami (solid3d,VFE_TECH,VFE_TECH_MIXED);
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,(Vobject*)matlfun);
vfe_Solid3DNumIntPnt (solid3d,SYS_ANALYSIS_STRUCTURAL,&nintpnt);
vfe_Solid3DNumDof (solid3d,SYS_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,SYS_ANALYSIS_STRUCTURAL,loc,tag);
/* Allocate plasticity histories for steps N and N+1 */
hist1 = (Vdouble*)malloc(numel*nhist*nintpnt*sizeof(Vdouble));
hist2 = (Vdouble*)malloc(numel*nhist*nintpnt*sizeof(Vdouble));
memset (hist1,0,numel*nhist*nintpnt*sizeof(Vdouble));
memset (hist2,0,numel*nhist*nintpnt*sizeof(Vdouble));
/* Create DofTab object */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vfs_DofTabElemUse (doftab,n,nix,ix);
}
for(n = 1; n <= numnp; n++) {
vis_ConnectNodeAssoc (connect,VIS_MISCID1,1,&n,&aid);
if(aid != 5) continue;
vfs_DofTabNodePerm (doftab,n,3,bctags);
}
vfs_DofTabProcess (doftab,&neq,&nre);
/* Create stiffness matrix */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
for(n = 1; n <= numel; n++) {
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DInitHist (solid3d);
}
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixProcess (stiff);
/* Create and form mass matrix (it does not change) */
mass = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (mass,SYSMATRIX_DIAG,neq);
vfs_SysMatrixPreProcess (mass);
vfs_SysMatrixProcess (mass);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,coords);
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DMassDiag (solid3d,coords,md);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysMatrixAssemDiag (mass,nedofs,lm,md);
}
/* Instance SolProc object */
solproc = vfs_SolProcBegin ();
/* Define it for structural analyis for "neq" dofs */
vfs_SolProcDef (solproc,SYS_ANALYSIS_STRUCTURAL,neq,nre);
/* Specify a transient, linear analysis for 100 units of time */
vfs_SolProcSetStrategy (solproc,SYS_SOL_TRANSIENT,SOLPROC_NONLINEAR,0.1);
vfs_SolProcInitTimeStep (solproc,1.e-5);
vfs_SolProcSetParamd (solproc,SOLPROC_MINTIMESTEP,1.e-6);
/* Infinite loop for reverse communication */
while(vfs_SolProcExecute (solproc,&action),
action != VFS_ACTION_COMPLETE &&
action != VFS_ACTION_INCOMPLETE) {
/* initial conditions - u=0., v= linear with x */
if(action == VFS_ACTION_INITIAL) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U0,&sysv);
vfs_SysVectorInit (sysv,0.);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_V0,&sysv);
vfs_SysVectorInit (sysv,0.);
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,1,&bctags[1],lm);
vis_ConnectCoordsdv (connect,1,&n,coords);
value = v_max * coords[0][0] / L;
vfs_SysVectorScatter (sysv,1,lm,&value);
}
/* factor mass matrix: needed for initial acceleration, if nonzero */
} else if(action == VFS_ACTION_FACTORMASS) {
vfs_SysMatrixFactor (mass);
/* compute element reactions */
} else if(action == VFS_ACTION_REACT) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&sysv);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SysVectorInit (react,0.);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,coords);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
memset (ul,0,nedofs*sizeof(Vdouble));
vfs_SysVectorGather (sysv,nedofs,lm,ul);
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DReact (solid3d,coords,ul,fl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
}
/* solve for initial acceleration */
} else if(action == VFS_ACTION_SOLVEMASS) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&sysv);
vfs_SysMatrixSolve (mass,react,sysv);
/* form predictor */
} else if(action == VFS_ACTION_PREDICTOR) {
/*
* Nothing to be done here. This action is reserved for problems
* where displacements, velocities, or accelerations at at nodes
* are given as boundary conditions. In this case, those values
* should be set here
*/
/* compute nonlinear element reactions and stiffness */
} else if(action == VFS_ACTION_REACTSTIFF) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&sysv);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SysVectorInit (react,0.);
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,coords);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
memset (ul,0,nedofs*sizeof(Vdouble));
vfs_SysVectorGather (sysv,nedofs,lm,ul);
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DReactStiff (solid3d,coords,ul,SYS_ON,fl,sl);
vfs_SysVectorAssem (react,nedofs,lm,fl);
vfs_SysMatrixAssem (stiff,nedofs,lm,sl);
}
vfs_SolProcGetDouble (solproc,SOLPROC_MASSFACTOR,&massfactor);
vfs_SysMatrixSum (stiff,massfactor,mass);
/* compute linear stiffness */
} else if(action == VFS_ACTION_STIFF) {
vfs_SysMatrixInit (stiff,0.);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,coords);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
memset (ul,0,nedofs*sizeof(Vdouble));
vfs_SysVectorGather (sysv,nedofs,lm,ul);
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DReactStiff (solid3d,coords,ul,SYS_ON,fl,sl);
vfs_SysMatrixAssem (stiff,nedofs,lm,sl);
}
vfs_SolProcGetDouble (solproc,SOLPROC_MASSFACTOR,&massfactor);
vfs_SysMatrixSum (stiff,massfactor,mass);
/* compute inertia */
} else if(action == VFS_ACTION_INERTIA) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_A1,&sysv);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_INERTIA,&react);
vfs_SysMatrixProd (mass,sysv,react);
/* factor stiffness */
} else if(action == VFS_ACTION_FACTOR) {
vfs_SysMatrixFactor (stiff);
vfs_SysMatrixGetInteger (stiff,SYSMATRIX_SINGULAR,&singular);
if(singular) {
printf("Singular matrix detected... aborted\n");
break;
}
/* solve linear system */
} else if(action == VFS_ACTION_SOLVE) {
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_REACT,&react);
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_DELTAINC,&sysv);
vfs_SysMatrixSolve (stiff,react,sysv);
/* update */
} else if(action == VFS_ACTION_UPDATE) {
/*
* Nothing to be done here. This action is reserved for problems
* with large rotations that need to be updated by the user by
* procedures such as an exponential map
*/
/* end step */
} else if(action == VFS_ACTION_ENDSTEP) {
/* print corner displacements */
vfs_SolProcGetSysVector (solproc,VFS_VECTOR_U1,&sysv);
vfs_DofTabNodeDof (doftab,outputnode,3,bctags,lm);
memset (ul,0,3*sizeof(Vdouble));
vfs_SysVectorGather (sysv,3,lm,ul);
vfs_SolProcGetDouble (solproc,SOLPROC_TIME,&time);
vfs_SolProcGetInteger (solproc,SOLPROC_STEPID,&stepid);
printf("stepid= %5d time= %e y-displacement= %e\n",
stepid,time,ul[1]);
for(n = 1; n <= numel; n++) {
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,coords);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
memset (ul,0,nedofs*sizeof(Vdouble));
vfs_SysVectorGather (sysv,nedofs,lm,ul);
vfe_Solid3DSetHistPtr (solid3d,hist1+(n-1)*nhist*nintpnt,
hist2+(n-1)*nhist*nintpnt);
vfe_Solid3DStrsStrn (solid3d,coords,ul,strs,strn);
}
/* swap history variables for next step */
tmphist = hist1;
hist1 = hist2;
hist2 = tmphist;
} else {
printf("must code action above\n");
exit(0);
}
}
if(action == VFS_ACTION_COMPLETE) {
printf("Normal termination\n");
} else {
printf("Terminated with errors\n");
}
vfs_SolProcEnd (solproc);
vfs_SysMatrixEnd (stiff);
vfs_SysMatrixEnd (mass);
vfe_MatlFunEnd (matlfun);
vfe_PlasMatEnd (plasmat);
vfs_DofTabEnd (doftab);
vis_ConnectEnd (connect);
vis_MapMeshEnd (mapmesh);
vfe_Solid3DEnd (solid3d);
free (hist1);
free (hist2);
return 0;
}
