brk8.c

Go to the documentation of this file.

/*
%M This modules contains the BRICK8 element sub-class methods and 
   definitions
%a Joao Luiz Elias Campos.
%d September 1st, 1998.
%r $Id: brk8.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:  23-Fev-06  Juan Pablo Ibañez
*    Modificada a funcão BRICK8ReadInitStress que passa a retornar um valor int.
* 
*
*  Modified:  Agosto-06  André Luis Muller
*    Criada a função _BRICK8ElementCenter que retorna as coordenadas do centro de
*    um elemento do tipo BRICK8.
*    Adicionado a função BRICK8StressStrain 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 BRICK8InterForce para possibilitar
*    a condição de tensões iniciais.
*
*    Modificação na função BRICK8New. Inicialização das tensões iniciais
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função BRICK8WriteGaussResults 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 BRICK8 element data definition
*/
typedef struct _brk8elm {
 int    matid;                /* Material identification          */
 int    tckid;                /* Element thickness identification */
 int    NumNodes;             /* Number of element nodes          */
 int    NumTensComp;          /* Number of tension components     */
 int    inc[8];               /* Element incidence                */
 double istr[8][6];           /* Element initial stress           */
 double effdef[8];            /* Effective deformation            */
 double prof;                 /* Element profile                  */
} sBrk8Data;


/*
%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[8][8] = {
                                 {  0.34150635, -0.09150635,
                                   -0.09150635, -0.52451905,
                                    0.77451905,  0.34150635,
                                    0.34150635, -0.09150635 },
                                 { -0.09150635,  0.34150635, 
                                   -0.52451905, -0.09150635,
                                    0.34150635,  0.77451905,
                                   -0.09150635,  0.34150635 },
                                 { -0.52451905, -0.09150635,
                                   -0.09150635,  0.34150635,
                                   -0.09150635,  0.34150635,
                                    0.34150635,  0.77451905 },
                                 { -0.09150635, -0.52451905,
                                    0.34150635, -0.09150635,
                                    0.34150635, -0.09150635,
                                    0.77451905,  0.34150635 },
                                 {  0.77451905,  0.34150635,
                                    0.34150635, -0.09150635,
                                    0.34150635, -0.09150635,
                                   -0.09150635, -0.52451905 },
                                 {  0.34150635,  0.77451905,
                                   -0.09150635,  0.34150635,
                                   -0.09150635,  0.34150635,
                                   -0.52451905, -0.09150635 },
                                 { -0.09150635,  0.34150635,
                                    0.34150635,  0.77451905,
                                   -0.52451905, -0.09150635,
                                   -0.09150635,  0.34150635 },
                                 {  0.34150635, -0.09150635,
                                    0.77451905,  0.34150635,
                                   -0.09150635, -0.52451905,
                                    0.34150635, -0.09150635 }
                                };

/*
%V Element edges index.
*/
static int _Edges[12][2] = { { 0, 1 }, { 1, 2 }, { 2, 3 }, { 3, 0 },
                             { 0, 4 }, { 1, 5 }, { 2, 6 }, { 3, 7 },
                             { 4, 5 }, { 5, 6 }, { 6, 7 }, { 7, 4 } };                           

static double GForce = 9.85;

/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _BRICK8GetCoord       ( sElement *, sCoord [] );
static void   _BRICK8ShapeFunc      ( double, double, double, double [] );
static void   _BRICK8Jacobian       ( double, double, double, sCoord [],
                                      double *, double [3][3] );
static void   _BRICK8BMatrix        ( double, double, double, sCoord [],
                                      double [6][24] );
static void   _BRICK8DerivShp       ( double, double, double, sDerivRST [] );
static void   _BRICK8GetDisplacement( sElement *, double *, double [] );
static double _BRICK8MinHeight      ( sCoord * );
static double _BRICK8CalcLength     ( sCoord *, sCoord * );
static void   _BRICK8ElementCenter  ( sElement *elm, sCoord coord[1] );

static int  _GetNol( sElement *elm, int noi, int noj, int nok, int type, double *v1x, double *v1y, double *v1z);
static void _Normal( sCoord coord[3], sCoord *normal );
static void _Q4ShapeFunc( double r, double s, double shape[4] );
static void _Q4Jacobian( double r, double s, sCoord coord[4], double *detjac,
                         double jac[2][2] );
