dkt.c

Go to the documentation of this file.

/*
%M This modules contains the DKT element sub-class methods and 
   definitions
%a André Luis Muller.
%d 08-2007.
OBS: a implementação desse elemento difere do elemento (discrete Kirchhoff triangle)
     original, afim de atender algumas especificações do programa Recon.

/*
** ------------------------------------------------------------------------
** 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"
#include "load.h"
#include "alg.h"


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

/*
%T DKT element data definition
*/
typedef struct _dktelm {
 int      matid;             /* Material identification          */
 int      tckid;             /* Element thickness identification */
 int      NumNodes;          /* Number of element nodes          */
 int      NumTensComp;       /* Number of tension components     */
 int      inc[3];            /* Element incidence                */
 double   istr[4];           /* Element initial stress           */
 double   effdef;            /* Effective deformation            */
 sTensor *iPress;            /* Element internal pressure        */
 double prof;                /* Element profile                  */
} sDKTData;

/*
%V Gauss point coordinates 
*/
static double _GPoint  = 0.3333333;

/*
** ------------------------------------------------------------------------
** Local functions:
*/
static void   _DKTGetCoord       ( sElement *, sCoord [] );
static double _DKTArea           ( sCoord [] );
static double _DKTMinHeight      ( sCoord [] );
static void   _DKTGetDisplacement( sElement *, double *, double [] );
static void   _DKTBMatrix        ( sElement *, double [6][9] );
static void   _DKTDerivShp       ( double, double, double, sDerivRST [] );

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

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

/* Get element coordinates
 */
 for( i = 0; i < 3; i++ ) {
  coord[i].x = elm->nodes[data->inc[i]-1].coord.x;
  coord[i].y = elm->nodes[data->inc[i]-1].coord.y;
  coord[i].z = elm->nodes[data->inc[i]-1].coord.z;
 }

} /* End of _DKTGetCoord */


/*
%F This function computes the element area
%i Element coordinates.
%o Element area.
*/
static double _DKTArea( sCoord coord[3] )
{
 double area = 0.0;
 double u[3],v[3],w[3];
/* Computes element area
 */
 u[0] = coord[1].x - coord[0].x;
 u[1] = coord[1].y - coord[0].y; 
 u[2] = coord[1].z - coord[0].z;
 v[0] = coord[2].x - coord[0].x; 
 v[1] = coord[2].y - coord[0].y; 
 v[2] = coord[2].z - coord[0].z;

 w[0] = u[1] * v[2] - u[2] * v[1];
 w[1] = u[2] * v[0] - u[0] * v[2];
 w[2] = u[0] * v[1] - u[1] * v[0];
 area = sqrt(w[0]*w[0] + w[1]*w[1] + w[2]*w[2]) / 2.0;

 return( area );

} /* End of _DKTArea */


/*
%F This function computes the minimum height for the element.
%i Element coordinate.
%o Minimum height
*/
static double _DKTMinHeight( sCoord coord[3] )
{
 int    i;
 double area;
 double bcorr;
 double bmax;
 double hmin;

/* Initializa variables
 */
 bmax = bcorr = hmin = area = 0.0;

/* Loop over the element edges to compute the maximum edge
 */
 for( i = 0; i < 3; i++ ) {
  bcorr = sqrt( ( coord[(i+1)%3].x - coord[i].x ) * ( coord[(i+1)%3].x - coord[i].x ) +
                ( coord[(i+1)%3].y - coord[i].y ) * ( coord[(i+1)%3].y - coord[i].y ) +
                ( coord[(i+1)%3].z - coord[i].z ) * ( coord[(i+1)%3].z - coord[i].z ));

  if( bcorr > bmax )
   bmax = bcorr;
 }

/* Computes the element area
 */
 area = _DKTArea( coord );

/* Computes the minimum height
 */
 hmin = (2.0 * area) / bmax;

 return( hmin );

} /* End of _DKTMinHight */


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

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

/* Fill element displacement vector
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  ue[3*i]   = U[3*(data->inc[i]-1)];
  ue[3*i+1] = U[3*(data->inc[i]-1)+1];
  ue[3*i+2] = U[3*(data->inc[i]-1)+2];
 }

} /* End of _DKTGetDisplacement */

