tetra4.c

Go to the documentation of this file.

/*
%M This modules contains the TETR4 element sub-class methods and 
   definitions
%a Joao Luiz Elias Campos.
%d February 20th, 2001.
%r $Id: tetra4.c,v 1.4 2004/06/24 18:32:35 joaoluiz Exp $
%w (C) COPYRIGHT 1995-1996, Eduardo Nobre Lages.
   (C) COPYRIGHT 1997-2001, 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:  27-Apr-2005    Alexandre A. Del Savio
*    Foram substituídas todas as alocações dinâmicas feitas com malloc 
*    por calloc.
*
*  Modified:  23-Fev-2006    Juan Pablo Ibañez
*    Modificada a funcão TETR4ReadInitStress que passa a retornar um valor int.
*
*  Modified:  Agosto-06  André Luis Muller
*    Criada a função _TETRA4ElementCenter que retorna as coordenadas do centro de
*    um elemento do tipo TETRA4.
*    Adicionado a função TETRA4StressStrain 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 TETRA4InterForce para possibilitar
*    a condição de tensões iniciais.
*
*    Modificação na função TETRA4New. Inicialização das tensões iniciais
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função TETRA4WriteGaussResults para possibilitar
*    a escolha das respostas a serem escritas no arquivo de resultados.
*
*/

/*
** ------------------------------------------------------------------------
** Global variables and symbols:
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.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 TETR4 element data definition
*/
typedef struct _tetr4elm
{
 int    matid;                /* Material identification          */
 int    tckid;                /* Element thickness identification */
 int    NumNodes;             /* Number of element nodes          */
 int    NumTensComp;          /* Number of tension components     */
 int    inc[4];               /* Element incidence                */
 double istr[6];              /* Element initial stress           */
 double effdef;               /* Effective deformation            */
 double prof;                 /* Element profile                  */
} sTetr4Data;

/*
%V Gauss point coordinates and its weights
*/
static double _GPoint  = 0.25000000;
static double _GWeight = 0.16666667;

static double GForce = 9.85;

/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _TETR4GetCoord       ( sElement *, sCoord [] );
static void   _TETR4ShapeFunc      ( double, double, double, double [] );
static void   _TETR4Jacobian       ( double, double, double, sCoord [],
                                     double *, double [3][3] );
static void   _TETR4BMatrix        ( double, double, double, sCoord [],
                                     double [6][12] );
static void   _TETR4DerivShp       ( double, double, double, sDerivRST [] );
static void   _TETR4GetDisplacement( sElement *, double *, double [] );
static double _TETR4MinHeight      ( sCoord * );
static double _TETR4Volume         ( sCoord * );
static void   _TETR4ElementCenter  ( sElement *elm, sCoord coord[1] );
static void   _TETR4FaceCenter( sCoord *coord0, sCoord *coord1, 
                                sCoord *coord2, sCoord *center );
                                
static double RigidCoeff( sElement *elm, int noi, int noj, int nok );
static double _Area( sCoord coord[3] );
static void _Normal( sCoord coord[3], sCoord *normal );  

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

/* Get element coordinates
 */
 coord[0].x = elm->nodes[noi-1].coord.x;
 coord[0].y = elm->nodes[noi-1].coord.y;
 coord[0].z = elm->nodes[noi-1].coord.z;
 coord[1].x = elm->nodes[noj-1].coord.x;
 coord[1].y = elm->nodes[noj-1].coord.y;
 coord[1].z = elm->nodes[noj-1].coord.z;
 coord[2].x = elm->nodes[nok-1].coord.x;
 coord[2].y = elm->nodes[nok-1].coord.y;
 coord[2].z = elm->nodes[nok-1].coord.z;

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

 return( area );

} /* End of RigidCoeff */

static double _Area( sCoord coord[3] )
{
 double area = 0.0;
 double u[3],v[3],w[3];

/* Computes element area
 */
 u[0] = coord[1].x - coord[0].x;
 u[1] = coord[1].y - coord[0].y; 
 u[2] = coord[1].z - coord[0].z;
 v[0] = coord[2].x - coord[0].x; 
 v[1] = coord[2].y - coord[0].y; 
 v[2] = coord[2].z - coord[0].z;

 w[0] = u[1] * v[2] - u[2] * v[1];
 w[1] = u[2] * v[0] - u[0] * v[2];
 w[2] = u[0] * v[1] - u[1] * v[0];
 area = sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]) / 2.0;
          
 return( area );

} /* End of _Area */



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

 /* Element element data
 */
 data = (sTetr4Data *)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;
  coord[0].z += NodeVector[node].coord.z;
 }
 /* Compute element center
 */
 coord[0].x /= data->NumNodes;
 coord[0].y /= data->NumNodes;
 coord[0].z /= data->NumNodes;
} /* End of _TETR4ElementCenter */

/*
%F This function computes the element volume.
*/
static double _TETR4Volume( sCoord *coord )
{
 double jac[3][3];
 double detjac, volume = 0.0;

/* Compute the jacobian matrix
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &detjac, jac );

/* Rigidity coefficient
 */
 volume = detjac * _GWeight;

 return volume;

} /* End of _TETR4Volume */

