line2.c

Go to the documentation of this file.

/*
 * $Id: line2.c,v 1.4 2004/06/24 18:32:35 joaoluiz Exp $
 * (C) COPYRIGHT 2004, 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:    01-Sep-05    Alexandre A. Del Savio
 *     Criada as funcoes LINE2SetConnect e LINE2UpdateConnect para
 *     serem usadas quando os segmentos de uma linha foram estendidos.
 *
 *  Modified:  23-Fev-06  Juan Pablo Ibañez
 *    Modificada a funcão LINE2ReadInitStress que passa a retornar 
 *    um valor int.
 *     
 *
 *  Modified:  14-Dez-06  Gisele Holtz
 *    Criada a função _LINE2Normalize e modificada a função 
 *    LINE2Jacobian para normalizar os vetores do sistema de 
 *    coordenadas local.
 */

/*
 * Global variables and symbols:
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "rio.h"
#include "load.h"
#include "elm.h"
#include "node.h"
#include "material.h"

/*
 * Local variables and symbols:
 */

typedef struct _line2elm {
 int matid;
 int tckid;
 int NumNodes;
 int NumTensComp;
 int inc[2];
 double istr[4];
 double effdef;
 sTensor *iPress;
} sLINE2Data;

/*
 * Local functions:
 */
static void _LINE2GetCoord(sElement *, sCoord []);
static double _LINE2Area(sCoord []);
static double _LINE2MinHeight(sCoord []);
static void _LINE2GetDisplacement(sElement *, double *, double []);
static void _LINE2BMatrix(sElement *, double [3][6]);
static double _LINE2Lenght(sCoord *);
static sCoord _LINE2Cross(sCoord *, sCoord *);
static void _LINE2Inverse(double [3][3], double [3][3]);
static void  _LINE2Normalize( sCoord *vec, sCoord *n );

static void  _LINE2Normalize( sCoord *vec, sCoord *n )
{
  double      len;

  len = _LINE2Lenght(vec);

  n->x = vec->x / len;
  n->y = vec->y / len;
  n->z = vec->z / len;
}

