intf.c

Go to the documentation of this file.

/*
%M This modules contains the INTERFACE element sub-class methods and 
   definitions
%a Joao Luiz Elias Campos.
%d September 2nd, 1998.
%r $Id: intf.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 INTERFACEReadInitStress que passa a retornar 
*    um valor int.
*
*  Modified:  07-06 André Luis Muller
*    Modificação na função INTERFACEWriteGaussResults 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 INTERFACE element data definition
*/
typedef struct _itfelm {
 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];              /* Element initial stress           */
 double   effdef;               /* Effective deformation            */
 sTensor *iPress;               /* Element internal pressure        */
} sItfData;


/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _INTERFACEGetCoord       ( sElement *, sCoord [] );
static double _INTERFACEArea           ( sCoord [] );
static void   _INTERFACEBMatrix        ( sElement *, double [2][8] );
static void   _INTERFACEGetDisplacement( sElement *, double *, double [] );


/*
%F
*/
static void _INTERFACEBMatrix( sElement *elm, double bm[2][8] )
{
 sCoord coord[2];
 double ct, st, l;

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

/* Compute element lenght
 */
 l = _INTERFACEArea( coord );

/* Compute derivatives
 */
 ct = (coord[1].x - coord[0].x) / l;
 st = (coord[1].y - coord[0].y) / l;

/* Mount B matrix
 */
 bm[0][0] = bm[0][2] = bm[1][1] = bm[1][3] = -ct / 2.0;
 bm[0][1] = bm[0][3] = bm[1][4] = bm[1][6] = -st / 2.0;
 bm[0][4] = bm[0][6] = bm[1][5] = bm[1][7] =  ct / 2.0;
 bm[0][5] = bm[0][7] = bm[1][0] = bm[1][2] =  st / 2.0;

} /* End of _INTERFACEBMatrix */


/*
%F
*/
static void _INTERFACEGetDisplacement( sElement *elm, double *U, double ue[4] )
{
 sItfData *data = 0L;
 int       i;

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


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

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

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


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

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

 return( area );

} /* End of _INTERFACEArea */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   INTERFACENew( int, int, int, int, sElement **, sElement **, sNode * );
static void   INTERFACEFree                  ( sElement * );
static void   INTERFACERead                  ( sElement * );
static int    INTERFACEReadInitStress        ( sElement * );
static void   INTERFACEMassMatrix            ( sElement *, double * );
static double INTERFACERigidCoeff            ( sElement * );
static void   INTERFACEConnect               ( sElement *, int * );
static void   INTERFACENumNodes              ( sElement *, int * );
static void   INTERFACEAssVector             ( sElement *, double *, double * );
static void   INTERFACEStressStrain          ( sElement *, double, double *,
                                               double *, sTensor *, sTensor * );
static void   INTERFACEInterForce            ( sElement *, sTensor *, 
                                                           double * );
static void   INTERFACEWriteGaussResult      ( sElement *, FILE *, FILE * );
static void   INTERFACEWriteGaussVectorResult( sElement *, int, FILE *, FILE * );
static void   INTERFACEWriteNodalResult      ( sElement *, FILE *, FILE * );
static void   INTERFACEWriteStress           ( sElement *, FILE *, double *, 
                                                                   double * );
static void   INTERFACESetPressure           ( sElement *, double );
static void   INTERFACESetInitStress         ( sElement *, sTensor * );


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

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

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

/* Fill INTERFACE element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 4;
 data->NumTensComp = 2;
 data->effdef      = 0.0;
 data->iPress      = 0L;

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


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

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

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

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

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

} /* End of INTERFACEFree */


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

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

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

} /* End of INTERFACERead */


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

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

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

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

 return 1;

} /* End of INTERFACEReadInitStress */


/*
%F This method computes the mass matrix for the INTERFACE element.
%i Element descriptor.
%o Element mass matrix.
*/
static void INTERFACEMassMatrix( sElement *elm, double *mass )
{
 int i;

/* Compute the mass matrix
 */
 for( i = 0; i < 8; i++ )
  mass[i] = 0.0;
 
} /* End of INTERFACEMassMatrix */


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

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

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

 return( area );

} /* End of INTERFACERigidCoeff */


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

 data = (sItfData *)elm->data;

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

} /* End of INTERFACEConnect */


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

 (*nnodes) = data->NumNodes;

} /* End of INTERFACENumNodes */