/* Calcula o centro de um triangulo */
static void _TETR4FaceCenter( sCoord *coord0, sCoord *coord1, 
                              sCoord *coord2, sCoord *center )
{
 sCoord med[2]; /* medianas do triângulo */
 double dx, dy, dz;

/* O ponto de interseção das três medianas de um triangulo é o
 * seu baricentro. O baricentro divide a mediana em dois segmentos.
 * O segmento que une o vértice ao baricentro vale o dobro do 
 * segmento que une o baricentro ao lado oposto deste vértice. 
 */

 med[0].x = coord0->x;
 med[0].y = coord0->y;
 med[0].z = coord0->z;
 med[1].x = ( coord1->x + coord2->x ) / 2.0;
 med[1].y = ( coord1->y + coord2->y ) / 2.0;
 med[1].z = ( coord1->z + coord2->z ) / 2.0;
 
 dx = ( med[1].x - med[0].x ) / 3.0;
 dy = ( med[1].y - med[0].y ) / 3.0;
 dz = ( med[1].z - med[0].z ) / 3.0;

 center->x = med[0].x + (2*dx);
 center->y = med[0].y + (2*dy);
 center->z = med[0].z + (2*dz);
 }

/*
%F This function computes the minimun height for the TETR4 element.
 */

#if 0
static double _TETR4MinHeight( sCoord *coord )
{
 double volume, bcorr, bmax, hmin;
 int    i;

/* Initialize variables
 */
 bmax = bcorr = hmin = volume = 0.0;

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

  if( bcorr > bmax )
   bmax = bcorr;
 }

/* Compute element volume
 */
 volume = _TETR4Volume( coord );

/* Compute the minimum height
 */
 hmin = (2.0 * volume) / bmax;

 return hmin;
} /* End of _TETR4MinHeight */

#else
static double _TETR4MinHeight( sCoord *coord )
{
 sCoord center;
 double hmin = 0.0;
 double h;
 int i;
 int first = 1;
 
 for( i = 0; i < 4; i++ ) 
 {
  _TETR4FaceCenter( &coord[i+1], &coord[(i+2)%4], &coord[(i+3)%4], &center );
  h = sqrt( ( (coord[i].x - center.x)*(coord[i].x - center.x) ) +
            ( (coord[i].y - center.y)*(coord[i].y - center.y) ) +
            ( (coord[i].z - center.z)*(coord[i].z - center.z) ) ); 
  if( first || ( h < hmin ) ) 
  {
   first = 0;
   hmin = h;
  }
 }
 return hmin;
}
#endif

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

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

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

} /* End of _TETR4GetDisplacement */


/*
%F This function computes the stress x strain matrix for the TETR4 element.
*/
static void _TETR4BMatrix( double r, double s, double t, sCoord coord[4], 
                                                         double bm[6][12] )
{
 double    detjac;
 double    jac[3][3], ijac[3][3];
 sDerivRST dshp[4];
 int       i;

/* Initialize B matrix
 */
 for( i = 0; i < 12; i++ )
  bm[0][i] = bm[1][i] = bm[2][i] = 
  bm[3][i] = bm[4][i] = bm[5][i] = 0.0;

/* Compute the element jacobian matrix
 */
 _TETR4Jacobian( r, s, t, coord, &detjac, jac );

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


/* Compute shape functions derivatives
 */
 _TETR4DerivShp( r, s, t, dshp );

/* Compute the B matrix coefficientes
 */
 for( i = 0; i < 4; i++ ) {
  bm[0][3*i]   = bm[3][3*i+1] = bm[4][3*i+2] = (ijac[0][0] * dshp[i].r) +
                                               (ijac[0][1] * dshp[i].s) + 
                                               (ijac[0][2] * dshp[i].t);
  bm[1][3*i+1] = bm[3][3*i]   = bm[5][3*i+2] = (ijac[1][0] * dshp[i].r) +
                                               (ijac[1][1] * dshp[i].s) + 
                                               (ijac[1][2] * dshp[i].t);
  bm[2][3*i+2] = bm[4][3*i]   = bm[5][3*i+1] = (ijac[2][0] * dshp[i].r) +
                                               (ijac[2][1] * dshp[i].s) + 
                                               (ijac[2][2] * dshp[i].t);
 }

} /* End of _TETR4BMatrix */


