elemswitch.cpp

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

      case spring_4 :
      case spring_5 :
      case spring_6 :
	break;
	
	/*
	  case planeelementlt:{
	  Pelt->nod_others (lcid,i,0,0);
	  break;
	  }
	*/

	/*
	  case planeelementrotlt:{
	  Perlt->nod_others (lcid,i,0,0);
	  break;
	  }
	*/
	
      case planeelementlq:{
	Pelq->nod_eqother_ip (lcid,i,0,0);
	break;
      }
	
      case planeelementqq:{
	Peqq->nod_eqother_ip (lcid,i,0,0);
	break;
      }
	/*
	  case planeelementrotlq:{
	  Perlq->nod_stresses (lcid,i,0,0);
	  break;
	  }
	  
	  case planeelementsubqt:{
	  Pesqt->nod_stresses (lcid,i,0,0);
	  break;
	  }
	  
	  case cctel:{
	  Cct->nod_stresses (lcid,i,0,0);
	  break;
	  }
	  case dktel:{
	  Dkt->nod_stresses (lcid,i,0,0);
	  break;
	  }
	  case dstel:{
	  Dst->nod_stresses (lcid,i,0,0);
	  break;
	  }
	  case q4plateel:{
	  Q4pl->nod_stresses (lcid,i,0,0);
	  break;
	  }*/
	
      case axisymmlt:{
	Asymlt->nod_eqother_ip (lcid,i);
	break;
      }
      case axisymmlq:{
	Asymlq->nod_eqother_ip (lcid,i);
	break;
      }
      case axisymmqq:{
	Asymqq->nod_eqother_ip (lcid,i);
	break;
      }
	
	/*
	  case shelltrelem:{
	  Shtr->nod_stresses (lcid,i);
	  break;
	  }
	*/
	
      case lineartet:{
	Ltet->nod_eqother_ip (lcid,i);
	break;
      }
      case linearhex:{
	Lhex->nod_eqother_ip (lcid,i,0,0);
	break;
      }
      case quadrhex:{
	Qhex->nod_eqother_ip (lcid,i,0,0);
	break;
      }
	
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function compute_nodeothers (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
    }
  }
  
  for (i=0;i<Mt->nn;i++){
    Mt->nodes[i].other_averageval (lcid);
  }
}





