elm.c

Go to the documentation of this file.

/*
%M This modules contains the element super class methods and definitions
%a Joao Luiz Elias Campos.
%d September 1st, 1998.
%r $Id: elm.c,v 1.2 2004/06/22 17:53:10 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.
*/

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

#include "load.h"
#include "elm.h"

/*
%V Element class, element list and number of mesh elements
*/
sElmClass   ElmClass[NumElmTypes];
sElement  **ElmList     = 0L;
int         NumElements = 0;

/*
%V Integration order and thickness
*/
sOrder *IntOrder     = 0L;
double *Thickness    = 0L;
int     NumIntOrder  = 0;
int     NumThickness = 0;

/*
%V Rezone variables
*/
int NumRezones = 0;
sRezone *Rezone = 0L;


/*
** ------------------------------------------------------------------------
** Local variables and symbols:
*/
#ifndef PI
#define PI 3.141592654
#endif
/*
** ------------------------------------------------------------------------
** Local functions:
*/
void T3Init(void);
void Q4Init(void);
void BRICK8Init(void);
void TETR4Init(void);
void INFINITEInit(void);
void INTERFACEInit(void);
void LINE2Init(void);
void DKTInit(void);

static int fn_tp(double *x, double *y);

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

/*
%F This function initializes the globally shared "ElmClass" variable. This
   variable holds the pointers to the necessary functions which implement 
   the element class methods.
*/
void ElementInit(void)
{
 int i;

/* Define the subclasses initialization methods
 */
 ElmClass[T3].init = T3Init;
 ElmClass[Q4].init = Q4Init;
 ElmClass[BRICK8].init = BRICK8Init;
 ElmClass[TETR4].init = TETR4Init;
 ElmClass[INTERFACE].init = INTERFACEInit;
 ElmClass[INFINITE].init = INFINITEInit;
 ElmClass[LINE2].init = LINE2Init;
 ElmClass[DKT].init = DKTInit; 

 for( i = 0; i < NumElmTypes; i++ )
  ElmInit(i);

/* Allocate memory for element list
 */
 if( ElmList != 0L )
  free(ElmList);

 ElmList = (sElement **)calloc( NumElements, sizeof(sElement *) );

} /* End of ElementInit */


/*
%F This function release memory used by element super-class.
*/
void ElementFree(void)
{
 int i;

/* Releases sub-classes
 */
 for( i = 0; i < NumElements; i++ )
  ElmFree(ElmList[i]);

/* Release element list and reset it
 */
 free(ElmList);

 ElmList = 0L;

} /* End of ElementFree */


/* 
%F This function builds the element adjacence for the line load
*/
void ElementBuildAdjacence( void )
{
 int       i, j, k;
 int       noi, noj, nnode;
 int       conn[4];
 sElement *elm = 0L;

 for( i = 0; i < NumLineForce; i++ )
 {
  LineForce[i].adj = 0;
  for( j = 0; j < NumElements; j++ )
  {
   elm   = ElmList[j];
   ElmNumNodes(elm,&nnode);
   ElmConnect(elm,conn);
   if( (LineForce[i].elem - 1) != j )
   {
    for( k = 0; k < nnode; k++ )
    {
     noi = conn[k];
     noj = conn[(k+1)%nnode];
     if( (noj == LineForce[i].noi) && (noi == LineForce[i].noj))
      LineForce[i].adj = j + 1;
    }
   }
  }
 }

} /* End of ElementBuildAdjacence */


