q4.c

Go to the documentation of this file.

/*
%M This modules contains the Q4 element sub-class methods and definitions
%a Joao Luiz Elias Campos.
%d September 7th, 1998.
%r $Id: q4.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:  23-Fev-2006    Juan Pablo Ibañez
*    Modificada a funcão Q4ReadInitStress que passa a retornar 
*    um valor int. 
*
*
*  Modified:  Agosto-06  André Luis Muller
*    Criada a função _Q4ElementCenter que retorna as coordenadas do centro de
*    um elemento do tipo Q4.
*    Adicionado a função Q4StressStrain 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 Q4InterForce para possibilitar
*    a condição de tensões iniciais.
*
*    Modificação na função Q4New. Inicialização das tensões iniciais
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função Q4WriteGaussResults 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 "xgplib.h"
#include "alg.h"

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

/*
%T Q4 element data definition
*/
typedef struct _q4elm {
 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[4][3];           /* Element initial stress           */
 double effdef[4];            /* Effective deformation            */
 double prof;                 /* Element profile                  */
} sQ4Data;


/*
%V Gauss point coordinates and its weights
*/
static double _GPoint[2]  = { -0.577350269189626, 0.577350269189626 };
static double _GWeight[2] = {  1.000000000000000, 1.000000000000000 };


/*
%V Transformation matrix
*/
static double _TRMatrix[4][4] = {
                                 {  1.116025404,  0.250000000,
                                    0.250000000, -0.616025404 },
                                 {  0.250000000,  1.116025404,
                                   -0.616025404,  0.250000000 },
                                 { -0.616025404,  0.250000000,
                                    0.250000000,  1.116025404 },
                                 {  0.250000000, -0.616025404,
                                    1.116025404,  0.250000000 }
                                };

/*
%V Element edges ids.
*/
static int _Edges[4][2] = { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 } };

static double GForce = 9.85;

/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _Q4GetCoord       ( sElement *, sCoord [] );
static void   _Q4ShapeFunc      ( double, double, double [] );
static void   _Q4Jacobian       ( double, double, sCoord [], double *,
                                  double [2][2] );
static void   _Q4BMatrix        ( double, double, sCoord [], double [3][8] );
static void   _Q4DerivShp       ( double, double, sDerivRST [] );
static void   _Q4GetDisplacement( sElement *, double *, double [] );
static double _Q4MinHeight      ( sCoord * );
static double _Q4CalcLength     ( sCoord *, sCoord * );
static void   _Q4ElementCenter  ( sElement *, sCoord [] );
static double _Q4PoroPressure   ( sElement *elm, double lambda );

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

 /* Element element data
 */
 data = (sQ4Data *)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 _Q4ElementCenter */

/*
%F This function computes the length of a given segment.
*/
static double _Q4CalcLength( sCoord *pt1, sCoord *pt2 )
{
 double dx = pt2->x - pt1->x;
 double dy = pt2->y - pt1->y;
 double dz = pt2->z - pt1->z;

 return sqrt( (dx * dx) + (dy * dy) + (dz * dz) );

} /* End of _BRICK8CalcLength */

/*
%F This function computes the minimun height for the Q4 element.
*/
static double _Q4MinHeight( sCoord *coord )
{
 double minlen;
 int    i;

/* For the first element edge.
 */
 minlen = _Q4CalcLength( &coord[_Edges[0][0]], &coord[_Edges[0][1]] );

/* For the remain element edge
 */
 for( i = 1; i < 4; i++ ) {
 /* Compute the edge length
  */
  double len = _Q4CalcLength( &coord[_Edges[i][0]], &coord[_Edges[i][1]] );

 /* Get the smallest.
  */
  if( minlen > len )
   minlen = len;
 }

 return minlen;

} /* End of _Q4MinHeight */


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

/* Element element data
 */
 data = (sQ4Data *)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 = 0.0;
 }

} /* End of _Q4GetCoord */


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

/* Get element descriptor
 */
 data = (sQ4Data *)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 _Q4GetDisplacement */


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

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

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

/* 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 shape functions derivatives
 */
 _Q4DerivShp( r, s, dshp );

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

} /* End of _Q4BMatrix */


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

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