static void _Q4DerivShp( double r, double s, sDerivRST dshp[4] );
static void EqLoad (sCoord F[4], sCoord coord[4], double v1x,
                    double v1y, double v1z);
static void GetCoordQ4(int i, sElement *elm, int noi, int noj,
                       int nok, int nol, sCoord coord[4]);
                    
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 */


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


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



static int _GetNol( sElement *elm, int noi, int noj, int nok, int type, double *v1x, double *v1y, double *v1z)
{
 sCoord coord[3];
 int    i,nol,node;
 sCoord normali;
 sCoord normalj;
 sBrk8Data *data = 0L;
 double p = (*v1x);
 sCoord vec;
 double prod;

/* Get element coordinates from first triangle
 */
 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;

 /*Get normal
 */
 _Normal( coord, &normali );
 
 /*Get nol from elm
 */
 data = (sBrk8Data *)elm->data;
 
 for( i = 0; i < data->NumNodes; i++ ) {
  node = data->inc[i];
  if((node!=noi)&&(node!=noj)&&(node!=nok)){
   coord[1].x = NodeVector[node-1].coord.x;
   coord[1].y = NodeVector[node-1].coord.y;
   coord[1].z = NodeVector[node-1].coord.z;
   _Normal( coord, &normalj );
   if((normali.x==normalj.x)&&
      (normali.y==normalj.y)&&
      (normali.z==normalj.z)){
    nol = node;
    break;
   }
  }
 }
 
 /* para garantir a direção da carga de pressão*/  
 
 for( i = 0; i < data->NumNodes; i++ ) {
  if((data->inc[i]!=noi)&&(data->inc[i]!=noj)
     &&(data->inc[i]!=nok)&&(data->inc[i]!=nol)){
   node = data->inc[i];
   break;
  }
 }
            
 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 * normali.x; 
 if( prod < 0.0)
  normali.x *= (-1.0);
  
 prod = vec.y * normali.y;
 if( prod < 0.0)
  normali.y *= (-1.0);
  
 prod = vec.z * normali.z;
 if( prod < 0.0)
  normali.z *= (-1.0);
    
 if(type == 1){
  *v1x = (p * (-normali.x));
  *v1y = (p * (-normali.y));
  *v1z = (p * (-normali.z));
 }
 
 return( nol );

} /* End of GetNol */


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

static void EqLoad (sCoord F[4], sCoord coord[4], double v1x,
                    double v1y, double v1z)
{
 int i,j,k;
 double shape[4];
 double jac[2][2], detjac;
 
 for( i = 0; i < 4; i++ ) {
  F[i].x = F[i].y = F[i].z = 0.0;
 }
 
 for( j = 0; j < 2; j++ ) {
  for( i = 0; i < 2; i++ ) {
   /* Calcula as funções de forma
   */
   _Q4ShapeFunc( _GPoint[i], _GPoint[j], shape );
   /* Calcula detjac
   */    
   _Q4Jacobian( _GPoint[i], _GPoint[j], coord, &detjac, jac );
   /* Calcula as forças nodais equivalentes
   */
   for( k = 0; k < 4; k++ ) {
    F[k].x+=shape[k]*v1x*fabs(detjac);
    F[k].y+=shape[k]*v1y*fabs(detjac);
    F[k].z+=shape[k]*v1z*fabs(detjac);
   }
  }
 }
} /* End EqLoad */



