t3.c

Go to the documentation of this file.

/*
%M This modules contains the T3 element sub-class methods and definitions
%a Joao Luiz Elias Campos.
%d September 1st, 1998.
%r $Id: t3.c,v 1.4 2004/06/24 18:32:35 joaoluiz Exp $
%w (C) COPYRIGHT 1995-1996, Eduardo Nobre Lages.
   (C) COPYRIGHT 1997-1999, Joao Luiz Elias Campos.
   All Rights Reserved
   Duplication of this program or any part thereof without the express
   written consent of the author is prohibited.
*
*  Modified:  28-Apr-2005    Alexandre A. Del Savio
*    Foram substituídas todas as alocações dinâmicas feitas com malloc 
*    por calloc.
*
*  Modified:  26-Jan-2005    Alexandre A. Del Savio
*    Criadas as funções T3Jacobian e T3Volume. E na funcão T3Init 
*    passou-se a fazer referência para estas funções, que
*    anteriormente eram nulas.
*
*  Modified:  Agosto-06  André Luis Muller
*    Criada a função _T3ElementCenter que retorna as coordenadas do centro de
*    um elemento do tipo T3.
*    Adicionado a função T3StressStrain o cálculo da poro pressão e o 
*    cálculo da tensão efetiva.
*
*  Modified:  Setembro-06 André Luis Muller
*    Modificação na função T3InterForce para possibilitar
*    a condição de tensões iniciais.
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função T3WriteGaussResults para possibilitar
*    a escolha das respostas a serem escritas no arquivo de resultados.
*/

/*
** ------------------------------------------------------------------------
** Global variables and symbols:
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "rio.h"
#include "load.h"
#include "elm.h"
#include "node.h"
#include "material.h"
#include "load.h"
#include "alg.h"


/*
** ------------------------------------------------------------------------
** Local variables and symbols:
*/

/*
%T T3 element data definition
*/
typedef struct _t3elm {
 int      matid;             /* Material identification          */
 int      tckid;             /* Element thickness identification */
 int      NumNodes;          /* Number of element nodes          */
 int      NumTensComp;       /* Number of tension components     */
 int      inc[3];            /* Element incidence                */
 double   istr[4];           /* Element initial stress           */
 double   effdef;            /* Effective deformation            */
 sTensor *iPress;            /* Element internal pressure        */
 double prof;                /* Element profile                  */
} sT3Data;


/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _T3GetCoord       ( sElement *, sCoord [] );
static double _T3Area           ( sCoord [] );
static double _T3MinHeight      ( sCoord [] );
static void   _T3GetDisplacement( sElement *, double *, double [] );
static void   _T3BMatrix        ( sElement *, double [3][6] );
static void   _T3ElementCenter  ( sElement *, sCoord [] );
static double _T3PoroPressure( sElement *elm, double lambda );

static double GForce = 9.85;

/*
%F This function gets the element coordinates.
%i Element descriptor.
%o Element coordinares.
*/
static void _T3GetCoord( sElement *elm, sCoord coord[3] )
{
 sT3Data *data = 0L;
 int      i;

/* Element element data
 */
 data = (sT3Data *)elm->data;

/* Get element coordinates
 */
 for( i = 0; i < 3; i++ ) {
  coord[i].x = elm->nodes[data->inc[i]-1].coord.x;
  coord[i].y = elm->nodes[data->inc[i]-1].coord.y;
  coord[i].z = 0.0;
 }

} /* End of _T3GetCoord */


/*
%F This function gets the element center.
%i Element descriptor.
%o Element center.
*/
static void _T3ElementCenter( sElement *elm, sCoord coord[1] )
{
 sT3Data *data = 0L;
 int      i, node;
 
 coord[0].x = 0.0;
 coord[0].y = 0.0;
 coord[0].z = 0.0;

 /* Element element data
 */
 data = (sT3Data *)elm->data;
 
 for( i = 0; i < data->NumNodes; i++ ) {
  node = data->inc[i]-1;
  coord[0].x += NodeVector[node].coord.x;
  coord[0].y += NodeVector[node].coord.y;
 }
 /* Compute element center
 */
 coord[0].x /= data->NumNodes;
 coord[0].y /= data->NumNodes;
} /* End of _T3ElementCenter */

/*
%F This function computes the element area
%i Element coordinates.
%o Element area.
*/
static double _T3Area( sCoord coord[3] )
{
 double area = 0.0;

/* Computes element area
 */
 area = ((coord[1].x - coord[0].x) * (coord[2].y - coord[0].y) -
         (coord[1].y - coord[0].y) * (coord[2].x - coord[0].x)) / 2.0;

 return( area );

} /* End of _T3Area */


/*
%F This function computes the minimum height for the element.
%i Element coordinate.
%o Minimum height
*/
static double _T3MinHeight( sCoord coord[3] )
{
 int    i;
 double area;
 double bcorr;
 double bmax;
 double hmin;

/* Initializa variables
 */
 bmax = bcorr = hmin = area = 0.0;

/* Loop over the element edges to compute the maximum edge
 */
 for( i = 0; i < 3; i++ ) {
  bcorr = sqrt( ( coord[(i+1)%3].x - coord[i].x ) * ( coord[(i+1)%3].x - coord[i].x ) +
                ( coord[(i+1)%3].y - coord[i].y ) * ( coord[(i+1)%3].y - coord[i].y ) );

  if( bcorr > bmax )
   bmax = bcorr;
 }

/* Computes the element area
 */
 area = _T3Area( coord );

/* Computes the minimum height
 */
 hmin = (2.0 * area) / bmax;

 return( hmin );

} /* End of _T3MinHight */