/* Compute the jacobian matrix
 */
 for( i = 0; i < 2; i++ )
  jac[i][0] = jac[i][1] = 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[1][0] += dmap[i].s * coord[i].x;
  jac[1][1] += dmap[i].s * coord[i].y;
 }

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

} /* End of _Q4Jacobian */


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

} /* End of _Q4ShapeFunc */


/*
%F This function computes the shape functions derivatives for the
   quadrilateral 4 node element.
*/
static void _Q4DerivShp( double r, double s, sDerivRST dshp[4] )
{
/* Compute shape functions derivatives
 */
 dshp[0].r = -0.25 * (1.0 - s);
 dshp[1].r =  0.25 * (1.0 - s);
 dshp[2].r =  0.25 * (1.0 + s);
 dshp[3].r = -0.25 * (1.0 + s);

 dshp[0].s = -0.25 * (1.0 - r);
 dshp[1].s = -0.25 * (1.0 + r);
 dshp[2].s =  0.25 * (1.0 + r);
 dshp[3].s =  0.25 * (1.0 - r);

} /* End of _INFINITEDerivShp */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   Q4New( int, int, int, int, sElement **, sElement **, sNode * );
static void   Q4Free                  ( sElement * );
static void   Q4Read                  ( sElement * );
static int    Q4ReadInitStress        ( sElement * );
static void   Q4ReadProfile           ( sElement * );
static void   Q4MassMatrix            ( sElement *, double * );

static double Q4RigidCoeff            ( sElement * );
static void   Q4Connect               ( sElement *, int * );
static void   Q4NumNodes              ( sElement *, int * );
static void   Q4Gravity               ( sElement *, double *, double *,
                                                    double * );
static void   Q4AssVector             ( sElement *, double *, double * );
static void   Q4TimeStep              ( sElement *, double * );
static void   Q4InterForce            ( sElement *, sTensor *, double * );
static void   Q4StressStrain          ( sElement *, double, double *, double *,
                                                    sTensor *, sTensor * );
static void   Q4WriteStress           ( sElement *, FILE *, double *,
                                                            double * );
static void   Q4WriteGaussResult      ( sElement *, FILE *, FILE * );
static void   Q4WriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   Q4WriteNodalResult      ( sElement *, FILE *, FILE * );
static void   Q4UpdateStress          ( sElement *, double, double *,
                                                            sTensor * );
static void   Q4PercForces            ( sElement *, double * );
static void   Q4SetPressure           ( sElement *, double );
static void   Q4SetInitStress         ( sElement *, sTensor * );
static void   Q4ViscoForce            ( sElement *, double, sTensor *,
                                        double * );
static void Q4KMatrix( sElement *, double [24][24] );                                                                         
static void Q4GetDof( sElement *, int [24], int * );

/*
%F This method allocates memory for a Q4 element and fills its data with 
   the specific data.
%i Element label (number)
%o Element descriptor.
*/
static void Q4New( int label, int matid, int intord, int tckid, 
                   sElement **elm, sElement **elist, sNode *nodes )
{
 sQ4Data *data = 0L;
 int      i;

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

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

/* Fill Q4 element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 4;
 data->NumTensComp = 4;
 data->prof        = 0.0;
 for( i = 0; i < 4; i++ )
  data->effdef[i] = 0.0;

 for( i =0; i < 4; i++){
  data->istr[i][0] = 0.0;
  data->istr[i][1] = 0.0;
  data->istr[i][2] = 0.0;
 }

/* Fill element descriptor
 */
 (*elm)->type    = Q4;
 (*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 Q4New */


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

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

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

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

} /* End of Q4Free */


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

/* Get the element data
 */
 data = (sQ4Data *)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 Q4Read */


/*
%F This method reads the Q4 initial stress
*/
static int Q4ReadInitStress( sElement *elm )
{
 int      i, numpg = 4;
 double   sxx, syy, sxy ,szz, sxz, syz;
 sQ4Data *data = 0L;

/* Get element data
 */
 data = (sQ4Data *)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 != 4 ) {
  printf( "\n\nNumber of gauss points must be equal to four.\n\n" );
  return 0;
 }

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

 return 1;

} /* End of Q4ReadInitStress */