/*
%F
*/
static void INTERFACEAssVector( sElement *elm, double *Gvec, double *vec )
{
 int       i;
 int       conn[4];
 sItfData *data;

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

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

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

} /* End of INTERFACEAssVector */


/*
%F
*/
static void INTERFACEStressStrain( sElement *elm, double dt, double *U,
                                                  double *yield, sTensor *stre,
                                                  sTensor *stra )
{
 sItfData  *data = 0L;
 sMaterial *mat  = 0L;
 double     cm[6][6];
 double     bm[2][8];
 double     ue[8];
 double     defe[2], str[3];
 double     pf = 1.0;
 int        i, j;

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

/* Get Material object
 */
 if( data->matid-1 >= 0 ){
  mat = MatList[data->matid-1];

 /* Update the material parameter
  */
  MatUpdateParameter( mat, data->effdef );
 
 /* Get material constitutive matrix
  */
  MatConstitutiveMatrix( mat, cm );
 
 /* Get element stress x strain matrix
  */
  _INTERFACEBMatrix( elm, bm );
 
 /* Get element displacement
  */
  _INTERFACEGetDisplacement( elm, U, ue );
 
 /* Compute the deformations
  */
  for( i = 0; i < 2; i++ ) {
   defe[i] = 0.0;
   for( j = 0; j < (2*data->NumNodes); j++ )
    defe[i] += bm[i][j] * ue[j];
  }

 /* Compute the element stress
  */
  for( i = 0; i < 2; i++ ) {
   str[i] = 0.0;
   for( j = 0; j < 2; j++ )
    str[i] += cm[i][j] * defe[j];
  }
  str[2] = 0.0;

 /* Correct stress based on the material model
  */
  MatUpdateStress( mat, dt, &pf, &data->effdef, str, defe );

 /* Set stress values 
  */
  stre->xx = str[0];
  stre->yy = str[1];
 
 /* Set stress values 
  */
  if( stra != 0L ) {
   stra->xx = defe[0];
   stra->yy = defe[1];
  }
 }
 else {
  pf = 0.0;
 /* Set yield parameter
 */
  if( yield != 0L )
   (*yield) = (pf > 0.0) ? 0.0 : 1.0;

  stre->xx = stre->yy = stre->xy = 0.0;
  stra->xx = stra->yy = stra->xy = 0.0;
 }
 
} /* End of INTERFACEStressStrain */


/*
%F This function computes the internal forces for the INTERFACE element
*/
static void INTERFACEInterForce( sElement *elm, sTensor *stress, 
                                                 double *intforce )
{
 int       i, j;
 double    coeff = 0.0;
 double    bm[2][8];
 double    str[2];
 sItfData *data = 0L;

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

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

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

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

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

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

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

} /* End of INTERFACEInterForce */


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

/* Check to see if the element was rezoned
 */
 if( elm->rezone == DONE ) {
  stress.xx = stress.xy = stress.xz =
              stress.yy = stress.yz =
                          stress.zz = 0.0;
  strd.xx = strd.xy = strd.xz =
            strd.yy = strd.yz =
                      strd.zz = 0.0;
  strain.xx = strain.xy = strain.xz =
              strain.yy = strain.yz =
                          strain.zz = 0.0;
  yield = 0.0;
 }
 else {
 /* Compute the element stress and strain
  */
  INTERFACEStressStrain( elm, 0.0, U, &yield, &stress, &strain );
 }

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

} /* End of INTERFACEWriteStress */


/*
%F
*/
static void INTERFACEWriteNodalResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor stress;
 sTensor strain;
 sTensor strd;
 double  yield;
 int     i;

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

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 for( i = 0; i < 6; i++ ) {
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t%+10.5e\t", stress.xx,
                                                        stress.yy,
                                                        stress.zz,
                                                        0.0 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", 0.0,
                                               0.0,
                                               0.0 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", 0.0,
                                               0.0,
                                               0.0 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", 0.0,
                                               0.0,
                                               0.0 );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", 0.0,
                                               0.0,
                                               0.0 );
  fprintf( out, "%+10.5e\n", yield );
 }

} /* End of INTERFACEWriteNodalResult */


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

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

