5. Examples

This chapter contains example VfxTools programs. The source code for all examples resides in devtools/src/vfx/exam.

5.1. Example 1, Attachment Point Load Distribution using GenLCase

This example illustrates using GenLCase to compute the equivalent concentrated loads on the nodes of a group of element faces due to a force and moment applied to a point in space. The point in space is referred to as an attachment point. The resulting equivalent nodal loads consist of translational forces only. The GenLCase object requires a Model object to be registered using vfx_GenLCaseSetObject(). Only the Connect object component of Model is required for the computaton of distributing loads. The function vfx_GenLCaseSetAttach() is required to be called once for each configuration of attachment point location and associated element face group. Any number of attachment point configurations may be defined, each assigned an identifying index. Once the attachment point configuration has been defined then load distributions may be computed repeatedly using vfx_GenLCaseAttach(). Note that the computed loads are added to the State object input to vfx_GenLCaseAttach(). The loads computed at the corner nodes of the faces of higher order (parabolic) elements may be counter-intuitive. In this example, using parabolic quadrilateral faces, the loads at the corner nodes are in the opposite direction and of smaller magnitude than the loads computed at the midside nodes. For triangular parabolic faces the loads at the corner nodes are zero.

#include <stdio.h>
#include "vfx/vfx.h"
/*
                  9----10----11----12----13
   y              |           |           |
   |              |           |           |
   --x            6           7           8
  /               |           |           |
 z                |           |           |
                  1-----2-----3-----4-----5      o (attachment point)
*/

#define MAX_ELEM  2
#define MAX_NODE 13

static Vint conn[MAX_ELEM][8] = {
   {1, 3,11, 9, 2, 7,10, 6},
   {3, 5,13,11, 4, 8,12, 7} };

static Vdouble coords[MAX_NODE][3] = {
   {0.,0.,0.},  {1.,0.,0.},  {2.,0.,0.},  {3.,0.,0.},  {4.,0.,0.},
   {0.,1.,0.},               {2.,1.,0.},               {4.,1.,0.},
   {0.,2.,0.},  {1.,2.,0.},  {2.,2.,0.},  {3.,2.,0.},  {4.,2.,0.} };

static Vdouble xp[3] = {
   5.,0.,0.};

static Vdouble force[2][3]  = { {1000.,0.,0.}, {500.,100.,  0.} };
static Vdouble moment[2][3] = { {   0.,0.,0.}, {  0.,  0.,200.} };

static void
print_LCase(vis_LCase *lcase, vis_Connect *connect, Vchar *title);

/*----------------------------------------------------------------------
                      Attachment Point Load Distribution using GenLCase
----------------------------------------------------------------------*/
int
main()
{
   Vint i;
   vis_Model *model;
   vis_Connect *connect;
   vis_Group *group;
   vis_LCase *lcase;
   vfx_GenLCase *genlcase;

                   /* create connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,MAX_NODE,MAX_ELEM);

                   /* set element node connectivity */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetTopology (connect,i+1,SYS_SHAPEQUAD,3,0,0);
      vis_ConnectSetElemNode (connect,i+1,conn[i]);
   }

                   /* set node coordinates */ 
   for(i = 0; i < MAX_NODE; i++) {
      vis_ConnectSetCoordsdv (connect,i+1,coords[i]);
   }
                   /* create model object */ 
   model = vis_ModelBegin();
   vis_ModelSetObject (model,VIS_CONNECT,connect);

                   /* create group of element faces/edges or nodes */ 
   group = vis_GroupBegin();

                   /* create genlcase object */ 
   genlcase = vfx_GenLCaseBegin();
   vfx_GenLCaseSetObject (genlcase,VIS_MODEL,model);

                   /* compute distributing forces on faces of elem 1 and 2 */ 
   vis_GroupDef (group,MAX_ELEM,SYS_ELEM,SYS_FACE);
   vis_GroupSetIndex (group,1,1);
   vis_GroupSetIndex (group,2,1);
   vfx_GenLCaseSetAttach (genlcase,1,xp,group);

                   /* create lcase object */ 
   lcase = vis_LCaseBegin();
                   /* compute two load distributions */
   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,1,force[0],moment[0],lcase);
   print_LCase (lcase,connect,"Distribute on faces, case 1");

   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,1,force[1],moment[1],lcase);
   print_LCase (lcase,connect,"Distribute on faces, case 2");

                   /* compute distributing forces on edges 1 and 2 of elem 2 */ 
   vis_GroupDef (group,MAX_ELEM,SYS_ELEM,SYS_EDGE);
   vis_GroupSetIndex (group,2,3);
   vfx_GenLCaseSetAttach (genlcase,2,xp,group);

                   /* compute two load distributions */
   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,2,force[0],moment[0],lcase);
   print_LCase (lcase,connect,"Distribute on edges, case 1");

   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,2,force[1],moment[1],lcase);
   print_LCase (lcase,connect,"Distribute on edges, case 2");

                   /* compute distributing forces on nodes */ 
   vis_GroupDef (group,MAX_NODE,SYS_NODE,SYS_NONE);
   vis_GroupSetIndex (group, 3,1);
   vis_GroupSetIndex (group, 4,1);
   vis_GroupSetIndex (group, 5,1);
   vis_GroupSetIndex (group, 8,1);
   vis_GroupSetIndex (group,13,1);
   vfx_GenLCaseSetAttach (genlcase,3,xp,group);

                   /* create 2 lcase objects */ 
   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,3,force[0],moment[0],lcase);
   print_LCase (lcase,connect,"Distribute on nodes, case 1");
   vis_LCaseClear(lcase);
   vfx_GenLCaseAttach (genlcase,3,force[1],moment[1],lcase);
   print_LCase (lcase,connect,"Distribute on nodes, case 2");

                   /* end objects */ 
   vis_ConnectEnd (connect);
   vis_ModelEnd (model);
   vis_GroupEnd (group);
   vis_LCaseEnd (lcase);
   vfx_GenLCaseEnd (genlcase);
   return 0;
}

static void
print_LCase(vis_LCase *lcase, vis_Connect *connect, Vchar *title)
{
   Vint i;
   Vint maxindex, flags;
   Vdouble v[3], x[3], f[3], m[3];

                   /* print forces and sum forces and moments to verify */ 
   printf("\n%s:\n\n",title);
   f[0] = f[1] = f[2] = 0.;
   m[0] = m[1] = m[2] = 0.;
   vis_LCaseConcMax (lcase,LCASE_FORCE,&maxindex);
   for(i = 1; i <= maxindex; i++) {
      vis_LCaseConcFlag (lcase,i,LCASE_FORCE,&flags);
      if(flags == 0)  continue;
      vis_LCaseConcdv (lcase,i,LCASE_FORCE,v);
      printf("i= %d, force= %e %e %e\n",i,v[0],v[1],v[2]);
      f[0] += v[0];
      f[1] += v[1];
      f[2] += v[2];
      vis_ConnectCoordsdv(connect,1,&i,&x);
      m[0] += v[2]*(x[1]-xp[1]) - v[1]*(x[2]-xp[2]);
      m[1] += v[0]*(x[2]-xp[2]) - v[2]*(x[0]-xp[0]);
      m[2] += v[1]*(x[0]-xp[0]) - v[0]*(x[1]-xp[1]);
   }
   printf("Total force              = %e %e %e\n",f[0],f[1],f[2]);
   printf("Total moment about attach= %e %e %e\n",m[0],m[1],m[2]);
}

5.2. Example 2, Statically Determinant Constraints using GenMCase

This example illustrates using GenMCase to generate a set of statically determinant multipoint constraints for a free body. The body is assumed to be a single, continuous, fully connected body. The GenMCase object requires a Model object to be registered using vfx_GenMCaseSetObject(). Only the Connect object component of Model is required for the computaton of determinant constraints. Use the function vfx_GenMCaseIRelief() to compute the multipoint constraints (MPC) and add them to the input MCase object. The resulting MPC’s are printed.

#include <stdio.h>
#include "vfx/vfx.h"
/*
                 4------5------6
                /.     /.     /|
  y            / .    / .    / |
  |           /  .   /  .   /  |
  --x        /   1../...2../...3
 /         10-----11-----12   /
z           |  .   |  .   |  /
            | .    | .    | /
            |.     |.     |/
            7------8------9
*/
#define MAX_ELEM  2
#define MAX_NODE 12

static Vint conn[MAX_ELEM][8] = {
   {1, 2, 5, 4, 7, 8,11,10},
   {2, 3, 6, 5, 8, 9,12,11} };

static Vdouble coords[MAX_NODE][3] = {
   {0.,0.,0.},  {1.,0.,0.},  {2.,0.,0.},
   {0.,1.,0.},  {1.,1.,0.},  {2.,1.,0.},
   {0.,0.,1.},  {1.,0.,1.},  {2.,0.,1.},
   {0.,1.,1.},  {1.,1.,1.},  {2.,1.,1.} };

/*----------------------------------------------------------------------
                      Statically Determinant Constraints using GenMCase
----------------------------------------------------------------------*/
int
main()
{
   Vint i;
   vis_Model *model;
   vis_Connect *connect;
   vis_MCase *mcase;
   vfx_GenMCase *genmcase;
   Vint maxindex, maxterms, maxrhs;
   Vint index, nterms;
   Vdouble r, c[3];
   Vint ig[3], ix[3];

                   /* create connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,MAX_NODE,MAX_ELEM);

                   /* set element node connectivity */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetTopology (connect,i+1,SYS_SHAPEHEX,2,0,0);
      vis_ConnectSetElemNode (connect,i+1,conn[i]);
   }

                   /* set node coordinates */ 
   for(i = 0; i < MAX_NODE; i++) {
      vis_ConnectSetCoordsdv (connect,i+1,coords[i]);
   }
                   /* create model object */ 
   model = vis_ModelBegin();
   vis_ModelSetObject (model,VIS_CONNECT,connect);

                   /* create genmcase object */ 
   genmcase = vfx_GenMCaseBegin();
   vfx_GenMCaseSetObject (genmcase,VIS_MODEL,model);

                   /* create mcase object */ 
   mcase = vis_MCaseBegin();
   vfx_GenMCaseIRelief (genmcase,NULL,mcase);

                   /* query MPC's */
   vis_MCaseMax (mcase,&maxindex,&maxterms,&maxrhs);
   printf("\nMultiPoint Constraints\n");
   for(index = 1; index <= maxindex; index++) {
      vis_MCaseNum (mcase,index,&nterms);
      if(nterms) {
         vis_MCaseMPCdv (mcase,index,&nterms,ix,ig,c,&r);
         printf("index= %d, nterms= %d, r= %f\n",index,nterms,r);
         for(i = 0; i < nterms; i++) {
            printf(" ix= %d, tag= %d, c= %f\n",ix[i],ig[i],c[i]);
         }
      }
   }
                   /* end objects */ 
   vis_ConnectEnd (connect);
   vis_ModelEnd (model);
   vis_MCaseEnd (mcase);
   vfx_GenMCaseEnd (genmcase);
   return 0;
}

5.3. Example 3, Solving Inertia Relief and Attachment Point Loading

This example illustrates using ProSolve with GenLCase and GenMCase to perform static linear analysis on a number of load cases. The analysis is typical of the techniques involved to solve for displacements and stresses given a series of attachment point loads generated from an up-stream mechanical motion analysis. The inertia relief feature of ProSolve is used to automatically compute inertia based loads to equilibrate the applied loads.

A simple 10x2x1 aluminum beam model is generated using 10 node tetrahedra. Two attachment points are located 1 unit from the center of each end face. Three load cases are solved for. The first load case is a force directed exactly along the axis of the beam. The remaining load cases are simply manufactured numbers for illustration only.

#include <stdlib.h>
#include "vfx/vfx.h"

                   /* number of attachment points */ 
#define NATTACH     2
                   /* number of load cases */ 
#define NLOADCASES  3

                   /* helper functions */ 
static void
createModel(vis_Model *model);
static void
printDisp(vis_State *state, Vint caseid);
static void
printStress(vis_State *state, Vint caseid);
                   /* monitor function */ 
static void
monitorfun(vfx_ProSolve *prosolve, Vobject *object);

                   /* attachment point locations */ 
static Vdouble xattach[NATTACH][3] = { {-1.,1.,0.5}, { 11.,1.,0.5 } };

                   /* attachment point forces and moments */ 
static Vdouble force[NLOADCASES][NATTACH][3] = 
                 { { { 1., 0., 0.}, { 0., 0., 0.} },
                   { {-2.,0.5, 3.}, { 5., 4., 3.} },
                   { { 1., 2., 7.}, {-9., 3., 0.} } };
static Vdouble moment[NLOADCASES][NATTACH][3] = 
                 { { { 0., 0., 0.}, { 0., 0., 0.} },
                   { {-4., 7., 8.}, { 5., 1., 6.} },
                   { { 2., 9., 2.}, {-3., 4., 2.} } };

                   /* solution result types */ 
static Vint res[2] = {SYS_RES_D,SYS_RES_S};

/*----------------------------------------------------------------------
  Solving multiple load cases with inertia relief and attachment points
----------------------------------------------------------------------*/
int
main()
{
   Vint          i, n;
   Vint          numel, numnp;
   vsy_HashTable *mchash,*lchash;
   vsy_List      *spropl;
   vsy_TextFun   *tftext,*tfhtml,*tftee;
   vsy_PlainText *plaintext;
   vsy_HTMLText  *htmltext;
   vsy_TextTee   *texttee;
   vis_Model     *model;
   vis_MCase     *mcase;
   vis_SProp     *sprop;
   vis_Group     *fgroup,*ngroup;
   vis_Connect   *connect;
   vis_LCase     *lcase;
   vis_GridFun   *gridfun;
   vis_State     *stated, *states;
   vis_RProp     *rpropd, *rprops;
   vfx_GenMCase  *genmcase;
   vfx_GenLCase  *genlcase;
   vfx_ProSolve  *prosolve;
   Vint init;

                   /* create model */
   model = vis_ModelBegin ();
   createModel (model);

                   /* get connect object and number of nodes and elements */ 
   vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
   vis_ConnectNumber (connect,SYS_ELEM,&numel);
   vis_ConnectNumber (connect,SYS_NODE,&numnp);

                   /* compute constraints and add MCase to Model */
   genmcase = vfx_GenMCaseBegin ();
   vfx_GenMCaseSetObject (genmcase,VIS_MODEL,model);
   mcase = vis_MCaseBegin ();
   vfx_GenMCaseIRelief (genmcase,NULL,mcase);
   vfx_GenMCaseEnd (genmcase);
   mchash = vsy_HashTableBegin ();
   vsy_HashTableInsert (mchash,1,(Vobject*)mcase);
   vis_ModelSetHashTable (model,VIS_MCASE,mchash);

                   /* define attachment point configurations in GenLCase */
   genlcase = vfx_GenLCaseBegin ();
   vfx_GenLCaseSetObject (genlcase,VIS_MODEL,model);
   fgroup = vis_GroupBegin ();
   vis_GroupDef (fgroup,numel,SYS_ELEM,SYS_FACE);
   ngroup = vis_GroupBegin ();
   vis_GroupDef (ngroup,numnp,SYS_NODE,SYS_NONE);

   for (i = 1; i <= NATTACH; i++) {
      vis_GroupClear (fgroup);
      vis_GroupClear (ngroup);
      vis_ConnectSetGroupParami (connect,CONNECT_ASSOCID,i);
      vis_ConnectSetGroupParami (connect,CONNECT_ASSOCTYPE,VIS_MISCID);
      vis_ConnectNodeGroup (connect,CONNECT_ASSOC,NULL,ngroup);
      vis_ConnectFaceGroup (connect,CONNECT_CONTAINED,ngroup,fgroup);
      vfx_GenLCaseSetAttach (genlcase,i,xattach[i-1],fgroup);
   }
   vis_GroupEnd (fgroup);
   vis_GroupEnd (ngroup);

                   /* build results states for displacement and stress */ 
   gridfun = vis_GridFunBegin ();
   vis_ConnectGridFun (connect,gridfun);
                   /* displacement at nodes */ 
   stated = vis_StateBegin ();
   vis_StatePre (stated,SYS_DOUBLE);
   rpropd = vis_RPropBegin ();
   vis_RPropDef (rpropd,SYS_NODE,SYS_NONE);
   vis_RPropSetType (rpropd,SYS_RES_D);
                   /* stress at element nodes */ 
   states = vis_StateBegin ();
   vis_StatePre (states,SYS_DOUBLE);
   vis_StateSetObject (states,VIS_GRIDFUN,gridfun);
   rprops = vis_RPropBegin ();
   vis_RPropDef (rprops,SYS_ELEM,SYS_NODE);
   vis_RPropSetType (rprops,SYS_RES_S);

                   /* initialize ProSolve */
   prosolve = vfx_ProSolveBegin ();
   init = 0;
   vfx_ProSolveSetFunction (prosolve,PROSOLVE_FUN_MONITOR,monitorfun,
                            (Vobject*)&init);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_PRINTLEVEL,0);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,SYS_OFF);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,SYS_OFF);

                   /* set console text and html object */ 
   plaintext = vsy_PlainTextBegin ();
   vsy_PlainTextConnectFile (plaintext,stdout);
   tftext = vsy_TextFunBegin ();
   vsy_PlainTextTextFun (plaintext,tftext);

   htmltext = vsy_HTMLTextBegin ();
   vsy_HTMLTextOpenFile (htmltext,"exam3.htm");
   tfhtml = vsy_TextFunBegin ();
   vsy_HTMLTextTextFun (htmltext,tfhtml);

                   /* link these formatters to a TextTee */
   texttee = vsy_TextTeeBegin ();
   vsy_TextTeeSetObject (texttee,VSY_TEXTFUN_1,tftext);
   vsy_TextTeeSetObject (texttee,VSY_TEXTFUN_2,tfhtml);
   tftee = vsy_TextFunBegin ();
   vsy_TextTeeTextFun (texttee,tftee);

                   /* set ProSolve Model and TextFun objects */ 
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);
   vfx_ProSolveSetObject (prosolve,VSY_TEXTFUN,tftee);
   vfx_ProSolveInit (prosolve);

                   /* execute each SProp */ 
   lchash = vsy_HashTableBegin ();
   vis_ModelSetHashTable (model,VIS_LCASE,lchash);
   spropl = vsy_ListBegin ();
   vis_ModelSetList (model,VIS_SPROP,spropl);
   for (n = 0; n < NLOADCASES; n++) {
      sprop = vis_SPropBegin ();
      vsy_ListInsert (spropl,n+1,(Vobject*)sprop);
      vis_SPropDef (sprop,SYS_SOL_STATIC);
      vis_SPropSetValuei (sprop,SPROP_INERTIARELIEF,SYS_ON);
      vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
      vis_SPropSetValuei (sprop,SPROP_LCASE,n+1);
      vis_SPropSetValuei (sprop,SPROP_MCASE,1);
      vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
      vis_SPropSetValuei (sprop,SPROP_RESPRINT_NUM,2);
      vis_SPropSetValueiv (sprop,SPROP_RESPRINT,res);
      lcase = vis_LCaseBegin ();
      vsy_HashTableInsert (lchash,n+1,(Vobject*)lcase);
      for (i = 0; i < NATTACH; i++) {
         vfx_GenLCaseAttach (genlcase,i+1,force[n][i],moment[n][i],lcase);
      }
      vfx_ProSolveRange (prosolve,n+1,n+1);
                   /* perform solution, this operation is compute intensive */ 
      vfx_ProSolveExec (prosolve);

                   /* check for abort */
      if(init == 2)  exit(0);

                   /* get displacement result and print */ 
      vfx_ProSolveState (prosolve,rpropd,stated);
      printDisp (stated,n+1);
                   /* get stress result and print */ 
      vfx_ProSolveState (prosolve,rprops,states);
      printStress (states,n+1);
   }
   vfx_GenLCaseEnd (genlcase);

                   /* delete output objects */ 
   vsy_PlainTextEnd (plaintext);
   vsy_HTMLTextEnd  (htmltext);
   vsy_TextTeeEnd   (texttee);
   vsy_TextFunEnd   (tftext);
   vsy_TextFunEnd   (tfhtml);
   vsy_TextFunEnd   (tftee);

                   /* delete Model object contents and object */ 
   vis_ModelDelete (model);
   vis_ModelEnd (model);

                   /* delete other objects and ProSolve object */ 
   vis_GridFunEnd (gridfun);
   vis_RPropEnd   (rpropd);
   vis_RPropEnd   (rprops);
   vis_StateEnd   (stated);
   vis_StateEnd   (states);
   vfx_ProSolveEnd (prosolve);
   return 0;
}