/*
%F This function computes the shape functions derivatives for the 
   _DKT element.
*/
static void _DKTDerivShp( double r, double s, double t, sDerivRST dshp[3] )
{
/* Compute shape functions derivatives
 */
 dshp[0].r =  1.0;
 dshp[1].r =  0.0;
 dshp[2].r = -1.0;
 

 dshp[0].s =  0.0;
 dshp[1].s =  1.0;
 dshp[2].s = -1.0;
 

 dshp[0].t =  0.0;
 dshp[1].t =  0.0;
 dshp[2].t =  0.0;
 

} /* End of _DKTDerivShp */


/*
%F This function computes the stress x strain matrix for the DKT element.
*/
static void _DKTBMatrix( sElement *elm, double bm[6][9] )
{
 int i;
 double detjac;
 double jac[3][3], ijac[3][3];
 sCoord coord[3];
 sDerivRST dshp[3];
 
 for( i = 0; i < 9; i++ )
  bm[0][i] = bm[1][i] = bm[2][i] = 
  bm[3][i] = bm[4][i] = bm[5][i] = 0.0;

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

/* Compute the element jacobian matrix
 */
 jac[0][0] = coord[0].x - coord[2].x;
 jac[0][1] = coord[0].y - coord[2].y;
 jac[0][2] = coord[0].z - coord[2].z;
 jac[1][0] = coord[1].x - coord[2].x;
 jac[1][1] = coord[1].y - coord[2].y;
 jac[1][2] = coord[1].z - coord[2].z;
 jac[2][0] = jac[2][1] = 0.0;
 jac[2][2] = 1.0;

/* Compute the matrix determinant
 */ 
 detjac = _DKTArea( coord );

/* Compute the inverse
 */
 ijac[0][0] =  (jac[1][1] * jac[2][2] - jac[1][2] * jac[2][1]) / detjac;
 ijac[0][1] = -(jac[0][1] * jac[2][2] - jac[0][2] * jac[2][1]) / detjac;
 ijac[0][2] =  (jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]) / detjac;
 ijac[1][0] = -(jac[1][0] * jac[2][2] - jac[1][2] * jac[2][0]) / detjac;
 ijac[1][1] =  (jac[0][0] * jac[2][2] - jac[0][2] * jac[2][0]) / detjac;
 ijac[1][2] = -(jac[0][0] * jac[1][2] - jac[0][2] * jac[1][0]) / detjac;
 ijac[2][0] =  (jac[1][0] * jac[2][1] - jac[1][1] * jac[2][0]) / detjac;
 ijac[2][1] = -(jac[0][0] * jac[2][1] - jac[0][1] * jac[2][0]) / detjac;
 ijac[2][2] =  (jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) / detjac;

/* Compute the B matrix coefficientes
 */
 _DKTDerivShp( _GPoint, _GPoint, _GPoint, dshp );
 
 for( i = 0; i < 3; i++ ) {
  bm[0][3*i]   = bm[3][3*i+1] = bm[4][3*i+2] = (ijac[0][0] * dshp[i].r) +
                                               (ijac[0][1] * dshp[i].s) + 
                                               (ijac[0][2] * dshp[i].t);
  bm[1][3*i+1] = bm[3][3*i]   = bm[5][3*i+2] = (ijac[1][0] * dshp[i].r) +
                                               (ijac[1][1] * dshp[i].s) + 
                                               (ijac[1][2] * dshp[i].t);
  bm[2][3*i+2] = bm[4][3*i]   = bm[5][3*i+1] = (ijac[2][0] * dshp[i].r) +
                                               (ijac[2][1] * dshp[i].s) + 
                                               (ijac[2][2] * dshp[i].t);
 }
 
} /* End of _DKTBMatrix */


/*
** ------------------------------------------------------------------------
** Subclass methods:
*/
static void   DKTNew( int, int, int, int, sElement **, sElement **, sNode * );
static void   DKTFree                  ( sElement * );
static void   DKTRead                  ( sElement * );
static void   DKTMassMatrix            ( sElement *, double * );

static double DKTRigidCoeff            ( sElement * );
static void   DKTConnect               ( sElement *, int * );
static void   DKTNumNodes              ( sElement *, int * );
static void   DKTGravity               ( sElement *, double *, double *, 
                                                    double * );
static void   DKTAssVector             ( sElement *, double *, double * );
static void   DKTStressStrain          ( sElement *, double, double *, double *,
                                                    sTensor *, sTensor * );
static void   DKTTimeStep              ( sElement *, double * );
static void   DKTInterForce            ( sElement *, sTensor *, double * );
static void   DKTWriteStress           ( sElement *, FILE *, double *,
                                                            double * );
