plelemqq.cpp

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

#include "plelemqq.h"
#include "plelemlq.h"
#include "global.h"
#include "globmat.h"
#include "genfile.h"
#include "adaptivity.h"
#include "gadaptivity.h"
#include "node.h"
#include "intpoints.h"
#include "element.h"
#include "loadcase.h"
#include <stdlib.h>
#include <math.h>



planeelemqq::planeelemqq (void)
{
  long i,j;
  
  //  number nodes on element
  nne=8;
  //  number of DOFs on element
  ndofe=16;
  //  number of strain/stress components
  tncomp=3;
  //  number of functions approximated
  napfun=2;
  //  order of numerical integration of mass matrix
  intordmm=4;
  //  number of edges on element
  ned=4;
  //  number of nodes on one edge
  nned=3;
  //  order of numerical integration on element edges (boundaries)
  intordb=3;
  
  //  number of blocks (parts of geometric matrix)
  nb=1;
  
  //  number of strain/stress components
  ncomp = new long [nb];
  ncomp[0]=3;
  
  //  cumulative number of components approximated
  cncomp = new long [nb];
  cncomp[0]=0;
  
  //  number of integration points
  //  order of numerical integration of stiffness matrix
  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]=9;
  
  //  total number of integration points
  tnip=0;
  for (i=0;i<nb;i++){
    for (j=0;j<nb;j++){
      tnip+=nip[i][j];
    }
  }
  
  intordsm[0][0]=3;
}

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

void planeelemqq::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];
    }
  }
}

/**
   procedure approximates function defined by nodal values
   
   @param xi,eta - natural coordinates on element
   @param nodval - nodal values
   
   JK
*/
double planeelemqq::approx (double xi,double eta,vector &nodval)
{
  double f;
  vector bf(nne);
  
  bf_quad_4_2d (bf.a,xi,eta);
  
  scprd (bf,nodval,f);
  
  return f;
}

/**
   function returns %matrix of approximation functions
   
   @param n - %matrix of approximation functions
   @param xi,eta - natural coordinates
   
   JK, 25.8.2001
*/
void planeelemqq::bf_matrix (matrix &n,double xi,double eta)
{
  long i,j,k;
  vector bf(nne);
  
  fillm (0.0,n);
  
  bf_quad_4_2d (bf.a,xi,eta);
  
  j=0;  k=1;
  for (i=0;i<nne;i++){
    n[0][j]=bf[i];
    n[1][k]=bf[i];
    j+=2;  k+=2;
  }
}


/**
   function assembles geometric %matrix (strain-displacement %matrix)
   
   @param gm - geometric %matrix
   @param x,y - array containing node coordinates
   @param xi,eta - natural coordinates
   @param jac - Jacobian
   
   JK, 9.7.2001
*/
void planeelemqq::geom_matrix (matrix &gm,vector &x,vector &y,double xi,double eta,double &jac)
{
  long i,i1,i2;
  vector dx(nne),dy(nne);
  
  dx_bf_quad_4_2d (dx.a,xi,eta);
  dy_bf_quad_4_2d (dy.a,xi,eta);
  
  derivatives_2d (dx,dy,jac,x,y,xi,eta);
  
  fillm (0.0,gm);
  
  i1=0;  i2=1;
  for (i=0;i<nne;i++){
    gm[0][i1]=dx[i];
    gm[1][i2]=dy[i];
    gm[2][i1]=dy[i];
    gm[2][i2]=dx[i];
    i1+=2;  i2+=2;
  }
}

/**
   function assembles transformation %matrix from local nodal coordinate
   system to the global coordinate system x_g = T x_l
   
   @param nodes - element nodes
   @param tmat - transformation %matrix
   
   JK,
*/
void planeelemqq::transf_matrix (ivector &nodes,matrix &tmat)
{
  long i,n,m;
  
  fillm (0.0,tmat);
  
  n=nodes.n;
  m=tmat.m;
  for (i=0;i<m;i++){
    tmat[i][i]=1.0;
  }
  
  for (i=0;i<n;i++){
    if (Mt->nodes[nodes[i]].transf>0){
      tmat[i*2][i*2]   = Mt->nodes[nodes[i]].e1[0];    tmat[i*2][i*2+1]   = Mt->nodes[nodes[i]].e2[0];
      tmat[i*2+1][i*2] = Mt->nodes[nodes[i]].e1[1];    tmat[i*2+1][i*2+1] = Mt->nodes[nodes[i]].e2[1];
    }
  }
}