/*
%F This function computes the jacobian matrix for the quadrilateral 
   brick 8 element.
*/
static void _TETR4Jacobian( double r, double s, double t, sCoord coord[4], 
                            double *detjac, double jac[3][3] )
{
 sDerivRST dmap[4];
 int       i;

/* Compute the map derivatives
 */
 _TETR4DerivShp( r, s, t, dmap );

/* Compute the jacobian matrix
 */
 for( i = 0; i < 3; i++ )
  jac[i][0] = jac[i][1] = jac[i][2] = 0.0;

 for( i = 0; i < 4; i++ ) {
  jac[0][0] += dmap[i].r * coord[i].x;
  jac[0][1] += dmap[i].r * coord[i].y;
  jac[0][2] += dmap[i].r * coord[i].z;
  jac[1][0] += dmap[i].s * coord[i].x;
  jac[1][1] += dmap[i].s * coord[i].y;
  jac[1][2] += dmap[i].s * coord[i].z;
  jac[2][0] += dmap[i].t * coord[i].x;
  jac[2][1] += dmap[i].t * coord[i].y;
  jac[2][2] += dmap[i].t * coord[i].z;
 }

/* Determiant
 */
 (*detjac) = jac[0][0] * (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1]) -
             jac[0][1] * (jac[1][0] * jac[2][2] - jac[1][2] * jac[2][0]) +
             jac[0][2] * (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]);

} /* End of _TETR4Jacobian */


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

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

/* Get element coordinates
 */
 for( i = 0; i < 4; 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 = elm->nodes[data->inc[i]-1].coord.z;
 }
 
} /* End of _TETR4GetCoord */


/*
%F This function computes the shape functions for the quadrilateral 
   brick 8 element.
*/
static void _TETR4ShapeFunc( double r, double s, double t, double shape[4] )
{
/* Compute shape functions
 */
 shape[0] = r;
 shape[1] = s;
 shape[2] = t;
 shape[3] = 1.0 - r - s - t;

} /* End of _TETR4ShapeFunc */


/*
%F This function computes the shape functions derivatives for the 
   quadrilateral brick 8 element.
*/
static void _TETR4DerivShp( double r, double s, double t, sDerivRST dshp[8] )
{
/* Compute shape functions derivatives
 */
 dshp[0].r =  1.0;
 dshp[1].r =  0.0;
 dshp[2].r =  0.0;
 dshp[3].r = -1.0;

 dshp[0].s =  0.0;
 dshp[1].s =  1.0;
 dshp[2].s =  0.0;
 dshp[3].s = -1.0;

 dshp[0].t =  0.0;
 dshp[1].t =  0.0;
 dshp[2].t =  1.0;
 dshp[3].t = -1.0;

} /* End of _TETR4DerivShp */

static void _Normal( sCoord coord[3], sCoord *normal )
{
 double size;
 int first;
 int curr;
 int prev;
 int    i;

 normal->x = 0.0;
 normal->y = 0.0;
 normal->z = 0.0;
 
 first = 0;

 for( i = 2; i < 3; i++ )
 {
  curr = i;
  prev = i-1;
  
  normal->x += (coord[prev].y - coord[first].y) * 
                 (coord[curr].z - coord[first].z) -
                 (coord[prev].z - coord[first].z) *
                 (coord[curr].y - coord[first].y);
                 
  normal->y += (coord[prev].x - coord[first].x) * 
                 (coord[curr].z - coord[first].z) +
                 (coord[prev].z - coord[first].z) *
                 (coord[curr].x - coord[first].x);
                 
  normal->z += (coord[prev].x - coord[first].x) * 
                 (coord[curr].y - coord[first].y) -
                 (coord[prev].y - coord[first].y) *
                 (coord[curr].x - coord[first].x);
 }

 size = sqrt( (normal->x) * (normal->x) + (normal->y) *
              (normal->y) + (normal->z) * (normal->z) );

 if( size > 0.0 ){
  normal->x /= size;
  normal->y /= size;
  normal->z /= size;
 }
 
} /* end of _Normal */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   TETR4New( int, int, int, int, sElement **, sElement **, sNode * );
static void   TETR4Free                  ( sElement * );
static void   TETR4Read                  ( sElement * );
static int    TETR4ReadInitStress        ( sElement * );
static void   TETR4ReadProfile           ( sElement * );
static void   TETR4MassMatrix            ( sElement *, double * );
static void   TETR4Gravity               ( sElement *, double *, double *,
                                                                 double * );
static double TETR4RigidCoeff            ( sElement * );
static void   TETR4Connect               ( sElement *, int * );
static void   TETR4NumNodes              ( sElement *, int * );
static void   TETR4AssVector             ( sElement *, double *, double * );
static void   TETR4TimeStep              ( sElement *, double * );
static void   TETR4InterForce            ( sElement *, sTensor *, double * );
static void   TETR4StressStrain          ( sElement *, double, double *,
                                                       double *, sTensor *, 
                                                       sTensor * );
static void   TETR4WriteStress           ( sElement *, FILE *, double *, 
                                                               double * );
static void   TETR4WriteGaussResult      ( sElement *, FILE *, FILE * );
static void   TETR4WriteNodalResult      ( sElement *, FILE *, FILE * );
static void   TETR4WriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   TETR4UpdateStress          ( sElement *, double, double *,
                                                               sTensor * );
static void   TETR4PercForces            ( sElement *, double * );
static void   TETR4SetPressure           ( sElement *, double );
static void   TETR4SetInitStress         ( sElement *, sTensor * );
static void   TETR4ViscoForce            ( sElement *, double, sTensor *, 
                                           double * );                                                                                     
