Solid2D 2D solid element. Solid3D 3D solid element. Shell3D 3D monocoque or composite shell element. Mem3D 3D membrane element. Beam3D 3D beam and shell stiffener element. Truss3D 3D truss or cable element. Inter2D 2D interface element. Inter3D 3D interface element. Bulk Volumetric mass or heat capacitance ConMass Concentrated or degree of freedom mass ConElas Grounded or degree of freedom elastic spring Gap Point to point gap element Spring Spring element MPC Multi-point constraint
LinMat Linear isotropic, orthotropic or anisotropic materials PlasMat Isotropic, rate independent, large strain plasticity HyperMat Mooney-Rivlin and Ogden materials ShellProp Shell element wall properties BeamProp Beam element section properties InterProp Interface and gap element properties MatlFun Define function pointers to material model
ElemTran Element transformations for rotated coordinate systems. NodeGeom Node accelerations, rigid body modes, inertia relief. Corot Corotational large rotation utilities
Each instance of a module is termed an object. Users can generate as many objects as needed. The internal state of an object may be changed at any time to reflect, for example, the element topology it is meant to accept. As an example, users may change the current topology accepted by a Solid2D object from a 4 node quadrilateral to a 6 node triangle.
VfeTools can be used in a straightforward manner by the host application. The general methodology for generating element stiffness and mass matrices can be described by the following steps:
vfe_Solid2D *solid2d; solid2d = vfe_Solid2DBegin (); vfe_Solid2DSetTopology (solid2d,SYS_SHAPEQUAD,2,0);
vfe_LinMat *linmat; Vdouble young = 10000000.0; Vdouble poisson = 0.3; Vdouble props[2]; linmat = vfe_LinMatBegin (); vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC); props[0] = young; props[1] = poisson; vfe_LinMatSetElasProp (linmat,props); vfe_LinMatSetDensity (linmat,0.0000133);
vfe_MatlFun *matlfun; matlfun = vfe_MatlFunBegin (); vfe_LinMatMatlFun (linmat,matlfun); vfe_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun);
Vint nedofs; Vint loc[VFE_MAXDOF]; Vint tag[VFE_MAXDOF]; vfe_Solid2DNumDof (solid2d,SYS_ANALYSIS_STRUCTURAL,&nedofs); vfe_Solid2DDofMap (solid2d,SYS_ANALYSIS_STRUCTURAL,loc,tag);
Vint loc[8] = { 1,1, 2,2, 3,3, 4,4 }; Vint tag[8] = { SYS_DOF_TX,SYS_DOF_TY, SYS_DOF_TX,SYS_DOF_TY, SYS_DOF_TX,SYS_DOF_TY, SYS_DOF_TX,SYS_DOF_TY };
-- -- | k11 | | k12 k22 | | k13 k23 k33 | | k14 k24 k34 k44 | -- --Element force and diagonal mass are output as vectors. The host application is responsible for the memory required to hold the output stiffness, mass, force, etc. Note that the defined constant VFE_MAXNODE is equal to the maximum number of nodes used by any element formulation in VfeTools (about 30).
Vdouble x[VFE_MAXNODE][3]; Vdouble stiff[VFE_MAXDOF*(VFE_MAXDOF)/2], mdiag[VFE_MAXDOF]; /* localize coordinates, x, from host application */ . . /* generate element stiffness and diagonal mass */ vfe_Solid2DStiff (solid2d,x,stiff); vfe_Solid2DMassDiag (solid2d,x mdiag);
A variety of element topologies are supported. Basic topologies include lines, triangles, quadrilaterals, tetrahedra, pyramids, wedges and hexahedra. VfeTools supports linear and quadratic element orders in general, however cubic 2D and 3D solid elements are supported. There are two distinct forms for each element topology:
Module_Name(number_of_nodes,technology)The number_of_nodes and technology fields are optional.
The module name defines the element type and dimensionality. For example 2D solid elements are implemented in the Solid2D module. The number of nodes serves as a convenient and unique notation for the element shape and polynomial order within any module. The technologies, such as "isoparametric", are abbreviated as follows: isoparametric, ISOP; uniform reduced, URED. For example, the notation for the 2D solid Serendipity parabolic quadrilateral using uniform reduced technology would be as follows:
Solid2D(8,URED)If an asterisk, *, appears in place of the number of nodes or technology, it is meant to apply for all numbers of nodes or all technologies.
Coordinates
Solid2D(*) x, y Solid3D(*) x, y, z
Degrees of freedom
Solid2D(*) tx, ty Solid3D(*) tx, ty, tz
Stresses and strains
Solid2D(*) sxx,syy,szz,sxy,syz,szx exx,eyy,ezz,gxy,gyz,gzx Solid3D(*) sxx,syy,szz,sxy,syz,szx exx,eyy,ezz,gxy,gyz,gzxIn all 2D analysis syz, szx, gyz and gzx are zero. In 2D plane stress szz is zero, in 2D plane strain ezz is zero.
Integration and polynomial order
Normal Reduced Solid2D(3) 1 1 linear triangle Solid2D(4) 2x2 1 linear quadrilateral Solid2D(6) 3 1 parabolic triangle Solid2D(8) 3x3 2x2 parabolic Serendipity quadrilateral Solid2D(9) 3x3 2x2 parabolic Lagrange quadrilateral Normal Reduced Solid3D(4) 1 1 linear tetrahedron Solid3D(5) 5 1 linear pyramid Solid3D(6) 1x2 1 linear wedge Solid3D(8) 2x2x2 1 linear hexahedron Solid3D(10) 4 1 parabolic tetrahedron Solid3D(13) 3x3x3 5 parabolic Serendipity pyramid Solid3D(14) 3x3x3 5 parabolic Lagrange pyramid Solid3D(15) 3x3 3x2 parabolic Serendipity wedge Solid3D(18) 3x3 3x2 parabolic Lagrange wedge Solid3D(20) 3x3x3 2x2x2 parabolic Serendipity hexahedron Solid3D(27) 3x3x3 2x2x2 parabolic Lagrange hexahedronGaussian quadrature is used for all quadrilateral and hexahedral elements and in the prismatic direction of wedge elements.
Coordinates
Shell3D(*) x, y, z
Degrees of freedom
Shell3D(*) tx, ty, tz, rx, ry, rz
Stress resultants and midsurface strains and curvatures in shell local coordinate system
Shell3D(*) Nxx,Nyy,Nxy,Mxx,Myy,Mxy,Qxz,Qyz exx,eyy,gxy,kxx,kyy,kxy,gxz,gyzThe number of integration points used by each shell element type appears below. The integration orders are for the shell surface integration. Gaussian quadrature is used for all quadrilateral elements.
Normal Reduced Shell3D(3) 1 1 linear triangle Shell3D(4) 2x2 1 linear quadrilateral Shell3D(6) 3 1 parabolic triangle Shell3D(8) 3x3 2x2 parabolic Serendipity quadrilateral Shell3D(9) 3x3 2x2 parabolic Lagrange quadrilateral
Coordinates
Beam3D(*) x, y, z
Degrees of freedom
Beam3D(*) tx, ty, tz, rx, ry, rz
Stress resultants and reference axis strains, curvatures and twist in beam local coordinate system
Beam3D(*) Nxx,Myy,Mzz,Txx,Qxy,Qzx exx,kyy,kzz,twist,gxy,gzx
The number of integration points used by each beam element type appears below. The integration orders are for the beam axis integration. Gaussian quadrature is used for all elements.
Normal Reduced Beam3D(2) 1 1 linear line Beam3D(3) 2 1 parabolic line
The value for the normalized tip displacement in the direction of the load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 3.0E-05 for extension, 0.1081 for in-plane shear, 0.4321 for out-of-plane shear and 0.03208 for twist.
Tip loading direction (a) Rectangular elements Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) Extension 0.995 0.997 0.999 0.999 In-plane shear 0.903 0.983 0.987 0.994 Out-of-plane shear 0.980 0.992 0.994 0.994 Twist 0.946 0.954 0.949 0.949 NASTRAN QUAD4 QUAD8 Extension 0.995 0.999 In-plane shear 0.904 0.987 Out-of-plane shear 0.986 0.991 Twist 0.941 0.950 (b) Trapezoidal elements Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) Extension 0.988 0.994 0.999 0.999 In-plane shear 0.043 0.891 0.965 0.989 Out-of-plane shear 0.971 0.979 0.995 0.994 Twist 0.927 0.947 0.950 0.952 NASTRAN QUAD4 QUAD8 Extension 0.996 0.999 In-plane shear 0.071 0.946 Out-of-plane shear 0.968 0.998 Twist 0.951 0.943 (c) Parallelogram elements Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) Extension 1.011 0.993 0.999 0.999 In-plane shear 0.017 0.888 0.987 0.987 Out-of-plane shear 0.980 0.986 0.994 0.995 Twist 0.855 0.956 0.949 0.949 NASTRAN QUAD4 QUAD8 Extension 0.996 0.999 In-plane shear 0.080 0.995 Out-of-plane shear 0.977 0.985 Twist 0.945 0.965
The value for the normalized tip displacement in direction of load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.08734 for the in-plane load and 0.5022 for the out-of-plane load.
Tip loading direction Normalized tip displacement in direction of load Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) In-plane 0.833 0.998 1.021 1.011 Out-of-plane 0.936 0.938 0.957 0.957 NASTRAN QUAD4 QUAD8 In-plane 0.833 1.007 Out-of-plane 0.951 0.971 Solid3D(8) Solid3D(10) Solid3D(15) Solid3D(20) Solid3D(27) In-plane 0.822 0.887 0.947 1.016 1.021 Out-of-plane 0.790 0.868 0.908 0.929 0.932 NASTRAN HEXA(8) HEX20(R) In-plane 0.880 1.006 Out-of-plane 0.849 0.959 Beam3D(2) Beam3D(3) In-plane 1.022 1.042 Out-of-plane 0.991 1.003
The value for the normalized tip displacement in direction of load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.005424 for the in-plane load and 0.001754 for the out-of-plane load.
Tip loading direction Normalized tip displacement in direction of load Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) In-plane 0.987 0.984 0.997 1.001 Out-of-plane 0.981 0.987 1.003 1.003 NASTRAN QUAD4 QUAD8 In-plane 0.985 0.998 Out-of-plane 0.993 0.998
The value for the normalized displacement at the center (intersection of two symmetry boundaries) of the plate is used as a measure of solution accuracy. We have assumed the solutions for the normalization of results to be the following for shell elements, multiply by 1.0e-06 for solid elements since the thickness for solid elements is 100 times that used for shell elements.
Boundary Aspect Displacement at center of plate supports ratio b/a Uniform pressure Concentrated force Simple 1.0 4.062 11.60 Simple 5.0 12.97 16.96 Clamped 1.0 1.26 5.60 Clamped 5.0 2.56 7.23
Simple support uniform load, aspect ratio = 1.0 Normalized lateral deflection at center Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 0.987 1.093 0.985 1.003 4 0.997 1.029 1.001 1.003 6 0.999 1.014 1.001 1.001 8 0.999 1.008 1.000 1.000 NASTRAN QUAD4 QUAD8 2 0.981 0.927 4 1.004 0.996 6 1.003 0.999 8 1.002 1.000 Simple support uniform load, aspect ratio = 5.0 Normalized lateral deflection at center Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 0.890 0.595 1.279 1.245 4 0.941 0.779 0.986 0.995 6 0.970 0.844 0.990 0.996 8 0.982 0.890 0.994 0.998 NASTRAN QUAD4 QUAD8 2 1.052 1.223 4 0.991 1.003 6 0.997 1.000 8 0.998 1.000 X Clamped support concentrated load, aspect ratio = 1.0 Normalized lateral deflection at center Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 0.881 0.001 1.116 1.116 4 0.971 0.879 0.902 1.011 6 0.988 0.966 0.990 1.004 8 0.994 0.987 0.984 1.003 NASTRAN QUAD4 QUAD8 2 0.934 1.076 4 1.010 0.969 6 1.012 0.992 8 1.010 0.997 Clamped support concentrated load, aspect ratio = 5.0 Normalized lateral deflection at center Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 0.319 0.001 0.423 0.423 4 0.841 0.401 0.942 0.950 6 0.922 0.556 0.957 0.988 8 0.954 0.672 0.969 0.994 NASTRAN QUAD4 QUAD8 2 0.519 0.542 4 0.863 0.754 6 0.940 0.932 8 0.972 0.975
The value for the midside vertical displacement is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.3024 .
Normalized vertical deflection at midpoint of free edge Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 1.387 0.872 0.832 0.863 4 1.048 0.863 0.955 0.961 6 1.017 0.970 0.987 0.988 8 1.010 0.992 0.993 0.993 10 1.009 0.997 0.994 0.994 NASTRAN QUAD4 QUAD8 2 1.376 1.021 4 1.050 0.984 6 1.018 1.002 8 1.008 0.997 10 1.004 0.996
The radial displacement under the load is used as a measure of solution accuracy. We have assumed the solution for the normalization of results to be 0.094 .
Normalized radial deflection at load point Shell3D(4) Shell3D(6) Shell3D(8) Shell3D(9) 2 0.909 0.910 1.333 1.774 4 0.999 0.860 0.646 0.889 8 0.990 1.012 0.961 0.989 12 0.990 1.019 0.987 0.991 NASTRAN QUAD4 QUAD8 2 0.972 0.025 4 1.024 0.121 8 1.005 0.823 12 0.998 0.992 ABAQUS S4 STRI65 S8R5 S9R5 2 0.444 0.229 0.796 0.669 4 0.925 0.916 0.951 0.944 8 0.974 0.967 0.982 0.982 12 0.980 0.974 0.982 0.982 Normalized radial deflection at load point Solid3D(8) Solid3D(20) Solid3D(27) 4 0.034 0.022 0.367 8 0.702 0.571 0.980 12 0.951 0.919 1.007 NASTRAN HEXA(8) HEX20(R) 4 0.039 0.162 8 0.730 0.776 12 0.955 0.972Table of Contents
F_reaction = F_applied + F_inertia
F_reaction = K u
F_inertia = - M a
The collection of all nodal and elemental unknowns is denoted as the degree of freedom map or "dof map" of the element. Typical degrees of freedom in structural mechanics are the Cartesian displacements or translations. Shell and beam elements require nodal rotational degrees of freedom in addition to the translations.
In heat transfer simulations, typical degrees of freedom can be the temperature and/or its gradient. As the analysis type and the element technology are selected along with the element topology, the number and type of degrees of freedom associated with the element are determined. Note that the number of degrees of freedom and dof map depends on the discipline (e.g., stress analysis, thermal analysis transfer, etc.), the technology utilized (e.g., isoparametric, enhanced strains, etc.), and the element topology and order (e.g. quadratic triangle, linear quadrilateral).
Every element in VfeTools provides a method to query the element for the number of degrees of freedom, as well as the dof map. The number of degrees of freedom, nedofs, is queried as follows:
vfe_Solid2DNumDof (solid2d,SYS_ANALYSIS_STRUCTURAL,&nedofs);while the dof map, vectors loc and tag (described below), are queried as follows:
vfe_Solid2DDofMap (solid2d,SYS_ANALYSIS_STRUCTURAL,loc,tag);The vector loc is dimensioned nedofs long, and contains for every degree of freedom in the element, the node index into the element connectivity associated with the degree of freedom. A node index of 1 indicates that the degree of freedom is associated with the first node in the element connection list. If the degree of freedom happens to be internal to the element, the node index is zero.
The array tag is also dimensioned nedofs, and describes the type of each of the degrees of freedom. The degree of freedom types are specified by a defined constant as follows:
Primitive materials are characterized by a material response which excludes any element geometry. Examples include stress - strain relations, heat flux - temperature gradient relations, etc. Within VfeTools all modules which implement primitive materials are given names with the string "Mat" appended. As an example, the LinMat module implements linear isotropic, orthotropic and anisotropic materials. Primitive materials include:
Element materials are material responses which include the effects of element properties. Examples include layered composites which include element thickness and offset information as well as primitive material definitions for each layer. An element material need not necessarily include a primitive material. For example an element material for an elastic spring may simply contain a single scalar spring constant. Within VfeTools all modules which implement element materials are given names with the string "Prop" appended. As an example, the ShellProp module implements layered composites for shell element formulations. Element materials include:
In the finite element method, the displacement field is usually explicitly assumed to be continuous between elements. However, under such an assumption, the computed stress field will generally be discontinuous between elements to some degree depending upon the accuracy of the finite element discretization. This fact can be exploited to compute an estimate of the error of a given discretization and suggest a new discretization size to obtain an improved solution.
An effective method for estimating error is to compute the strain energy error on an element by element basis and compare it to the total strain energy. This computation requires an integration over the element of the difference between the exact and computed nodal stresses. The computed nodal stresses are determined by the element given the displacement field. An approximation to the exact nodal stresses may be obtained by averaging the computed stresses between adjacent elements where the stresses are expected to be continuous. This approximation process is termed "stress recovery" and the method suggested above for recovering stress is one of many possible methods.
A number of element modules in VfeTools provide methods to compute the element energy and energy error. Supported element types include solid, shell, beam and truss element types. For structural analysis the StrsAdapt function is provided. This function computes the element strain energy, setot and strain energy error, seerr. In addition, useful quantities such as the characteristic element length, h, effective polynomial order, p, and dimension of the element, d, are returned.
vfe_Solid2DStrsAdapt (solid2d,x,u,strss,setot,seerr,&h,&p,&d);The element characteristic length is the edge length of an element with the same area or volume assuming equal length, straight edges and maximum minimum angles between edges. The element polynomial order is approximately the polynomial order of displacement interpolation along an element edge. The element dimension is the number of local element dimensions spanned by the element connectivity. For example, 3D shells and 2D solids have an element dimensionality of 2, 3D solids have an element dimensionality of 3.
A similar function is available in heat transfer analysis for computing total heat energy and heat energy error given the nodal temperatures and recovered nodal flux.
vfe_Solid2DHFlxAdapt (solid2d,x,u,hflxs,hetot,heerr,&h,&p,&d);To obtain an estimate of the global error the elemental energies and energy errors are computed using the functions above and summed for all applicable elements. Let E represent the total strain energy and e the energy error obtained in this manner, then the total solution error err may be obtained as:
err = sqrt(e/E)A solution is considered converged when the solution error is less than some predefined tolerance, tol, usually 5%. It is interesting to note that is has been shown that if seerr is equal for all elements, then the finite element model using a given number of elements is the optimal one. This fact suggests a adaptation process in which a new characteristic element length is computed for each element to equilibrate the solution error and reduce it to the convergence tolerance.
First an estimate of the number of elements in the adapted mesh is determined which yields a solution error equal to the convergence tolerance. This is done by assuming that the mesh size will be uniformly scaled by a scaling factor, f. Since the energy error, seerr, in any element is of the order h**p where h and p are the element characteristic length and polynomial order, for a mesh characterized by p, the energy error in the new mesh scaled by f is related to the current energy error by
e_new = (f**p)*eDividing by E and solving for f yields
f = (tol/err)**(1/p)The number of elements in the new mesh is given by
N_new = N/(f**d)where N is the current number of elements and d is the element dimensionality in the mesh. Assuming that the element energy error is equilibrated in the new mesh, the target element energy error is given as
seerr_tar = tol*E/sqrt(N_new)Although the target element energy error is computed assuming a uniform scaling of the current mesh, it is more realistic and efficient to adjust the mesh size locally according to the current contribution of each element to the solution error. This suggests that the new characteristic element size for each element be computed as
h_new = (seerr_tar/seerr)**(1/p)*hNew element sizes computed in this fashion may be given to an automatic mesh generator to generate an entirely new mesh which will be refined or derefined in a very localized way.
To address performance issues from a strict operation count, cache size and processor architecture point of view, VfeTools modules avoid zero arithmetic.
VfeTools is basically platform independent and as a result no user defined C preprocessor directives are required to compile VfeTools on any supported platform. However for performance reasons certain platform specific compiler options may be necessary to notify the compiler that pointers are not "aliased", etc.
It is suggested that during the development cycle that the source be conditionally compiled with error checking by defining VKI_CHECK_ERRS as described in base library.
For example, on SGI systems, create a directory SGI under lib to hold the final devtools.a archive file. Then from the devtools/src/vfe directory compile using
creating .o files in the vfe directory. As mentioned above it is suggested to use the appropriate compiler option to specify that pointers are not aliased plus additional optimizations as follows:cc -c -O2 -I.. *.c
cc -c -O2 -OPT:ieee_arith=3 -OPT:roundoff=3 -OPT:alias=restrict *.c
cc -c -fast -xO3 -xrestrict *.c
The object files may be deleted after the devtools.a archive is successfully created.ar rs ../../lib/SGI/devtools.a *.o
To compile the base source, change directory to devtools/src/base and compile using
creating .o files in the base directory. To add these objects to the previously created devtools.a archive issuecc -c -O2 *.c
To compile the vis source, change directory to devtools/src/vis and compile usingar rs ../../lib/SGI/devtools.a *.o
creating .o files in the vis directory. To add these objects to the previously created devtools.a archive issuecc -c -O2 *.c
Again object files may deleted at this time. At this point the devtools.a archive contains all vfe, vis, and base objects. Place the devtools.a archive immediately before the graphics subsystem libraries in the load line. A devtools.a archive must be built for each computer platform using the methodology outlined above.ar rs ../../lib/SGI/devtools.a *.o
#include <stdio.h> #include "base/base.h" #include "vfe/vfe.h" static Vdouble x[3][3] = { {0.,0.,0.}, {1.,0.,0.}, {0.,1.,0.} }; /*---------------------------------------------------------------------- Generate Stiffness Matrix for 3 Node Triangle ----------------------------------------------------------------------*/ int main() { Vint i,j; Vdouble prop[2]; Vdouble kl[21]; vfe_Solid2D *solid2d; vfe_MatlFun *matlfun; vfe_LinMat *linmat; /* Create material object and set properties */ linmat = vfe_LinMatBegin (); vfe_LinMatDef (linmat,SYS_MAT_ISOTROPIC); prop[0] = 1.0; prop[1] = .3; vfe_LinMatSetElasProp (linmat,prop); /* Create MatlFun object and load material functions */ matlfun = vfe_MatlFunBegin (); vfe_LinMatMatlFun (linmat, matlfun); /* Create 2D solid element formulation */ solid2d = vfe_Solid2DBegin (); vfe_Solid2DSetTopology (solid2d,SYS_SHAPETRI,2,0); /* Register material functions */ vfe_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun); /* Form and print stiffness matrix */ vfe_Solid2DStiff (solid2d,x,kl); printf("Lower triangle of stiffness matrix\n"); for(i = 0; i < 6; i++) { for(j = 0; j <= i; j++) { printf("%12.4e ",kl[i*(i+1)/2+j]); } printf("\n"); } /* Delete objects */ vfe_Solid2DEnd (solid2d); vfe_MatlFunEnd (matlfun); vfe_LinMatEnd (linmat); return 0; }
The output of this example program appears below. Only the lower half of the matrix is printed.
Lower triangle of stiffness matrix 7.4176e-01 3.5714e-01 7.4176e-01 -5.4945e-01 -1.6484e-01 5.4945e-01 -1.9231e-01 -1.9231e-01 0.0000e+00 1.9231e-01 -1.9231e-01 -1.9231e-01 0.0000e+00 1.9231e-01 1.9231e-01 -1.6484e-01 -5.4945e-01 1.6484e-01 0.0000e+00 0.0000e+00 5.4945e-01
#include <stdio.h> #include "base/base.h" #include "vfe/vfe.h" static Vdouble x[3][3] = { {0.,0.,0.}, {1.,0.,0.}, {0.,1.,0.} }; /*---------------------------------------------------------------------- Generate Stiffness Matrix for 3 Node Triangle ----------------------------------------------------------------------*/ int main() { Vint i,j; Vdouble prop[2]; Vdouble kl[21]; vfe_Solid2D *solid2d; vfe_MatlFun *matlfun; vfe_LinMat *linmat; /* Create material object and set properties */ linmat = new vfe_LinMat (); linmat->Def (SYS_MAT_ISOTROPIC); prop[0] = 1.0; prop[1] = .3; linmat->SetElasProp (prop); /* Create MatlFun object and load material functions */ matlfun = new vfe_MatlFun (); linmat->MatlFun (matlfun); /* Create 2D solid element formulation */ solid2d = new vfe_Solid2D (); solid2d->SetTopology (SYS_SHAPETRI,2,0); /* Register material functions */ solid2d->SetObject (VFE_MATLFUN,matlfun); /* Form and print stiffness matrix */ solid2d->Stiff (x,kl); printf("Lower triangle of stiffness matrix\n"); for(i = 0; i < 6; ++i ) { for(j = 0; j <= i; ++j) { printf("%12.4e ",kl[i*(i+1)/2+j]); } printf("\n"); } /* Delete objects */ delete solid2d; delete matlfun; delete linmat; return 0; }Table of Contents
C----------------------------------------------------------------------- C Generate Stiffness Matrix for 3 Node Triangle C----------------------------------------------------------------------- program main include 'base/fortran/base.inc' include 'vfe/fortran/vfe.inc' integer i,j double precision prop(2) double precision kl(21) double precision solid2d double precision matlfun double precision linmat double precision x(3,3) data x / 0.,0.,0., 1.,0.,0., 0.,1.,0. / C Create material object and set properties call vfef_LinMatBegin (linmat) call vfef_LinMatDef (linmat,SYS_MAT_ISOTROPIC) prop(1) = 1.0 prop(2) = .3 call vfef_LinMatSetElasProp (linmat,prop) C Create MatlFun object and load material functions call vfef_MatlFunBegin (matlfun) call vfef_LinMatMatlFun (linmat, matlfun) C Create 2D solid element formulation call vfef_Solid2DBegin (solid2d) call vfef_Solid2DSetTopology (solid2d,SYS_SHAPETRI,2,0) C Register material functions call vfef_Solid2DSetObject (solid2d,VFE_MATLFUN,matlfun) C Form and print stiffness matrix call vfef_Solid2DStiff (solid2d,x,kl) write(6,*) "Lower triangle of stiffness matrix" do i = 0,5 write(6,1000) (kl(i*(i+1)/2+j+1),j=0,i) enddo C Delete objects call vfef_Solid2DEnd (solid2d) call vfef_MatlFunEnd (matlfun) call vfef_LinMatEnd (linmat) C 1000 format(1x,6e12.4) stop endTable of Contents
using System; using System.Runtime.InteropServices; using System.Reflection; using System.Text; using DevTools; public class intro1 { public static double [] x = { 0.0,0.0,0.0, 1.0,0.0,0.0, 0.0,1.0,0.0 }; /*---------------------------------------------------------------------- Generate Stiffness Matrix for 3 Node Triangle ----------------------------------------------------------------------*/ public static void Main() { int i = 0,j = 0; double [] prop = new double[2]; double [] kl = new double[21]; IntPtr solid2d; IntPtr matlfun; IntPtr linmat; /* Create material object and set properties */ linmat = vfe.LinMatBegin (); vfe.LinMatDef (linmat,vsy.SYS_MAT_ISOTROPIC); prop[0] = 1.0; prop[1] = 0.3; vfe.LinMatSetElasProp (linmat,prop); /* Create MatlFun object and load material functions */ matlfun = vfe.MatlFunBegin (); vfe.LinMatMatlFun (linmat, matlfun); /* Create 2D solid element formulation */ solid2d = vfe.Solid2DBegin (); vfe.Solid2DSetTopology (solid2d,vsy.SYS_SHAPETRI,2,0); /* Register material functions */ vfe.Solid2DSetObject (solid2d,vfe.VFE_MATLFUN,matlfun); /* Form and print stiffness matrix */ vfe.Solid2DStiff (solid2d,x,kl); Console.Write("Lower triangle of stiffness matrix\n"); for(i = 0; i < 6; i++) { for(j = 0; j <= i; j++) { Console.Write("{0} ",kl[i*(i+1)/2+j]); } Console.Write("\n"); } /* Delete objects */ vfe.Solid2DEnd (solid2d); vfe.MatlFunEnd (matlfun); vfe.LinMatEnd (linmat); } }