mechtop.cpp

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

      }
      destrv (nod);
    }
  }
}

/**
   function marks elements necessary for reaction evaluation and
   alocates array containing reactions
*/
void mechtop::comreac (void)
{
  comreacnod ();
  comreacelem ();
}

/**
   function computes nodal values with help of least square problem solution
   
   @param nodes - array containing nodes of element
   @param nx,ny,nz - arrays containing natural coordinates of nodes
   @param lhs - array of coefficients of linear functions (solution of least square problem)
   @param dim - problem dimension
   @param ncomp - number of stored components
   @param pos - position of first component in nodal strain array

   10.5.2002
*/
void mechtop::strain_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
				   double *lhs,long dim,long fi,long ncomp,long lcid)
{
  long i,j;
  vector eps;
  
  allocv (ncomp,eps);

  if (dim==1){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	eps[j] = lhs[j*2]+lhs[j*2+1]*nx[i];
      }
      nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
    }
  }
  if (dim==2){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	eps[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
      }
      nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
    }
  }
  if (dim==3){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	eps[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
      }
      nodes[nod[i]].storestrain (lcid,fi,ncomp,eps);
    }
  }
  
  destrv (eps);
}

/**
   function computes nodal values with help of least square problem solution
   
   10.5.2002
*/
void mechtop::stress_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
				   double *lhs,long dim,long fi,long ncomp,long lcid)
{
  long i,j;
  vector sig;
  
  allocv (ncomp,sig);
  if (dim==2){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	//Mt->nodes[nodes[i]].stress[pos+j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
	sig[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
      }
      nodes[nod[i]].storestress (lcid,fi,ncomp,sig);

    }
  }
  if (dim==3){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	//Mt->nodes[nodes[i]].stress[pos+j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
	sig[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
      }
      nodes[nod[i]].storestress (lcid,fi,ncomp,sig);
    }
  }
  destrv (sig);
}

/**
   function computes nodal values with help of least square problem solution
   
   10.5.2002
*/
void mechtop::other_nodal_values (ivector &nod,vector &nx,vector &ny,vector &nz,
				   double *lhs,long dim,long fi,long ncomp,long lcid)
{
  long i,j;
  vector other;
  
  allocv (ncomp,other);
  if (dim==1){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	other[j] = lhs[j*2]+lhs[j*2+1]*nx[i];
      }
      nodes[nod[i]].storeother (lcid,fi,ncomp,other);
    }
  }
  if (dim==2){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	other[j] = lhs[j*3]+lhs[j*3+1]*nx[i]+lhs[j*3+2]*ny[i];
      }
      nodes[nod[i]].storeother (lcid,fi,ncomp,other);

    }
  }
  if (dim==3){
    for (i=0;i<nod.n;i++){
      for (j=0;j<ncomp;j++){
	other[j] = lhs[j*4]+lhs[j*4+1]*nx[i]+lhs[j*4+2]*ny[i]+lhs[j*4+3]*nz[i];
      }
      nodes[nod[i]].storeother (lcid,fi,ncomp,other);
    }
  }
  destrv (other);
}