/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 1 );
 
 for( j = 0; j < Config.numgaussresults; j++){
   if((strcmp( Config.gaussresults[j], "Yield" ) == 0))
    fprintf( out, "%+10.5e\t", yield );
   else if((strcmp( Config.gaussresults[j], "Rf" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "Sxx" ) == 0))
    fprintf( out, "%+10.5e\t", stress.xx );
   else if((strcmp( Config.gaussresults[j], "Syy" ) == 0))
    fprintf( out, "%+10.5e\t", stress.yy );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", stress.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", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SDxx" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SDyy" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "S1" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "S3" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   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", 0.0 );
   else if((strcmp( Config.gaussresults[j], "SD3" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "Exx" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "Eyy" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "E1" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "E3" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "Ev" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
   else if((strcmp( Config.gaussresults[j], "K0" ) == 0))
    fprintf( out, "%+10.5e\t", 0.0 );
 }                                                    

} /* End of INTERFACEWriteGaussResult */


/*
%F
*/
static void INTERFACEWriteGaussVectorResult( sElement *elm, int version,
                                             FILE *out, FILE *tmp )
{
 sTensor stress;
 sTensor strain;
 double  yield;

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

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

 if( version > 1 ) {
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
  fprintf( out, "%e\t%e\n", 0.0, 0.0 );
 }
 else {
  fprintf( out, "%e\t%e\t%e\t", 0.0, 0.0, 0.0 );
  fprintf( out, "%e\t%e\t%e\t", 0.0, 0.0, 0.0 );
  fprintf( out, "%e\t%e\t%e\t", 0.0, 0.0, 0.0 );
  fprintf( out, "%e\t%e\t%e\n", 0.0, 0.0, 0.0 );
 }

} /* End of INTERFACEWriteGaussVectorResult */


/*
%F
*/
static void INTERFACESetPressure( sElement *elm, double pot )
{
 sCoord    coord[2];
 sItfData *data = 0L;
 double    elev;

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

/* Compute element CG elevation
 */
 elev = 0.5 * (coord[0].y + coord[1].y);

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

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

} /* End of INTERFACESetPressure */

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

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

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

} /* End of INTERFACESetInitStress */

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

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

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

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

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

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

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

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

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

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

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

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

} /* End of INTERFACEViscoForce */

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

/*
%F This function initializes the sub-class "INTERFACE".
*/
void INTERFACEInit( void );
void INTERFACEInit( void )
{
/* Initialize INTERFACE element sub-class
 */
 ElmClass[INTERFACE].new               = INTERFACENew;
 ElmClass[INTERFACE].free              = INTERFACEFree;
 ElmClass[INTERFACE].read              = INTERFACERead;
 ElmClass[INTERFACE].readinitstr       = INTERFACEReadInitStress;
 ElmClass[INTERFACE].mass              = INTERFACEMassMatrix;
 ElmClass[INTERFACE].rigidcoeff        = INTERFACERigidCoeff;
 ElmClass[INTERFACE].load              = 0L;
 ElmClass[INTERFACE].connect           = INTERFACEConnect;
 ElmClass[INTERFACE].numnodes          = INTERFACENumNodes;
 ElmClass[INTERFACE].gravity           = 0L;
 ElmClass[INTERFACE].assvector         = INTERFACEAssVector;
 ElmClass[INTERFACE].strstrain         = INTERFACEStressStrain;
 ElmClass[INTERFACE].timestep          = 0L;
 ElmClass[INTERFACE].intforce          = INTERFACEInterForce;
 ElmClass[INTERFACE].writestr          = INTERFACEWriteStress;
 ElmClass[INTERFACE].writendlresult    = INTERFACEWriteNodalResult;
 ElmClass[INTERFACE].writegauresult    = INTERFACEWriteGaussResult;
 ElmClass[INTERFACE].writegauvecresult = INTERFACEWriteGaussVectorResult;
 ElmClass[INTERFACE].updatestress      = 0L;
 ElmClass[INTERFACE].percforce         = 0L;
 ElmClass[INTERFACE].setpressure       = INTERFACESetPressure;
 ElmClass[INTERFACE].setinitstress     = INTERFACESetInitStress;
 ElmClass[INTERFACE].viscoforce        = INTERFACEViscoForce;
 ElmClass[INTERFACE].jacobian          = 0L;
 ElmClass[INTERFACE].volume            = 0L;
 ElmClass[INTERFACE].KMatrix           = 0L;
 ElmClass[INTERFACE].GetDof            = 0L;

} /* End of INTERFACEInit */


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

---