mechtop.cpp

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

  case beamg3d:{             ndofe=Beam3dg->ndofe;            break; }
  case subsoilbeam:{         ndofe=Sbeam->ndofe;              break; }
  case beam2dsp:{            ndofe=Spbeam2d->ndofe;           break; }
  case spring_1:{            ndofe=Spring->give_ndofe(eid);   break; }
  case spring_2:{            ndofe=Spring->give_ndofe(eid);   break; }
  case spring_3:{            ndofe=Spring->give_ndofe(eid);   break; }
  case spring_4:{            ndofe=Spring->give_ndofe(eid);   break; }
  case spring_5:{            ndofe=Spring->give_ndofe(eid);   break; }
  case spring_6:{            ndofe=Spring->give_ndofe(eid);   break; }
  case planeelementlt:{      ndofe=Pelt->ndofe;               break; }
  case planeelementqt:{      ndofe=Peqt->ndofe;               break; }
  case planeelementrotlt:{   ndofe=Perlt->ndofe;              break; }
  case planeelementlq:{      ndofe=Pelq->ndofe;               break; }
  case planeelementqq:{      ndofe=Peqq->ndofe;               break; }
  case planeelementrotlq:{   ndofe=Perlq->ndofe;              break; }
  case planeelementsubqt:{   ndofe=Pesqt->ndofe;              break; }
  case planequadcontact:{    ndofe=Pqcon->ndofe;              break; }
  case cctel:{               ndofe=Cct->ndofe;                break; }
  case dktel:{               ndofe=Dkt->ndofe;                break; }
  case dstel:{               ndofe=Dst->ndofe;                break; }
  case q4plateel:{           ndofe=Q4pl->ndofe;               break; }
  case subsoilplatetr:{      ndofe=Spltr->ndofe;              break; }
  case subsoilplateq:{       ndofe=Splq->ndofe;               break; }
  case axisymmlt:{           ndofe=Asymlt->ndofe;             break; }
  case axisymmlq:{           ndofe=Asymlq->ndofe;             break; }
  case axisymmqq:{           ndofe=Asymqq->ndofe;             break; }
  case shelltrelem:{         ndofe=Shtr->ndofe;               break; }
  case shellqelem:{          ndofe=Shq->ndofe;                break; }
  case lineartet:{           ndofe=Ltet->ndofe;               break; }
  case linearhex:{           ndofe=Lhex->ndofe;               break; }
  case quadrhex:{            ndofe=Qhex->ndofe;               break; }
  case linearwed:{           ndofe=Lwed->ndofe;               break; }
  case quadrwed:{            ndofe=Qwed->ndofe;               break; }
  case particleelem:{        ndofe=Pelem->ndofe;              break; }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function give_ndofe (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return ndofe;
}

/**
   function returns number of DOFs of node
   
   @param nid - node id
*/
long mechtop::give_ndofn (long nid)
{
  return Gtm->give_ndofn (nid);
}

/**
   function returns appropriate DOF of node
   
   @param nid - node id
   @param n - number of required DOF
*/
long mechtop::give_dof (long nid,long n)
{
  return Gtm->give_dof (nid,n);
}

/**
   function returns number of nodes of the element
   
   @param eid - element id
*/
long mechtop::give_nne (long eid)
{
  long nne;
  elemtype te;
  
  te = give_elem_type (eid);
  nne = give_nne_inner (te);
  return nne;
}

