plquadcontact.cpp

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

#include "plquadcontact.h"
#include "global.h"
#include "element.h"
#include "node.h"
#include "globmat.h"

plquadcontact::plquadcontact (void)
{
  long i,j;

  nne=4;  ndofe=8;
  nb=1;  tncomp=2;

  ncomp = new long [nb];
  ncomp[0]=2;

  nip = new long* [nb];
  //intordsm = new long* [nb];
  for (i=0;i<nb;i++){
    nip[i] = new long [nb];
    //intordsm[i] = new long [nb];
  }
  
  nip[0][0]=2;
  
  tnip=0;
  for (i=0;i<nb;i++){
    for (j=0;j<nb;j++){
      tnip+=nip[i][j];
    }
  }
}

plquadcontact::~plquadcontact (void)
{
  long i;
  
  for (i=0;i<nb;i++){
    delete [] nip[i];
    //delete [] intordsm[i];
  }
  delete [] nip;
  //delete [] intordsm;

  delete [] ncomp;
  
}

/**
   function initializes basic data on elements
   
   @param eid - element id
   
   JK, 11.6.2006
*/
void plquadcontact::eleminit (long eid)
{

  long ii,jj;

  Mt->elements[eid].nb=nb;
  //Mt->elements[eid].intordsm = new long* [nb];
  Mt->elements[eid].nip = new long* [nb];

  for (ii=0;ii<nb;ii++){
    //Mt->elements[eid].intordsm[ii] = new long [nb];
    Mt->elements[eid].nip[ii] = new long [nb];
    for (jj=0;jj<nb;jj++){
      //Mt->elements[eid].intordsm[ii][jj]=intordsm[ii][jj];
      Mt->elements[eid].nip[ii][jj]=nip[ii][jj];
    }
  }

}

/**
   function assembles transformation %matrix x_g = T x_l
   
   @param nodes - array containing node numbers
   @param tmat - transformation %matrix
   
   JK, 11.6.2006
*/
void plquadcontact::transf_matrix (ivector &nod,matrix &tmat)
{
  long i,n,m;

  fillm (0.0,tmat);

  n=nod.n;
  m=tmat.m;
  for (i=0;i<m;i++){
    tmat[i][i]=1.0;
  }
  
  for (i=0;i<n;i++){
    if (Mt->nodes[nod[i]].transf>0){
      tmat[i*2][i*2]   = Mt->nodes[nod[i]].e1[0];    tmat[i*2][i*2+1]   = Mt->nodes[nod[i]].e2[0];
      tmat[i*2+1][i*2] = Mt->nodes[nod[i]].e1[1];    tmat[i*2+1][i*2+1] = Mt->nodes[nod[i]].e2[1];
    }
  }
}
















/**
   function computes stiffness %matrix of plane rectangular
   finite element for contact problems
   
   @param eid - number of element
   @param ri,ci - row and column indices
   @param sm - stiffness %matrix
   @param x,y - vectors of nodal coordinates
   
   JK, 11.6.2006
*/
void plquadcontact::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
{
  long ipp;
  matrix d(2,2);
  
  fillm (0.0,sm);
  
  //  number of appropriate integration point
  ipp=Mt->elements[eid].ipp[ri][ci];

  //  matrix of stiffness of the material
  Mm->matstiff (d,ipp);
  
  sm[0][0]=d[0][0];
  sm[0][2]=0.0-d[0][0];

  sm[1][1]=d[1][1];
  sm[1][3]=0.0-d[1][1];

  sm[2][0]=0.0-d[0][0];
  sm[2][2]=d[0][0];
  
  sm[3][1]=0.0-d[1][1];
  sm[3][3]=d[1][1];
  


  sm[4][4]=d[0][0];
  sm[4][6]=0.0-d[0][0];

  sm[5][5]=d[1][1];
  sm[5][7]=0.0-d[1][1];

  sm[6][4]=0.0-d[0][0];
  sm[6][6]=d[0][0];

  sm[7][5]=0.0-d[1][1];
  sm[7][7]=d[1][1];

}