/*----------------------------------------------------------------------
             Function to monitor ProSolve execution
----------------------------------------------------------------------*/
static void
monitorfun(vfx_ProSolve *prosolve, Vobject *object)
{
   Vint  phase;
   Vint *init = (Vint*)object;
   Vlong fm;

   vfx_ProSolveGetInteger (prosolve,PROSOLVE_PHASE,&phase);
   if(*init == 1 && phase != SYS_PHASE_SYMFACTOR) {
      vfx_ProSolveGetLong (prosolve,PROSOLVE_FACTOR_NONZ,&fm);

      /* flag abort has been called */
      *init = 2;
      vfx_ProSolveAbort (prosolve);
   }

   if(phase == SYS_PHASE_MASS) {
   } else if(phase == SYS_PHASE_SYMFACTOR) {
      printf("Phase: SYS_PHASE_SYMFACTOR\n");
      *init = 1;
   } else if(phase == SYS_PHASE_INITCOND) {
      printf("Phase: SYS_PHASE_INITCOND\n");
   } else if(phase == SYS_PHASE_PREDICTOR) {
      printf("Phase: SYS_PHASE_PREDICTOR\n");
   } else if(phase == SYS_PHASE_STIFF) {
      printf("Phase: SYS_PHASE_STIFF\n");
   } else if(phase == SYS_PHASE_FACTOR) {
      printf("Phase: SYS_PHASE_FACTOR\n");
   } else if(phase == SYS_PHASE_UPDATE) {
      printf("Phase: SYS_PHASE_UPDATE\n");
   } else if(phase == SYS_PHASE_REACT) {
      printf("Phase: SYS_PHASE_REACT\n");
   } else if(phase == SYS_PHASE_REACTSTIFF) {
      printf("Phase: SYS_PHASE_REACTSTIFF\n");
   } else if(phase == SYS_PHASE_SOLVE) {
      printf("Phase: SYS_PHASE_SOLVE\n");
   } else if(phase == SYS_PHASE_GEOMSTIFF) {
      printf("Phase: SYS_PHASE_GEOMSTIFF\n");
   } else if(phase == SYS_PHASE_BUCKLING) {
      printf("Phase: SYS_PHASE_BUCKLING\n");
   } else if(phase == SYS_PHASE_VIBRATION) {
      printf("Phase: SYS_PHASE_VIBRATION\n");
   } else if(phase == SYS_PHASE_OUTPUT) {
      printf("Phase: SYS_PHASE_OUTPUT\n");
   }
}

/*----------------------------------------------------------------------
             Function to create a simple beam model
----------------------------------------------------------------------*/
static void
createModel(vis_Model *model)
{
   Vint i;
   Vint numel;
   vis_MapMesh   *mapmesh;
   vis_Connect   *connect;
   vis_MProp     *mprop;
   vis_EProp     *eprop;
   vsy_HashTable *mphash,*ephash;
   Vint    patch[8] = {1,2,3,4,5,6,7,8};
   Vdouble x[8][3] = { {0.,0.,0.}, {10.,0.,0.}, {10.,2.,0.}, {0.,2.,0.},
                       {0.,0.,1.}, {10.,0.,1.}, {10.,2.,1.}, {0.,2.,1.} };

                   /* generate mesh */ 
   mapmesh = vis_MapMeshBegin ();
   vis_MapMeshDef (mapmesh,8,1);
   for (i = 0; i < 8; i++) {
      vis_MapMeshSetPoint (mapmesh,i+1,x[i]);
   }
   vis_MapMeshSetPatch (mapmesh,1,VIS_SHAPEHEX,2,0,0,patch);
   vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID,1,SYS_FACE,5,1);
   vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID,1,SYS_FACE,6,2);
   vis_MapMeshSetParami (mapmesh,VIS_MESH_MAXI,3);
   vis_MapMeshSetPatchParami (mapmesh,1,VIS_MESH_SHAPE,VIS_SHAPETET);
   vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMI,10);
   vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMJ,2);
   vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMK,2);
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,0,0);
   vis_MapMeshGenerate (mapmesh,connect);
   vis_ConnectKernel (connect,0);
   vis_MapMeshEnd (mapmesh);
   vis_ConnectNumber (connect,SYS_ELEM,&numel);
   for (i = 1; i <= numel; i++) {
      vis_ConnectSetElemAssoc (connect,VIS_PROPID,i,1);
   }
   vis_ModelSetObject (model,VIS_CONNECT,connect);

                   /* generate property */ 
   ephash = vsy_HashTableBegin ();
   eprop = vis_EPropBegin ();
   vis_EPropDef (eprop,VIS_ELEM_SOLID);
   vis_EPropSetValuei (eprop,EPROP_MID,1);
   vsy_HashTableInsert (ephash,1,eprop);
   vis_ModelSetHashTable (model,VIS_EPROP,ephash);

                   /* generate material */ 
   mphash = vsy_HashTableBegin ();
   mprop = vis_MPropBegin ();
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   vis_MPropSetValued (mprop,MPROP_E,1.e+7);
   vis_MPropSetValued (mprop,MPROP_NU,.3);
   vis_MPropSetValued (mprop,MPROP_DENSITY,0.0000133);
   vsy_HashTableInsert (mphash,1,(Vobject*)mprop);
   vis_ModelSetHashTable (model,VIS_MPROP,mphash);
}

/*----------------------------------------------------------------------
             Function to print displacements
----------------------------------------------------------------------*/
static void
printDisp(vis_State *state, Vint caseid)
{
   Vint n,nument,enttype,subtype,datatype;
   Vdouble d[6];

   vis_StateInq (state,&nument,&enttype,&subtype,&datatype);
   printf("Displacements for load case %d:\n",caseid);

   for (n = 1; n <= nument; n++) {
      vis_StateDatadv (state,1,&n,d);
      printf("%8d %11.4e %11.4e %11.4e\n",n,d[0],d[1],d[2]);
   }
}

/*----------------------------------------------------------------------
             Function to print von mises stresses
----------------------------------------------------------------------*/
static void
printStress(vis_State *state, Vint caseid)
{
   Vint i,n,nument,enttype,subtype,datatype;
   Vdouble d[10];

   vis_StateInq (state,&nument,&enttype,&subtype,&datatype);
   vis_StateSetDerive (state,VIS_TENSOR_VONMISES);
   printf("Von Mises stresses for load case %d:\n",caseid);

   for (n = 1; n <= nument; n++) {
      vis_StateDatadv (state,1,&n,d);
      printf("Element = %d\n",n);
      for (i = 0; i < 10; i++) {
         printf("   %11.4e\n",d[i]);
      }
   }
}

5.4. Example 4, Design Analysis Finite Element Solver

This example illustrates using ProSolve to implement a complete design analysis finite element solver.

#include <stdlib.h>
#include "base/base.h"
#include "base/vututil.h"
#include "vfx/prosolve.h"
#include "vfx/logfile.h"
#include "vfx/prorender.h"

static void
vki_Term(void);
static void
render_model(vis_Model *model, vsy_TextFun *textfun, Vchar *htm);
static void
render_state(vfx_ProSolve *prosolve, vis_Model *model,
             vsy_TextFun *textfun, Vchar *htm);
static void
checksingular(vfx_ProSolve *prosolve);

vfx_ProSolve  *prosolve  = NULL;
vdm_DataFun   *datafun   = NULL;
vdm_NatLib    *natlib    = NULL;
vdm_NASFil    *nasfil    = NULL;
vdm_ABAFil    *abafil    = NULL;
vdm_ANSFil    *ansfil    = NULL;
vdm_LMan      *lman      = NULL;
vis_Model     *model     = NULL;

vsy_HTMLText  *htmltext  = NULL;
vsy_PlainText *plaintext = NULL;
vsy_PlainText *tty       = NULL;
vsy_TextFun   *tf0       = NULL;
vsy_TextFun   *tf1       = NULL;
vsy_TextFun   *tf2       = NULL;
vsy_TextFun   *tf3       = NULL;
vsy_TextTee   *texttee   = NULL;

vsy_TextFun   *textfun   = NULL;
vfx_LogFile   *logfile   = NULL;

/*----------------------------------------------------------------------
                      Design Analysis Finite Element Solver
----------------------------------------------------------------------*/
int
main(int argc, char *argv[])
{
   Vint  i;
   Vint  filetype;
   Vchar job[256],dat[256],log[256],out[256],foc[256],htm[256];
   Vchar res[256], ver[256], moc[256];
   Vint  sol=0,ord=PROSOLVE_MATRIXORDER_METIS;
   Vint  usr=0,prs=1,prl=1,prh=0,mem=256,ooc=0,lmp=0, lbs=7;
   Vint  mpc=VFX_MPCTYPE_TRANSFORM,prt=1,nas=0,llv=2;
   Vdouble cf1 = 5., cf2 = 1.7;
   Vint  sng=0;
   Vint  savemodel=0;
   Vint  restype;
   Vint svm=1, svr=1;
   Vdouble utl=5.e-7,ftl=5.e-3;
   Vint slt=SYS_SOLVERTYPE_MFP;
   Vint npc=0;
   Vint egs=SYS_EIGEN_LANCZOS;

                   /* initialize file name strings */
   sprintf(job,"chexa");
   dat[0] = '\0';
   log[0] = '\0';
   res[0] = '\0';
   out[0] = '\0';
   foc[0] = '\0';
   moc[0] = '\0';
   htm[0] = '\0';
   res[0] = '\0';

                   /* parse command line options */ 
   for(i = 1; i < argc;) {

                   /* job name */
      if(strcmp(argv[i],"-job") == 0) {
         sprintf(job,"%s",argv[i+1]);

                   /* equation solver type */ 
      } else if(strcmp(argv[i],"-sol") == 0) {
         sscanf(argv[i+1],"%d",&sol);
                   /* equation ordering type */ 
      } else if(strcmp(argv[i],"-ord") == 0) {
         sscanf(argv[i+1],"%d",&ord);
                   /* user flag */ 
      } else if(strcmp(argv[i],"-usr") == 0) {
         sscanf(argv[i+1],"%d",&usr);
                   /* log level */
      } else if(strcmp(argv[i],"-llv") == 0) {
         sscanf(argv[i+1],"%d",&llv);
                   /* printout level */
      } else if(strcmp(argv[i],"-prt") == 0) {
         sscanf(argv[i+1],"%d",&prt);
                   /* out-of-core megabytes */
      } else if(strcmp(argv[i],"-mem") == 0) {
         sscanf(argv[i+1],"%d",&mem);
                   /* HTML file name */
      } else if(strcmp(argv[i],"-htm") == 0) {
         sprintf(htm,"%s",argv[i+1]);
                   /* input data file name */
      } else if(strcmp(argv[i],"-dat") == 0) {
         sprintf(dat,"%s",argv[i+1]);
                   /* output res file name */
      } else if(strcmp(argv[i],"-res") == 0) {
         sprintf(res,"%s",argv[i+1]);
                   /* printout file name */
      } else if(strcmp(argv[i],"-out") == 0) {
         sprintf(out,"%s",argv[i+1]);
                   /* log file name */
      } else if(strcmp(argv[i],"-log") == 0) {
         sprintf(log,"%s",argv[i+1]);
                   /* out-of-core stiffness factor file */
      } else if(strcmp(argv[i],"-foc") == 0) {
         sprintf(foc,"%s",argv[i+1]);
                   /* out-of-core mass factor file */
      } else if(strcmp(argv[i],"-moc") == 0) {
         sprintf(moc,"%s",argv[i+1]);

                   /* save model, save results */
      } else if(strcmp(argv[i],"-svm") == 0) {
         sscanf(argv[i+1],"%d",&svm);
      } else if(strcmp(argv[i],"-svr") == 0) {
         sscanf(argv[i+1],"%d",&svr);
                   /* set low memory preconditioner */
      } else if(strcmp(argv[i],"-lmp") == 0) {
         sscanf(argv[i+1],"%d",&lmp);
                   /* Lanczos block size */
      } else if(strcmp(argv[i],"-lbs") == 0) {
         sscanf(argv[i+1],"%d",&lbs);
                   /* set out-of-core flag*/
      } else if(strcmp(argv[i],"-ooc") == 0) {
         sscanf(argv[i+1],"%d",&ooc);
                   /* external solver */
                   /* stdout log output */
      } else if(strcmp(argv[i],"-prs") == 0 || 
                strcmp(argv[i],"-std") == 0) {
         sscanf(argv[i+1],"%d",&prs);
                   /* plain text log file output */
      } else if(strcmp(argv[i],"-prl") == 0) {
         sscanf(argv[i+1],"%d",&prl);
                   /* HTML log file output */
      } else if(strcmp(argv[i],"-prh") == 0) {
         sscanf(argv[i+1],"%d",&prh);
                   /* set number of processors */
      } else if(strcmp(argv[i],"-npc") == 0) {
         sscanf(argv[i+1],"%d",&npc);
                   /* set nastran style results calculations */
      } else if(strcmp(argv[i],"-nas") == 0) {
         sscanf(argv[i+1],"%d",&nas);
                   /* set mpc approach */
      } else if(strcmp(argv[i],"-mpc") == 0) {
         sscanf(argv[i+1],"%d",&mpc);
         if(mpc != VFX_MPCTYPE_PENALTY && mpc != VFX_MPCTYPE_LAGMULT &&
            mpc != VFX_MPCTYPE_AUGMENT && mpc != VFX_MPCTYPE_TRANSFORM) {
            printf("Invalid value of \"mpc\". Reset to TRANSFORM\n");
            printf("Valid values are: 1 (penalty)\n");
            printf("                  2 (lagrange multiplier)\n");
            printf("                  3 (augmented lagrangian)\n");
            printf("                  4 (transform)\n");
            mpc = VFX_MPCTYPE_TRANSFORM;
         }
                   /* amls cf1 */
      } else if(strcmp(argv[i],"-cf1") == 0) {
         sscanf(argv[i+1],"%le",&cf1);
                   /* amls cf2 */
      } else if(strcmp(argv[i],"-cf2") == 0) {
         sscanf(argv[i+1],"%le",&cf2);
                   /* u tolerance */
      } else if(strcmp(argv[i],"-utl") == 0) {
         sscanf(argv[i+1],"%le",&utl);
                   /* f tolerance */
      } else if(strcmp(argv[i],"-ftl") == 0) {
         sscanf(argv[i+1],"%le",&ftl);
                   /* f tolerance */
      } else if(strcmp(argv[i],"-sng") == 0) {
         sscanf(argv[i+1],"%d",&sng);
      } else if(strcmp(argv[i],"-slt") == 0) {
         sscanf(argv[i+1],"%d",&slt);
         if(slt != SYS_SOLVERTYPE_LL && slt != SYS_SOLVERTYPE_MF &&
            slt != SYS_SOLVERTYPE_MFP && slt != SYS_SOLVERTYPE_MUMPS &&
            slt != SYS_SOLVERTYPE_PARDISO) {
            printf("Invalid value of \"slt\". Reset to SYS_SOLVERTYPE_LL\n");
            slt = SYS_SOLVERTYPE_LL;
         }
      } else if(strcmp(argv[i],"-egs") == 0) {
         sscanf(argv[i+1],"%d",&egs);
         if(egs != SYS_EIGEN_LANCZOS && egs != SYS_EIGEN_SUBSPACE &&
            egs != SYS_EIGEN_AMLS) {
            printf("Invalid value of \"egs\". Reset to SYS_EIGEN_LANCZOS\n");
            egs = SYS_EIGEN_LANCZOS;
         }
                   /* unrecognized */ 
      } else {
         printf("Unrecognized switch: %s\n",argv[i]);
         return 1;  
      }
      i += 2;
   }
                   /* check for invalid solver type for out-of-core */
   if(ooc && slt == SYS_SOLVERTYPE_LL) {
      printf("Out-of-core solver cannot be used with SYS_SOLVERTYPE_LL\n");
      return 1;
   }
                   /* set file names, defaults to job name */ 
   if(dat[0] == '\0') {
      sprintf(dat,"%s.dat",job);
   }
   if(res[0] == '\0') {
      sprintf(res,"%s.vdm",job);
   }
   if(out[0] == '\0') {
      sprintf(out,"%s.out",job);
   }
   if(log[0] == '\0') {
      sprintf(log,"%s.log",job);
   }
   if(foc[0] == '\0') {
      sprintf(foc,"%s.foc",job);
   }
   if(moc[0] == '\0') {
      sprintf(moc,"%s.moc",job);
   }
   if(htm[0] == '\0') {
      sprintf(htm,"%s.htm",job);
   }
                   /* text stream to the log files */
   logfile = vfx_LogFileBegin ();
   texttee = vsy_TextTeeBegin ();
   textfun = vsy_TextFunBegin ();
   vsy_TextTeeTextFun (texttee,textfun);
   vfx_LogFileSetObject(logfile,VSY_TEXTFUN,textfun);
                   /* set error handler to stream to log file */
   vfx_LogFileSetError (logfile);
   tf0 = vsy_TextFunBegin ();
   vsy_TextFunAPI(tf0,TEXTFUN_APIRETURN);
   vsy_TextTeeSetObject (texttee,VSY_TEXTFUN,tf0);
                   /* plain text file */ 
   if(prl) {
      plaintext = vsy_PlainTextBegin ();
      tf1 = vsy_TextFunBegin ();
      vsy_PlainTextTextFun (plaintext,tf1);
      vsy_TextFunOpenFile (tf1,log);
      if(vsy_TextFunError(tf1)) {
         return 1;
      }
      vsy_TextTeeSetObject (texttee,VSY_TEXTFUN_1,tf1);
   }
                   /* HTML text file */ 
   if(prh) {
      htmltext = vsy_HTMLTextBegin ();
      tf2 = vsy_TextFunBegin ();
      vsy_HTMLTextTextFun (htmltext,tf2);
      vsy_TextFunOpenFile(tf2,htm);
      if(vsy_TextFunError(tf2)) {
         return 1;
      }
      vsy_TextTeeSetObject (texttee,VSY_TEXTFUN_2,tf2);
   }
                   /* plain text to standard output */ 
   if(prs) {
      tty = vsy_PlainTextBegin ();
      tf3 = vsy_TextFunBegin ();
      vsy_PlainTextTextFun (tty,tf3);
      vsy_TextFunConnectFile (tf3,stdout);
      vsy_TextTeeSetObject (texttee,VSY_TEXTFUN_3,tf3);
   }

                   /* top of log file */ 
   vsy_TextFunString (textfun,"Log File");
   vsy_TextFunFormAction (textfun,SYS_TEXT_LINEBREAK);
                   /* output DevTools version */
   sprintf(ver,"Version %d",SYS_VERSION);
   vsy_TextFunString (textfun,ver);
   vsy_TextFunFormAction (textfun,SYS_TEXT_LINEBREAK);

                   /* log file, echo input parameters */
   vfx_LogFileTableInit (logfile,"Command Summary");
   vfx_LogFileTableEntryc (logfile,"(-job) job name",job);
   vfx_LogFileTableEntryc (logfile,"(-dat) data file",dat);
   vfx_LogFileTableEntryc (logfile,"(-res) res file",res);
   vfx_LogFileTableEntryc (logfile,"(-out) out file",out);
   vfx_LogFileTableEntryc (logfile,"(-log) log file",log);
   vfx_LogFileTableEntryc (logfile,"(-htm) html file",htm);
   vfx_LogFileTableEntryc (logfile,"(-foc) foc file",foc);
   vfx_LogFileTableEntryc (logfile,"(-moc) moc file",moc);
   vfx_LogFileTableEntryi (logfile,"(-prs) log stdout flag",prs);
   vfx_LogFileTableEntryi (logfile,"(-prl) log file flag",prl);
   vfx_LogFileTableEntryi (logfile,"(-prh) log html flag",prh);
   vfx_LogFileTableEntryi (logfile,"(-sol) solver",sol);
   vfx_LogFileTableEntryi (logfile,"(-ord) matrix ordering",ord);
   vfx_LogFileTableEntryi (logfile,"(-mem) ooc factor memory (MB)",mem);
   vfx_LogFileTableEntryi (logfile,"(-prt) print level",prt);
   vfx_LogFileTableEntryi (logfile,"(-lmp) low memory preconditioner",lmp);
   vfx_LogFileTableEntryi (logfile,"(-lbs) Lanczos block size",lbs);
   vfx_LogFileTableEntryi (logfile,"(-npc) Processors",npc);
   vfx_LogFileTableEntryi (logfile,"(-mpc) MPC approach",mpc);
   vfx_LogFileTableEntryi (logfile,"(-ooc) Out-of-core flag",ooc);
   vfx_LogFileTableEntryi (logfile,"(-slt) Solver type",slt);
   vfx_LogFileTableEntryi (logfile,"(-egs) Eigensolver type",egs);
   if(sol == 5) {
      vfx_LogFileTableEntryd (logfile,"(-utl) U tolerance",utl);
   }
   if(sol == 5 || sol == 6) {
      vfx_LogFileTableEntryd (logfile,"(-ftl) F tolerance",ftl);
   }
   vfx_LogFileTableEntryd (logfile,"(-cf1) AMLS factor 1",cf1);
   vfx_LogFileTableEntryd (logfile,"(-cf2) AMLS factor 2",cf2);
   vfx_LogFileTableTerm (logfile);

                   /* Access input file */
   datafun = vdm_DataFunBegin ();
                   /* Native .vdm file */
   if(strstr(dat,".vdm") != NULL) {
      filetype = VDM_NATIVE;
      natlib = vdm_NatLibBegin ();
      vdm_NatLibDataFun (natlib,datafun);
      if(strcmp(res,dat)!=0) {
         savemodel = 1;
      }
                   /* ABAQUS inp file */
   } else if(strstr(dat,".inp") != NULL) {
      filetype = VDM_ABAQUS_INPUT;
      abafil = vdm_ABAFilBegin ();
      vdm_ABAFilDataFun (abafil,datafun);
      savemodel = 1;
                   /* ANSYS cdb file */
   } else if(strstr(dat,".cdb") != NULL ||
             strstr(dat,".ans") != NULL) {
      filetype = VDM_ANSYS_INPUT;
      ansfil = vdm_ANSFilBegin ();
      vdm_ANSFilDataFun (ansfil,datafun);
      savemodel = 1;
                   /* NASTRAN Bulk Data file */
   } else if(strstr(dat,".bdf") != NULL ||
             strstr(dat,".dat") != NULL ||
             strstr(dat,".nas") != NULL ||
             strstr(dat,".fem") != NULL ||
             strstr(dat,".DAT") != NULL) {
      filetype = VDM_NASTRAN_BULKDATA;
      nasfil = vdm_NASFilBegin ();
      vdm_NASFilDataFun (nasfil,datafun);
      savemodel = 1;
   } else {
      vfx_LogFileTextString (logfile,"Error: Unrecognized input file type");
      vki_Term();
      return 1;
   }
   vdm_DataFunSetConvention (datafun,VDM_CONVENTION_DOUBLE);
   vdm_DataFunOpen (datafun,0,dat,filetype);
   if(vdm_DataFunError(datafun)) {
      vfx_LogFileTextString (logfile,"Error: Unable to open input data file");
      vki_Term();
      return 1;
   }
                   /* create lman as helper object to read model */
   lman = vdm_LManBegin ();
   vdm_LManSetObject (lman,VDM_DATAFUN,datafun);

                   /* load model */
   model = vis_ModelBegin ();
   vdm_LManLoadModel (lman,model);
   if(vdm_LManError(lman)) {
      vfx_LogFileTextString (logfile,"Error: Unable to load model information");
      vki_Term();
      return 1;
   }
                   /* close input file */ 
   vdm_LManEnd (lman);
   lman = NULL;
   vdm_DataFunClose (datafun);
   vdm_DataFunEnd (datafun);
   datafun = NULL;

                   /* delete input file objects */ 
   if(filetype == VDM_NATIVE) {
      vdm_NatLibEnd (natlib);
      natlib = NULL;
   } else if(filetype == VDM_ABAQUS_INPUT) {
      vdm_ABAFilEnd (abafil);
      abafil = NULL;
   } else if(filetype == VDM_ANSYS_INPUT) {
      vdm_ANSFilEnd (ansfil);
      ansfil = NULL;
   } else if(filetype == VDM_NASTRAN_BULKDATA) {
      vdm_NASFilEnd (nasfil);
      nasfil = NULL;
   }
                   /* process output file */
   if(strlen(res)==0) {
      restype = VDM_NATIVE;
   } else if(strstr(res,".op2")!=NULL) {
      restype = VDM_NASTRAN_OUTPUT2;
   } else if(strstr(res,".fil")!=NULL) {
      restype = VDM_ABAQUS_FIL;
   } else if(strstr(res,".vdm")!=NULL) {
      restype = VDM_NATIVE;
   } else if(strstr(res,".odb")!=NULL) {
      restype = VDM_ABAQUS_ODB;
   } else if(strstr(res,".vh5")!=NULL) {
      restype = VDM_NATIVE_HDF5;
   }
                   /* render model to html based log file */
   if(prh) {
      render_model (model,textfun,job);
   }
                   /* initialize ProSolve */
   prosolve = vfx_ProSolveBegin ();
   if(svm) {
      vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,savemodel);
      vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,svr);
   } else {
      vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,0);
      vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,0);
   }
   vfx_ProSolveSetParamd (prosolve,PROSOLVE_AMLS_CUTFACT1,cf1);
   vfx_ProSolveSetParamd (prosolve,PROSOLVE_AMLS_CUTFACT2,cf2);
   vfx_ProSolveSetObject (prosolve,VSY_TEXTFUN,textfun);

                   /* set parameters */
   vfx_ProSolveSetParami (prosolve,PROSOLVE_ITEREXPLEVEL,3);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_ITERMEMLEVEL,lmp);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_LANCZOS_BLOCKSIZE,lbs);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SOLVER,sol);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_FACTOR_OOC,ooc);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_UNROLL,1);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_MATRIXORDER,ord);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_USERFLAG,usr);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_PRINTLEVEL,prt);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_LOGLEVEL,llv);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_OOCMEM_FACTOR,mem);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_MPCTYPE,mpc);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_RESTYPE,restype);
   vfx_ProSolveSetParamc (prosolve,PROSOLVE_RESFILE,res);
   vfx_ProSolveSetParamc (prosolve,PROSOLVE_OUTFILE,out);
   vfx_ProSolveSetParamc (prosolve,PROSOLVE_KFACTOR_FILE,foc);
   vfx_ProSolveSetParamc (prosolve,PROSOLVE_MFACTOR_FILE,moc);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_NASTRANOUTPUT,nas);
   vfx_ProSolveSetParamd (prosolve,PROSOLVE_ITER_UTOL,utl);
   vfx_ProSolveSetParamd (prosolve,PROSOLVE_ITER_FTOL,ftl);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SOLVERTYPE,slt);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_EIGEN_METHOD,egs);
   if(npc) {
      vfx_ProSolveSetNumThreads (prosolve,npc);
   }

   /* set model */
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);

   /* check for singularity */
   if(sng) {
      checksingular (prosolve);
   }

   /* execute */
   vfx_ProSolveExec (prosolve);
   if(vfx_ProSolveError(prosolve)) {
      vki_Term();
      return 1;
   }
                   /* render results to html based log file */
   if(prh) {
      render_state (prosolve,model,textfun,job);
   }
                   /* close log file text function interfaces */ 
   if(prl) {
      vsy_TextFunCloseFile (tf1);
   }
   if(prh) {
      vsy_TextFunCloseFile (tf2);
   }
   if(prs) {
      vsy_TextFunDisconnectFile (tf3);
   }
   vki_Term();
   return 0;
}