/**
   function returns number nodes on elements
   
   @param te - element type
   
   27.2.2002
*/
long mechtop::give_nne_inner (elemtype te)
{
  long nne=0;

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

/**
   function returns number of integration points
   
   @param eid - element id
   @param ri,ci - row and column indices
   
   30.12.2001
*/
long mechtop::give_nip (long eid,long ri,long ci)
{
  long nip=0;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){

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

/**
   function returns number of integration points
   
   @param eid - element id
   
   21.6.2004, JK
*/
long mechtop::give_tnip (long eid)
{
  long tnip = -1;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){
  case bar2d:   {           tnip=Bar2d->tnip;      break;  }
  case bar3d:   {           tnip=Bar3d->tnip;      break;  }
  case barq2d:  {           tnip=Barq2d->tnip;     break;  }
  case barq3d:  {           tnip=Barq3d->tnip;     break;  }
  case beam2d:{             tnip=Beam2d->tnip;     break;  }
  case beam3d:{             tnip=Beam3d->tnip;     break;  }
  case subsoilbeam:{        tnip=Sbeam->tnip;      break;  }
  case beamg3d:{            tnip=Beam3dg->tnip;    break;  }
  case beam2dsp:{           tnip=Spbeam2d->tnip;   break;  }
  case spring_1:{           tnip=Spring->tnip;     break;  }
  case spring_2:{           tnip=Spring->tnip;     break;  }
  case spring_3:{           tnip=Spring->tnip;     break;  }
  case spring_4:{           tnip=Spring->tnip;     break;  }
  case spring_5:{           tnip=Spring->tnip;     break;  }
  case spring_6:{           tnip=Spring->tnip;     break;  }
  case planeelementlt:{     tnip=Pelt->tnip;       break;  }
  case planeelementqt:{     tnip=Peqt->tnip;       break;  }
  case planeelementrotlt:{  tnip=Perlt->tnip;      break;  }
  case planeelementlq:{     tnip=Pelq->tnip;       break;  }
  case planeelementqq:{     tnip=Peqq->tnip;       break;  }
  case planeelementrotlq:{  tnip=Perlq->tnip;      break;  }
  case planeelementsubqt:{  tnip=Pesqt->tnip;      break;  }
  case planequadcontact:{   tnip=Pqcon->tnip;      break;  }
  case cctel:{              tnip=Cct->tnip;        break;  }
  case dktel:{              tnip=Dkt->tnip;        break;  }
  case dstel:{              tnip=Dst->tnip;        break;  }
  case q4plateel:{          tnip=Q4pl->tnip;       break;  }
  case subsoilplatetr:{     tnip=Spltr->tnip;      break;  }
  case subsoilplateq:{      tnip=Splq->tnip;       break;  }
  case axisymmlt:{          tnip=Asymlt->tnip;     break;  }
  case axisymmlq:{          tnip=Asymlq->tnip;     break;  }
  case axisymmqq:{          tnip=Asymqq->tnip;     break;  }
  case shelltrelem:{        tnip=Shtr->tnip;       break;  }
  case shellqelem:{         tnip=Shq->tnip;        break;  }
  case lineartet:{          tnip=Ltet->tnip;       break;  }
  case linearhex:{          tnip=Lhex->tnip;       break;  }
  case quadrhex:{           tnip=Qhex->tnip;       break;  }
  case linearwed:{          tnip=Lwed->tnip;       break;  }
  case quadrwed:{           tnip=Qwed->tnip;       break;  }
  case particleelem:{       tnip=1;                break;  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function give_tnip (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return tnip;
}

/**
   function returns total number of integration points on element
   
   @param eid - element id
   
   30.12.2001
*/
long mechtop::give_totnip (long eid)
{
  if (eid < Mt->ne-1)
    return elements[eid+1].ipp[0][0] - elements[eid].ipp[0][0];

  return Mm->tnip-elements[eid].ipp[0][0];
}

/**
   function returns number of strain/stress components
   
   @param eid - element id
   
   30.12.2001
*/
long mechtop::give_ncomp (long eid)
{
  long ncomp=0;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){
  case bar2d:{
    //ncomp=4;
    ncomp=Bar2d->tncomp;
    break;  }
  case bar3d:{              ncomp=Bar3d->tncomp;      break;  }
  case barq2d:{             ncomp=Barq2d->tncomp;     break;  }
  case barq3d:{             ncomp=Barq3d->tncomp;     break;  }
  case beam2d:{             ncomp=Beam2d->tncomp;     break;  }
  case beam3d:{             ncomp=Beam3d->tncomp;     break;  }
  case beamg3d:{            ncomp=Beam3dg->tncomp;    break;  }
  case beam2dsp:{           ncomp=Spbeam2d->tncomp;   break;  }
  case subsoilbeam:{        ncomp=Sbeam->tncomp;      break;  }
  case spring_1:{           ncomp=Spring->tncomp;     break;  }
  case spring_2:{           ncomp=Spring->tncomp;     break;  }
  case spring_3:{           ncomp=Spring->tncomp;     break;  }
  case spring_4:{           ncomp=Spring->tncomp;     break;  }
  case spring_5:{           ncomp=Spring->tncomp;     break;  }
  case spring_6:{           ncomp=Spring->tncomp;     break;  }
  case planeelementlt:{     ncomp=Pelt->tncomp;       break;  }
  case planeelementqt:{     ncomp=Peqt->tncomp;       break;  }
  case planeelementrotlt:{  ncomp=Perlt->tncomp;      break;  }
  case planeelementlq:{     ncomp=Pelq->tncomp;       break;  }
  case planeelementqq:{     ncomp=Peqq->tncomp;       break;  }
  case planeelementrotlq:{  ncomp=Perlq->tncomp;      break;  }
  case planeelementsubqt:{  ncomp=Pesqt->tncomp;      break;  }
  case planequadcontact:{   ncomp=Pqcon->tncomp;      break;  }
  case cctel:{              ncomp=Cct->tncomp;        break;  }
  case dktel:{              ncomp=Dkt->tncomp;        break;  }
  case dstel:{              ncomp=Dst->tncomp;        break;  }
  case q4plateel:{          ncomp=Q4pl->tncomp;       break;  }
  case subsoilplatetr:{     ncomp=Spltr->tncomp;      break;  }
  case subsoilplateq:{      ncomp=Splq->tncomp;       break;  }
  case axisymmlt:{          ncomp=Asymlt->tncomp;     break;  }
  case axisymmlq:{          ncomp=Asymlq->tncomp;     break;  }
  case axisymmqq:{          ncomp=Asymqq->tncomp;     break;  }
  case shelltrelem:{        ncomp=Shtr->tncomp;       break;  }
  case shellqelem:{         ncomp=Shq->tncomp;        break;  }
  case lineartet:{          ncomp=Ltet->tncomp;       break;  }
  case linearhex:{          ncomp=Lhex->tncomp;       break;  }
  case quadrhex:{           ncomp=Qhex->tncomp;       break;  }
  case linearwed:{          ncomp=Lwed->tncomp;       break;  }
  case quadrwed:{           ncomp=Qwed->tncomp;       break;  }
  case particleelem:{       ncomp=Pelem->ndofe;       break;  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function give_ncomp (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return ncomp;
}

/**
   function returns number of edges of the element
   
   @param eid - element id
   
   30.12.2001
*/
long mechtop::give_ned (long eid)
{
  long ned=0;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){
  case bar3d:{              ned=0;            break;  }
  case beam3d:{             ned=0;            break;  }
  case planeelementlt:{     ned=Pelt->ned;    break;  }
  case planeelementqt:{     ned=Peqt->ned;    break;  }
  case planeelementrotlt:{  ned=Perlt->ned;   break;  }
  case planeelementlq:{     ned=Pelq->ned;    break;  }
  case planeelementqq:{     ned=Peqq->ned;    break;  }
  case planeelementrotlq:{  ned=Perlq->ned;   break;  }
  case planeelementsubqt:{  ned=Pesqt->ned;   break;  }
  case cctel:{              ned=Cct->ned;     break;  }
  case dktel:{              ned=Dkt->ned;     break;  }
  case dstel:{              ned=Dst->ned;     break;  }
  case q4plateel:{          ned=Q4pl->ned;    break;  }
  case subsoilplatetr:{     ned=Spltr->ned;   break;  }
  case subsoilplateq:{      ned=Splq->ned;    break;  }
  case axisymmlt:{          ned=Asymlt->ned;  break;  }
  case axisymmlq:{          ned=Asymlq->ned;  break;  }
  case axisymmqq:{          ned=Asymqq->ned;  break;  }
  case shelltrelem:{        ned=Shtr->ned;    break;  }
  case shellqelem:{         ned=Shq->ned;     break;  }
  case lineartet:{          ned=Ltet->ned;    break;  }
  case linearhex:{          ned=Lhex->ned;    break;  }
  case quadrhex:{           ned=Qhex->ned;    break;  }
  case linearwed:{          ned=Lwed->ned;    break;  }
  case quadrwed:{           ned=Qwed->ned;    break;  }
  case particleelem:{                         break;  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function give_ned (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  return ned;
}

/**
   function returns number of nodes on one element edge
   
   @param eid - number of element
   
   30.12.2001
*/
long mechtop::give_nned (long eid)
{
  long nned=0;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){
  case bar3d:{              nned=0;             break;  }
  case beam3d:{             nned=0;             break;  }