static void TETR4KMatrix( sElement *, double [24][24] );                                                                         
static void TETR4GetDof( sElement *, int [24], int * );

/*
%F This function computes the equivalent load applied to the element 
   face.
*/
static void TETR4Load( 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;
 double v1x = (*q1x);
 double v1y = (*q1y);
 double v1z = (*q1z);
 double v2x = (*q2x);
 double v2y = (*q2y);
 double v2z = (*q2z);
 double area;
 sCoord normal;
 sCoord coord[3]; 
 sTetr4Data *data = 0L;
 int i, node; 
 sCoord vec;
 double prod;
 
 *nol = -1;

/* Compute the equivalent forces
 */
 if( ltype == UNIFORM ) {
  /* 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)) +
          ((elm->nodes[noi-1].coord.z - elm->nodes[noj-1].coord.z)  *
           (elm->nodes[noi-1].coord.z - elm->nodes[noj-1].coord.z)));
           
  v1x *= l / 2.0;
  v1y *= l / 2.0;
  v1z *= l / 2.0;

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

  (*q2x) = v1x;
  (*q2y) = v1y;
  (*q2z) = v1z;
 }
 if( ltype == AREAUNIFORM ) {
        
  area = RigidCoeff(elm,noi,noj,nok);
  
  (*q1x) = 1.0/3.0*v1x*area;
  (*q1y) = 1.0/3.0*v1y*area;
  (*q1z) = 1.0/3.0*v1z*area;
  (*q2x) = 1.0/3.0*v1x*area;
  (*q2y) = 1.0/3.0*v1y*area;
  (*q2z) = 1.0/3.0*v1z*area;
  (*q3x) = 1.0/3.0*v1x*area;
  (*q3y) = 1.0/3.0*v1y*area;
  (*q3z) = 1.0/3.0*v1z*area;
 }
 if( ltype == AREAPRESSURE ) {
        
  area = RigidCoeff(elm,noi,noj,nok);
  
  coord[0].x = elm->nodes[noi-1].coord.x;
  coord[0].y = elm->nodes[noi-1].coord.y;
  coord[0].z = elm->nodes[noi-1].coord.z;
  coord[1].x = elm->nodes[noj-1].coord.x;
  coord[1].y = elm->nodes[noj-1].coord.y;
  coord[1].z = elm->nodes[noj-1].coord.z;
  coord[2].x = elm->nodes[nok-1].coord.x;
  coord[2].y = elm->nodes[nok-1].coord.y;
  coord[2].z = elm->nodes[nok-1].coord.z;
  
  _Normal( coord, &normal );
  
  /* para garantir a direção da carga*/
  
  data = (sTetr4Data *)elm->data;
  for( i = 0; i < data->NumNodes; i++ ) {
   if((data->inc[i]!=noi)&&(data->inc[i]!=noj)&&(data->inc[i]!=nok)){
    node = data->inc[i];
   }
  }
          
  vec.x = (elm->nodes[noi-1].coord.x - NodeVector[node-1].coord.x);
  vec.y = (elm->nodes[noi-1].coord.y - NodeVector[node-1].coord.y);
  vec.z = (elm->nodes[noi-1].coord.z - NodeVector[node-1].coord.z); 
  
  prod = vec.x * normal.x + vec.y * normal.y + vec.z * normal.z; 
  if( prod < 0.0){
   normal.x *= (-1.0);
   normal.y *= (-1.0);
   normal.z *= (-1.0);
  }
  
  (*q1x) = 1.0/3.0 * v1x * (-normal.x) * area;
  (*q1y) = 1.0/3.0 * v1x * (-normal.y) * area;
  (*q1z) = 1.0/3.0 * v1x * (-normal.z) * area;
  (*q2x) = 1.0/3.0 * v1x * (-normal.x) * area;
  (*q2y) = 1.0/3.0 * v1x * (-normal.y) * area;
  (*q2z) = 1.0/3.0 * v1x * (-normal.z) * area;
  (*q3x) = 1.0/3.0 * v1x * (-normal.x) * area;
  (*q3y) = 1.0/3.0 * v1x * (-normal.y) * area;
  (*q3z) = 1.0/3.0 * v1x * (-normal.z) * area;
 }                                   
                                     
} /* End of TETR4Load */             
                                     
                                     
/*                                   
%F This method allocates memory for a TETR4 element and fills its data 
   with the specific data.
%i Element label (number)
%o Element descriptor.
*/
static void TETR4New( int label, int matid, int intord, int tckid, 
                      sElement **elm, sElement **elist, sNode *nodes )
{
 sTetr4Data *data = 0L;

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

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

/* Fill TETR4 element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 4;
 data->NumTensComp = 6;
 data->effdef      = 0.0;
 data->prof        = 0.0;
 data->istr[0]     = 0.0;
 data->istr[1]     = 0.0;
 data->istr[2]     = 0.0;
 data->istr[3]     = 0.0;
 data->istr[4]     = 0.0;
 data->istr[5]     = 0.0;
 
/* Fill element descriptor
 */
 (*elm)->type    = TETR4;
 (*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 TETR4New */


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

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

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

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

} /* End of TETR4Free */


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

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

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

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

} /* End of TETR4Read */


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