/**
   function interpolates nodal values to integration points on one element
   
   @param eid - element id
   @param nodval - %vector of nodal values
   @param ipval - %vector of values at integration points
   
   JK, 29.11.2006
*/
void elem_intpointval (long eid,vector &nodval,vector &ipval)
{
  elemtype te;

  te = Mt->give_elem_type (eid);

  switch (te){
  case bar2d:{
    Bar2d->intpointval (eid,nodval,ipval);
    break;
  }
  case bar3d:{
    Bar3d->intpointval (eid,nodval,ipval);
    break;
  }
  case barq2d:{
    Barq2d->intpointval (eid,nodval,ipval);
    break;
  }
    
  case spring_1:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
  case spring_2:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
  case spring_3:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
  case spring_4:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
  case spring_5:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
  case spring_6:{
    Spring->intpointval (eid,nodval,ipval);
    break;
  }
    
  case planeelementlq:{
    Pelq->intpointval (eid,nodval,ipval);
    break;
  }
  case planeelementqq:{
    Peqq->intpointval (eid,nodval,ipval);
    break;
  }
    
  case axisymmlq:{
    Asymlq->intpointval (eid,nodval,ipval);
    break;
  }
  case axisymmqq:{
    Asymqq->intpointval (eid,nodval,ipval);
    break;
  }
    
  case linearhex:{
    Lhex->intpointval (eid,nodval,ipval);
    break;
  }
  case quadrhex:{
    Qhex->intpointval (eid,nodval,ipval);
    break;
  }
    
    /*
      case planeelementrotlq:{
      Perlq->res_internal_forces (lcid,i,ifor);
      break;
      }
      
      case lineartet:{
      Ltet->res_internal_forces (lcid,i,ifor);
      break;
      }
    */
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function elem_intpointval (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
}

/**
   function interpolates nodal values into integration points
   
   @param gv - array containing values at all nodes of the mesh
   @param nmq - type of non-mechanical quantity
   
   21.6.2004, JK
*/
void intpointval (double *gv,nonmechquant nmq,double scale)
{
  long i,j,nne,nip,ipid;
  elemtype te;
  ivector nodes;
  vector nodval,ipval;
  double temp;

  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      nne = Mt->give_nne (i);
      nip = Mt->give_tnip (i);
      
      allocv (nne,nodes);
      allocv (nne,nodval);
      allocv (nip,ipval);
      
      Mt->give_elemnodes (i,nodes);
      
      for (j=0;j<nne;j++){
	nodval[j]=gv[nodes[j]]*scale;
      }
      
      
      //  interpolation of nodal values to integration points
      elem_intpointval (i,nodval,ipval);
      
      //added by TKr 4.11.2007
      //averaged values for element with linear approximation functions 
      /* if((nmq == temperature) || (nmq == moist)){
	 te = Mt->give_elem_type (i);
	 switch (te){
	 case bar2d:
	 case bar3d:
	 case planeelementlt:
	 case planeelementrotlt:
	 case planeelementlq:
	 case axisymmlt:
	 case axisymmlq:
	 case lineartet:
	 case linearhex:{
	 temp = 0.0;
	 for (j=0;j<nip;j++){
	 temp += ipval[j];
	 }
	 for (j=0;j<nip;j++){
	 ipval[j] = temp/nip;
	 } 
	 }
	 default:{
	 }
	 }
	 }
      */

      //  number of the first integration point
      ipid=Mt->elements[i].ipp[0][0];
      
      for (j=0;j<nip;j++){
	switch (nmq){
	case temperature:{
	  Mm->tempr[ipid]=ipval[j];
	  ipid++;
	  break;
	}
	case moist:{
	  Mm->moist[ipid]=ipval[j];
	  ipid++;
	  break;
	}
	default:{
	  fprintf (stderr,"\n\n unknown type of non-mechanical quantity is required in function intpointval (file %s, line %d).\n",__FILE__,__LINE__);
	}
	}
	
      }
      
      destrv (nodval);
      destrv (nodes);
      destrv (ipval);
    }
  }
}


/**
   function interpolates nodal values into integration points
   
   @param gv - array containing values at all nodes of the mesh
   @param nmq - type of non-mechanical quantity
   
   21.6.2004, JK
*/
void intpointval2 (double *gv,nonmechquant nmq)
{
  long i,j,nne,nip,ipp;
  elemtype te;
  ivector nodes;
  vector nodval,ipval;
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      te = Mt->give_elem_type (i);
      nne = Mt->give_nne (i);
      nip = Mt->give_tnip (i);
      
      allocv (nne,nodes);
      allocv (nne,nodval);
      allocv (nip,ipval);
      
      Mt->give_elemnodes (i,nodes);
      
      for (j=0;j<nne;j++){
	nodval[j]=gv[nodes[j]];
      }
      
      switch (te){
      case barq2d:{
	Barq2d->intpointval2 (i,nodval,ipval);
	break;
      }
      case barq3d:{
	Barq3d->intpointval2 (i,nodval,ipval);
	break;
      }
	
	
      case spring_1:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
      case spring_2:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
      case spring_3:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
      case spring_4:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
      case spring_5:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
      case spring_6:{
	Spring->intpointval (i,nodval,ipval);
	break;
      }
	
	/*
	  case planeelementlq:{
	  Pelq->intpointval (i,nodval,ipval);
	  break;
	  }
	*/
	
      case planeelementqq:{
	Peqq->intpointval2 (i,nodval,ipval);
	break;
      }
      case axisymmqq:{
	Asymqq->intpointval2 (i,nodval,ipval);
	break;
      }
	
      case quadrhex:{
	Qhex->intpointval2 (i,nodval,ipval);
	break;
      }
	
	/*
	  case planeelementrotlq:{
	  Perlq->res_internal_forces (lcid,i,ifor);
	  break;
	  }
	  
	  case lineartet:{
	  Ltet->res_internal_forces (lcid,i,ifor);
	  break;
	  }
	  case linearhex:{
	  Lhex->res_internal_forces (lcid,i,ifor);
	  break;
	  }
	  
	*/
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function intpointval (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
      
      //  number of the first integration point
      ipp=Mt->elements[i].ipp[0][0];
      
      for (j=0;j<nip;j++){
	switch (nmq){
	case temperature:{
	  Mm->tempr[ipp]=ipval[j];
	  ipp++;
	  break;
	}
	case moist:{
	  Mm->moist[ipp]=ipval[j];
	  ipp++;
	  break;
	}
	default:{
	  fprintf (stderr,"\n\n unknown type of non-mechanical quantity is required in function intpointval (file %s, line %d).\n",__FILE__,__LINE__);
	}
	}
	
      }
      
      destrv (nodval);
      destrv (nodes);
      destrv (ipval);
    }
  }
}