static void   DKTWriteNodalResult      ( sElement *, FILE *, FILE * );
static void   DKTWriteGaussResult      ( sElement *, FILE *, FILE * );
static void   DKTWriteGaussVectorResult( sElement *, int, FILE *, FILE * );


static void   DKTSetInitStress         ( sElement *, sTensor * );

static void   DKTVolume                ( sElement *, double * );

static void DKTKMatrix( sElement *, double [24][24] );                                                                         
static void DKTGetDof( sElement *, int [24], int * );
/*
%F This method allocates memory for a DKT element and fills its data with 
   the specific data.
%i Element label (number)
%o Element descriptor.
*/
static void DKTNew( int label, int matid, int intord, int tckid, 
                   sElement **elm, sElement **elist, sNode *nodes )
{
 sDKTData *data = 0L;

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

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

/* Fill DKT element data
 */
 data->matid       = matid;
 data->tckid       = tckid;
 data->NumNodes    = 3;
 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    = DKT;
 (*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 DKTNew */


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

/* Get DKT element data
 */
 data = (sDKTData *)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 DKTFree */


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

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

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

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

} /* End of DKTRead */


/*
%F This method reads the DKT initial stress
*/
static int DKTReadInitStress( sElement *elm )
{
 int      numpg = 1;
 double   sxx, syy, sxy ,szz, sxz, syz;
 sDKTData *data = 0L;

/* Get element data
 */
 data = (sDKTData *)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 %lf %lf %lf", &sxx, &syy, &szz, &sxy, &sxz, &syz);
 data->istr[0] = sxx;
 data->istr[1] = syy;
 data->istr[2] = sxy;

 return 1;
} /* End of DKTReadInitStress */


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

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

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

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

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

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


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

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

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

 return( area );

} /* End of DKTRigidCoeff */


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

 data = (sDKTData *)elm->data;

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

} /* End of DKTConnect */


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

 (*nnodes) = data->NumNodes;

} /* End of DKTNumNodes */


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

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

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

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

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

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

} /* End of DKTGravity */


/*
%F
*/
static void DKTAssVector( sElement *elm, double *GMatrix, double *matrix )
{
 int      i;
 int      conn[3];
 sDKTData *data;

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

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

/* Assemble matrix
 */
 for( i = 0; i < data->NumNodes; i++ ) {
  GMatrix[3*(conn[i]-1)]   += matrix[3*i];
  GMatrix[3*(conn[i]-1)+1] += matrix[3*i+1];
  GMatrix[3*(conn[i]-1)+2] += matrix[3*i+2];
 }
 
} /* End of DKTAssVector */