/*----------------------------------------------------------------------
                      delete objects
----------------------------------------------------------------------*/
static void
vki_Term(void)
{
   if(prosolve ) vfx_ProSolveEnd(prosolve);
   if(datafun  ) vdm_DataFunEnd(datafun);
   if(natlib   ) vdm_NatLibEnd(natlib);
   if(nasfil   ) vdm_NASFilEnd(nasfil);
   if(abafil   ) vdm_ABAFilEnd(abafil);
   if(ansfil   ) vdm_ANSFilEnd(ansfil);
   if(lman     ) vdm_LManEnd(lman);
   if(model    ) {
      vis_ModelDelete (model);
      vis_ModelEnd(model);
   }
   if(htmltext ) vsy_HTMLTextEnd(htmltext);
   if(plaintext) vsy_PlainTextEnd(plaintext);
   if(tty      ) vsy_PlainTextEnd(tty);
   if(texttee  ) vsy_TextTeeEnd(texttee);
   if(tf0      ) vsy_TextFunEnd(tf0);
   if(tf1      ) vsy_TextFunEnd(tf1);
   if(tf2      ) vsy_TextFunEnd(tf2);
   if(tf3      ) vsy_TextFunEnd(tf3);
   if(textfun  ) vsy_TextFunEnd(textfun);
   if(logfile  ) vfx_LogFileEnd(logfile);
}

/*----------------------------------------------------------------------
                      render model and results to log file
----------------------------------------------------------------------*/
static void
render_model(vis_Model *model, vsy_TextFun *textfun, Vchar *htm)
{
   Vint i;
   vfx_ProRender *prorender;
   Vchar gifb[256], gifc[256];
   Vint glen;
                   /* create file names for image */ 
   strcpy(gifc,htm);
   glen = (Vint)strlen(gifc);
   for(i = glen-1; i >= 0; i--) {
      if(htm[i] == '/'  ||  htm[i] == '\\') {
         break;
      } else if(i == 0) {
         i = -1;
         break;
      }
   }
   strcpy(gifb,&gifc[i+1]);
   strcat(gifb,"_mod.gif");
   strcpy(&gifc[i+1],gifb);

                   /* instance renderer */ 
   prorender = vfx_ProRenderBegin();
   vfx_ProRenderSetDimension (prorender,240,180);
   vfx_ProRenderSetObject (prorender,VIS_MODEL,model);

                   /* generate model retained graphics data */ 
   vfx_ProRenderModel (prorender);

                   /* render graphics data to image file */ 
   vfx_ProRenderSetParamc (prorender,PRORENDER_OUTFILE,gifc);
   vfx_ProRenderSetLookAt (prorender,1.,1.,1.,0.,0.,0.,0.,0.,1.);
   vfx_ProRenderDraw (prorender);
   vfx_ProRenderEnd (prorender);

                   /* place image in log file */ 
   vsy_TextFunHeadingLevel (textfun,4);
   vsy_TextFunFormInit (textfun,SYS_TEXT_HEADING);
   vsy_TextFunString (textfun,"Initial Model Rendering");
   vsy_TextFunFormTerm (textfun,SYS_TEXT_HEADING);
   vsy_TextFunImageFile (textfun,gifb);
   vsy_TextFunFormAction (textfun,SYS_TEXT_LINEBREAK);
}

static void
render_state(vfx_ProSolve *prosolve, vis_Model *model,
             vsy_TextFun *textfun, Vchar *htm)
{
   Vint i;
   vfx_ProRender *prorender;
   vis_State *state;
   vis_Connect *connect;
   vis_GridFun *gridfun;
   vis_RProp *rprop;
   Vchar gifb[256],gifc[256], gifd[256], gife[256];
   Vint glen;
   Vint analysistype;
   Vint numnp,numel;

   rprop = vis_RPropBegin();
   state = vis_StateBegin();

   vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
   vis_ConnectNumber (connect,SYS_NODE,&numnp);
   vis_ConnectNumber (connect,SYS_ELEM,&numel);

   vfx_LogFileTableInit (logfile,"Final Results Rendering");

                   /* determine analysis type */ 
   vfx_ProSolveGetInteger (prosolve,PROSOLVE_ANALYSIS,&analysistype);
   if(analysistype == SYS_ANALYSIS_STRUCTURAL) {
      vfx_LogFileTableEntry3 (logfile,LOGFILE_STRING,"Displacement Magnitude",
                                      LOGFILE_SPACE,"  ",
                                      LOGFILE_STRING,"Stress VonMises");
   } else {
      vfx_LogFileTableEntry3 (logfile,LOGFILE_STRING,"Temperature",
                                      LOGFILE_SPACE,"  ",
                                      LOGFILE_STRING,"Flux Magnitude");
   }
                   /* create file names for image */ 
   strcpy(gifc,htm);
   strcpy(gife,htm);
   glen = (Vint)strlen(htm);
   for(i = glen-1; i >= 0; i--) {
      if(htm[i] == '/'  ||  htm[i] == '\\') {
         break;
      } else if(i == 0) {
         i = -1;
         break;
      }
   }
   strcpy(gifb,&gifc[i+1]);
   strcat(gifb,"_dis.gif");
   strcpy(&gifc[i+1],gifb);
   strcpy(gifd,&gife[i+1]);
   strcat(gifd,"_str.gif");
   strcpy(&gife[i+1],gifd);
  
   prorender = vfx_ProRenderBegin();
   vfx_ProRenderSetDimension (prorender,240,180);
   vfx_ProRenderSetObject (prorender,VIS_MODEL,model);

                   /* get displacement or temperature result */ 
   if(analysistype == SYS_ANALYSIS_STRUCTURAL) {
      vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
      vis_RPropSetType (rprop,SYS_RES_D);
   } else {
      vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
      vis_RPropSetType (rprop,SYS_RES_TEMP);
   }
   vfx_ProSolveState (prosolve,rprop,state);
   if(vfx_ProSolveError(prosolve)) {
      vfx_LogFileTableTerm (logfile);
      vis_RPropEnd (rprop);
      vis_StateEnd (state);
      return;
   }
                   /* generate and draw */ 
   vfx_ProRenderSetObject (prorender,VIS_STATE,state);
   vfx_ProRenderState (prorender);
   vfx_ProRenderSetParamc (prorender,PRORENDER_OUTFILE,gifc);
   vfx_ProRenderSetLookAt (prorender,1.,1.,1.,0.,0.,0.,0.,0.,1.);
   vfx_ProRenderDraw (prorender);

                   /* get stress or flux */ 
   if(analysistype == SYS_ANALYSIS_STRUCTURAL) {
      vis_RPropDef (rprop,SYS_ELEM,SYS_NODE);
      vis_RPropSetType (rprop,SYS_RES_S);
   } else {
      vis_RPropDef (rprop,SYS_ELEM,SYS_NODE);
      vis_RPropSetType (rprop,SYS_RES_HEAT_FLUX);
   }
                   /* register GridFun object in State */ 
   gridfun = vis_GridFunBegin ();
   vis_ConnectGridFun (connect,gridfun);
   vis_StateSetObject (state,VIS_GRIDFUN,gridfun);

   vfx_ProSolveState (prosolve,rprop,state);
   if(vfx_ProSolveError(prosolve)) {
      vfx_LogFileTableTerm (logfile);
      vis_RPropEnd (rprop);
      vis_StateEnd (state);
      vis_GridFunEnd (gridfun);
      return;
   }
                   /* generate and draw */ 
   vfx_ProRenderSetObject (prorender,VIS_STATE,state);
   vfx_ProRenderState (prorender);
   vfx_ProRenderSetParamc (prorender,PRORENDER_OUTFILE,gife);
   vfx_ProRenderSetLookAt (prorender,1.,1.,1.,0.,0.,0.,0.,0.,1.);
   vfx_ProRenderDraw (prorender);
   vfx_ProRenderEnd (prorender);

                   /* place images in log file */ 
   vfx_LogFileTableEntry3 (logfile,LOGFILE_IMAGE,gifb,
                                   LOGFILE_SPACE,"  ",
                                   LOGFILE_IMAGE,gifd);
   vfx_LogFileTableTerm (logfile);
   vis_RPropEnd (rprop);
   vis_StateEnd (state);
   vis_GridFunEnd (gridfun);
}

/*----------------------------------------------------------------------
                      setup SProp for singularity detection
----------------------------------------------------------------------*/
static void
checksingular(vfx_ProSolve *prosolve)
{
   vis_Model *model;
   vis_SProp *sprop_old, *sprop_new;
   vsy_List  *spl;
   Vint flag, analysis, rcaseid, mcaseid;

                   /* extract original SProp */ 
   vfx_ProSolveGetObject (prosolve,VIS_MODEL,(Vobject**)&model);
   vis_ModelGetList (model,VIS_SPROP,&spl);
   vsy_ListRef (spl,1,(Vobject**)&sprop_old);

                   /* create and populate new SProp */
   sprop_new = vis_SPropBegin ();
   vis_SPropDef (sprop_new,SYS_SOL_SINGULAR);
   vsy_ListInsert (spl,1,(Vobject*)sprop_new);
   vis_SPropValueFlag (sprop_old,SPROP_ANALYSIS,&flag);
   vis_SPropSetValuei (sprop_new,SPROP_EIGEN_NUM,12);
                   /* set analysis type and save singular modes */ 
   if(flag) {
      vis_SPropValueInteger (sprop_old,SPROP_ANALYSIS,&analysis);
      vis_SPropSetValuei (sprop_new,SPROP_ANALYSIS,analysis);
      vis_SPropSetValuei (sprop_new,SPROP_RESFILE_NUM,1);
      if(analysis == SYS_ANALYSIS_STRUCTURAL) {
         vis_SPropSetValuei (sprop_new,SPROP_RESFILE,SYS_RES_D);
      } else {
         printf("Singular analysis only valid for SYS_ANALYSIS_STRUCTURAL\n");
         return;
      }
   }
                   /* need restraint case */
   vis_SPropValueFlag (sprop_old,SPROP_RCASE,&flag);
   if(flag) {
      vis_SPropValueInteger (sprop_old,SPROP_RCASE,&rcaseid);
      vis_SPropSetValuei (sprop_new,SPROP_RCASE,rcaseid);
   }
                   /* need MPC case */
   vis_SPropValueFlag (sprop_old,SPROP_MCASE,&flag);
   if(flag) {
      vis_SPropValueInteger (sprop_old,SPROP_MCASE,&mcaseid);
      vis_SPropSetValuei (sprop_new,SPROP_MCASE,mcaseid);
   }
                   /* set ProSolve to only run this SProp */
   vfx_ProSolveRange (prosolve,1,1);
                   /* delete unused SProp */
   vis_SPropEnd (sprop_old);
}

