The material properties are isotropic linear elastic and characteristic of aluminum. A LinMat object is instanced and defined to contain an isotropic material. The Young's modulus, Poisson's ratio and density of the material are set. A MatlFun object is instanced and filled with the interface to the LinMat object using vfe_LinMatMatlFun. A Solid3D element formulation object is instanced and configured to generate element matrices and vectors for isoparametric parabolic Serendipity tetrahedra (10 nodes). The MatlFun object is registered with Solid3D using vfe_Solid3DSetObject. This completes the connection between the element formulation and the material model. Element quantities such as stiffness, mass, etc. can now be computed.
The element stiffness matrix is generated using vfe_Solid3DStiff The lower triangle of the matrix is returned in array kl. All isoparametric element formulations supported by Solid3D use 3 translations per node. The total number of degrees of freedom for the 10 node tetrahedron is 30. The output arrays for matrices and vectors have been dimensioned exactly to accommodate the output for this element. Both the consistent and diagonal mass matrices are generated.
Two varieties of distributed loads are computed; a constant pressure and a variable force per unit area directed along the global x axis. Each load is applied to face 3 of the tetrahedron. Note that the function vfe_Solid3DDistLoad which is used to compute consistent nodal loads, requires the element node coordinate array, x, and element node load magnitudes, v, even though only the nodes on face 3 will be referenced. The output load vector, f, contains the computed consistent loads for each element degree of freedom.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[10][3] = {
{0.,0.,0.}, {2.,0.,0.}, {0.,2.,0.}, {0.,0.,2.},
{1.,0.,0.}, {1.,1.,0.}, {0.,1.,0.},
{0.,0.,1.}, {1.,0.,1.}, {0.,1.,1.} };
/*----------------------------------------------------------------------
Generate Stiffness Matrix for 10 Node Tetrahedron
----------------------------------------------------------------------*/
int
main()
{
Vint i, j;
Vdouble prop[2];
Vdouble kl[465];
Vdouble mc[465], md[30];
Vdouble v[10], f[30];
vfe_Solid3D *solid3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 10000000.;
prop[1] = .3;
vfe_LinMatSetElasProp (linmat,prop);
vfe_LinMatSetDensity (linmat,.0000133);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create 3D solid element formulation */
solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetParami (solid3d,VFE_TECH,VFE_TECH_ISOP);
vfe_Solid3DSetTopology (solid3d,SYS_SHAPETET,3,0,0);
/* Register material functions */
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
/* Form stiffness matrix */
vfe_Solid3DStiff (solid3d,x,kl);
/* Form consistent mass matrix */
vfe_Solid3DMass (solid3d,x,mc);
/* Form diagonal mass matrix */
vfe_Solid3DMassDiag (solid3d,x,md);
printf("Diagonal mass matrix\n");
for(i = 0; i < 10; i++) {
for(j = 0; j < 3; j++) {
printf("%12.4e ",md[3*i+j]);
}
printf("\n");
}
/* Form load vector for constant pressure */
for(i = 0; i < 10; i++) v[i] = 1.;
vfe_Solid3DDistLoad (solid3d,x,SYS_FACE,3,VFE_DISTLOAD_PRES,v,f);
printf("Constant pressure load vector\n");
for(i = 0; i < 10; i++) {
for(j = 0; j < 3; j++) {
printf("%12.4e ",f[3*i+j]);
}
printf("\n");
}
/* Delete objects */
vfe_Solid3DEnd (solid3d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
return 0;
}
#include <stdlib.h>
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[10][3] = {
{0.,0.,0.}, {2.,0.,0.}, {0.,2.,0.}, {0.,0.,2.},
{1.,0.,0.}, {1.,1.,0.}, {0.,1.,0.},
{0.,0.,1.}, {1.,0.,1.}, {0.,1.,1.} };
/*----------------------------------------------------------------------
Query for Degree of Freedom Map for 10 Node Tetrahedron
----------------------------------------------------------------------*/
int
main()
{
Vint i;
Vdouble prop[2];
Vdouble *kl;
Vint nedofs;
Vint loc[VFE_MAXDOF], tag[VFE_MAXDOF];
vfe_Solid3D *solid3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 10000000.;
prop[1] = .3;
vfe_LinMatSetElasProp (linmat,prop);
vfe_LinMatSetDensity (linmat,.0000133);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create 3D solid element formulation */
solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetParami (solid3d,VFE_TECH,VFE_TECH_ISOP);
vfe_Solid3DSetTopology (solid3d,SYS_SHAPETET,3,0,0);
/* Register material functions */
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
/* Query for element degree of freedom map */
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
/* Set tag names */
printf("Degree of freedom map\n");
for(i = 0; i < nedofs; i++) {
if(tag[i]==SYS_DOF_TX) printf("loc[%2d]= %2d, tag[%2d]= TX\n",i,loc[i],i);
if(tag[i]==SYS_DOF_TY) printf("loc[%2d]= %2d, tag[%2d]= TY\n",i,loc[i],i);
if(tag[i]==SYS_DOF_TZ) printf("loc[%2d]= %2d, tag[%2d]= TZ\n",i,loc[i],i);
}
/* Allocate storage for element stiffness */
kl = (Vdouble*) malloc(sizeof(Vdouble)*(nedofs*(nedofs+1)/2));
/* Form stiffness matrix */
vfe_Solid3DStiff (solid3d,x,kl);
/* Free allocated storage for element stiffness */
free(kl);
/* Delete objects */
vfe_Solid3DEnd (solid3d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
return 0;
}
#include <stdio.h>
#include <stdlib.h>
#include "base/base.h"
#include "vfe/vfe.h"
static void
gat_vector(Vdouble x[][3], Vint np, Vint con[4], Vdouble xe[][3]);
static void
set_scalar(Vdouble s, Vint np, Vdouble se[]);
static Vint is[2][6] = {
{ 1, 7, 9, 4, 8, 5 },
{ 9, 3, 1, 6, 2, 5 } };
/*----------------------------------------------------------------------
2D Plane Stress Cantilever
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, k;
Vdouble prop[2];
Vdouble x[33][3];
Vint ix[10][6];
Vdouble xe[6][3], te[6];
Vint nelems;
Vint nedofs;
Vdouble *kl;
vfe_Solid2D *solid2d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
/* Generate 3 by 11 mesh of nodes */
k = 0;
for(j = 0; j < 11; j++) {
for(i = 0; i < 3; i++) {
x[k][0] = .2 * j;
x[k][1] = .1 * i;
x[k][2] = 0.;
k++;
}
}
/* set number of elements */
nelems = 10;
/* Generate 10 triangular elements */
k = 0;
for(j = 0; j < nelems; j += 2) {
for(i = 0; i < 6; i++) {
ix[j][i] = k+is[0][i];
ix[j+1][i] = k+is[1][i];
}
k += 6;
}
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 10000000.;
prop[1] = .3;
vfe_LinMatSetElasProp (linmat,prop);
vfe_LinMatSetDensity (linmat,.000254);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create 2D solid element formulation */
solid2d = vfe_Solid2DBegin ();
vfe_Solid2DSetParami (solid2d,VFE_TECH,VFE_TECH_ISOP);
vfe_Solid2DSetTopology (solid2d,SYS_SHAPETRI,3,0);
/* Register material functions */
vfe_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun);
/* Determine number of degrees of freedom */
vfe_Solid2DNumDof (solid2d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
/* Allocate memory for element stiffnesses */
kl = (Vdouble*) malloc(sizeof(Vdouble)*(nedofs*(nedofs+1)/2));
/* Form stiffness matrix element by element */
vfe_Solid2DSetPropPtr (solid2d,VFE_PROP_DEPTH,te);
for(i = 0; i < nelems; i++) {
set_scalar (.5,6,te);
gat_vector (x,6,ix[i],xe);
vfe_Solid2DStiff (solid2d,xe,kl);
/* print first and last element */
if(i == 0 || i == nelems-1) {
printf("Lower triangle of stiffness matrix\n");
for(k = 0; k < nedofs; k++) {
for(j = 0; j <= k; j++) {
printf("%12.4e ",kl[k*(k+1)/2+j]);
}
printf("\n");
}
}
}
/* Free memory */
free(kl);
/* Delete objects */
vfe_Solid2DEnd (solid2d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
return 0;
}
/*----------------------------------------------------------------------
Utility function to gather element vectors
----------------------------------------------------------------------*/
static void
gat_vector(Vdouble x[][3], Vint np, Vint con[4], Vdouble xe[][3])
{
int i;
for(i = 0; i < np; i++) {
xe[i][0] = x[con[i]-1][0];
xe[i][1] = x[con[i]-1][1];
xe[i][2] = x[con[i]-1][2];
}
}
/*----------------------------------------------------------------------
Utility function to set element scalars
----------------------------------------------------------------------*/
static void
set_scalar(Vdouble s, Vint np, Vdouble se[])
{
int i;
for(i = 0; i < np; i++) {
se[i] = s;
}
}
As in Example 1, the material properties are isotropic linear elastic and characteristic of aluminum. A LinMat object is instanced and defined to contain an isotropic material. The Young's modulus, Poisson's ratio and density of the material are set. A MatlFun object is instanced and filled with the interface to the LinMat object using vfe_LinMatMatlFun. At this point this example departs from Example 1 in that a ShellProp object must be used to specify the shell wall (through the shell thickness) construction. In this case a monocoque (single layer) shell wall is specified using vfe_ShellPropDef and vfe_ShellPropSetMonocoque. For linear material properties it is only necessary to use 2 point Gaussian quadrature to integrate through the shell thickness. A monocoque shell wall may only consist of a single material which is specified to be material 1. The MatlFun object associated with the linear material properties is registered with ShellProp as material 1 using vfe_ShellPropSetMatlFun. An additional MatlFun object will now be required to connect the ShellProp object to the Shell3D element formulation object. Therefore another MatlFun object is instanced and filled with the interface to the ShellProp object using vfe_ShellPropMatlFun.
A Shell3D element formulation object is instanced and configured to generate element matrices and vectors for linear quadrilateral shells (4 nodes) using an assumed natural strain technology. A preintegrated formulation (using moments and curvatures) is chosen.
The MatlFun object, matlfunprop, is registered with Shell3D using vfe_Shell3DSetObject. This completes the connections between the element formulation, the element properties and the material model. The total thickness at each element node is set to .1 and is entered into the array th. The pointer to this array is registered with the Shell3D object using vfe_Shell3DSetPropPtr. Any other element nodal properties such as offsets or shell normals would be specified in a similar manner. Element quantities such as stiffness, mass, etc. can now be computed.
The element stiffness matrix is generated using vfe_Shell3DStiff The lower triangle of the matrix is returned in array kl. All element formulations supported by Shell3D use 3 translations and 3 rotations per node. The total number of degrees of freedom for the 4 node quadrilateral shell is 24. The output arrays for matrices and vectors have been dimensioned exactly to accommodate the output for this element. Both the consistent and diagonal mass matrices are generated. A constant pressure load is computed using vfe_Shell3DDistLoad.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[4][3] = {
{0.,0.,0.}, {1.,0.,0.}, {1.,1.,0.}, {0.,1.,0.} };
/*----------------------------------------------------------------------
Generate Stiffness Matrix for 4 Node Quadrilateral Shell
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, k;
Vdouble prop[2];
Vdouble sprop[2];
Vdouble kl[300];
Vdouble mc[300], md[24];
Vdouble v[4], f[24];
vfe_Shell3D *shell3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
vfe_MatlFun *matlfunprop;
vfe_ShellProp *shellprop;
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 10000000.;
prop[1] = .3;
vfe_LinMatSetElasProp (linmat,prop);
vfe_LinMatSetDensity (linmat,.0000133);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create shell property object */
shellprop = vfe_ShellPropBegin ();
vfe_ShellPropDef (shellprop,SHELLPROP_MONOCOQUE);
sprop[0] = .1;
sprop[1] = 5./6.;
vfe_ShellPropSetMonocoque (shellprop,sprop,2,SYS_RULE_GAUSS,0.,1);
vfe_ShellPropSetMatlFun (shellprop,1,matlfun);
/* Create MatlFun object and load shell property functions */
matlfunprop = vfe_MatlFunBegin ();
vfe_ShellPropMatlFun (shellprop,matlfunprop);
/* Create 3D shell element formulation */
shell3d = vfe_Shell3DBegin ();
vfe_Shell3DSetParami (shell3d,VFE_TECH,VFE_TECH_ANS);
vfe_Shell3DSetParami (shell3d,VFE_ROTINERTIA,SYS_ON);
vfe_Shell3DSetTopology (shell3d,SYS_SHAPEQUAD,2,0);
/* Register material functions */
vfe_Shell3DSetObject (shell3d,VFE_MATLFUN,matlfunprop);
/* Form stiffness matrix */
vfe_Shell3DStiff (shell3d,x,kl);
printf("Diagonal of stiffness matrix\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
k = 6*i+j+1;
printf("%12.4e ",kl[k*(k+1)/2-1]);
}
printf("\n");
}
/* Form consistent mass matrix */
vfe_Shell3DMass (shell3d,x,mc);
/* Form diagonal mass matrix */
vfe_Shell3DMassDiag (shell3d,x,md);
printf("Diagonal mass matrix\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",md[6*i+j]);
}
printf("\n");
}
/* Form load vector for constant pressure */
for(i = 0; i < 4; i++) v[i] = 1.;
vfe_Shell3DDistLoad (shell3d,x,SYS_FACE,0,VFE_DISTLOAD_PRES,v,f);
printf("Constant pressure load vector\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",f[6*i+j]);
}
printf("\n");
}
/* Delete objects */
vfe_Shell3DEnd (shell3d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
vfe_MatlFunEnd (matlfunprop);
vfe_ShellPropEnd (shellprop);
return 0;
}
Only a single MatlFun object is required. All material functions are registered using vfe_MatlFunSet. Note that the material functions registered in the MatlFun object provide for a user specified pointer (object) to be passed as the first argument. In this example the pointer argument is used to pass the element thickness, effective shear factor and material modulii so that the shell wall matrix can be computed.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[4][3] = {
{0.,0.,0.}, {1.,0.,0.}, {1.,1.,0.}, {0.,1.,0.} };
typedef struct {
Vdouble t; /* thickness */
Vdouble shearfactor; /* shear factor */
Vdouble E; /* Young's modulus */
Vdouble nu; /* Poisson's ratio */
Vdouble rho; /* density */
} MaterialProperties;
/* function to compute shell wall matrix */
static void
exam_DConstruct(Vobject *obj, Vdouble dmatrix[8][8])
{
MaterialProperties *mat = (MaterialProperties*)obj;
Vdouble tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7;
Vdouble E = mat->E;
Vdouble t = mat->t;
Vdouble nu = mat->nu;
Vdouble shearfactor = mat->shearfactor;
Vint i,j;
tmp1 = E*t*t*t/(12.*(1.-nu*nu));
tmp2 = nu*tmp1;
tmp3 = 0.5*(1-nu)*tmp1;
tmp4 = E*t/(1.-nu*nu);
tmp5 = nu*tmp4;
tmp6 = 0.5*(1.-nu)*tmp4;
tmp7 = shearfactor*E*t/(2.*(1+nu));
for(i = 0; i < 8; ++i) {
for(j = 0; j < 8; ++j) {
dmatrix[i][j] = 0.;
}
}
dmatrix[0][0] = tmp4;
dmatrix[0][1] = tmp5;
dmatrix[1][0] = tmp5;
dmatrix[1][1] = tmp4;
dmatrix[2][2] = tmp6;
dmatrix[3][3] = tmp1;
dmatrix[3][4] = tmp2;
dmatrix[4][3] = tmp2;
dmatrix[4][4] = tmp1;
dmatrix[5][5] = tmp3;
dmatrix[6][6] = tmp7;
dmatrix[7][7] = tmp7;
}
/*----------------------------------------------------------------------
Material functions for laminate description
----------------------------------------------------------------------*/
static void
exam_DCond(Vobject *obj, Vint cond[])
{
/* cond contains all zeros for 3D preintegrated shell */
}
static void
exam_DInfo(Vobject *obj, Vint type, Vint *info)
{
if(type == VFE_MATLELEM) {
*info = SYS_ELEM_SHELL;
} else if(type == VFE_MATLFULL) {
*info = SYS_OFF;
} else if(type == VFE_STRAINTYPE) {
*info = VFE_SMALLSTRAIN;
}
}
static void
exam_DMatrix(Vobject *obj, Vdouble dm[])
{
Vint i,j,k;
Vdouble dmatrix[8][8];
exam_DConstruct (obj,dmatrix);
k = 0;
for(i = 0; i < 8; i++) {
for(j = 0; j <= i; j++) {
dm[k++] = dmatrix[i][j];
}
}
}
static void
exam_Stress(Vobject *obj, Vdouble dfrm[], Vdouble strs[])
{
Vint i, j;
Vdouble dmatrix[8][8];
exam_DConstruct (obj,dmatrix);
for(i = 0; i < 8; i++) {
strs[i] = 0.;
for(j = 0; j < 8; j++) {
strs[i] += dmatrix[i][j] * dfrm[j];
}
}
}
static void
exam_StressDMatrix(Vobject *obj,
Vdouble dfrm[], Vint dflag, Vdouble strs[], Vdouble dm[])
{
exam_Stress (obj,dfrm,strs);
/* assume DMatrix is not a function of deformation */
if(dflag) {
exam_DMatrix (obj,dm);
}
}
static void
exam_Modulus(Vobject *obj, Vdouble *modulus)
{
MaterialProperties *mat = (MaterialProperties*)obj;
Vdouble t = mat->t;
Vdouble dmatrix[8][8];
exam_DConstruct (obj,dmatrix);
*modulus = dmatrix[0][0]*t*t/12.;
}
static void
exam_Density(Vobject *obj, Vdouble density[])
{
MaterialProperties *mat = (MaterialProperties*)obj;
Vdouble t = mat->t;
Vdouble rho = mat->rho;
density[0] = rho*t;
density[1] = 0.;
density[2] = rho*t*t*t/12.;
}
static void
exam_Prop(Vobject *obj, Vint ptype, Vdouble prop)
{
/* only called if vfe_Shell3DSetPropPtr is called */
}
static void
exam_NumHist(Vobject *obj, Vint *numhist)
{
*numhist = 0;
}
static void
exam_ElemIntPnt(Vobject *obj,
Vint elemintpnt, Vdouble r[3], Vdouble w, Vdouble dj)
{
/* use this function to compute any user defined results */
}
static void
exam_NumStressStrain(Vobject *obj, Vint *num)
{
/* use this function to compute layered stress and strain */
*num = 0;
}
/*----------------------------------------------------------------------
Generate Stiffness Matrix for 4 Node Quadrilateral Shell
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, k;
Vdouble kl[300];
Vdouble mc[300], md[24];
Vdouble v[4], f[24];
vfe_Shell3D *shell3d;
vfe_MatlFun *mf;
MaterialProperties mat;
/* shell and material properties */
mat.t = 0.1;
mat.shearfactor = 5./6.;
mat.E = 1.e+7;
mat.nu = 0.3;
mat.rho = 0.0000133;
/* Create MatlFun object and load material functions */
mf = vfe_MatlFunBegin ();
vfe_MatlFunSet (mf,MATLFUN_DCOND ,(void(*)(void))exam_DCond);
vfe_MatlFunSet (mf,MATLFUN_DINFO ,(void(*)(void))exam_DInfo);
vfe_MatlFunSet (mf,MATLFUN_DMATRIX ,(void(*)(void))exam_DMatrix);
vfe_MatlFunSet (mf,MATLFUN_STRESS ,(void(*)(void))exam_Stress);
vfe_MatlFunSet (mf,MATLFUN_STRESSDMATRIX ,(void(*)(void))exam_StressDMatrix);
vfe_MatlFunSet (mf,MATLFUN_MODULUS ,(void(*)(void))exam_Modulus);
vfe_MatlFunSet (mf,MATLFUN_DENSITY ,(void(*)(void))exam_Density);
vfe_MatlFunSet (mf,MATLFUN_PROP ,(void(*)(void))exam_Prop);
vfe_MatlFunSet (mf,MATLFUN_NUMHIST ,(void(*)(void))exam_NumHist);
vfe_MatlFunSet (mf,MATLFUN_ELEMINTPNT ,(void(*)(void))exam_ElemIntPnt);
vfe_MatlFunSet (mf,MATLFUN_NUMSTRESSSTRAIN,(void(*)(void))exam_NumStressStrain);
vfe_MatlFunSetObj(mf,(Vobject*)&mat);
/* Create 3D shell element formulation */
shell3d = vfe_Shell3DBegin ();
vfe_Shell3DSetParami (shell3d,VFE_TECH,VFE_TECH_ANS);
vfe_Shell3DSetParami (shell3d,VFE_ROTINERTIA,SYS_ON);
vfe_Shell3DSetTopology (shell3d,SYS_SHAPEQUAD,2,0);
/* Register material functions */
vfe_Shell3DSetObject (shell3d,VFE_MATLFUN,mf);
/* Form stiffness matrix */
vfe_Shell3DStiff (shell3d,x,kl);
printf("Diagonal of stiffness matrix\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
k = 6*i+j+1;
printf("%12.4e ",kl[k*(k+1)/2-1]);
}
printf("\n");
}
/* Form consistent mass matrix */
vfe_Shell3DMass (shell3d,x,mc);
/* Form diagonal mass matrix */
vfe_Shell3DMassDiag (shell3d,x,md);
printf("Diagonal mass matrix\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",md[6*i+j]);
}
printf("\n");
}
/* Form load vector for constant pressure */
for(i = 0; i < 4; i++) v[i] = 1.;
vfe_Shell3DDistLoad (shell3d,x,SYS_FACE,0,VFE_DISTLOAD_PRES,v,f);
printf("Constant pressure load vector\n");
for(i = 0; i < 4; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",f[6*i+j]);
}
printf("\n");
}
/* Delete objects */
vfe_Shell3DEnd (shell3d);
vfe_MatlFunEnd (mf);
return 0;
}
The procedure for defining the element properties and material properties is similar to Example 3 for shell elements. A BeamProp object must be used to specify the beam section construction. In this case the BeamProp object is configured to manage general section area properties using vfe_BeamPropDef and vfe_BeamPropSetSection. A beam section may only consist of a single material which is specified to be material 1. The MatlFun object associated with the linear material properties is registered with BeamProp as material 1 using vfe_BeamPropSetMatlFun. An additional MatlFun object will now be required to connect the BeamProp object to the Beam3D element formulation object. Therefore another MatlFun object is instanced and filled with the interface to the BeamProp object using vfe_BeamPropMatlFun.
A Beam3D element formulation object is instanced and configured to generate element matrices and vectors for parabolic beams (3 nodes) using an isoparametric technology. A preintegrated formulation (using moments and curvatures) is chosen.
The MatlFun object, matlfunprop, is registered with Beam3D using vfe_Beam3DSetObject. This completes the connections between the element formulation, the element properties and the material model. The area properties for a constant 1. by .2 rectangular beam section are specified at the element nodes. The area, moments of inertia about the y and z section axes and the torsional constant are set using vfe_Beam3DSetPropPtr Any other element nodal properties such as offsets or beam node normals would be specified in a similar manner. Element quantities such as stiffness, mass, etc. can now be computed.
The element stiffness matrix is generated using vfe_Beam3DStiff The lower triangle of the matrix is returned in array kl. All element formulations supported by Beam3D use 3 translations and 3 rotations per node. The total number of degrees of freedom for the 3 node beam is 18. The output arrays for matrices and vectors have been dimensioned exactly to accommodate the output for this element. Both the consistent and diagonal mass matrices are generated. A constant line load is computed using vfe_Beam3DDistLoad.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[4][3] = {
{0.,0.,0.}, {1.,0.,0.}, {2.,0.,0.} };
/*----------------------------------------------------------------------
Generate Stiffness Matrix for 3 Node Beam
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, k;
Vdouble prop[2];
Vdouble bprop[10];
Vdouble kl[171];
Vdouble mc[171], md[18];
Vdouble v[3][3], f[18];
vfe_Beam3D *beam3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
vfe_MatlFun *matlfunprop;
vfe_BeamProp *beamprop;
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 10000000.;
prop[1] = .3;
vfe_LinMatSetElasProp (linmat,prop);
/*
vfe_LinMatSetDensity (linmat,.0000133);
*/
vfe_LinMatSetDensity (linmat,1.33);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create shell property object */
beamprop = vfe_BeamPropBegin ();
vfe_BeamPropDef (beamprop,BEAMPROP_SECTION);
bprop[0] = 2.000000e-02;
bprop[1] = 6.666667e-05;
bprop[2] = 1.666667e-05;
bprop[3] = 0.;
bprop[4] = 4.575848e-05;
bprop[5] = 5./6.;
bprop[6] = 5./6.;
bprop[7] = 0.;
bprop[8] = 0.;
bprop[9] = 0.;
vfe_BeamPropSetSection (beamprop,bprop,1);
vfe_BeamPropSetMatlFun (beamprop,1,matlfun);
/* Create MatlFun object and load shell property functions */
matlfunprop = vfe_MatlFunBegin ();
vfe_BeamPropMatlFun (beamprop,matlfunprop);
/* Create 3D shell element formulation */
beam3d = vfe_Beam3DBegin ();
vfe_Beam3DSetParami (beam3d,VFE_TECH,VFE_TECH_ISOP);
vfe_Beam3DSetParami (beam3d,VFE_ROTINERTIA,SYS_ON);
vfe_Beam3DSetTopology (beam3d,3);
/* Register material functions */
vfe_Beam3DSetObject (beam3d,VFE_MATLFUN,matlfunprop);
/* Form stiffness matrix */
vfe_Beam3DStiff (beam3d,x,kl);
printf("Diagonal of stiffness matrix\n");
for(i = 0; i < 3; i++) {
for(j = 0; j < 6; j++) {
k = 6*i+j+1;
printf("%12.4e ",kl[k*(k+1)/2-1]);
}
printf("\n");
}
/* Form consistent mass matrix */
vfe_Beam3DMass (beam3d,x,mc);
/* Form diagonal mass matrix */
vfe_Beam3DMassDiag (beam3d,x,md);
printf("Diagonal mass matrix\n");
for(i = 0; i < 3; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",md[6*i+j]);
}
printf("\n");
}
/* Form load vector for constant pressure */
for(i = 0; i < 3; i++) {
v[i][0] = 1.;
v[i][1] = 0.;
v[i][2] = 0.;
}
vfe_Beam3DDistLoad (beam3d,x,SYS_EDGE,0,VFE_DISTLOAD_TRAC,(Vdouble*)v,f);
printf("Constant traction load vector\n");
for(i = 0; i < 3; i++) {
for(j = 0; j < 6; j++) {
printf("%12.4e ",f[6*i+j]);
}
printf("\n");
}
/* Delete objects */
vfe_Beam3DEnd (beam3d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
vfe_MatlFunEnd (matlfunprop);
vfe_BeamPropEnd (beamprop);
return 0;
}
In this example the constraint is satisfied using the penalty function technique. The multi-point constraint is represented by a properly constructed stiffness matrix which when assembled into the global stiffness matrix will ensure that the constraint is satisfied approximately. The greater the magnitude of the stiffness coefficients used to satisfy the constraint the more accurate the approximation of the constraint will become. The magnitude of the stiffness coefficients is controlled by the penalty parameter which is specified using vfe_MPCSetParamd. This parameter is generally chosen to be 4 or 5 orders of magnitude larger than the stiffnesses of the elastic elements to which it is attached. The stiffness of constraint is generated for each link using vfe_MPCStiff. The number of degrees of freedom in the link stiffness is returned using vfe_MPCNumDof. In practice, the link stiffness should be assembled into the global system stiffness matrix using the same methodology as any finite element such as a solid, shell or beam.
#include <stdio.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vint itwo[6] = {2,2,2,2,2,2};
static Vdouble x[3][2][3] = {
{ {0.,0.,1.}, {0.,0.,0.} },
{ {0.,0.,1.}, {1.,0.,0.} },
{ {0.,0.,1.}, {0.,1.,0.} } };
/*----------------------------------------------------------------------
Rigid Body MPC Using Penalty Function
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n;
Vint loc[6], tag[6], tagdep[6];
Vdouble kl[12*(12+1)/2];
vfe_MPC *mpc;
Vint nedofs;
/* Create MPC object */
mpc = vfe_MPCBegin ();
/* specify rigid body type constraint */
vfe_MPCDef (mpc,MPC_KINE);
/* set penalty factor */
vfe_MPCSetParamd (mpc,MPC_PENALTY,1.e+12);
/* set RBE info */
loc[0] = 1;
loc[1] = 1;
loc[2] = 1;
loc[3] = 1;
loc[4] = 1;
loc[5] = 1;
tag[0] = SYS_DOF_TX;
tag[1] = SYS_DOF_TY;
tag[2] = SYS_DOF_TZ;
tag[3] = SYS_DOF_RX;
tag[4] = SYS_DOF_RY;
tag[5] = SYS_DOF_RZ;
tagdep[0] = SYS_DOF_TX;
tagdep[1] = SYS_DOF_TY;
tagdep[2] = SYS_DOF_TZ;
tagdep[3] = SYS_DOF_RX;
tagdep[4] = SYS_DOF_RY;
tagdep[5] = SYS_DOF_RZ;
vfe_MPCSetKine (mpc,loc,tag,6,itwo,tagdep);
/* return total number of element dofs */
vfe_MPCNumDof (mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
/* Form stiffness matrix */
for(n = 0; n < 3; n++) {
vfe_MPCStiff (mpc,x[n],kl);
printf("Lower triangle of stiffness matrix, element= %d\n",n+1);
for(i = 0; i < nedofs; ++i ) {
for(j = 0; j <= i; ++j) {
printf("%9.2e ",kl[i*(i+1)/2+j]);
}
printf("\n");
}
}
/* Delete object */
vfe_MPCEnd (mpc);
return 0;
}
A Young's modulus of 206.9 and a Poisson's ratio of 0.29 are set. The exponential law is chosen as the hardening curve, with a hardening modulus of 16.93 and a hardening exponent of 0.265. The yield stress is chosen as 0.45. The PlasMat module requires that pointers be set for the history variables. The number of required history variables is obtained with a call to vfe_MatlFunNumHist. Two history variables are allocated: one to store the converged results, and one for the output of the current solution. The arrays hist1 and hist2 are allocated for this purpose, and are set in the PlasMat object. These values must be internally initialized with a call to vfe_MatlFunInitHist. Users are advised to always use the PlasMat history initialization function, as not all history variables can be initialized to zero.
Two pointers are declared, "curold" and "curnew". These two pointers will facilitate the updating of the internal variables by a simple swap of the arrays pointed to by these pointers. The material is cycled through a strain history that takes two axial strains through the following path: 0 -> ebar -> -ebar -> ebar, where "ebar" is 10 times the elastic strain as the stress reaches the yield value, or ebar = 10*yield/Elastic modulus
An important point is the fact that as the call to vfe_MatlFunStressDMatrix is made and the results converged, the converged history variables must be set as the "old" history variables. This is done via the pointer swap mentioned above.
#include <stdlib.h>
#include <stdio.h>
#include "vfe/vfe.h"
/*----------------------------------------------------------------------
Infinitesimal elasto-plasticity under bi-axial stretching
----------------------------------------------------------------------*/
int
main()
{
Vint m;
vfe_PlasMat *plasmat;
vfe_MatlFun *matlfun;
Vdouble elasprop[2];
Vdouble strain[6], stress[6];
Vdouble yield, hardmod, hardexp, ebar;
Vdouble *hist1, *hist2, *curold, *curnew, *tmp;
Vint nhist;
/* Create material object and set properties */
plasmat = vfe_PlasMatBegin ();
matlfun = vfe_MatlFunBegin ();
vfe_PlasMatMatlFun (plasmat,matlfun);
vfe_PlasMatSetParami (plasmat,VFE_STRAINTYPE,VFE_SMALLSTRAIN);
elasprop[0] = 206.9;
elasprop[1] = 0.29;
vfe_PlasMatSetElasProp (plasmat,elasprop);
vfe_PlasMatDef (plasmat,PLASMAT_EXPLAW);
yield = 0.45;
hardmod = 16.93;
hardexp = 0.265;
vfe_PlasMatSetExponent (plasmat,yield,hardmod,hardexp);
/* Query number of history variables */
vfe_MatlFunNumHist (matlfun,&nhist);
/* Allocate history variables for single point */
hist1 = (Vdouble*)malloc(nhist*sizeof(Vdouble));
hist2 = (Vdouble*)malloc(nhist*sizeof(Vdouble));
/* Set pointer to the current old and new histories */
curold = hist1;
curnew = hist2;
/* Set history pointers and initialize history */
vfe_MatlFunHistPtr (matlfun,curold,curnew);
vfe_MatlFunInitHist (matlfun);
/* Set strain limits */
ebar = 10.*yield/elasprop[0];
/* Loop over increments (including initialization) */
printf(" Strain Stress\n");
for(m = 0; m <= 100; m++) {
/* Compute uniform strain for this load increment */
if(m <= 20) {
strain[0] = m*ebar/20.;
strain[1] = strain[0];
} else if(m > 20 && m <= 60) {
strain[0] = ebar - (m-20)*ebar/20.;
strain[1] = strain[0];
} else if(m > 60 && m <= 100) {
strain[0] = -ebar + (m-60)*ebar/20.;
strain[1] = strain[0];
}
strain[2] = strain[3] = strain[4] = strain[5] = 0.;
/* Compute stress */
vfe_MatlFunStress (matlfun,strain,stress);
printf("%12.4e %12.4e\n",strain[0],stress[0]);
/* Swap history */
tmp = curnew;
curnew = curold;
curold = tmp;
vfe_MatlFunHistPtr (matlfun,curold,curnew);
}
/* free history variables */
free(hist1);
free(hist2);
/* Delete object */
vfe_PlasMatEnd (plasmat);
vfe_MatlFunEnd (matlfun);
return 0;
}
The deformation is defined as an array of 9 numbers, with entries given in the following order: F_xx, F_yx, F_zx, F_xy, F_yy, F_zy, F_xz, F_yz, F_zz. In the example, the entries F_yx and F_xy are set to zero, which will result in no shear stress in the plane of the plate/shell.
The condensation flags indicate that only the zz stress component is to be eliminated. The transverse shear components are assumed present, and are treated elastically, assuming a linear relationship between the transverse shear components and its corresponding infinitesimal strains. A shear constant G is internally computed for this purpose.
The example shows how the Mooney-Rivlin material - and other hyperelastic materials - are generalized. Users are required to enter a bulk modulus to enforce the incompressibility condition. A large value of the bulk modulus enforces this condition in a penalized form.
#include <stdio.h>
#include "vfe/vfe.h"
/* define deformation gradient matrix */
static Vdouble F[9] = { 2.900, 0.000, 0.029,
0.000, 1.984, 0.029,
0.029, 0.029, 1.000 };
/*----------------------------------------------------------------------
Mooney-Rivlin material in plane stress
----------------------------------------------------------------------*/
int
main()
{
Vint i, j, n;
vfe_HyperMat *hypermat;
vfe_MatlFun *matlfun;
Vdouble bulk, mooney[2], stress[6], dm[21];
Vint cond[6] = {0,0,1,0,0,0};
/* Create material object and set properties */
hypermat = vfe_HyperMatBegin ();
matlfun = vfe_MatlFunBegin ();
vfe_HyperMatMatlFun (hypermat,matlfun);
vfe_HyperMatDef (hypermat,HYPERMAT_MOONEY);
bulk = 7.5e+8;
mooney[0] = 37500.;
mooney[1] = 375.;
vfe_HyperMatSetMooney (hypermat,bulk,mooney);
/* Set condensation flag for plane stress */
vfe_MatlFunDCond (matlfun,cond);
/* Compute Cauchy stress and tangent matrix */
vfe_MatlFunStressDMatrix (matlfun,F,SYS_ON,stress,dm);
printf("Deformation Gradient:\n");
for(i = 0; i < 3; i++) {
for(j = 0; j < 3; j++) {
printf("%12.4e ",F[3*j+i]);
}
printf("\n");
}
printf("\n");
printf("Cauchy Stress:\n");
for(i = 0; i < 6; i++) {
printf("%12.4e\n",stress[i]);
}
printf("\n");
printf("Tangent Matrix:\n");
for(n = 0, j = 0; j < 6; j++) {
for(i = 0; i <= j; i++) {
printf("%12.4e ",dm[n++]);
}
printf("\n");
}
/* Delete object */
vfe_HyperMatEnd (hypermat);
vfe_MatlFunEnd (matlfun);
return 0;
}
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "base/base.h"
#include "vfe/vfe.h"
static Vdouble x[8][3] = {
{0.,0.,0.}, {1.,0.,0.}, {1.,1.,0.}, {0.,1.,0.},
{0.,0.,1.}, {1.,0.,1.}, {1.,1.,1.}, {0.,1.,1.} };
/*----------------------------------------------------------------------
Corotation of a Solid3D element
----------------------------------------------------------------------*/
int
main()
{
Vint i,j,n;
vfe_Solid3D *solid3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
vfe_Corot *corot;
Vdouble prop[2];
Vdouble x1[8][3], x2[8][3];
Vint nedofs, loc[24], tag[24];
Vdouble u[24];
Vdouble tm[3][3];
Vdouble rf[24],ke[300];
Vdouble rot[3][3], angle;
/* Create material object and set properties */
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
prop[0] = 1.;
prop[1] = 0.;
vfe_LinMatSetElasProp (linmat,prop);
/* Create MatlFun object and load material functions */
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat, matlfun);
/* Create 3D solid element formulation */
solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetParami (solid3d,VFE_TECH,VFE_TECH_MIXED);
vfe_Solid3DSetTopology (solid3d,SYS_SHAPEHEX,2,0,0);
/* Register material functions */
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
/* Determine number of degrees of freedom */
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
/* Create corotational object */
corot = vfe_CorotBegin ();
vfe_CorotSetTopology (corot,SYS_SHAPEHEX,2,0,0);
vfe_CorotSetMap(corot,nedofs,loc,tag);
/* Shear in x-y plane */
for(n = 0; n < 8; n++) {
x1[n][0] = x[n][0] + 0.1*x[n][1];
x1[n][1] = x[n][1];
x1[n][2] = x[n][2];
}
/* Add a 60-degree rotation in x-y plane */
angle = acos(-1.)/3;
rot[0][0] = cos(angle);
rot[0][1] = -sin(angle);
rot[0][2] = 0.;
rot[1][1] = cos(angle);
rot[1][0] = sin(angle);
rot[1][2] = 0.;
rot[2][0] = 0.;
rot[2][1] = 0.;
rot[2][2] = 1.;
for(n = 0; n < 8; n++) {
for(i = 0; i < 3; i++) {
x2[n][i] = 0.;
for(j = 0; j < 3; j++) {
x2[n][i] += rot[i][j]*x1[n][j];
}
}
}
/* Compute element displacements */
for(n = 0; n < 8; n++) {
for(i = 0; i < 3; i++) {
u[3*n + i] = x2[n][i] - x[n][i];
}
}
/* Pull-back deformation to reference frame */
vfe_CorotPull (corot,x,u,tm);
printf("\nNodal displacements in corotated frame:\n");
printf("Node ux uy uz\n");
for(n = 0; n < 8; n++) {
printf("%4d %15.8f %15.8f %15.8f\n",n+1,u[3*n],u[3*n+1],u[3*n+2]);
}
/* Form reaction force and stiffness in reference frame */
vfe_Solid3DReactStiff (solid3d,x,u,SYS_ON,rf,ke);
printf("\nReaction forces in inertial frame:\n");
printf("Node fx fy fz\n");
for(n = 0; n < 8; n++) {
printf("%4d %15.8f %15.8f %15.8f\n",n+1,rf[3*n],rf[3*n+1],rf[3*n+2]);
}
/* Push reaction force and stiffness to corotated frame */
vfe_CorotPush (corot,SYS_ON,rf,SYS_ON,ke);
printf("\nReaction forces in corotated frame:\n");
printf("Node fx fy fz\n");
for(n = 0; n < 8; n++) {
printf("%4d %15.8f %15.8f %15.8f\n",n+1,rf[3*n],rf[3*n+1],rf[3*n+2]);
}
/* Delete objects */
vfe_Solid3DEnd (solid3d);
vfe_CorotEnd (corot);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
return 0;
}
#include "base/base.h"
#include "vfe/vfe.h"
/*----------------------------------------------------------------------
Rigid Body Mode Detection using NodeGeom
----------------------------------------------------------------------*/
int
main()
{
Vint i;
vfe_NodeGeom *nodegeom;
Vdouble x[2][3];
Vint ntags, loc[4], tag[4];
Vint ntra, nrot;
Vdouble tra[3][3], rot[3][6], off[3][3];
Vint mode;
/* Create NodeGeom object */
nodegeom = vfe_NodeGeomBegin ();
/* two node coordintes at which 4 constraints applied */
x[0][0] = 1.;
x[0][1] = 0.;
x[0][2] = 0.;
x[1][0] = 1.;
x[1][1] = 1.;
x[1][2] = 0.;
/* Ty and Rx at first node, Tx and Ty at second node */
ntags = 4;
loc[0] = 1;
loc[1] = 1;
loc[2] = 2;
loc[3] = 2;
tag[0] = SYS_DOF_TY;
tag[1] = SYS_DOF_RX;
tag[2] = SYS_DOF_TX;
tag[3] = SYS_DOF_TY;
/* set center to first node as an example */
vfe_NodeGeomSetCenter (nodegeom,x[0]);
/* compute unconstrained rigid body modes */
vfe_NodeGeomFreeMode (nodegeom,x,0,NULL,ntags,loc,tag,
&ntra,tra,&nrot,rot,off);
mode = 0;
/* print unconstrained modes */
printf("Translational modes\n");
for(i = 0; i < ntra; i++) {
mode += 1;
printf("mode %d, vector= %e %e %e\n",mode,
tra[i][0],tra[i][1],tra[i][2]);
}
printf("Rotational modes\n");
for(i = 0; i < nrot; i++) {
mode += 1;
printf("mode %d, vector= %e %e %e %e %e %e\n",mode,
rot[i][0],rot[i][1],rot[i][2], rot[i][3],rot[i][4],rot[i][5]);
printf(" offset= %e %e %e\n",
off[i][0],off[i][1],off[i][2]);
}
/* Delete object */
vfe_NodeGeomEnd (nodegeom);
return 0;
}
The solid element stiffness and load computations are performed using the Solid3D module. It is initialized for a 10 node tetrahedral element formulation. The linear elastic material is modelled using a LinMat module. The MatlFun module provides the abstract material function interface to be registered with the Solid3D module.
The global degree of freedom table is generated by looping through all elements and declaring their degree of freedom usage with a DofTab object. All permanent restraints are then defined and the global degree of freedom table may be generated using vfs_DofTabProcess. The system stiffness matrix, applied force vector and displacement solution vectors are represented by a SysMatrix and two SysVector objects respectively. Three optional solution strategies are implemented, in-core direct sparse solver, out of core direct sparse solver and iterative solver. In practice one of these three technologies would be selected considering element type, problem size and computer resources, in this example the in-core sparse solver is selected. Note that if the iterative solver is selected using the "Fast" preconditioner, then the DofTab object must be given the element topology and node coordinate information as shown.
The structure of the assembled stiffness matrix is determined by the pre-assembly process and then the actual element stiffness matrix computation and assembly can occur. The function vfs_SysMatrixProcess must be called to perform any processing required before factorization such as degree of freedom reordering to optimize a sparse matrix factorization. The function vfs_SysMatrixFactor is then called to factorize the matrix. If a direct solver is used, a complete factorization is performed, if the iterative solver is used, a preconditioner is computed.
The equivalent nodal loads due to pressure are computed and assembled into the force SysVector. The displacement under the force is computed using vfs_SysMatrixSolve and the nodal displacements are printed. The applied pressure places the block under uniform compression stress equal in magnitude to the applied pressure. Finally the stress and strain in each element is computed and printed.
#include <stdlib.h>
#include "base/basedefs.h"
#include "vfe/vfe.h"
#include "vfs/vfs.h"
/* node coordinates for 27 nodes */
static Vdouble xg[27][3] = { {0.0,0.0,0.0}, {0.5,0.0,0.0}, {1.0,0.0,0.0},
{0.0,1.0,0.0}, {0.5,1.0,0.0}, {1.0,1.0,0.0},
{0.0,2.0,0.0}, {0.5,2.0,0.0}, {1.0,2.0,0.0},
{0.0,0.0,1.5}, {0.5,0.0,1.5}, {1.0,0.0,1.5},
{0.0,1.0,1.5}, {0.5,1.0,1.5}, {1.0,1.0,1.5},
{0.0,2.0,1.5}, {0.5,2.0,1.5}, {1.0,2.0,1.5},
{0.0,0.0,3.0}, {0.5,0.0,3.0}, {1.0,0.0,3.0},
{0.0,1.0,3.0}, {0.5,1.0,3.0}, {1.0,1.0,3.0},
{0.0,2.0,3.0}, {0.5,2.0,3.0}, {1.0,2.0,3.0} };
/* element connectivity for 6 10-node tetrahedra */
static Vint ix[6][10] = { { 1, 9, 7,25, 5, 8, 4,13,17,16},
{ 1, 9,25,19, 5,17,13,10,14,22},
{ 9,25,19,27,17,22,14,18,26,23},
{ 1, 3, 9,19, 2, 6, 5,10,11,14},
{ 3, 9,19,27, 6,14,11,15,18,23},
{ 3,27,19,21,15,23,11,12,24,20} };
/* permanent node restraints in x, y and z */
static Vint xbc[9] = { 1, 4, 7,10,13,16,19,22,25};
static Vint ybc[9] = { 1, 2, 3,10,11,12,19,20,21};
static Vint zbc[9] = { 1, 2, 3, 4, 5, 6, 7, 8, 9};
/* pressure loads on second face of elements 5 and 6 */
static Vint lface[2][2] = { {5,2}, {6,2} };
static Vdouble pres[10] = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
/* material properties */
static Vdouble eprop[2] = { 1.e+6, 0.3 };
/* nodal degree of freedom types */
static Vint doftags[3] = {SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ};
/*----------------------------------------------------------------------
Linear Stress Analysis Using 10 Node Tetrahedra
----------------------------------------------------------------------*/
int
main()
{
/* VfeTools objects */
vfe_Solid3D *solid3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
/* VfsTools objects */
vfs_DofTab *doftab;
vfs_SysMatrix *stiff;
vfs_SysVector *disp, *force;
Vint i, k, m, n, no;
Vint nloadfaces, nxbc, nybc, nzbc;
Vint nix, numel, numnp;
Vint method;
Vint neq, nre;
Vint nedofs, loc[VFE_MAXDOF], tag[VFE_MAXDOF], lm[VFE_MAXDOF];
Vdouble ke[VFE_MAXDOF*(VFE_MAXDOF+1)/2], f[VFE_MAXDOF], u[VFE_MAXDOF];
Vdouble x[VFE_MAXNODE][3];
Vdouble strs[VFE_MAXNODE*6], strn[VFE_MAXNODE*6];
/* set problem parameters */
nloadfaces = 2;
nxbc = 9;
nybc = 9;
nzbc = 9;
nix = 10;
numel = 6;
numnp = 27;
/* element initializations */
solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetTopology (solid3d,SYS_SHAPETET,3,0,0);
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
/* material initializations */
matlfun = vfe_MatlFunBegin ();
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vfe_LinMatSetElasProp (linmat,eprop);
vfe_LinMatMatlFun (linmat,matlfun);
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
/* select solver method */
method = 1;
/* declare element degree of freedom usage */
doftab = vfs_DofTabBegin ();
vfs_DofTabDef (doftab,numnp,numel);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
for(m = 1; m <= numel; m++) {
vfs_DofTabSetTopology (doftab,m,SYS_SHAPETET,3,0,0);
vfs_DofTabElemUse (doftab,m,nix,ix[m-1]);
}
if(method == 3) {
vfs_DofTabSetParami (doftab,DOFTAB_REDUCE,DOFTAB_REDUCE_ITER);
for(n = 1; n <= numnp; n++) {
vfs_DofTabSetCoords (doftab,n,xg[n-1]);
}
}
/* set permanent restraints */
for(i = 0; i < nxbc; i++) {
vfs_DofTabNodePerm (doftab,xbc[i],1,&doftags[0]);
}
for(i = 0; i < nybc; i++) {
vfs_DofTabNodePerm (doftab,ybc[i],1,&doftags[1]);
}
for(i = 0; i < nzbc; i++) {
vfs_DofTabNodePerm (doftab,zbc[i],1,&doftags[2]);
}
vfs_DofTabProcess (doftab,&neq,&nre);
printf("Number of equations: %d\n",neq);
printf("Number of reactions: %d\n",nre);
/* define system matrix and vectors */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
disp = vfs_SysVectorBegin ();
vfs_SysVectorDef (disp,neq);
force = vfs_SysVectorBegin ();
vfs_SysVectorDef (force,neq);
/* in-core sparse solver */
if(method == 1) {
printf("In core sparse solver\n");
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetParami (stiff,SYSMATRIX_SOLVERTYPE,
SYSMATRIX_SOLVERTYPE_MFP);
/* out-of-core sparse solver */
} else if(method == 2) {
printf("Out of core sparse solver\n");
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetParami (stiff,SYSMATRIX_TECH,SYSMATRIX_TECH_SPARSE);
vfs_SysMatrixSetParami (stiff,SYSMATRIX_SOLVERTYPE,
SYSMATRIX_SOLVERTYPE_MFP);
vfs_SysMatrixSetParami (stiff,SYSMATRIX_FACTOR_OOC,SYS_ON);
vfs_SysMatrixSetOocFile (stiff,SYSMATRIX_KFACTOR_FILE,"outofcore.ooc");
/* iterative solver */
} else if(method == 3) {
printf("Iterative solver\n");
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_ITER,neq);
vfs_SysMatrixSetParami (stiff,SYSMATRIX_SOLVERTYPE,
SYSMATRIX_SOLVERTYPE_MFP);
}
/* preprocess stiffness matrix */
vfs_SysMatrixPreProcess (stiff);
vfs_SysMatrixInit (stiff,0.);
/* assemble stiffness matrix */
for(m = 1; m <= numel; m++) {
for(i = 0; i < nix; i++) {
n = ix[m-1][i];
x[i][0] = xg[n-1][0];
x[i][1] = xg[n-1][1];
x[i][2] = xg[n-1][2];
}
vfe_Solid3DStiff (solid3d,x,ke);
vfs_DofTabElemDof (doftab,m,nix,ix[m-1],&nedofs,lm);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
/* factor stiffness matrix */
vfs_SysMatrixProcess (stiff);
vfs_SysMatrixFactor (stiff);
/* form distributed load on element faces and assemble */
vfs_SysVectorInit (force,0.0);
for(k = 0; k < nloadfaces; k++) {
m = lface[k][0];
no = lface[k][1];
for(i = 0; i < nix; i++) {
n = ix[m-1][i];
x[i][0] = xg[n-1][0];
x[i][1] = xg[n-1][1];
x[i][2] = xg[n-1][2];
}
vfe_Solid3DDistLoad (solid3d,x,SYS_FACE,no,VFE_DISTLOAD_PRES,pres,f);
vfs_DofTabElemDof (doftab,m,nix,ix[m-1],&nedofs,lm);
vfs_SysVectorAssem (force,nedofs,lm,f);
}
/* solve for displacements */
vfs_SysMatrixSolve (stiff,force,disp);
/* print displacements */
printf("Displacements:\n");
printf(" Node ux uy uz\n");
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,3,doftags,lm);
for(i = 0; i < 3; i++) u[i] = 0.;
vfs_SysVectorGather (disp,3,lm,u);
printf(" %4d %8.1e %8.1e %8.1e\n",n,u[0],u[1],u[2]);
}
printf("\n");
/* compute and print stresses and strains */
for(m = 1; m <= numel; m++) {
for(i = 0; i < nix; i++) {
n = ix[m-1][i];
x[i][0] = xg[n-1][0];
x[i][1] = xg[n-1][1];
x[i][2] = xg[n-1][2];
}
vfs_DofTabElemDof (doftab,m,nix,ix[m-1],&nedofs,lm);
for(i = 0; i < nedofs; i++) u[i] = 0.;
vfs_SysVectorGather (disp,nedofs,lm,u);
vfe_Solid3DStrsStrn (solid3d,x,u,strs,strn);
printf("Element %2d\n\n",m);
printf(" Node exx eyy ezz exy eyz ezx\n");
for(i = 0; i < nix; i++) {
printf("%6d %8.1e %8.1e %8.1e %8.1e %8.1e %8.1e\n",
i+1,strn[6*i ],strn[6*i+1],strn[6*i+2],
strn[6*i+3],strn[6*i+4],strn[6*i+5]);
}
printf("\n");
printf(" Node sxx syy szz sxy syz szx\n");
for(i = 0; i < nix; i++) {
printf("%6d %8.1e %8.1e %8.1e %8.1e %8.1e %8.1e\n",
i+1,strs[6*i ],strs[6*i+1],strs[6*i+2],
strs[6*i+3],strs[6*i+4],strs[6*i+5]);
}
printf("\n");
}
/* delete objects */
vfe_Solid3DEnd (solid3d);
vfe_MatlFunEnd (matlfun);
vfe_LinMatEnd (linmat);
vfs_SysVectorEnd (disp);
vfs_SysVectorEnd (force);
vfs_SysMatrixEnd (stiff);
vfs_DofTabEnd (doftab);
return 0;
}
An initial .vdm file containing a finite element model may be generated using VisTools example vis/exam/exam49vdm.c. The resulting .vdm file may be visualized using VisTools example vis/exam/exam30edev.c.
#include <stdlib.h>
#include "base/base.h"
#include "vfe/vfe.h"
#include "vfs/vfs.h"
#include "vis/vis.h"
#include "vdm/vdm.h"
/* useful macro for parsing group activity flags */
#define VIEW_FLAG(flags,ind) (((flags) >> ((ind)-1)) & 1)
static void
LoadMatl(vsy_HashTable *matlfunh, vis_MProp *mprop, Vint mid);
static void
SolidProp(Vint m, vis_Connect *connect, vfe_Solid3D *solid3d,
vsy_HashTable *matlfunh, vsy_HashTable *eproph,
Vint *nix, Vint ix[], Vdouble x[][3], vfe_MatlFun **matlfun);
static void
SolidTopo(Vint m, vis_Connect *connect, vfe_Solid3D *solid3d,
Vint *nedofs, Vint loc[], Vint tag[]);
/* nodal degree of freedom types */
static Vint doftags[3] = {SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ};
/*----------------------------------------------------------------------
Linear Stress Analysis Using VisTools and VdmTools
----------------------------------------------------------------------*/
int
main(int argc, char **argv)
{
/* VfeTools objects */
vfe_Solid3D *solid3d;
vfe_MatlFun *matlfun;
vfe_LinMat *linmat;
/* VfsTools objects */
vfs_DofTab *doftab;
vfs_SysMatrix *stiff;
vfs_SysVector *disp, *forc;
/* VisTools objects */
vis_Model *model;
vis_Connect *connect;
vis_GridFun *gridfun;
vis_MProp *mprop;
vis_SProp *sprop;
vis_RCase *rcase;
vis_LCase *lcase;
vis_State *state;
vis_Group *group;
vis_RProp *rprop;
/* VdmTools objects */
vdm_DataFun *datafun;
vdm_NatLib *natlib;
vdm_LMan *lman;
/* Base object */
vsy_HashTable *mproph;
vsy_HashTable *eproph;
vsy_HashTable *matlfunh;
vsy_HashTable *rcaseh;
vsy_HashTable *lcaseh;
vsy_List *spropl;
Vchar path[256];
Vdouble x[VFE_MAXNODE][3];
Vdouble ke[VFE_MAXDOF*(VFE_MAXDOF+1)/2];
Vdouble r[VFE_MAXDOF], f[3*VFE_MAXNODE], u[VFE_MAXDOF];
Vdouble strs[6*VFE_MAXNODE], strn[6*VFE_MAXNODE];
Vdouble etemp[VFE_MAXNODE], *temp;
Vdouble value;
Vint i, j, l, m, n;
Vint numel, numnp;
Vint neq, nre;
Vint nfaces, no, flags, ntypes, itype[4];
Vint nedofs, loc[VFE_MAXDOF], tag[VFE_MAXDOF], lm[VFE_MAXDOF];
Vint nix, ix[VFE_MAXNODE];
Vint featype;
Vint shape, maxi, maxj, maxk;
Vint ntags, rtag[3], type, master;
Vint mid, sid;
Vint rcaseidold, rcaseid;
Vint maxlcaseids, *lcaseids, nlcases;
Vint itrac, tstrain;
Vint ierr;
/* check for proper number of arguments */
if(argc < 2) {
fprintf(stderr,"Usage: %s input/outputfile\n",argv[0]);
fprintf(stderr," input/outputfile is blank, exam10adev.vdm is assumed\n");
strcpy(path,"exam10adev.vdm");
} else {
strcpy(path,argv[1]);
}
/* open native input file */
datafun = vdm_DataFunBegin ();
natlib = vdm_NatLibBegin ();
vdm_NatLibDataFun (natlib,datafun);
vdm_DataFunSetStatus (datafun,VDM_STATUS_ADD);
vdm_DataFunSetConvention (datafun,VDM_CONVENTION_DOUBLE);
vdm_DataFunOpen (datafun,0,path,VDM_NATIVE);
ierr = vdm_DataFunError(datafun);
if(ierr) {
fprintf(stderr,"Error: opening file %s\n",path);
exit(1);
}
/* create lman as helper object */
lman = vdm_LManBegin ();
vdm_LManSetObject (lman,VDM_DATAFUN,datafun);
/* load model */
model = vis_ModelBegin ();
vdm_LManLoadModel (lman,model);
/* extract Connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
printf("Number of nodes= %10d\n",numnp);
printf("Number of elements= %10d\n",numel);
/* check for 3D solid elements only */
for(m = 1; m <= numel; m++) {
vis_ConnectElemAssoc (connect,VIS_FEATYPE,1,&m,&featype);
vis_ConnectTopology (connect,m,&shape,&maxi,&maxj,&maxk);
if(featype != SYS_ELEM_SOLID ||
(shape != SYS_SHAPETET &&
shape != SYS_SHAPEWED &&
shape != SYS_SHAPEHEX)) {
printf("Error: Unsupported element type, exiting\n");
exit(0);
}
}
/* element initializations */
solid3d = vfe_Solid3DBegin ();
for(i = 0; i < VFE_MAXNODE; i++) {
etemp[i] = 0.;
}
vfe_Solid3DSetPropPtr (solid3d,VFE_PROP_TEMPERATURE,etemp);
/* material initializations */
vis_ModelGetHashTable (model,VIS_MPROP,&mproph);
matlfunh = vsy_HashTableBegin ();
vsy_HashTableInitIter (mproph);
while(vsy_HashTableNextIter(mproph,&mid,(Vobject**)&mprop), mprop) {
LoadMatl (matlfunh,mprop,mid);
}
vis_ModelGetHashTable (model,VIS_EPROP,&eproph);
/* instance VfsTools objects */
doftab = vfs_DofTabBegin ();
stiff = vfs_SysMatrixBegin ();
disp = vfs_SysVectorBegin ();
forc = vfs_SysVectorBegin ();
/* create auxiliary group for distributed loads */
group = vis_GroupBegin ();
/* create output state */
state = vis_StateBegin ();
gridfun = vis_GridFunBegin ();
vis_ConnectGridFun (connect,gridfun);
vis_StateSetObject (state,VIS_GRIDFUN,gridfun);
rprop = vis_RPropBegin ();
/* allocate temperature array */
temp = (Vdouble*)malloc(numnp*sizeof(Vdouble));
/* extract load case information */
vis_ModelGetHashTable (model,VIS_LCASE,&lcaseh);
vsy_HashTableCount (lcaseh,&maxlcaseids);
lcaseids = (Vint*)malloc(maxlcaseids*sizeof(Vint));
/* extract restraint case information */
vis_ModelGetHashTable (model,VIS_RCASE,&rcaseh);
rcaseidold = -1;
/* loop over solution properties */
vis_ModelGetList (model,VIS_SPROP,&spropl);
vsy_ListInitIter (spropl);
while(vsy_ListNextIter(spropl,&sid,(Vobject**)&sprop),sprop) {
/* generate degree of freedom table */
vis_SPropValueInteger (sprop,SPROP_RCASE,&rcaseid);
/* check whether restraint set has changed */
if(rcaseid != rcaseidold) {
vfs_DofTabDef (doftab,numnp,numel);
for(m = 1; m <= numel; m++) {
vis_ConnectElemNode (connect,m,&nix,ix);
SolidTopo (m,connect,solid3d,&nedofs,loc,tag);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,m,nix,ix);
}
/* process zero boundary conditions */
if(rcaseid) {
vsy_HashTableLookup (rcaseh,rcaseid,(Vobject**)&rcase);
for(n = 1; n <= numnp; n++) {
vis_RCaseSPCTag (rcase,n,&ntags,rtag);
if(ntags) {
for(i = 0; i < ntags; i++) {
vis_RCaseSPCdv (rcase,n,rtag[i],&type,&value,&master);
if(type == RCASE_FIXED) {
vfs_DofTabNodePerm (doftab,n,1,&rtag[i]);
}
}
}
}
}
vfs_DofTabProcess (doftab,&neq,&nre);
printf("Number of equations= %10d\n",neq);
printf("Number of reactions= %10d\n",nre);
/* define sparse system matrix and vectors */
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysVectorDef (disp,neq);
vfs_SysVectorDef (forc,neq);
/* pre-process stiffness matrix */
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (stiff);
/* initialize matrix */
vfs_SysMatrixInit (stiff,0.);
/* form and assemble stiffness matrix */
printf("Form and assemble stiffness matrix\n");
for(m = 1; m <= numel; m++) {
SolidProp (m,connect,solid3d,matlfunh,eproph,&nix,ix,x,&matlfun);
SolidTopo (m,connect,solid3d,&nedofs,loc,tag);
vfs_DofTabElemDof (doftab,m,nix,ix,&nedofs,lm);
vfe_Solid3DStiff (solid3d,x,ke);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
/* process stiffness matrix */
printf("Process stiffness matrix\n");
vfs_SysMatrixProcess (stiff);
/* factor stiffness */
printf("Factor stiffness matrix\n");
vfs_SysMatrixFactor (stiff);
}
rcaseidold = rcaseid;
/* loop over load cases */
printf("Generating load for solution sequence= %3d\n",sid);
vis_SPropValueInteger (sprop,SPROP_LCASE_NUM,&nlcases);
vis_SPropValueInteger (sprop,SPROP_LCASE,lcaseids);
/* loop over load cases to be added */
vfs_SysVectorInit (forc,0.);
itrac = 0;
tstrain = 0;
vis_SPropValueInteger (sprop,SPROP_THERMALSTRAIN,&tstrain);
for(l = 0; l < nlcases; l++) {
/* extract load case */
vsy_HashTableLookup (lcaseh,lcaseids[l],(Vobject**)&lcase);
/* face based distributed loads */
vis_GroupDef (group,numel,SYS_ELEM,SYS_FACE);
vis_LCaseFaceGroup (lcase,NULL,group);
vis_GroupInitIndex (group);
while(vis_GroupNextIndex (group,&m,&flags),m) {
SolidProp (m,connect,solid3d,matlfunh,eproph,&nix,ix,x,&matlfun);
SolidTopo (m,connect,solid3d,&nedofs,loc,tag);
vfs_DofTabElemDof (doftab,m,nix,ix,&nedofs,lm);
vis_ConnectElemNum (connect,SYS_FACE,m,&nfaces);
for(no = 1; no <= nfaces; no++) {
if(VIEW_FLAG(flags,no) == 0) continue;
vis_LCaseDistType (lcase,SYS_FACE,m,no,&ntypes,itype);
for(i = 0; i < ntypes; i++) {
if(itype[i] == LCASE_TRAC) {
itrac = 1;
for(j = 0; j < 3*VFE_MAXNODE; j++) {
f[j] = 0.;
}
vis_LCaseDistdv (lcase,SYS_FACE,m,no,itype[i],1,f);
vfe_Solid3DDistLoad (solid3d,x,SYS_FACE,no,itype[i],f,r);
vfs_SysVectorAssem (forc,nedofs,lm,r);
}
}
}
}
/* node based thermal loads */
if(tstrain) {
for(i = 0; i < numnp; i++) {
temp[i] = 0.;
}
vis_GroupDef (group,numnp,SYS_NODE,SYS_NONE);
vis_LCaseNodeGroup (lcase,NULL,group);
vis_GroupInitIndex (group);
while (vis_GroupNextIndex (group,&n,&flags), n) {
vis_LCaseConcType (lcase,n,&ntypes,itype);
for(i = 0; i < ntypes; i++) {
if(itype[i] == LCASE_TEMP) {
vis_LCaseConcdv (lcase,n,LCASE_TEMP,&temp[n-1]);
}
}
}
for(m = 1; m <= numel; m++) {
SolidProp (m,connect,solid3d,matlfunh,eproph,&nix,ix,x,&matlfun);
vfe_MatlFunGetObj (matlfun,(Vobject**)&linmat);
vfe_LinMatSetParami (linmat,VFE_THERMALSTRAIN,SYS_ON);
SolidTopo (m,connect,solid3d,&nedofs,loc,tag);
for(i = 0; i < nix; i++) {
etemp[i] = temp[ix[i]-1];
}
for(i = 0; i < nedofs; i++) {
u[i] = 0.;
}
vfs_DofTabElemDof (doftab,m,nix,ix,&nedofs,lm);
vfe_Solid3DReact (solid3d,x,u,r);
vfs_SysVectorAssem (forc,nedofs,lm,r);
}
}
}
if(itrac) {
printf("Distributed traction detected\n");
}
if(tstrain) {
printf("Thermal strains detected\n");
}
/* solve for displacements */
printf("Solve system\n");
vfs_SysMatrixSolve (stiff,forc,disp);
/* save displacements */
printf("Save results\n");
vis_StateDef (state,numnp,SYS_NODE,SYS_NONE,SYS_VECTOR);
for(n = 1; n <= numnp; n++) {
u[0] = u[1] = u[2] = 0.;
vfs_DofTabNodeDof (doftab,n,3,doftags,lm);
vfs_SysVectorGather (disp,3,lm,u);
vis_StateSetDatadv (state,n,u);
}
vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
vis_RPropSetValuec (rprop,RPROP_TITLE,"exam10adev");
vis_RPropSetType (rprop,SYS_RES_D);
vis_RPropSetIds (rprop,sid,1,0);
vdm_LManSaveState (lman,state,rprop);
/* save stresses */
vis_StateDef (state,numel,SYS_ELEM,SYS_NODE,SYS_TENSOR);
vis_StateSetObject (state,VIS_GRIDFUN,gridfun);
vis_RPropDef (rprop,SYS_ELEM,SYS_NODE);
vis_RPropSetValuec (rprop,RPROP_TITLE,"exam10adev");
vis_RPropSetType (rprop,SYS_RES_S);
vis_RPropSetIds (rprop,sid,1,0);
for(m = 1; m <= numel; m++) {
SolidProp (m,connect,solid3d,matlfunh,eproph,&nix,ix,x,&matlfun);
SolidTopo (m,connect,solid3d,&nedofs,loc,tag);
for(i = 0; i < nedofs; i++) {
u[i] = 0.;
}
vfs_DofTabElemDof (doftab,m,nix,ix,&nedofs,lm);
vfs_SysVectorGather (disp,nedofs,lm,u);
vfe_Solid3DStrsStrn (solid3d,x,u,strs,strn);
vis_StateSetDatadv (state,m,strs);
}
vdm_LManSaveState (lman,state,rprop);
}
vdm_DataFunClose (datafun);
/* delete objects */
vdm_DataFunEnd (datafun);
vdm_NatLibEnd (natlib);
vdm_LManEnd (lman);
vis_ModelDelete(model);
vis_ModelEnd (model);
vfe_Solid3DEnd (solid3d);
vsy_HashTableInitIter (matlfunh);
while (vsy_HashTableNextIter (matlfunh,&mid,(Vobject**)&matlfun), matlfun) {
vfe_MatlFunGetObj (matlfun,(Vobject**)&linmat);
vfe_LinMatEnd (linmat);
}
vsy_HashTableForEach (matlfunh,(void(*)(Vobject*))vfe_MatlFunEnd);
vsy_HashTableEnd (matlfunh);
vfs_DofTabEnd (doftab);
vfs_SysMatrixEnd (stiff);
vfs_SysVectorEnd (disp);
vfs_SysVectorEnd (forc);
vis_GroupEnd (group);
vis_StateEnd (state);
vis_GridFunEnd (gridfun);
vis_RPropEnd (rprop);
free(lcaseids);
if(tstrain) {
free(temp);
}
return 0;
}
/*----------------------------------------------------------------------
load material
----------------------------------------------------------------------*/
static void
LoadMatl(vsy_HashTable *matlfunh, vis_MProp *mprop, Vint mid)
{
Vint i;
Vint youngflag, poissonflag;
Vint ntypes,mattype[MPROP_MAX];
Vdouble tref, alpha, values[2];
vfe_LinMat *linmat;
vfe_MatlFun *matlfun;
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vis_MPropValueType (mprop,&ntypes,mattype);
youngflag = poissonflag = 0;
for(i = 0; i < ntypes; i++) {
if(mattype[i] == MPROP_E) {
vis_MPropValueDouble (mprop,MPROP_E,&values[0]);
youngflag = 1;
} else if(mattype[i] == MPROP_NU) {
vis_MPropValueDouble (mprop,MPROP_NU,&values[1]);
poissonflag = 1;
} else if(mattype[i] == MPROP_TREF) {
vis_MPropValueDouble (mprop,MPROP_TREF,&tref);
vfe_LinMatSetRefTemp (linmat,tref);
} else if(mattype[i] == MPROP_A) {
vis_MPropValueDouble (mprop,MPROP_A,&alpha);
vfe_LinMatSetThermExp (linmat,&alpha);
}
}
/* Both Young's modulus and Poisson's ratio must be defined */
if(youngflag && poissonflag) {
vfe_LinMatSetElasProp (linmat,values);
}
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat,matlfun);
/* orthotropic material */
vsy_HashTableInsert (matlfunh,mid,matlfun);
}
/*----------------------------------------------------------------------
solid element properties
----------------------------------------------------------------------*/
static void
SolidProp(Vint m, vis_Connect *connect, vfe_Solid3D *solid3d,
vsy_HashTable *matlfunh, vsy_HashTable *eproph,
Vint *nix, Vint ix[], Vdouble x[][3], vfe_MatlFun **matlfun)
{
Vint pid, mid;
vis_EProp *eprop;
vis_ConnectElemNode (connect,m,nix,ix);
vis_ConnectCoordsdv (connect,*nix,ix,x);
vis_ConnectElemAssoc (connect,VIS_PROPID,1,&m,&pid);
if(eproph) vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
/* use material id */
if(eprop == NULL) {
vis_ConnectElemAssoc (connect,VIS_MATLID,1,&m,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)matlfun);
/* use material id from element property */
} else {
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)matlfun);
}
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,*matlfun);
}
/*----------------------------------------------------------------------
solid element topology and degrees of freedom
----------------------------------------------------------------------*/
static void
SolidTopo(Vint m, vis_Connect *connect, vfe_Solid3D *solid3d,
Vint *nedofs, Vint loc[], Vint tag[])
{
Vint shape, maxi, maxj, maxk;
Vint featech;
/* set element topology */
vis_ConnectTopology (connect,m,&shape,&maxi,&maxj,&maxk);
vfe_Solid3DSetTopology (solid3d,shape,maxi,maxj,maxk);
/* set element technology */
vis_ConnectElemAssoc (connect,VIS_FEATECH,1,&m,&featech);
/* not specified so set a reasonable technology */
if(featech == 0) {
if(shape == SYS_SHAPETET) {
featech = SYS_TECH_ISOP;
} else if(shape == SYS_SHAPEWED) {
if(maxi == 3) {
featech = SYS_TECH_ISOP;
} else {
featech = SYS_TECH_MIXED;
}
} else if(shape == SYS_SHAPEHEX) {
featech = SYS_TECH_MIXED;
}
}
vfe_Solid3DSetParami (solid3d,VFE_TECH,featech);
/* query degree of freedom information */
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
}
#include <stdio.h>
#include "vfe/vfe.h"
static Vdouble x3d[10][3] = {
{ 0., 0., 0.}, {10., 0., 0.}, { 0.,10., 0.}, { 0., 0.,10.},
{ 5., 0., 0.}, { 5., 5., 0.}, { 0., 5., 0.},
{ 0., 0., 5.}, { 5., 0., 5.}, { 0., 5., 5.} };
static Vint nix3d = 10;
static Vdouble x2d[4][3] = {
{ 0., 0., 0. }, { 10., 0., 0.}, { 10., 10., 0.}, { 0., 10., 0.}};
static Vdouble x2dxz[4][3] = {
{ 0., 0., 0. }, { 10., 0., 0.}, { 10., 0., 10.}, { 0., 0., 10.}};
static Vint nix2d = 4;
/*----------------------------------------------------------------------
Generate Gradients of Shape Functions
----------------------------------------------------------------------*/
int
main() {
vfe_Solid3D *solid3d;
vfe_Solid2D *solid2d;
Vdouble phr3d[64*64][3], phx[64*64][3], h[64*64], dj[64], phr2d[16*16][2];
Vint node, ipt, ind;
/* instance Solid3D */
printf("Solid3D:\n");
solid3d = vfe_Solid3DBegin ();
/* set solid3d topology */
vfe_Solid3DSetTopology (solid3d,SYS_SHAPETET,3,0,0);
/* compute gradients of shape functions at element nodes */
vfe_Solid3DShapeGrad (solid3d,x3d,SYS_OFF,h,phr3d,dj,phx);
for(ipt = 0; ipt < nix3d; ++ipt) {
printf("Point location = node %d, dj= %e\n",ipt+1, dj[ipt]);
if(dj[ipt] == 0.) continue;
for(node = 0; node < nix3d; ++node) {
ind = nix3d*ipt+node;
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[ind],phx[ind][0],phx[ind][1],phx[ind][2]);
}
}
/* compute gradients of shape functions at centroid */
vfe_Solid3DShapeGrad (solid3d,x3d,SYS_ON,h,phr3d,dj,phx);
printf("Point location = centroid, dj= %e\n",dj[0]);
if(dj[0] != 0.) {
for(node = 0; node < nix3d; ++node) {
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[node],phx[node][0],phx[node][1],phx[node][2]);
}
}
/* delete Solid3D */
vfe_Solid3DEnd (solid3d);
/* instance Solid2D */
printf("\nSolid2D in x-y plane:\n");
solid2d = vfe_Solid2DBegin ();
/* set solid2d topology */
vfe_Solid2DSetTopology (solid2d,SYS_SHAPEQUAD,2,0);
/* compute gradients of shape functions at element nodes */
vfe_Solid2DShapeGrad (solid2d,x2d,SYS_OFF,h,phr2d,dj,phx);
for(ipt = 0; ipt < nix2d; ++ipt) {
printf("Point location = node %d, dj= %e\n",ipt+1, dj[ipt]);
if(dj[ipt] == 0.) continue;
for(node = 0; node < nix2d; ++node) {
ind = nix2d*ipt+node;
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[ind],phx[ind][0],phx[ind][1],phx[ind][2]);
}
}
/* compute gradients of shape functions at centroid */
vfe_Solid2DShapeGrad (solid2d,x2d,SYS_ON,h,phr2d,dj,phx);
printf("Point location = centroid, dj= %e\n",dj[0]);
if(dj[0] != 0.) {
for(node = 0; node < nix2d; ++node) {
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[node],phx[node][0],phx[node][1],phx[node][2]);
}
}
/* assume 2D element on x-z plane instead */
printf("\nSolid2D in x-z plane:\n");
vfe_Solid2DSetParami (solid2d,VFE_2DXZ,SYS_ON);
/* compute gradients of shape functions at element nodes */
vfe_Solid2DShapeGrad (solid2d,x2dxz,SYS_OFF,h,phr2d,dj,phx);
for(ipt = 0; ipt < nix2d; ++ipt) {
printf("Point location = node %d, dj= %e\n",ipt+1, dj[ipt]);
if(dj[ipt] == 0.) continue;
for(node = 0; node < nix2d; ++node) {
ind = nix2d*ipt+node;
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[ind],phx[ind][0],phx[ind][1],phx[ind][2]);
}
}
/* compute gradients of shape functions at centroid */
vfe_Solid2DShapeGrad (solid2d,x2dxz,SYS_ON,h,phr2d,dj,phx);
printf("Point location = centroid, dj= %e\n",dj[0]);
if(dj[0] != 0.) {
for(node = 0; node < nix2d; ++node) {
printf(" node= %2d, %11.4e %11.4e %11.4e %11.4e\n",node+1,
h[node],phx[node][0],phx[node][1],phx[node][2]);
}
}
/* delete Solid2D */
vfe_Solid2DEnd (solid2d);
return 0;
}
The main program creates a slender cylindrical shell section spanning a 90 degree arc. Two off-section attachment points are created at either end of the shell section at the center of curvature. The attachment point at one end is linked to the shell nodes at the corresponding end using NASTRAN RBE2 like constraints. The attachment point at the other end is linked to the corresponding end using a NASTRAN RBE3 like constraint. The model geometry including element, material and solution properties are entered into the appropriate VisTools modules and registered with a single Model object.
A structure, FlexBody, is introduced which contains all instances of VfeTools modules required to support the finite element calculations of the element types, material models, etc. contained in the Model object. The FlexBodyInit function traverses the input Model object and initializes the FlexBody structure with the appropriate VfeTools objects. This FlexBody is passed to all subsequent functions to perform the degree of freedom table generation, FlexBodyDof, the flexible body mass and stiffness computations, FlexBodyGen, and various stress computations FlexBodyStr and FlexBodyStrElem. Finally all instances of VfeTools objects in FlexBody are destroyed using FlexBodyTerm.
The FlexBodyGen function performs the actual superelement mass and stiffness calculation. VfsTools objects are instanced for the required full system stiffness and mass matrices and temporary vectors. The solution procedure contained in the SProp object must be stored as solution sequence 1 in the solution property List object. The FlexBodyGen function only supports superelement creation, SYS_SOL_SUPERELEMENT.
The FlexBodyStr function illustrates stress computation for a system vector of displacements. In the example, two different uses of the FlexBodyStr function are illustrated. First an example set of generalized displacements are generated and expanded to the full system using vfs_SysVectorExpand. The stress results associated with this displacement vector are then computed and printed. Second the stresses associated with each component mode are computed.
The FlexBodyStr function consists of a loop over all elements in the model. The displacements for each element are gathered from the full system solution vector and the element node stresses and strains are computed. For shell and beam elements, stress resultants (forces and moments) are computed. Stress at the shell midsurface and beam centroid are computed by dividing the stress resultants by the shell thickness and beam area respectively. These stress quantities are then transformed to global coordinates by querying the elements for the direction cosine matrices at each element node of the coordinate axes in which the stresses are expressed. The element node stresses in global coordinates are then entered into an element node State object. When all element stresses have been computed, the stresses are processed to yield a unique stress tensor at each node point using vis_StateMap.
The FlexBodyStrElem function illustrates computing the stress matrix which relates a set of generalized coordinates for a superelement to the element node stresses for a specific element, in this case for element number 1. The procedure is similar to the stress computation procedure in function FlexBodyStr. A loop is entered over each superelement component mode. The displacements from each mode are gathered and the component stresses and strains are computed and form a column of the stress-displacement matrix.
#include <stdlib.h>
#include <math.h>
#include "base/base.h"
#include "vis/vis.h"
#include "vis/vismesh.h"
#include "vfs/vfs.h"
#include "vfe/vfe.h"
static Vint itwos[6] = {2,2,2,2,2,2};
typedef struct FlexBody {
vfe_Solid3D *solid3d;
vfe_Shell3D *shell3d;
vfe_Beam3D *beam3d;
vfe_MPC *mpc;
vsy_HashTable *matlfunh;
vfe_ShellProp *shellprop;
vfe_MatlFun *matlfunshell;
vfe_BeamProp *beamprop;
vfe_MatlFun *matlfunbeam;
Vdouble thickness[9];
Vdouble area[3], iyy[3], izz[3], iyz[3], jtr[3];
} FlexBody;
static Vint
FlexBodyInit(FlexBody *fb, vis_Model *model, Vint verbose);
static void
FlexBodyTerm(FlexBody *fb);
static Vint
FlexBodyDof(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab, Vint verbose);
static Vint
FlexBodyGen(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysMatrix *redM, vfs_SysMatrix *redK,
vfs_SysVector ***modes, vfs_SysVector *nodalMass,
Vint *pNumModes, Vint *pNrd, Vint verbose);
static void
FlexBodyStr(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector *sysV, vis_State *statent, Vint verbose);
static void
FlexBodyStrElem(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector ***modes,
Vint Nrd, Vint nelement, Vdouble sm[], Vint verbose);
static void
transStress3D(Vdouble tm[3][3], Vdouble t[6][6]);
static void
matvecprod6(Vdouble t[6][6], Vdouble a[6], Vdouble b[6]);
/*----------------------------------------------------------------------
Generate Flexible Body Using Component Mode Synthesis
----------------------------------------------------------------------*/
int
main()
{
/* FlexBody object */
FlexBody fb;
/* Modelling objects */
vis_MapMesh *mapmesh;
vis_Model *model;
vis_RCase *rcase;
vis_EProp *eprop;
vis_MProp *mprop;
vis_SProp *sprop;
vis_Connect *connect;
vsy_HashTable *mproph, *eproph, *rcaseh;
vsy_List *spropl;
/* Solution objects */
vfs_DofTab *doftab;
vfs_SysMatrix *redK, *redM;
vfs_SysVector **modes, *nodalMass;
/* Reduced solution vector */
vfs_SysVector *redV;
/* Full solution vector */
vfs_SysVector *sysV;
/* Result objects */
vis_State *stated, *statent;
vis_RProp *rprop;
/* 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;
Vint nelemcurve = 8;
Vint nelemwidth = 40;
/* number of component modes requested */
Vint nummodes = 10;
Vint Nrd, NumModes;
Vint neq, nre;
Vint verbose = 1;
/* set patch point connectivity */
Vint pix[8] = {1,2,3,4,5,6,7,8};
/* 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;
Vint *ix;
Vdouble val;
Vdouble *sm;
Vdouble u[VFE_MAXDOF];
Vint lm[VFE_MAXDOF];
Vint writeunvflag = 1;
Vdouble dzero = 0.;
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 Connect object */
connect = vis_ConnectBegin ();
vis_ConnectPre (connect,SYS_DOUBLE);
vis_ConnectDef (connect,0,0);
/* 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);
/* 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 as multiple elements */
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,MPC_KINE);
}
}
/* add rbe3 as a single element */
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,MPC_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);
/* write SDRC Universal file of basic model */
if(writeunvflag) {
vis_ConnectWrite (connect,SYS_SDRC_UNIVERSAL,"exam20dev.unv");
}
/* hash table of restraint cases */
rcaseh = vsy_HashTableBegin();
vsy_HashTableDef (rcaseh,10);
/* Suppress freedoms at rbe nodes */
rcase = vis_RCaseBegin();
for(i = 0; i < 6; i++) {
vis_RCaseSetSPCdv (rcase,rbe2node,tags[i],RCASE_APPLIED,&dzero,0);
vis_RCaseSetSPCdv (rcase,rbe3node,tags[i],RCASE_APPLIED,&dzero,0);
}
vsy_HashTableInsert (rcaseh,1,rcase);
/* hash table of element properties */
eproph = vsy_HashTableBegin();
vsy_HashTableDef (eproph,10);
/* 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 */
eprop = vis_EPropBegin();
vis_EPropDef (eprop,SYS_ELEM_RIGID);
vis_EPropSetValued (eprop,EPROP_PENALTY,1.e+9);
vis_EPropSetValuei (eprop,EPROP_NUMDIST,nedge);
vsy_HashTableInsert (eproph,2,eprop);
eprop = vis_EPropBegin();
vis_EPropDef (eprop,SYS_ELEM_RIGID);
vis_EPropSetValued (eprop,EPROP_PENALTY,1.e+9);
vis_EPropSetValuei (eprop,EPROP_NUMDIST,nedge);
vsy_HashTableInsert (eproph,3,eprop);
/* hash table of material properties */
mproph = vsy_HashTableBegin();
vsy_HashTableDef (mproph,10);
/* 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);
/* list of solution properties */
spropl = vsy_ListBegin();
vsy_ListDef (spropl,10);
/* solution properties for superelement generation */
sprop = vis_SPropBegin();
vis_SPropDef (sprop,SYS_SOL_SUPERELEMENT);
vis_SPropSetValuei (sprop,SPROP_EIGEN_NUM,nummodes);
vis_SPropSetValuei (sprop,SPROP_RCASE,1);
vis_SPropSetValuei (sprop,SPROP_MASSDIAG,SYS_ON);
vsy_ListInsert (spropl,1,(Vobject*)sprop);
/* load Model object */
model = vis_ModelBegin();
vis_ModelSetObject (model,VIS_CONNECT,connect);
vis_ModelSetHashTable (model,VIS_EPROP,eproph);
vis_ModelSetHashTable (model,VIS_MPROP,mproph);
vis_ModelSetHashTable (model,VIS_RCASE,rcaseh);
vis_ModelSetList (model,VIS_SPROP,spropl);
/* instance solver objects */
doftab = vfs_DofTabBegin();
nodalMass = vfs_SysVectorBegin();
redM = vfs_SysMatrixBegin();
redK = vfs_SysMatrixBegin();
/* initialize VfeTools filters */
FlexBodyInit (&fb,model,verbose);
/* form degree of freedom table */
FlexBodyDof (&fb,model,doftab,verbose);
/* generate flexible body */
FlexBodyGen (&fb,model,doftab,redM,redK,&modes,nodalMass,
&NumModes,&Nrd,verbose);
/* instance solution vector and fill with some numbers */
redV = vfs_SysVectorBegin ();
vfs_SysVectorDef (redV,Nrd);
for(i = 0; i < Nrd; ++i) {
val = .001*(i+1);
vfs_SysVectorScatter (redV,1,&i,&val);
}
/* instance result objects */
vfs_DofTabProcess (doftab,&neq,&nre);
sysV = vfs_SysVectorBegin ();
vfs_SysVectorDef (sysV,neq);
/* expand to full system */
vfs_SysVectorExpand (redV,Nrd,modes,sysV);
/* load node translations and rotations in state object */
stated = vis_StateBegin();
vis_StateDef (stated,numnp,SYS_NODE,SYS_NONE,VIS_SIXDOF);
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,6,tags,lm);
vfs_SysVectorGather (sysV,6,lm,u);
vis_StateSetDatadv (stated,n,u);
}
/* compute stress */
statent = vis_StateBegin();
vis_StateDef (statent,numnp,SYS_NODE,SYS_NONE,VIS_TENSOR);
FlexBodyStr (&fb,model,doftab,sysV,statent,verbose);
/* compute mode stress shapes */
if(verbose) printf("Begin computing stress modes\n");
for(i = 0; i < Nrd; ++i) {
vis_StateClear (stated);
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,6,tags,lm);
vfs_SysVectorGather (modes[i],6,lm,u);
vis_StateSetDatadv (stated,n,u);
}
vis_StateClear (statent);
if(verbose) printf(" Computing mode= %d\n",i+1);
FlexBodyStr (&fb,model,doftab,modes[i],statent,0);
/* write SDRC Universal file of displacement and stress */
if(writeunvflag) {
rprop = vis_RPropBegin();
vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
vis_RPropSetIds (rprop,i+1,0,0);
vis_RPropSetType (rprop,SYS_RES_D);
vis_StateWrite (stated,rprop,SYS_SDRC_UNIVERSAL,"exam20dev.unv");
vis_RPropSetType (rprop,SYS_RES_S);
vis_StateWrite (statent,rprop,SYS_SDRC_UNIVERSAL,"exam20dev.unv");
vis_RPropEnd (rprop);
}
}
if(verbose) printf("done!\n");
/* compute stress matrix for element number 1 */
n = 1;
sm = (Vdouble*)malloc(8*27*Nrd*sizeof(Vdouble));
FlexBodyStrElem (&fb,model,doftab,&modes,Nrd,n,sm,verbose);
free(sm);
/* terminate VfeTools filters */
FlexBodyTerm (&fb);
/* delete solution vectors */
vfs_SysVectorEnd (redV);
vfs_SysVectorEnd (sysV);
/* delete result objects */
vis_StateEnd (stated);
vis_StateEnd (statent);
/* delete solver objects */
for(i = 0; i < Nrd; ++i) {
vfs_SysVectorEnd(modes[i]);
}
free(modes);
vfs_SysMatrixEnd (redK);
vfs_SysMatrixEnd (redM);
vfs_SysVectorEnd (nodalMass);
vfs_DofTabEnd (doftab);
/* delete model objects */
vsy_HashTableForEach (rcaseh,(void(*)(Vobject*))vis_RCaseEnd);
vsy_HashTableForEach (eproph,(void(*)(Vobject*))vis_EPropEnd);
vsy_HashTableForEach (mproph,(void(*)(Vobject*))vis_MPropEnd);
vsy_ListForEach (spropl,(void(*)(Vobject*))vis_SPropEnd);
vsy_HashTableEnd (mproph);
vsy_HashTableEnd (eproph);
vsy_HashTableEnd (rcaseh);
vsy_ListEnd (spropl);
vis_ConnectEnd (connect);
vis_ModelEnd (model);
return 0;
}
/*----------------------------------------------------------------------
Initialize VfeTools objects
----------------------------------------------------------------------*/
static Vint
FlexBodyInit(FlexBody *fb, vis_Model *model, Vint verbose)
{
Vint i;
Vdouble shellprops[2] = {1., 5./6.};
Vdouble beamprops[10] = {1., 1., 1., 1., 1., 5./6., 5./6., 0., 0., 0.};
vsy_HashTable *mproph;
Vint nprops;
vis_MProp *mprop;
Vint mtype, mid;
vfe_LinMat *linmat;
vfe_MatlFun *matlfun;
Vint ntypes, mattype[MPROP_MAX];
Vint densityflag, youngflag, poissonflag;
Vdouble density, young, poisson, elasprop[2];
/* instance element formulations */
fb->solid3d = vfe_Solid3DBegin ();
fb->shell3d = vfe_Shell3DBegin ();
fb->beam3d = vfe_Beam3DBegin ();
fb->mpc = vfe_MPCBegin ();
/* shell properties */
fb->shellprop = vfe_ShellPropBegin ();
vfe_ShellPropDef (fb->shellprop,SHELLPROP_MONOCOQUE);
vfe_ShellPropSetMonocoque (fb->shellprop,shellprops,2,SYS_RULE_GAUSS,0.,1);
fb->matlfunshell = vfe_MatlFunBegin ();
vfe_ShellPropMatlFun (fb->shellprop,fb->matlfunshell);
vfe_Shell3DSetObject (fb->shell3d,VFE_MATLFUN,fb->matlfunshell);
/* shell thickness */
vfe_Shell3DSetPropPtr (fb->shell3d,VFE_PROP_THICKNESS,fb->thickness);
/* beam properties */
fb->beamprop = vfe_BeamPropBegin ();
vfe_BeamPropDef (fb->beamprop,BEAMPROP_SECTION);
vfe_BeamPropSetSection (fb->beamprop,beamprops,1);
fb->matlfunbeam = vfe_MatlFunBegin ();
vfe_BeamPropMatlFun (fb->beamprop,fb->matlfunbeam);
vfe_Beam3DSetObject (fb->beam3d,VFE_MATLFUN,fb->matlfunbeam);
/* beam sections */
vfe_Beam3DSetPropPtr (fb->beam3d,VFE_PROP_AREA,fb->area);
vfe_Beam3DSetPropPtr (fb->beam3d,VFE_PROP_IYY,fb->iyy);
vfe_Beam3DSetPropPtr (fb->beam3d,VFE_PROP_IZZ,fb->izz);
vfe_Beam3DSetPropPtr (fb->beam3d,VFE_PROP_IYZ,fb->iyz);
vfe_Beam3DSetPropPtr (fb->beam3d,VFE_PROP_J,fb->jtr);
/* extract material property hash table */
vis_ModelGetHashTable (model,VIS_MPROP,&mproph);
if(!mproph) {
if(verbose) printf("No material properties in model\n");
return 1;
}
/* create hashtable of material interfaces */
fb->matlfunh = vsy_HashTableBegin ();
vsy_HashTableCount (mproph,&nprops);
vsy_HashTableDef (fb->matlfunh,nprops);
/* loop through material properties */
vsy_HashTableInitIter (mproph);
while(vsy_HashTableNextIter (mproph,&mid,(Vobject**)&mprop),mprop) {
vis_MPropInq (mprop,&mtype);
/* only isotropic, non temperature dependent materials */
if(mtype == SYS_MAT_ISOTROPIC) {
linmat = vfe_LinMatBegin ();
vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC);
vis_MPropValueType (mprop,&ntypes,mattype);
densityflag = youngflag = poissonflag = 0;
for(i = 0; i < ntypes; i++) {
if(mattype[i] == MPROP_DENSITY) {
vis_MPropValueDouble (mprop,MPROP_DENSITY,&density);
densityflag = 1;
} else if(mattype[i] == MPROP_E) {
vis_MPropValueDouble (mprop,MPROP_E,&young);
youngflag = 1;
} else if(mattype[i] == MPROP_NU) {
vis_MPropValueDouble (mprop,MPROP_NU,&poisson);
poissonflag = 1;
}
}
if(densityflag) {
vfe_LinMatSetDensity (linmat,density);
}
/* Both Young's modulus and Poisson's ratio */
if(youngflag && poissonflag) {
elasprop[0] = young;
elasprop[1] = poisson;
vfe_LinMatSetElasProp (linmat,elasprop);
}
matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (linmat,matlfun);
}
vsy_HashTableInsert (fb->matlfunh,mid,matlfun);
}
return 0;
}
/*----------------------------------------------------------------------
Terminate VfeTools objects
----------------------------------------------------------------------*/
static void
FlexBodyTerm(FlexBody *fb)
{
Vint mid;
vfe_LinMat *linmat;
vfe_MatlFun *matlfun;
/* free material functions */
vsy_HashTableInitIter (fb->matlfunh);
while(vsy_HashTableNextIter(fb->matlfunh,&mid,(Vobject**)&matlfun),matlfun) {
vfe_MatlFunGetObj (matlfun,(Vobject**)&linmat);
vfe_LinMatEnd (linmat);
}
vsy_HashTableForEach (fb->matlfunh,(void(*)(Vobject*))vfe_MatlFunEnd);
vsy_HashTableEnd (fb->matlfunh);
vfe_MatlFunEnd (fb->matlfunbeam);
vfe_MatlFunEnd (fb->matlfunshell);
vfe_BeamPropEnd (fb->beamprop);
vfe_ShellPropEnd (fb->shellprop);
vfe_MPCEnd (fb->mpc);
vfe_Beam3DEnd (fb->beam3d);
vfe_Shell3DEnd (fb->shell3d);
vfe_Solid3DEnd (fb->solid3d);
}
/*----------------------------------------------------------------------
Generate Degree of freedom table
----------------------------------------------------------------------*/
static Vint
FlexBodyDof(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab, Vint verbose)
{
/* localized FlexBody objects */
vfe_Solid3D *solid3d;
vfe_Shell3D *shell3d;
vfe_Beam3D *beam3d;
vfe_MPC *mpc;
/* VisTools objects */
vis_Connect *connect;
vis_RCase *rcase;
vis_SProp *sprop;
/* Base objects */
vsy_HashTable *rcaseh;
vsy_List *spropl;
Vint i, n;
Vint tags[6] = {
SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TZ,
SYS_DOF_RX, SYS_DOF_RY, SYS_DOF_RZ };
Vint locone[6] = {1,1,1,1,1,1};
Vint numel, numnp;
Vint shape, maxi, maxj, maxk;
Vint featype, feaspec, featech;
Vint maxnedofs, nedofs, *loc, *tag, ntags;
Vint nix, maxnix, *ix;
Vint rcaseid;
Vint constrtype, constrmaster;
Vdouble constrval;
Vint neq, nre;
if(verbose) printf("Begin formation of global dof table\n");
/* extract number of nodes and elements */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* localize VfeTools objects */
solid3d = fb->solid3d;
shell3d = fb->shell3d;
beam3d = fb->beam3d;
mpc = fb->mpc;
/* allocate arrays with maximum number of element nodes */
vis_ConnectMaxElemNode (connect,&maxnix);
ix = (Vint*)malloc(maxnix*sizeof(Vint));
/* maximum size of non-MPC element */
maxnedofs = VFE_MAXDOF;
loc = (Vint*)malloc(VFE_MAXDOF*sizeof(Vint));
tag = (Vint*)malloc(VFE_MAXDOF*sizeof(Vint));
/* prepare DofTab object */
vfs_DofTabDef (doftab,numnp,numel);
/* loop over all elements to declare element use */
for(n = 1; n <= numel; n++) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectElemAssoc (connect,VIS_FEATYPE,1,&n,&featype);
vis_ConnectElemAssoc (connect,VIS_FEATECH,1,&n,&featech);
vis_ConnectElemAssoc (connect,VIS_FEASPEC,1,&n,&feaspec);
if(featype == SYS_ELEM_SOLID) {
vfe_Solid3DSetTopology (solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DSetParami (solid3d,VFE_TECH,featech);
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
} else if(featype == SYS_ELEM_SHELL) {
vfe_Shell3DSetTopology (shell3d,shape,maxi,maxj);
vfe_Shell3DSetParami (shell3d,VFE_TECH,featech);
vfe_Shell3DNumDof (shell3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Shell3DDofMap (shell3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
} else if(featype == SYS_ELEM_BEAM) {
vfe_Beam3DSetTopology (beam3d,maxi);
vfe_Beam3DSetParami (beam3d,VFE_TECH,featech);
vfe_Beam3DNumDof (beam3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Beam3DDofMap (beam3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
} else if(featype == SYS_ELEM_RIGID) {
vfe_MPCDef (mpc,feaspec);
/* RBE3 style */
if(feaspec == MPC_DIST) {
vfe_MPCSetDist (mpc,MPC_DIST_3D,nix-1,SYS_OFF,NULL);
/* RBE2 style */
} else if(feaspec == MPC_KINE) {
vfe_MPCSetKine (mpc,locone,tags,6,itwos,tags);
}
vfe_MPCNumDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
/* large MPCs may need to increase storage */
if(nedofs > maxnedofs) {
loc = (Vint*)realloc(loc,nedofs*sizeof(Vint));
tag = (Vint*)realloc(tag,nedofs*sizeof(Vint));
maxnedofs = nedofs;
}
vfe_MPCDofMap (mpc,VFE_ANALYSIS_STRUCTURAL,loc,tag);
}
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,n,nix,ix);
}
/* Check for permanent node restraints */
vis_ModelGetList (model,VIS_SPROP,&spropl);
vsy_ListRef (spropl,1,(Vobject**)&sprop);
vis_SPropValueInteger (sprop,SPROP_RCASE,&rcaseid);
vis_ModelGetHashTable (model,VIS_RCASE,&rcaseh);
vsy_HashTableLookup (rcaseh,rcaseid,(Vobject**)&rcase);
for(n = 1; n <= numnp; n++) {
vis_RCaseSPCTag (rcase,n,&ntags,tag);
for(i = 0; i < ntags; i++) {
vis_RCaseSPCdv (rcase,n,tag[i],&constrtype,&constrval,&constrmaster);
if(constrtype == RCASE_FIXED) {
vfs_DofTabNodePerm (doftab,n,1,&tag[i]);
}
}
}
/* finalize equation count */
vfs_DofTabProcess (doftab,&neq,&nre);
if(verbose) {
printf("Number of equations = %d\n",neq);
if(neq == 0) {
printf("No equations\n");
return 1;
}
}
/* free memory */
free(ix);
free(loc);
free(tag);
if(verbose) printf("done!\n");
return 0;
}
/*----------------------------------------------------------------------
Generate flexible body
----------------------------------------------------------------------*/
static Vint
FlexBodyGen(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysMatrix *redM, vfs_SysMatrix *redK,
vfs_SysVector ***modes, vfs_SysVector *nodalMass,
Vint *pNumModes, Vint *pNrd, Vint verbose)
{
/* localized FlexBody objects */
vfe_Solid3D *solid3d;
vfe_Shell3D *shell3d;
vfe_Beam3D *beam3d;
vfe_MPC *mpc;
vfe_ShellProp *shellprop;
vfe_BeamProp *beamprop;
vsy_HashTable *matlfunh;
/* temporary instances */
vfe_MatlFun *matlfun;
/* VfsTools objects */
vfs_SysMatrix *sysk, *sysm;
vfs_SysVector *temp, *rede, **redv;
/* VisTools objects */
vis_Connect *connect;
vis_RCase *rcase;
vis_EProp *eprop;
vis_SProp *sprop;
/* Base objects */
vsy_HashTable *rcaseh;
vsy_HashTable *eproph;
vsy_List *spropl;
/* auxiliary variables */
Vint i, j, n;
Vint tags[6] = {
SYS_DOF_TX, SYS_DOF_TY, SYS_DOF_TZ,
SYS_DOF_RX, SYS_DOF_RY, SYS_DOF_RZ };
Vint locone[6] = {1,1,1,1,1,1};
Vint dmFlag;
Vint numel, numnp, numRests;
Vint shape, maxi, maxj, maxk;
Vint featype, feaspec, featech, massflag, pid, mid;
Vint maxnedofs, nedofs, *loc, *tag, ntags;
Vint nix, maxnix, *ix;
Vint wflag, ndist, numModes, rcaseid;
Vint stype;
Vint constrtype, constrmaster;
Vint neq, nre;
Vdouble *weights;
Vdouble constrval, penalty;
Vdouble done = 1.;
Vdouble *ke, *me, *mc, (*xl)[3];
Vint *lm, *lmrest;
Vint neigen, nrd;
Vdouble *eigenval;
Vdouble value;
/* Retrieve list of solution properties */
vis_ModelGetList (model,VIS_SPROP,&spropl);
if(spropl == NULL) {
if(verbose) printf("No solution properties\n");
return 1;
}
/* assume first solution sequence */
vsy_ListRef (spropl,1,(Vobject**)&sprop);
if(sprop == NULL) {
if(verbose) printf("No solution properties\n");
return 1;
}
/* get superelement solution properties */
vis_SPropInq (sprop,&stype);
if(stype == SYS_SOL_SUPERELEMENT) {
dmFlag = 1;
vis_SPropValueInteger (sprop,SPROP_MASSDIAG,&dmFlag);
numModes = 1;
vis_SPropValueInteger (sprop,SPROP_EIGEN_NUM,&numModes);
rcaseid = 1;
vis_SPropValueInteger (sprop,SPROP_RCASE,&rcaseid);
} else {
if(verbose) printf("Only superelement solution supported\n");
return 1;
}
/* extract number of nodes and elements */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* localize VfeTools objects */
solid3d = fb->solid3d;
shell3d = fb->shell3d;
beam3d = fb->beam3d;
mpc = fb->mpc;
shellprop = fb->shellprop;
beamprop = fb->beamprop;
matlfunh = fb->matlfunh;
/* extract element property hash table */
vis_ModelGetHashTable (model,VIS_EPROP,&eproph);
if(!eproph) {
if(verbose) printf("No element properties in model\n");
return 1;
}
/* extract restraint case hash table */
vis_ModelGetHashTable (model,VIS_RCASE,&rcaseh);
if(!rcaseh) {
if(verbose) printf("No restraint case in model\n");
return 1;
}
vsy_HashTableLookup (rcaseh,rcaseid,(Vobject**)&rcase);
if(rcase == NULL) {
if(verbose) printf("RCase %d not found\n",rcaseid);
return 1;
}
/* allocate arrays with maximum number of element nodes */
vis_ConnectMaxElemNode (connect,&maxnix);
ix = (Vint*)malloc(maxnix*sizeof(Vint));
weights = (Vdouble*)malloc(maxnix*sizeof(Vdouble));
/* query for number of equations, reactions */
vfs_DofTabProcess (doftab,&neq,&nre);
/* allocate arrays with maximum size element dof map */
vfs_DofTabMaxElemDof (doftab,&maxnedofs);
loc = (Vint*)malloc(maxnedofs*sizeof(Vint));
tag = (Vint*)malloc(maxnedofs*sizeof(Vint));
/* Generate storage for largest element matrix */
ke = (Vdouble*)malloc((maxnedofs*(maxnedofs+2)/2)*sizeof(Vdouble));
me = (Vdouble*)malloc(maxnedofs*sizeof(Vdouble));
mc = (Vdouble*)malloc((maxnedofs*(maxnedofs+2)/2)*sizeof(Vdouble));
lm = (Vint*)malloc(maxnedofs*sizeof(Vint));
xl = (Vdouble(*)[3])malloc(3*maxnix*sizeof(Vdouble));
/* form stiffness and mass */
sysk = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (sysk,SYSMATRIX_SYMM_SPARSE,neq);
/* create either lumped or consistent mass object */
sysm = vfs_SysMatrixBegin();
if(dmFlag) {
vfs_SysMatrixDef (sysm,SYSMATRIX_DIAG,neq);
} else {
vfs_SysMatrixDef (sysm,SYSMATRIX_SYMM_SPARSE,neq);
}
/* preprocess and initialize stiffness and mass */
if(verbose) printf("Begin preprocessing stiffness and mass\n");
vfs_SysMatrixSetObject (sysk,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixPreProcess (sysk);
vfs_SysMatrixPreProcess (sysm);
vfs_SysMatrixInit (sysk,0.);
vfs_SysMatrixInit (sysm,0.);
if(verbose) printf("done!\n");
/* perpare vector for nodal mass assembly */
vfs_SysVectorDef (nodalMass,numnp);
/* form and assemble k and m */
if(verbose) printf("Begin assembling stiffness and mass\n");
vfs_SysMatrixInit (sysk,0.);
vfs_SysMatrixInit (sysm,0.);
for(n = 1; n <= numel; n++) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,xl);
vis_ConnectElemAssoc (connect,VIS_FEATYPE,1,&n,&featype);
vis_ConnectElemAssoc (connect,VIS_FEATECH,1,&n,&featech);
vis_ConnectElemAssoc (connect,VIS_FEASPEC,1,&n,&feaspec);
vis_ConnectElemAssoc (connect,VIS_PROPID,1,&n,&pid);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
massflag = 0;
if(featype == SYS_ELEM_SOLID) {
massflag = 1;
vfe_Solid3DSetTopology (solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DSetParami (solid3d,VFE_TECH,featech);
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
vfe_Solid3DStiff (solid3d,xl,ke);
if(!dmFlag) {
vfe_Solid3DMass (solid3d,xl,mc);
}
vfe_Solid3DMassDiag (solid3d,xl,me);
} else if(featype == SYS_ELEM_SHELL) {
massflag = 1;
vfe_Shell3DSetTopology (shell3d,shape,maxi,maxj);
vfe_Shell3DSetParami (shell3d,VFE_TECH,featech);
vfe_Shell3DNumDof (shell3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Shell3DDofMap (shell3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_ShellPropSetMatlFun (shellprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_THICKNESS,shape,maxi,maxj,fb->thickness);
vfe_Shell3DStiff (shell3d,xl,ke);
if(!dmFlag) {
vfe_Shell3DMass (shell3d,xl,mc);
}
vfe_Shell3DMassDiag (shell3d,xl,me);
} else if(featype == SYS_ELEM_BEAM) {
massflag = 1;
vfe_Beam3DSetTopology (beam3d,maxi);
vfe_Beam3DSetParami (beam3d,VFE_TECH,featech);
vfe_Beam3DNumDof (beam3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Beam3DDofMap (beam3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_BeamPropSetMatlFun (beamprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_AREA,shape,maxi,maxj,fb->area);
vis_EPropEvaldv (eprop,EPROP_IYY, shape,maxi,maxj,fb->iyy);
vis_EPropEvaldv (eprop,EPROP_IZZ, shape,maxi,maxj,fb->izz);
vis_EPropEvaldv (eprop,EPROP_IYZ, shape,maxi,maxj,fb->iyz);
vis_EPropEvaldv (eprop,EPROP_J, shape,maxi,maxj,fb->jtr);
vfe_Beam3DStiff (beam3d,xl,ke);
if(!dmFlag) {
vfe_Beam3DMass (beam3d,xl,mc);
}
vfe_Beam3DMassDiag (beam3d,xl,me);
} else if(featype == SYS_ELEM_RIGID) {
vfe_MPCDef (mpc,feaspec);
if(feaspec == MPC_DIST) {
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_NUMDIST,&ndist);
vis_EPropValueFlag (eprop,EPROP_WEIGHTS,&wflag);
if(wflag) {
vis_EPropValueDouble (eprop,EPROP_WEIGHTS,weights);
vfe_MPCSetDist (mpc,MPC_DIST_3D,ndist,wflag,weights);
} else {
vfe_MPCSetDist (mpc,MPC_DIST_3D,ndist,wflag,NULL);
}
} else if(feaspec == MPC_KINE) {
vfe_MPCSetKine (mpc,locone,tags,6,itwos,tags);
}
vfe_MPCNumDof (mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCDofMap (mpc,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueDouble (eprop,EPROP_PENALTY,&penalty);
if(penalty != 0.) {
vfe_MPCSetParamd (mpc,MPC_PENALTY,penalty);
}
vfe_MPCStiff (mpc,xl,ke);
}
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysMatrixAssem (sysk,nedofs,lm,ke);
if(massflag) {
if(dmFlag) {
vfs_SysMatrixAssemDiag (sysm,nedofs,lm,me);
} else {
vfs_SysMatrixAssem (sysm,nedofs,lm,mc);
}
for(j = 0; j < nedofs; j++) {
if(tag[j] == SYS_DOF_TX) {
vfs_SysVectorAssem (nodalMass,1,&ix[loc[j]-1],&me[j]);
}
}
}
}
if(verbose) printf("done!\n");
/* define attachment restraints */
/* count number of restraints */
numRests = 0;
for(n = 1; n <= numnp; n++) {
vis_RCaseSPCTag (rcase,n,&ntags,tag);
for(i = 0; i < ntags; i++) {
vis_RCaseSPCdv (rcase,n,tag[i],&constrtype,&constrval,&constrmaster);
if(constrtype == RCASE_APPLIED) {
numRests++;
}
}
}
if(numRests == 0) {
if(verbose) printf("No restraints defined\n");
return 1;
}
/* determine dof of restraints */
lmrest = (Vint*)malloc(numRests*sizeof(Vint));
numRests = 0;
for(n = 1; n <= numnp; n++) {
vis_RCaseSPCTag (rcase,n,&ntags,tag);
for(i = 0; i < ntags; i++) {
vis_RCaseSPCdv (rcase,n,tag[i],&constrtype,&constrval,&constrmaster);
if(constrtype == RCASE_APPLIED) {
vfs_DofTabNodeDof (doftab,n,1,&tag[i],&lmrest[numRests]);
numRests++;
}
}
}
/* form numRests force vectors due to unit displacements */
/* use modes vectors here to store the force vectors */
/* temp is a scratch vector */
temp = vfs_SysVectorBegin();
vfs_SysVectorDef (temp,neq);
/* a "1" is inserted at the unit displacement */
*modes = (vfs_SysVector**)malloc((numModes+numRests)*sizeof(vfs_SysVector*));
eigenval = (Vdouble*)malloc(numModes*sizeof(Vdouble));
for(i = 0; i < numRests; i++) {
vfs_SysVectorInit (temp,0.);
vfs_SysVectorScatter (temp,1, &lmrest[i], &done);
(*modes)[i] = vfs_SysVectorBegin();
vfs_SysVectorDef ((*modes)[i],neq);
vfs_SysMatrixProd (sysk,temp,(*modes)[i]);
vfs_SysVectorScale ((*modes)[i],(*modes)[i], -done);
}
/* set restraints as "temporary" */
vfs_SysMatrixSetTemp (sysk, numRests, lmrest);
/* process and factorize */
if(verbose) printf("Begin factorizing\n");
vfs_SysMatrixProcess (sysk);
vfs_SysMatrixFactor (sysk);
/* check for factorization problems */
if(vfs_SysMatrixError(sysk) != 0) {
if(verbose) printf("Factorization error encountered\n");
return 1;
}
if(verbose) printf("done!\n");
/* form numRests static modes */
/* a "1" is reinserted at the unit displacement */
if(verbose) printf("Begin calculating static modes\n");
for(i = 0; i < numRests; i++) {
vfs_SysVectorMove (temp, (*modes)[i]);
vfs_SysMatrixSolve (sysk, temp, (*modes)[i]);
vfs_SysVectorScatter ((*modes)[i], 1, &lmrest[i], &done);
}
if(verbose) printf("done!\n");
/* form numModes vibration modes */
if(verbose) printf("Begin calculating vibration modes\n");
vfs_SysMatrixSetParami (sysk,SYSMATRIX_EIGEN,SYSMATRIX_EIGEN_VIBE);
vfs_SysMatrixSetParami (sysk,SYSMATRIX_EIGEN_NUM,numModes);
vfs_SysMatrixSetParami (sysk,SYSMATRIX_EIGEN_TYPE,SYS_EIGEN_NEAREST);
vfs_SysMatrixEigen (sysk,sysm,&neigen,eigenval);
if(verbose) printf("done!\n");
if(verbose) {
printf("Eigenvalues of full system\n");
for(i = 0; i < neigen; i++) {
printf("eigenvalue[%d]= %f\n",i+1,eigenval[i]);
}
}
/* Warning, reduce number of component modes */
if(neigen != numModes) {
if(verbose) {
printf("Warning: number of eigenvalues less than requested\n");
}
numModes = neigen;
}
*pNumModes = numModes;
/* retrieve vibration modes */
for(i = 0; i < numModes; i++) {
(*modes)[numRests+i] = vfs_SysVectorBegin();
vfs_SysVectorDef ((*modes)[numRests+i],neq);
vfs_SysMatrixEigenMode (sysk,i+1,(*modes)[numRests+i]);
}
/* free up internal eigenvector storage */
vfs_SysMatrixEigenFree (sysk);
/* the total number of dof in the reduced system */
if(verbose) printf("Begin reduction\n");
nrd = numRests + numModes;
*pNrd = nrd;
/* turn off restraints */
vfs_SysMatrixSetTemp (sysk,0,NULL);
/* create symmetric full matrices for reduced system */
/* create reduced stiffness */
/* Assuming redK already allocated & initialized */
vfs_SysMatrixDef (redK,SYSMATRIX_SYMM_FULL,nrd);
vfs_SysMatrixPreProcess (redK);
vfs_SysMatrixReduce (sysk,nrd,numRests,numModes,lmrest,*modes,redK);
/* create reduced mass */
/* Assuming redM already allocated & initialized */
vfs_SysMatrixDef (redM,SYSMATRIX_SYMM_FULL,nrd);
vfs_SysMatrixPreProcess (redM);
vfs_SysMatrixReduce (sysm,nrd,0,numModes,NULL,*modes,redM);
if(verbose) printf("done!\n");
/* compute all eigenvalues of reduced system */
/* create reduced vector for eigenvalues */
rede = vfs_SysVectorBegin();
vfs_SysVectorDef (rede,nrd);
/* Create nrd reduced vectors for eigenvectors */
redv = (vfs_SysVector**)malloc(nrd*sizeof(vfs_SysVector*));
for(i = 0; i < nrd; i++) {
redv[i] = vfs_SysVectorBegin();
vfs_SysVectorDef (redv[i],nrd);
}
vfs_SysMatrixSetParami (redK, SYSMATRIX_EIGEN, SYSMATRIX_EIGEN_VIBE);
vfs_SysMatrixEigenAll (redK, redM, rede, redv);
/* print eigenvalues */
if(verbose) {
printf("Eigenvalues of reduced system\n");
for(i = 1; i <= nrd; i++) {
vfs_SysVectorGather (rede,1,&i,&value);
printf("eigenvalue[%d]= %e\n",i,value);
}
}
/* free objects */
for(i = 0; i < nrd; i++) {
vfs_SysVectorEnd (redv[i]);
}
free(redv);
vfs_SysVectorEnd (rede);
free(eigenval);
vfs_SysVectorEnd (temp);
free(lmrest);
vfs_SysMatrixEnd (sysm);
vfs_SysMatrixEnd (sysk);
free(xl);
free(lm);
free(mc);
free(me);
free(ke);
free(tag);
free(loc);
free(ix);
free(weights);
if(verbose) printf("Model solved\n");
return 0;
}
/*----------------------------------------------------------------------
Compute stress
----------------------------------------------------------------------*/
static void
FlexBodyStr(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector *sysV, vis_State *statent, Vint verbose)
{
/* localized FlexBody objects */
vfe_Solid3D *solid3d;
vfe_Shell3D *shell3d;
vfe_Beam3D *beam3d;
vfe_ShellProp *shellprop;
vfe_BeamProp *beamprop;
vsy_HashTable *matlfunh;
/* localized Model objects */
vis_Connect *connect;
vsy_HashTable *eproph;
vfe_MatlFun *matlfun;
vis_EProp *eprop;
/* temporary objects for stress mapping */
vis_State *stateent;
vis_GridFun *gridfun;
Vint i, j, n;
Vint numel, numnp;
Vint shape, maxi, maxj, maxk;
Vint featype, feaspec, featech, pid, mid;
Vint nix, ix[27];
Vdouble xl[27][3];
Vint nedofs, loc[VFE_MAXDOF], tag[VFE_MAXDOF], lm[VFE_MAXDOF];
Vdouble u[VFE_MAXDOF];
Vdouble strs[27*8], strn[27*8];
Vdouble lstrs[6], gstrs[9][6], tm[9][3][3], stm[6][6];
Vdouble z,z2,t,t2,factor1,factor2;
/* prepare to compute stresses */
/* extract number of nodes and elements */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* temporary objects for element node stress */
gridfun = vis_GridFunBegin();
vis_ConnectGridFun (connect,gridfun);
stateent = vis_StateBegin();
vis_StateDef (stateent,numel,SYS_ELEM,SYS_NODE,VIS_TENSOR);
vis_StateSetObject (stateent,VIS_GRIDFUN,gridfun);
/* localize VfeTools objects */
solid3d = fb->solid3d;
shell3d = fb->shell3d;
beam3d = fb->beam3d;
shellprop = fb->shellprop;
beamprop = fb->beamprop;
matlfunh = fb->matlfunh;
/* extract element property hash table */
vis_ModelGetHashTable (model,VIS_EPROP,&eproph);
/* compute stresses */
if(verbose) printf("Begin stress computation\n");
for(n = 1; n <= numel; n++) {
vis_ConnectElemAssoc (connect,VIS_FEATYPE,1,&n,&featype);
if(featype == SYS_ELEM_RIGID) {
continue;
}
if(verbose) {
printf("element= %d\n",n);
}
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,xl);
vis_ConnectElemAssoc (connect,VIS_FEATECH,1,&n,&featech);
vis_ConnectElemAssoc (connect,VIS_FEASPEC,1,&n,&feaspec);
vis_ConnectElemAssoc (connect,VIS_PROPID,1,&n,&pid);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
if(featype == SYS_ELEM_SOLID) {
vfe_Solid3DSetTopology (solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DSetParami (solid3d,VFE_TECH,featech);
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysVectorGather (sysV,nedofs,lm,u);
vfe_Solid3DStrsStrn (solid3d,xl,u,strs,strn);
if(verbose) {
for(i = 0; i < nix; i++) {
printf(" node= %d\n",ix[i]);
j = i*6;
printf(" sxx,syy,szz= %f %f %f\n",strs[j+0],strs[j+1],strs[j+2]);
printf(" sxy,syz,szx= %f %f %f\n",strs[j+3],strs[j+4],strs[j+5]);
}
}
vis_StateSetDatadv (stateent,n,strs);
} else if(featype == SYS_ELEM_SHELL) {
vfe_Shell3DSetTopology (shell3d,shape,maxi,maxj);
vfe_Shell3DSetParami (shell3d,VFE_TECH,featech);
vfe_Shell3DNumDof (shell3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Shell3DDofMap (shell3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_ShellPropSetMatlFun (shellprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_THICKNESS,shape,maxi,maxj,fb->thickness);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysVectorGather (sysV,nedofs,lm,u);
vfe_Shell3DStrsStrn (shell3d,xl,u,strs,strn);
vfe_Shell3DDirCos (shell3d,xl,u,tm);
if(verbose) {
for(i = 0; i < nix; i++) {
printf(" node= %d\n",ix[i]);
j = i*8;
printf(" Nxx,Nyy,Nxy= %f %f %f\n",strs[j+0],strs[j+1],strs[j+2]);
printf(" Mxx,Myy,Mxy= %f %f %f\n",strs[j+3],strs[j+4],strs[j+5]);
printf(" Qxz,Qyz= %f %f\n",strs[j+6],strs[j+7]);
}
}
/*
* The variable "z" is the distance from the mid-surface.
* if z == 0., then the point is at the mid-surface
* if z == +thickness/2., then the point is on the top surface
* if z == -thickness/2., then the point is on the bottom surface
*
* The example shows the stress on the top surface, as z is set to
* +thickness/2. Users can extract stresses at any point by using
* the formulae below
*/
for(i = 0; i < nix; i++) {
j = i*8;
t = fb->thickness[i];
t2 = t*t;
z = t/2.;
z2 = z*z;
factor1 = 12.*z/t2;
factor2 = 1.5*(1. - 4.*z2/t2)/t;
lstrs[0] = strs[j+0] / t + factor1*strs[j+3];
lstrs[1] = strs[j+1] / t + factor1*strs[j+4];
lstrs[2] = 0.;
lstrs[3] = strs[j+2] / t + factor1*strs[j+5];
lstrs[4] = factor2*strs[j+7];
lstrs[5] = factor2*strs[j+6];
transStress3D (tm[i],stm);
matvecprod6 (stm,lstrs,gstrs[i]);
}
vis_StateSetDatadv (stateent,n,(Vdouble*)gstrs);
} else if(featype == SYS_ELEM_BEAM) {
vfe_Beam3DSetTopology (beam3d,maxi);
vfe_Beam3DSetParami (beam3d,VFE_TECH,featech);
vfe_Beam3DNumDof (beam3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Beam3DDofMap (beam3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_BeamPropSetMatlFun (beamprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_AREA,shape,maxi,maxj,fb->area);
vis_EPropEvaldv (eprop,EPROP_IYY, shape,maxi,maxj,fb->iyy);
vis_EPropEvaldv (eprop,EPROP_IZZ, shape,maxi,maxj,fb->izz);
vis_EPropEvaldv (eprop,EPROP_IYZ, shape,maxi,maxj,fb->iyz);
vis_EPropEvaldv (eprop,EPROP_J, shape,maxi,maxj,fb->jtr);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfs_SysVectorGather (sysV,nedofs,lm,u);
vfe_Beam3DStrsStrn (beam3d,xl,u,strs,strn);
vfe_Beam3DDirCos (beam3d,xl,u,tm);
if(verbose) {
for(i = 0; i < nix; i++) {
printf(" node= %d\n",ix[i]);
j = i*6;
printf(" N,Myy,Mzz= %f %f %f\n",strs[j+0],strs[j+1],strs[j+2]);
printf(" T,Qxy,Qzx= %f %f %f\n",strs[j+3],strs[j+4],strs[j+5]);
}
}
for(i = 0; i < nix; i++) {
j = i*6;
lstrs[0] = strs[j+0] / fb->area[i];
lstrs[1] = 0.;
lstrs[2] = 0.;
lstrs[3] = strs[j+4] / fb->area[i];
lstrs[4] = 0.;
lstrs[5] = strs[j+5] / fb->area[i];
transStress3D (tm[i],stm);
matvecprod6 (stm,lstrs,gstrs[i]);
}
vis_StateSetDatadv (stateent,n,(Vdouble*)gstrs);
}
}
if(verbose) printf("done!\n");
/* map element node stresses to nodes */
vis_StateMap (statent,stateent,NULL);
/* end temporary objects */
vis_StateEnd (stateent);
vis_GridFunEnd (gridfun);
}
/*----------------------------------------------------------------------
Compute stress matrix for a single element
----------------------------------------------------------------------*/
static void
FlexBodyStrElem(FlexBody *fb, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector ***modes,
Vint Nrd, Vint nelement, Vdouble sm[], Vint verbose)
{
/* localized FlexBody objects */
vfe_Solid3D *solid3d;
vfe_Shell3D *shell3d;
vfe_Beam3D *beam3d;
vfe_ShellProp *shellprop;
vfe_BeamProp *beamprop;
vsy_HashTable *matlfunh;
/* localized Model objects */
vis_Connect *connect;
vsy_HashTable *eproph;
vfe_MatlFun *matlfun;
vis_EProp *eprop;
Vint i, j, k, n;
Vint shape, maxi, maxj, maxk;
Vint featype, feaspec, featech, pid, mid;
Vint nix, ix[27];
Vdouble xl[27][3];
Vint nedofs, loc[VFE_MAXDOF], tag[VFE_MAXDOF], lm[VFE_MAXDOF];
Vdouble u[VFE_MAXDOF];
Vint nst;
Vdouble strn[27*8];
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
/* localize VfeTools objects */
solid3d = fb->solid3d;
shell3d = fb->shell3d;
beam3d = fb->beam3d;
shellprop = fb->shellprop;
beamprop = fb->beamprop;
matlfunh = fb->matlfunh;
/* extract element property hash table */
vis_ModelGetHashTable (model,VIS_EPROP,&eproph);
/* compute element stresses */
n = nelement;
/* gather element information */
vis_ConnectElemAssoc (connect,VIS_FEATYPE,1,&n,&featype);
/* skip rigid type elements */
if(featype == SYS_ELEM_RIGID) return;
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,xl);
vis_ConnectElemAssoc (connect,VIS_FEATECH,1,&n,&featech);
vis_ConnectElemAssoc (connect,VIS_FEASPEC,1,&n,&feaspec);
vis_ConnectElemAssoc (connect,VIS_PROPID,1,&n,&pid);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
/* set up element for stress computation */
if(featype == SYS_ELEM_SOLID) {
vfe_Solid3DSetTopology (solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DSetParami (solid3d,VFE_TECH,featech);
vfe_Solid3DNumDof (solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_Solid3DSetObject (solid3d,VFE_MATLFUN,matlfun);
} else if(featype == SYS_ELEM_SHELL) {
vfe_Shell3DSetTopology (shell3d,shape,maxi,maxj);
vfe_Shell3DSetParami (shell3d,VFE_TECH,featech);
vfe_Shell3DNumDof (shell3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Shell3DDofMap (shell3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_ShellPropSetMatlFun (shellprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_THICKNESS,shape,maxi,maxj,fb->thickness);
} else if(featype == SYS_ELEM_BEAM) {
vfe_Beam3DSetTopology (beam3d,maxi);
vfe_Beam3DSetParami (beam3d,VFE_TECH,featech);
vfe_Beam3DNumDof (beam3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Beam3DDofMap (beam3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vsy_HashTableLookup (eproph,pid,(Vobject**)&eprop);
vis_EPropValueInteger (eprop,EPROP_MID,&mid);
vsy_HashTableLookup (matlfunh,mid,(Vobject**)&matlfun);
vfe_BeamPropSetMatlFun (beamprop,1,matlfun);
vis_EPropEvaldv (eprop,EPROP_AREA,shape,maxi,maxj,fb->area);
vis_EPropEvaldv (eprop,EPROP_IYY, shape,maxi,maxj,fb->iyy);
vis_EPropEvaldv (eprop,EPROP_IZZ, shape,maxi,maxj,fb->izz);
vis_EPropEvaldv (eprop,EPROP_IYZ, shape,maxi,maxj,fb->iyz);
vis_EPropEvaldv (eprop,EPROP_J, shape,maxi,maxj,fb->jtr);
}
/* get element degree of freedoms */
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
/* loop through modes */
if(verbose) printf("Stress Matrix\n");
for(i = 0; i < Nrd; i++) {
vfs_SysVectorGather ((*modes)[i],nedofs,lm,u);
if(featype == SYS_ELEM_SOLID) {
nst = 6;
j = nst*nix*i;
vfe_Solid3DStrsStrn (solid3d,xl,u,&sm[j],strn);
} else if(featype == SYS_ELEM_SHELL) {
nst = 8;
j = nst*nix*i;
vfe_Shell3DStrsStrn (shell3d,xl,u,&sm[j],strn);
} else if(featype == SYS_ELEM_BEAM) {
nst = 6;
j = nst*nix*i;
vfe_Beam3DStrsStrn (beam3d,xl,u,&sm[j],strn);
}
if(verbose) {
for(k = 0; k < nst*nix; k++) {
printf("%f ",sm[j+k]);
}
printf("\n");
}
}
}
/*----------------------------------------------------------------------
stress rotation utilities
----------------------------------------------------------------------*/
static void
transStress3D(Vdouble tm[3][3], Vdouble t[6][6])
{
/* form 6x6 transformation matrix for stress */
t[0][0] = tm[0][0]*tm[0][0];
t[0][1] = tm[0][1]*tm[0][1];
t[0][2] = tm[0][2]*tm[0][2];
t[0][3] = 2.*tm[0][0]*tm[0][1];
t[0][4] = 2.*tm[0][1]*tm[0][2];
t[0][5] = 2.*tm[2][0]*tm[0][0];
t[1][0] = tm[1][0]*tm[1][0];
t[1][1] = tm[1][1]*tm[1][1];
t[1][2] = tm[1][2]*tm[1][2];
t[1][3] = 2.*tm[1][0]*tm[1][1];
t[1][4] = 2.*tm[1][1]*tm[1][2];
t[1][5] = 2.*tm[1][2]*tm[1][0];
t[2][0] = tm[2][0]*tm[2][0];
t[2][1] = tm[2][1]*tm[2][1];
t[2][2] = tm[2][2]*tm[2][2];
t[2][3] = 2.*tm[2][0]*tm[2][1];
t[2][4] = 2.*tm[2][1]*tm[2][2];
t[2][5] = 2.*tm[2][2]*tm[2][0];
t[3][0] = tm[0][0]*tm[1][0];
t[3][1] = tm[0][1]*tm[1][1];
t[3][2] = tm[0][2]*tm[1][2];
t[3][3] = tm[0][0]*tm[1][1] + tm[0][1]*tm[1][0];
t[3][4] = tm[0][1]*tm[1][2] + tm[0][2]*tm[1][1];
t[3][5] = tm[0][2]*tm[1][0] + tm[0][0]*tm[1][2];
t[4][0] = tm[1][0]*tm[2][0];
t[4][1] = tm[1][1]*tm[2][1];
t[4][2] = tm[1][2]*tm[2][2];
t[4][3] = tm[1][0]*tm[2][1] + tm[1][1]*tm[2][0];
t[4][4] = tm[1][1]*tm[2][2] + tm[1][2]*tm[2][1];
t[4][5] = tm[1][2]*tm[2][0] + tm[1][0]*tm[2][2];
t[5][0] = tm[0][0]*tm[2][0];
t[5][1] = tm[0][1]*tm[2][1];
t[5][2] = tm[0][2]*tm[2][2];
t[5][3] = tm[0][0]*tm[2][1] + tm[0][1]*tm[2][0];
t[5][4] = tm[0][1]*tm[2][2] + tm[0][2]*tm[2][1];
t[5][5] = tm[0][2]*tm[2][0] + tm[0][0]*tm[2][2];
}
static void
matvecprod6(Vdouble t[6][6], Vdouble a[6], Vdouble b[6])
{
Vint i,j;
/* a is input local stress, b is output global stress */
for(i = 0; i < 6; i++) {
b[i] = 0.;
for(j = 0; j < 6; j++) {
b[i] += t[i][j]*a[j];
}
}
}
The main program creates a 5 x 5 x 10 rectangular part using 10 node tetrahedral solid elements. The part consists of a single isotropic, temperature independent material. There are assumed to be two attachment points, NATTACH, at which forces and moments are defined for two times, NLOADS, in the load history. These attachment point locations, and forces and moments are contained in the arrays xlex and fmex. The model geometry including element, material and solution properties are entered into the appropriate VisTools modules and registered with a single Model object. Since the distributed forces plus the inertial forces are in balance, the part may be constrained by any set of statically determinant constraints. The function, formMCase, automatically computes these 6 constraints at 3 widely separated nodes and implements them as MPC's which are stored in an MCase object.
A structure, IRelief, is introduced which contains all instances of VfeTools modules required to support the finite element calculations of the element types, material models, etc. contained in the Model object. In addition, there is auxiliary data which is used for performing the load distribution and inertia relief calculations. The IReliefInit function initializes the IRelief structure with the appropriate VisTools, VfeTools and VfsTools objects. This IRelief structure is passed as the first argument to all subsequent functions to generate the global degree of freedom table, IReliefDof, form and factor the stiffness matrix, IReliefStiff, and compute distributed load matrices relating forces and moments at attachment points to nodes on the finite element model, IReliefAttachMatrix and IReliefLoadMatrix. Then, for each set of loads at attachment points, the distributed forces and balancing inertial forces are computed using IReliefWorkLCase and IReliefFormLCase, the displacements are computed under the total load using IReliefSolve and stresses are computed using IReliefOutput. Finally all objects in the IRelief structure are deleted by calling IReliefTerm.
The technique for computing inertia relief consists of a number of steps. First the forces and moments at the attachment points must be distributed to the nodes on an associated face of the finite element model. These nodes have been marked as element associations during the mesh generation procedure. The attachment points themselves do not actually appear in the finite element mesh. The function IReliefAttachMatrix uses a MPC object to compute the load distribution for each attachment point. The distribution coefficients are stored in the IRelief structure. These coefficients need only be computed once since they are dependent only upon the geometry of each of the attachment points and their associated nodes. Second the matrix which is used to generate the inertial forces which balance the distributed loads is computed using the function IReliefLoadMatrix. The balancing inertial forces are a linear combination of six force vectors which are equal to the mass times acceleration vectors which are equivalent to the six rigid body motion vectors. This matrix is dependent only upon the geometry of the finite element model and only needs to be computed once.
Then a loop over load sets at the attachment points is entered. The forces and moments at the attachment points are distributed to the associated nodes on the finite element model using the function IReliefWorkLCase. This function uses an intermediate working LCase object to hold the distributed forces at the nodes. It may be useful for visualization purposes to use this working LCase object to plot the distributed loads upon the part. The function IReliefFormLCase transfers the loads from the working LCase to the final LCase object to which the balancing inertial forces are also added. The displacements are solved for using the function IReliefSolve. Then stresses are computed and placed in State objects along with the displacements and applied forces by the function IReliefOutput. The State objects may be used for visualization. In this example they are written out in SDRC Universal File format appended to the previously written Connect object. Finally all objects instanced in the IRelief structure are freed using IReliefTerm.
#include <stdlib.h>
#include <math.h>
#include "base/base.h"
#include "vis/vis.h"
#include "vis/vismesh.h"
#include "vfs/vfs.h"
#include "vfe/vfe.h"
typedef struct IRelief {
vis_Group *group;
vis_LCase *lcasew;
vfe_NodeGeom *nodegeom;
vfe_Solid3D *solid3d;
vfe_LinMat *linmat;
vfe_MatlFun *matlfun;
vfe_MPC *mpc;
vfs_SysMatrix *am;
vfs_SysVector *fm;
vfs_SysVector *sm;
Vdouble (**cdist)[6];
Vint *ndist;
} IRelief;
static void
formMCase(vis_Model *model, vis_MCase *mcase, Vint verbose);
static void
IReliefInit(IRelief *ir, vis_Model *model, Vint verbose);
static void
IReliefTerm(IRelief *ir, Vint nattach, Vint verbose);
static void
IReliefDof(IRelief *ir, vis_Model *model, vfs_DofTab *doftab, Vint verbose);
static void
IReliefOutput(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysVector *disp, vfs_SysVector *force,
vis_State *states, vis_State *stated, vis_State *statef,
Vint verbose);
static void
IReliefSolve(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysMatrix *stiff,
vfs_SysVector *force, vfs_SysVector *disp, Vint verbose);
static void
IReliefStiff(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysMatrix *stiff, vfs_SysMatrix *mass,
Vint diag, Vint verbose);
static void
IReliefAttachMatrix(IRelief *ir, vis_Model *model,
Vint nattach, Vdouble xlex[][3]);
static void
IReliefLoadMatrix(IRelief *ir, vis_Model *model, vfs_DofTab *doftab,
vfs_SysMatrix *mass, vfs_SysVector **inert, Vint verbose);
static void
IReliefWorkLCase(IRelief *ir, vis_Model *model,
Vint nattach, Vdouble fmex[][6]);
static void
IReliefFormLCase(IRelief *ir, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector **inert, Vint verbose);
/*----------------------------------------------------------------------
Generate Model for Inertia Relief
----------------------------------------------------------------------*/
int
main(int argc, char *argv[])
{
/* Define number of attachment points */
#define NATTACH 2
/* Define number of load sets */
#define NLOADS 2
/* IRelief object */
IRelief ir;
/* Modeling objects */
vis_MapMesh *mapmesh;
vis_Model *model;
vis_EProp *eprop;
vis_MProp *mprop;
vis_SProp *sprop;
vis_Connect *connect;
vis_LCase *lcase;
vis_MCase *mcase;
vsy_HashTable *mproph, *eproph, *lcaseh, *mcaseh;
vsy_List *spropl;
vis_GridFun *gridfun;
/* Solution objects */
vfs_DofTab *doftab;
vfs_SysMatrix *stiff, *mass;
vfs_SysVector *disp, *force, *inert[6];
/* Result objects */
vis_State *states, *stated, *statef;
vis_RProp *rprop;
/* Data used for model generation */
Vint writeunvflag = 1;
Vint verbose = 1;
Vint diag,i,numnp,numel,neq,nre,iload;
Vdouble density = 7.82e+3,
young = 2.068e+11,
poisson = 0.29;
Vint pix[8] = {1,2,3,4,5,6,7,8};
Vdouble p[8][3] =
{ { 0.,0.,0. }, {10.,0.,0. }, {10.,5.,0.}, { 0.,5.,0.},
{ 0.,0.,5. }, {10.,0.,5. }, {10.,5.,5.}, { 0.,5.,5.} };
Vdouble xlex[NATTACH][3] =
{ {-1.,2.5,2.5 }, {11.,2.5,2.5} };
Vdouble fmex[NLOADS][NATTACH][6] =
{ { {1.,0.,0.,0.,0.,0.}, { 0.,0.,0.,0.,0.,0.} },
{ {0.,1.,0.,0.,0.,0.}, {-1.,0.,0.,2.,0.,0.} } };
/* Model Generation */
connect = vis_ConnectBegin ();
vis_ConnectPre (connect,SYS_DOUBLE);
vis_ConnectDef (connect,0,0);
gridfun = vis_GridFunBegin ();
vis_ConnectGridFun (connect,gridfun);
mapmesh = vis_MapMeshBegin ();
vis_MapMeshDef (mapmesh,8,1);
for(i = 1; i <= 8; i++) {
vis_MapMeshSetPoint (mapmesh,i,p[i-1]);
}
vis_MapMeshSetPatch (mapmesh,1,VIS_SHAPEHEX,2,0,0,pix);
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,3);
vis_MapMeshSetPatchParami (mapmesh,1,MAPMESH_NUMELEMK,3);
vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID1,1,SYS_FACE,5,1);
vis_MapMeshSetPatchAssoc (mapmesh,VIS_MISCID2,1,SYS_FACE,6,2);
vis_MapMeshGenerate (mapmesh,connect);
vis_ConnectKernel(connect,0);
vis_MapMeshEnd (mapmesh);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* mark elements as solids */
for(i = 1; i <= numel; i++) {
vis_ConnectSetElemAssoc (connect,VIS_FEATYPE,i,SYS_ELEM_SOLID);
vis_ConnectSetElemAssoc (connect,VIS_PROPID,i,1);
}
/* solid element properties */
eproph = vsy_HashTableBegin ();
vsy_HashTableDef (eproph,1);
eprop = vis_EPropBegin ();
vis_EPropDef (eprop,SYS_ELEM_SOLID);
vis_EPropSetValuei (eprop,EPROP_MID,1);
vsy_HashTableInsert (eproph,1,eprop);
/* isotropic material properties */
mproph = vsy_HashTableBegin ();
vsy_HashTableDef (mproph,1);
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 for example */
lcaseh = vsy_HashTableBegin ();
vsy_HashTableDef (lcaseh,1);
lcase = vis_LCaseBegin ();
vsy_HashTableInsert (lcaseh,1,lcase);
/* MCase */
mcaseh = vsy_HashTableBegin ();
vsy_HashTableDef (mcaseh,1);
mcase = vis_MCaseBegin ();
vsy_HashTableInsert (mcaseh,1,mcase);
/* list of solution properties */
spropl = vsy_ListBegin ();
vsy_ListDef (spropl,1);
/* solution properties for static solution */
sprop = vis_SPropBegin ();
vis_SPropDef (sprop,SYS_SOL_STATIC);
vis_SPropSetValuei (sprop,SPROP_MASSDIAG,SYS_OFF);
vis_SPropSetValuei (sprop,SPROP_LCASE_NUM,1);
vis_SPropSetValuei (sprop,SPROP_LCASE,1);
vis_SPropSetValued (sprop,SPROP_LCASE_FACTOR,1.);
vsy_ListInsert (spropl,1,(Vobject*)sprop);
/* load Model object */
model = vis_ModelBegin ();
vis_ModelSetObject (model,VIS_CONNECT,connect);
vis_ModelSetHashTable (model,VIS_EPROP,eproph);
vis_ModelSetHashTable (model,VIS_MPROP,mproph);
vis_ModelSetHashTable (model,VIS_LCASE,lcaseh);
vis_ModelSetHashTable (model,VIS_MCASE,mcaseh);
vis_ModelSetList (model,VIS_SPROP,spropl);
/* write element connectivity and coordinates */
if(writeunvflag) {
vis_ConnectWrite (connect,SYS_SDRC_UNIVERSAL,"exam21dev.unv");
}
/* form statically determinant constraints as MPCs */
formMCase (model,mcase,verbose);
/* Inertia Relief Procedure */
/* extract SProp */
vis_ModelGetList (model,VIS_SPROP,&spropl);
vsy_ListRef (spropl,1,(Vobject**)&sprop);
/* extract diagonal matrix flag */
diag = 0;
vis_SPropValueInteger (sprop,SPROP_MASSDIAG,&diag);
/* initialize struct */
IReliefInit (&ir,model,verbose);
/* form degree of freedom table */
doftab = vfs_DofTabBegin ();
IReliefDof (&ir,model,doftab,verbose);
vfs_DofTabProcess (doftab,&neq,&nre);
/* instance other VfsTools objects */
/* force and displacement */
force = vfs_SysVectorBegin ();
vfs_SysVectorDef (force,neq);
disp = vfs_SysVectorBegin ();
vfs_SysVectorDef (disp,neq);
/* stiffness and mass matrices */
stiff = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (stiff,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (stiff,VFS_DOFTAB,(Vobject*)doftab);
vfs_SysMatrixInit (stiff,0.);
mass = vfs_SysMatrixBegin ();
if(diag) {
vfs_SysMatrixDef (mass,SYSMATRIX_DIAG,neq);
} else {
vfs_SysMatrixDef (mass,SYSMATRIX_SYMM_SPARSE,neq);
vfs_SysMatrixSetObject (mass,VFS_DOFTAB,(Vobject*)doftab);
}
vfs_SysMatrixInit (mass,0.);
/* displacement, force and stress results states */
stated = vis_StateBegin ();
vis_StateDef (stated,numnp,SYS_NODE,SYS_NONE,SYS_VECTOR);
vis_StateSetObject (stated,VIS_GRIDFUN,gridfun);
statef = vis_StateBegin ();
vis_StateDef (statef,numnp,SYS_NODE,SYS_NONE,SYS_VECTOR);
vis_StateSetObject (statef,VIS_GRIDFUN,gridfun);
states = vis_StateBegin ();
vis_StateDef (states,numel,SYS_ELEM,SYS_NODE,SYS_TENSOR);
vis_StateSetObject (states,VIS_GRIDFUN,gridfun);
/* form and factorize stiffness and form mass */
IReliefStiff (&ir,model,doftab,stiff,mass,diag,verbose);
/* working vectors for inertia relief calculation */
for(i = 0; i < 6; i++) {
inert[i] = vfs_SysVectorBegin ();
vfs_SysVectorDef (inert[i],neq);
}
/* form distributed load matrices */
IReliefAttachMatrix (&ir,model,NATTACH,xlex);
/* generate load matrix */
IReliefLoadMatrix (&ir,model,doftab,mass,inert,verbose);
/* loop over number of load cases */
for(iload = 0; iload < NLOADS; iload++) {
/* rebuild working LCase */
IReliefWorkLCase (&ir,model,NATTACH,fmex[iload]);
/* rebuild applied LCase */
IReliefFormLCase (&ir,model,doftab,inert,verbose);
/* solve for displacement */
IReliefSolve (&ir,model,doftab,stiff,force,disp,verbose);
/* compute stress and fill results states */
IReliefOutput (&ir,model,doftab,disp,force,states,stated,statef,verbose);
/* write displacment and stress result states */
if(writeunvflag) {
rprop = vis_RPropBegin();
vis_RPropDef (rprop,SYS_NODE,SYS_NONE);
vis_RPropSetIds (rprop,iload+1,0,0);
vis_RPropSetType (rprop,SYS_RES_D);
vis_StateWrite (stated,rprop,SYS_SDRC_UNIVERSAL,"exam21dev.unv");
vis_RPropSetType (rprop,SYS_RES_XF);
vis_StateWrite (statef,rprop,SYS_SDRC_UNIVERSAL,"exam21dev.unv");
vis_RPropDef (rprop,SYS_ELEM,SYS_NODE);
vis_RPropSetIds (rprop,iload+1,0,0);
vis_RPropSetType (rprop,SYS_RES_S);
vis_StateWrite (states,rprop,SYS_SDRC_UNIVERSAL,"exam21dev.unv");
vis_RPropEnd (rprop);
}
}
IReliefTerm (&ir,NATTACH,verbose);
/* delete objects */
vis_ConnectEnd (connect);
vis_GridFunEnd (gridfun);
vsy_HashTableEnd (eproph);
vis_EPropEnd (eprop);
vsy_HashTableEnd (mproph);
vis_MPropEnd (mprop);
vsy_HashTableEnd (lcaseh);
vis_LCaseEnd (lcase);
vsy_HashTableEnd (mcaseh);
vis_MCaseEnd (mcase);
vsy_ListEnd (spropl);
vis_SPropEnd (sprop);
vis_ModelEnd (model);
vfs_DofTabEnd (doftab);
vfs_SysVectorEnd (force);
vfs_SysVectorEnd (disp);
vfs_SysMatrixEnd (stiff);
vfs_SysMatrixEnd (mass);
vis_StateEnd (stated);
vis_StateEnd (statef);
vis_StateEnd (states);
for(i = 0; i < 6; i++) {
vfs_SysVectorEnd (inert[i]);
}
return 0;
}
/*----------------------------------------------------------------------
Form MPC's
----------------------------------------------------------------------*/
static void
formMCase(vis_Model *model, vis_MCase *mcase, Vint verbose)
{
vis_Connect *connect;
Vint i,j,n,numnp,loc[3],tags[3],nedofs,nddofs;
Vdouble xmin = 1.e+20, xmax = -1.e+20;
Vdouble ymin = 1.e+20, ymax = -1.e+20;
Vdouble zmin = 1.e+20, zmax = -1.e+20;
Vint i1,n1,n2,n3,limits[6],ix[3],tag;
Vdouble dist,distmax,x[3][3],climits[6][3],unit[3],proj[3][3],v[3],w[3];
Vdouble coeff[3],vunit[3],c;
vfe_MPC *mpc;
if(verbose) printf("Begin MPC Generation\n");
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* find nodes with min/max coordinates */
for(n = 1; n <= numnp; n++) {
vis_ConnectCoordsdv (connect,1,&n,x);
if(x[0][0] > xmax) {
xmax = x[0][0];
limits[0] = n;
climits[0][0] = x[0][0];
climits[0][1] = x[0][1];
climits[0][2] = x[0][2];
}
if(x[0][1] > ymax) {
ymax = x[0][1];
limits[1] = n;
climits[1][0] = x[0][0];
climits[1][1] = x[0][1];
climits[1][2] = x[0][2];
}
if(x[0][2] > zmax) {
zmax = x[0][2];
limits[2] = n;
climits[2][0] = x[0][0];
climits[2][1] = x[0][1];
climits[2][2] = x[0][2];
}
if(x[0][0] < xmin) {
xmin = x[0][0];
limits[3] = n;
climits[3][0] = x[0][0];
climits[3][1] = x[0][1];
climits[3][2] = x[0][2];
}
if(x[0][1] < ymin) {
ymin = x[0][1];
limits[4] = n;
climits[4][0] = x[0][0];
climits[4][1] = x[0][1];
climits[4][2] = x[0][2];
}
if(x[0][2] < zmin) {
zmin = x[0][2];
limits[5] = n;
climits[5][0] = x[0][0];
climits[5][1] = x[0][1];
climits[5][2] = x[0][2];
}
}
/* from these 6, find the 2 most far apart */
distmax = 0.;
for(i = 0; i < 6; i++) {
for(j = i+1; j < 6; j++) {
unit[0] = climits[i][0] - climits[j][0];
unit[1] = climits[i][1] - climits[j][1];
unit[2] = climits[i][2] - climits[j][2];
dist = unit[0]*unit[0] + unit[1]*unit[1] + unit[2]*unit[2];
if(dist > distmax) {
distmax = dist;
vunit[0] = unit[0];
vunit[1] = unit[1];
vunit[2] = unit[2];
n1 = limits[i];
n2 = limits[j];
i1 = i;
}
}
}
/* locate node in model furthest from line n1-n2 */
distmax = 1./sqrt(distmax);
vunit[0] *= distmax;
vunit[1] *= distmax;
vunit[2] *= distmax;
proj[0][0] = 1. - vunit[0]*vunit[0];
proj[1][0] = - vunit[1]*vunit[0];
proj[2][0] = - vunit[2]*vunit[0];
proj[0][1] = - vunit[0]*vunit[1];
proj[1][1] = 1. - vunit[1]*vunit[1];
proj[2][1] = - vunit[2]*vunit[1];
proj[0][2] = - vunit[0]*vunit[2];
proj[1][2] = - vunit[1]*vunit[2];
proj[2][2] = 1. - vunit[2]*vunit[2];
distmax = 0.;
for(n = 1; n <= numnp; n++) {
vis_ConnectCoordsdv (connect,1,&n,x);
v[0] = climits[i1][0] - x[0][0];
v[1] = climits[i1][1] - x[0][1];
v[2] = climits[i1][2] - x[0][2];
w[0] = proj[0][0]*v[0] + proj[0][1]*v[1] + proj[0][2]*v[2];
w[1] = proj[1][0]*v[0] + proj[1][1]*v[1] + proj[1][2]*v[2];
w[2] = proj[2][0]*v[0] + proj[2][1]*v[1] + proj[2][2]*v[2];
dist = v[0]*w[0] + v[1]*w[1] + v[2]*w[2];
if(dist > distmax) {
n3 = n;
distmax = dist;
}
}
/* initialize MCase and auxiliary MPC object */
vis_MCaseClear (mcase);
mpc = vfe_MPCBegin ();
/* add ground mpc to model */
vfe_MPCDef (mpc,MPC_GROUND);
vis_ConnectCoordsdv (connect,1,&n1,x);
tag = SYS_DOF_TX;
vfe_MPCSetGround (mpc,tag);
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vis_MCaseSetMPCdv (mcase,1,1,&n1,&tag,coeff,c);
tag = SYS_DOF_TY;
vfe_MPCSetGround (mpc,tag);
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vis_MCaseSetMPCdv (mcase,2,1,&n1,&tag,coeff,c);
tag = SYS_DOF_TZ;
vfe_MPCSetGround (mpc,tag);
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vis_MCaseSetMPCdv (mcase,3,1,&n1,&tag,coeff,0.);
/* plane/tang mpc */
vfe_MPCDef (mpc,MPC_PLANE);
vfe_MPCSetPlane (mpc,MPC_PLANE_TANG3D);
vfe_MPCNumDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCNumDepDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nddofs);
ix[0] = n2;
ix[1] = n1;
ix[2] = n3;
vis_ConnectCoordsdv (connect,3,ix,x);
ix[0] = n2;
ix[1] = n2;
ix[2] = n2;
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vfe_MPCDofMap (mpc,VFE_ANALYSIS_STRUCTURAL,loc,tags);
vis_MCaseSetMPCdv (mcase,4,nedofs-nddofs+1,ix,tags,coeff,0.);
/* first plane/norm mpc */
vfe_MPCDef (mpc,MPC_PLANE);
vfe_MPCSetPlane (mpc,MPC_PLANE_NORM3D);
vfe_MPCNumDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCNumDepDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nddofs);
ix[0] = n2;
ix[1] = n1;
ix[2] = n3;
vis_ConnectCoordsdv (connect,3,ix,x);
ix[0] = n2;
ix[1] = n2;
ix[2] = n2;
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vfe_MPCDofMap (mpc,VFE_ANALYSIS_STRUCTURAL,loc,tags);
vis_MCaseSetMPCdv (mcase,5,nedofs-nddofs+1,ix,tags,coeff,0.);
/* second plane/norm mpc */
vfe_MPCDef (mpc,MPC_PLANE);
vfe_MPCSetPlane (mpc,MPC_PLANE_NORM3D);
vfe_MPCNumDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCNumDepDof(mpc,VFE_ANALYSIS_STRUCTURAL,&nddofs);
ix[0] = n3;
ix[1] = n2;
ix[2] = n1;
vis_ConnectCoordsdv (connect,3,ix,x);
ix[0] = n3;
ix[1] = n3;
ix[2] = n3;
vfe_MPCPrimCoeff (mpc,1,x,coeff,&c);
vfe_MPCDofMap (mpc,VFE_ANALYSIS_STRUCTURAL,loc,tags);
vis_MCaseSetMPCdv (mcase,6,nedofs-nddofs+1,ix,tags,coeff,0.);
vfe_MPCEnd (mpc);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Initialize VfeTools objects
----------------------------------------------------------------------*/
static void
IReliefInit(IRelief *ir, vis_Model *model, Vint verbose)
{
Vint i;
vsy_HashTable *mproph;
vis_MProp *mprop;
vis_Connect *connect;
Vint numnp,nprops,mtype, mid;
Vint ntypes, mattype[MPROP_MAX];
Vint densityflag, youngflag, poissonflag;
Vdouble density, young, poisson, elasprop[2];
if(verbose) printf("Begin initialization of VfeTools objects\n");
/* extract numnp */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* instance solid element formulation */
ir->solid3d = vfe_Solid3DBegin ();
vfe_Solid3DSetParami (ir->solid3d,VFE_TECH,VFE_TECH_ISOP);
/* instance MPC for MCase */
ir->mpc = vfe_MPCBegin ();
vfe_MPCDef (ir->mpc,MPC_STAN);
/* extract material property hash table */
vis_ModelGetHashTable (model,VIS_MPROP,&mproph);
/* create hashtable of material interfaces */
vsy_HashTableCount (mproph,&nprops);
/* loop through material properties */
vsy_HashTableInitIter (mproph);
while(vsy_HashTableNextIter (mproph,&mid,(Vobject**)&mprop),mprop) {
vis_MPropInq (mprop,&mtype);
/* only isotropic, non temperature dependent materials */
if(mtype == SYS_MAT_ISOTROPIC) {
ir->linmat = vfe_LinMatBegin ();
vfe_LinMatDef (ir->linmat,SYS_MAT_ISOTROPIC);
vis_MPropValueType (mprop,&ntypes,mattype);
densityflag = youngflag = poissonflag = 0;
for(i = 0; i < ntypes; i++) {
if(mattype[i] == MPROP_DENSITY) {
vis_MPropValueDouble (mprop,MPROP_DENSITY,&density);
densityflag = 1;
} else if(mattype[i] == MPROP_E) {
vis_MPropValueDouble (mprop,MPROP_E,&young);
youngflag = 1;
} else if(mattype[i] == MPROP_NU) {
vis_MPropValueDouble (mprop,MPROP_NU,&poisson);
poissonflag = 1;
}
}
if(densityflag) {
vfe_LinMatSetDensity (ir->linmat,density);
}
/* Both Young's modulus and Poisson's ratio */
if(youngflag && poissonflag) {
elasprop[0] = young;
elasprop[1] = poisson;
vfe_LinMatSetElasProp (ir->linmat,elasprop);
}
ir->matlfun = vfe_MatlFunBegin ();
vfe_LinMatMatlFun (ir->linmat,ir->matlfun);
vfe_Solid3DSetObject (ir->solid3d,VFE_MATLFUN,ir->matlfun);
}
}
/* create 6x6 auxiliary matrix and vectors */
ir->am = vfs_SysMatrixBegin ();
vfs_SysMatrixDef (ir->am,SYSMATRIX_SYMM_FULL,6);
ir->fm = vfs_SysVectorBegin ();
vfs_SysVectorDef (ir->fm,6);
ir->sm = vfs_SysVectorBegin ();
vfs_SysVectorDef (ir->sm,6);
/* create auxiliary group and nodegeom objects */
ir->group = vis_GroupBegin ();
vis_GroupDef (ir->group,numnp,SYS_NODE,SYS_NONE);
ir->nodegeom = vfe_NodeGeomBegin ();
/* working load case */
ir->lcasew = vis_LCaseBegin ();
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Terminate VfeTools objects
----------------------------------------------------------------------*/
static void
IReliefTerm(IRelief *ir, Vint nattach, Vint verbose)
{
Vint j;
if(verbose) printf("Begin termination\n");
/* free material functions */
vfe_LinMatEnd (ir->linmat);
vfe_MatlFunEnd (ir->matlfun);
/* free VfeTools objects */
vfe_MPCEnd (ir->mpc);
vfe_Solid3DEnd (ir->solid3d);
/* free auxiliary VfsTools objects */
vfs_SysMatrixEnd (ir->am);
vfs_SysVectorEnd (ir->fm);
vfs_SysVectorEnd (ir->sm);
vis_GroupEnd (ir->group);
vfe_NodeGeomEnd (ir->nodegeom);
vis_LCaseEnd (ir->lcasew);
for(j = 0; j < nattach; j++) {
free(ir->cdist[j]);
}
free(ir->cdist);
free(ir->ndist);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Generate Degree of freedom table
----------------------------------------------------------------------*/
static void
IReliefDof(IRelief *ir, vis_Model *model, vfs_DofTab *doftab, Vint verbose)
{
/* VisTools objects */
vis_Connect *connect;
Vint n;
Vint numel, numnp, neq, nre;
Vint nterms;
Vint nedofs, loc[VFE_MAXDOF], tag[VFE_MAXDOF];
Vint shape, maxi, maxj, maxk;
Vint nix, ix[VFE_MAXNODE];
Vint maxterms, maxrhs, maxindex, icase;
Vdouble r, coeff[3];
vis_MCase *mcase;
vsy_HashTable *mcaseh;
if(verbose) printf("Begin formation of global dof table\n");
/* extract number of nodes and elements */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
vis_ModelGetHashTable (model,VIS_MCASE,&mcaseh);
vsy_HashTableLookup (mcaseh,1,(Vobject**)&mcase);
vis_MCaseMax (mcase,&maxindex,&maxterms,&maxrhs);
/* prepare DofTab object */
vfs_DofTabDef (doftab,numnp,numel+maxindex);
/* loop over all elements to declare element use */
for(n = 1; n <= numel; n++) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vfe_Solid3DSetTopology (ir->solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DNumDof (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,n,nix,ix);
}
/* loop over MPC's to declare element use */
for(icase = 1; icase <= maxindex; icase++) {
vis_MCaseMPCdv (mcase,icase,&nterms,ix,tag,coeff,&r);
vfe_MPCSetStan (ir->mpc,nterms,tag,coeff,r);
vfe_MPCNumDof (ir->mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCDofMap (ir->mpc,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vfs_DofTabSetMap (doftab,nedofs,loc,tag);
vfs_DofTabElemUse (doftab,numel+icase,nix,ix);
}
/* finalize equation count */
vfs_DofTabProcess (doftab,&neq,&nre);
if(verbose) printf("number of equations= %d\n",neq);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Form output
----------------------------------------------------------------------*/
static void
IReliefOutput(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysVector *disp, vfs_SysVector *force,
vis_State *states, vis_State *stated, vis_State *statef,
Vint verbose)
{
Vint i;
Vint n, numel, numnp;
Vint nedofs, loc[VFE_MAXDOF],tag[VFE_MAXDOF],lm[VFE_MAXDOF];
Vint nix, ix[VFE_MAXNODE];
Vint shape, maxi, maxj, maxk;
Vint tags[3] = {SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ};
Vdouble x[VFE_MAXNODE][3], u[VFE_MAXDOF];
Vdouble strs[6*VFE_MAXNODE], strn[6*VFE_MAXNODE];
vis_Connect *connect;
if(verbose) printf("Begin results computation\n");
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* compute and output element stresses */
vis_StateClear (states);
for(n = 1; n <= numel; n++) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,x);
vfe_Solid3DSetTopology (ir->solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DNumDof (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
for(i = 0; i < nedofs; i++) u[i] = 0.;
vfs_SysVectorGather (disp,nedofs,lm,u);
vfe_Solid3DStrsStrn (ir->solid3d,x,u,strs,strn);
vis_StateSetDatadv (states,n,strs);
}
/* output displacement */
vis_StateClear (stated);
for(n = 1; n <= numnp; n++) {
u[0] = u[1] = u[2] = 0.;
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
vfs_SysVectorGather (disp,3,lm,u);
vis_StateSetDatadv (stated,n,u);
}
/* output applied force */
vis_StateClear (statef);
for(n = 1; n <= numnp; n++) {
u[0] = u[1] = u[2] = 0.;
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
vfs_SysVectorGather (force,3,lm,u);
vis_StateSetDatadv (statef,n,u);
}
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Solve for displacements
----------------------------------------------------------------------*/
static void
IReliefSolve(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysMatrix *stiff,
vfs_SysVector *force, vfs_SysVector *disp, Vint verbose)
{
Vint i;
vsy_HashTable *lcaseh;
vsy_List *spropl;
vis_SProp *sprop;
vis_Connect *connect;
vis_LCase *lcase;
Vdouble factor;
Vdouble f[3],x[3];
Vint numnp,n,flags,ntypes,type[4];
Vint lm[3];
Vint tags[3] = { SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ };
if(verbose) printf("Begin load formation\n");
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* extract SProp object */
vis_ModelGetList (model,VIS_SPROP,&spropl);
vsy_ListRef (spropl,1,(Vobject**)&sprop);
vis_SPropValueDouble (sprop,SPROP_LCASE_FACTOR,&factor);
/* extract initial load case */
vis_ModelGetHashTable (model,VIS_LCASE,&lcaseh);
vsy_HashTableLookup (lcaseh,1,(Vobject**)&lcase);
/* initialize force vector */
vfs_SysVectorInit (force,0.0);
/* add force from point loads */
vis_LCaseNodeGroup (lcase,NULL,ir->group);
vis_GroupInitIndex (ir->group);
while (vis_GroupNextIndex (ir->group,&n,&flags), n) {
vis_ConnectCoordsdv (connect,1,&n,&x);
vis_LCaseConcType (lcase,n,&ntypes,type);
for(i = 0; i < ntypes; i++) {
if(type[i] == LCASE_FORCE) {
vis_LCaseConcdv (lcase,n,LCASE_FORCE,f);
f[0] *= factor;
f[1] *= factor;
f[2] *= factor;
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
vfs_SysVectorAssem (force,3,lm,f);
}
}
}
/* solve linear system */
vfs_SysMatrixSolve (stiff,force,disp);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Form and factorize stiffness, and form mass
----------------------------------------------------------------------*/
static void
IReliefStiff(IRelief *ir, vis_Model *model,
vfs_DofTab *doftab, vfs_SysMatrix *stiff, vfs_SysMatrix *mass,
Vint diag, Vint verbose)
{
vis_Connect *connect;
vis_MCase *mcase;
vsy_HashTable *mcaseh;
Vint n,numel,numnp;
Vint nedofs,lm[VFE_MAXDOF],loc[VFE_MAXDOF],tag[VFE_MAXDOF];
Vint nix, ix[VFE_MAXNODE];
Vint shape, maxi, maxj, maxk;
Vint maxindex,maxterms,maxrhs,icase,nterms;
Vdouble x[VFE_MAXNODE][3];
Vdouble ke[VFE_MAXDOF*(VFE_MAXDOF+1)/2];
Vdouble ms[VFE_MAXDOF*(VFE_MAXDOF+1)/2];
Vdouble coeff[3],r;
if(verbose) printf("Begin stiffness\n");
/* extract number of nodes and elements */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* loop over all elements to form stiffness and mass */
for(n = 1; n <= numel; n++) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,x);
vfe_Solid3DSetTopology (ir->solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DNumDof (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_Solid3DDofMap (ir->solid3d,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vfs_DofTabElemDof (doftab,n,nix,ix,&nedofs,lm);
vfe_Solid3DStiff (ir->solid3d,x,ke);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
if(diag) {
vfe_Solid3DMassDiag (ir->solid3d,x,ms);
vfs_SysMatrixAssemDiag (mass,nedofs,lm,ms);
} else {
vfe_Solid3DMass (ir->solid3d,x,ms);
vfs_SysMatrixAssem (mass,nedofs,lm,ms);
}
}
/* assemble MPC's */
vis_ModelGetHashTable (model,VIS_MCASE,&mcaseh);
vsy_HashTableLookup (mcaseh,1,(Vobject**)&mcase);
vis_MCaseMax (mcase,&maxindex,&maxterms,&maxrhs);
for(icase = 1; icase <= maxindex; icase++) {
vis_MCaseMPCdv (mcase,icase,&nterms,ix,tag,coeff,&r);
vis_ConnectCoordsdv (connect,nterms,ix,x);
vfe_MPCSetStan (ir->mpc,nterms,tag,coeff,r);
vfe_MPCNumDof (ir->mpc,VFE_ANALYSIS_STRUCTURAL,&nedofs);
vfe_MPCDofMap (ir->mpc,VFE_ANALYSIS_STRUCTURAL,loc,tag);
vfs_DofTabElemDof (doftab,numel+icase,nterms,ix,&nedofs,lm);
vfe_MPCStiff (ir->mpc,x,ke);
vfs_SysMatrixAssem (stiff,nedofs,lm,ke);
}
/* process and factorize */
vfs_SysMatrixProcess (stiff);
vfs_SysMatrixFactor (stiff);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Form attachment point matrices
----------------------------------------------------------------------*/
static void
IReliefAttachMatrix(IRelief *ir, vis_Model *model,
Vint nattach, Vdouble xlex[][3])
{
vis_Connect *connect;
Vint i,j,k,n,numnp,numel,nchild,flags,numdofs,numdepdofs;
vfe_MPC *mpcdist;
Vdouble (*xdist)[3],*coeff,c;
#define WEIGHTS
#ifdef WEIGHTS
Vdouble xl[VFE_MAXNODE][3];
Vdouble ones[VFE_MAXNODE] = {1.,1.,1.,1.,1.,1.,1.,1.,1.,1.};
Vdouble load[VFE_MAXNODE];
vis_Group *facegroup;
vfs_SysVector *nodal;
Vint shape,maxi,maxj,maxk,ix[VFE_MAXNODE],nix,dof;
Vdouble *weight;
#endif
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
/* instance MPC for distributing load */
mpcdist = vfe_MPCBegin ();
vfe_MPCDef (mpcdist,MPC_DIST);
/* allocate array of matrices */
ir->cdist = (Vdouble(**)[6])malloc(nattach*sizeof(Vdouble(**)[6]));
ir->ndist = (Vint*)malloc(nattach*sizeof(Vint));
#ifdef WEIGHTS
facegroup = vis_GroupBegin ();
vis_GroupDef (facegroup,numel,SYS_ELEM,SYS_FACE);
nodal = vfs_SysVectorBegin ();
vfs_SysVectorDef (nodal,numnp);
#endif
/* loop through points of applied load */
for(j = 0; j < nattach; j++) {
/* distribute forces and moments to surrounding nodes */
vis_GroupClear (ir->group);
vis_ConnectSetGroupParami (connect,CONNECT_ASSOCTYPE,VIS_MISCID1+j);
vis_ConnectSetGroupParami (connect,CONNECT_ASSOCID,j+1);
vis_ConnectNodeGroup (connect,CONNECT_ASSOC,NULL,ir->group);
vis_GroupCount (ir->group,&ir->ndist[j],&nchild);
#ifdef WEIGHTS
vis_GroupClear (facegroup);
vis_ConnectFaceGroup (connect,CONNECT_CONTAINED,ir->group,facegroup);
vfs_SysVectorInit (nodal,0.);
vis_GroupInitIndex (facegroup);
while(vis_GroupNextIndex (facegroup,&n,&flags), n) {
for(k = 1, dof = 1; k <= 4; k++, dof *= 2) {
if(flags & dof) {
vis_ConnectTopology (connect,n,&shape,&maxi,&maxj,&maxk);
vis_ConnectElemNode (connect,n,&nix,ix);
vis_ConnectCoordsdv (connect,nix,ix,xl);
vfe_Solid3DSetTopology (ir->solid3d,shape,maxi,maxj,maxk);
vfe_Solid3DDistHeat (ir->solid3d,xl,SYS_FACE,k,ones,load);
vfs_SysVectorAssem (nodal,nix,ix,load);
}
}
}
weight = (Vdouble*)malloc(ir->ndist[j]*sizeof(Vdouble));
vis_GroupInitIndex (ir->group);
for(k = 0; k < ir->ndist[j]; k++) {
vis_GroupNextIndex (ir->group,&n,&flags);
vfs_SysVectorGather (nodal,1,&n,&weight[k]);
if(weight[k] < 0.) weight[k] = -weight[k];
}
#endif
/* allocate space for attachment info */
ir->cdist[j] = (Vdouble(*)[6])malloc(ir->ndist[j]*3*6*sizeof(Vdouble));
xdist = (Vdouble(*)[3])malloc((1+ir->ndist[j])*3*sizeof(Vdouble));
#ifdef WEIGHTS
vfe_MPCSetDist (mpcdist,MPC_DIST_3D,ir->ndist[j],1,weight);
#else
vfe_MPCSetDist (mpcdist,MPC_DIST_3D,ir->ndist[j],SYS_OFF,NULL);
#endif
/* set first point */
xdist[0][0] = xlex[j][0];
xdist[0][1] = xlex[j][1];
xdist[0][2] = xlex[j][2];
/* set distributed points */
vis_GroupInitIndex (ir->group);
for(k = 1; k <= ir->ndist[j]; k++) {
vis_GroupNextIndex (ir->group,&n,&flags);
vis_ConnectCoordsdv (connect,1,&n,(Vdouble(*)[3])xdist[k]);
}
vfe_MPCNumDof (mpcdist,SYS_ANALYSIS_STRUCTURAL,&numdofs);
vfe_MPCNumDepDof (mpcdist,SYS_ANALYSIS_STRUCTURAL,&numdepdofs);
coeff = (Vdouble*)malloc(numdofs*sizeof(Vdouble));
for(n=0, k=0; k
vfe_MPCPrimCoeff (mpcdist,k+1,xdist,coeff,&c);
for(i=numdepdofs; i
ir->cdist[j][i-numdepdofs][k] = coeff[i-numdepdofs+1];
}
}
free(coeff);
/* free space used for attachment info */
free(xdist);
#ifdef WEIGHTS
free(weight);
#endif
}
vfe_MPCEnd (mpcdist);
#ifdef WEIGHTS
vis_GroupEnd (facegroup);
vfs_SysVectorEnd (nodal);
#endif
}
/*----------------------------------------------------------------------
Form inertia relief load
----------------------------------------------------------------------*/
static void
IReliefLoadMatrix(IRelief *ir, vis_Model *model, vfs_DofTab *doftab,
vfs_SysMatrix *mass, vfs_SysVector **inert, Vint verbose)
{
vis_Connect *connect;
vfs_SysVector *accel;
Vdouble force[3],moment[3],x[3],matrix[6][6],matrixlt[21],fmx[3];
Vdouble u[3];
Vint numnp,n,neq,type;
Vint lmf[6] = {1,2,3,4,5,6};
Vint i,j,lm[3];
Vint tags[3] = { SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ };
if(verbose) printf("Begin inertia relief load formation\n");
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
/* form six rigid body accelerations and forces */
vfs_SysMatrixInq (mass,&type,&neq);
accel = vfs_SysVectorBegin ();
vfs_SysVectorDef (accel,neq);
for(i = 1; i <= 6; i++) {
vfs_SysVectorInit (accel,0.);
for(n = 1; n <= numnp; n++) {
vis_ConnectCoordsdv (connect,1,&n,&x);
vfe_NodeGeomRigid (ir->nodegeom,i,x,3,tags,u);
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
vfs_SysVectorScatter (accel,3,lm,u);
}
vfs_SysMatrixProd (mass,accel,inert[i-1]);
}
vfs_SysVectorEnd (accel);
/* form small linear system */
for(i = 0; i < 6; i++) {
for(j = 0; j < 6; j++) {
matrix[i][j] = 0.;
}
}
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
for(i = 0; i < 6; i++) {
for(j = 0; j < 3; j++) fmx[j] = 0.;
vfs_SysVectorGather (inert[i],3,lm,fmx);
vis_ConnectCoordsdv (connect,1,&n,&x);
vfe_NodeGeomForceMoment (ir->nodegeom,x,3,tags,fmx,force,moment);
matrix[i][0] += force[0];
matrix[i][1] += force[1];
matrix[i][2] += force[2];
matrix[i][3] += moment[0];
matrix[i][4] += moment[1];
matrix[i][5] += moment[2];
}
}
for(n = 0, j = 0; n < 6; n++) {
for(i = 0; i <= n; i++) {
matrixlt[j++] = matrix[i][n];
}
}
vfs_SysMatrixPreProcess (ir->am);
vfs_SysMatrixInit (ir->am,0.);
vfs_SysMatrixAssem (ir->am,6,lmf,matrixlt);
vfs_SysMatrixProcess (ir->am);
vfs_SysMatrixFactor (ir->am);
if(verbose) printf("done!\n");
}
/*----------------------------------------------------------------------
Distribute load from attachment points
----------------------------------------------------------------------*/
static void
IReliefWorkLCase(IRelief *ir, vis_Model *model,
Vint nattach, Vdouble fmex[][6])
{
vis_Connect *connect;
Vint i,j,k,n,flags;
Vdouble f[3];
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
/* rebuild working LCase */
vis_LCaseClear (ir->lcasew);
for(j = 0; j < nattach; j++) {
vis_GroupClear (ir->group);
vis_ConnectSetGroupParami (connect,CONNECT_ASSOCTYPE,VIS_MISCID1+j);
vis_ConnectSetGroupParami (connect,CONNECT_ASSOCID,j+1);
vis_ConnectNodeGroup (connect,CONNECT_ASSOC,NULL,ir->group);
vis_GroupInitIndex (ir->group);
for(k = 0; k < ir->ndist[j]; k++) {
f[0] = 0.; f[1] = 0.; f[2] = 0.;
for(i = 0; i < 6; i++) {
f[0] += ir->cdist[j][3*k ][i] * fmex[j][i];
f[1] += ir->cdist[j][3*k+1][i] * fmex[j][i];
f[2] += ir->cdist[j][3*k+2][i] * fmex[j][i];
}
vis_GroupNextIndex (ir->group,&n,&flags);
vis_LCaseAddConcdv (ir->lcasew,n,LCASE_FORCE,f);
}
}
}
/*----------------------------------------------------------------------
Form inertia relief load
----------------------------------------------------------------------*/
static void
IReliefFormLCase(IRelief *ir, vis_Model *model, vfs_DofTab *doftab,
vfs_SysVector **inert, Vint verbose)
{
vsy_HashTable *lcaseh;
vsy_List *spropl;
vis_SProp *sprop;
vis_Connect *connect;
vis_LCase *lcase;
Vdouble ftotal[6],coeff[6],factor;
Vdouble f[3],aforce[3],moment[3],x[3];
Vint numel,numnp,n,flags,ntypes,type[4];
Vint i,lm[3],lmf[6]={1,2,3,4,5,6};
Vint tags[6] = { SYS_DOF_TX,SYS_DOF_TY,SYS_DOF_TZ };
if(verbose) printf("Begin load formation\n");
/* extract connect object */
vis_ModelGetObject (model,VIS_CONNECT,(Vobject**)&connect);
vis_ConnectNumber (connect,SYS_NODE,&numnp);
vis_ConnectNumber (connect,SYS_ELEM,&numel);
/* extract SProp object for load factor */
vis_ModelGetList (model,VIS_SPROP,&spropl);
vsy_ListRef (spropl,1,(Vobject**)&sprop);
vis_SPropValueDouble (sprop,SPROP_LCASE_FACTOR,&factor);
/* extract load case */
vis_ModelGetHashTable (model,VIS_LCASE,&lcaseh);
vsy_HashTableLookup (lcaseh,1,(Vobject**)&lcase);
/* initialize LCase */
vis_LCaseClear (lcase);
/* add force and moment from point loads */
vfs_SysVectorInit (ir->fm,0.0);
for(i = 0; i < 6; i++) {
ftotal[i] = 0.;
}
vis_LCaseNodeGroup (ir->lcasew,NULL,ir->group);
vis_GroupInitIndex (ir->group);
while (vis_GroupNextIndex (ir->group,&n,&flags), n) {
vis_ConnectCoordsdv (connect,1,&n,&x);
vis_LCaseConcType (ir->lcasew,n,&ntypes,type);
for(i = 0; i < ntypes; i++) {
if(type[i] == LCASE_FORCE) {
/* extract given concentrated load */
vis_LCaseConcdv (ir->lcasew,n,LCASE_FORCE,f);
/* add given concentrated load to working LCase */
vis_LCaseAddConcdv (lcase,n,LCASE_FORCE,f);
/* continue computation of balancing load */
f[0] *= factor;
f[1] *= factor;
f[2] *= factor;
vfe_NodeGeomForceMoment (ir->nodegeom,x,3,tags,f,aforce,moment);
ftotal[0] += aforce[0];
ftotal[1] += aforce[1];
ftotal[2] += aforce[2];
ftotal[3] += moment[0];
ftotal[4] += moment[1];
ftotal[5] += moment[2];
}
}
}
vfs_SysVectorAssem (ir->fm,6,lmf,ftotal);
if(verbose) {
printf("Applied force: %e, %e, %e\n",ftotal[0],ftotal[1],ftotal[2]);
printf("Applied moment: %e, %e, %e\n",ftotal[3],ftotal[4],ftotal[5]);
}
vfs_SysMatrixSolve (ir->am,ir->fm,ir->sm);
vfs_SysVectorGather (ir->sm,6,lmf,coeff);
/* update LCase */
for(i = 0; i < 6; i++) {
for(n = 1; n <= numnp; n++) {
vfs_DofTabNodeDof (doftab,n,3,tags,lm);
f[0] = f[1] = f[2] = 0.;
vfs_SysVectorGather (inert[i],3,lm,f);
f[0] *= -coeff[i];
f[1] *= -coeff[i];
f[2] *= -coeff[i];
vis_LCaseAddConcdv (lcase,n,LCASE_FORCE,f);
}
}
if(verbose) printf("done!\n");
}