/* Get element data
 */
 data = (sTetr4Data *)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] = szz;
 data->istr[3] = sxy;
 data->istr[4] = sxz;
 data->istr[5] = syz;

 return 1;

} /* End of TETR4ReadInitStress */

/*
%F This method reads the TETRA4 profile
*/
static void TETRA4ReadProfile( sElement *elm )
{
 double prof;
 sTetr4Data *data = 0L;

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

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

} /* End of TETRA4ReadProfile */

/*
%F This method computes the mass matrix for the TETR4 element.
%i Element descriptor.
%o Element mass matrix.
*/
static void TETR4MassMatrix( sElement *elm, double *mass )
{
 sTetr4Data *data = 0L;
 sCoord      coord[4];
 double      shape[4], jac[3][3];
 double      dens, detjac, coeff;
 int         i, j;

/* Initialize mass vector
 */
 for( i = 0; i < 12; i++ )
  mass[i] = 0.0;

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

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

/* Compute shape functions
 */
 _TETR4ShapeFunc( _GPoint, _GPoint, _GPoint, shape );

/* Compute the jacobian matrix
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &detjac, jac );

/* Rigidity coefficient
 */ 
 coeff = dens * detjac * _GWeight;

/* Compute mass
 */
 for( i = 0; i < 4; i++ ) {
  for( j = 0; j < 4; j++ ) {
   mass[3*i]   += coeff * shape[i] * shape[j];
   mass[3*i+1] += coeff * shape[i] * shape[j];
   mass[3*i+2] += coeff * shape[i] * shape[j];
  }
 }

} /* End of TETR4MassMatrix */


/*
%F This functions computes the element rigidity coefficient.
%i Element descriptor.
%o Element rigidity coefficient.
*/
static double TETR4RigidCoeff( sElement *elm )
{
 sCoord coord[4];
 double jac[3][3];
 double detjac, volume = 0.0;

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

/* Compute the jacobian matrix
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &detjac, jac );

/* Rigidity coefficient
 */
 volume = detjac * _GWeight;

 return volume;

} /* End of TETR4RigidCoeff */


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

 data = (sTetr4Data *)elm->data;

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

} /* End of TETR4Connect */


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

 (*nnodes) = data->NumNodes;

} /* End of TETR4NumNodes */


/*
%F
*/
static void TETR4AssVector( sElement *elm, double *GMatrix, double *matrix )
{
 int         i;
 int         conn[4];
 sTetr4Data *data;

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

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

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

} /* End of TETR4AssVector */


/*
%F
*/

static void TETR4StressStrain( sElement *elm, double dt, double *U,
                                              double *yield, sTensor *stre, 
                                              sTensor *stra )
{
 sTetr4Data *data = 0L;
 sMaterial  *mat  = 0L;
 sCoord      coord[4];
 double      cm[6][6];
 double      bm[6][12];
 double      ue[12];
 double      def[6], str[6];
 double      pf = 1.0;
 int         i, j;
 double     poropressure;
 double     lambda, gamma;
 double     distance;
 
 
/* Get element descriptor
 */
 data = (sTetr4Data *)elm->data;
 
/* Get Material object
 */
 mat = MatList[data->matid-1];
 
/* Get element coordinates
 */
 _TETR4GetCoord( elm, coord );

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

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

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

/* Get element stress x strain matrix
 */
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, bm );

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

/* Compute the element stress
 */
 for( i = 0; i < 6; i++ ) {
  str[i] = 0.0;
  for( j = 0; j < 6; j++ )
   str[i] += cm[i][j] * def[j];
  str[i] += data->istr[i];
 }
 
/* Get element lambda
 */
 lambda = MatList[data->matid-1]->Lambda;
 
 if(lambda>0.0){
 
  /* Get element gamma
  */
  MatDensity( mat, &gamma );
 
  /* Get element prof
  */
  distance = data->prof;
 
  /* Compute the poropressure
  */ 
  poropressure = lambda*gamma*GForce*distance*Config.loadfactor;
  
  /* Compute efective stress
  */
  for( i = 0; i < 3; 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 < 3; 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->zz = str[2];
 stre->xy = str[3];
 stre->yz = str[4];
 stre->xz = str[5];

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

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

} /* End of TETR4StressStrain */


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

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

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

/* Compute element B matrix
 */
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, bm );