/*
%F This method reads the Q4 profile
*/
static void Q4ReadProfile( sElement *elm )
{
 double prof;
 sQ4Data *data = 0L;

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

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

} /* End of Q4ReadProfile */

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

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

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

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

/* Calculate the element node contribution for the element mass
 */
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Compute shape functions
   */
   _Q4ShapeFunc( _GPoint[i], _GPoint[j], shape );

  /* Compute the jacobian matrix
   */
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &detjac, jac );

  /* Rigidity coefficient
   */
   coeff = dens * detjac * _GWeight[i] * _GWeight[j];

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

  }
 }
 
} /* End of Q4MassMatrix */


/*
%F This functions computes the element rigidity coefficient.
%i Element descriptor.
%o Element rigidity coefficient.
*/
static double Q4RigidCoeff( sElement *elm )
{
 sCoord coord[4];
 double shape[4], jac[2][2];
 double detjac, area = 0.0;
 int    i, j;

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

/* Calculate the element node contribution for the element area
 */
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Compute shape functions
   */
   _Q4ShapeFunc( _GPoint[i], _GPoint[j], shape );

  /* Compute the jacobian matrix
   */
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &detjac, jac );

  /* Rigidity coefficient
   */
   area += detjac * _GWeight[i] * _GWeight[j];
  }
 }

 return( area );

} /* End of Q4RigidCoeff */


/*
%F This function computes the equivalent load applied to the element 
   face.
*/
static void Q4Load( 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 Q4Load */


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

 data = (sQ4Data *)elm->data;

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

} /* End of Q4Connect */


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

 (*nnodes) = data->NumNodes;

} /* End of Q4NumNodes */


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

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

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

/* Compute element area
 */
 area = Q4RigidCoeff( elm );

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

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

} /* End of Q4Gravity */


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

/* Get the element data
 */
 data = (sQ4Data *)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 Q4AssVector */


/*
%F
*/


static void Q4StressStrain( sElement *elm, double dt, double *U, double *yield,
                                           sTensor *stre, sTensor *stra )
{
 sQ4Data    *data = 0L;
 sMaterial  *mat  = 0L;
 sCoord      coord[4];
 double      cm[6][6];
 double      bm[3][8];
 double      ue[8];
 double      def[3], str[5];
 double      pf = 1.0;
 int         i, j, k, l, npg;
 double     poropressure;
 double     lambda;
 double     nu;

/* Initialize stress and strain vector
 */
 for( i = 0; i < 4; i++ ) {
  stre[i].xx = stre[i].xy = stre[i].xz = 
               stre[i].yy = stre[i].yz = 
                            stre[i].zz = 0.0;
  if( stra != 0L ) {
   stra[i].xx = stra[i].xy = stra[i].xz = 
                stra[i].yy = stra[i].yz = 
                             stra[i].zz = 0.0;
  }
 }

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

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

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

/* Get element displacement
 */
 _Q4GetDisplacement( elm, U, ue );
 
/* Get poisson coeff.
 */
 MatNuParameter( mat, &nu );


/* Get element lambda
 */
 lambda = MatList[data->matid-1]->Lambda;
 
 if(lambda>0.0){
  poropressure = _Q4PoroPressure(elm,lambda );
 }
/* Compute stress and strain for each gauss point
 */
 npg = 0;
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Update the material parameter
   */
   MatUpdateParameter( mat, data->effdef[npg] );

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

  /* Get element stress x strain matrix
   */
   _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

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

  /* Compute the element stress
   */
   for( k = 0; k < 3; k++ ) {
    str[k] = 0.0;
    for( l = 0; l < 3; l++ )
     str[k] += cm[k][l] * def[l];
    str[k] += data->istr[npg][k];
   }
   str[3] = str[4] = 0.0;

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

  /* Set stress values
   */
   stre[npg].xx = str[0];
   stre[npg].yy = str[1];
   stre[npg].xy = str[2];
   stre[npg].rf = str[4];
   
 /* Compute off plane stress
  */
   ElementOffPlaneStress( nu, &stre[npg] );

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

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

  /* Update gauss point number
   */
   npg++;

  }
 }

} /* End of Q4StressStrain */

