14. Examples
This chapter contains example VfeTools programs. The source code for all examples resides in CeetronSAM/src/sam/vfe/exam.
14.1. Example 1, Generate Stiffness, Mass and Loads for a 3D 10 Node Tetrahedron
This example illustrates forming the stiffness matrix, consistent and diagonal mass matrices and the load vectors for two distributed loads for a 3D 10 node tetrahedron. Three objects are required; a LinMat object to hold the material properties, a Solid3D object to form the element stiffness, mass, etc. and a MatlFun object to provide the interface between the element formulation, Solid3D, and the underlying material model, LinMat.
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 "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.2. Example 1a, Query for Element Degrees of Freedom
This example uses the 10 node tetrahedron introduced in Example 1 and illustrates the procedure for querying for the number and type of degrees of freedom, called the dof “map”, used by the element.
#include <stdlib.h>
#include <stdio.h>
#include "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.3. Example 2, Generate Stiffness, Loads for a 2D Plane Stress Cantilever Beam
This example illustrates generating the stiffness and loads for a 2D plane stress cantilever beam consisting of 10 6 node triangles.
#include <stdio.h>
#include <stdlib.h>
#include "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
}
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
}
Figure 9-1, 2D Plane Stress Cantilever Beam Mesh
14.4. Example 3, Generate Stiffness, Mass and Loads for a 3D 4 Node Shell
This example illustrates forming the stiffness matrix, consistent and diagonal mass matrices and the load vector for a pressure load for a 3D 4 node quadrilateral shell element. Five objects are required; a LinMat object to hold the material properties, a ShellProp object to hold information related to the element properties. a Shell3D object to form the element stiffness, mass, etc. and two MatlFun objects to provide the interfaces between the element formulation, Shell3D, the element properties, ShellProp,and the underlying material model, LinMat.
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 "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.5. Example 3a, Using a Custom MatlFun Object with Example 3
This example illustrates using a custom MatlFun (material functions) object to interface a user defined material model to an element formulation. This is a modification of Example 3 in which the material model for a 4 node quadrilateral shell element was implemented using the ShellProp and LinMat material models. The custom material model is designed to reproduce exactly the material model used in Example 3, so that the computed stiffness and mass matrices are identical.
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 "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.6. Example 4, Generate Stiffness, Mass and Loads for a 3D 3 Node Beam
This example illustrates forming the stiffness matrix, consistent and diagonal mass matrices and the load vector for a line load for a 3D 3 node beam element. Five objects are required; a LinMat object to hold the material properties, a BeamProp object to hold information related to the element properties. a Beam3D object to form the element stiffness, mass, etc. and two MatlFun objects to provide the interfaces between the element formulation, Beam3D, the element properties, BeamProp,and the underlying material model, LinMat.
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 "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.7. Example 5, Rigid Body MPC Using Penalty Function
This example illustrates generating multi-point constraints using
the MPC module.
There are three multi-point constraints in this example which form a tripod
and which could be used to transmit forces and moments applied to the
single node at the tip of the tripod to the 3 nodes at the base
of the tripod. This would be useful to attach a beam element to the face of
a tetrahedral solid element.
The type of constraint to be modeled is defined using vfe_MPCDef().
The rigid body type constraint is defined which can generally relate
the translations and rotations of one node to the translations and rotations
of another node which is not coincident with the first node.
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 "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.8. Example 6, Infinitesimal Elasto-plasticity under Biaxial Stretching
This example illustrates using the PlasMat module to extract the stress originating from a cyclic plastic deformation under a bi-axial state of strain. In the example, all calls to a MatlFun function assume that the PlasMat object has been registered in the MatlFun object.
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 "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.9. Example 7, Mooney-Rivlin Material in Plane Stress
This example illustrates using the HyperMat module to extract the Cauchy stress and the tangent material matrix to be used in a non-linear finite element program. It is assumed that a plane stress state is required, as is the case for plates and shells. The following operations are performed: instance a HyperMat object; set its material properties, including the bulk modulus; set condensation flags for the plane stress condition; and extract the stress and material tangent matrix for a given deformation gradient.
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 "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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};
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.10. Example 8, Corotational Frame for a 3D Solid Element
This example illustrates the use of a Corot object to compute a corotational
frame for 3D solid elements.
A corotational frame may be used to isolate the effects of large rigid body
motions, in particular large rotations, from an element formulation.
A Corot object is instanced and applied to the formation of stiffness
and reaction force computed using a Solid3D object. The setup of the
Solid3D object for computing stiffness and force follows standard practice.
The Corot object is first used to transform the element coordinates to an
initial corotational frame embedded in the element in its undeformed
configuration. This computation is performed using vfe_CorotPull().
A set of element displacements representing a small shearing and a large 60
degree rotation is manufactured and transformed to the corotational system
representing the deformed state using vfe_CorotPush(). These transformed
displacements along with the transformed element coordinates computed
by vfe_CorotPull() are now used as input to the standard
vfe_Solid3DReactStiff() function to compute stiffness and reactions in the
corotational frame. These element stiffness and force vectors must now
be transformed back to the global system using vfe_CorotPush() for
normal assembly, factorization and solution.
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.11. Example 9, Rigid Body Mode Detection using NodeGeom
This example illustrates the use of a NodeGeom object to compute unconstrained rigid body modes. A set of four constraints are applied to two nodes. A translational constraint in the y direction and a rotational constraint about the x axis are applied at the first node located at (1.,0.,0.), Translational constraints in the x direction and the y direction are applied at the second node located at (1.,1.,0.), Since there is no translational constraint in the z direction there will be a translational free mode in the z direction. Translational constraints can prevent rotational rigid body modes if there are two translational constraints in the same direction which are offset in space in a direction perpendicular to their direction. The two translational constraints in the y direction are offset parallel to their direction and therefore to not prevent rotation. As a result, even though there are 4 constraints, they only prevent 3 rigid body motions. Two rotational rigid body motions remain, one about the y axis passing through the specified center point and a second about the z axis passing through a point offset by 1. from the center point in the y direction, ie. the location of node 2.
#include "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.12. Example 10, Linear Stress Analysis Using 10 Node Tetrahedra
This example illustrates the use of VfeTools and VfsTools to develop the basic framework for linear static analysis using 10 node tetrahedral solid elements. The finite element model is a 1. x 2. x 3. block with a uniform 3 x 3 x 3 grid of nodes and connected by 6 10-node tetrahedral solid elements. Permanent restraints are applied so that the x = 0., y = 0. and z = 0. planes are symmetry planes. Uniform pressure load is applied to the element faces at the x = 1. plane. The material properties are linear elastic with Young’s modulus = 1.e+6 and Poisson’s ratio = .3 . The finite element model data structures appear as simple static arrays in this example.
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.
Figure 10-1, Finite Element Model of 10-Node Tetrahedra
#include <stdlib.h>
#include "sam/base/basedefs.h"
#include "sam/vfe/vfe.h"
#include "sam/vfs/vfs.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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, (Vchar*)"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;
}
14.13. Example 10a, Linear Stress Analysis Using 3D Solid Elements
This example is a generalization of Example 10 in which a finite element model definition is loaded from a VdmTools native .vdm file into a VisTools Model object hierarchy. The solution sequences defined in the finite element model are solved and the resulting displacements and stresses are saved on the input .vdm file. The finite element model may contain any mixture of 3D solid elements. The load cases may consist of distributed tractions on element faces and thermal loads. This example could be extended in a straight forward way to include concentrated nodal forces, body forces, etc.
An initial .vdm file containing a finite element model may be generated using VisTools example vis/exam/exam49vdm.cpp. The resulting .vdm file may be visualized using VisTools example vis/exam/exam30edev.cpp.
#include <stdlib.h>
#include "sam/base/base.h"
#include "sam/vfe/vfe.h"
#include "sam/vfs/vfs.h"
#include "sam/vis/vis.h"
#include "sam/vdm/vdm.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.h"
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], flagno;
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 = NULL, nlcases;
Vint itrac, tstrain = 0;
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]);
}
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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++) {
vis_GroupGetEntFlag(group, m, no, &flagno);
if (flagno == 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, (Vchar*)"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, (Vchar*)"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 = NULL;
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);
}
14.14. Example 11, Generate Gradients of Shape Functions
This example illustrates computing shape functions and gradients
at element nodes and element centroid in solid elements
using vfe_Solid2DShapeGrad() and vfe_Solid3DShapeGrad().
#include <stdio.h>
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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;
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* 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;
}
14.15. Example 20, Generate Flexible Body Using Component Mode Synthesis
This extensive example illustrates the use of VfeTools, VfsTools, VisTools and Base modules to build a full scale finite element solver for generating superelements using component mode synthesis. Superelements built in this fashion can be used in multi-body simulation systems to model flexible components. VisTools and Base modules are utilized as data structures to hold a complete description of the finite element model and the solution procedure to be applied to the model. VfeTools and VfsTools modules are used to implement the element formulations, material models and numerical analysis procedures used to create the superelement mass and stiffness matrices. The example includes support for solid, shell, beam and rigid type elements even though only shell and rigid elements are explicitly created.
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.
Figure 20-1, Stress Contour on Deflected Shape of Mode 5
#include <stdlib.h>
#include <math.h>
#include "sam/base/base.h"
#include "sam/vis/vis.h"
#include "sam/msh/msh.h"
#include "sam/vfs/vfs.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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 */
msh_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 = NULL, *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 = 0, 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;
}
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* create MapMesh object */
mapmesh = msh_MapMeshBegin();
msh_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++) {
msh_MapMeshSetPoint(mapmesh, i + 1, coords[i]);
}
/* define patch connectivity */
msh_MapMeshSetParami(mapmesh, VIS_MESH_MAXI, 2);
msh_MapMeshSetPatch(mapmesh, 1, VIS_SHAPEQUAD, 4, 2, 0, pix);
msh_MapMeshSetPatchParami(mapmesh, 1, VIS_MESH_SHAPE, VIS_SHAPEQUAD);
msh_MapMeshSetPatchParami(mapmesh, 1, MAPMESH_NUMELEMI, nelemcurve);
msh_MapMeshSetPatchParami(mapmesh, 1, MAPMESH_NUMELEMJ, nelemwidth);
/* set assoc of 100 at z=0 and 200 at z=-width */
/* use these to identify RBE connections */
msh_MapMeshSetPatchAssoc(mapmesh, VIS_MISCID, 1, SYS_EDGE, 1, 100);
msh_MapMeshSetPatchAssoc(mapmesh, VIS_MISCID, 1, SYS_EDGE, 3, 200);
/* generate shell mesh */
msh_MapMeshGenerate(mapmesh, connect);
/* delete MapMesh objects */
msh_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 = NULL;
Vint nprops;
vis_MProp* mprop = NULL;
Vint mtype, mid;
vfe_LinMat* linmat = NULL;
vfe_MatlFun* matlfun = NULL;
Vint ntypes, mattype[MPROP_MAX];
Vint densityflag, youngflag, poissonflag;
Vdouble density = 0, young = 0, poisson = 0, 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 = 0, *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*)vut_MemoryRealloc(loc, nedofs * sizeof(Vint));
tag = (Vint*)vut_MemoryRealloc(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 = 0, 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 = 0;
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];
}
}
}
14.16. Example 21, Inertia Relief and Linear Stress Analysis
This extensive example illustrates the use of VfeTools, VfsTools, VisTools and Base modules to compute body forces using inertia relief and perform linear stress analysis. The example is designed to support the case in which a part in a multi body dynamic simulation is subjected to a time history of forces and moments at its attachment points. It is desired to distribute the forces and moments at the attachment points at one or more times in the load history to corresponding sets of nodes on associated faces of the finite element model. These distributed loads are then balanced by additional inertial forces which are subsequently applied to the finite element model to compute stress.
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.
Figure 21-1a, Effective Stress Contour on Deflected Shape of Load 1
Figure 21-1b, Effective Stress Contour on Deflected Shape of Load 2
#include <stdlib.h>
#include <math.h>
#include "sam/base/base.h"
#include "sam/vis/vis.h"
#include "sam/msh/msh.h"
#include "sam/vfs/vfs.h"
#include "sam/vfe/vfe.h"
#include "sam/base/license.h"
#include "sam/CEETRON_SAM_license.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 */
msh_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.}}};
vsy_LicenseValidate(CEETRON_SAM_LICENSE);
/* Model Generation */
connect = vis_ConnectBegin();
vis_ConnectPre(connect, SYS_DOUBLE);
vis_ConnectDef(connect, 0, 0);
gridfun = vis_GridFunBegin();
vis_ConnectGridFun(connect, gridfun);
mapmesh = msh_MapMeshBegin();
msh_MapMeshDef(mapmesh, 8, 1);
for (i = 1; i <= 8; i++) {
msh_MapMeshSetPoint(mapmesh, i, p[i - 1]);
}
msh_MapMeshSetPatch(mapmesh, 1, VIS_SHAPEHEX, 2, 0, 0, pix);
msh_MapMeshSetParami(mapmesh, VIS_MESH_MAXI, 3);
msh_MapMeshSetPatchParami(mapmesh, 1, VIS_MESH_SHAPE, VIS_SHAPETET);
msh_MapMeshSetPatchParami(mapmesh, 1, MAPMESH_NUMELEMI, 10);
msh_MapMeshSetPatchParami(mapmesh, 1, MAPMESH_NUMELEMJ, 3);
msh_MapMeshSetPatchParami(mapmesh, 1, MAPMESH_NUMELEMK, 3);
msh_MapMeshSetPatchAssoc(mapmesh, VIS_MISCID1, 1, SYS_FACE, 5, 1);
msh_MapMeshSetPatchAssoc(mapmesh, VIS_MISCID2, 1, SYS_FACE, 6, 2);
msh_MapMeshGenerate(mapmesh, connect);
vis_ConnectKernel(connect, 0);
msh_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 = 0, n1 = 0, n2 = 0, n3 = 0, 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] = {0}, 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 = 0, young = 0, poisson = 0, 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 = 0, 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 */
void* temp = malloc(nattach * sizeof(Vdouble(*)[6]));
ir->cdist = (Vdouble(**)[6])temp;
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 < numdepdofs; ++k) {
vfe_MPCPrimCoeff(mpcdist, k + 1, xdist, coeff, &c);
for (i = numdepdofs; i < numdofs; ++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");
}
