inf.c

Go to the documentation of this file.

/*
%M This modules contains the INFINITE element sub-class methods and 
   definitions
%a Joao Luiz Elias Campos.
%d September 1st, 1998.
%r $Id: inf.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 INFINITEReadInitStress que passa a retornar um valor int.
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função INFINITEWriteGaussResults 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 "alg.h"

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


/*
%T INFINITE sub-class definition
*/
typedef enum _inftype
{
 INF_T3,                         /* Infinite element with 3 nodes    */
 INF_Q4,                         /* Infinite element with 4 nodes    */
 NumInfTypes                     /* Number of infinite sub-classes   */
} eInfType;


/*
%T INFINITE element sub-class definition
*/
typedef struct _infclass
{
 void (*mapfunc) ( double, double, double [] );
 void (*mapderiv)( double, double, sDerivRST [] );
 void (*shpderiv)( double, double, sDerivRST [] );
} sInfClass;


/*
%T INFINITE element data definition
*/
typedef struct _infelm
{
 eInfType type;                  /* Infinite element type            */
 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            */
} sInfData;


/*
%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.683012702,  0.683012702, 
                                   -0.183012702, -0.183012702 },
                                 { -0.183012702, -0.183012702, 
                                    0.683012702,  0.683012702 },
                                 {  0.250000000, -0.616025404, 
                                    1.116025404,  0.250000000 }
                                };


/*
%V Infiniti element sub-class
*/
static sInfClass InfClass[NumInfTypes];


/*
%K Macros definitions
*/
#define INFINITEMappFunc(d,r,s,map)                     \
        if(InfClass[d->type].mapfunc != 0L)             \
         (*InfClass[d->type].mapfunc)( r, s, map )

#define INFINITEDerivMap(d,r,s,dmap)                    \
        if(InfClass[d->type].mapderiv != 0L)            \
         (*InfClass[d->type].mapderiv)( r, s, dmap )

#define INFINITEDerivShp(d,r,s,dshp)                    \
        if(InfClass[d->type].shpderiv != 0L)            \
         (*InfClass[d->type].shpderiv)( r, s, dshp )


/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _INFINITEGetCoord       ( sElement *, sCoord [] );
static void   _INFINITEJacobian       ( sInfData *, double, double, sCoord [],
                                        double *, double [2][2] );
static void   _INFINITEBMatrix        ( sInfData *, double, double, sCoord [], 
                                        double [3][8] );
static void   _INFINITEGetDisplacement( sElement *, double *, double [] );
static double _INFINITEMinHeight      ( sInfData *, sCoord * );


/*
%F This function computes the minimun height for the INFINITE element.
*/
static double _INFINITEMinHeight( sInfData *data, sCoord *coord )
{
 double jac[2][2], len[2];
 double detjac, sum, minlen = 1.0e+8;
 int    i, j, k, l;

/* For each gauss point
 */
 for( j = 0; j < 2; j++ )
 {
  for( i = 0; i < 2; i++ )
  {
  /* Compute the element jacobiam matrix
   */
   _INFINITEJacobian( data, _GPoint[i], _GPoint[j], coord, &detjac, jac );

  /* Compute the characteristic lenght
   */
   for( k = 0; k < 2; k++ )
   {
    sum = 0.0;
    for( l = 0; l < 2; l++ )
     sum += jac[l][k] * jac[k][l];
    len[k] = 2.0 * sum;
   }

  /* Get the minimun element lenght
   */
   if( minlen > len[0] ) minlen = len[0];
   if( minlen > len[1] ) minlen = len[1];
  }
 }

 return( minlen );

} /* End of _INFINITEMinHeight */


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

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