static void GetCoordQ4(int i, sElement *elm, int noi, int noj,
                       int nok, int nol, sCoord coord[4])
{
 switch(i){
 
  case(0):
   coord[0].x = elm->nodes[noi-1].coord.x;
   coord[0].y = elm->nodes[noi-1].coord.y;
   coord[1].x = elm->nodes[nok-1].coord.x;
   coord[1].y = elm->nodes[nok-1].coord.y;
   coord[2].x = elm->nodes[nol-1].coord.x;
   coord[2].y = elm->nodes[nol-1].coord.y;
   coord[3].x = elm->nodes[noj-1].coord.x;
   coord[3].y = elm->nodes[noj-1].coord.y;
  break;
  case(1):
   coord[0].x = elm->nodes[noi-1].coord.x;
   coord[0].y = elm->nodes[noi-1].coord.z;
   coord[1].x = elm->nodes[nok-1].coord.x;
   coord[1].y = elm->nodes[nok-1].coord.z;
   coord[2].x = elm->nodes[nol-1].coord.x;
   coord[2].y = elm->nodes[nol-1].coord.z;
   coord[3].x = elm->nodes[noj-1].coord.x;
   coord[3].y = elm->nodes[noj-1].coord.z;
  break;
  case(2):
   coord[0].x = elm->nodes[noi-1].coord.y;
   coord[0].y = elm->nodes[noi-1].coord.z;
   coord[1].x = elm->nodes[nok-1].coord.y;
   coord[1].y = elm->nodes[nok-1].coord.z;
   coord[2].x = elm->nodes[nol-1].coord.y;
   coord[2].y = elm->nodes[nol-1].coord.z;
   coord[3].x = elm->nodes[noj-1].coord.y;
   coord[3].y = elm->nodes[noj-1].coord.z;
  break;
 }
} /* End GetCoordQ4*/

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

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