/*
%F This function get the element nodal displacement from the global 
   vector.
*/
static void _T3GetDisplacement( sElement *elm, double *U, double ue[6] )
{
 sT3Data *data = 0L;
 int      i;

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Fill element displacement vector
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  ue[2*i]   = U[2*(data->inc[i]-1)];
  ue[2*i+1] = U[2*(data->inc[i]-1)+1];
 }

} /* End of _T3GetDisplacement */


/*
%F This function computes the stress x strain matrix for the T3 element.
*/
static void _T3BMatrix( sElement *elm, double bm[3][6] )
{
 double detjac;
 double jac[2][2], ijac[2][2];
 sCoord coord[3];

/* Get element coordinate
 */
 _T3GetCoord( elm, coord );

/* Compute the element jacobian matrix
 */
 jac[0][0] = coord[1].x - coord[0].x;
 jac[0][1] = coord[1].y - coord[0].y;
 jac[1][0] = coord[2].x - coord[0].x;
 jac[1][1] = coord[2].y - coord[0].y;

/* Compute the matrix determinant
 */ 
 detjac = (jac[0][0] * jac[1][1]) - (jac[0][1] * jac[1][0]);

/* Compute the inverse
 */
 ijac[0][0] =   jac[1][1] / detjac;
 ijac[1][0] = - jac[1][0] / detjac;
 ijac[0][1] = - jac[0][1] / detjac;
 ijac[1][1] =   jac[0][0] / detjac;

/* Compute the B matrix coefficientes
 */
 bm[0][0] = bm[2][1] = - ijac[0][0] - ijac[0][1];
 bm[1][1] = bm[2][0] = - ijac[1][0] - ijac[1][1];
 bm[0][2] = bm[2][3] =   ijac[0][0];
 bm[1][3] = bm[2][2] =   ijac[1][0];
 bm[0][4] = bm[2][5] =   ijac[0][1];
 bm[1][5] = bm[2][4] =   ijac[1][1];

 bm[0][1] = bm[0][3] = bm[0][5] =
 bm[1][0] = bm[1][2] = bm[1][4] = 0.0;

} /* End of _T3BMatrix */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   T3New( int, int, int, int, sElement **, sElement **, sNode * );
static void   T3Free                  ( sElement * );
static void   T3Read                  ( sElement * );
static int    T3ReadInitStress        ( sElement * );
static void   T3ReadProfile           ( sElement * );
static void   T3MassMatrix            ( sElement *, double * );

static double T3RigidCoeff            ( sElement * );
static void   T3Connect               ( sElement *, int * );
static void   T3NumNodes              ( sElement *, int * );
static void   T3Gravity               ( sElement *, double *, double *, 
                                                    double * );
static void   T3AssVector             ( sElement *, double *, double * );
static void   T3StressStrain          ( sElement *, double, double *, double *,
                                                    sTensor *, sTensor * );
static void   T3TimeStep              ( sElement *, double * );
static void   T3InterForce            ( sElement *, sTensor *, double * );
static void   T3WriteStress           ( sElement *, FILE *, double *,
                                                            double * );
static void   T3WriteNodalResult      ( sElement *, FILE *, FILE * );
static void   T3WriteGaussResult      ( sElement *, FILE *, FILE * );
static void   T3WriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   T3UpdateStress          ( sElement *, double, double *,
                                                    sTensor * );
static void   T3PercForces            ( sElement *, double * );
static void   T3SetPressure           ( sElement *, double );
static void   T3SetInitStress         ( sElement *, sTensor * );
static void   T3ViscoForce            ( sElement *, double, sTensor *,
                                        double * );
static void   T3Jacobian              ( sElement *, double [3][3], 
                                        double [3][3] );
static void   T3Volume                ( sElement *, double * );

static void T3KMatrix( sElement *, double [24][24] );                                                                         
static void T3GetDof( sElement *, int [24], int * );
/*
%F This method allocates memory for a T3 element and fills its data with 
   the specific data.
%i Element label (number)
%o Element descriptor.
*/
static void T3New( int label, int matid, int intord, int tckid, 
                   sElement **elm, sElement **elist, sNode *nodes )
{
 sT3Data *data = 0L;

/* Get memory for the element descriptor
 */
 (*elm) = (sElement *)calloc(1, sizeof(sElement));

/* Get memory for the T3 element data
 */
 data = (sT3Data *)calloc(1, sizeof(sT3Data));

/* Fill T3 element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 3;
 data->NumTensComp = 4;
 data->effdef      = 0.0;
 data->istr[0]     = 0.0;
 data->istr[1]     = 0.0;
 data->istr[2]     = 0.0;
 data->iPress      = 0L;
 data->prof        = 0.0;

/* Fill element descriptor
 */
 (*elm)->type    = T3;
 (*elm)->rezone  = NONE;
 (*elm)->display = DSP_NO;
 (*elm)->curve   = 0;
 (*elm)->label   = label;
 (*elm)->data    = (void *)data;
 (*elm)->nodes   = nodes;

/* Add to the element list
 */
 elist[label-1] = (*elm);

} /* End of T3New */


/*
%F This method frees the T3 element data for a given element.
%i Element descriptor.
*/
static void T3Free( sElement *elm )
{
 sT3Data *data = 0L;

/* Get T3 element data
 */
 data = (sT3Data *)elm->data;

/* Release internal pressure
 */
 if( data->iPress != 0L )
  free( data->iPress );

/* Release allocated memory
 */
 free( data );

/* Reset element data
 */
 elm->data = 0L;

} /* End of T3Free */


