mechtop.cpp

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

/**
   function returns node numbers on required element edge
   
   @param eid - element id
   @param edid - edge id
   @param nodes - array containing edge nodes
   
   5.3.2002
*/
void mechtop::give_edge_nodes (long eid,long edid,long *nodes)
{
  long nne;
  elemtype te;
  ivector nod;
  
  te = give_elem_type (eid);
  nne = give_nne (eid);
  allocv (nne,nod);
  give_elemnodes (eid,nod);
  
  switch (te){
  case planeelementlt:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[0]; }
    break;
  }
  case planeelementqt:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1];  nodes[2]=nod[3]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2];  nodes[2]=nod[4]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[0];  nodes[2]=nod[5]; }
    break;
  }
  case planeelementrotlt:{
    if (edid==0){  nodes[0]=nod[1];  nodes[1]=nod[2]; }
    if (edid==1){  nodes[0]=nod[2];  nodes[1]=nod[0]; }
    if (edid==2){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    break;
  }
  case planeelementlq:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[3]; }
    if (edid==3){  nodes[0]=nod[3];  nodes[1]=nod[0]; }
    break;
  }
  case planeelementqq:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[4];  nodes[2]=nod[1]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[5];  nodes[2]=nod[2]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[6];  nodes[2]=nod[3]; }
    if (edid==3){  nodes[0]=nod[3];  nodes[1]=nod[7];  nodes[2]=nod[0]; }
    break;
  }
  case planeelementrotlq:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[3]; }
    if (edid==3){  nodes[0]=nod[3];  nodes[1]=nod[0]; }
    break;
  }
  case planeelementsubqt:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1];  nodes[2]=nod[3]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2];  nodes[2]=nod[4]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[0];  nodes[2]=nod[5]; }
    break;
  }
  case axisymmlt:{
    if (edid==0){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==1){  nodes[0]=nod[1];  nodes[1]=nod[2]; }
    if (edid==2){  nodes[0]=nod[2];  nodes[1]=nod[0]; }
    break;
  }
  case linearhex:{
    if (edid==0){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==1){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==2){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==3){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==4){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==5){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==6){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==7){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==8){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==9){   nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==10){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    if (edid==11){  nodes[0]=nod[0];  nodes[1]=nod[1]; }
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function give_edge_nodes (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  destrv (nod);
}



/**
   function returns numbers of nodes defining apropriate element
   
   @param eid - element id
   @param nodes - array containing nodes
   
   22.7.2001
*/
void mechtop::give_elemnodes (long eid,ivector &nodes)
{
  Gtm->give_nodes (eid,nodes);
}

/**
   function extracts all code numbers of actual element
   
   @param eid - element id
   @param cn - array containing all code numbers on element
   
   25.6.2001
*/
void mechtop::give_code_numbers (long eid,long *cn)
{
  Gtm->give_code_numbers (eid,cn);
}

/**
   function extracts code numbers of actual node
   
   @param nid - node id
   @param cn - array containing code numbers on node
   
   7.12.2002
*/
void mechtop::give_node_code_numbers (long nid,long *cn)
{
  Gtm->give_node_code_numbers (nid,cn);
}

/**
   function extracts code numbers of multipliers on actual node
   
   @param nid - node id
   @param lid - layer id
   @param cn - array containing all code numbers on node
   
   7.12.2002
*/
void mechtop::give_mult_code_numbers (long nid,long lid,long *cn)
{
  Gtm->give_mult_code_numbers (nid,lid,cn);
}

/**
   function returns node coordinates of the appropriate element
   
   @param x,y - vectors containing node coordinates
   @param eid - element id
   
   19.7.2001
*/
void mechtop::give_node_coord1d (vector &x,long eid)
{
  Gtm->give_node_coord1d (x,eid);
}

/**
   function returns node coordinates of the appropriate element
   
   @param x,y - vectors containing node coordinates
   @param eid - element id
   
   19.7.2001
*/
void mechtop::give_node_coord2d (vector &x,vector &y,long eid)
{
  Gtm->give_node_coord2d (x,y,eid);
}

/**
   function returns node coordinates of the appropriate element
   
   @param x,z - vectors containing coordinates
   @param eid - element id
   
   19.7.2001
*/
void mechtop::give_node_coord2dxz (vector &x,vector &z,long eid)
{
  Gtm->give_node_coord2dxz (x,z,eid);
}

/**
   function returns node coordinates of the appropriate element
   
   @param x,y,z - vectors containing coordinates
   @param eid - element id
   
   19.7.2001
*/
void mechtop::give_node_coord3d (vector &x,vector &y,vector &z,long eid)
{
  Gtm->give_node_coord3d (x,y,z,eid);
}

/**
   function returns length of the appropriate element

   @param eid - element id

   19.7.2001
*/
double mechtop::give_length (long eid)
{
  long nn = give_nne(eid);
  double l;
  if (nn > 2)
  {
    fprintf(stderr, "\n\nInvalid number of node on element is detected");
    fprintf(stderr, "\n function give_length(), (file %s, line %d)", __FILE__, __LINE__);
  }
  vector x(nn), y(nn), z(nn);
  Gtm->give_node_coord3d (x, y, z, eid);
  l = sqrt(sqr(x[1]-x[0]) + sqr(y[1]-y[0]) + sqr(z[1]-z[0]));
  return l;
}

/**
   function returns area of the appropriate element

   @param eid - element id

   19.7.2001
*/
double mechtop::give_area (long eid)
{
  long nn = give_nne(eid);
  double a = 0.0;
  if (nn < 2)
  {
    fprintf(stderr, "\n\nInvalid number of node on element is detected");
    fprintf(stderr, "\n function give_area(), (file %s, line %d)", __FILE__, __LINE__);
  }
  vector x(nn), y(nn), z(nn);
  Gtm->give_node_coord3d (x, y, z, eid);
  switch (Mt->give_elem_type(eid))
  {
    case planeelementlt:
    case planeelementqt:
    case planeelementrotlt:
      a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
      break;
    case planeelementlq:
    case planeelementqq:
    case planeelementrotlq:
      a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
      break;
    case cctel:
    case dktel:
    case dstel:
      a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]))/2.0;
      break;
    case q4plateel:
      a = ((x[1]-x[0])*(y[2]-y[0])-(x[2]-x[0])*(y[1]-y[0]) + (x[2]-x[0])*(y[3]-y[0])-(x[3]-x[0])*(y[2]-y[0]) ) / 2.0;
      break;
    default :
      fprintf(stderr, "\n\nInvalide element type is required in function give_area (file %s, line %d).\n ", __FILE__, __LINE__);
  }
  return a;
}