static void _LINE2Inverse(double m[3][3], double d[3][3])
{
 double det = m[0][0] * ((m[1][1] * m[2][2]) - (m[1][2] * m[2][1])) -
              m[0][1] * ((m[1][0] * m[2][2]) - (m[1][2] * m[2][0])) +
              m[0][2] * ((m[1][0] * m[2][1]) - (m[1][1] * m[2][0]));

 d[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
 d[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]) / det;
 d[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
 d[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]) / det;
 d[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
 d[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]) / det;
 d[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
 d[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]) / det;
 d[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;

} /* End of _LINE2Inverse */

static double _LINE2Lenght(sCoord *v)
{
 return sqrt((v->x * v->x) + (v->y * v->y) + (v->z * v->z));

} /* End of _LINE2Lenght */

static sCoord _LINE2Cross(sCoord *u, sCoord *v)
{
 sCoord w;
 w.x = (u->y * v->z) - (u->z * v->y);
 w.y = (u->z * v->x) - (u->x * v->z);
 w.z = (u->x * v->y) - (u->y * v->x);
 return w;

} /* End of _LINE2Cross */

static void _LINE2GetCoord(sElement *elm, sCoord coord[2])
{
 sLINE2Data *data = 0L;

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

/* Get element coordinates
 */
 coord[0].x = elm->nodes[data->inc[0]-1].coord.x;
 coord[0].y = elm->nodes[data->inc[0]-1].coord.y;
 coord[0].z = elm->nodes[data->inc[0]-1].coord.z;
 coord[1].x = elm->nodes[data->inc[1]-1].coord.x;
 coord[1].y = elm->nodes[data->inc[1]-1].coord.y;
 coord[1].z = elm->nodes[data->inc[1]-1].coord.z;

} /* End of _LINE2GetCoord */

static double _LINE2Area(sCoord coord[2])
{
 double dx = coord[1].x - coord[0].x;
 double dy = coord[1].y - coord[0].y;
 double dz = coord[1].z - coord[0].z;
 return sqrt((dx * dx) + (dy * dy) + (dz * dz));

} /* End of _LINE2Area */

static double _LINE2MinHeight(sCoord coord[2])
{
 return _LINE2Area(coord);

} /* End of _LINE2MinHight */

static void _LINE2GetDisplacement(sElement *elm, double *U, double ue[6])
{
 sLINE2Data *data = 0L;

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

/* Fill element displacement vector
 */
 ue[0] = U[2*(data->inc[0]-1)];
 ue[1] = U[2*(data->inc[0]-1)+1];
 ue[2] = U[2*(data->inc[1]-1)];
 ue[3] = U[2*(data->inc[1]-1)+1];

} /* End of _LINE2GetDisplacement */

static void _LINE2BMatrix(sElement *elm, double bm[3][6])
{
 double jac[3][3], ijac[3][3];
 double len;
 int i;
 sCoord n1, n2, n3;
 sCoord ex = {1.0, 0.0, 0.0};
 sCoord ey = {0.0, 1.0, 0.0};
 sCoord ez = {0.0, 0.0, 1.0};
 sCoord coord[2], dxyz[2];

 _LINE2GetCoord(elm, coord);

 n1.x = 0.5 * (coord[1].x - coord[0].x);
 n1.y = 0.5 * (coord[1].y - coord[0].y);
 n1.z = 0.5 * (coord[1].z - coord[0].z);

 len = _LINE2Lenght(&n1);
 n2 = _LINE2Cross(&ex, &n1);
 if( (_LINE2Lenght(&n2)/len) < 1.0e-5 ) {
  n2 = _LINE2Cross(&ey, &n1);
  if( (_LINE2Lenght(&n2)/len) < 1.0e-5 ) {
   n2 = _LINE2Cross(&ez, &n1);
  }
 }
 n3 = _LINE2Cross(&n1, &n2);

 jac[0][0] = n1.x;
 jac[0][1] = n1.y;
 jac[0][2] = n1.z;
 jac[1][0] = n2.x;
 jac[1][1] = n2.y;
 jac[1][2] = n2.z;
 jac[2][0] = n3.x;
 jac[2][1] = n3.y;
 jac[2][2] = n3.z;

 _LINE2Inverse(jac, ijac);

 dxyz[0].x = -0.5 * ijac[0][0];
 dxyz[0].y = -0.5 * ijac[1][0];
 dxyz[0].z = -0.5 * ijac[2][0];
 dxyz[1].x =  0.5 * ijac[0][0];
 dxyz[1].y =  0.5 * ijac[1][0];
 dxyz[1].z =  0.5 * ijac[2][0];

 for( i = 0; i < 2; i++ ) {
  bm[0][2*i] = bm[3][2*i+1] = dxyz[i].x;
  bm[1][2*i+1] = bm[3][2*i] = dxyz[i].y;
  bm[2][2*i+1] = bm[2][2*i] = 0.0;
  bm[0][2*i+1] = bm[1][2*i] = 0.0;
 }
 
} /* End of _LINE2BMatrix */

/*
 * Subclass methods:
 */
static void LINE2New(int, int, int, int, sElement **, sElement **, sNode * );
static void LINE2Free(sElement * );
static void LINE2Read(sElement * );
static int  LINE2ReadInitStress(sElement * );
static void LINE2MassMatrix(sElement *, double * );
static double LINE2RigidCoeff(sElement *);
static void LINE2Connect(sElement *, int *);
static void LINE2SetConnect(sElement *, int *);
static void LINE2UpdateConnect(sElement *);
static void LINE2NumNodes(sElement *, int *);
static void LINE2Gravity(sElement *, double *, double *, double *);
static void LINE2AssVector(sElement *, double *, double *);
static void LINE2StressStrain(sElement *, double, double *, double *,
                              sTensor *, sTensor *);
static void LINE2TimeStep(sElement *, double *);
static void LINE2InterForce(sElement *, sTensor *, double *);
static void LINE2WriteStress(sElement *, FILE *, double *, double *);
static void LINE2WriteNodalResult(sElement *, FILE *, FILE *);
static void LINE2WriteGaussResult(sElement *, FILE *, FILE *);
static void LINE2WriteGaussVectorResult(sElement *, int, FILE *, FILE *);
static void LINE2UpdateStress(sElement *, double, double *, sTensor *);
static void LINE2PercForces(sElement *, double *);
static void LINE2SetPressure(sElement *, double);
static void LINE2SetInitStress(sElement *, sTensor *);
static void LINE2ViscoForce(sElement *, double, sTensor *, double *);
static void LINE2Jacobian(sElement *, double [3][3], double [3][3]);
static void LINE2Volume(sElement *, double *);

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

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

/* Get memory for the LINE2 element data
 */
 data = (sLINE2Data *)malloc( sizeof(sLINE2Data) );

/* Fill LINE2 element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 2;
 data->NumTensComp = 4;
 data->effdef      = 0.0;
 data->istr[0]     = 0.0;
 data->istr[1]     = 0.0;
 data->istr[2]     = 0.0;
 data->iPress      = 0L;

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


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

/* Get LINE2 element data
 */
 data = (sLINE2Data *)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 LINE2Free */


/*
%F This method reads the LINE2 element information.
%i Element descriptor.
*/
static void LINE2Read( sElement *elm )
{
 sLINE2Data *data = 0L;
 int      i1, i2;

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

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

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

} /* End of LINE2Read */


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

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


/*
%F This method computes the mass matrix for the LINE2 element.
%i Element descriptor.
%o Element mass matrix.
*/
static void LINE2MassMatrix( sElement *elm, double *mass )
{
 sLINE2Data *data = 0L;
 sCoord   coord[2];
 double   area, ncont, dens;
 int      i, matid;

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

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

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

/* Calculate the element node contribution for the element mass
 */
 ncont = (area * dens) / 2.0;

/* Compute the mass matrix
 */
 for( i = 0; i < 6; i++ )
  mass[i] = ncont;
 
} /* End of LINE2MassMatrix */


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

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

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

 return( area );

} /* End of LINE2RigidCoeff */

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

 data = (sLINE2Data *)elm->data;

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

} /* End of LINE2Connect */

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

 data = (sLINE2Data *)elm->data;

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

} /* End of LINE2SetConnect */