/*
%F This method reads the T3 element information.
%i Element descriptor.
*/
static void T3Read( sElement *elm )
{
 sT3Data *data = 0L;
 int      i1, i2, i3;

/* Get the element data
 */
 data = (sT3Data *)elm->data;

/* Read the element incidence
 */
 fscanf( nf, "%d %d %d", &i1, &i2, &i3 );

/* Fill element information
 */
 data->inc[0] = i1;
 data->inc[1] = i2;
 data->inc[2] = i3;

} /* End of T3Read */


/*
%F This method reads the T3 initial stress
*/
static int T3ReadInitStress( sElement *elm )
{
 int      numpg = 1;
 double   sxx, syy, sxy ,szz, sxz, syz;
 sT3Data *data = 0L;

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Read the number of integration points and check to see if is compatible 
 * whith the element type.
 
 /*fscanf( nf, "%d", &numpg );*/
 if( numpg != 1 ) {
  printf( "\n\nNumber of gauss points must be equal to one.\n\n" );
  return 0;
 }
 

/* Read the initial stress information
 */
 fscanf( nf, "%lf %lf %lf %lf %lf %lf", &sxx, &syy, &szz, &sxy, &sxz, &syz);
 data->istr[0] = sxx;
 data->istr[1] = syy;
 data->istr[2] = sxy;

 return 1;
} /* End of T3ReadInitStress */

/*
%F This method reads the T3 profile
*/
static void T3ReadProfile( sElement *elm )
{
 double prof;
 sT3Data *data = 0L;

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Read the element profile
 */
 fscanf( nf, "%lf ", &prof);
 data->prof = prof;

} /* End of T3ReadProfile */


/*
%F This method computes the mass matrix for the T3 element.
%i Element descriptor.
%o Element mass matrix.
*/
static void T3MassMatrix( sElement *elm, double *mass )
{
 sT3Data *data = 0L;
 sCoord   coord[3];
 double   area, ncont, dens;
 int      i, matid;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Compute element area
 */
 area = _T3Area( coord );

/* Get material density
 */
 data  = (sT3Data *)elm->data;
 matid = data->matid;
 MatDensity( MatList[matid-1], &dens );

/* Calculate the element node contribution for the element mass
 */
 ncont = (area * dens) / 3.0;

/* Compute the mass matrix
 */
 for( i = 0; i < 6; i++ )
  mass[i] = ncont;
 
} /* End of T3MassMatrix */


/*
%F This functions computes the element rigidity coefficient.
%i Element descriptor.
%o Element rigidity coefficient.
*/
static double T3RigidCoeff( sElement *elm )
{
 sCoord coord[3];
 double area = 0.0;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Compute the element area
 */
 area = _T3Area( coord );

 return( area );

} /* End of T3RigidCoeff */


/*
%F This function computes the equivalent load applied to the element 
   face.
*/
static void T3Load( sElement *elm, eLoadType ltype, int key,
                    int noi, int noj, int nok, int *nol,
                    double *q1x, double *q1y, double *q1z,
                    double *q2x, double *q2y, double *q2z,
                    double *q3x, double *q3y, double *q3z,
                    double *q4x, double *q4y, double *q4z )
{
 double l, aux, teta;
 double v1x = (*q1x);
 double v1y = (*q1y);
 double v2x = (*q2x);
 double v2y = (*q2y);

/* Check for the local system
 */ 
 if( key ) {
  teta = atan2((elm->nodes[noj-1].coord.y - elm->nodes[noi-1].coord.y),
               (elm->nodes[noj-1].coord.x - elm->nodes[noi-1].coord.x));

  aux = (v1x * cos(teta)) - (v1y * sin(teta));
  v1y = (v1x * sin(teta)) + (v1y * cos(teta));
  v1x = aux;

  aux = (v2x * cos(teta)) - (v2y * sin(teta));
  v2y = (v2x * sin(teta)) + (v2y * cos(teta));
  v2x = aux;
 }

/* Compute the element edge lenght
 */
 l = sqrt(((elm->nodes[noi-1].coord.x - elm->nodes[noj-1].coord.x)  *
           (elm->nodes[noi-1].coord.x - elm->nodes[noj-1].coord.x)) +
          ((elm->nodes[noi-1].coord.y - elm->nodes[noj-1].coord.y)  *
           (elm->nodes[noi-1].coord.y - elm->nodes[noj-1].coord.y)));

/* Compute the equivalent forces
 */
 if( ltype == UNIFORM ) {
  v1x *= l / 2.0;
  v1y *= l / 2.0;

  (*q1x) = v1x;
  (*q1y) = v1y;

  (*q2x) = v1x;
  (*q2y) = v1y;
 }
 else {
  (*q1x) = ((2.0 * v1x * l) / 6.0) + ((v2x * l) / 6.0);
  (*q1y) = ((2.0 * v1y * l) / 6.0) + ((v2y * l) / 6.0);
 
  (*q2x) = ((2.0 * v2x * l) / 6.0) + ((v1x * l) / 6.0);
  (*q2y) = ((2.0 * v2y * l) / 6.0) + ((v1y * l) / 6.0);
 }

} /* End of T3Load */


/*
%F
*/
static void T3Connect( sElement *elm, int *conn )
{
 int      i;
 sT3Data *data = 0L; 

 data = (sT3Data *)elm->data;

 for( i = 0; i < data->NumNodes; i++ )
  conn[i] = data->inc[i];

} /* End of T3Connect */


/*
%F
*/
static void T3NumNodes( sElement *elm, int *nnodes )
{
 sT3Data *data = (sT3Data *)elm->data;

 (*nnodes) = data->NumNodes;

} /* End of T3NumNodes */