/* Compute rigidity coefficient
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &coeff, jac );

 coeff *= _GWeight;

/* Check for the initial stress
 */
 if( stress == 0L ) {
  for( i = 0; i < (3*data->NumNodes); i++ ) {
   intforce[i] = 0.0;
   for( j = 0; j < 6; 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->zz;
 str[3] = stress->xy;
 str[4] = stress->yz;
 str[5] = stress->xz;

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

} /* End of TETR4InterForce */


/*
%F
*/
static void TETR4WriteStress( sElement *elm, FILE *out, double *U, double *V )
{
 sTensor  stress, strain, strd;
 sPTensor pstr;
 double   yield;

/* Check to see if the element was rezoned
 */
 if( elm->rezone == DONE ) {

 /* Set the stress and strain to 0
  */
  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
  */
  TETR4StressStrain( elm, 0.0, U, &yield, &stress, &strain );
 }

/* Compute desviatory stress
 */
 PrincipalTensor( &stress, &pstr);

 strd.xx = stress.xx - ((pstr.dir1 + pstr.dir2 + pstr.dir3)/3.0);
 strd.yy = stress.yy - ((pstr.dir1 + pstr.dir2 + pstr.dir3)/3.0);
 strd.zz = stress.zz - ((pstr.dir1 + pstr.dir2 + pstr.dir3)/3.0);

/* Write element stress and strain
 */
 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 TETR4WriteStress */


/*
%F
*/
static void TETR4WriteNodalResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  str, strd, def;
 sPTensor pstress, pstrain;
 double   yield;
 int      i;

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

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 PrincipalTensor( &str, &pstress );
 PrincipalTensor( &def, &pstrain );
 for( i = 0; i < 4; i++ ) {
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", str.xx, str.yy, str.zz );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", str.xy, str.xz, str.yz );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", pstress.dir1,
                                               pstress.dir2,
                                               pstress.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", strd.xx, strd.yy, strd.zz );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", def.xx, def.yy, def.zz );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", def.xy, def.xz, def.yz );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", pstrain.dir1,
                                               pstrain.dir2,
                                               pstrain.dir3 );
  fprintf( out, "%+10.5e\n", yield );
 }

} /* End of TETR4WriteNodalResult */


/*
%F
*/
static void TETR4WriteGaussResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  stress, strain, strd;
 sPTensor pstress, pstrain;
 double   yield;
 int      j;

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 1 );
 PrincipalTensor( &stress, &pstress );
 PrincipalTensor( &strain, &pstrain );
 
 for( j = 0; j < Config.numgaussresults; j++){
   if((strcmp( Config.gaussresults[j], "Sxx" ) == 0))
    fprintf( out, "%+10.5e\t", stress.xx );
   else if((strcmp( Config.gaussresults[j], "Syy" ) == 0))
    fprintf( out, "%+10.5e\t", stress.yy );
   else if((strcmp( Config.gaussresults[j], "Szz" ) == 0))
    fprintf( out, "%+10.5e\t", stress.zz );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", stress.xy );
   else if((strcmp( Config.gaussresults[j], "Sxz" ) == 0))
    fprintf( out, "%+10.5e\t", stress.xz );
   else if((strcmp( Config.gaussresults[j], "Syz" ) == 0))
    fprintf( out, "%+10.5e\t", stress.yz );
   else if((strcmp( Config.gaussresults[j], "S1" ) == 0))
    fprintf( out, "%+10.5e\t", pstress.dir1 );
   else if((strcmp( Config.gaussresults[j], "S2" ) == 0))
    fprintf( out, "%+10.5e\t", pstress.dir2 );
   else if((strcmp( Config.gaussresults[j], "S3" ) == 0))
    fprintf( out, "%+10.5e\t", pstress.dir3 );
   else if((strcmp( Config.gaussresults[j], "Sdx" ) == 0))
    fprintf( out, "%+10.5e\t", strd.xx );
   else if((strcmp( Config.gaussresults[j], "Sdy" ) == 0))
    fprintf( out, "%+10.5e\t", strd.yy );
   else if((strcmp( Config.gaussresults[j], "Sdz" ) == 0))
    fprintf( out, "%+10.5e\t", strd.zz );
   else if((strcmp( Config.gaussresults[j], "Exx" ) == 0))
    fprintf( out, "%+10.5e\t", strain.xx );
   else if((strcmp( Config.gaussresults[j], "Eyy" ) == 0))
    fprintf( out, "%+10.5e\t", strain.yy );
   else if((strcmp( Config.gaussresults[j], "Ezz" ) == 0))
    fprintf( out, "%+10.5e\t", strain.zz );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", strain.xy );
   else if((strcmp( Config.gaussresults[j], "Exz" ) == 0))
    fprintf( out, "%+10.5e\t", strain.xz );
   else if((strcmp( Config.gaussresults[j], "Eyz" ) == 0))
    fprintf( out, "%+10.5e\t", strain.yz );
   else if((strcmp( Config.gaussresults[j], "E1" ) == 0))
    fprintf( out, "%+10.5e\t", pstrain.dir1 );
   else if((strcmp( Config.gaussresults[j], "E2" ) == 0))
    fprintf( out, "%+10.5e\t", pstrain.dir2 );
   else if((strcmp( Config.gaussresults[j], "E3" ) == 0))
    fprintf( out, "%+10.5e\t", pstrain.dir3 );
   else if((strcmp( Config.gaussresults[j], "Ev" ) == 0))
    fprintf( out, "%+10.5e\t", (pstrain.dir1+pstrain.dir2+pstrain.dir3) );
  }

} /* End of TETR4WriteGaussResult */