/*
%F
*/
static void LINE2UpdateConnect( sElement *elm )
{
 int      i;
 sLINE2Data *data = 0L;

 data = (sLINE2Data *)elm->data;

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

} /* End of LINE2UpdateConnect */

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

 (*nnodes) = data->NumNodes;

} /* End of LINE2NumNodes */


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

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

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

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

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

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

} /* End of LINE2Gravity */


/*
%F
*/
static void LINE2AssVector( sElement *elm, double *GMatrix, double *matrix )
{
 int      i;
 int      conn[2];
 sLINE2Data *data;

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


/*
%F
*/
static void LINE2StressStrain( sElement *elm, double dt, double *U,
                               double *yield, sTensor *stre, sTensor *stra )
{
 sLINE2Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     cm[6][6];
 double     bm[3][6];
 double     ue[6];
 double     def[3], str[4];
 double     nu;
 double     pf = 1.0;
 int        i, j;

/* Initialize stress and strain variables.
 */
 stre = (sTensor *)memset( (void *)stre, 0, sizeof(sTensor) );
 if( stra != 0L )
  stra = (sTensor *)memset( (void *)stra, 0, sizeof(sTensor) );

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

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

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

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

/* Get poisson coeff.
 */
 MatNuParameter( mat, &nu );

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

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

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

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

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

/* Set stress values
 */
 stre->xx = str[0];
 stre->yy = str[1];
 stre->xy = str[2];
 stre->zz = str[3];

/* Compute off plane stress
 */
 ElementOffPlaneStress( nu, stre );

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

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

} /* End of LINE2StressStrain */


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

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

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

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


/*
%F This function computes the internal forces for the LINE2 element
*/
static void LINE2InterForce( sElement *elm, sTensor *stress, double *intforce )
{
 int      i, j;
 double   coeff = 0.0;
 double   bm[3][6], str[3];
 sLINE2Data *data = 0L;

/* Get element descriptor
 */ 
 data = (sLINE2Data *)elm->data;
 
/* Compute element rigidity coefficient
 */ 
 coeff = ElmRigidCoeff( elm );

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

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

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

/* Compute element internal forces
 */
 for( i = 0; i < (2*data->NumNodes); i++ ) {
  intforce[i] = 0.0;
  for( j = 0; j < 3; 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 LINE2InterForce */


/*
%F
*/
static void LINE2WriteStress( sElement *elm, FILE *out, double *U, double *V )
{
 sTensor  stress;
 sTensor  strain;
 sTensor  strd = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
 sPTensor pstr = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
 double   yield = 0.0;
 
/* 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
  */
  LINE2StressStrain( elm, 0.0, U, &yield, &stress, &strain );
 }

/* Compute desviatory stress
 */
 PrincipalTensor( &stress, &pstr );
 
 strd.xx = stress.xx - ((pstr.dir1 + pstr.dir3 + pstr.dir2)/3.0);
 strd.yy = stress.yy - ((pstr.dir1 + pstr.dir3 + pstr.dir2)/3.0);

/* Write element stress, strain and yield parameter
 */
 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 LINE2WriteStress */


/*
%F
*/
static void LINE2WriteNodalResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  sts, sta;
 sPTensor psts, psta;
 sTensor  strd;
 double   yield;
 int      i;

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

/* Computes the principal tensors
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 for( i = 0; i < 2; i++ ) {
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t%+10.5e\t", sts.xx,
                                                        sts.yy,
                                                        sts.zz,
                                                        sts.xy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psts.dir1,
                                               psts.dir2,
                                               psts.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\t", strd.xx,
                                      strd.yy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", sta.xx,
                                               sta.yy,
                                               sta.xy );
  fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psta.dir1,
                                               psta.dir2,
                                               psta.dir3 );
  fprintf( out, "%+10.5e\n", yield );
 }

} /* End of LINE2WriteNodalResult */


/*
%F
*/
static void LINE2WriteGaussResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  sts, sta;
 sPTensor psts, psta;
 sTensor  strd;
 double   yield;
 double   d1, d2, d3, smises;

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

/* Computes the principal stress
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );

/* Compute Von Mises stress.
 */
 d1 = (psts.dir1 - psts.dir2)*(psts.dir1 - psts.dir2);
 d2 = (psts.dir1 - psts.dir3)*(psts.dir1 - psts.dir3);
 d3 = (psts.dir2 - psts.dir3)*(psts.dir2 - psts.dir3);
 smises = sqrt( 0.5*(d1 + d2 + d3) );
 
/* Write results on output file
 */
 fprintf( out, "%d\t%d\n", elm->label, 1 );
 fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t%+10.5e\t", sts.xx,
                                                       sts.yy,
                                                       sts.zz,
                                                       sts.xy );
 fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psts.dir1,
                                              psts.dir2,
                                              psts.dir3 );
 fprintf( out, "%+10.5e\t%+10.5e\t", strd.xx,
                                     strd.yy );
 fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", sta.xx,
                                              sta.yy,
                                              sta.xy );
 fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\t", psta.dir1,
                                              psta.dir2,
                                              psta.dir3 );
 fprintf( out, "%+10.5e\t", (psta.dir1+psta.dir2+psta.dir3) );
 fprintf( out, "%+10.5e\t", yield );
 fprintf( out, "%+10.5e\n", smises );

} /* End of LINE2WriteGaussResult */