/*
%F
*/
static void T3Gravity( sElement *elm, double *qx, double *qy, double *qz )
{
 sT3Data *data = 0;
 sCoord   coord[3];
 double   area, mass, gamma;
 int      matid;

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Compute element area
 */
 area = _T3Area( coord );

/* Get material id, and the material density property
 */
 matid = data->matid;
 MatDensity( MatList[matid-1], &gamma );
 
/* Compute the gravitational load
 */
 mass = (gamma * area) / 3.0;

 (*qx) = mass * GravForce.gx;
 (*qy) = mass * GravForce.gy;
 (*qz) = mass * GravForce.gz;

} /* End of T3Gravity */


/*
%F
*/
static void T3AssVector( sElement *elm, double *GMatrix, double *matrix )
{
 int      i;
 int      conn[3];
 sT3Data *data;

/* Get the element data
 */
 data = (sT3Data *)elm->data;

/* Get element connectivity
 */
 ElmConnect( elm, conn );

/* Assemble matrix
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  GMatrix[2*(conn[i]-1)]   += matrix[2*i];
  GMatrix[2*(conn[i]-1)+1] += matrix[2*i+1];
 }
 
} /* End of T3AssVector */


/*
%F
*/


static void T3StressStrain( sElement *elm, double dt, double *U, double *yield, 
                                           sTensor *stre, sTensor *stra )
{
 sT3Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     cm[6][6];
 double     bm[3][6];
 double     ue[6];
 double     def[3], str[5];
 double     nu;
 double     pf = 1.0;
 int        i, j;
 double     poropressure;
 double     lambda;
 
/* Initialize stress and strain variables.
 */
 stre = (sTensor *)memset( (void *)stre, 0, sizeof(sTensor) );
 if( stra != 0L )
  stra = (sTensor *)memset( (void *)stra, 0, sizeof(sTensor) );

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Get Material object
 */
 mat = MatList[data->matid-1];

/* Update the material parameter
 */
 MatUpdateParameter( mat, data->effdef );

/* Get material constitutive matrix
 */
 MatConstitutiveMatrix( mat, cm );

/* Get poisson coeff.
 */
 MatNuParameter( mat, &nu );

/* Get element stress x strain matrix
 */
 _T3BMatrix( elm, bm );

/* Get element displacement
 */
 _T3GetDisplacement( elm, U, ue );

/* Compute the deformations
 */
 for( i = 0; i < 3; i++ ) {
  def[i] = 0.0;
  for( j = 0; j < (2*data->NumNodes); j++ )
   def[i] += bm[i][j] * ue[j];
 }

/* Compute the element stress
 */
 for( i = 0; i < 3; i++ ) {
  str[i] = 0.0;
  for( j = 0; j < 3; j++ )
   str[i] += cm[i][j] * def[j];
  str[i] += data->istr[i];
 }
 str[3] = str[4] = 0.0;
 
 /* Get element lambda
 */
 lambda = MatList[data->matid-1]->Lambda;
 
 if(lambda>0.0){
 
  poropressure = _T3PoroPressure(elm,lambda);

  /* Compute efective stress
  */
  for( i = 0; i < 2; i++ ) {
   str[i] += poropressure;
  }

  /* Correct stress based on the material model
  */
  MatUpdateStress( mat, dt, &pf, &data->effdef, str, def );
 
  /* Compute total stress
  */
  for( i = 0; i < 2; i++ ) {
   str[i] -= poropressure;
  }
 }
 else{
  /* Correct stress based on the material model
  */
  MatUpdateStress( mat, dt, &pf, &data->effdef, str, def );
 }

/* Set stress values
 */
 stre->xx = str[0];
 stre->yy = str[1];
 stre->xy = str[2];
 stre->zz = str[3];
 stre->rf = str[4];

/* Compute off plane stress
 */
 ElementOffPlaneStress( nu, stre );

/* Set strain values
 */
 if( stra != 0L ) {
  stra->xx = def[0];
  stra->yy = def[1];
  stra->xy = def[2];
 }

/* Set yield parameter
 */
 if( yield != 0L )
  (*yield) = (pf > 0.0) ? 0.0 : 1.0;

} /* End of T3StressStrain */

static double _T3PoroPressure( sElement *elm, double lambda )
{
 sT3Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     poropressure;
 double     gamma;
 double     distance;
 
 /* Get element descriptor
 */
 data = (sT3Data *)elm->data;

 /* Get Material object
 */
 mat = MatList[data->matid-1];
 
 /* Get element gamma
  */
 MatDensity( mat, &gamma );
 
 /* Get element prof
  */
 distance = data->prof;
 
 /* Compute the poropressure
  */ 
 poropressure = lambda*gamma*GForce*distance*Config.loadfactor;
 
 return (poropressure);
} /* End of _T3PoroPressure */



/*
%F
*/
static void T3TimeStep( sElement *elm, double *dt )
{
 sCoord     coord[3];
 double     hmin, dtime;
 sT3Data   *data = 0L;
 sMaterial *mat  = 0L;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Compute the element minimum height
 */
 hmin = _T3MinHeight( coord );

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Get material object
 */
 mat = MatList[data->matid-1];

/* Compute the time step according the material type
 */
 MatTimeStep( mat, &dtime );

/* Compute the time step
 */
 if( mat->type == KELVIN ) (*dt) = dtime;
 else                      (*dt) = hmin * dtime;

} /* End of T3TimeStep */


/*
%F This function computes the internal forces for the T3 element
*/
static void T3InterForce( sElement *elm, sTensor *stress, double *intforce )
{
 int      i, j;
 double   coeff = 0.0;
 double   bm[3][6], str[3];
 sT3Data *data = 0L;

/* Get element descriptor
 */ 
 data = (sT3Data *)elm->data;

/* Compute element rigidity coefficient
 */ 
 coeff = ElmRigidCoeff( elm );

/* Compute element B matrix
 */
 _T3BMatrix( elm, bm );

/* Check for the initial stress
 */
 if( stress == 0L ) {
  for( i = 0; i < (2*data->NumNodes); i++ ) {
   intforce[i] = 0.0;
   for( j = 0; j < 3; j++ ) 
    intforce[i] += bm[j][i] * data->istr[j];
   intforce[i] *= coeff;
  }
  return;
 }

/* Get element stress
 */
 str[0] = stress->xx;
 str[1] = stress->yy;
 str[2] = stress->xy;

/* Compute element internal forces
 */
 for( i = 0; i < (2*data->NumNodes); i++ ) {
  intforce[i] = 0.0;
  for( j = 0; j < 3; j++ ) 
   intforce[i] += bm[j][i] * str[j];
  intforce[i] *= coeff;
 }
 for( i = 0; i < (2*data->NumNodes); i++ ) 
  intforce[i] *= -1.0;

/* Compute element internal pressure
 */
 if( data->iPress != 0L ) {
  for( i = 0; i < (2*data->NumNodes); i++ ) {
   intforce[i] += coeff * bm[0][i] * data->iPress->xx;
   intforce[i] += coeff * bm[1][i] * data->iPress->yy;
  }
 }

} /* End of T3InterForce */


/*
%F
*/
static void T3WriteStress( sElement *elm, FILE *out, double *U, double *V )
{
 sTensor  stress;
 sTensor  strain;
 sTensor  strd = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
 double   yield = 0.0;
 double   iso;
 
/* Check to see if the element was rezoned
 */
 if( elm->rezone == DONE ) {
  stress.xx = stress.xy = stress.xz = 
              stress.yy = stress.yz = 
                          stress.zz = 0.0;
  strd.xx = strd.xy = strd.xz = 
            strd.yy = strd.yz = 
                      strd.zz = 0.0;
  strain.xx = strain.xy = strain.xz = 
              strain.yy = strain.yz = 
                          strain.zz = 0.0;
  yield = 0.0;
 }
 else {
 /* Compute the element stress and strain
  */
  T3StressStrain( elm, 0.0, U, &yield, &stress, &strain );
 }

/* Compute desviatory stress
 */
 iso = (stress.xx + stress.yy + stress.zz)/3.0;
 strd.xx = stress.xx - iso;
 strd.yy = stress.yy - iso;

/* Write element stress, strain and yield parameter
 */
 fwrite( &stress, sizeof(sTensor), 1, out );
 fwrite( &strain, sizeof(sTensor), 1, out );
 fwrite( &strd, sizeof(sTensor), 1, out );
 fwrite( &yield, sizeof(double), 1, out );

} /* End of T3WriteStress */


/*
%F
*/
static void T3WriteNodalResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  sts, sta;
 sPTensor psts, psta;
 sTensor  strd;
 double   yield;
 int      i;
 double   iso;

/* Read results from temporary file
 */
 fread( &sts, sizeof(sTensor), 1, tmp );
 fread( &sta, sizeof(sTensor), 1, tmp );
 fread( &strd, sizeof(sTensor), 1, tmp );
 fread( &yield, sizeof(double), 1, tmp );

/* Computes the principal tensors
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );
 
 /* Computes the isotropic stress
 */
 iso = (sts.xx + sts.yy + sts.zz)/3.0;

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 for( i = 0; i < 3; i++ ) {
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t%+10.5e\t", sts.xx,
                                                        sts.yy,
                                                        sts.zz,
                                                        sts.xy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psts.dir1,
                                               psts.dir2,
                                               psts.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", strd.xx,
                                               strd.yy,
                                               iso );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", sta.xx,
                                               sta.yy,
                                               sta.xy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psta.dir1,
                                               psta.dir2,
                                               psta.dir3 );
  fprintf( out, "%+10.5e\n", yield );
 }

} /* End of T3WriteNodalResult */