/*
%F This function computes the length of a given segment.
*/
static double _BRICK8CalcLength( 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 BRICK8 element.
*/
static double _BRICK8MinHeight( sCoord *coord )
{
 double minlen;
 int    i;

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

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

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

 return minlen;

} /* End of _BRICK8MinHeight */


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

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


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

/* Initialize B matrix
 */
 for( i = 0; i < 24; 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
 */
 _BRICK8Jacobian( 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
 */
 _BRICK8DerivShp( r, s, t, dshp );

/* Compute the B matrix coefficientes
 */
 for( i = 0; i < 8; 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 _BRICK8BMatrix */


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

/* Compute the map derivatives
 */
 _BRICK8DerivShp( 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 < 8; 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 _BRICK8Jacobian */


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

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

/* Get element coordinates
 */
 for( i = 0; i < 8; 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 _BRICK8GetCoord */


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

} /* End of _BRICK8ShapeFunc */


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

 dshp[0].s = -0.125 * (1.0 - r) * (1.0 + t);
 dshp[1].s = -0.125 * (1.0 + r) * (1.0 + t);
 dshp[2].s =  0.125 * (1.0 + r) * (1.0 + t);
 dshp[3].s =  0.125 * (1.0 - r) * (1.0 + t);
 dshp[4].s = -0.125 * (1.0 - r) * (1.0 - t);
 dshp[5].s = -0.125 * (1.0 + r) * (1.0 - t);
 dshp[6].s =  0.125 * (1.0 + r) * (1.0 - t);
 dshp[7].s =  0.125 * (1.0 - r) * (1.0 - t);

 dshp[0].t =  0.125 * (1.0 - r) * (1.0 - s);
 dshp[1].t =  0.125 * (1.0 + r) * (1.0 - s);
 dshp[2].t =  0.125 * (1.0 + r) * (1.0 + s);
 dshp[3].t =  0.125 * (1.0 - r) * (1.0 + s);
 dshp[4].t = -0.125 * (1.0 - r) * (1.0 - s);
 dshp[5].t = -0.125 * (1.0 + r) * (1.0 - s);
 dshp[6].t = -0.125 * (1.0 + r) * (1.0 + s);
 dshp[7].t = -0.125 * (1.0 - r) * (1.0 + s);

} /* End of _BRICK8DerivShp */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   BRICK8New( int, int, int, int, sElement **, sElement **, sNode * );
static void   BRICK8Free                  ( sElement * );
static void   BRICK8Read                  ( sElement * );
static int    BRICK8ReadInitStress        ( sElement * );
static void   BRICK8ReadProfile           ( sElement * );
static void   BRICK8MassMatrix            ( sElement *, double * );
static void   BRICK8Gravity               ( sElement *, double *, double *,
                                                                  double * );
static double BRICK8RigidCoeff            ( sElement * );
static void   BRICK8Connect               ( sElement *, int * );
static void   BRICK8NumNodes              ( sElement *, int * );
static void   BRICK8AssVector             ( sElement *, double *, double * );
static void   BRICK8TimeStep              ( sElement *, double * );
static void   BRICK8InterForce            ( sElement *, sTensor *, double * );
static void   BRICK8StressStrain          ( sElement *, double, double *,
                                                        double *, sTensor *, 
                                                        sTensor * );
static void   BRICK8WriteStress           ( sElement *, FILE *, double *, 
                                                                double * );
static void   BRICK8WriteGaussResult      ( sElement *, FILE *, FILE * );
static void   BRICK8WriteNodalResult      ( sElement *, FILE *, FILE * );
static void   BRICK8WriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   BRICK8UpdateStress          ( sElement *, double, double *,
                                                                sTensor * );
static void   BRICK8PercForces            ( sElement *, double * );
static void   BRICK8SetPressure           ( sElement *, double );
static void   BRICK8SetInitStress         ( sElement *, sTensor * );
static void   BRICK8ViscoForce            ( sElement *, double, sTensor *, 
                                            double * );              
static void BRICK8KMatrix( sElement *, double [24][24] );                                                                         
static void BRICK8GetDof( sElement *, int [24], int * );
                       
/*
%F This function computes the equivalent load applied to the element 
   face.
*/
static void BRICK8Load( 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);
 int    i,type;
 sCoord coord[4], F[4];
 
 *nol = -1;

/* Calcula as forças nodais equivalentes
 */
 if( ltype == UNIFORM ) {
        
  /* Calcula o comprimento da aresta                                    
  */                                                                   
  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) || (ltype == AREAPRESSURE) ) {
 
  (*q1x) = (*q1y) = (*q1z) = 0.0;
  (*q2x) = (*q2y) = (*q2z) = 0.0;
  (*q3x) = (*q3y) = (*q3z) = 0.0;
  (*q4x) = (*q4y) = (*q4z) = 0.0;
  
  if( ltype == AREAUNIFORM ){
   type = 0;
  }
  else if(ltype == AREAPRESSURE ){
   type = 1;
  }
  /* Pega o nol e corrige carga se type = 1
  */    
  *nol = _GetNol(elm,noi,noj,nok,type,&v1x,&v1y,&v1z);
 
  for( i = 0; i < 3; i++ ){
   /* Pega as coordenadas da face
   */
   GetCoordQ4(i,elm,noi,noj,nok,*nol,coord);
   /* Calcula forças nodais equivalentes
   */
   EqLoad(F,coord,v1x,v1y,v1z);
  
   (*q1x) += F[0].x; 
   (*q1y) += F[0].y;
   (*q1z) += F[0].z;
   (*q2x) += F[1].x; 
   (*q2y) += F[1].y;
   (*q2z) += F[1].z;
   (*q3x) += F[2].x; 
   (*q3y) += F[2].y;
   (*q3z) += F[2].z;
   (*q4x) += F[3].x; 
   (*q4y) += F[3].y;
   (*q4z) += F[3].z;
  }
 }
}                       

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

/* Get memory for the element descriptor
 */
 (*elm) = (sElement *)malloc( sizeof(sElement) );

/* Get memory for the BRICK8 element data
 */
 data = (sBrk8Data *)malloc( sizeof(sBrk8Data) );

/* Fill BRICK8 element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 8;
 data->NumTensComp = 6;
 data->prof = 0.0;
 for( i = 0; i < 8; i++ )
  data->effdef[i] = 0.0;
 for( i = 0; i < 8; i++ ){
  data->istr[i][0] = 0.0;
  data->istr[i][1] = 0.0;
  data->istr[i][2] = 0.0;
  data->istr[i][3] = 0.0;
  data->istr[i][4] = 0.0;
  data->istr[i][5] = 0.0;
 }

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

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

} /* End of BRICK8New */


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

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

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

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

} /* End of BRICK8Free */


/*
%F This method reads the BRICK8 element information.
%i Element descriptor.
*/
static void BRICK8Read( sElement *elm )
{
 sBrk8Data *data = 0L;
 int      i1, i2, i3, i4, i5, i6, i7, i8;

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

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

/* Fill element information
 */
 data->inc[0] = i1;
 data->inc[1] = i2;
 data->inc[2] = i3;
 data->inc[3] = i4;
 data->inc[4] = i5;
 data->inc[5] = i6;
 data->inc[6] = i7;
 data->inc[7] = i8;

} /* End of BRICK8Read */


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

/* Get element data
 */
 data = (sBrk8Data *)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 != 8 ) {
  printf( "\n\nNumber of gauss points must be equal to eight.\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] = szz;
  data->istr[i][3] = sxy;
  data->istr[i][4] = sxz;
  data->istr[i][5] = syz;
 }

 return 1;

} /* End of BRICK8ReadInitStress */