/*
%F This function computes the stress x strain matrix for the INFINITE element.
*/
static void _INFINITEBMatrix( sInfData *data, 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
 */
 _INFINITEJacobian( data, 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
 */
 INFINITEDerivShp( data, r, s, dshp );

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


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

/* Compute the map derivatives
 */
 INFINITEDerivMap( data, 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 < data->NumNodes; 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 _INFINITEJacobian */


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

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

/* Get element coordinates
 */
 for( i = 0; i < data->NumNodes; 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 _INFINITEGetCoord */


/*
** ------------------------------------------------------------------------
** Infinite subclasses methods:
*/
static void INFT3MappFunc  ( double, double, double [] );
static void INFT3MappDeriv ( double, double, sDerivRST [] );
static void INFT3ShapeDeriv( double, double, sDerivRST [] );
static void INFQ4MappFunc  ( double, double, double [] );
static void INFQ4MappDeriv ( double, double, sDerivRST [] );
static void INFQ4ShapeDeriv( double, double, sDerivRST [] );


/*
%F This function computes the mapp functions for the triangular 
   infinite element.
*/
static void INFT3MappFunc( double r, double s, double shape[4] )
{
/* Compute map functions
 */
 shape[0] = (r * s + 3.0 * (-1.0 - r - s)) / ((1.0 - r) * (1.0 - s));
 shape[1] = 2.0 * (1.0 + r) / ((1.0 - r) * (1.0 - s ));
 shape[2] = 2.0 * (1.0 + s) / ((1.0 - r) * (1.0 - s ));

} /* End of INFT3MappFunc */


/*
%F This function computes the mapp functions derivatives for the 
   triangular infinite element.
*/
static void INFT3MappDeriv( double r, double s, sDerivRST dshp[4] )
{
/* Compute map functions derivatives
 */
 dshp[0].r =  2.0 * (s + 3.0) / ((r - 1.0) * (r - 1.0) * (s - 1.0));
 dshp[1].r = -4.0 / ((r - 1.0) * (r - 1.0) * (s - 1.0));
 dshp[2].r = -2.0 * (s + 1.0) / ((r - 1.0) * (r - 1.0) * (s - 1.0));

 dshp[0].s =  2.0 * (r + 3.0) / ((r - 1.0) * (s - 1.0) * (s - 1.0));
 dshp[1].s = -2.0 * (r + 1.0) / ((r - 1.0) * (s - 1.0) * (s - 1.0));
 dshp[2].s = -4.0 / ((r - 1.0) * (s - 1.0) * (s - 1.0));

} /* End of INFT3MappDeriv */


/*
%F This function computes the shape functions derivatives for the 
   triangular infinite element.
*/
static void INFT3ShapeDeriv( double r, double s, sDerivRST dshp[4] )
{
/* Compute shape functions derivatives
 */
 dshp[0].r = (1.0 - s) * (2.0 * r + s) / 4.0;
 dshp[1].r = r * (s - 1.0);
 dshp[2].r = (s * s - 1.0) / 2.0;

 dshp[0].s = (1.0 - r) * (2.0 * s + r) / 4.0;
 dshp[1].s = (r * r - 1.0) / 2.0;
 dshp[2].s = s * (r - 1.0);

} /* End of INFT3ShapeDeriv */


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

} /* End of INFQ4MappFunc */


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

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

} /* End of INFQ4MappDeriv */


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

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

} /* End of INFQ4ShapeDeriv */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   INFINITENew( int, int, int, int, sElement **, sElement **, sNode * );
static void   INFINITEFree                  ( sElement * );
static void   INFINITERead                  ( sElement * );
static int    INFINITEReadInitStress        ( sElement * );
static double INFINITERigidCoeff            ( sElement * );
static void   INFINITEMassMatrix            ( sElement *, double * );
static void   INFINITEConnect               ( sElement *, int * );
static void   INFINITENumNodes              ( sElement *, int * );
static void   INFINITEAssVector             ( sElement *, double *, double * );
static void   INFINITETimeStep              ( sElement *, double * );
static void   INFINITEInterForce            ( sElement *, sTensor *, double [] );
static void   INFINITEStressStrain          ( sElement *, double, double *, 
                                              double *, sTensor *, sTensor * );
static void   INFINITEWriteStress           ( sElement *, FILE *, double *, 
                                                                  double * );
static void   INFINITEWriteGaussResult      ( sElement *, FILE *, FILE * );
static void   INFINITEWriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   INFINITEUpdateStress          ( sElement *, double, double *, 
                                                                  sTensor * );
static void   INFINITESetPressure           ( sElement *, double );
static void   INFINITESetInitStress         ( sElement *, sTensor * );
static void   INFINITEViscoForce            ( sElement *, double, sTensor *,
                                              double * );


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

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

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

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

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


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

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

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

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

} /* End of INFINITEFree */


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

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