/**
   function computes stiffness %matrix of quadrilateral element
   for contact problems
   
   @param lcid - load case id
   @param eid - element id
   @param sm - stiffness %matrix
   
   JK, 11.6.2006
*/
void plquadcontact::res_stiffness_matrix (long lcid,long eid,matrix &sm)
{
  long transf;
  ivector nodes(nne);
  vector x(nne),y(nne);
  Mt->give_node_coord2d (x,y,eid);
  Mt->give_elemnodes (eid,nodes);
  
  stiffness_matrix (eid,0,0,sm,x,y);
  
  //  transformation of stiffness matrix
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
}



















/**
   function computes strains at integration points
   
   @param lcid - load case id
   @param eid - element id
   
   JK, 11.6.2006
*/
void plquadcontact::res_mainip_strains (long lcid,long eid)
{
  vector aux,x(nne),y(nne),r(ndofe);
  ivector cn(ndofe),nodes(nne);
  matrix tmat;
  
  //  coordinates of nodes of element
  Mt->give_node_coord2d (x,y,eid);
  //  node numbers of element
  Mt->give_elemnodes (eid,nodes);
  //  code numbers (numbers of DOFs)
  Mt->give_code_numbers (eid,cn.a);
  //  nodal displacements
  eldispl (lcid,eid,r.a,cn.a,ndofe);
  
  //  transformation of displacement vector
  long transf = Mt->locsystems (nodes);
  if (transf>0){
    allocv (ndofe,aux);
    allocm (ndofe,ndofe,tmat);
    transf_matrix (nodes,tmat);
    //locglobtransf (aux,r,tmat);
    lgvectortransf (aux,r,tmat);
    copyv (aux,r);
    destrv (aux);
    destrm (tmat);
  }
  
  //
  mainip_strains (lcid,eid,0,0,x,y,r);
  
}

/**
   function computes strains at integration points of element
   function is used in geometrically linear problems
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   @param ii - number of block
   @param x,y - arrays with node coordinates
   @param r - %vector of nodal displacements
   
   JK, 11.6.2006
*/
void plquadcontact::mainip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
{
  long ipp;
  vector displdiff(2);
  
  ipp=Mt->elements[eid].ipp[ri][ci];
  
  displdiff[0]=r[0]-r[2];
  displdiff[1]=r[1]-r[3];
  
  Mm->storestrain (lcid,ipp,displdiff);
  ipp++;
  
  displdiff[0]=r[6]-r[4];
  displdiff[1]=r[7]-r[5];
  
  Mm->storestrain (lcid,ipp,displdiff);
  
}












/**
   function computes internal forces
   
   @param lcid - number of load case
   @param eid - element id
   @param ri,ci - row and column indices
   @param ifor - vector of internal forces
   @param x,y - vectors of nodal coordinates
   
   JK, 11.6.2006
*/
void plquadcontact::internal_forces (long lcid,long eid,long ri,long ci,vector &ifor,vector &x,vector &y)
{
  long ipp;
  vector nodfor (2);

  //  number of apropriate integration point
  ipp=Mt->elements[eid].ipp[ri][ci];
  
  if (Mp->strcomp==1)
    Mm->computenlstresses (ipp);
  
  Mm->givestress (lcid,ipp,nodfor);
  
  ifor[0]=nodfor[0];
  ifor[1]=nodfor[1];
  ifor[2]=0.0-nodfor[0];
  ifor[3]=0.0-nodfor[1];

  if (Mp->strcomp==1)
    Mm->computenlstresses (ipp);
  
  Mm->givestress (lcid,ipp,nodfor);
  
  ifor[4]=0.0-nodfor[0];
  ifor[5]=0.0-nodfor[1];
  ifor[6]=nodfor[0];
  ifor[7]=nodfor[1];

}


/**
   function computes internal forces
   
   @param lcid - load case id
   @param eid - element id
   @param ifor - %vector of internal forces
   
   JK, 11.6.2006
*/
void plquadcontact::res_internal_forces (long lcid,long eid,vector &ifor)
{
  vector x(nne),y(nne);

  Mt->give_node_coord2d (x,y,eid);

  internal_forces (lcid,eid,0,0,ifor,x,y);
}