5.5. Example 5, Generate Software Rendered Images

This example illustrates using ProRender to generate software rendered image files. The example uses standard VdmTools procedures to load a Model object from a file input on the command line.

The Model object is registered with ProRender using vfx_ProRenderSetObject(). An internal graphics object is generated using vfx_ProRenderModel(). This graphics object is then rendered twice using different viewpoints. Each rendering requires a different output file name to be specified. The viewpoint is specified using vfx_ProRenderSetLookAt(). The actual rendering occurs using vfx_ProRenderDraw().

#include "base/base.h"
#include "vfx/vfx.h"
#include "vfx/prorender.h"

/*----------------------------------------------------------------------
                      Generate software rendered images
----------------------------------------------------------------------*/
int
main(int argc, char *argv[])
{
   vdm_LMan *lman;
   Vint ierr;
                   /* libraries */
   vdm_SDRCLib *sdrclib;
   vdm_RASLib  *raslib;
   vdm_ABALib  *abalib;
   vdm_ANSLib  *anslib;
   vdm_NASLib  *naslib;
   vdm_NatLib  *natlib;
   vdm_NASFil  *nasfil;
   vdm_STLFil  *stlfil;
   vdm_DataFun *datafun;
   char inputfile[256];
   Vint filetype;

   vis_Model *model;
   vfx_ProRender *prorender;

                   /* determine input file */
   if(argc < 2) {
      fprintf(stderr,"Usage: %s pathname\n",argv[0]);
      fprintf(stderr," chexa.dat is assumed\n");
      strcpy(inputfile,"chexa.dat");
   } else {
      strcpy(inputfile,argv[1]);
   }
                   /* create data function object */
   datafun = vdm_DataFunBegin ();

                   /* determine first file type from file extension */
   if(strstr(inputfile,".bdf") != NULL  ||
      strstr(inputfile,".dat") != NULL) {
      filetype = VDM_NASTRAN_BULKDATA;
      nasfil = vdm_NASFilBegin();
      vdm_NASFilDataFun (nasfil,datafun);
   } else if(strstr(inputfile,".op2") != NULL) {
      filetype = VDM_NASTRAN_OUTPUT2;
      naslib = vdm_NASLibBegin();
      vdm_NASLibDataFun (naslib,datafun);
   } else if(strstr(inputfile,".rst") != NULL) {
      filetype = VDM_ANSYS_RESULT;
      anslib = vdm_ANSLibBegin();
      vdm_ANSLibDataFun (anslib,datafun);
   } else if(strstr(inputfile,".fil") != NULL) {
      filetype = VDM_ABAQUS_FIL;
      abalib = vdm_ABALibBegin();
      vdm_ABALibDataFun (abalib,datafun);
   } else if(strstr(inputfile,".odb") != NULL) {
      filetype = VDM_ABAQUS_ODB;
      abalib = vdm_ABALibBegin();
      vdm_ABALibDataFun (abalib,datafun);
   } else if(strstr(inputfile,".neu") != NULL) {
      filetype = VDM_MECHANICA_STUDY;
      raslib = vdm_RASLibBegin();
      vdm_RASLibDataFun (raslib,datafun);
   } else if(strstr(inputfile,".unv") != NULL  ||
             strstr(inputfile,".bun") != NULL) {
      filetype = VDM_SDRC_UNIVERSAL;
      sdrclib = vdm_SDRCLibBegin();
      vdm_SDRCLibDataFun (sdrclib,datafun);
   } else if(strstr(inputfile,".stl") != NULL) {
      filetype = VDM_STL;
      stlfil = vdm_STLFilBegin();
      vdm_STLFilDataFun (stlfil,datafun);
   } else if(strstr(inputfile,".STL") != NULL) {
      filetype = VDM_STLBIN;
      stlfil = vdm_STLFilBegin();
      vdm_STLFilDataFun (stlfil,datafun);
   } else if(strstr(inputfile,".vdm") != NULL) {
      filetype = VDM_NATIVE;
      natlib = vdm_NatLibBegin();
      vdm_NatLibDataFun (natlib,datafun);
   } else {
      fprintf(stderr,"Error: Bad input file %s\n",inputfile);
      exit(0);
   }
                   /* open library device */
   vdm_DataFunOpen (datafun,0,inputfile,filetype);

                   /* check for error */
   ierr = vdm_DataFunError(datafun);
   if(ierr) {
      fprintf(stderr,"Error: opening file %s\n",inputfile);
      exit(0);
   }
                   /* instance model object for finite element model */
   lman = vdm_LManBegin ();
   vdm_LManSetObject (lman,VDM_DATAFUN,datafun);
   model = vis_ModelBegin ();
   vdm_LManLoadModel (lman,model);
   if(vdm_LManError(lman)) {
      fprintf(stderr,"Error: Unable to load model information\n");
      exit(0);
   }

   prorender = vfx_ProRenderBegin();
   vfx_ProRenderSetObject (prorender,VIS_MODEL,model);
                   /* set edge type before generating internal display lists */ 
   vfx_ProRenderSetParami (prorender,PRORENDER_EDGETYPE,CONNECT_UNIQUE);
   vfx_ProRenderModel (prorender);
   vfx_ProRenderSetParamc (prorender,PRORENDER_OUTFILE,"exam5_1.gif");
   vfx_ProRenderSetLookAt (prorender,1.,1.,1.,0.,0.,0.,0.,0.,1.);
   vfx_ProRenderDraw (prorender);

   vfx_ProRenderSetParamc (prorender,PRORENDER_OUTFILE,"exam5_2.gif");
   vfx_ProRenderSetLookAt (prorender,0.,0.,1.,0.,0.,0.,0.,1.,0.);
   vfx_ProRenderDraw (prorender);

                   /* destroy objects */ 
   vfx_ProRenderEnd (prorender);

                   /* close library device and delete interface */
   vdm_DataFunClose (datafun);
   vdm_DataFunEnd (datafun);
   if(filetype == VDM_NASTRAN_BULKDATA) {
      vdm_NASFilEnd (nasfil);
   } else if(filetype == VDM_NASTRAN_OUTPUT2) {
      vdm_NASLibEnd (naslib);
   } else if(filetype == VDM_ANSYS_RESULT) {
      vdm_ANSLibEnd (anslib);
   } else if(filetype == VDM_ABAQUS_FIL) {
      vdm_ABALibEnd (abalib);
   } else if(filetype == VDM_MECHANICA_STUDY) {
      vdm_RASLibEnd (raslib);
   } else if(filetype == VDM_SDRC_UNIVERSAL) {
      vdm_SDRCLibEnd (sdrclib);
   } else if(filetype == VDM_STL ||
             filetype == VDM_STLBIN) {
      vdm_STLFilEnd (stlfil);
   } else if(filetype == VDM_NATIVE) {
      vdm_NatLibEnd (natlib);
   }
   vdm_LManEnd (lman);
   vis_ModelDelete (model);
   vis_ModelEnd (model);
   return 0;
}

5.6. Example 6, Flexible Body Generation

This example illustrates using ProSolve to generate a flexible body using a superelement solution procedure. The superelement solution generates a set of modes and forms a reduced stiffness and mass matrix using the modes to project the full finite element degree of freedom system to the reduced system. A finite element model is created in the example function CreateModel. Two types of rigid elements are used to connect two attachment points to the adjacent finite element model, 1) kinematic couplings (NASTRAN RBE2) and 2) distributing couplings (NASTRAN RBE3). The computed modes consist of 12 attachment point modes, 10 vibration modes and 2 load modes. The solution procedure is defined in an SProp object. Various SProp properties are set to define the solution process.

The attachment point modes are defined as applied restraints in an RCase object. The RCase object is installed in a HashTable with an identifier of 1. The SPROP_RCASE solution property is set to 1 to point to this RCase. The vibration modes are defined using the SPROP_EIGEN_NUM solution property. The eigenmodes nearest to zero are requested.

The 2 load modes are defined using two State objects which are installed in a HashTable with identifiers 1 and 2 respectively. The SPROP_LCASE_NUM and SPROP_LCASE properties point to these State objects. SPROP_LCASE_MODE is enabled to indicate that the load cases are to be used as modes, and SPROP_LCASE_UNIT is enabled to indicated that each State is to be treated as an independent mode rather than to have them all added into a single load mode. Specific results are requested using the SPROP_RESFILE_NUM and SPROP_RESFILE solution properties, In this case displacements and stresses are requested. The final Model object encapsulates all aspects of the finite element model definition as well as the solution procedure.

The ProSolve object is instanced and parameters are set to suppress printing to the output file and to enable saving the model and final results to a VDM_NATIVE (.vdm) file. The Model object is registered with ProSolve using vfx_ProSolveSetObject(). The function vfx_ProSolveInit() initializes the solution process. The function vfx_ProSolveExec() performs the actual finite element solution. The displacement and stress modes must be queried during the solution process using a monitor function monitorfun. The final reduced stiffness and mass may be queried after vfx_ProSolveExec() returns using vfx_ProSolveSysMatrix(). The pointer to the monitor function is installed using vfx_ProSolveSetFunction(). This function provides for a user defined object to be propagated to the monitor function. In this example, a pointer to a user defined structure, UserObj, is used to provide a mechanism to retrieve the computed displacement and stress modes as they are computed during the execution of vfx_ProSolveExec().

The monitor function is called regularly during the execution of the vfx_ProSolveExec() function. Results are only available during the output phase. Query for the current phase using vfx_ProSolveGetInteger() in the monitor function. If in the output phase, then query for the current case identifier and step. For single step procedures such as superelement generation, the step is interpreted as a mode number. In the output phase the monitor function will be called once for each mode. The displacement and stress results are obtained using vfx_ProSolveState().

When the function vfx_ProSolveExec() returns, the resulting reduced (modal) stiffness and mass matrices may be queried using vfx_ProSolveSysMatrix(). This function returns a pointer to the SysMatrix object, owned by ProSolve, which contains the desired reduced matrix. In this example we compute the eigenvalues and eigenvectors of the reduced system using vfs_SysMatrixEigenAll(). There should be six rigid body modes recovered by this procedure. Note that since the SysMatrix objects are owned by ProSolve they should not be explicitly deleted using vfs_SysMatrixEnd() by the application.

#include <stdlib.h>
#include <math.h>
#include "vfx/vfx.h"
#include "vis/vismesh.h"

static void
CreateModel(vis_Model *model);
static void
monitorfun(vfx_ProSolve *prosolve, Vobject *object);

                   /* define user object for monitor function */ 
typedef struct {
   vis_State *stated;
   vis_State *statent;
   vis_RProp *rpropd;
   vis_RProp *rpropent;
   Vint       numnp;
   Vint       numel;
} UserObj;

/*----------------------------------------------------------------------
                      Generate Flexible Body Using Component Mode Synthesis
----------------------------------------------------------------------*/
int
main()
{
   vfx_ProSolve  *prosolve;
   vis_Model     *model;
   vis_RProp     *rpropred;
   vsy_TextFun   *textfun;
   vsy_PlainText *plaintext;
   vis_Connect   *connect;
   vis_GridFun   *gridfun;
   UserObj        user;
   vfs_SysMatrix *redk, *redm;
   vfs_SysVector *rede, **redv;
   Vint           n;
   Vint           type, neq;
   Vdouble        value;

                   /* Create Model object */
   model = vis_ModelBegin ();
   CreateModel (model);

                   /* Create ProSolve object */
   prosolve = vfx_ProSolveBegin ();
   vfx_ProSolveSetParami (prosolve,PROSOLVE_PRINTLEVEL,0);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,SYS_ON);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,SYS_ON);
   vfx_ProSolveSetParamc (prosolve,PROSOLVE_RESFILE,"exam6.vdm");
   vfx_ProSolveSetParami (prosolve,PROSOLVE_EIGEN_METHOD,SYS_EIGEN_LANCZOS);
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);
   vfx_ProSolveSetFunction (prosolve,PROSOLVE_FUN_MONITOR,monitorfun,
                            (Vobject*)&user);

                   /* Create log file output objects */
   textfun = vsy_TextFunBegin ();
   plaintext = vsy_PlainTextBegin ();
   vsy_PlainTextTextFun (plaintext,textfun);
   vsy_TextFunConnectFile (textfun,stdout);
   vfx_ProSolveSetObject (prosolve,VSY_TEXTFUN,textfun);

                   /* Create state retrieval objects */
   user.stated   = vis_StateBegin ();
   user.rpropd   = vis_RPropBegin ();
   user.statent  = vis_StateBegin ();
   user.rpropent = vis_RPropBegin ();
   vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
   vis_ModelGetObject (model,VIS_GRIDFUN,(Vobject**)&gridfun);
   vis_StateSetObject (user.statent,VIS_GRIDFUN,gridfun);
   vis_ConnectNumber (connect,SYS_NODE,&user.numnp);
   vis_ConnectNumber (connect,SYS_ELEM,&user.numel);

                   /* Generate flexible body */ 
   vfx_ProSolveInit (prosolve);
   vfx_ProSolveExec (prosolve);

                   /* Retrieve reduced stiffness and mass */ 
   rpropred = vis_RPropBegin();
   vis_RPropDef (rpropred,SYS_DOF,SYS_NONE);

   vis_RPropSetType (rpropred,SYS_RES_K_MAT);
   redk = vfs_SysMatrixBegin ();
   vfx_ProSolveSysMatrix (prosolve,rpropred,redk);

   vis_RPropSetType (rpropred,SYS_RES_M_MAT);
   redm = vfs_SysMatrixBegin ();
   vfx_ProSolveSysMatrix (prosolve,rpropred,redm);
   vfs_SysMatrixInq (redk,&type,&neq);

                   /* Extract eigenvalues of reduced system */ 
   vfs_SysMatrixSetParami (redk,SYSMATRIX_EIGEN,SYSMATRIX_EIGEN_VIBE);

                   /* vector of eigenvalues */ 
   rede = vfs_SysVectorBegin ();
   vfs_SysVectorDef (rede,neq);
                   /* array of eigenvectors */ 
   redv = (vfs_SysVector**)malloc(neq*sizeof(vfs_SysVector*));
   for(n = 0; n < neq; n++) {
      redv[n] = vfs_SysVectorBegin ();
      vfs_SysVectorDef (redv[n],neq);
   }
                   /* compute all eigenvalues */ 
   vfs_SysMatrixEigenAll (redk,redm,rede,redv);
                   /* print eigenvalues */ 
   printf("Eigenvalues of reduced system\n");
   for(n = 1; n <= neq; n++) {
      vfs_SysVectorGather (rede,1,&n,&value);
      printf("eigenvalue[%d]= %e\n",n,value);
   }
                   /* Delete model objects */ 
   vis_ModelDelete (model);
   vis_ModelEnd (model);
                   /* Delete other objects */ 
   vfs_SysVectorEnd (rede);
   for(n = 0; n < neq; n++) {
      vfs_SysVectorEnd (redv[n]);
   }
   vfs_SysMatrixEnd (redk);
   vfs_SysMatrixEnd (redm);
   free(redv);
   vis_StateEnd (user.stated);
   vis_StateEnd (user.statent);
   vis_RPropEnd (user.rpropd);
   vis_RPropEnd (user.rpropent);
   vis_RPropEnd (rpropred);
   vsy_PlainTextEnd (plaintext);
   vsy_TextFunEnd (textfun);
   vfx_ProSolveEnd (prosolve);
   return 0;
}

/*----------------------------------------------------------------------
             Function to monitor ProSolve execution
----------------------------------------------------------------------*/
static void
monitorfun(vfx_ProSolve *prosolve, Vobject *object)
{
   Vint  phase,step,caseid;
   UserObj *user = (UserObj*)object;

                   /* Retrieve ProSolve phase status */ 
   vfx_ProSolveGetInteger (prosolve,PROSOLVE_PHASE,&phase);

                   /* Retrieve displacement and stress modes at output phase */ 
                   /* Note that the states are overwritten for each mode.
                      In general the user should modify this procedure to
                      retain the states at each mode by providing a separate
                      state for each mode, writing to a file, etc. */
   if(phase == SYS_PHASE_OUTPUT) {
                   /* Retrieve mode number */ 
      vfx_ProSolveGetInteger (prosolve,PROSOLVE_CASE,&caseid);
      vfx_ProSolveGetInteger (prosolve,PROSOLVE_STEP,&step);

                   /* Retrieve mode shape */ 
      vis_RPropDef (user->rpropd,SYS_NODE,SYS_NONE);
      vis_RPropSetType (user->rpropd,SYS_RES_D);
      vis_RPropSetIds (user->rpropd,caseid,step,0);
      vfx_ProSolveState (prosolve,user->rpropd,user->stated);

                   /* Retrieve mode element node stresses */ 
      vis_RPropDef (user->rpropent,SYS_ELEM,SYS_NODE);
      vis_RPropSetIds (user->rpropent,caseid,step,0);
      vis_RPropSetType (user->rpropent,SYS_RES_S);
      vfx_ProSolveState (prosolve,user->rpropent,user->statent);
   }
}