static double _Q4PoroPressure( sElement *elm, double lambda )
{
 sQ4Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     poropressure;
 double     gamma;
 double     distance;
 
 /* Get element descriptor
 */
 data = (sQ4Data *)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 _Q4PoroPressure */

/*
%F This function computes the internal forces for the Q4 element
*/
static void Q4InterForce( sElement *elm, sTensor *stress, double *intforce )
{
 int      i, j, k, l, npg;
 double   coeff;
 double   bm[3][8], str[3], jac[2][2], iforce[8];
 sCoord   coord[4];
 sQ4Data *data = 0L;

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

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

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

/* Compute internal force for each gauss point
 */
/* Check for the initial stress
 */
 if( stress == 0L ) {
  npg = 0;
  for( j = 0; j < 2; j++ ) {
   for( i = 0; i < 2; i++ ) {
   /* Compute element B matrix
   */
    _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

   /* Compute rigidity coefficient
   */
    _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &coeff, jac );

    coeff *= _GWeight[i] * _GWeight[j];

    for( k = 0; k < (2*data->NumNodes); k++ ) {
     iforce[k] = 0.0;
     for( l = 0; l < 3; l++ )
      iforce[k] += bm[l][k] * data->istr[npg][l];
     iforce[k] *= coeff;
    }
    /* Add in the element vector
    */
    for( k = 0; k < (2*data->NumNodes); k++ )
     intforce[k] += iforce[k];
    npg++;
   }
  }
  return;
 }
 npg = 0;
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Compute element B matrix
  */
   _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

  /* Compute rigidity coefficient
  */
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &coeff, jac );

   coeff *= _GWeight[i] * _GWeight[j];

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

  /* Compute element internal forces for the current gauss point
   */
   for( k = 0; k < (2*data->NumNodes); k++ ) {
    iforce[k] = 0.0;
    for( l = 0; l < 3; l++ )
     iforce[k] += bm[l][k] * str[l];
    iforce[k] *= coeff;
   } 

  /* Add in the element vector
   */
   for( k = 0; k < (2*data->NumNodes); k++ )
    intforce[k] += iforce[k];

  /* Update gauss point number
   */
   npg++;

  }
 }

 for( i = 0; i < (2*data->NumNodes); i++ )
  intforce[i] *= -1.0;
 
} /* End of Q4InterForce */


/*
%F
*/
static void Q4WriteStress( sElement *elm, FILE *out, double *U, double *V )
{
 int     i;
 sTensor  stress[4];
 sTensor  strain[4];
 sTensor  strd[4];
 double   yield[4];
 double   iso;

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

 /* Set the stress and strain vectors to 0
  */
  for( i = 0; i < 4; i++ ) {
   stress[i].xx = stress[i].xy = stress[i].xz = 
                  stress[i].yy = stress[i].yz = 
                                 stress[i].zz = 0.0;
   strd[i].xx = strd[i].xy = strd[i].xz = 
                strd[i].yy = strd[i].yz = 
                             strd[i].zz = 0.0;
   strain[i].xx = strain[i].xy = strain[i].xz = 
                  strain[i].yy = strain[i].yz = 
                                 strain[i].zz = 0.0;
   yield[i] = 0.0;
  }
 }
 else {
 /* Compute the element stress and strain
  */
  Q4StressStrain( elm, 0.0, U, yield, stress, strain );
 }

/* Compute desviatory stress
 */
 for( i = 0; i < 4; i++ ) {
  iso = (stress[i].xx + stress[i].yy + stress[i].zz)/3.0;
  strd[i].xx = stress[i].xx - iso;
  strd[i].yy = stress[i].yy - iso;
 }

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

} /* End of Q4WriteStress */


/*
%F
*/
static void Q4WriteNodalResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  strgau[4],strndl[4];
 sTensor  strdgau[4],strdndl[4];
 sTensor  defgau[4],defndl[4];
 sPTensor pstress, pstrain;
 double   yieldg[4], yieldn[4];
 int      i, j;
 double   iso[4];

/* Read results from temporary file
 */
 fread( strgau, sizeof(sTensor), 4, tmp );
 fread( defgau, sizeof(sTensor), 4, tmp );
 fread( strdgau, sizeof(sTensor), 4, tmp );
 fread( yieldg, sizeof(double), 4, tmp );