/*
%F This function computes the principal components for a given tensor.
*/
void PrincipalTensor( sTensor *tensor, sPTensor *ptensor )
{
 double sxx = tensor->xx;
 double syy = tensor->yy;
 double sxy = tensor->xy;
 double szz = tensor->zz;
 double sxz = tensor->xz;
 double syz = tensor->yz;
 double tmed, ratio, ra;
 double ang1, ang2;
 double len1, len2;
 double tp[3], de1, de2;
 double I1, I2, I3;
 double q, r, phi;
 double ZERO = 1.0E-10;

 if( NDof == 2) {
   if( sxy >= ZERO ){
    ra = sxx * sxx - 2.0 * sxx * syy + syy * syy + 4 * sxy * sxy;
    if(ra > 0.0)
     ratio = sqrt(ra);
    else ratio = 0.0;
    tmed  = 0.5 * (sxx + syy);
    tp[0] = tmed - 0.5 * ratio;
    tp[1] = tmed + 0.5 * ratio;
    de1 = syy - sxx - ratio;
    de2 = syy - sxx + ratio;
  
    /* artificio utilizado para
     evitar erro numerico
     na determinacao dos cossenos
     diretores
    */
     if(de1 == 0.0)
      de1 = ZERO;
     if(de2 == 0.0)
      de2 = ZERO;
    /*
     determina as tensoes principais
     maxima e minima
    */
     if(tp[0] <= tp[1]){
      ptensor->dir1 = tp[0];
      ptensor->dir3 = tp[1];
      ang1 = -2.0 * sxy / de1;
      ang2 = -2.0 * sxy / de2;
     }
   else{
    ptensor->dir1 = tp[1];
    ptensor->dir3 = tp[0];
    ang1 = -2.0 * sxy / de2;
    ang2 = -2.0 * sxy / de1;
   }
  /*
   calcula os cossenos diretores
  */
   len1 = sqrt(ang1*ang1+1.0);
   ptensor->cos1x = ang1 / len1;
   ptensor->cos1y = 1.0 / len1;
 
   len2 = sqrt(ang2*ang2+1.0);
   ptensor->cos3x = ang2 / len2;
   ptensor->cos3y = 1.0 / len2;
  }
  else{
  /* nesse caso as tensoes principais
   sao iguais as tensoes  sxx e syy
   bastando apenas verificar os valores
   maximo e minimo
  */ 
   if(syy <= sxx){
    ptensor->dir1 = syy;
    ptensor->dir3 = sxx;
    ptensor->cos1x = 0.0;
    ptensor->cos1y = 1.0;
    ptensor->cos3x = 1.0;
    ptensor->cos3y = 0.0;
   }
   else{
    ptensor->dir1 = sxx;
    ptensor->dir3 = syy;
    ptensor->cos1x = 1.0;
    ptensor->cos1y = 0.0;
    ptensor->cos3x = 0.0;
    ptensor->cos3y = 1.0;
   }
  }
  ptensor->dir2 = szz;
  ptensor->cos1z = ptensor->cos3z = 
  ptensor->cos2x = ptensor->cos2y = 
  ptensor->cos2z = 0.0;
 }
 else{
  /*
   calcula os invariantes de tensoes
  */
  I1 = sxx + syy + szz; 
  I2 = sxx * syy - sxy * sxy + sxx * szz -
       sxz * sxz + syy * szz - syz * syz;
  I3 = sxx * syy * szz - sxx * syz * syz -
       sxy * sxy * szz + 2.0 * sxy * sxz *
       syz - sxz * sxz * syy;   
  /* 
   calcula as raizes do polinomio caracteristico
   determinando as tensoes principais
  */
  q = (3.0 * I2 - (I1 * I1))/9.0;
  if(q <= 0.0){
   r = (9.0 * I1 * I2 - 27.0 * I3 - 2.0 * I1 * I1 * I1)/54.0;
   
   if((q*q*q+r*r)>=0.0){
    ptensor->dir1 = 0.0;
    ptensor->dir2 = 0.0;
    ptensor->dir3 = 0.0;
   }
   else{
    ang1 = 120.0 * PI/180.0;
    ang2 = 240.0 * PI/180.0;
    phi = acos(r/sqrt(-(q * q * q)));
   
    tp[0] = 2 * sqrt(-q)*cos(1.0/3.0 * phi) + 1.0/3.0 * I1;
    tp[1] = 2 * sqrt(-q)*cos(1.0/3.0 * phi + ang1) + 1.0/3.0 * I1;
    tp[2] = 2 * sqrt(-q)*cos(1.0/3.0 * phi + ang2) + 1.0/3.0 * I1;
   
    qsort( tp, 3, sizeof(double), (void *) fn_tp ); 
   
    ptensor->dir1 = tp[0];
    ptensor->dir2 = tp[1];
    ptensor->dir3 = tp[2];
   }
  }
  else{
   ptensor->dir1 = ptensor->dir2 = ptensor->dir3 = 0.0;
  }
  /*
   fazer o calculo das dir. principais
  */
  ptensor->cos1z = ptensor->cos2z = 
  ptensor->cos3z = ptensor->cos1x = 
  ptensor->cos2x = ptensor->cos3x =
  ptensor->cos1y = ptensor->cos2y = 
  ptensor->cos3y = 0.0;
     
 }

} /* End of PrincipalTensor */

static int fn_tp(double *x, double *y)
{
 if (*x > *y)
  return 1;
 else if (*x == *y)
  return 0;
 return -1;
}

/*
%F Compute the off plane stress for the plane strain analysis
*/
void ElementOffPlaneStress( double nu, sTensor *str )
{
 sPTensor pstr;

/* Compute the principal stress
 */
 PrincipalTensor( str, &pstr );

/* Compute the off plane stress
 */
 str->zz = nu * (pstr.dir1 + pstr.dir3);

} /* End of ElementOffPlaneStress */


/*
%F This function initializes a given tensor.
*/
void ElementInitTensor( sTensor *tensor )
{
/* Initialize tensor
 */
 tensor->xx = tensor->xy = tensor->xz =
              tensor->yy = tensor->yz =
                           tensor->zz = 0.0;

} /* End of ElementInitTensor */


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

---