/*
%F
*/
static void DKTStressStrain( sElement *elm, double dt, double *U, double *yield, 
                                           sTensor *stre, sTensor *stra )
{
 sDKTData   *data = 0L;
 sMaterial *mat  = 0L;
 double     cm[6][6];
 double     bm[6][9];
 double     ue[9];
 double     def[6], str[6];
 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 = (sDKTData *)elm->data;

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

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

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

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

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

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

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

/* Compute the element stress
 */
 for( i = 0; i < 6; i++ ) {
  str[i] = 0.0;
  for( j = 0; j < 6; j++ )
   str[i] += cm[i][j] * def[j];
  str[i] += data->istr[i];
 }
 /* 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->zz = str[2];
 stre->xy = str[3];
 stre->yz = str[4];
 stre->xz = str[5];

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

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

} /* End of DKTtressStrain */



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

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

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

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


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

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

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

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

/* Check for the initial stress
 */
 if( stress == 0L ) {
  for( i = 0; i < (3*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->zz;
 str[3] = stress->xy;
 str[4] = stress->yz;
 str[5] = stress->xz;

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


} /* End of DKTInterForce */


/*
%F
*/
static void DKTWriteStress( 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 };
 double   yield = 0.0;
 double   iso;
 
/* 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
  */
  DKTStressStrain( elm, 0.0, U, &yield, &stress, &strain );
 }

/* Compute desviatory stress
 */
 iso = (stress.xx + stress.yy + stress.zz)/3.0;
 strd.xx = stress.xx - iso;
 strd.yy = stress.yy - iso;

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


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

/* 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 );
 
 /* Computes the isotropic stress
 */
 iso = (sts.xx + sts.yy + sts.zz)/3.0;

/* Write results on output file
 */
 fprintf( out, "%d\n", elm->label );
 for( i = 0; i < 3; 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%+10.5e\t", strd.xx,
                                               strd.yy,
                                               iso );
  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 DKTWriteNodalResult */


/*
%F
*/
static void DKTWriteGaussResult( sElement *elm, FILE *out, FILE *tmp )
{
 sTensor  sts, sta;
 sPTensor psts, psta, epsts, dpsts;
 sTensor  strd, estr;
 double   yield;
 double   iso, piso;
 double   poropressure = 0.0;
 sMaterial *mat  = 0L;
 sDKTData   *data = 0L;
 int      j;

/* 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 );
 
/* Computes the isotropic stress
 */
 iso = (sts.xx + sts.yy + sts.zz)/3.0;
 
/* Compute desviatory principal stress
 */
 piso = (psts.dir1 + psts.dir2 + psts.dir3)/3.0;
 
 dpsts.dir1 = psts.dir1 - piso;
 dpsts.dir3 = psts.dir3 - piso;
 
 /* Get element descriptor
 */
 data = (sDKTData *)elm->data;

/* Get Material object
 */
 mat = MatList[data->matid-1];
 
 estr.xx = sts.xx;
 estr.yy = sts.yy;
 estr.xy = sts.xy;
 
/* Computes the principal efetive stress
 */
 PrincipalTensor( &estr, &epsts );  
 
/* 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], "Sxx" ) == 0))
    fprintf( out, "%+10.5e\t", sts.xx );
   else if((strcmp( Config.gaussresults[j], "Syy" ) == 0))
    fprintf( out, "%+10.5e\t", sts.yy );
   else if((strcmp( Config.gaussresults[j], "Szz" ) == 0))
    fprintf( out, "%+10.5e\t", sts.zz );
   else if((strcmp( Config.gaussresults[j], "Sxy" ) == 0))
    fprintf( out, "%+10.5e\t", sts.xy );
   else if((strcmp( Config.gaussresults[j], "Sxz" ) == 0))
    fprintf( out, "%+10.5e\t", sts.xz );
   else if((strcmp( Config.gaussresults[j], "Syz" ) == 0))
    fprintf( out, "%+10.5e\t", sts.yz );
   else if((strcmp( Config.gaussresults[j], "S1" ) == 0))
    fprintf( out, "%+10.5e\t", psts.dir1 );
   else if((strcmp( Config.gaussresults[j], "S2" ) == 0))
    fprintf( out, "%+10.5e\t", psts.dir2 );
   else if((strcmp( Config.gaussresults[j], "S3" ) == 0))
    fprintf( out, "%+10.5e\t", psts.dir3 );
   else if((strcmp( Config.gaussresults[j], "Sdx" ) == 0))
    fprintf( out, "%+10.5e\t", strd.xx );
   else if((strcmp( Config.gaussresults[j], "Sdy" ) == 0))
    fprintf( out, "%+10.5e\t", strd.yy );
   else if((strcmp( Config.gaussresults[j], "Sdz" ) == 0))
    fprintf( out, "%+10.5e\t", strd.zz );
   else if((strcmp( Config.gaussresults[j], "Exx" ) == 0))
    fprintf( out, "%+10.5e\t", sta.xx );
   else if((strcmp( Config.gaussresults[j], "Eyy" ) == 0))
    fprintf( out, "%+10.5e\t", sta.yy );
   else if((strcmp( Config.gaussresults[j], "Ezz" ) == 0))
    fprintf( out, "%+10.5e\t", sta.zz );
   else if((strcmp( Config.gaussresults[j], "Exy" ) == 0))
    fprintf( out, "%+10.5e\t", sta.xy );
   else if((strcmp( Config.gaussresults[j], "Exz" ) == 0))
    fprintf( out, "%+10.5e\t", sta.xz );
   else if((strcmp( Config.gaussresults[j], "Eyz" ) == 0))
    fprintf( out, "%+10.5e\t", sta.yz );
   else if((strcmp( Config.gaussresults[j], "E1" ) == 0))
    fprintf( out, "%+10.5e\t", psta.dir1 );
   else if((strcmp( Config.gaussresults[j], "E2" ) == 0))
    fprintf( out, "%+10.5e\t", psta.dir2 );
   else if((strcmp( Config.gaussresults[j], "E3" ) == 0))
    fprintf( out, "%+10.5e\t", psta.dir3 );
   else if((strcmp( Config.gaussresults[j], "Ev" ) == 0))
    fprintf( out, "%+10.5e\t", (psta.dir1+psta.dir2+psta.dir3) );
 }                               
 
} /* End of DKTWriteGaussResult */


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

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

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

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

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

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

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

} /* End of DKTVolume */

static void DKTKMatrix( sElement *elm, double k[24][24] )
{         
 int    i, j, l;           
 sCoord  coord[3];                  
 double detjac;            
 double C[6][6];          
 double B[6][9];          
 double BTC[9][6];        
 double ke[9][9]; 
 sDKTData  *data = 0L;
 sMaterial  *mat  = 0L;  

 /* Inicia matriz k = 0
 */
 for( i = 0; i < 24; i++ ) {
  for( j = 0; j < 24; j++ ) {
   k[i][j] = 0.0;
  }
 }
 
 /* Pega descrição do elemento
 */
 data = (sDKTData *)elm->data;
 
/* Pega material
 */
 mat = MatList[data->matid-1];
 
 /* Pega as coordenadas nodais
 */
 _DKTGetCoord(elm, coord);
 
 /* Pega matriz constitutiva [C]*/
 //MatPSMatrix( mat, C );
 MatConstitutiveMatrix( mat, C ); 

 /* Calcula matriz B*/
 _DKTBMatrix(elm, B );
 
 /* Calcula a área do elemento */
 detjac = _DKTArea( coord);
 
 /*Calcula [B] transposto * [C] */
 for( i = 0; i < 9; i++ ) {
  for( j = 0; j < 6; j++)  {
   BTC[i][j] =0.0;
   for(l=0; l<6; l++){
    BTC[i][j] += B[l][i]*C[l][j];
   }
  }
 }
 /* Calcula [B] transposto * [C] * [B] * |jac|*/
 for( i = 0; i < 9; i++ ) {
  for( j = 0; j < 9; j++)  {
   ke[i][j] =0.0;
   for(l=0; l<6; l++){
    ke[i][j] += BTC[i][l]*B[l][j]*detjac;
   }
  }
 }
 /* Calcula matriz de rigidez do elemento*/
 for( i = 0; i < 9; i++ ) {
  for( j = 0; j < 9; j++)  {
   k[i][j] += ke[i][j];
  }
 }
 
} /* DKTKMatrix */

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

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

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

} /* End of DKTGetDof */

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

/* Get element descriptor
 */
 data = (sDKTData *)elm->data;
 
 for( i = 0; i < data->NumNodes; i++ ) {
  inc[i] = data->inc[i]-1;
 }
 
 *index = 3;
} /* End of DKTGetInc */

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