static void
CreateModel(vis_Model *model)
{
   vis_MapMesh   *mapmesh;
   vis_RCase     *rcase;
   vis_EProp     *eprop;
   vis_MProp     *mprop;
   vis_SProp     *sprop;
   vis_LCase     *lcase;
   vis_Connect   *connect;
   vis_GridFun   *gridfun;
   vsy_HashTable *mproph, *eproph, *rcaseh, *lcaseh;
   vsy_List      *spropl;
   Vdouble        pres[4],force[3];
   Vint           lcases[2];
   Vint           restypes[4];
   Vdouble        zero = 0.;

                   /* geometry parameters */
   Vdouble radius = 10.;
   Vdouble theta  = 90.;
   Vdouble width  = 125.;
   Vdouble thick  = 2.;

                   /* material properties */ 
   Vdouble young   = 210000.;
   Vdouble poisson = 0.3;
   Vdouble density = 7.5e-09;

                   /* set patch point connectivity */
   Vint pix[8] = {1,2,3,4,5,6,7,8};
   Vint nelemcurve =  8;
   Vint nelemwidth = 40;
   
                   /* number of component modes requested */ 
   Vint nummodes = 10;

                   /* other parameters */
   Vint i, m, n;
   Vint tags[6] = {
      SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TZ,
      SYS_DOF_RX, SYS_DOF_RY, SYS_DOF_RZ };
   Vdouble rbe2center[3];
   Vdouble rbe3center[3];
   Vint rbe2node;
   Vint rbe3node;
   Vdouble coords[8][3];
   Vdouble dtheta;
   Vint numnp, numel, nedge, aid, numshell;
   Vint *ix;
   Vint dep[2],indep[2];

   dtheta = theta*(3.14159265/180.)/3.;

                   /* generate 4 points at each end of cylinder patch */
   for(i = 0; i < 4; i++) {
      coords[i][0] = radius*(1.-cos(dtheta*i));
      coords[i][1] = radius*sin(dtheta*i);
      coords[i][2] = 0.;

      coords[i+4][0] = coords[i][0];
      coords[i+4][1] = coords[i][1];
      coords[i+4][2] = -width;
   }
                   /* create MapMesh object */
   mapmesh = vis_MapMeshBegin ();
   vis_MapMeshDef (mapmesh,8,1);

                   /* create GridFun object */ 
   gridfun = vis_GridFunBegin ();

                   /* create Connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,0,0);
                   /* set up GridFun */ 
   vis_ConnectGridFun (connect,gridfun);

                   /* define points */ 
   for(i = 0; i < 8; i++) {
      vis_MapMeshSetPoint (mapmesh,i+1,coords[i]);
   }
                   /* define patch connectivity */ 
   vis_MapMeshSetParami (mapmesh,VIS_MESH_MAXI,2);
   vis_MapMeshSetPatch (mapmesh,1,VIS_SHAPEQUAD,4,2,0,pix);
   vis_MapMeshSetPatchParami (mapmesh,1,VIS_MESH_SHAPE,VIS_SHAPEQUAD);
   vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMI,nelemcurve);
   vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMJ,nelemwidth);

                   /* set assoc of 100 at z=0 and 200 at z=-width */ 
                   /* use these to identify RBE connections */
   vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID,1,SYS_EDGE,1,100);
   vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID,1,SYS_EDGE,3,200);

                   /* generate shell mesh */
   vis_MapMeshGenerate (mapmesh,connect);
                   /* delete MapMesh objects */ 
   vis_MapMeshEnd (mapmesh);

                   /* extract number of nodes and elements in shell */ 
   vis_ConnectNumber (connect,SYS_NODE,&numnp);
   vis_ConnectNumber (connect,SYS_ELEM,&numel);
   printf("Number of shell nodes    = %5d\n",numnp);
   printf("Number of shell elements = %5d\n",numel);
   numshell = numel;

                   /* set useful shell element associations and property id */ 
   for(n = 1; n <= numel; n++) {
      vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,n,SYS_ELEM_SHELL);
      vis_ConnectSetElemAssoc (connect,VIS_FEATECH,n,SYS_TECH_ANS);
      vis_ConnectSetElemAssoc (connect,VIS_PROPID,n,1);
   }

                   /* add center nodes for attachment points */
   rbe2center[0] = coords[3][0];
   rbe2center[1] = 0.;
   rbe2center[2] = coords[3][2];
   numnp += 1;
   rbe2node = numnp;
   vis_ConnectSetCoordsdv (connect,rbe2node,rbe2center);
   rbe3center[0] = coords[7][0];
   rbe3center[1] = 0.;
   rbe3center[2] = coords[7][2];
   numnp += 1;
   rbe3node = numnp;
   vis_ConnectSetCoordsdv (connect,rbe3node,rbe3center);

                   /* count number of nodes on edge */
   nedge = 0;
   for(i = 1; i <= numnp; i++) {
      vis_ConnectNodeAssoc (connect,VIS_MISCID,1,&i,&aid);
      if(aid == 100) {
        nedge++;
      }
   }
   ix = (Vint*)malloc((nedge+1)*sizeof(Vint));

                   /* add rbe2 (kinematic coupling) as multiple elements */
                   /* rigid property id 2 */
   ix[0] = rbe2node;
   for(i = 1; i <= numnp; i++) {
      vis_ConnectNodeAssoc (connect,VIS_MISCID,1,&i,&aid);
      if(aid == 100) {
         ix[1] = i;
         numel += 1;
         vis_ConnectSetTopology (connect,numel,SYS_SHAPEPOINT,2,0,0);
         vis_ConnectSetElemNode (connect,numel,ix);
         vis_ConnectSetElemAssoc (connect,VIS_PROPID,numel,2);
         vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,numel,SYS_ELEM_RIGID);
         vis_ConnectSetElemAssoc (connect,VIS_FEASPEC,numel,SYS_RIGID_KINE);
      }
   }

                   /* add rbe3 (distributing coupling) as a single element */
                   /* rigid property id 3 */
   ix[0] = rbe3node;
   m = 0;
   for(i = 1; i <= numnp; i++) {
      vis_ConnectNodeAssoc (connect,VIS_MISCID,1,&i,&aid);
      if(aid == 200) {
         m++;
         ix[m] = i;
      }
   }
   numel += 1;
   vis_ConnectSetTopology (connect,numel,SYS_SHAPEPOINT,nedge+1,0,0);
   vis_ConnectSetElemNode (connect,numel,ix);
   vis_ConnectSetElemAssoc (connect,VIS_PROPID,numel,3);
   vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,numel,SYS_ELEM_RIGID);
   vis_ConnectSetElemAssoc (connect,VIS_FEASPEC,numel,SYS_RIGID_DIST);
   free(ix);

                   /* finite element nodes and elements complete */
   printf("Total number of nodes    = %5d\n",numnp);
   printf("Total number of elements = %5d\n",numel);

                   /* hash table of restraint cases */ 
   rcaseh = vsy_HashTableBegin();

                   /* Suppress freedoms at rbe nodes */
   rcase = vis_RCaseBegin();
   for(i = 0; i < 6; i++) {
      vis_RCaseSetSPCdv (rcase,rbe2node,tags[i],RCASE_APPLIED,&zero,0);
      vis_RCaseSetSPCdv (rcase,rbe3node,tags[i],RCASE_APPLIED,&zero,0);
   }
   vsy_HashTableInsert (rcaseh,1,rcase);

                   /* hash table of element properties */ 
   eproph = vsy_HashTableBegin();

                   /* shell element properties */ 
   eprop = vis_EPropBegin();
   vis_EPropDef (eprop,SYS_ELEM_SHELL);
   vis_EPropSetValued (eprop,EPROP_THICKNESS,thick);
   vis_EPropSetValuei (eprop,EPROP_MID,1);
   vsy_HashTableInsert (eproph,1,eprop);

                   /* rigid element properties */ 
                   /* rbe2 property */ 
   eprop = vis_EPropBegin();
   vis_EPropDef (eprop,SYS_ELEM_RIGID);
   vis_EPropSetValued (eprop,EPROP_PENALTY,1.e+10);
                   /* 63 is binary 111111, indicating 3 translations
                      and 3 rotations */ 
   vis_EPropSetValuei (eprop,EPROP_DOFFLAG_NUM,2);
   dep[0] = 0;
   dep[1] = 63;
   indep[0] = 63;
   indep[1] = 0;
   vis_EPropSetValueiv (eprop,EPROP_DOFFLAG_DEP,dep);
   vis_EPropSetValueiv (eprop,EPROP_DOFFLAG_IND,indep);
   vsy_HashTableInsert (eproph,2,eprop);

                   /* rbe3 property */ 
   eprop = vis_EPropBegin();
   vis_EPropDef (eprop,SYS_ELEM_RIGID);
   vis_EPropSetValued (eprop,EPROP_PENALTY,1.e+10);
   vis_EPropSetValuei (eprop,EPROP_NUMDIST,nedge);
   vsy_HashTableInsert (eproph,3,eprop);

                   /* hash table of material properties */ 
   mproph = vsy_HashTableBegin();

                   /* isotropic material properties */ 
   mprop = vis_MPropBegin();
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   vis_MPropSetValued (mprop,MPROP_DENSITY,density);
   vis_MPropSetValued (mprop,MPROP_E,young);
   vis_MPropSetValued (mprop,MPROP_NU,poisson);
   vsy_HashTableInsert (mproph,1,mprop);

                   /* load case 1: fz point load at node 50 */
   lcaseh = vsy_HashTableBegin ();
   lcase = vis_LCaseBegin ();
   force[0] = 0.;
   force[1] = 0.;
   force[2] = 1.e+6;
   vis_LCaseAddConcdv(lcase,150,LCASE_FORCE,force);
   vsy_HashTableInsert (lcaseh,1,lcase);

                   /* load case 2: pressure */
   lcase = vis_LCaseBegin ();
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,gridfun);
   pres[0] = 1.e+3;
   pres[1] = 1.e+3;
   pres[2] = 1.e+3;
   pres[3] = 1.e+3;
   for(n = 1; n <= numshell; n++) {
      vis_LCaseAddDistdv(lcase,SYS_FACE,n,1,LCASE_PRES,pres);
   }
   vsy_HashTableInsert (lcaseh,2,lcase);

                   /* list of solution properties */ 
   spropl = vsy_ListBegin();

                   /* solution properties for superelement generation */ 
   sprop = vis_SPropBegin();
   vis_SPropDef (sprop,SYS_SOL_SUPERELEMENT);
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
   vis_SPropSetValuei (sprop,SPROP_EIGEN_TYPE,SYS_EIGEN_NEAREST);
   vis_SPropSetValued (sprop,SPROP_EIGEN_SHIFT,0.);
   vis_SPropSetValuei (sprop,SPROP_EIGEN_NUM,nummodes);
   vis_SPropSetValuei (sprop,SPROP_RCASE,1);
   vis_SPropSetValuei (sprop,SPROP_MASSDIAG,SYS_ON);
                   /* two results, displacement and stress */ 
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,4);
   restypes[0] = SYS_RES_D;
   restypes[1] = SYS_RES_S;
   restypes[2] = SYS_RES_K_MAT;
   restypes[3] = SYS_RES_M_MAT;
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,restypes);
                   /* two load cases 1 and 2 */
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,2);
   lcases[0] = 1;
   lcases[1] = 2;
   vis_SPropSetValueiv (sprop,SPROP_LCASE,lcases);
   vis_SPropSetValuei  (sprop,SPROP_LCASE_UNIT,SYS_ON);
   vis_SPropSetValuei  (sprop,SPROP_LCASE_MODE,SYS_ON);
   vsy_ListInsert (spropl,1,(Vobject*)sprop);

                   /* load Model object */ 
   vis_ModelSetObject (model,VIS_CONNECT,connect);
   vis_ModelSetObject (model,VIS_GRIDFUN,gridfun);
   vis_ModelSetHashTable (model,VIS_EPROP,eproph);
   vis_ModelSetHashTable (model,VIS_MPROP,mproph);
   vis_ModelSetHashTable (model,VIS_RCASE,rcaseh);
   vis_ModelSetHashTable (model,VIS_LCASE,lcaseh);
   vis_ModelSetList (model,VIS_SPROP,spropl);
}

5.7. Example 7, Error Analysis and Mesh Refinement

This example illustrates computing error estimates of linear static analysis results and performing H-mesh refinement until analysis error norms are reduced to a user specified tolerance. The major modules used are the VisTools SurfMesh and TetMesh modules for performing the meshing operations, ElemChk for computing error estimates and new element sizes, and the VfxTools ProSolve module for performing the finite element analysis. The model is an L-shaped solid fixed at one end and loaded at the other. This results in a high stress (theoretically infinite) at the crease of the L-shape. The process involves generating an initial surface mesh and volume mesh of the solid and then entering a refinement loop in which an analysis is performed and error estimates are computed. If the error is below the user specified tolerance, the loop exits. Otherwise, new element sizes are predicted, the volume mesh is refined and the loop is traversed again.

The command line arguments to the example are the tolerance and the element size. If no arguments are specified, the tolerance is set to .05 (5 percent) and the initial element size is the width of the vertical leg of the L-shape (2.).

#include <stdlib.h>
#include "vfx/vfx.h"

                    /* geometry */
/*
   (fixed face)
      -----
      |   |
      |   |
    D |   |                 depth = E
      |   |
      |   |
      ------------------------------
      |   |                        |
    C |   |                        | (unit traction face in vertical direction)
      |   |                        |
      ------------------------------
        A             B
*/

#define A   2.0
#define B   8.0
#define C   1.0
#define D   6.0
#define E   1.0

#define VIEW_FLAG(flags,ind) (((flags) >> ((ind)-1)) & 1)

static void
generate_surfmesh(void);
static void
object_setup(Vdouble tol);
static void
load_bc(void);
static void
generate_tetmesh(void);
static void
refine_tetmesh(Vdouble tol, Vint *iter);

static vis_SurfMesh  *surfmesh    = NULL;
static vis_TetMesh   *tetmesh     = NULL;
static vis_Connect   *surfconnect = NULL;
static vis_Connect   *tetconnect  = NULL;
static vis_Model     *model       = NULL;
static vfx_ProSolve  *prosolve    = NULL;
static vsy_HashTable *lch         = NULL;
static vsy_HashTable *rch         = NULL;
static vsy_HashTable *eph         = NULL;
static vsy_HashTable *mph         = NULL;
static vsy_List      *spl         = NULL;
static vis_EProp     *eprop       = NULL;
static vis_MProp     *mprop       = NULL;
static vis_SProp     *sprop       = NULL;
static vis_LCase     *lcase       = NULL;
static vis_RCase     *rcase       = NULL;
static vis_Group     *ngroup      = NULL;
static vis_Group     *egroup      = NULL;
static vis_Group     *fgroup      = NULL;
static vsy_TextFun   *textfun     = NULL;
static vsy_PlainText *plaintext   = NULL;
static vis_GridFun   *gridfun     = NULL;
static vfe_MatlFun   *matlfun     = NULL;
static vis_State     *errorstate  = NULL;
static vis_State     *enerstate   = NULL;
static vis_State     *sesize      = NULL;
static vis_State     *tesize      = NULL;
static vis_RProp     *errorrprop  = NULL;
static vis_RProp     *enerrprop   = NULL;
static vis_IdTran    *idtran      = NULL;
static vis_ElemChk   *elemchk     = NULL;

/*----------------------------------------------------------------------
                Mesh refinement and error analysis
----------------------------------------------------------------------*/

static Vint structural = 1;
static Vint linearelem = 0;

int
main(int argc, char *argv[])
{
   Vdouble tol, edgelength, minedgelength;
   Vint    iter = 1;
   Vchar   name[80];

                    /* read error tolerance */
   if(argc < 2) {
      printf("Usage: %s tolerance\n",argv[0]);
      tol = .05;
      printf("assume tolerance= %f\n",tol);
   } else {
      sscanf(argv[1],"%le",&tol);
   }
                    /* read initial element size */
   if(argc >= 3) {
      sscanf(argv[2],"%le",&edgelength);
   } else {
      edgelength = A;
      printf("assume edgelength= %f\n",edgelength);
   }
                   /* object setup */
   object_setup (tol);

                   /* set target lengths */
   minedgelength = 0.1*edgelength;
   vis_SurfMeshSetParamd (surfmesh,VIS_MESH_EDGELENGTH,edgelength);
   vis_SurfMeshSetParamd (surfmesh,VIS_MESH_MINEDGELENGTH,minedgelength);
   vis_TetMeshSetParamd (tetmesh,VIS_MESH_EDGELENGTH,edgelength);

                   /* generate initial surface mesh */
   generate_surfmesh ();

                   /* loop until tolerance is achieved */
   do {
                   /* generate tet mesh */
      if(iter == 1) {
         generate_tetmesh ();
      }
                   /* add loads and boundary conditions and set pids */
      load_bc ();

                   /* solve */
      sprintf(name,"exam7_%d.vdm",iter);
      vfx_ProSolveInit (prosolve);
      vfx_ProSolveSetParamc (prosolve,PROSOLVE_RESFILE,name);
      vfx_ProSolveExec (prosolve);

                   /* refine tetmesh using error estimation */
                   /* iter = 0, means convergence obtained */
      refine_tetmesh (tol,&iter);
   } while(iter);

                   /* clean up */
   vis_ModelDelete (model);
   vis_ModelEnd (model);
   vis_SurfMeshEnd (surfmesh);
   vis_TetMeshEnd (tetmesh);
   vis_ConnectEnd (surfconnect);
   vfx_ProSolveEnd (prosolve);
   vis_GroupEnd (ngroup);
   vis_GroupEnd (egroup);
   vis_GroupEnd (fgroup);
   vsy_TextFunEnd (textfun);
   vsy_PlainTextEnd (plaintext);
   vis_GridFunEnd (gridfun);
   vfe_MatlFunEnd (matlfun);
   vis_StateEnd (errorstate);
   vis_StateEnd (enerstate);
   vis_StateEnd (sesize);
   vis_StateEnd (tesize);
   vis_RPropEnd (errorrprop);
   vis_RPropEnd (enerrprop);
   vis_IdTranEnd (idtran);
   vis_ElemChkEnd (elemchk);
   return 0;
}