/*
%F
*/
static void T3WriteGaussResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  sts, sta;
 sPTensor psts, psta, epsts, dpsts;
 sTensor  strd, estr;
 double   yield;
 double   iso, piso;
 double   lambda;
 double   poropressure = 0.0;
 sMaterial *mat  = 0L;
 sT3Data   *data = 0L;
 int      j;

/* Read results from temporary file
 */
 fread( &sts, sizeof(sTensor), 1, tmp );
 fread( &sta, sizeof(sTensor), 1, tmp );
 fread( &strd, sizeof(sTensor), 1, tmp );
 fread( &yield, sizeof(double), 1, tmp );

/* Computes the principal stress
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );
 
/* Computes the isotropic stress
 */
 iso = (sts.xx + sts.yy + sts.zz)/3.0;
 
/* Compute desviatory principal stress
 */
 piso = (psts.dir1 + psts.dir2 + psts.dir3)/3.0;
 
 dpsts.dir1 = psts.dir1 - piso;
 dpsts.dir3 = psts.dir3 - piso;
 
 /* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Get Material object
 */
 mat = MatList[data->matid-1];
 
 /* Get element lambda
 */
 lambda = MatList[data->matid-1]->Lambda;
 
 if(lambda>0.0){
  poropressure = _T3PoroPressure(elm,lambda);
  /* Compute efective stress
  */
  estr.xx = sts.xx + poropressure;
  estr.yy = sts.yy + poropressure; 
 }
 else{
  estr.xx = sts.xx;
  estr.yy = sts.yy;
 }
 estr.xy = sts.xy;
 
/* Computes the principal efetive stress
 */
 PrincipalTensor( &estr, &epsts );  
 