/*
%F
*/
static void LINE2WriteGaussVectorResult( sElement *elm, int version, FILE *out,
                                      FILE *tmp )
{
 sPTensor psts, psta;
 sTensor  sts, sta;
 sTensor  strd;
 double   yield;

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

/* Computes the principal stress
 */
 PrincipalTensor( &sts, &psts );
 PrincipalTensor( &sta, &psta );

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

 if( version > 1 ) {
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir1, psta.dir1 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir2, psta.dir2 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.dir3, psta.dir3 );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1x, psta.cos1x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2x, psta.cos2x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3x, psta.cos3x );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1y, psta.cos1y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2y, psta.cos2y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3y, psta.cos3y );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos1z, psta.cos1z );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos2z, psta.cos2z );
  fprintf( out, "%+10.5e\t%+10.5e\n", psts.cos3z, psta.cos3z );
 }
 else {
  if( psts.dir1 > 0.0 )
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", (psts.dir1*psts.cos1x),
                                                (psts.dir1*psts.cos1y), 0.0 );
  else
   fprintf( out, "%+10.5e\t%+10.5e\t%+10.5e\n", 0.0, 0.0, 0.0 );

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

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

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

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

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

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

} /* End of LINE2WriteGaussVectorResult */


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

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

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

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

/* Compute the displacement increment
 */
 for( i = 0; i < 2*data->NumNodes; i++ )
  ve[i] *= dtime;

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

/* Compute the element stress
 */
 sxx = stre->xx;
 syy = stre->yy;

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

} /* End of LINE2UpdateStress */


/*
%F
*/
static void LINE2PercForces( sElement *elm, double *pforce )
{
 sCoord     coord[2];
 sLINE2Data   *data = 0L;
 sMaterial *mat  = 0L;
 double     bm[3][6], phi[2];
 sTensor    grad;
 int        i, k;

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

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

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

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

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

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

/* Compute element internal forces
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  pforce[2*i]   = grad.xx;
  pforce[2*i+1] = grad.yy;
 }

} /* End of LINE2PercForces */


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

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

/* Compute element CG elevation
 */
 elev = 0.5 * (coord[0].y + coord[1].y);
 
/* Get element data
 */
 data = (sLINE2Data *)elm->data;

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

} /* End of LINE2SetPressure */

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

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

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

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

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

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

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

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

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