/*----------------------------------------------------------------------
                      generate geometry
----------------------------------------------------------------------*/
static void
generate_surfmesh(void)
{
   Vdouble xgeom[16][3];
   Vint eflag[28][3]  =
      {  {1,0,1},{0,0,0},{0,0,0},{1,1,0},{1,0,0},{1,1,0},
         {1,0,1},{0,0,0},{1,0,0},{0,1,1},{0,0,1},{0,1,1},
         {1,1,0},{0,1,1},{1,1,0},{0,1,1},
         {1,0,0},{0,1,1},{1,1,0},{0,1,0},
         {1,1,0},{0,1,0},{1,0,0},{0,1,1},
         {1,0,1},{1,1,0},{1,0,1},{1,1,0}  };
   Vint ixgeom[28][3] =
      {  { 1, 4, 2},{ 4, 5, 2},{ 2, 5, 3},{ 5, 6, 3},{ 4, 7, 5},{ 7, 8, 5},
         { 9,10,12},{12,10,13},{10,11,13},{13,11,14},{12,13,15},{15,13,16},
         {11, 3, 6},{11, 6,14},{13, 5, 8},{13, 8,16},
         { 9,12, 4},{ 9, 4, 1},{12,15, 7},{12, 7, 4},
         {10, 9, 1},{10, 1, 2},{11,10, 2},{11, 2, 3},
         {15,16, 7},{16, 8, 7},{13,14, 5},{14, 6, 5}  };
   Vdouble normal[28][3][3] =
      {  {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  {-1.,0.,0.},{-1.,0.,0.},{-1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  { 1.,0.,0.},{ 1.,0.,0.},{ 1.,0.,0.}  },
         {  {0., 1.,0.},{0., 1.,0.},{0., 1.,0.}  },
         {  {0., 1.,0.},{0., 1.,0.},{0., 1.,0.}  },
         {  {0., 1.,0.},{0., 1.,0.},{0., 1.,0.}  },
         {  {0., 1.,0.},{0., 1.,0.},{0., 1.,0.}  },
         {  {0.,-1.,0.},{0.,-1.,0.},{0.,-1.,0.}  },
         {  {0.,-1.,0.},{0.,-1.,0.},{0.,-1.,0.}  },
         {  {0.,-1.,0.},{0.,-1.,0.},{0.,-1.,0.}  },
         {  {0.,-1.,0.},{0.,-1.,0.},{0.,-1.,0.}  },
         {  {0.,0.,-1.},{0.,0.,-1.},{0.,0.,-1.}  },
         {  {0.,0.,-1.},{0.,0.,-1.},{0.,0.,-1.}  },
         {  {0.,0.,-1.},{0.,0.,-1.},{0.,0.,-1.}  },
         {  {0.,0.,-1.},{0.,0.,-1.},{0.,0.,-1.}  },
         {  {0.,0., 1.},{0.,0., 1.},{0.,0., 1.}  },
         {  {0.,0., 1.},{0.,0., 1.},{0.,0., 1.}  },
         {  {0.,0., 1.},{0.,0., 1.},{0.,0., 1.}  },
         {  {0.,0., 1.},{0.,0., 1.},{0.,0., 1.}  }  };
   Vint n;
   Vint numnp,numel;

                   /* back face points */
   xgeom[ 0][0] = 0.0; xgeom[ 0][1] = 0.0; xgeom[ 0][2] = 0.0;
   xgeom[ 1][0] = 0.0; xgeom[ 1][1] = A;   xgeom[ 1][2] = 0.0;
   xgeom[ 2][0] = 0.0; xgeom[ 2][1] = A+B; xgeom[ 2][2] = 0.0;
   xgeom[ 3][0] = 0.0; xgeom[ 3][1] = 0.0; xgeom[ 3][2] = C;  
   xgeom[ 4][0] = 0.0; xgeom[ 4][1] = A;   xgeom[ 4][2] = C;  
   xgeom[ 5][0] = 0.0; xgeom[ 5][1] = A+B; xgeom[ 5][2] = C;    
   xgeom[ 6][0] = 0.0; xgeom[ 6][1] = 0.0; xgeom[ 6][2] = C+D;
   xgeom[ 7][0] = 0.0; xgeom[ 7][1] = A;   xgeom[ 7][2] = C+D;

                   /* front face points */
   xgeom[ 8][0] = E;   xgeom[ 8][1] = 0.0; xgeom[ 8][2] = 0.0;
   xgeom[ 9][0] = E;   xgeom[ 9][1] = A;   xgeom[ 9][2] = 0.0;
   xgeom[10][0] = E;   xgeom[10][1] = A+B; xgeom[10][2] = 0.0;
   xgeom[11][0] = E;   xgeom[11][1] = 0.0; xgeom[11][2] = C;  
   xgeom[12][0] = E;   xgeom[12][1] = A;   xgeom[12][2] = C;  
   xgeom[13][0] = E;   xgeom[13][1] = A+B; xgeom[13][2] = C;    
   xgeom[14][0] = E;   xgeom[14][1] = 0.0; xgeom[14][2] = C+D;
   xgeom[15][0] = E;   xgeom[15][1] = A;   xgeom[15][2] = C+D;

                   /* set mesh parameters */
   vis_SurfMeshDef (surfmesh,16,28);
   if(linearelem) {
      vis_SurfMeshSetParami (surfmesh,VIS_MESH_MAXI,2);
   } else {
      vis_SurfMeshSetParami (surfmesh,VIS_MESH_MAXI,3);
   }
   vis_SurfMeshSetParami (surfmesh,VIS_MESH_SHAPE,VIS_SHAPETRI);

                   /* set surface geometry */
   for(n = 1; n <= 16; n++) {
      vis_SurfMeshSetPoint (surfmesh,n,xgeom[n-1],1);
   }
   for(n = 1; n <= 28; n++) {
      vis_SurfMeshSetTri (surfmesh,n,ixgeom[n-1],eflag[n-1]);
      vis_SurfMeshSetTriNorm (surfmesh,n,normal[n-1]);
   }

                   /* set association of 1 for fixed boundary */
   vis_SurfMeshSetTriAssoc (surfmesh,VIS_GEOFACE,25,SYS_FACE,0,1);
   vis_SurfMeshSetTriAssoc (surfmesh,VIS_GEOFACE,26,SYS_FACE,0,1);

                   /* set association of 2 for load boundary */
   vis_SurfMeshSetTriAssoc (surfmesh,VIS_GEOFACE,13,SYS_FACE,0,2);
   vis_SurfMeshSetTriAssoc (surfmesh,VIS_GEOFACE,14,SYS_FACE,0,2);

                /* generate surface mesh */
   vis_ConnectDef (surfconnect,0,0);
 
   vis_SurfMeshGenerate (surfmesh,surfconnect);
   vis_ConnectNumber (surfconnect,SYS_NODE,&numnp);
   vis_ConnectNumber (surfconnect,SYS_ELEM,&numel);
   printf("Original surface mesh: %d nodes, %d elements\n",numnp,numel);
}

/*----------------------------------------------------------------------
                      setup objects
----------------------------------------------------------------------*/
static void
object_setup(Vdouble tol)
{
   Vint numres, result[4];

                   /* instance all objects */
   textfun     = vsy_TextFunBegin ();
   plaintext   = vsy_PlainTextBegin ();
   surfmesh    = vis_SurfMeshBegin ();
   tetmesh     = vis_TetMeshBegin ();
   surfconnect = vis_ConnectBegin ();
   tetconnect  = vis_ConnectBegin ();
   prosolve    = vfx_ProSolveBegin ();
   model       = vis_ModelBegin ();
   lch         = vsy_HashTableBegin ();
   rch         = vsy_HashTableBegin ();
   eph         = vsy_HashTableBegin ();
   mph         = vsy_HashTableBegin ();
   spl         = vsy_ListBegin ();
   eprop       = vis_EPropBegin ();
   mprop       = vis_MPropBegin ();
   sprop       = vis_SPropBegin ();
   rcase       = vis_RCaseBegin ();
   lcase       = vis_LCaseBegin ();
   ngroup      = vis_GroupBegin ();
   egroup      = vis_GroupBegin ();
   fgroup      = vis_GroupBegin ();
   gridfun     = vis_GridFunBegin ();
   matlfun     = vfe_MatlFunBegin ();
   errorstate  = vis_StateBegin ();
   enerstate   = vis_StateBegin ();
   sesize      = vis_StateBegin ();
   tesize      = vis_StateBegin ();
   errorrprop  = vis_RPropBegin ();
   enerrprop   = vis_RPropBegin ();
   idtran      = vis_IdTranBegin ();
   elemchk     = vis_ElemChkBegin ();

                   /* setup textfun */
   vsy_PlainTextTextFun (plaintext,textfun);
   vsy_TextFunConnectFile (textfun,stdout);

                   /* setup gridfun */
   vis_ConnectGridFun (tetconnect,gridfun);

                   /* setup elemchk */
   vis_ElemChkSetObject (elemchk,VIS_GRIDFUN,gridfun);
   vis_ElemChkSetParamd (elemchk,ELEMCHK_ADAPTTOL,tol);

                   /* setup lcase */
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,gridfun);

                   /* setup prosolve */
   vfx_ProSolveSetParami (prosolve,PROSOLVE_PRINTLEVEL,0);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,SYS_ON);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,SYS_ON);

   vfx_ProSolveSetParami (prosolve,PROSOLVE_SOLVER,PROSOLVE_SYMM_ITER);
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);
   vfx_ProSolveSetObject (prosolve,VSY_TEXTFUN,textfun);
   vfx_ProSolveRange (prosolve,1,1);

                   /* setup sprop */
   vis_SPropDef (sprop,SYS_SOL_STATIC);
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
   vis_SPropSetValuei (sprop,SPROP_LCASE,1);
   vis_SPropSetValued (sprop,SPROP_LCASE_FACTOR,1.0);
   vis_SPropSetValuei (sprop,SPROP_RCASE,1);
   vis_SPropSetValued (sprop,SPROP_RCASE_FACTOR,1.0);
   if(structural) {
      vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
   } else {
      vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_THERMAL);
   }
                   /* these results are necessary for error estimation */ 
   numres = 4;
   if(structural) {
      result[0] = SYS_RES_D;
      result[1] = SYS_RES_S;
      result[2] = SYS_RES_SE_ERROR;
      result[3] = SYS_RES_SE;
   } else {
      result[0] = SYS_RES_TEMP;
      result[1] = SYS_RES_HEAT_FLUX;
      result[2] = SYS_RES_TE_ERROR;
      result[3] = SYS_RES_TE;
   }
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,numres);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,result);

                   /* setup eprop */
   vis_EPropDef (eprop,VIS_ELEM_SOLID);
   vis_EPropSetValuei (eprop,EPROP_MID,1);

                   /* setup mprop */
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   if(structural) {
      vis_MPropSetValued (mprop,MPROP_E,1.);
      vis_MPropSetValued (mprop,MPROP_NU,.0);
   } else {
      vis_MPropSetValued (mprop,MPROP_K,1.);
   }

                   /* setup rprop */
   vis_RPropDef (errorrprop,SYS_ELEM,SYS_NONE);
   vis_RPropSetIds (errorrprop,1,1,0);
   vis_RPropDef (enerrprop,SYS_ELEM,SYS_NONE);
   vis_RPropSetIds (enerrprop,1,1,0);
   if(structural) {
      vis_RPropSetType (errorrprop,SYS_RES_SE_ERROR);
      vis_RPropSetType (enerrprop,SYS_RES_SE);
   } else {
      vis_RPropSetType (errorrprop,SYS_RES_TE_ERROR);
      vis_RPropSetType (enerrprop,SYS_RES_TE);
   }

                   /* setup states */
   vis_StatePre (sesize,SYS_DOUBLE);
   vis_StatePre (errorstate,SYS_DOUBLE);
   vis_StatePre (enerstate,SYS_DOUBLE);
   vis_StatePre (tesize,SYS_DOUBLE);
   vis_StateSetObject (tesize,VIS_GRIDFUN,gridfun);

                   /* create model data structure */
   vis_ModelSetObject (model,VIS_CONNECT,tetconnect);
   vis_ModelSetHashTable (model,VIS_LCASE,lch);
   vis_ModelSetHashTable (model,VIS_RCASE,rch);
   vis_ModelSetHashTable (model,VIS_EPROP,eph);
   vis_ModelSetHashTable (model,VIS_MPROP,mph);
   vis_ModelSetList (model,VIS_SPROP,spl);
   vsy_HashTableInsert (lch,1,(Vobject*)lcase);
   vsy_HashTableInsert (rch,1,(Vobject*)rcase);
   vsy_HashTableInsert (eph,1,(Vobject*)eprop);
   vsy_HashTableInsert (mph,1,(Vobject*)mprop);
   vsy_ListInsert (spl,1,(Vobject*)sprop);
}

/*----------------------------------------------------------------------
                      loads and boundary conditions
----------------------------------------------------------------------*/
static void
load_bc(void)
{
   Vint n;
   Vint numnp, numel, no, flags;
   Vdouble heat[6] = {
      1.,1.,1.,1.,1.,1.};
   Vdouble ztraction[18] = {
      0.,0.,1.,0.,0.,1.,0.,0.,1.,0.,0.,1.,0.,0.,1.,0.,0.,1.};

                   /* initialize lcase and rcase */
   vis_RCaseClear (rcase);
   vis_LCaseClear (lcase);

   vis_ConnectNumber (tetconnect,SYS_NODE,&numnp);
   vis_ConnectNumber (tetconnect,SYS_ELEM,&numel);
   vis_ConnectKernel (tetconnect,0);

   vis_GroupDef (ngroup,numnp,SYS_NODE,SYS_NONE);
   vis_GroupDef (egroup,numel,SYS_ELEM,SYS_FACE);

   vis_ConnectSetGroupParami (tetconnect,CONNECT_ASSOCTYPE,VIS_GEOFACE);

                   /* set boundary conditions */
                   /* extract unique nodes on face 1 */
   vis_ConnectSetGroupParami (tetconnect,CONNECT_ASSOCID,1);
   vis_ConnectFaceGroup (tetconnect,CONNECT_ASSOC,NULL,egroup);
   vis_ConnectNodeGroup (tetconnect,CONNECT_UNIQUE,egroup,ngroup);
   vis_GroupInitIndex (ngroup);
   while(vis_GroupNextIndex (ngroup,&n,&flags), n) {
      if(structural) {
         vis_RCaseSetSPCdv (rcase,n,SYS_DOF_TX,RCASE_FIXED,NULL,0);
         vis_RCaseSetSPCdv (rcase,n,SYS_DOF_TY,RCASE_FIXED,NULL,0);
         vis_RCaseSetSPCdv (rcase,n,SYS_DOF_TZ,RCASE_FIXED,NULL,0);
      } else {
         vis_RCaseSetSPCdv (rcase,n,SYS_DOF_TEMP,RCASE_FIXED,NULL,0);
      }
   }
                   /* set loads */
                   /* extract element faces on face 2 */
   vis_ConnectSetGroupParami (tetconnect,CONNECT_ASSOCID,2);
   vis_GroupClear (egroup);
   vis_ConnectFaceGroup (tetconnect,CONNECT_ASSOC,NULL,egroup);
   vis_GroupInitIndex (egroup);
   while(vis_GroupNextIndex (egroup,&n,&flags), n) {
      for(no = 1; no <= 4; no++) {
         if(VIEW_FLAG(flags,no) == 0) continue;
         if(structural) {
            vis_LCaseAddDistdv (lcase,SYS_FACE,n,no,LCASE_TRAC,ztraction);
         } else {
            vis_LCaseAddDistdv (lcase,SYS_FACE,n,no,LCASE_HEATFLUX,heat);
         }
      }
   }
                   /* set pids */
   for(n = 1; n <= numel; n++) {
      vis_ConnectSetElemAssoc (tetconnect,VIS_PROPID,n,1);
      vis_ConnectSetElemAssoc (tetconnect,VIS_FEATYPE,n,SYS_ELEM_SOLID);
   }
}

/*----------------------------------------------------------------------
                      generate mesh
----------------------------------------------------------------------*/
static void
generate_tetmesh(void)
{
   Vint numel,numnp;

                /* generate tetrahedral mesh */
   vis_TetMeshConnect (tetmesh,surfconnect);

   vis_ConnectDef (tetconnect,0,0);
   vis_TetMeshGenerate (tetmesh,tetconnect);
   vis_ConnectNumber (tetconnect,SYS_NODE,&numnp);
   vis_ConnectNumber (tetconnect,SYS_ELEM,&numel);
   printf("Generate tet mesh: %d nodes, %d elements\n",numnp,numel);
   vis_ConnectKernel (tetconnect,0);
}

/*----------------------------------------------------------------------
                      refine mesh
----------------------------------------------------------------------*/
static void
refine_tetmesh(Vdouble tol, Vint *iter)
{
   Vdouble err,setot,seerr;
   Vint numnp, numel;

                   /* extract error and energy states */
   vfx_ProSolveState (prosolve,errorrprop,errorstate);
   vfx_ProSolveState (prosolve,enerrprop,enerstate);
                   /* sum total strain energy and error */ 
   vis_StateSumdv (errorstate,NULL,&seerr);
   vis_StateSumdv (enerstate,NULL,&setot);
   setot += seerr;
   if(setot != 0.) {
      err = sqrt(seerr/setot);
   } else {
      err = 0.;
   }
                   /* check convergence */ 
   printf("Error: %e\n",err);
   if(err < tol) {
      *iter = 0;
      return;
   }
                   /* get new mesh sizes */
   vis_ElemChkAdapt (elemchk,enerstate,errorstate,NULL,tesize);
                   /* refine tetmesh */ 
   *iter += 1;
   vis_ConnectDef (tetconnect,0,0);
   vis_TetMeshRefine (tetmesh,tesize,tetconnect);
   vis_ConnectNumber (tetconnect,SYS_NODE,&numnp);
   vis_ConnectNumber (tetconnect,SYS_ELEM,&numel);
   printf("Refined tetmesh: %d nodes, %d elements\n",numnp,numel);
   if(*iter == 10) *iter = 0;
}

5.8. Example 8, Thermal-Structural Analysis

This example illustrates using ProSolve to perform a thermal-structural analysis.

A Model object is initially generated where structural and thermal properties - including the coefficient of thermal expansion - are defined for all elements. Two SProp objects are created: the first performs a thermal analysis, while the second takes the temperatures computed in the first to perform a structural analysis. It is important to note that the SProp object for the structural analysis contains the flag SPROP_THERMALSTRAIN turned on. In addition, a reference temperature is set as a load case with the SProp parameter SPROP_LCASE_TREF.

Once the thermal analysis is completed the temperature State is retrieved and the temperature values are transfered to a concentrated load case of type LCASE_TEMP. Likewise, the reference temperatures are set in another State.

The changing of analysis type - from SYS_ANALYSIS_THERMAL to SYS_ANALYSIS_STRUCTURAL - requires ProSolve to be reinitialized with vfx_ProSolveInit().

At the end of the structural analysis displacements, stresses, and strains are retrieved and printed.

#include <stdio.h>
#include "vfx/vfx.h"
/*
                 4------5------6
                /.     /.     /|
  y            / .    / .    / |
  |           /  .   /  .   /  |
  --x        /   1../...2../...3
 /         10-----11-----12   /
z           |  .   |  .   |  /
            | .    | .    | /
            |.     |.     |/
            7------8------9
*/
#define MAX_ELEM  2
#define MAX_NODE 12

static Vint conn[MAX_ELEM][8] = {
   {1, 2, 5, 4, 7, 8,11,10},
   {2, 3, 6, 5, 8, 9,12,11} };

static Vdouble coords[MAX_NODE][3] = {
   {0.,0.,0.},  {1.,0.,0.},  {2.,0.,0.},
   {0.,1.,0.},  {1.,1.,0.},  {2.,1.,0.},
   {0.,0.,1.},  {1.,0.,1.},  {2.,0.,1.},
   {0.,1.,1.},  {1.,1.,1.},  {2.,1.,1.} };