/* Write results on output file
 */
 
 fprintf( out, "%d\t%d\n", elm->label, 1 );
 
 for( j = 0; j < Config.numgaussresults; j++){
   if((strcmp( Config.gaussresults[j], "Yield" ) == 0))
    fprintf( out, "%+10.5e\t", yield );
   else if((strcmp( Config.gaussresults[j], "Rf" ) == 0))
    fprintf( out, "%+10.5e\t", sts.rf );
   else if((strcmp( Config.gaussresults[j], "Sxx" ) == 0))
    fprintf( out, "%+10.5e\t", sts.xx );
   else if((strcmp( Config.gaussresults[j], "Syy" ) == 0))
    fprintf( out, "%+10.5e\t", sts.yy );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", sts.xy );
   else if((strcmp( Config.gaussresults[j], "Pf" ) == 0))
    fprintf( out, "%+10.5e\t", poropressure );
   else if((strcmp( Config.gaussresults[j], "SExx" ) == 0))
    fprintf( out, "%+10.5e\t", estr.xx );
   else if((strcmp( Config.gaussresults[j], "SEyy" ) == 0))
    fprintf( out, "%+10.5e\t", estr.yy );
   else if((strcmp( Config.gaussresults[j], "Siso" ) == 0))
    fprintf( out, "%+10.5e\t", iso );
   else if((strcmp( Config.gaussresults[j], "SDxx" ) == 0))
    fprintf( out, "%+10.5e\t", strd.xx );
   else if((strcmp( Config.gaussresults[j], "SDyy" ) == 0))
    fprintf( out, "%+10.5e\t", strd.yy );
   else if((strcmp( Config.gaussresults[j], "S1" ) == 0))
    fprintf( out, "%+10.5e\t", psts.dir1 );
   else if((strcmp( Config.gaussresults[j], "S3" ) == 0))
    fprintf( out, "%+10.5e\t", psts.dir3 );
   else if((strcmp( Config.gaussresults[j], "SE1" ) == 0))
    fprintf( out, "%+10.5e\t", epsts.dir1 );
   else if((strcmp( Config.gaussresults[j], "SE3" ) == 0))
    fprintf( out, "%+10.5e\t", epsts.dir3 );
   else if((strcmp( Config.gaussresults[j], "SD1" ) == 0))
    fprintf( out, "%+10.5e\t", dpsts.dir1 );
   else if((strcmp( Config.gaussresults[j], "SD3" ) == 0))
    fprintf( out, "%+10.5e\t", dpsts.dir3 );
   else if((strcmp( Config.gaussresults[j], "Exx" ) == 0))
    fprintf( out, "%+10.5e\t", sta.xx );
   else if((strcmp( Config.gaussresults[j], "Eyy" ) == 0))
    fprintf( out, "%+10.5e\t", sta.yy );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", sta.xy );
   else if((strcmp( Config.gaussresults[j], "E1" ) == 0))
    fprintf( out, "%+10.5e\t", psta.dir1 );
   else if((strcmp( Config.gaussresults[j], "E3" ) == 0))
    fprintf( out, "%+10.5e\t", psta.dir3 );
   else if((strcmp( Config.gaussresults[j], "Ev" ) == 0))
    fprintf( out, "%+10.5e\t", (psta.dir1+psta.dir2+psta.dir3) );
   else if((strcmp( Config.gaussresults[j], "K0" ) == 0)){
    if(estr.yy != 0.0){
     fprintf( out, "%+10.5e\t", (estr.xx/estr.yy) );
    }
    else{
     fprintf( out, "%+10.5e\t", 0.0 );
    }
   }
 }                               
 
} /* End of T3WriteGaussResult */