/* Compute viscos stress
 */
 for( i = 0; i < 3; i++ ) {
  strv[i] = 0.0;
  for( j = 0; j < 3; 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 < 3; j++ )
   vforce[i] += bm[j][i] * strv[j];
  vforce[i] *= coeff;
 } 

} /* End of LINE2ViscoForce */


static void LINE2Jacobian(sElement *elm, double jac[3][3], double ijac[3][3])
{
 double len;
 sCoord n1, n2, n3;
 sCoord ex = {1.0, 0.0, 0.0};
 sCoord ey = {0.0, 1.0, 0.0};
 sCoord ez = {0.0, 0.0, 1.0};
 sCoord coord[2];

 _LINE2GetCoord(elm, coord);

 n1.x = (coord[1].x - coord[0].x);
 n1.y = (coord[1].y - coord[0].y);
 n1.z = (coord[1].z - coord[0].z);

 len = _LINE2Lenght(&n1);

 /* normalize n1 */
 _LINE2Normalize( &n1, &n1 );

 n2 = _LINE2Cross(&ex, &n1);
 if( (_LINE2Lenght(&n2)/len) < 1.0e-5 ) {
  n2 = _LINE2Cross(&ey, &n1);
  if( (_LINE2Lenght(&n2)/len) < 1.0e-5 ) {
   n2 = _LINE2Cross(&ez, &n1);
  }
 }

 /* normalize n2 */
 _LINE2Normalize( &n2, &n2 );

 n3 = _LINE2Cross(&n1, &n2);

 jac[0][0] = n1.x;
 jac[0][1] = n1.y;
 jac[0][2] = n1.z;
 jac[1][0] = n2.x;
 jac[1][1] = n2.y;
 jac[1][2] = n2.z;
 jac[2][0] = n3.x;
 jac[2][1] = n3.y;
 jac[2][2] = n3.z;

 _LINE2Inverse(jac, ijac);

} /* End of LINE2Jacobian */

static void LINE2Volume(sElement *elm, double *v)
{
 sCoord coord[2];

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

/* Compute element area
 */
 (*v) = _LINE2Area(coord);

} /* End of LINE2Volume */

/*
 * Public functions:
 */

/*
%F This function initializes the sub-class "LINE2".
*/
void LINE2Init(void);
void LINE2Init(void)
{
/* Initialize LINE2 element sub-class
 */
 ElmClass[LINE2].new = LINE2New;
 ElmClass[LINE2].free = LINE2Free;
 ElmClass[LINE2].read = LINE2Read;
 ElmClass[LINE2].readinitstr = LINE2ReadInitStress;
 ElmClass[LINE2].mass = LINE2MassMatrix;
 ElmClass[LINE2].rigidcoeff = LINE2RigidCoeff;
 ElmClass[LINE2].load = NULL;
 ElmClass[LINE2].connect = LINE2Connect;
 ElmClass[LINE2].setconnect = LINE2SetConnect;
 ElmClass[LINE2].updateconnect = LINE2UpdateConnect;
 ElmClass[LINE2].numnodes = LINE2NumNodes;
 ElmClass[LINE2].gravity = LINE2Gravity;
 ElmClass[LINE2].assvector = LINE2AssVector;
 ElmClass[LINE2].strstrain = LINE2StressStrain;
 ElmClass[LINE2].timestep = LINE2TimeStep;
 ElmClass[LINE2].intforce = LINE2InterForce;
 ElmClass[LINE2].writestr = LINE2WriteStress;
 ElmClass[LINE2].writendlresult = LINE2WriteNodalResult;
 ElmClass[LINE2].writegauresult = LINE2WriteGaussResult;
 ElmClass[LINE2].writegauvecresult = LINE2WriteGaussVectorResult;
 ElmClass[LINE2].updatestress = LINE2UpdateStress;
 ElmClass[LINE2].percforce = LINE2PercForces;
 ElmClass[LINE2].setpressure = LINE2SetPressure;
 ElmClass[LINE2].setinitstress = LINE2SetInitStress;
 ElmClass[LINE2].viscoforce = LINE2ViscoForce;
 ElmClass[LINE2].jacobian = LINE2Jacobian;
 ElmClass[LINE2].volume = LINE2Volume;
 ElmClass[LINE2].KMatrix           = 0L;
 ElmClass[LINE2].GetDof            = 0L;

} /* End of LINE2Init */

/*
 * End of File.
 */

---