/*----------------------------------------------------------------------
                      Thermal-Structural Analysis
----------------------------------------------------------------------*/
int
main()
{
   vis_Model *model;
   vis_GridFun *gridfun;
   vfx_ProSolve *prosolve;
   vsy_List *spl;
   vsy_HashTable *lch, *rch, *eph, *mph;
   vis_Connect *connect;
   vis_SProp *sprop;
   vis_LCase *lcase;
   vis_RCase *rcase;
   vis_EProp *eprop;
   vis_MProp *mprop;
   vis_RProp *rprop;
   vis_State *state;
   Vint i, n;
   Vint lcaset, rcaset, lcases, lcasetref, rcases;
   Vfloat tref      = 100.;
   Vfloat ptheat    = 5.e+3;
   Vfloat fixedtemp = 100.;
   Vfloat young     = 2.e+7;
   Vfloat poisson   = 0.3;
   Vfloat thermcond =  50.;
   Vfloat thermexp  = 1.e-6;
   Vfloat temperature;
   Vint   pid       =   7;
   Vint   mid       =   3;
   Vint   resultt    = SYS_RES_TEMP;
   Vint   results[3] = {SYS_RES_D, SYS_RES_S, SYS_RES_E};
   Vfloat data[100];

                   /* create and fill model object */
   model = vis_ModelBegin ();
   gridfun = vis_GridFunBegin ();

                   /* create connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,MAX_NODE,MAX_ELEM);

   for(i = 0; i < MAX_ELEM; i++) {
                   /* set element connectivity */
      vis_ConnectSetTopology (connect,i+1,SYS_SHAPEHEX,2,0,0);
      vis_ConnectSetElemNode (connect,i+1,conn[i]);

                   /* assign property id */
      vis_ConnectSetElemAssoc (connect,VIS_PROPID,i+1,pid);
   }
                   /* set nodal coordinates */
   for(i = 0; i < MAX_NODE; i++) {
      vis_ConnectSetCoordsdv (connect,i+1,coords[i]);
   }
   vis_ModelSetObject (model,VIS_CONNECT,connect);
   vis_ConnectGridFun (connect,gridfun);

                   /* define thermal and structural material properties */
   eph = vsy_HashTableBegin ();
   vis_ModelSetHashTable (model,VIS_EPROP,eph);
   eprop = vis_EPropBegin ();
   vsy_HashTableInsert (eph,pid,(Vobject*)eprop);
   vis_EPropDef (eprop,SYS_ELEM_SOLID);
   vis_EPropSetValuei (eprop,EPROP_MID,mid);

   mph = vsy_HashTableBegin ();
   vis_ModelSetHashTable (model,VIS_MPROP,mph);
   mprop = vis_MPropBegin ();
   vsy_HashTableInsert (mph,mid,(Vobject*)mprop);
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   vis_MPropSetValuef (mprop,MPROP_E,young);
   vis_MPropSetValuef (mprop,MPROP_NU,poisson);
   vis_MPropSetValuef (mprop,MPROP_K,thermcond);
   vis_MPropSetValuef (mprop,MPROP_A,thermexp);

                   /* create list of sprops, rcase, and lcase */
   spl = vsy_ListBegin ();
   vis_ModelSetList (model,VIS_SPROP,spl);
   lch = vsy_HashTableBegin ();
   vis_ModelSetHashTable (model,VIS_LCASE,lch);
   rch = vsy_HashTableBegin ();
   vis_ModelSetHashTable (model,VIS_RCASE,rch);

                   /* create SProp for thermal problem */
   lcaset = 10;
   rcaset = 20;
   sprop = vis_SPropBegin ();
   vsy_ListInsert (spl,1,(Vobject*)sprop);
   vis_SPropDef (sprop,SYS_SOL_STATIC);
   vis_SPropSetValuec (sprop,SPROP_TITLE,"Thermal Analysis");
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_THERMAL);
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
   vis_SPropSetValuei (sprop,SPROP_LCASE,lcaset);
   vis_SPropSetValuei (sprop,SPROP_RCASE,rcaset);
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,1);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,&resultt);

                   /* create thermal loads */
   lcase = vis_LCaseBegin ();
   vsy_HashTableInsert (lch,lcaset,(Vobject*)lcase);
   vis_LCaseAddConc (lcase, 7,LCASE_HEAT,&ptheat);
   vis_LCaseAddConc (lcase, 8,LCASE_HEAT,&ptheat);
   vis_LCaseAddConc (lcase, 9,LCASE_HEAT,&ptheat);
   vis_LCaseAddConc (lcase,10,LCASE_HEAT,&ptheat);
   vis_LCaseAddConc (lcase,11,LCASE_HEAT,&ptheat);
   vis_LCaseAddConc (lcase,12,LCASE_HEAT,&ptheat);

                   /* create thermal restraints at 100 degrees */
   rcase = vis_RCaseBegin ();
   vsy_HashTableInsert (rch,rcaset,(Vobject*)rcase);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TEMP,RCASE_APPLIED,&fixedtemp,0);
   vis_RCaseSetSPC (rcase,3,SYS_DOF_TEMP,RCASE_APPLIED,&fixedtemp,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TEMP,RCASE_APPLIED,&fixedtemp,0);
   vis_RCaseSetSPC (rcase,6,SYS_DOF_TEMP,RCASE_APPLIED,&fixedtemp,0);

                   /* create SProp for structural problem */
   lcases    = 11;
   lcasetref = 12;
   rcases    = 21;
   sprop = vis_SPropBegin ();
   vsy_ListInsert (spl,2,(Vobject*)sprop);
   vis_SPropDef (sprop,SYS_SOL_STATIC);
   vis_SPropSetValuec (sprop,SPROP_TITLE,"Structural Analysis");
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
   vis_SPropSetValuei (sprop,SPROP_LCASE,lcases);
   vis_SPropSetValuei (sprop,SPROP_RCASE,rcases);
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,3);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,results);
   vis_SPropSetValuei (sprop,SPROP_THERMALSTRAIN,SYS_ON);

                   /*
                    * Assign a LCase for reference temperature
                    * (reference temperature defaults to zero if undefined)
                    */
   vis_SPropSetValuei (sprop,SPROP_LCASE_TREF,lcasetref);

                   /* create structural restraints */
   rcase = vis_RCaseBegin ();
   vsy_HashTableInsert (rch,rcases,(Vobject*)rcase);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TY,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TZ,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TZ,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,6,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,7,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   printf("Finished building model\n");

                   /* instance and init ProSolve */
   prosolve = vfx_ProSolveBegin ();
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVEMODEL,SYS_OFF);
   vfx_ProSolveSetParami (prosolve,PROSOLVE_SAVERES,SYS_OFF);
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,(Vobject*)model);
   vfx_ProSolveInit (prosolve);

                   /* run ProSolve for thermal problem */
   vfx_ProSolveRange (prosolve,1,1);
   vfx_ProSolveExec (prosolve);
   printf("Finished running thermal analysis\n");

                   /* extract temperatures for thermal load */
   lcase = vis_LCaseBegin ();
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,(Vobject*)gridfun);
   vsy_HashTableInsert (lch,lcases,(Vobject*)lcase);
   rprop = vis_RPropBegin ();
   state = vis_StateBegin ();
   vis_StateSetObject (state,VIS_GRIDFUN,(Vobject*)gridfun);
   vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
   vis_RPropSetType (rprop,SYS_RES_TEMP);
   vis_StateDef (state,MAX_NODE,SYS_NODE,SYS_NONE,SYS_SCALAR);
   vfx_ProSolveState (prosolve,rprop,state);
   printf("Temperatures:\n");
   for(n = 1; n <= MAX_NODE; ++n) {
      vis_StateData (state,1,&n,&temperature);
      printf("%5d %12.4e\n",n,temperature);
      vis_LCaseAddConc (lcase,n,LCASE_TEMP,&temperature);
   }
   printf("Finished transfering nodal temperatures\n");

                   /* create LCase for reference temperature if nonzero */
   lcase = vis_LCaseBegin ();
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,(Vobject*)gridfun);
   vsy_HashTableInsert (lch,lcasetref,(Vobject*)lcase);
   for(n = 1; n <= MAX_NODE; ++n) {
      vis_LCaseAddConc (lcase,n,LCASE_TEMP,&tref);
   }
   printf("Finished setting reference temperature\n");
                   /* initialize ProSolve for different analysis */
   vfx_ProSolveInit (prosolve);

                   /* run ProSolve for structural problem */
   vfx_ProSolveRange (prosolve,2,2);
   vfx_ProSolveExec (prosolve);
   printf("Finished running structural analysis\n");

                   /* retrieve structural solution */
   vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
   vis_RPropSetType (rprop,SYS_RES_D);
   vis_StateDef (state,MAX_NODE,SYS_NODE,SYS_NONE,SYS_SIXDOF);
   vfx_ProSolveState (prosolve,rprop,state);
   printf("Displacements\n");
   for(n = 1; n <= MAX_NODE; ++n) {
      vis_StateData (state,1,&n,data);
      printf("%5d %12.4e %12.4e %12.4e\n",n,data[0],data[1],data[2]);
   }
   printf("\n");

   vis_RPropDef (rprop,SYS_ELEM,SYS_NODE);
   vis_RPropSetType (rprop,SYS_RES_S);
   vis_StateDef (state,MAX_ELEM,SYS_ELEM,SYS_NODE,SYS_TENSOR);
   vfx_ProSolveState (prosolve,rprop,state);
   printf("Stresses\n");
   for(n = 1; n <= MAX_ELEM; ++n) {
      vis_StateData (state,1,&n,data);
      printf("%5d\n",n);
      for(i = 0; i < 8; ++i) {
         printf("   %5d %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n",
                n,data[6*i  ],data[6*i+1],data[6*i+2],
                  data[6*i+3],data[6*i+4],data[6*i+5]);
      }
   }
   printf("\n");

   vis_RPropSetType (rprop,SYS_RES_E);
   vis_StateDef (state,MAX_ELEM,SYS_ELEM,SYS_NODE,SYS_TENSOR);
   vfx_ProSolveState (prosolve,rprop,state);
   printf("Strains\n");
   for(n = 1; n <= MAX_ELEM; ++n) {
      vis_StateData (state,1,&n,data);
      printf("%5d\n",n);
      for(i = 0; i < 8; ++i) {
         printf("   %5d %12.4e %12.4e %12.4e %12.4e %12.4e %12.4e\n",
                n,data[6*i  ],data[6*i+1],data[6*i+2],
                  data[6*i+3],data[6*i+4],data[6*i+5]);
      }
   }
   printf("\n");

                   /* end objects */ 
   vis_RPropEnd (rprop);
   vis_StateEnd (state);
   vis_ModelDelete (model);
   vis_ModelEnd (model);
   vis_GridFunEnd (gridfun);
   vfx_ProSolveEnd (prosolve);
   return 0;
}

5.9. Example 10, Buckling Analysis

This example illustrates using ProSolve to perform a buckling analysis.

Like the thermal-structural analysis in Example 8, two SProp objects are required: the first performs a static analysis to compute the pre-stress condition; the second computes the resulting buckling modes.

Note that unlike Example 8 there is no need to call vfx_ProSolveInit() between the two SProp analysis. As a result a single vfx_ProSolveExec() command is requested using a vfx_ProSolveRange() from 1 to 2.

ProSolve only holds output for a single step. As a result a monitor function is introduced to retrieve State results when more than one step is computed. The example requests two buckling modes to demonstrate how the monitor function enables not only the retrieval of results for the static and the buckling SProp s, but also for the retrieval of multiple modes from the buckling analysis. An identical procedure could be used if, instead of buckling, a vibration analysis were requested.

#include "vis/vis.h"
#include "vfx/vfx.h"
/*
                 4------5------6
       r        /.     /.     /.
       e       / .    / .    / .
       s      /  .   /  .   /  .
       t     /   1../...2../...3
       r    10-----11-----12  /    +--- pressure 50.
       a    |  .   |  .   |  /
       i    | .    | .    | /
       t    |.     |.     |/
            7------8------9

         y
         |
         --- x
       /
      z
*/
#define MAX_ELEM  2
#define MAX_NODE 12

                   /* element connectivity */ 
static Vint conn[MAX_ELEM][8] = {
   { 1, 2, 5, 4, 7, 8,11,10 },
   { 2, 3, 6, 5, 8, 9,12,11 } };

                   /* element topologies */ 
static Vint shap[MAX_ELEM] = { VIS_SHAPEHEX, VIS_SHAPEHEX };
static Vint maxi[MAX_ELEM] = { 2, 2 };
                   /* element types */ 
static Vint featype[MAX_ELEM] = { VIS_ELEM_SOLID, VIS_ELEM_SOLID };

                   /* element property ids */ 
static Vint pid[MAX_ELEM] = { 1, 1 };
                   /* node coordinates */ 
static Vdouble coords[MAX_NODE][3] = {
   { 0.,0.,0.}, { 1.,0.,0.}, { 2.,0.,0.},
   { 0.,1.,0.}, { 1.,1.,0.}, { 2.,1.,0.},
   { 0.,0.,2.}, { 1.,0.,2.}, { 2.,0.,2.},
   { 0.,1.,2.}, { 1.,1.,2.}, { 2.,1.,2.} };

static void
mymonitor(vfx_ProSolve *prosolve, Vobject *o)
{
   Vint phase, numres;
   vis_RProp   *rprop;
   vis_State   *state;
   vis_GridFun *gridfun;
   vis_Connect *connect;
   vis_Model   *model;
   Vint n, enttype, subtype, id1, id2, id3, nument, datatype;
   Vint i, j, k, size, flag, objtype, dtyp, nval, num;
   Vfloat data[100];
   Vint idata[100];
   Vchar cdata[100];
   Vchar name[65];

   vfx_ProSolveGetInteger (prosolve,PROSOLVE_PHASE,&phase);
   if(phase == SYS_PHASE_OUTPUT) {

      vfx_ProSolveGetInteger (prosolve,PROSOLVE_NUMRES,&numres);
      rprop = vis_RPropBegin ();
      state = vis_StateBegin ();
      gridfun = vis_GridFunBegin ();
      vfx_ProSolveGetObject (prosolve,VIS_MODEL,(Vobject**)&model);
      vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
      vis_ConnectGridFun (connect,gridfun);
      vis_StateSetObject (state,VIS_GRIDFUN,gridfun);

      for(n = 1; n <= numres; n++) {
         printf("\nnew result:\n");
         vfx_ProSolveFillResRProp (prosolve,n,rprop,&flag,&objtype);
         vis_RPropInq (rprop,&enttype,&subtype);
         vis_RPropGetIds (rprop,&id1,&id2,&id3);
         printf("   id1= %d, id2= %d\n",id1,id2);
         if(enttype == SYS_NODE) {
            printf("   nodal result\n");
         } else {
            if(subtype == SYS_NODE) {
               printf("   element-nodal result\n");
            } else {
               printf("   element centroidal result\n");
            }
         }
         vfx_ProSolveState (prosolve,rprop,state);
         vis_StateInq (state,&nument,&enttype,&subtype,&datatype);
         for(i = 0; i < RPROP_MAX; i++) {
            vis_RPropValueFlag (rprop,i,&flag);
            if(flag == 0)  continue;

            vis_RPropValueName (rprop,i,name);
            printf("   %s=",name);
            vis_RPropValueParams (rprop,i,&nval,&dtyp);
            if(dtyp == SYS_CHAR) {
               vis_RPropValueString (rprop,i,cdata);
               printf(" %s\n",cdata);
            } else if(dtyp == SYS_INTEGER) {
               vis_RPropValueInteger (rprop,i,idata);
               for(j = 0; j < nval; j++) {
                  printf(" %d",idata[j]);
               }
               printf("\n");
            } else if(dtyp == SYS_FLOAT) {
               vis_RPropValueFloat (rprop,i,data);
               for(j = 0; j < nval; j++) {
                  printf(" % 12.5e",data[j]);
               }
               printf("\n");
            }
         }
         for(i = 1; i <= nument; i++) {
            vis_StateData (state,1,&i,data);
            if(datatype == SYS_SCALAR) {
               size = 1;
            } else if(datatype == SYS_VECTOR) {
               size = 3;
            } else if(datatype == SYS_SIXDOF) {
               size = 6;
            } else if(datatype == SYS_TENSOR) {
               size = 6;
            }
            num = 1;
            if(enttype == SYS_NODE) {
               printf("   node= %3d",i);
            } else if(enttype == SYS_ELEM) {
               printf("   elem= %3d\n",i);
               if(subtype == SYS_NODE) {
                  vis_ConnectElemNum (connect,SYS_NODE,i,&num);
               }
            }
            for(k = 0; k < num; k++) {
               if(num > 1) {
                  printf("      node= %2d, ",k+1);
               } else {
                  if(enttype == SYS_ELEM) {
                     printf("      centroidal value= ");
                  } else {
                     printf("      value= ");
                  }
               }
               for(j = 0; j < size; j++) {
                  printf(" % 12.5e",data[j]);
               }
               printf("\n");
            }
         }
      }
      vis_RPropEnd (rprop);
      vis_StateEnd (state);
      vis_GridFunEnd (gridfun);
   }
}

/*----------------------------------------------------------------------
                      Buckling Analysis
----------------------------------------------------------------------*/
int
main(int argc, char *argv[])
{
   vfx_ProSolve *prosolve;
   Vint i;
   vis_Model *model;
   vis_Connect *connect;
   vis_GridFun *gf;
   vis_MProp *mprop;
   vis_EProp *eprop;
   vis_SProp *sprop;
   vis_RCase *rcase;
   vis_LCase *lcase;
   vsy_HashTable *ephash, *mphash, *rchash, *lchash;
   vsy_List *splist;
   Vdouble pres[4];
   Vint res[2];

                   /* create connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,MAX_NODE,MAX_ELEM);

                   /* set node coordinates and user ids */ 
   for(i = 0; i < MAX_NODE; i++) {
      vis_ConnectSetCoordsdv (connect,i+1,coords[i]);
   }
                   /* set element topologies */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetTopology (connect,i+1,shap[i],maxi[i],0,0);
   }
                   /* set element node connectivity and associations */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetElemNode (connect,i+1,conn[i]);
      vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,i+1,featype[i]);
      vis_ConnectSetElemAssoc (connect,VIS_PROPID,i+1,pid[i]);
   }

                   /* create model object hierarchy */ 
   model = vis_ModelBegin ();

                   /* hashtables of element and material properties */
   ephash = vsy_HashTableBegin ();
   mphash = vsy_HashTableBegin ();

                   /* hashtables of restraint and load cases */
   rchash = vsy_HashTableBegin ();
   lchash = vsy_HashTableBegin ();

                   /* list of solution procedures */ 
   splist = vsy_ListBegin ();

                   /* material 1 */ 
   mprop = vis_MPropBegin ();
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   vis_MPropSetValued (mprop,MPROP_E,1.e+7);
   vis_MPropSetValued (mprop,MPROP_NU,.3);
   vsy_HashTableInsert (mphash,1,mprop);

                   /* solid property 1 */ 
   eprop = vis_EPropBegin ();
   vis_EPropDef (eprop,VIS_ELEM_SOLID);
   vis_EPropSetValuei (eprop,EPROP_MID,1);
   vsy_HashTableInsert (ephash,1,eprop);

                   /* GridFun for use with LCase objects with element loads */ 
   gf = vis_GridFunBegin ();
   vis_ConnectGridFun (connect,gf);

                   /* restraint case 1 */ 
                   /* fix nodes 1, 4, 7, 10 in x */
                   /* fix node 4 in z, node 7 in y */
   rcase = vis_RCaseBegin ();
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TY,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TZ,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TZ,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,7,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,7,SYS_DOF_TY,RCASE_FIXED,NULL,0);
   vis_RCaseSetSPC (rcase,10,SYS_DOF_TX,RCASE_FIXED,NULL,0);
   vsy_HashTableInsert (rchash,1,rcase);

                   /* load case 1 */ 
                   /* uniform pressure on face 6 element 2 */ 
                   /* register GridFun object with LCase */
   lcase = vis_LCaseBegin ();
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,gf);
   pres[0] = 50.;
   pres[1] = 50.;
   pres[2] = 50.;
   pres[3] = 50.;
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,6,LCASE_PRES,pres);
   vsy_HashTableInsert (lchash,1,lcase);

                   /* solution step 1 - static analysis */
   sprop = vis_SPropBegin ();
   vis_SPropDef (sprop,SYS_SOL_STATIC);
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
   vis_SPropSetValuei (sprop,SPROP_CASEID,1);
   vis_SPropSetValuei (sprop,SPROP_RCASE,1);
   vis_SPropSetValued (sprop,SPROP_RCASE_FACTOR,1.);
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
   vis_SPropSetValuei (sprop,SPROP_LCASE,1);
   vis_SPropSetValued (sprop,SPROP_LCASE_FACTOR,1.);
                   /* output requests, displacement, stress */ 
   res[0] = SYS_RES_D;
   res[1] = SYS_RES_S;
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,2);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,res);
   vsy_ListInsert (splist,1,sprop);

                   /* solution step 2 - buckling analysis */
   sprop = vis_SPropBegin ();
   vis_SPropDef (sprop,SYS_SOL_BUCKLING);
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_STRUCTURAL);
   vis_SPropSetValuei (sprop,SPROP_CASEID,2);
   vis_SPropSetValuei (sprop,SPROP_EIGEN_NUM,2);
   vis_SPropSetValuei (sprop,SPROP_EIGEN_TYPE,SYS_EIGEN_NEAREST);
   vis_SPropSetValued (sprop,SPROP_EIGEN_SHIFT,0.);

                   /* output requests, displacement, stress */ 
   res[0] = SYS_RES_D;
   res[1] = SYS_RES_S;
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,2);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,res);
   vsy_ListInsert (splist,2,sprop);

                   /* register Connect in Model */ 
   vis_ModelSetObject (model,VIS_CONNECT,connect);

                   /* register property hashtables in Model */ 
   vis_ModelSetHashTable (model,VIS_MPROP,mphash);
   vis_ModelSetHashTable (model,VIS_EPROP,ephash);

                   /* register case hashtables in Model */ 
   vis_ModelSetHashTable (model,VIS_RCASE,rchash);
   vis_ModelSetHashTable (model,VIS_LCASE,lchash);

                   /* register solution property list in Model */ 
   vis_ModelSetList (model,VIS_SPROP,splist);

                   /* create and configure ProSolve object */
   prosolve = vfx_ProSolveBegin ();
   vfx_ProSolveSetFunction (prosolve,PROSOLVE_FUN_MONITOR,mymonitor,NULL);
   vfx_ProSolveRange (prosolve,1,2);
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);

                   /* solve */
   vfx_ProSolveExec (prosolve);

                   /* end objects */ 
   vis_GridFunEnd (gf);
                   /* Use convenience routine ModelDelete to destroy
                      all objects registered in Model */
   vis_ModelDelete (model);
                   /* finally destroy Model object */ 
   vis_ModelEnd (model);
   vfx_ProSolveEnd (prosolve);
   return 0;
}