/**
   function interpolates arbitrary quantity defined by its
   nodal values to required inner point
   function returns required value
   
   @param eid - element id
   @param nodval - array containing nodal values
   @param coord - array containing coordinates
   
   JK, 23.2.2002
*/
double interpolate (long eid,double *nodval,double *coord)
{
  long i,nne;
  double val;
  elemtype te;
  vector areacoord(3),volcoord(4),nv;
  
  te = Mt->give_elem_type (eid);
  nne = Mt->give_nne (eid);
  allocv (nne,nv);
  for (i=0;i<nne;i++){
    nv[i]=nodval[i];
  }

  switch (te){
  case planeelementlt:{
    areacoord[0]=coord[0];  areacoord[1]=coord[1];  areacoord[2]=1.0-areacoord[0]-areacoord[1];
    val=Pelt->approx (areacoord,nv);
    break;
  }
  case planeelementqt:{
    areacoord[0]=coord[0];  areacoord[1]=coord[1];  areacoord[2]=1.0-areacoord[0]-areacoord[1];
    val=Peqt->approx (areacoord[0],areacoord[1],nv);
    break;
  }
  case planeelementrotlt:{
    areacoord[0]=coord[0];  areacoord[1]=coord[1];  areacoord[2]=1.0-areacoord[0]-areacoord[1];
    val=Perlt->approx (areacoord,nv);
    break;
  }
  case planeelementlq:{
    val=Pelq->approx (coord[0],coord[1],nv);
    break;
  }
  case planeelementqq:{
    val=Peqq->approx (coord[0],coord[1],nv);
    break;
  }
  case planeelementrotlq:{
    val=Perlq->approx (coord[0],coord[1],nv);
    break;
  }
  case lineartet:{
    volcoord[0]=coord[0];  volcoord[1]=coord[1];  volcoord[2]=coord[2];  volcoord[3]=1.0-volcoord[0]-volcoord[1]-volcoord[2];
    val=Ltet->approx (volcoord,nv);
    break;
  }
  case linearhex:{
    val=Lhex->approx (coord[0],coord[1],coord[2],nv);
    break;
  }
  case quadrhex:{
    val=Qhex->approx (coord[0],coord[1],coord[2],nv);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function interpolate (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  destrv (nv);
  
  return val;
}


/**
   function computes volume appropriate to integration point
   
   2.3.2004, JK
*/
void ipvolume ()
{
  long i;
  elemtype te;
  
  if (Mm->ipv == NULL)
    delete [] Mm->ipv;
  
  Mm->ipv = new double [Mm->tnip];
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      te = Mt->give_elem_type (i);
      
      switch (te){
      case planeelementlt:{
	Pelt->ipvolume (i,0,0);
	break;
      }
      case planeelementqt:{
	Peqt->ipvolume (i,0,0);
	break;
      }
      case planeelementrotlt:{
	Perlt->ipvolume (i,0,0);
	break;
      }
      case planeelementlq:{
	Pelq->ipvolume (i,0,0);
	break;
      }
      case planeelementqq:{
	Peqq->ipvolume (i,0,0);
	break;
      }
      case planeelementrotlq:{
	Perlq->ipvolume (i,0,0);
	break;
      }
      case lineartet:{
	Ltet->ipvolume (i,0,0);
	break;
      }
      case linearhex:{
	Lhex->ipvolume (i,0,0);
	break;
      }
      case quadrhex:{
	Qhex->ipvolume (i,0,0);
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function ipvolume (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
    }
    
  }
}