/**
   function computes stiffness %matrix of plane stress rectangular
   finite element with biquadratic approximation functions
   
   this function is used in plane stress/strain elements (function is called
   by function res_stiffness_matrix) and shell elements

   @param eid - number of element
   @param ri,ci - row and column indices
   @param sm - stiffness %matrix
   @param x,y - node coordinates
   
   JK, 25.8.2001
*/
void planeelemqq::stiffness_matrix (long eid,long ri,long ci,matrix &sm,vector &x,vector &y)
{
  long i,j,ii,jj,ipp;
  double xi,eta,jac,thick;
  ivector nodes(nne);
  vector w,gp,t(nne);
  matrix d(tncomp,tncomp),gm(tncomp,ndofe);
  
  //  element nodes
  Mt->give_elemnodes (eid,nodes);
  //  thickness of the element
  Mc->give_thickness (eid,nodes,t);
  
  fillm (0.0,sm);
  
  for (ii=0;ii<nb;ii++){
    for (jj=0;jj<nb;jj++){
      if (intordsm[ii][jj]==0)  continue;
      allocv (intordsm[ii][jj],w);
      allocv (intordsm[ii][jj],gp);
      gauss_points (gp.a,w.a,intordsm[ii][jj]);
      
      ipp=Mt->elements[eid].ipp[ri+ii][ci+jj];
      
      for (i=0;i<intordsm[ii][jj];i++){
	xi=gp[i];
	for (j=0;j<intordsm[ii][jj];j++){
	  eta=gp[j];
	  
	  //  geometric matrix
	  geom_matrix (gm,x,y,xi,eta,jac);
	  
	  //  matrix of material stiffness
	  Mm->matstiff (d,ipp);
	  
	  thick = approx (xi,eta,t);
	  
	  jac*=thick*w[i]*w[j];
	  
	  //  contribution to the stiffness matrix of the element
	  bdbj (sm.a,gm.a,d.a,jac,gm.m,gm.n);
	  
	  ipp++;
	}
      }
      destrv (gp);  destrv (w);
    }
  }
  
}

/**
   function assembles stiffness %matrix of plane stress rectangular
   finite element with biquadratic approximation functions
   
   @param eid - element id
   @param sm - stiffness %matrix
   
   JK
*/
void planeelemqq::res_stiffness_matrix (long eid,matrix &sm)
{
  long transf;
  ivector nodes(nne);
  vector x(nne),y(nne);
  matrix tmat (ndofe,ndofe);
  
  Mt->give_node_coord2d (x,y,eid);
  
  stiffness_matrix (eid,0,0,sm,x,y);

  //  transformation of stiffness matrix
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tmat);
    glmatrixtransf (sm,tmat);
  }
}

/**
   function computes mass %matrix of the plane stress rectangular
   finite element with biquadratic approximation functions
   
   this function is used in plane stress/strain elements (function is called
   by function res_mass_matrix) and shell elements

   @param eid - number of element
   @param mm - mass %matrix
   @param x,y - node coordinates
   
   JK, 25.8.2001
*/
void planeelemqq::mass_matrix (long eid,matrix &mm,vector &x,vector &y)
{
  long i,j;
  double jac,xi,eta,w1,w2,thick,rho;
  ivector nodes(nne);
  vector w(intordmm),gp(intordmm),t(nne),dens(nne);
  matrix n(napfun,ndofe);
  
  //  element nodes
  Mt->give_elemnodes (eid,nodes);
  //  thickness of the element
  Mc->give_thickness (eid,nodes,t);
  //  density of material (defined at nodes or on element)
  Mc->give_density (eid,nodes,dens);

  gauss_points (gp.a,w.a,intordmm);
  
  fillm (0.0,mm);

  for (i=0;i<intordmm;i++){
    xi=gp[i];  w1=w[i];
    for (j=0;j<intordmm;j++){
      eta=gp[j];  w2=w[i];
      jac_2d (jac,x,y,xi,eta);
      
      //  matrix of approximation functions
      bf_matrix (n,xi,eta);
      
      //  thickness at integration point
      thick = approx (xi,eta,t);
      //  density at integration point
      rho = approx (xi,eta,dens);
      jac*=w1*w2*thick*rho;
      
      nnj (mm.a,n.a,jac,n.m,n.n);
    }
  }
  
}

/**
   function assembles mass %matrix of plane stress rectangular
   finite element with biquadratic approximation functions
   
   @param eid - element id
   @param mm - mass %matrix
   
   JK
*/
void planeelemqq::res_mass_matrix (long eid,matrix &mm)
{
  long transf;
  ivector nodes(nne);
  vector x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);
  mass_matrix (eid,mm,x,y);
  
  //  transformation of mass matrix
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    glmatrixtransf (mm,tmat);
  }
}

/**
   function computes load %matrix of the plane stress rectangular
   finite element with biquadratic approximation functions
   load vector is obtained after premultiplying load %matrix
   by nodal load values
   
   this function is used in plane stress/strain elements (function is called
   by function res_load_matrix) and shell elements

   @param eid - number of element
   @param lm - load %matrix
   @param x,y - node coordinates

   JK, 25.8.2001
*/
void planeelemqq::load_matrix (long eid,matrix &lm,vector &x,vector &y)
{
  long i,j;
  double jac,xi,eta,w1,w2,thick;
  ivector nodes(nne);
  vector w(intordmm),gp(intordmm),t(nne);
  matrix n(napfun,ndofe);
  
  //  element nodes
  Mt->give_elemnodes (eid,nodes);
  //  thickness of the element
  Mc->give_thickness (eid,nodes,t);

  gauss_points (gp.a,w.a,intordmm);
  
  fillm (0.0,lm);

  for (i=0;i<intordmm;i++){
    xi=gp[i];  w1=w[i];
    for (j=0;j<intordmm;j++){
      eta=gp[j];  w2=w[j];

      jac_2d (jac,x,y,xi,eta);

      //  matrix of approximation functions
      bf_matrix (n,xi,eta);
      
      //  thickness at integration point
      thick = approx (xi,eta,t);
      jac*=w1*w2*thick;
      
      nnj (lm.a,n.a,jac,n.m,n.n);
    }
  }
  
}