/* Compute the nodal results
 */
 for( i = 0; i < 4; i++ ) {
  strndl[i].xx = strndl[i].xy = strndl[i].xz =
                 strndl[i].yy = strndl[i].yz =
                                strndl[i].zz = 0.0;
  strdndl[i].xx = strdndl[i].xy = strdndl[i].xz =
                  strdndl[i].yy = strdndl[i].yz =
                                  strdndl[i].zz = 0.0;
  defndl[i].xx = defndl[i].xy = defndl[i].xz = 
                 defndl[i].yy = defndl[i].yz =
                                defndl[i].zz = 0.0;
  yieldn[i] = 0.0;

  for( j = 0; j < 4; j++ ) {
   strndl[i].xx  += _TRMatrix[i][j] * strgau[j].xx;
   strndl[i].yy  += _TRMatrix[i][j] * strgau[j].yy;
   strndl[i].zz  += _TRMatrix[i][j] * strgau[j].zz;
   strndl[i].xy  += _TRMatrix[i][j] * strgau[j].xy;
   strdndl[i].xx += _TRMatrix[i][j] * strdgau[j].xx;
   strdndl[i].yy += _TRMatrix[i][j] * strdgau[j].yy;
   defndl[i].xx  += _TRMatrix[i][j] * defgau[j].xx;
   defndl[i].yy  += _TRMatrix[i][j] * defgau[j].yy;
   defndl[i].xy  += _TRMatrix[i][j] * defgau[j].xy;
   yieldn[i]     += _TRMatrix[i][j] * yieldg[j];
   iso[i] = (strndl[i].xx + strndl[i].yy + strndl[i].zz)/3.0;
  }
 }

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 for( i = 0; i < 4; i++ ) {
  PrincipalTensor( &strndl[i], &pstress );
  PrincipalTensor( &defndl[i], &pstrain );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t%+10.5e\t", strndl[i].xx,
                                                        strndl[i].yy,
                                                        strndl[i].zz,
                                                        strndl[i].xy );
  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", strdndl[i].xx,
                                               strdndl[i].yy,
                                               iso[i]);
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", defndl[i].xx,
                                               defndl[i].yy,
                                               defndl[i].xy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", pstrain.dir1,
                                               pstrain.dir2,
                                               pstrain.dir3 );
  fprintf( out, "%+10.5e\n", yieldn[i] );
 }

} /* End of Q4WriteNodalResult */


/*
%F
*/
static void Q4WriteGaussResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  stress[4];
 sTensor  strain[4];
 sTensor  strd[4],estr[4];
 sPTensor pstress, pstrain, epsts, dpsts;
 double   yield[4];
 int      i,j;
 double   iso, piso;
 double   lambda;
 double   poropressure = 0.0;
 sMaterial *mat  = 0L;
 sQ4Data   *data = 0L;

/* Read results from temporary file
 */
 fread( stress, sizeof(sTensor), 4, tmp );
 fread( strain, sizeof(sTensor), 4, tmp );
 fread( strd, sizeof(sTensor), 4, tmp );
 fread( yield, sizeof(double), 4, tmp );
 
 /* Get element descriptor
 */
 data = (sQ4Data *)elm->data;