/*
%F This method reads the BRICK8 profile
*/
static void BRICK8ReadProfile( sElement *elm )
{
 double prof;
 sBrk8Data *data = 0L;

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

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

} /* End of BRICK8ReadProfile */

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

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

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

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

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

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

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

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

   }
  }
 }
 
} /* End of BRICK8MassMatrix */


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

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

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

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

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

 return( volume );

} /* End of BRICK8RigidCoeff */


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

 data = (sBrk8Data *)elm->data;

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

} /* End of BRICK8Connect */


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

 (*nnodes) = data->NumNodes;

} /* End of BRICK8NumNodes */


/*
%F
*/
static void BRICK8AssVector( sElement *elm, double *GMatrix, double *matrix )
{
 int        i;
 int        conn[8];
 sBrk8Data *data;

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


/*
%F
*/


static void BRICK8StressStrain( sElement *elm, double dt, double *U,
                                               double *yield, sTensor *stre, 
                                               sTensor *stra )
{
 sBrk8Data  *data = 0L;
 sMaterial  *mat  = 0L;
 sCoord      coord[8];
 double      cm[6][6];
 double      bm[6][24];
 double      ue[24];
 double      def[6], str[6];
 double      pf = 1.0;
 int         i, j, k, l, m, npg;
 double     poropressure;
 double     lambda, gamma;
 double     distance;
 
 
/* Initialize stress and strain vector
 */
 for( i = 0; i < 8; 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 = (sBrk8Data *)elm->data;
 
/* Get Material object
 */
 mat = MatList[data->matid-1];
 
/* Get element coordinates
 */
 _BRICK8GetCoord( elm, coord );

/* Get element displacement
 */
 _BRICK8GetDisplacement( elm, U, ue );
 
/* 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 stress and strain for each gauss point
 */
 npg = 0;
 for( k = 0; k < 2; k++ ) {
  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
    */
    _BRICK8BMatrix( _GPoint[i], _GPoint[j], _GPoint[k], coord, bm );

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

   /* Compute the element stress
    */
    for( m = 0; m < 6; m++ ) {
     str[m] = 0.0;
     for( l = 0; l < 6; l++ )
      str[m] += cm[m][l] * def[l];
     str[m] += data->istr[npg][m];
    }
    
    if(lambda>0.0){
     /* Compute the poropressure
     */ 
     poropressure = lambda*gamma*GForce*distance*Config.loadfactor;
         /* 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; i++ ) {
      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].zz = str[2];
    stre[npg].xy = str[3];
    stre[npg].xz = str[4];
    stre[npg].yz = str[5];

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

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

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

   }
  }
 }

} /* End of BRICK8StressStrain */

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

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

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

/* Initialize element internal forces vector
 */
 for( i = 0; i < (3*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( k = 0; k < 2; k++ ) {
   for( j = 0; j < 2; j++ ) {
    for( i = 0; i < 2; i++ ) {
    /* Compute element B matrix
    */
     _BRICK8BMatrix( _GPoint[i], _GPoint[j], _GPoint[k], coord, bm );

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

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

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

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

   /* Get element stress
    */
    str[0] = stress[npg].xx;
    str[1] = stress[npg].yy;
    str[2] = stress[npg].zz;
    str[3] = stress[npg].xy;
    str[4] = stress[npg].xz;
    str[5] = stress[npg].yz;

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

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

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

   }
  }
  for( i = 0; i < (3*data->NumNodes); i++ )
   intforce[i] *= -1.0;
 }

} /* End of BRICK8InterForce */


/*
%F
*/
static void BRICK8WriteStress( sElement *elm, FILE *out, double *U, double *V )
{
 sTensor  stress[8];
 sTensor  strain[8];
 sTensor  strd[8];
 double   yield[8];
 int      i;
 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 < 8; 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
  */
  BRICK8StressStrain( elm, 0.0, U, yield, stress, strain );
 }

/* Compute desviatory stress
 */
 for( i = 0; i < 8; 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), 8, out );
 fwrite( strain, sizeof(sTensor), 8, out );
 fwrite( strd, sizeof(sTensor), 8, out );
 fwrite( yield, sizeof(double), 8, out );

} /* End of BRICK8WriteStress */


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

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

/* Compute the nodal results
 */
 for( i = 0; i < 8; 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 < 8; 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;
   strndl[i].xz  += _TRMatrix[i][j] * strgau[j].xz;
   strndl[i].yz  += _TRMatrix[i][j] * strgau[j].yz;
   strdndl[i].xx += _TRMatrix[i][j] * strdgau[j].xx;
   strdndl[i].yy += _TRMatrix[i][j] * strdgau[j].yy;
   strdndl[i].zz += _TRMatrix[i][j] * strdgau[j].zz;
   defndl[i].xx  += _TRMatrix[i][j] * defgau[j].xx;
   defndl[i].yy  += _TRMatrix[i][j] * defgau[j].yy;
   defndl[i].zz  += _TRMatrix[i][j] * defgau[j].zz;
   defndl[i].xy  += _TRMatrix[i][j] * defgau[j].xy;
   defndl[i].xz  += _TRMatrix[i][j] * defgau[j].xz;
   defndl[i].yz  += _TRMatrix[i][j] * defgau[j].yz;
   yieldn[i]     += _TRMatrix[i][j] * yieldg[j];
  }
 }

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


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

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 8 );
 for( i = 0; i < 8; i++ ) {
  PrincipalTensor( &stress[i], &pstress );
  PrincipalTensor( &strain[i], &pstrain );
  
  for( j = 0; j < Config.numgaussresults; j++){
   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], "Szz" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].zz );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].xy );
   else if((strcmp( Config.gaussresults[j], "Sxz" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].xz );
   else if((strcmp( Config.gaussresults[j], "Syz" ) == 0))
    fprintf( out, "%+10.5e\t", stress[i].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[i].xx );
   else if((strcmp( Config.gaussresults[j], "Sdy" ) == 0))
    fprintf( out, "%+10.5e\t", strd[i].yy );
   else if((strcmp( Config.gaussresults[j], "Sdz" ) == 0))
    fprintf( out, "%+10.5e\t", strd[i].zz );
   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], "Ezz" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].zz );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].xy );
   else if((strcmp( Config.gaussresults[j], "Exz" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].xz );
   else if((strcmp( Config.gaussresults[j], "Eyz" ) == 0))
    fprintf( out, "%+10.5e\t", strain[i].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 BRICK8WriteGaussResult */


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

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 8 );
 for( i = 0; i < 8; 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 BRICK8WriteGaussVectorResult */

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

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

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

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

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

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

   /* Compute rotation factor
    */
    rot = 0.0;
    for( l = 0; l < 24; 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 BRICK8UpdateStress */


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

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

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

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


/*
%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 BRICK8Gravity( sElement *elm, double *qx, double *qy, double *qz )
{
 sBrk8Data *data = 0;
 sCoord     coord[8];
 double     volume, mass, gamma;
 int        matid;

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

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

/* Compute element area
 */
 volume = BRICK8RigidCoeff( elm );

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

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

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

} /* End of BRICK8Gravity */


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

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

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

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

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

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

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

   /* Compute gradient
    */
    grad.xx = grad.yy = grad.zz = 0.0;
    for( k = 0; k < 8; 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 < 8; 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 BRICK8PercForce */


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

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

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

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

} /* End of BRICK8SetPressure */

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

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

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

} /* End of BRICK8SetInitStress */

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

/* Get element descriptor
 */
 data = (sBrk8Data *)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
 */
 _BRICK8GetCoord( elm, coord );

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

/* Compute visco force for each gauss point
 */
 npg = 0;
 for( k = 0; k < 2; k++ ) {
  for( j = 0; j < 2; j++ ) {
   for( i = 0; i < 2; i++ ) {
   /* Get element stress x strain matrix
    */
    _BRICK8BMatrix( _GPoint[i], _GPoint[j], _GPoint[k], coord, bm );

   /* Get element jacobiam
    */
    _BRICK8Jacobian( _GPoint[i], _GPoint[j], _GPoint[k], coord, &coeff, jac );

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

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

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

   /* 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 < (2*data->NumNodes); m++ ) {
     force[m] = 0.0;
     for( l = 0; l < 6; l++ )
      force[m] += bm[l][m] * strv[l];
     force[m] *= coeff;
    } 

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

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

   }
  }
 }

} /* End of BRICK8ViscoForce */


/*============================ BRICK8Kmatrix ===========================*/

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

 for( t = 0; t < 2; t++ ) {
  for( s = 0; s < 2; s++ ) {
   for( r = 0; r < 2; r++ ) {
    /* Calcula matriz B*/
    _BRICK8BMatrix( _GPoint[r], _GPoint[s], _GPoint[t], coord, B );
    /* Calcula o determinante do jacobiano */
    _BRICK8Jacobian( _GPoint[r], _GPoint[s], _GPoint[t], coord, &detjac, jac );
 
    /* Calcula [B] transposto * [C] */
    for( i = 0; i < 24; 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 < 24; i++ ) {
     for( j = 0; j < 24; j++)  {
      ke[i][j] =0.0;
      for(l=0; l<6; l++){
       ke[i][j] += BTC[i][l]*B[l][j]*detjac;
      }
     }
    }
    /* Calcula matriz de rigidez do elemento*/
    for( i = 0; i < 24; i++ ) {
     for( j = 0; j < 24; j++)  {
      k[i][j] += ke[i][j];
     }
    }
   }
  }
 }
} /* BRICK8KMatrix */

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

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

} /* End of _BRICK8GetDof */

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

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

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

} /* End of BRICK8GetInc */


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

/*
%F This function initializes the sub-class "BRICK8".
*/
void BRICK8Init( void );
void BRICK8Init( void )
{
/* Initialize BRICK8 element sub-class
 */
 ElmClass[BRICK8].new               = BRICK8New;
 ElmClass[BRICK8].free              = BRICK8Free;
 ElmClass[BRICK8].read              = BRICK8Read;
 ElmClass[BRICK8].readinitstr       = BRICK8ReadInitStress;
 ElmClass[BRICK8].readprofile       = BRICK8ReadProfile;
 ElmClass[BRICK8].mass              = BRICK8MassMatrix;
 ElmClass[BRICK8].rigidcoeff        = BRICK8RigidCoeff;
 ElmClass[BRICK8].load              = BRICK8Load;
 ElmClass[BRICK8].connect           = BRICK8Connect;
 ElmClass[BRICK8].numnodes          = BRICK8NumNodes;
 ElmClass[BRICK8].gravity           = BRICK8Gravity;
 ElmClass[BRICK8].assvector         = BRICK8AssVector;
 ElmClass[BRICK8].strstrain         = BRICK8StressStrain;
 ElmClass[BRICK8].timestep          = BRICK8TimeStep;
 ElmClass[BRICK8].intforce          = BRICK8InterForce;
 ElmClass[BRICK8].writestr          = BRICK8WriteStress;
 ElmClass[BRICK8].writendlresult    = BRICK8WriteNodalResult;
 ElmClass[BRICK8].writegauresult    = BRICK8WriteGaussResult;
 ElmClass[BRICK8].writegauvecresult = BRICK8WriteGaussVectorResult;
 ElmClass[BRICK8].updatestress      = BRICK8UpdateStress;
 ElmClass[BRICK8].percforce         = BRICK8PercForces;
 ElmClass[BRICK8].setpressure       = BRICK8SetPressure;
 ElmClass[BRICK8].setinitstress     = BRICK8SetInitStress;
 ElmClass[BRICK8].viscoforce        = BRICK8ViscoForce;
 ElmClass[BRICK8].jacobian          = 0L;
 ElmClass[BRICK8].volume            = 0L;
 ElmClass[BRICK8].KMatrix           = BRICK8KMatrix;
 ElmClass[BRICK8].GetDof            = BRICK8GetDof;
 ElmClass[BRICK8].GetInc            = BRICK8GetInc;

} /* End of BRICK8Init */


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

---