/*
%F
*/
static void T3WriteGaussVectorResult( sElement *elm, int version, FILE *out,
                                      FILE *tmp )
{
 sPTensor psts, psta;
 sTensor  sts, sta;
 sTensor  strd;
 double   yield;

/* Read results from temporary file
 */
 fread( &sts, sizeof(sTensor), 1, tmp );
 fread( &sta, sizeof(sTensor), 1, tmp );
 fread( &strd, sizeof(sTensor), 1, tmp );
 fread( &yield, sizeof(double), 1, tmp );

/* Computes the principal stress
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 1 );

 if( version > 1 ) {
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir1, psta.dir1 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir2, psta.dir2 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir3, psta.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1x, psta.cos1x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2x, psta.cos2x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3x, psta.cos3x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1y, psta.cos1y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2y, psta.cos2y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3y, psta.cos3y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1z, psta.cos1z );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2z, psta.cos2z );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3z, psta.cos3z );
 }
 else {
  if( psts.dir1 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psts.dir1*psts.cos1x),
                                                (psts.dir1*psts.cos1y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psts.dir1 < 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psts.dir1*psts.cos1x),
                                                (psts.dir1*psts.cos1y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psts.dir3 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psts.dir3*psts.cos1y),
                                                (psts.dir3*psts.cos3y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psts.dir3 < 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psts.dir3*psts.cos3x),
                                                (psts.dir3*psts.cos3y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psta.dir1 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psta.dir1*psta.cos1x),
                                                (psta.dir1*psta.cos1y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );
 
  if( psta.dir1 < 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psta.dir1*psta.cos1x),
                                                (psta.dir1*psta.cos1y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psta.dir3 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psta.dir3*psta.cos3x),
                                                (psta.dir3*psta.cos3y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

  if( psta.dir3 < 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psta.dir3*psta.cos3x),
                                                (psta.dir3*psta.cos3y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );
 }

} /* End of T3WriteGaussVectorResult */


/*
%F This function update the element stress to reflect the effects of 
   the larges displacements.
*/
static void T3UpdateStress( sElement *elm, double dtime, double *V, 
                                                         sTensor *stre )
{
 sT3Data *data = 0L;
 double   bm[3][6];
 double   ve[6];
 double   sxx, syy, rot;
 int      i;

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Get element stress x strain matrix
 */
 _T3BMatrix( elm, bm );

/* Get element displacement increment
 */
 _T3GetDisplacement( elm, V, ve );

/* Compute the displacement increment
 */
 for( i = 0; i < 6; i++ )
  ve[i] *= dtime;

/* Compute rotation factor
 */
 rot = 0.0;
 for( i = 0; i < 6; i++ )
  rot += bm[2][i] * ve[i] * pow( -1.0, (double)(i+1) );
 rot *= 0.5;

/* Compute the element stress
 */
 sxx = stre->xx;
 syy = stre->yy;

 stre->xx -= 2.0 * stre->xy * rot;
 stre->yy += 2.0 * stre->xy * rot;
 stre->xy += (sxx - syy) * rot;

} /* End of T3UpdateStress */


/*
%F
*/
static void T3PercForces( sElement *elm, double *pforce )
{
 sCoord     coord[3];
 sT3Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     bm[3][6], phi[3];
 sTensor    grad;
 int        i, k;

/* Initialize percolation forces vector
 */
 for( i = 0; i < 6; i++ )
  pforce[i] = 0.0;

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Get material object
 */
 mat = MatList[data->matid-1];

/* Get element pressure
 */
 for( i = 0; i < data->NumNodes; i++ )
  phi[i] = elm->nodes[data->inc[i]-1].dof.psi + coord[i].y;

/* Get element stress x strain matrix
 */
 _T3BMatrix( elm, bm );

/* Compute gradient
 */
 grad.xx = grad.yy = 0.0;
 for( k = 0; k < 3; k++ ) {
  grad.xx -= bm[0][2*k] * phi[k];
  grad.yy -= bm[1][2*k+1] * phi[k];
 }

/* Compute element internal forces
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  pforce[2*i]   = grad.xx;
  pforce[2*i+1] = grad.yy;
 }

} /* End of T3PercForces */


/*
%F
*/
static void T3SetPressure( sElement *elm, double pot )
{
 sCoord   coord[3];
 sT3Data *data = 0L;
 double   elev;
 double   ncg = 0.33333333333333333;

/* Get element coordinates
 */
 _T3GetCoord( elm, coord );

/* Compute element CG elevation
 */
 elev = ncg * (coord[0].y + coord[1].y + coord[2].y);
 
/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Compute and set the pressure for the element CG
 */
 data->iPress = (sTensor *)calloc(1, sizeof(sTensor));
 ElementInitTensor( data->iPress );
 data->iPress->xx = 
 data->iPress->yy = pot - elev;

} /* End of T3SetPressure */

/*
%F
*/
static void T3SetInitStress( sElement *elm, sTensor *istress )
{
 sT3Data *data = 0L;

/* Get element data
 */
 data = (sT3Data *)elm->data;

/* Set element initial stress
 */
 data->istr[0] = istress->xx;
 data->istr[1] = istress->yy;
 data->istr[2] = istress->xy;
 
} /* End of T3SetInitStress */

/*
%F
*/
static void T3ViscoForce( sElement *elm, double timeStep, sTensor *stress,
                          double *vforce )
{
 int        i, j;
 double     coeff;
 double     bm[3][6], cm[6][6];
 double     stav[3], strv[4];
 sTensor    vsta;
 sT3Data   *data = 0L;
 sMaterial *mat  = 0L;

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Get material object
 */
 mat = MatList[data->matid-1];

/* Initialize element visco forces vector
 */
 for( i = 0; i < (2*data->NumNodes); i++ )
  vforce[i] = 0.0;

/* Check to see if the material is a visco one
 */
 if( !MaterialIsVisco( mat ) )
  return;

/* Compute element rigidity coefficient
 */ 
 coeff = ElmRigidCoeff( elm );

/* Compute element B matrix
 */
 _T3BMatrix( elm, bm );

/* Get material constitutive matrix
 */
 MatConstitutiveMatrix( mat, cm );

/* Compute visco stress
 */
 MatViscoStrain( mat, timeStep, stress, &vsta );

/* Get viscos strain components
 */
 stav[0] = vsta.xx;
 stav[1] = vsta.yy;
 stav[2] = vsta.xy;

/* Compute viscos stress
 */
 for( i = 0; i < 3; i++ ) {
  strv[i] = 0.0;
  for( j = 0; j < 3; j++ )
   strv[i] += cm[i][j] * stav[j];
 }
 strv[3] = 0.0;

/* Compute element visco forces
 */
 for( i = 0; i < (2*data->NumNodes); i++ ) {
  vforce[i] = 0.0;
  for( j = 0; j < 3; j++ )
   vforce[i] += bm[j][i] * strv[j];
  vforce[i] *= coeff;
 } 

} /* End of T3ViscoForce */

static void T3Jacobian( sElement *elm, double jac[3][3], double ijac[3][3] )
{
  double detjac;
  sCoord coord[3];

  /* Get element coordinate
   */
  _T3GetCoord( elm, coord );

  /* Compute the element jacobian matrix
   */
  jac[0][0] = coord[1].x - coord[0].x;
  jac[0][1] = coord[1].y - coord[0].y;
  jac[1][0] = coord[2].x - coord[0].x;
  jac[1][1] = coord[2].y - coord[0].y;
  jac[0][2] = 0.0;
  jac[1][2] = 0.0;
  jac[2][0] = 0.0;
  jac[2][1] = 0.0;
  jac[2][2] = 0.0;

  /* Compute the matrix determinant
   */ 
  detjac = (jac[0][0] * jac[1][1]) - (jac[0][1] * jac[1][0]);

  /* Compute the inverse
   */
  ijac[0][0] =   jac[1][1] / detjac;
  ijac[1][0] = - jac[1][0] / detjac;
  ijac[0][1] = - jac[0][1] / detjac;
  ijac[1][1] =   jac[0][0] / detjac;
  ijac[0][2] = 0.0;
  ijac[1][2] = 0.0;
  ijac[2][0] = 0.0;
  ijac[2][1] = 0.0;
  ijac[2][2] = 0.0;

} /* End of T3Jacobian */


static void T3Volume( sElement *elm, double *v )
{
  sCoord coord[2];

  /* Get element coordinates
   */
  _T3GetCoord(elm, coord);

  /* Compute element area
   */
  (*v) = _T3Area(coord);

} /* End of T3Volume */

static void T3KMatrix( sElement *elm, double k[24][24] )
{         
 int    i, j, l;           
 sCoord  coord[3];                  
 double detjac;            
 double C[6][6];          
 double B[3][6];          
 double BTC[6][3];        
 double ke[6][6]; 
 sT3Data  *data = 0L;
 sMaterial  *mat  = 0L;       

 /* Inicia matriz k = 0
 */
 for( i = 0; i < 24; i++ ) {
  for( j = 0; j < 24; j++ ) {
   k[i][j] = 0.0;
  }
 }
 
 /* Pega descrição do elemento
 */
 data = (sT3Data *)elm->data;
 
/* Pega material
 */
 mat = MatList[data->matid-1];
 
 /* Pega as coordenadas nodais
 */
 _T3GetCoord(elm, coord);
 
 /* Pega matriz constitutiva [C]*/
 MatConstitutiveMatrix( mat, C );

 /* Calcula matriz B*/
 _T3BMatrix(elm, B );
 
 /* Calcula a área do elemento */
 detjac = _T3Area( coord);
 
 /*Calcula [B] transposto * [C] */
 for( i = 0; i < 6; i++ ) {
  for( j = 0; j < 3; j++)  {
   BTC[i][j] =0.0;
   for(l=0; l<3; l++){
    BTC[i][j] += B[l][i]*C[l][j];
   }
  }
 }
 /* Calcula [B] transposto * [C] * [B] * |jac|*/
 for( i = 0; i < 6; i++ ) {
  for( j = 0; j < 6; j++)  {
   ke[i][j] =0.0;
   for(l=0; l<3; l++){
    ke[i][j] += BTC[i][l]*B[l][j]*detjac;
   }
  }
 }
 /* Calcula matriz de rigidez do elemento*/
 for( i = 0; i < 6; i++ ) {
  for( j = 0; j < 6; j++)  {
   k[i][j] += ke[i][j];
  }
 }
 
} /* T3KMatrix */

static void T3GetDof( sElement *elm, int u[24], int *index )
{
 sT3Data *data = 0L;
 int        i;

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;

/* Fill element dof
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  u[2*i]   = 2*(data->inc[i]-1);
  u[2*i+1] = 2*(data->inc[i]-1)+1;
 }
 
 *index = 6;

} /* End of T3GetDof */

static void T3GetInc( sElement *elm, int inc[8], int *index )
{
 sT3Data *data = 0L;
 int        i;

/* Get element descriptor
 */
 data = (sT3Data *)elm->data;
 
 for( i = 0; i < data->NumNodes; i++ ) {
  inc[i] = data->inc[i]-1;
 }
 
 *index = 3;
} /* End of T3GetInc */

/*
** ------------------------------------------------------------------------
** Public functions:
*/

/*
%F This function initializes the sub-class "T3".
*/
void T3Init( void );
void T3Init( void )
{
/* Initialize T3 element sub-class
 */
 ElmClass[T3].new               = T3New;
 ElmClass[T3].free              = T3Free;
 ElmClass[T3].read              = T3Read;
 ElmClass[T3].readinitstr       = T3ReadInitStress;
 ElmClass[T3].readprofile       = T3ReadProfile;
 ElmClass[T3].mass              = T3MassMatrix;
 ElmClass[T3].rigidcoeff        = T3RigidCoeff;
 ElmClass[T3].load              = T3Load;
 ElmClass[T3].connect           = T3Connect;
 ElmClass[T3].numnodes          = T3NumNodes;
 ElmClass[T3].gravity           = T3Gravity;
 ElmClass[T3].assvector         = T3AssVector;
 ElmClass[T3].strstrain         = T3StressStrain;
 ElmClass[T3].timestep          = T3TimeStep;
 ElmClass[T3].intforce          = T3InterForce;
 ElmClass[T3].writestr          = T3WriteStress;
 ElmClass[T3].writendlresult    = T3WriteNodalResult;
 ElmClass[T3].writegauresult    = T3WriteGaussResult;
 ElmClass[T3].writegauvecresult = T3WriteGaussVectorResult;
 ElmClass[T3].updatestress      = T3UpdateStress;
 ElmClass[T3].percforce         = T3PercForces;
 ElmClass[T3].setpressure       = T3SetPressure;
 ElmClass[T3].setinitstress     = T3SetInitStress;
 ElmClass[T3].viscoforce        = T3ViscoForce;
 ElmClass[T3].jacobian          = T3Jacobian;
 ElmClass[T3].volume            = T3Volume;
 ElmClass[T3].KMatrix           = T3KMatrix;
 ElmClass[T3].GetDof            = T3GetDof;
 ElmClass[T3].GetInc            = T3GetInc;
 
} /* End of T3Init */


/* =======================================================  End of File  */

---