/* Get Material object
 */
 mat = MatList[data->matid-1];
 
 /* Get element lambda
 */
 lambda = MatList[data->matid-1]->Lambda;
 
 if(lambda>0.0){
  poropressure = _Q4PoroPressure(elm,lambda);
  /* Compute efective stress
  */
  for( i = 0; i < 4; i++ ) {
   estr[i].xx = stress[i].xx + poropressure;
   estr[i].yy = stress[i].yy + poropressure;
   estr[i].xy = stress[i].xy;
  } 
 }
 else{
  for( i = 0; i < 4; i++ ) {
   estr[i].xx = stress[i].xx;
   estr[i].yy = stress[i].yy;
   estr[i].xy = stress[i].xy;
  }
 }

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 4 );
 for( i = 0; i < 4; i++ ) {
  /* Compute principal stress and strain
  */
  PrincipalTensor( &stress[i], &pstress );
  PrincipalTensor( &strain[i], &pstrain );
  iso = (stress[i].xx + stress[i].yy + stress[i].zz)/3.0;
  /* Compute desviatory principal stress
  */
  piso = (pstress.dir1 + pstress.dir2 + pstress.dir3)/3.0;
  dpsts.dir1 = pstress.dir1 - piso;
  dpsts.dir3 = pstress.dir3 - piso;
  /* Computes the principal efetive stress
  */
  PrincipalTensor( &estr[i], &epsts );
  
  /* Write results on output file
  */
  
  for( j = 0; j < Config.numgaussresults; j++){
   if((strcmp( Config.gaussresults[j], "Yield" ) == 0))
    fprintf( out, "%+10.5e\t", yield[i] );
   else if((strcmp( Config.gaussresults[j], "Rf" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].rf );
   else if((strcmp( Config.gaussresults[j], "Sxx" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].xx );
   else if((strcmp( Config.gaussresults[j], "Syy" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].yy );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].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[i].xx );
   else if((strcmp( Config.gaussresults[j], "SEyy" ) == 0))
    fprintf( out, "%+10.5e\t", estr[i].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[i].xx );
   else if((strcmp( Config.gaussresults[j], "SDyy" ) == 0))
    fprintf( out, "%+10.5e\t", strd[i].yy );
   else if((strcmp( Config.gaussresults[j], "S1" ) == 0))
    fprintf( out, "%+10.5e\t", pstress.dir1 );
   else if((strcmp( Config.gaussresults[j], "S3" ) == 0))
    fprintf( out, "%+10.5e\t", pstress.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", strain[i].xx );
   else if((strcmp( Config.gaussresults[j], "Eyy" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].yy );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].xy );
   else if((strcmp( Config.gaussresults[j], "E1" ) == 0))
    fprintf( out, "%+10.5e\t", pstrain.dir1 );
   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) );
   else if((strcmp( Config.gaussresults[j], "K0" ) == 0)){
    if(estr[i].yy != 0.0){
     fprintf( out, "%+10.5e\t", (estr[i].xx/estr[i].yy) );
    }
    else{
     fprintf( out, "%+10.5e\t", 0.0 );
    }
   }
  }
 }

} /* End of Q4WriteGaussResult */


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

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 4 );
 for( i = 0; i < 4; i++ ) {
  PrincipalTensor( &stress[i], &pstress );
  PrincipalTensor( &strain[i], &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 Q4WriteGaussVectorResult */


/*
%F This function update the element stress to reflect the effects of
   the larges displacements.
*/
static void Q4UpdateStress( sElement *elm, double dtime, double *V,
                                                        sTensor *stre )
{
 sQ4Data *data = 0L;
 sCoord    coord[4];
 double    bm[3][8];
 double    ve[8];
 double    sxx, syy, rot;
 int       i, j, k, npg;

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

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

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

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

/* Update the element stress for each gauss point
 */
 npg = 0;
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Get element stress x strain matrix
   */
   _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

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

  /* Compute the element stress
   */
   sxx = stre[npg].xx;
   syy = stre[npg].yy;

   stre[npg].xx -= 2.0 * stre[npg].xy * rot;
   stre[npg].yy += 2.0 * stre[npg].xy * rot;
   stre[npg].xy += (sxx - syy) * rot;

  /* Update gauss point number
   */
   npg++;

  }
 }

} /* End of Q4UpdateStress */


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

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

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

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


/*
%F
*/
static void Q4PercForces( sElement *elm, double *pforce )
{
 sCoord     coord[4];
 sQ4Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     bm[3][8], phi[4], jac[2][2], shape[4];
 double     detjac;
 sTensor    grad;
 int        i, j, k;

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

/* Get element data
 */
 data = (sQ4Data *)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 < 8; i++ )
  pforce[i] = 0.0;

/* Compute percolation forces
 */
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Get element stress x strain matrix
   */
   _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

  /* Get element jacobiam
   */
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &detjac, jac );

   detjac *= _GWeight[i] * _GWeight[j] * mat->Rhof;

  /* Get element shape functions
   */
   _Q4ShapeFunc( _GPoint[i], _GPoint[j], shape );

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

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

  }
 }

} /* End of Q4PercForces */


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

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

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

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

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

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

} /* End of Q4SetInitStress */

