element.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

#include <string.h>
#include "element.h"
#include "global.h"
#include "elemhead.h"
#include "gtopology.h"
#include "intpoints.h"


element::element (void)
{
  tm = NULL;  idm = NULL; nm = 0;
  te = (elemtype) 0;
  prescdispl=0;  presctemp=0;  react=0;
  nb=0;  ipp=NULL;  intordsm=NULL;
  crst = (crsectype) 0;  idcs=0;
  ncomp=0;
  initdispl=NULL;
}

element::~element (void)
{
  long i;
  delete [] tm; delete [] idm;
  delete [] initdispl;
  for (i=0; i < nb; i++)
  {
    delete [] ipp[i];
    delete [] nip[i];
    delete [] intordsm[i];
  }
  delete [] nip;
  delete [] ipp;
  delete [] intordsm;
}

void element::read (XFILE *in,long eid)
{
  long nne,ndofe;
  gelemtype get;
  
  xfscanf (in,"%d",(int*)&te);
  
  switch (te){
  case bar2d:{                if (Bar2d==NULL)     Bar2d = new barel2d;           break;  }
  case bar3d:{                if (Bar3d==NULL)     Bar3d = new barel3d;           break;  }
  case barq2d:{               if (Barq2d==NULL)    Barq2d = new barelq2d;         break;  }
  case barq3d:{               if (Barq3d==NULL)    Barq3d = new barelq3d;         break;  }
  case beam2d:{               if (Beam2d==NULL)    Beam2d = new beamel2d;         break;  }
  case beam3d:{               if (Beam3d==NULL)    Beam3d = new beamel3d;         break;  }
  case beamg3d:{              if (Beam3dg==NULL)   Beam3dg = new beamgen3d;       break;  }
  case subsoilbeam:{          if (Sbeam==NULL)     Sbeam = new soilbeam;          break;  }
  case beam2dsp:{             if (Spbeam2d==NULL)  Spbeam2d = new beam2dspec;     break;  }
  case spring_1:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case spring_2:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case spring_3:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case spring_4:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case spring_5:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case spring_6:{             if (Spring==NULL)    Spring = new springel;         break;  }
  case planeelementlt:{       if (Pelt==NULL)      Pelt = new planeelemlt;        break;  }
  case planeelementqt:{       if (Peqt==NULL)      Peqt = new planeelemqt;        break;  }
  case planeelementrotlt:{    if (Perlt==NULL)     Perlt = new planeelemrotlt;    break;  }
  case planeelementlq:{       if (Pelq==NULL)      Pelq = new planeelemlq;        break;  }
  case planeelementqq:{       if (Peqq==NULL)      Peqq = new planeelemqq;        break;  }
  case planeelementrotlq:{    if (Perlq==NULL)     Perlq = new planeelemrotlq;    break;  }
  case planeelementsubqt:{    if (Pesqt==NULL)     Pesqt = new planeelemsubqt;    break;  }
  case planequadcontact:{     if (Pqcon==NULL)     Pqcon = new plquadcontact;     break;  }
  case cctel:{                if (Cct==NULL)       Cct = new cctelem;             break;  }
  case dktel:{                if (Dkt==NULL)       Dkt = new dktelem;             break;  }
  case dstel:{                if (Dst==NULL)       Dst = new dstelem;             break;  }
  case q4plateel:{            if (Q4pl==NULL)      Q4pl = new q4plate;            break;  }
  case subsoilplatetr:{       if (Spltr==NULL)     Spltr = new soilplatetr;       break;  }
  case subsoilplateq:{        if (Splq==NULL)      Splq = new soilplateq;         break;  }
  case axisymmlt:{            if (Asymlt==NULL)    Asymlt = new axisymlt;         break;  }
  case axisymmlq:{            if (Asymlq==NULL)    Asymlq = new axisymlq;         break;  }
  case axisymmqq:{            if (Asymqq==NULL)    Asymqq = new axisymqq;         break;  }
  case shelltrelem:{          if (Shtr==NULL)      Shtr = new shelltr;            break;  }
  case shellqelem:{           if (Shq==NULL)       Shq = new shellq;              break;  }
  case lineartet:{            if (Ltet==NULL)      Ltet = new lintet;             break;  }
  case linearhex:{            if (Lhex==NULL)      Lhex = new linhex;             break;  }
  case quadrhex:{             if (Qhex==NULL)      Qhex = new quadhex;            break;  }
  case linearwed:{            if (Lwed==NULL)      Lwed = new linwedge;           break;  }
  case quadrwed:{             if (Qwed==NULL)      Qwed = new quadwedge;          break;  }
  case particleelem:{
    if (Pelem==NULL){
      long nne,dim;
      xfscanf (in,"%ld %ld",&nne,&dim);
      Pelem = new elemparticle (nne,dim);
    }
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function");
    fprintf (stderr,"\n element::read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  switch (te){
  case planeelementlt:{     xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementqt:{     xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementrotlt:{  xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementlq:{     xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementqq:{     xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementrotlq:{  xfscanf (in,"%d",(int*)&ssst);  break;  }
  case planeelementsubqt:{  xfscanf (in,"%d",(int*)&ssst);  break;  }
  default:{
    //fprintf (stderr,"\n\n unknown element type is required in function");
    //fprintf (stderr,"\n element::read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

  switch (te){
  case bar2d:{                Bar2d->eleminit (eid);                     break;  }
  case bar3d:{                Bar3d->eleminit (eid);                     break;  }
  case barq2d:{               Barq2d->eleminit (eid);                    break;  }
  case barq3d:{               Barq3d->eleminit (eid);                    break;  }
  case beam2d:{               Beam2d->eleminit (eid);                    break;  }
  case beam3d:{               Beam3d->eleminit (eid);                    break;  }
  case beamg3d:{              Beam3dg->eleminit (eid);                   break;  }
  case subsoilbeam:{          Sbeam->eleminit (eid);                     break;  }
  case beam2dsp:{             Spbeam2d->eleminit (eid);                  break;  }
  case spring_1:{             Spring->eleminit (eid);                    break;  }
  case spring_2:{             Spring->eleminit (eid);                    break;  }
  case spring_3:{             Spring->eleminit (eid);                    break;  }
  case spring_4:{             Spring->eleminit (eid);                    break;  }
  case spring_5:{             Spring->eleminit (eid);                    break;  }
  case spring_6:{             Spring->eleminit (eid);                    break;  }
  case planeelementlt:{       Pelt->eleminit (eid);      get=lintriag;   break;  }
  case planeelementqt:{       Peqt->eleminit (eid);      get=quadtriag;  break;  }
  case planeelementrotlt:{    Perlt->eleminit (eid);     get=lintriag;   break;  }
  case planeelementlq:{       Pelq->eleminit (eid);      get=linquad;    break;  }
  case planeelementqq:{       Peqq->eleminit (eid);      get=quadquad;   break;  }
  case planeelementrotlq:{    Perlq->eleminit (eid);     get=linquad;    break;  }
  case planeelementsubqt:{    Pesqt->eleminit (eid);                     break;  }
  case planequadcontact:{     Pqcon->eleminit (eid);                     break;  }
  case cctel:{                Cct->eleminit (eid);       get=lintriag;   break;  }
  case dktel:{                Dkt->eleminit (eid);                       break;  }
  case dstel:{                Dst->eleminit (eid);                       break;  }
  case q4plateel:{            Q4pl->eleminit (eid);                      break;  }
  case subsoilplatetr:{       Spltr->eleminit (eid);                     break;  }
  case subsoilplateq:{        Splq->eleminit (eid);                      break;  }
  case axisymmlt:{            Asymlt->eleminit (eid);    get=lintriag;   break;  }
  case axisymmlq:{            Asymlq->eleminit (eid);    get=linquad;    break;  }
  case axisymmqq:{            Asymqq->eleminit (eid);    get=quadquad;   break;  }
  case shelltrelem:{          Shtr->eleminit (eid);                      break;  }
  case shellqelem:{           Shq->eleminit (eid);                       break;  }
  case lineartet:{            Ltet->eleminit (eid);      get=lintetra;   break;  }
  case linearhex:{            Lhex->eleminit (eid);      get=linhexa;    break;  }
  case quadrhex:{             Qhex->eleminit (eid);      get=quadhexa;   break;  }
  case linearwed:{            Lwed->eleminit (eid);                      break;  }
  case quadrwed:{             Qwed->eleminit (eid);                      break;  }
  case particleelem:{         Pelem->eleminit (eid);                     break;  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function");
    fprintf (stderr,"\n element::read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

  ndofe=Mt->give_ndofe (eid);
  nne = Mt->give_nne (eid);
  
  if (Mp->tprob != growing_mech_structure){
    Gtm->gelements[eid].read (in,nne,ndofe,get);
  }
  else{
    Gtm->gelements[eid].read_gf (in,nne,ndofe);
  }
  
  xfscanf (in,"%d",(int*)&crst);
  if (crst!=0){
    xfscanf (in,"%ld",&idcs);
    idcs--;
  }
  
  readmat (in);
  
}

void element::readmat (XFILE *in)
{
  long i, idx, idxt;

  xfscanf (in, "%ld", &nm);
  tm  = new mattype[nm];
  idm = new long[nm];
  memset (tm,  0, sizeof(*tm)*nm);
  memset (idm, 0, sizeof(*idm)*nm);
  idx = 0;
  idxt = nm-1;
  for (i = 0; i < nm; i++)
  {
    xfscanf (in, "%d %ld", (int*)(tm+idx), idm+idx);
    idm[idx]--;
    if (tm[idx] == therisodilat || tm[idx] == therisodilattime){
      //Mp->temperature=1;
      tm[idxt] = tm[idx];
      idm[idxt] = idm[idx];
      idxt--;
      idx--;
    }
    
    idx++;
  }
}

/**
   function allocates array for initial displacements
   initial displacements are used in problems with growing number of elements
   
   @praram ndofe - number of DOFs on element
   
   JK, 3.3.2006
*/
void element::alloc_initdispl (long ndofe)
{
  initdispl = new double [ndofe];
}

/**
   function defines initial displacements
   initial displacements are used in problems with growing number of elements
   
   @param r - array containing initial displacements
   @param ndofe - number of DOFs on element
   
   JK, 3.3.2006
*/
void element::initdisplacement (double *r,long ndofe)
{
  long i;
  
  for (i=0;i<ndofe;i++){
    initdispl[i]=r[i];
  }
  
}

/**
   function subtracts initial displacements from total displacements
   function is used in problems with changing number of elements
   
   @param r - total displacement
   @param ndofe - number of DOFs on element
   
   5.3.2006, JK
*/
void element::subtrinitdispl (double *r,long ndofe)
{
  long i;
  
  for (i=0;i<ndofe;i++){
    r[i]-=initdispl[i];
  }
}

/**
   function allocates array used for problems with
   changing nodes, elements and DOFs
   
   @param eid - element id
   
   7.11.2006, JK
*/
void element::alloc_growstr (long eid)
{
  long i,ndofe;
  
  ndofe=Mt->give_ndofe (eid);

  initdispl = new double [ndofe];

  for (i=0;i<ndofe;i++){
    initdispl[i]=0.0;
  }

}