/**
  The function converts code numbers stored at nodes to code numbers
  stored on elements. It is used for mixed 2D mesh of plane elements,
  bars and springs.

  Author : Tomas Koudelka - [e-mail]
*/
void mechtop::store_code_num_elem(void)
{
  long i, j, nne, ndofe, ndofn;
  elemtype te;
  ivector nod;
  long *cn;
  for(i=0; i < ne; i++)
  {
    te    = give_elem_type(i);
    nne   = give_nne(i);
    ndofe = give_ndofe(i);
    allocv(nne, nod);
    give_elemnodes(i, nod);
    switch(te)
    {
      case beam2d:
        break;
      case spring_1:
      {
        ndofn = give_ndofn(nod[0]);
        if (ndofn < 1)
        {
          fprintf(stderr, "\n\n Direction of spring_1 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[0] = cn[0];
          delete [] cn; 
        }
        break;
      }
      case spring_2:
      {
        ndofn = give_ndofn(nod[0]);
        if (ndofn < 2)
        {
          fprintf(stderr, "\n\n Direction of spring_2 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[1] = cn[1];
          delete [] cn; 
        }
        break;
      }
      case spring_3:
      {
        ndofn = give_ndofn(nod[0]);
        if (ndofn < 1)
        {
          fprintf(stderr, "\n\n Direction of spring_3 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[2] = cn[2];
          delete [] cn; 
        }
        break;
      }
      case spring_4:
      {
        ndofn = give_ndofn(nod[0]);
        if (ndofn < 1)
        {
          fprintf(stderr, "\n\n Direction of spring_4 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[3] = cn[3];
          delete [] cn; 
        }
        break;
      }
      case spring_5:
      {
        ndofn = give_ndofn(nod[0]);
        if (ndofn < 1)
        {
          fprintf(stderr, "\n\n Direction of spring_5 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[4] = cn[4];
          delete [] cn; 
        }
        break;
      }
      case spring_6: 
      {
        if (ndofn < 1)
        {
          fprintf(stderr, "\n\n Direction of spring_6 does not correspond to number\n");
          fprintf(stderr, " of dofs(%ld) in node %ld\n", ndofn, nod[0]);
          abort();
        }
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofn];
          cn = new long[ndofn];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofn);
          Mt->give_node_code_numbers(nod[0], cn);
          Gtm->gelements[i].cn[5] = cn[5];
          delete [] cn; 
        }
        break;
      }
      case bar2d :
      case barq2d:
      case planeelementlt:
      case planeelementqt:
      case planeelementlq:
      case planeelementqq:
        if (Gtm->gelements[i].cne == 0)
        {
          Gtm->gelements[i].cne = 1;
          Gtm->gelements[i].cn = new long[ndofe];
          memset(Gtm->gelements[i].cn, 0, sizeof(*Gtm->gelements[i].cn)*ndofe);
          for (j=0; j<nne; j++)
          {
            ndofn = Gtm->give_ndofn(nod[j]);
            cn = new long[ndofn];
            Mt->give_node_code_numbers(nod[j], cn);
            Gtm->gelements[i].cn[2*j]   = cn[0];
            Gtm->gelements[i].cn[2*j+1] = cn[1];
            delete [] cn; 
          }
        }
        break;
      default:
        fprintf(stderr, "\n\n Unknown type of element is required in function mechtop::store_code_num_elem\n"); 
        fprintf(stderr, " (file %s, line : %d)\n\n", __FILE__, __LINE__); 
    }
  }
}


/**
   function allocates arrays defined on nodes
   
   arrays strain, stress, other, ncontr_strain, ncontr_stress, ncontr_other are allocated on each node
   
   JK, 28.11.2006
*/
void mechtop::alloc_nodes_arrays (void)
{
  long i,j,nne,ncomp, ncompo;
  ivector nod;
  
  if (Mp->strainpos==2 || Mp->strainpos==3 || Mp->stresspos==2 || Mp->stresspos==3 || Mp->otherpos==2){
    for (i=0;i<ne;i++){
      if (Gtm->leso[i]==1){
	//  number of nodes on element
	nne = give_nne (i);
	allocv (nne,nod);
	give_elemnodes (i,nod);
	//  number of strain/stress components
	ncomp=give_ncomp (i);
	//  number of components in array other
	ncompo=Mm->ip[elements[i].ipp[0][0]].ncompeqother;     
	
	for (j=0;j<nne;j++){
	  if (nodes[nod[j]].strain==NULL){
	    //  the main loop is over elements because ncomp is defined on elements
	    //  this condition skips already allocated nodes
	    //  nodes are called several times
	    nodes[nod[j]].alloc(ncomp,ncompo,Mb->nlc);
	  }
	  else{
	    if (ncompo > nodes[nod[j]].ncompother){
	      //  if different number of components are used, the maximum number has to be used
	      nodes[nod[j]].realloc(ncompo,Mb->nlc);
	    }
	  }
 	}
	destrv (nod);
      }
    }
  }
}

/**
   function allocates nodes and elements for preprocessor

   @param nn - number of nodes
   @param ne - number of elements

   13.3.2002
*/
void mechtop::alloc_prep (long nn,long ne)
{
  nodes = new node [nn];
  elements = new element [ne];
}

/**
   function allocates array for meaning of DOFs
   
   1.2.2005, JK
*/
void mechtop::alloc_meaning ()
{
  long i;
  
  for (i=0;i<nn;i++){
    nodes[i].alloc_meaning (i);
  }
}

/**
   function allocates all necessary arrays for problems with changing
   number of nodes, elements and DOFs
   
   7.11.2006, JK
*/
void mechtop::alloc_growstr ()
{
  long i;
  
  for (i=0;i<nn;i++){
    nodes[i].alloc_growstr (i);
  }
  for (i=0;i<ne;i++){
    elements[i].alloc_