/**
   function returns volume of the appropriate element

   @param eid - element id

   19.7.2001
*/
double mechtop::give_volume (long eid)
{
  long nn = give_nne(eid);
  double v = 0.0;
  if (nn < 2)
  {
    fprintf(stderr, "\n\n Invalid number of node on element is detected function give_volume (file %s, line %d).\n", __FILE__, __LINE__);
  }
  vector x(nn), y(nn), z(nn);
  Gtm->give_node_coord3d (x, y, z, eid);
  switch (Mt->give_elem_type(eid))
  {
    case lineartet:
    case linearhex:
    case quadrhex:
      jac_3d (v, x, y, z, 0.0, 0.0, 0.0);
      break;
    default :
      fprintf(stderr, "\n\n Invalide element type is required in function give_volume (file %s, line %d).\n", __FILE__, __LINE__);
  }
  return v;
}

/**
   function returns maximum number of components of strain/stress array at nodes via 
   parameter max_ncompstrn and maximum number of components of strain/stress array on 
   elements via parameter max_ncompstre
   
   @param max_ncompstrn - maximum number of components of strain/stress array at nodes
   @param max_ncompstre - maximum number of components of strain/stress array on elements
*/
void mechtop::give_maxncompstr (long &max_ncompstrn, long &max_ncompstre)
{
  long i;
  max_ncompstrn = max_ncompstre = 0;

  for (i=0;i<ne;i++)
  {
    if (max_ncompstre < Mm->ip[elements[i].ipp[0][0]].ncompstr)
      max_ncompstre = Mm->ip[elements[i].ipp[0][0]].ncompstr;
  }      
  for (i=0;i<nn;i++)
  {
    if (max_ncompstrn < nodes[i].ncompstr)
      max_ncompstrn = nodes[i].ncompstr;
  }      
}

/**
   function returns maximum number of components of other array at nodes via 
   parameter max_nncompo and maximum number of components of eqother array on 
   elements via parameter max_encompo
   
   @param max_encompo - maximum number of components of other array at nodes
   @param max_nncompo - maximum number of components of eqother array on elements
*/
void mechtop::give_maxncompo (long &max_nncompo, long &max_encompo)
{
  long i;
  max_nncompo = max_encompo = 0;

  for (i=0;i<ne;i++)
  {
    if (max_encompo < Mm->ip[elements[i].ipp[0][0]].ncompeqother)
      max_encompo = Mm->ip[elements[i].ipp[0][0]].ncompeqother;
  }      
  for (i=0;i<nn;i++)
  {
    if (max_nncompo < nodes[i].ncompother)
      max_nncompo = nodes[i].ncompother;
  }      
}


/**
   function returns number of local coordinate systems on element
   
   @param nod - array containing node numbers
*/
long mechtop::locsystems (ivector &nod)
{
  long i,j,n;
  n=nod.n;
  j=0;
  for (i=0;i<n;i++){
    j+=nodes[nod[i]].transf;
  }
  return j;
}

/**
   function alocates array containing reactions
*/
void mechtop::comreacnod (void)
{
  long i,j,ndofn;
  
  for (i=0;i<nn;i++){
    ndofn = Gtm->give_ndofn (i);
    for (j=0;j<ndofn;j++){
      if (Gtm->give_dof (i,j)<1){
	nodes[i].react=1;
	if (nodes[i].r==NULL)
	  nodes[i].r = new double [ndofn];
	memset (nodes[i].r,0,ndofn*sizeof(double));
	break;
      }
    }
  }
}

/**
   function marks elements necessary for reaction evaluation
*/
void mechtop::comreacelem (void)
{
  long i,j,k,nne;
  ivector nod;
  
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      nne=Gtm->give_nne (i);
      allocv (nne,nod);
      Gtm->give_nodes (i,nod);
      for (j=0;j<nne;j++){
	k=nod[j];
	if (nodes[k].react==1){
	  elements[i].react=1;
	  break;
	}