/**
   function assembles load %matrix of plane stress rectangular
   finite element with biquadratic approximation functions
   
   @param eid - element id
   @param lm - load %matrix
   
   JK
*/
void planeelemqq::res_load_matrix (long eid,matrix &lm)
{
  long transf;
  ivector nodes(nne);
  vector x(nne),y(nne);
  
  Mt->give_node_coord2d (x,y,eid);

  load_matrix (eid,lm,x,y);

  //  transformation of load matrix
  //  (in the case of nodal coordinate systems)
  Mt->give_elemnodes (eid,nodes);
  transf = Mt->locsystems (nodes);
  if (transf>0){
    matrix tmat (ndofe,ndofe);
    transf_matrix (nodes,tmat);
    glmatrixtransf (lm,tmat);
  }

}

/**
   function computes strains at integration points
   
   @param lcid - load case id
   @param eid - element id
   
   JK, modified 23.11.2006
*/
void planeelemqq::res_ip_strains (long lcid,long eid)
{
  long transf;
  vector x(nne),y(nne),r(ndofe),aux(ndofe);
  ivector cn(ndofe),nodes(nne);
  matrix tmat(ndofe,ndofe);
  
  Mt->give_node_coord2d (x,y,eid);
  Mt->give_elemnodes (eid,nodes);
  Mt->give_code_numbers (eid,cn.a);
  eldispl (lcid,eid,r.a,cn.a,ndofe);
  
  //  transformation of displacement vector
  //  (in the case of nodal coordinate systems)
  transf = Mt->locsystems (nodes);
  if (transf>0){
    transf_matrix (nodes,tmat);
    lgvectortransf (aux,r,tmat);
    copyv (aux,r);
  }
  
  ip_strains (lcid,eid,0,0,x,y,r);
}

/**
   function computes strains at integration points of element
   
   this function is used in plane stress/strain elements (function is called
   by function res_ip_strains) and shell elements

   @param lcid - load case id
   @param eid - element id
   @param ri - row index
   @param ci - column index
   @param x,y - node coordinates
   @param r - nodal displacements
   
   10.5.2002, JK, modified 23.11.2006
*/
void planeelemqq::ip_strains (long lcid,long eid,long ri,long ci,vector &x,vector &y,vector &r)
{
  long i,j,ipp;
  double xi,eta,jac;
  vector gp,w,eps;
  matrix gm;

  allocv (intordsm[0][0],gp);
  allocv (intordsm[0][0],w);
  allocv (ncomp[0],eps);
  allocm (ncomp[0],ndofe,gm);
  
  gauss_points (gp.a,w.a,intordsm[0][0]);
  
  ipp=Mt->elements[eid].ipp[ri][ci];
  for (i=0;i<intordsm[0][0];i++){
    xi=gp[i];
    for (j=0;j<intordsm[0][0];j++){
      eta=gp[j];
      
      //  geometric matrix (strain-displacement matrix)
      geom_matrix (gm,x,y,xi,eta,jac);
      //  strain computation
      mxv (gm,r,eps);
      
      Mm->storestrain (lcid,ipp,eps);
      ipp++;
    }
  }
  
  destrm (gm);  destrv (eps);  destrv (w);  destrv (gp);
}

/**
   function assembles strains at nodes of element
   strains are obtained from the nearest integration points

   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices (default value is 0, nonzero values are used in shell elements)

   10.5.2002
*/
void planeelemqq::nod_strains_ip (long lcid,long eid,long ri,long ci)
{
  long i,j;
  ivector ipnum(nne),nod(nne);
  vector eps(tncomp);
  
  //  numbers of integration points closest to nodes
  nodipnum (eid,ri,ci,ipnum);
  
  //  node numbers of the element
  Mt->give_elemnodes (eid,nod);
  
  for (i=0;i<nne;i++){
    //  strains at the closest integration point
    Mm->givestrain (lcid,ipnum[i],eps);
    
    //  storage of strains to the node
    j=nod[i];
    Mt->nodes[j].storestrain (lcid,0,eps);
  }
  
}

/**
   function computes strains at required position on elements
   function possibly transforms strains to required coordinate systems
   
   @param lcid - load case id
   @param eid - element id
   @param ri,ci - row and column indices
   
   JK, 23.11.2006
*/
void planeelemqq::strains (long lcid,long eid,long ri,long ci)
{
  vector coord,eps;
  
  switch (Mm->stra.tape[eid]){
  case nowhere:{
    break;
  }
  case intpts:{
    //res_ip_strains (lcid,eid,ri,ci);
    break;
  }
  case enodes:{
    nod_strains_ip (lcid,eid,ri,ci);