/*
%F
*/
static void TETR4WriteGaussVectorResult( sElement *elm, int version,
                                         FILE *out, FILE *tmp )
{
 sTensor  stress, strain, strd;
 sPTensor pstress, pstrain;
 double   yield;

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 1 );
 PrincipalTensor( &stress, &pstress );
 PrincipalTensor( &strain, &pstrain );
 if( version > 1 ) {
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.dir1, pstrain.dir1 );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.dir2, pstrain.dir2 );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.dir3, pstrain.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos1x, pstrain.cos1x );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos2x, pstrain.cos2x );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos3x, pstrain.cos3x );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos1y, pstrain.cos1y );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos2y, pstrain.cos2y );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos3y, pstrain.cos3y );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos1z, pstrain.cos1z );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos2z, pstrain.cos2z );
  fprintf( out, "%+10.5e\t%+10.5e\n", pstress.cos3z, pstrain.cos3z );
 }
 else {
  if( pstress.dir1 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (pstress.dir1*pstress.cos1x),
                                                (pstress.dir1*pstress.cos1y),
                                                 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

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

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

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

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

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

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

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

} /* End of TETR4WriteGaussVectorResult */

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

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

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

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

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

/* Get element stress x strain matrix
 */
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, bm );

/* Compute rotation factor
 */
 rot = 0.0;
 for( i = 0; i < 12; 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;
 szz = stre->zz;

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

} /* End of TETR4UpdateStress */


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

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

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

/* Get element data
 */
 data = (sTetr4Data *)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 TETR4TimeStep */


/*
%F This function computes the load due to a gravitational field. The 
   gravitational load is given by the product between the element 
   volume and the gravity field components.
%i Element descriptor.
%o Load in X direction.
%o Load in Y direction.
%o Load in Z direction.
*/
static void TETR4Gravity( sElement *elm, double *qx, double *qy, double *qz )
{
 sTetr4Data *data = 0;
 sCoord      coord[4];
 double      volume, mass, gamma;
 int         matid;

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

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

/* Compute element area
 */
 volume = _TETR4Volume( coord );

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

/* Compute the gravitational load
 */
 mass = (gamma * volume) / 4.0;

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

} /* End of TETR4Gravity */


/*
%F This function computes the percolation forces for the TETR4 element.
*/
static void TETR4PercForces( sElement *elm, double *pforce )
{
 sCoord      coord[4];
 sTetr4Data *data = 0L;
 sMaterial  *mat  = 0L;
 double      bm[6][12], phi[4], jac[3][3], shape[4];
 double      detjac;
 sTensor     grad;
 int         i, k;

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

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

/* 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;

 for( i = 0; i < 12; i++ )
  pforce[i] = 0.0;

/* Get element stress x strain matrix
 */
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, bm );

/* Get element jacobiam
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &detjac, jac );

 detjac *= _GWeight * mat->Rhof;

/* Get element shape functions
 */
 _TETR4ShapeFunc( _GPoint, _GPoint, _GPoint, shape );

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

 for( k = 0; k < 4; k++ ) {
  pforce[3*k]   += detjac * shape[k] * grad.xx;
  pforce[3*k+1] += detjac * shape[k] * grad.yy;
  pforce[3*k+2] += detjac * shape[k] * grad.zz;
 }

} /* End of TETR4PercForce */


/*
%F
*/
static void TETR4SetPressure( sElement *elm, double pot )
{
 sCoord      coord[4];
 sTetr4Data *data = 0L;
 int         i, id;

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

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

/* Compute and set the pressure for the element node
 */
 for( i = 0; i < 4; i++ ) {
  id = data->inc[i] - 1;
  elm->nodes[id].dof.psi = pot - coord[i].y;
 }

} /* End of TETR4SetPressure */

/*
%F
*/
static void TETR4SetInitStress( sElement *elm, sTensor *istress )
{
 sTetr4Data *data = 0L;
 int         i;

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

/* Set initial stress
 */
 for( i = 0; i < 4; i++ ) {
  data->istr[0] = istress->xx;
  data->istr[1] = istress->yy;
  data->istr[2] = istress->zz;
  data->istr[3] = istress->xy;
  data->istr[4] = istress->yz;
  data->istr[5] = istress->xz;
 }

} /* End of TETR4SetInitStress */

/*
%F
*/
static void TETR4ViscoForce( sElement *elm, double timeStep, sTensor *stress,
                             double *vforce )
{
 int         i, l, m;
 double      coeff;
 double      bm[6][12], jac[3][3], cm[6][6];
 double      stav[6], strv[6];
 sTensor     vsta;
 sCoord      coord[4];
 sTetr4Data *data = 0L;
 sMaterial  *mat  = 0L;

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

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

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

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

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

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

/* Get element stress x strain matrix
 */
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, bm );

/* Get element jacobiam
 */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &coeff, jac );

 coeff *= _GWeight;

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