/*
%F This function initializes the sub-class "DKT".
*/
void DKTInit( void );
void DKTInit( void )
{
/* Initialize DKT element sub-class
 */
 ElmClass[DKT].new               = DKTNew;
 ElmClass[DKT].free              = DKTFree;
 ElmClass[DKT].read              = DKTRead;
 ElmClass[DKT].readinitstr       = DKTReadInitStress;
 ElmClass[DKT].readprofile       = 0L;
 ElmClass[DKT].mass              = DKTMassMatrix;
 ElmClass[DKT].rigidcoeff        = DKTRigidCoeff;
 ElmClass[DKT].load              = 0L;
 ElmClass[DKT].connect           = DKTConnect;
 ElmClass[DKT].numnodes          = DKTNumNodes;
 ElmClass[DKT].gravity           = DKTGravity;
 ElmClass[DKT].assvector         = DKTAssVector;
 ElmClass[DKT].strstrain         = DKTStressStrain;
 ElmClass[DKT].timestep          = DKTTimeStep;
 ElmClass[DKT].intforce          = DKTInterForce;
 ElmClass[DKT].writestr          = DKTWriteStress;
 ElmClass[DKT].writendlresult    = DKTWriteNodalResult;
 ElmClass[DKT].writegauresult    = DKTWriteGaussResult;
 ElmClass[DKT].writegauvecresult = DKTWriteGaussVectorResult;
 ElmClass[DKT].updatestress      = 0L;
 ElmClass[DKT].percforce         = 0L;
 ElmClass[DKT].setpressure       = 0L;
 ElmClass[DKT].setinitstress     = DKTSetInitStress;
 ElmClass[DKT].viscoforce        = 0L;
 ElmClass[DKT].jacobian          = 0L;
 ElmClass[DKT].volume            = DKTVolume;
 ElmClass[DKT].KMatrix           = DKTKMatrix;
 ElmClass[DKT].GetDof            = DKTGetDof;
 ElmClass[DKT].GetInc            = DKTGetInc;
 
} /* End of DKTInit */


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

---