/* Read the element incidence
 */
 fscanf( nf, "%d %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;

/* If necessary, update the infinite element type
 */
 if( i4 == 0 )
 {
  data->type     = INF_T3;
  data->NumNodes = 3;
 }

} /* End of INFINITERead */


/*
%F This method reads the INIFINIT initial stress
*/
static int INFINITEReadInitStress( sElement *elm )
{
 int       i, numpg;
 double    sxx, syy, sxy;
 sInfData *data = 0L;

/* Get element data
 */
 data = (sInfData *)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
 */
 for( i = 0; i < numpg; i++ )
 {
  fscanf( nf, "%lf %lf %lf", &sxx, &syy, &sxy );
  data->istr[i][0] = sxx;
  data->istr[i][1] = syy;
  data->istr[i][2] = sxy;
 }

 return 1;

} /* End of INFINITEReadInitStress */


/*
%F This method computes the mass matrix for the INFINITE element.
%i Element descriptor.
%o Element mass matrix.
*/
static void INFINITEMassMatrix( sElement *elm, double *mass )
{
 sInfData *data = 0L;
 sCoord    coord[4];
 double    mapp[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
 */
 _INFINITEGetCoord( elm, coord );

/* Get material density
 */
 data = (sInfData *)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 and mapp functions
   */
   INFINITEMappFunc( data, _GPoint[i], _GPoint[j], mapp );

  /* Compute the jacobian matrix
   */
   _INFINITEJacobian( data, _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 * mapp[k] * mapp[l];
     mass[2*k+1] += coeff * mapp[k] * mapp[l];
    }
   }

  }
 }
 
 for( i = 0; i < 8; i++ )
  mass[i] = fabs(mass[i]);

} /* End of INFINITEMassMatrix */


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

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

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

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

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

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

 return( area );

} /* End of Q4RigidCoeff */


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

 data = (sInfData *)elm->data;

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

} /* End of INFINITEConnect */


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

 (*nnodes) = data->NumNodes;

} /* End of INFINITENumNodes */


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

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


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

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

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

/* Compute the element minimum height
 */
 hmin = _INFINITEMinHeight( data, coord );

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