/* Get viscos strain components
 */
 stav[0] = vsta.xx;
 stav[1] = vsta.yy;
 stav[2] = vsta.zz;
 stav[3] = vsta.xy;
 stav[4] = vsta.yz;
 stav[5] = vsta.xz;

/* Compute viscos stress
 */
 for( m = 0; m < 6; m++ ) {
  strv[m] = 0.0;
  for( l = 0; l < 6; l++ )
   strv[m] += cm[m][l] * stav[l];
 }

/* Compute element visco forces for the current gauss point
 */
 for( m = 0; m < (3*data->NumNodes); m++ ) {
  vforce[m] = 0.0;
  for( l = 0; l < 6; l++ )
   vforce[m] += bm[l][m] * strv[l];
  vforce[m] *= coeff;
 } 

} /* End of TETR4ViscoForce */

/*============================ TETR4Kmatrix ===========================*/

static void TETR4KMatrix( sElement *elm, double k[24][24] )
{           
 int    i, j, l;           
 sCoord  coord[4];          
 double jac[3][3];         
 double detjac;            
 double C[6][6];          
 double B[6][12];          
 double BTC[12][6];        
 double ke[24][24]; 
 sTetr4Data  *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 = (sTetr4Data *)elm->data;
 
/* Pega material
 */
 mat = MatList[data->matid-1];
 
 /* Pega as coordenadas nodais
 */
 _TETR4GetCoord(elm, coord);
 
 /* Pega matriz constitutiva [C]*/
 MatConstitutiveMatrix( mat, C );

 /* Calcula matriz B*/
 _TETR4BMatrix( _GPoint, _GPoint, _GPoint, coord, B );
 /* Calcula o determinante do jacobiano */
 _TETR4Jacobian( _GPoint, _GPoint, _GPoint, coord, &detjac, jac );

 /* Calcula [B] transposto * [C] */
 for( i = 0; i < 12; i++ ) {
  for( j = 0; j < 6; j++)  {
   BTC[i][j] =0.0;
   for(l=0; l<6; l++){
    BTC[i][j] += B[l][i]*C[l][j];
   }
  }
 }
 /* Calcula [B] transposto * [C] * [B] * |jac|*/
 for( i = 0; i < 12; i++ ) {
  for( j = 0; j < 12; j++)  {
   ke[i][j] =0.0;
   for(l=0; l<6; l++){
    ke[i][j] += BTC[i][l]*B[l][j]*detjac*_GWeight;
   }
  }
 }
 /* Calcula matriz de rigidez do elemento*/
 for( i = 0; i < 24; i++ ) {
  for( j = 0; j < 24; j++)  {
   k[i][j] += ke[i][j];
  }
 }
   
} /* TETR4KMatrix */

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

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

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

} /* End of _TETR4GetDof */

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

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

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

} /* End of TETR4GetInc */


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

/*
%F This function initializes the sub-class "TETR4".
*/
void TETR4Init( void );
void TETR4Init( void )
{
/* Initialize TETR4 element sub-class
 */
 ElmClass[TETR4].new               = TETR4New;
 ElmClass[TETR4].free              = TETR4Free;
 ElmClass[TETR4].read              = TETR4Read;
 ElmClass[TETR4].readinitstr       = TETR4ReadInitStress;
 ElmClass[TETR4].readprofile       = TETRA4ReadProfile;
 ElmClass[TETR4].mass              = TETR4MassMatrix;
 ElmClass[TETR4].rigidcoeff        = TETR4RigidCoeff;
 ElmClass[TETR4].load              = TETR4Load; 
 ElmClass[TETR4].connect           = TETR4Connect;
 ElmClass[TETR4].numnodes          = TETR4NumNodes;
 ElmClass[TETR4].gravity           = TETR4Gravity;
 ElmClass[TETR4].assvector         = TETR4AssVector;
 ElmClass[TETR4].strstrain         = TETR4StressStrain;
 ElmClass[TETR4].timestep          = TETR4TimeStep;
 ElmClass[TETR4].intforce          = TETR4InterForce;
 ElmClass[TETR4].writestr          = TETR4WriteStress;
 ElmClass[TETR4].writendlresult    = TETR4WriteNodalResult;
 ElmClass[TETR4].writegauresult    = TETR4WriteGaussResult;
 ElmClass[TETR4].writegauvecresult = TETR4WriteGaussVectorResult;
 ElmClass[TETR4].updatestress      = TETR4UpdateStress;
 ElmClass[TETR4].percforce         = TETR4PercForces;
 ElmClass[TETR4].setpressure       = TETR4SetPressure;
 ElmClass[TETR4].setinitstress     = TETR4SetInitStress;
 ElmClass[TETR4].viscoforce        = TETR4ViscoForce;
 ElmClass[TETR4].jacobian          = 0L;
 ElmClass[TETR4].volume            = 0L;
 ElmClass[TETR4].KMatrix           = TETR4KMatrix;
 ElmClass[TETR4].GetDof            = TETR4GetDof;
 ElmClass[TETR4].GetInc            = TETR4GetInc;
} /* End of TETR4Init */


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

---