/*
%F
*/
static void Q4ViscoForce( sElement *elm, double timeStep, sTensor *stress,
                          double *vforce )
{
 int        i, j, k, l, npg;
 double     coeff;
 double     bm[3][8], jac[2][2], cm[6][6];
 double     stav[3], strv[4], force[8];
 sTensor    vsta;
 sCoord     coord[4];
 sQ4Data   *data = 0L;
 sMaterial *mat  = 0L;

/* Get element descriptor
 */
 data = (sQ4Data *)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;

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

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

/* Compute visco force for each gauss point
 */
 npg = 0;
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
  /* Compute element B matrix
   */
   _Q4BMatrix( _GPoint[i], _GPoint[j], coord, bm );

  /* Compute rigidity coefficient
   */
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &coeff, jac );

   coeff *= _GWeight[i] * _GWeight[j];

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

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

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

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

  /* Add in the element vector
   */
   for( k = 0; k < (2*data->NumNodes); k++ )
    vforce[k] += force[k];

  /* Update gauss point number
   */
   npg++;

  }
 }

} /* End of Q4ViscoForce */

/*============================ Q4KMatrix ===========================*/

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

 for( s = 0; s < 2; s++ ) {
  for( r = 0; r < 2; r++ ) {
   /* Calcula matriz B*/
   _Q4BMatrix( _GPoint[r], _GPoint[s], coord, B );
   /* Calcula o determinante do jacobiano */
   _Q4Jacobian( _GPoint[r], _GPoint[s],  coord, &detjac, jac );
 
   /* Calcula [B] transposto * [C] */
   for( i = 0; i < 8; 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 < 8; i++ ) {
    for( j = 0; j < 8; 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 < 8; i++ ) {
    for( j = 0; j < 8; j++)  {
     k[i][j] += ke[i][j];
    }
   }
  }
 }
} /* Q4KMatrix */

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

/* Get element descriptor
 */
 data = (sQ4Data *)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 = 8;

} /* End of _Q4GetDof */

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

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

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

} /* End of Q4GetInc */


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

/*
%F This function initializes the sub-class "Q4".
*/
void Q4Init( void );
void Q4Init( void )
{
/* Initialize Q4 element sub-class
 */
 ElmClass[Q4].new               = Q4New;
 ElmClass[Q4].free              = Q4Free;
 ElmClass[Q4].read              = Q4Read;
 ElmClass[Q4].readinitstr       = Q4ReadInitStress;
 ElmClass[Q4].readprofile       = Q4ReadProfile;
 ElmClass[Q4].mass              = Q4MassMatrix;
 ElmClass[Q4].rigidcoeff        = Q4RigidCoeff;
 ElmClass[Q4].load              = Q4Load;
 ElmClass[Q4].connect           = Q4Connect;
 ElmClass[Q4].numnodes          = Q4NumNodes;
 ElmClass[Q4].gravity           = Q4Gravity;
 ElmClass[Q4].assvector         = Q4AssVector;
 ElmClass[Q4].strstrain         = Q4StressStrain;
 ElmClass[Q4].timestep          = Q4TimeStep;
 ElmClass[Q4].intforce          = Q4InterForce;
 ElmClass[Q4].writestr          = Q4WriteStress;
 ElmClass[Q4].writendlresult    = Q4WriteNodalResult;
 ElmClass[Q4].writegauresult    = Q4WriteGaussResult;
 ElmClass[Q4].writegauvecresult = Q4WriteGaussVectorResult;
 ElmClass[Q4].updatestress      = Q4UpdateStress;
 ElmClass[Q4].percforce         = Q4PercForces;
 ElmClass[Q4].setpressure       = Q4SetPressure;
 ElmClass[Q4].setinitstress     = Q4SetInitStress;
 ElmClass[Q4].viscoforce        = Q4ViscoForce;
 ElmClass[Q4].jacobian          = 0L;
 ElmClass[Q4].volume            = 0L;
 ElmClass[Q4].KMatrix           = Q4KMatrix;
 ElmClass[Q4].GetDof            = Q4GetDof;
 ElmClass[Q4].GetInc            = Q4GetInc;
 
} /* End of Q4Init */


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

---