/*
%F
*/
static void INFINITEStressStrain( sElement *elm, double dt, double *U,
                                                 double *yield, sTensor *stre,
                                                 sTensor *stra )
{
 sInfData   *data = 0L;
 sMaterial  *mat  = 0L;
 sCoord      coord[4];
 double      cm[6][6];
 double      bm[3][8];
 double      ue[8];
 double      def[3], str[4];
 double      pf = 1.0;
 int         i, j, k, l, npg;
 
/* Get element descriptor
 */
 data = (sInfData *)elm->data;
 
/* Get Material object
 */
 mat = MatList[data->matid-1];
 
/* Get element coordinates
 */
 _INFINITEGetCoord( elm, coord );

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

/* 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
   */
   _INFINITEBMatrix( data, _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[3] = 0.0;

  /* 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].zz = str[3];

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


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

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

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

/* Initialize element internal force vector
 */
 for( i = 0; i < 8; i++ )
  intforce[i] = 0.0;

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

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

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

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

  /* 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 < 8; k++ )
    intforce[k] += iforce[k];

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

  }
 }

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


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

/* 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
  */
  INFINITEStressStrain( 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 INFINITEWriteStress */


/*
%F
*/
static void INFINITEWriteNodalResult( 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;

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

/* 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,
                                               0.0 );
  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 INFINITEWriteNodalResult */


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

/* 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 );
  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;
  
  /* 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", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SExx" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SEyy" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   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", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SE3" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   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))
    fprintf( out, "%+10.5e\t", 0.0 );
  }
 }

} /* End of INFINITEWriteGaussResult */


/*
%F
*/
static void INFINITEWriteGaussVectorResult( 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.cos1x),
                                                  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.cos1x),
                                                  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 INFINITEWriteGaussVectorResult */


/*
%F This function update the element stress to reflect the effects of
   the larges displacements.
*/
static void INFINITEUpdateStress( sElement *elm, double dtime, double *V,
                                    sTensor *stre )
{
 sInfData *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 = (sInfData *)elm->data;

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

/* Get element displacement increment
 */
 _INFINITEGetDisplacement( 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
   */
   _INFINITEBMatrix( data, _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 INFINITEUpdateStress */


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

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

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

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

} /* End of INFINITESetPressure */

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

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

/*
%F
*/
static void INFINITEViscoForce( 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];
 sInfData  *data = 0L;
 sMaterial *mat  = 0L;

/* Get element descriptor
 */
 data = (sInfData *)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
 */
 _INFINITEGetCoord( 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
   */
   _INFINITEBMatrix( data, _GPoint[i], _GPoint[j], coord, bm );

  /* Compute rigidity coefficient
   */
   _INFINITEJacobian( data, _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 INFINITEViscoForce */

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

/*
%F This function initializes the sub-class "INFINITE".
*/
void INFINITEInit( void );
void INFINITEInit( void )
{
/* Initialize INFINITE element sub-class
 */
 ElmClass[INFINITE].new               = INFINITENew;
 ElmClass[INFINITE].free              = INFINITEFree;
 ElmClass[INFINITE].read              = INFINITERead;
 ElmClass[INFINITE].readinitstr       = INFINITEReadInitStress;
 ElmClass[INFINITE].mass              = INFINITEMassMatrix;
 ElmClass[INFINITE].rigidcoeff        = INFINITERigidCoeff;
 ElmClass[INFINITE].load              = 0L;
 ElmClass[INFINITE].connect           = INFINITEConnect;
 ElmClass[INFINITE].numnodes          = INFINITENumNodes;
 ElmClass[INFINITE].gravity           = 0L;
 ElmClass[INFINITE].assvector         = INFINITEAssVector;
 ElmClass[INFINITE].strstrain         = INFINITEStressStrain;
 ElmClass[INFINITE].timestep          = 0L;
 ElmClass[INFINITE].intforce          = INFINITEInterForce;
 ElmClass[INFINITE].writestr          = INFINITEWriteStress;
 ElmClass[INFINITE].writendlresult    = INFINITEWriteNodalResult;
 ElmClass[INFINITE].writegauresult    = INFINITEWriteGaussResult;
 ElmClass[INFINITE].writegauvecresult = INFINITEWriteGaussVectorResult;
 ElmClass[INFINITE].updatestress      = INFINITEUpdateStress;
 ElmClass[INFINITE].percforce         = 0L;
 ElmClass[INFINITE].setpressure       = INFINITESetPressure;
 ElmClass[INFINITE].setinitstress     = INFINITESetInitStress;
 ElmClass[INFINITE].viscoforce        = INFINITEViscoForce;
 ElmClass[INFINITE].jacobian          = 0L;
 ElmClass[INFINITE].volume            = 0L;
 ElmClass[INFINITE].KMatrix           = 0L;
 ElmClass[INFINITE].GetDof            = 0L;

/* Define the Infinite subclasses initialization methods
 */
 InfClass[INF_T3].mapfunc  = INFT3MappFunc;
 InfClass[INF_T3].mapderiv = INFT3MappDeriv;
 InfClass[INF_T3].shpderiv = INFT3ShapeDeriv;

 InfClass[INF_Q4].mapfunc  = INFQ4MappFunc;
 InfClass[INF_Q4].mapderiv = INFQ4MappDeriv;
 InfClass[INF_Q4].shpderiv = INFQ4ShapeDeriv;

} /* End of INFINITEInit */


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

---