5.10. Example 11, Thermal Analysis

This example illustrates using ProSolve to perform a thermal analysis.

Boundary conditions include applied temperature, applied heat flux, convection and volumetric heat generation. A single SProp object is needed to specify the solution step.

#include "vis/vis.h"
#include "vfx/vfx.h"
/*

                 4------5------6
       a        /.     /.     /.
       p       / .    / .    / .
       p      /  .   /  .   /  .
       l     /   1../...2../...3
       i    10-----11-----12  /    +--- convection h = 1., T = 100.
       e    |  .   |  .   |  /     +--- F=5 = heat flux on all side faces
       d    | .    | .    | /      +--- T=10 = applied temperature on left face
            |.     |.     |/       +--- volumetric heat V=50 on both elements
       T    7------8------9

         y
         |
         --- x
       /
      z
*/
#define MAX_ELEM  2
#define MAX_NODE 12

                   /* element connectivity */ 
static Vint conn[MAX_ELEM][8] = {
   { 1, 2, 5, 4, 7, 8,11,10 },
   { 2, 3, 6, 5, 8, 9,12,11 } };

                   /* element topologies */ 
static Vint shap[MAX_ELEM] = { VIS_SHAPEHEX, VIS_SHAPEHEX };
static Vint maxi[MAX_ELEM] = { 2, 2 };
                   /* element types */ 
static Vint featype[MAX_ELEM] = { VIS_ELEM_SOLID, VIS_ELEM_SOLID };

                   /* element property ids */ 
static Vint pid[MAX_ELEM] = { 1, 1 };
                   /* node coordinates */ 
static Vdouble coords[MAX_NODE][3] = {
   { 0.,0.,0.}, { 1.,0.,0.}, { 2.,0.,0.},
   { 0.,1.,0.}, { 1.,1.,0.}, { 2.,1.,0.},
   { 0.,0.,2.}, { 1.,0.,2.}, { 2.,0.,2.},
   { 0.,1.,2.}, { 1.,1.,2.}, { 2.,1.,2.} };

static void
mymonitor(vfx_ProSolve *prosolve, Vobject *o)
{
   Vint phase, numres;
   vis_RProp   *rprop;
   vis_State   *state;
   vis_GridFun *gridfun;
   vis_Connect *connect;
   vis_Model   *model;
   Vint n, enttype, subtype, id1, id2, id3, nument, datatype;
   Vint i, j, k, size, flag, objtype, dtyp, nval, num;
   Vfloat data[100];
   Vint idata[100];
   Vchar cdata[100];
   Vchar name[65];

   vfx_ProSolveGetInteger (prosolve,PROSOLVE_PHASE,&phase);
   if(phase == SYS_PHASE_OUTPUT) {
      vfx_ProSolveGetInteger (prosolve,PROSOLVE_NUMRES,&numres);
      rprop = vis_RPropBegin ();
      state = vis_StateBegin ();
      gridfun = vis_GridFunBegin ();
      vfx_ProSolveGetObject (prosolve,VIS_MODEL,(Vobject**)&model);
      vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
      vis_ConnectGridFun (connect,gridfun);
      vis_StateSetObject (state,VIS_GRIDFUN,gridfun);

      for(n = 1; n <= numres; n++) {
         printf("\nnew result:\n");
         vfx_ProSolveFillResRProp (prosolve,n,rprop,&objtype,&flag);
         vis_RPropInq (rprop,&enttype,&subtype);
         vis_RPropGetIds (rprop,&id1,&id2,&id3);
         printf("   id1= %d, id2= %d\n",id1,id2);
         if(enttype == SYS_NODE) {
            printf("   nodal result\n");
         } else {
            if(subtype == SYS_NODE) {
               printf("   element-nodal result\n");
            } else {
               printf("   element centroidal result\n");
            }
         }
         vfx_ProSolveState (prosolve,rprop,state);
         vis_StateInq (state,&nument,&enttype,&subtype,&datatype);
         for(i = 0; i < RPROP_MAX; i++) {
            vis_RPropValueFlag (rprop,i,&flag);
            if(flag == 0)  continue;

            vis_RPropValueName (rprop,i,name);
            printf("   %s=",name);
            vis_RPropValueParams (rprop,i,&nval,&dtyp);
            if(dtyp == SYS_CHAR) {
               vis_RPropValueString (rprop,i,cdata);
               printf(" %s\n",cdata);
            } else if(dtyp == SYS_INTEGER) {
               vis_RPropValueInteger (rprop,i,idata);
               for(j = 0; j < nval; j++) {
                  printf(" %d",idata[j]);
               }
               printf("\n");
            } else if(dtyp == SYS_FLOAT) {
               vis_RPropValueFloat (rprop,i,data);
               for(j = 0; j < nval; j++) {
                  printf(" % 12.5e",data[j]);
               }
               printf("\n");
            }
         }
         for(i = 1; i <= nument; i++) {
            vis_StateData (state,1,&i,data);
            if(datatype == SYS_SCALAR) {
               size = 1;
            } else if(datatype == SYS_VECTOR) {
               size = 3;
            } else if(datatype == SYS_SIXDOF) {
               size = 6;
            } else if(datatype == SYS_TENSOR) {
               size = 6;
            }
            num = 1;
            if(enttype == SYS_NODE) {
               printf("   node= %3d",i);
            } else if(enttype == SYS_ELEM) {
               printf("   elem= %3d\n",i);
               if(subtype == SYS_NODE) {
                  vis_ConnectElemNum (connect,SYS_NODE,i,&num);
               }
            }
            for(k = 0; k < num; k++) {
               if(num > 1) {
                  printf("      node= %2d, ",k+1);
               } else {
                  if(enttype == SYS_ELEM) {
                     printf("      centroidal value= ");
                  } else {
                     printf("      value= ");
                  }
               }
               for(j = 0; j < size; j++) {
                  printf(" % 12.5e",data[j]);
               }
               printf("\n");
            }
         }
      }
      vis_RPropEnd (rprop);
      vis_StateEnd (state);
      vis_GridFunEnd (gridfun);
   }
}

/*----------------------------------------------------------------------
                      Thermal Analysis
----------------------------------------------------------------------*/
int
main(int argc, char *argv[])
{
   vfx_ProSolve *prosolve;
   Vint i;
   vis_Model *model;
   vis_Connect *connect;
   vis_GridFun *gf;
   vis_MProp *mprop;
   vis_EProp *eprop;
   vis_SProp *sprop;
   vis_RCase *rcase;
   vis_LCase *lcase;
   vsy_HashTable *ephash, *mphash, *rchash, *lchash;
   vsy_List *splist;
   Vdouble film[4][2], faceheat[4], volheat[8];
   Vint res[3];
   Vfloat ten = 10.;

                   /* create connect object */
   connect = vis_ConnectBegin ();
   vis_ConnectPre (connect,SYS_DOUBLE);
   vis_ConnectDef (connect,MAX_NODE,MAX_ELEM);

                   /* set node coordinates and user ids */ 
   for(i = 0; i < MAX_NODE; i++) {
      vis_ConnectSetCoordsdv (connect,i+1,coords[i]);
   }
                   /* set element topologies */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetTopology (connect,i+1,shap[i],maxi[i],0,0);
   }
                   /* set element node connectivity and associations */ 
   for(i = 0; i < MAX_ELEM; i++) {
      vis_ConnectSetElemNode (connect,i+1,conn[i]);
      vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,i+1,featype[i]);
      vis_ConnectSetElemAssoc (connect,VIS_PROPID,i+1,pid[i]);
   }

                   /* create model object hierarchy */ 
   model = vis_ModelBegin ();

                   /* hashtables of element and material properties */
   ephash = vsy_HashTableBegin ();
   mphash = vsy_HashTableBegin ();

                   /* hashtables of restraint and load cases */
   rchash = vsy_HashTableBegin ();
   lchash = vsy_HashTableBegin ();

                   /* list of solution procedures */ 
   splist = vsy_ListBegin ();

                   /* material 1 */ 
   mprop = vis_MPropBegin ();
   vis_MPropDef (mprop,SYS_MAT_ISOTROPIC);
   vis_MPropSetValued (mprop,MPROP_K,1.);
   vsy_HashTableInsert (mphash,1,mprop);

                   /* solid property 1 */ 
   eprop = vis_EPropBegin ();
   vis_EPropDef (eprop,VIS_ELEM_SOLID);
   vis_EPropSetValuei (eprop,EPROP_MID,1);
   vsy_HashTableInsert (ephash,1,eprop);

                   /* GridFun for use with LCase objects with element loads */ 
   gf = vis_GridFunBegin ();
   vis_ConnectGridFun (connect,gf);

                   /* restraint case 1 */ 
                   /* fix nodes 1, 4, 7, 10 in x */
                   /* fix node 4 in z, node 7 in y */
   rcase = vis_RCaseBegin ();
   vis_RCaseSetSPC (rcase,1,SYS_DOF_TEMP,RCASE_APPLIED,&ten,0);
   vis_RCaseSetSPC (rcase,4,SYS_DOF_TEMP,RCASE_APPLIED,&ten,0);
   vis_RCaseSetSPC (rcase,7,SYS_DOF_TEMP,RCASE_APPLIED,&ten,0);
   vis_RCaseSetSPC (rcase,10,SYS_DOF_TEMP,RCASE_APPLIED,&ten,0);
   vsy_HashTableInsert (rchash,1,rcase);

                   /* load case 1 */ 
                   /* register GridFun object with LCase */
   lcase = vis_LCaseBegin ();
   vis_LCaseSetObject (lcase,VIS_GRIDFUN,gf);

                   /* convection on face 6 element 2 */ 
   film[0][0] = 100.;
   film[1][0] = 100.;
   film[2][0] = 100.;
   film[3][0] = 100.;
   film[0][1] = 1.;
   film[1][1] = 1.;
   film[2][1] = 1.;
   film[3][1] = 1.;
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,6,LCASE_HEATCONV,(Vdouble*)film);

                   /* distributed heat on side faces of elements 1 and 2 */ 
   faceheat[0] = 5.;
   faceheat[1] = 5.;
   faceheat[2] = 5.;
   faceheat[3] = 5.;
   vis_LCaseSetDistdv (lcase,SYS_FACE,1,1,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,1,2,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,1,3,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,1,4,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,1,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,2,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,3,LCASE_HEATFLUX,faceheat);
   vis_LCaseSetDistdv (lcase,SYS_FACE,2,4,LCASE_HEATFLUX,faceheat);

                   /* volumetric heat on elements 1 and 2 */ 
   volheat[0] = 50.;
   volheat[1] = 50.;
   volheat[2] = 50.;
   volheat[3] = 50.;
   volheat[4] = 50.;
   volheat[5] = 50.;
   volheat[6] = 50.;
   volheat[7] = 50.;
   vis_LCaseSetElemdv (lcase,1,LCASE_HEATELEM,volheat);
   vis_LCaseSetElemdv (lcase,2,LCASE_HEATELEM,volheat);
   vsy_HashTableInsert (lchash,1,lcase);

                   /* static (steady state) thermal analysis */
   sprop = vis_SPropBegin ();
   vis_SPropDef (sprop,SYS_SOL_STATIC);
   vis_SPropSetValuei (sprop,SPROP_ANALYSIS,SYS_ANALYSIS_THERMAL);
   vis_SPropSetValuei (sprop,SPROP_CASEID,1);
   vis_SPropSetValuei (sprop,SPROP_RCASE,1);
   vis_SPropSetValued (sprop,SPROP_RCASE_FACTOR,1.);
   vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
   vis_SPropSetValuei (sprop,SPROP_LCASE,1);
   vis_SPropSetValued (sprop,SPROP_LCASE_FACTOR,1.);
                   /* output requests, displacement, stress */ 
   res[0] = SYS_RES_TEMP;
   res[1] = SYS_RES_TEMP_GRAD;
   res[2] = SYS_RES_HEAT_FLUX;
   vis_SPropSetValuei (sprop,SPROP_RESFILE_NUM,3);
   vis_SPropSetValueiv (sprop,SPROP_RESFILE,res);
   vsy_ListInsert (splist,1,sprop);

                   /* register Connect in Model */ 
   vis_ModelSetObject (model,VIS_CONNECT,connect);

                   /* register property hashtables in Model */ 
   vis_ModelSetHashTable (model,VIS_MPROP,mphash);
   vis_ModelSetHashTable (model,VIS_EPROP,ephash);

                   /* register case hashtables in Model */ 
   vis_ModelSetHashTable (model,VIS_RCASE,rchash);
   vis_ModelSetHashTable (model,VIS_LCASE,lchash);

                   /* register solution property list in Model */ 
   vis_ModelSetList (model,VIS_SPROP,splist);

                   /* create and configure ProSolve object */
   prosolve = vfx_ProSolveBegin ();
   vfx_ProSolveSetFunction (prosolve,PROSOLVE_FUN_MONITOR,mymonitor,NULL);
   vfx_ProSolveRange (prosolve,1,1);
   vfx_ProSolveSetObject (prosolve,VIS_MODEL,model);

                   /* solve */
   vfx_ProSolveExec (prosolve);

                   /* end objects */ 
   vis_GridFunEnd (gf);
                   /* Use convenience routine ModelDelete to destroy
                      all objects registered in Model */
   vis_ModelDelete (model);
                   /* finally destroy Model object */ 
   vis_ModelEnd (model);
   vfx_ProSolveEnd (prosolve);
   return 0;
}

5.11. Example 13, Frequency Response Analysis

This example illustrates using ProSolve to perform a mode based frequency response analysis.

A mode based frequency response analysis requires two solution procedures. The first solution step is essentially a superelement analysis to compute a set of vibration modes and attachment modes. The attachment modes represent degrees of freedom which are meant to have excitations applied. A reduced mass and stiffness matrix are computed in this first step to be used in the second solution procedure.

The second solution step is the frequency repsonse analysis step in which the exitations at the attachment points are defined and the frequencies at which responses are to be calculated are defined.

Note that the monitor function ignores output during the first solution step.

#include "vis/vis.h"
#include "vfx/vfx.h"
/*
                 4------5------6
       r        /.     /.     /.
       e       / .    / .    / .
       s      /  .   /  .   /  .
       t     /   1../...2../...3
       r    10-----11-----12  /    +--- pressure 50.
       a    |  .   |  .   |  /
       i    | .    | .    | /
       t    |.     |.     |/
            7------8------9

         y
         |
         --- x
       /
      z
*/
#define MAX_ELEM  2
#define MAX_NODE 12

static Vfloat indep1[2] = {1.,10000.};
static Vfloat dep1[2]   = {1.,1000.};
static Vfloat indep2[2] = {1.,10000.};
static Vfloat dep2[2]   = {1.,2000.};

static Vfloat zero = 0.;
static Vfloat ones[2] = {1.,1.};

                   /* element connectivity */ 
static Vint conn[MAX_ELEM][8] = {
   { 1, 2, 5, 4, 7, 8,11,10 },
   { 2, 3, 6, 5, 8, 9,12,11 } };

                   /* element topologies */ 
static Vint shap[MAX_ELEM] = { VIS_SHAPEHEX, VIS_SHAPEHEX };
static Vint maxi[MAX_ELEM] = { 2, 2 };
                   /* element types */ 
static Vint featype[MAX_ELEM] = { VIS_ELEM_SOLID, VIS_ELEM_SOLID };

                   /* element property ids */ 
static Vint pid[MAX_ELEM] = { 1, 1 };
                   /* node coordinates */ 
static Vdouble coords[MAX_NODE][3] = {
   { 0.,0.,0.}, { 1.,0.,0.}, { 2.,0.,0.},
   { 0.,1.,0.}, { 1.,1.,0.}, { 2.,1.,0.},
   { 0.,0.,2.}, { 1.,0.,2.}, { 2.,0.,2.},
   { 0.,1.,2.}, { 1.,1.,2.}, { 2.,1.,2.} };

/*----------------------------------------------------------------------
                      Monitor function
----------------------------------------------------------------------*/
static void
monfun(vfx_ProSolve *p, Vobject *obj)
{
   vis_Model *model;
   vis_State *state;
   vsy_List  *spl;
   vis_SProp *sprop;
   vis_RProp *rprop;
   vis_Connect *connect;
   vis_GridFun *gridfun;
   Vint phase, id, caseid, numfile, n, *restype, complexflag;
   Vint numel, numnp;

                   /* extract and print phase and percent done */ 
   vfx_ProSolveGetInteger (p,PROSOLVE_PHASE,&phase);

                   /* extract multiple output */ 
   if(phase == SYS_PHASE_OUTPUT) {
      /* check for frequency response solution property only */
      vfx_ProSolveGetInteger (p,PROSOLVE_CASE,&caseid);
      if(caseid == 1)  return;

      vfx_ProSolveGetInteger (p,PROSOLVE_STEP,&id);
      printf("   output available for caseid= %d, step= %d\n",caseid,id);
                   /* extract output list for this subcase */
      vfx_ProSolveGetObject (p,VIS_MODEL,(Vobject**)&model);
      vis_ModelGetList (model,VIS_SPROP,&spl);
      vsy_ListRef (spl,caseid,(Vobject**)&sprop);
      vis_SPropValueInteger (sprop,SPROP_RESFILE_NUM,&numfile);
      if(numfile) {
                   /* extract number of nodes and elements */
         vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
         gridfun = vis_GridFunBegin ();
         vis_ConnectGridFun (connect,gridfun);
         state = vis_StateBegin ();
         vis_StateSetObject (state,VIS_GRIDFUN,(Vobject*)gridfun);
         vis_ConnectNumber (connect,SYS_NODE,&numnp);
         vis_ConnectNumber (connect,SYS_ELEM,&n