elemswitch.cpp

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

#include "elemswitch.h"
#include "genfile.h"
#include "global.h"
#include "elemhead.h"
#include "node.h"
#include "element.h"
#include "intpoints.h"

// *******************
// *******************
// ****  STRAINS  ****
// *******************
// *******************

/**
   function computes strains at integration points
   
   all components at all integration points are evaluated
   
   @param lcid - load case id
   
   JK, 28.11.2006
*/
void compute_ipstrains (long lcid)
{
  long i,ne,te;
  
  //  strain array cleaning
  Mm->cleanstrain ();
  
  //  number of elements
  ne=Mt->ne;
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      //  only elements switched on are taken into account
      
      //  type of element
      te = Mt->give_elem_type (i);
      
      switch (te){
	
      case bar2d:{
	Bar2d->res_ip_strains (lcid,i);
	break;
      }
      case bar3d:{
	Bar3d->res_ip_strains (lcid,i);
	break;
      }
      case barq2d:{
	Barq2d->res_mainip_strains (lcid,i);
	break;
      }
      case barq3d:{
	Barq3d->res_mainip_strains (lcid,i);
	break;
      }
	
      case beam2d:{
	Beam2d->nodal_displ (lcid,i);
	break;
      }
      case beam3d:{
	Beam3d->nodal_displ (lcid,i);
	break;
      }
      case beamg3d:{
	Beam3dg->nodal_displ (i,lcid);
	break;
      }
      case subsoilbeam:{
	Sbeam->strains (lcid, i, 0, 0);
	break;
      }
      case beam2dsp:{
	Spbeam2d->res_mainip_strains (lcid,i);
	break;
      }
	
      case spring_1:{
	Spring->strains(i, lcid);
	break;
      }
      case spring_2:{
	Spring->strains(i, lcid);
	break;
      }
      case spring_3:{
	Spring->strains(i, lcid);
	break;
      }
      case spring_4:{
	Spring->strains(i, lcid);
	break;
      }
      case spring_5:{
	Spring->strains(i, lcid);
	break;
      }
      case spring_6:{
	Spring->strains(i, lcid);
	break;
      }
	
	
      case planeelementlt:{
	Pelt->res_ip_strains (lcid,i);
	break;
      }
      case planeelementqt:{
	Peqt->res_mainip_strains (lcid,i);
	break;
      }
	
      case planeelementrotlt:{
	Perlt->res_mainip_strains (lcid,i);
	break;
      }
	
      case planeelementlq:{
	Pelq->res_ip_strains (lcid,i);
	break;
      }
      case planeelementqq:{
	Peqq->res_ip_strains (lcid,i);
	break;
      }
      case planeelementrotlq:{
	Perlq->res_mainip_strains (lcid,i);
	break;
      }
	
      case planeelementsubqt:{
	Pesqt->res_mainip_strains (lcid,i);
	break;
      }
      case planequadcontact:{
	Pqcon->res_mainip_strains (lcid,i);
	break;
      }
	
      case cctel:{
	Cct->res_mainip_strains (lcid,i);
	break;
      }
      case dktel:{
	Dkt->res_mainip_strains (lcid,i);
	break;
      }
      case dstel:{
	Dst->res_mainip_strains (lcid,i);
	break;
      }
      case q4plateel:{
	Q4pl->res_ip_strains (lcid,i);
	break;
      }
	
      case subsoilplatetr:{
	break;
      }
	
      case subsoilplateq:{
	Splq->res_mainip_strains (lcid,i);
	break;
      }
	
	
      case shelltrelem:{
	Shtr->res_mainip_strains (lcid,i);
	break;
      }
	
	
      case axisymmlt:{
	Asymlt->res_mainip_strains (lcid,i);
	break;
      }
      case axisymmlq:{
	Asymlq->res_allip_strains (lcid,i);
	break;
      }
      case axisymmqq:{
	Asymqq->res_allip_strains (lcid,i);
	break;
      }
	
	
      case lineartet:{
	Ltet->res_mainip_strains (lcid,i);
	break;
      }
      case linearhex:{
	Lhex->res_ip_strains (lcid,i);
	break;
      }
      case quadrhex:{
	Qhex->res_ip_strains (lcid,i);
	break;
      }
	
      case particleelem:{
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function compute_ipstrains (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
    }
  }

  //  function modifies strains
  //  eigenstrains are subtracted from the total strains
  //  it works only for the first load case
  Mm->eigstrmod ();
  
}




/**
   function computes strains at element nodes
   nodal strains are averaged
   in case of bar and beam elements be carefull if averaging is reasonable
   
   @param lcid - load case id
   
   JK, 28.11.2006
*/
void compute_nodestrains (long lcid)
{
  long i;
  elemtype te;
  
  //  for open dx
  Mp->detnodstrain = 1;

  for (i=0;i<Mt->nn;i++){
    Mt->nodes[i].ncontr_strain[lcid]=0;
    Mt->nodes[i].nullstrain (lcid);
  }
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      te = Mt->give_elem_type (i);
      
      switch (te){
	
      case bar2d:{
	//Bar2d->nod_strains_ip (lcid, i, 0, 0);
	break;
      }
      case barq2d:{
	Barq2d->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case bar3d:{
	Bar3d->nod_strains_ip (lcid, i, 0, 0);
	break;
      }
      case barq3d:{
	Barq3d->nod_strains_ip (lcid,i,0,0);
	break;
      }
	
      case beam2d:{
	Beam2d->nodal_displ (lcid,i);
	break;
      }
      case beam3d:{
	Beam3d->nodal_displ (lcid,i);
	break;
      }

	
      case spring_1 :
      case spring_2 :
      case spring_3 :
      case spring_4 :
      case spring_5 :
      case spring_6 :
	break;
      case planeelementlt:{
	Pelt->nod_strains (lcid,i,0,0);
	break;
      }
      case planeelementqt:{
	Peqt->nod_strains (lcid,i,0,0);
	break;
      }
      case planeelementrotlt:{
	Perlt->nod_strains (lcid,i,0,0);
	break;
      }
	
      case planeelementlq:{
	Pelq->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case planeelementqq:{
	Peqq->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case planeelementrotlq:{
	//Perlq->nod_strains (lcid,i,0,0);
	break;
      }
	
      case planeelementsubqt:{
	Pesqt->nod_strains (lcid,i,0,0);
	break;
      }
	
      case cctel:{
	//Cct->nod_strains (lcid,i,0,0);
	break;
      }
      case dktel:{
	//Dkt->nod_strains (lcid,i,0,0);
	break;
      }
      case dstel:{
	//Dst->nod_strains (lcid,i,0,0);
	break;
      }
      case q4plateel:{
	//Q4pl->nod_strains (lcid,i,0,0);
	break;
      }
	
      case subsoilplatetr:{
	break;
      }
	
      case subsoilplateq:{
	//Splq->nod_strains (lcid,i,0,0);
	break;
      }
	
      case axisymmlt:{
	Asymlt->nod_strains_ip (lcid,i);
	break;
      }
      case axisymmlq:{
	Asymlq->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case axisymmqq:{
	Asymqq->nod_strains_ip (lcid,i);
	break;
      }
	
      case shelltrelem:{
	Shtr->nod_strains (lcid,i);
	break;
      }
	
      case lineartet:{
	Ltet->nod_strains (lcid,i,0,0);
	break;
      }
      case linearhex:{
	Lhex->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case quadrhex:{
	Qhex->nod_strains_ip (lcid,i,0,0);
	break;
      }
      case particleelem:{
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function compute_nodestrains (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
    }
  }

  //  averaging of nodal strains
  for (i=0;i<Mt->nn;i++){
    //if (Gtm->lnso[i]==1){
    Mt->nodes[i].strain_averageval (lcid);
    //}
  }
}


/**
   function computes strains at strain points
   
   @param lcid - load case id
   
   JK
*/
void computestrains (long lcid)
{
  
  //  function computes strains at integration points
  compute_ipstrains (lcid);
  
  /*
  long j,nip,ipp,ncomp,k;
  fprintf (Out,"\n\n\n\n kontrola napeti a sil \n\n");
  for (i=0;i<Mt->ne;i++){
    fprintf (Out,"\n prvek %ld  ",i+1);
    ipp=Mt->elements[i].ipp[0][0];
    fprintf (Out,"ipp  %ld  ",ipp);
    nip=Mt->give_tnip (i);
    fprintf (Out,"nip  %ld  ",nip);
    for (j=0;j<nip;j++){
      fprintf (Out,"\n int. bod %ld  ",j+1);
      ncomp=Mm->ip[ipp].ncompstr;
      for (k=0;k<ncomp;k++){
	fprintf (Out,"  %lf",Mm->ip[ipp].strain[k]);
      }
      ipp++;
    }
  }
  */
  

  /*
  //  number of elements
  ne=Mt->ne;
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      
      //  type of element
      te = Mt->give_elem_type (i);
      
      switch (te){
      case bar2d:{
	Bar2d->strains (lcid,i,0,0);
	break;
      }
      case bar3d:{
	Bar3d->strains (lcid,i,0,0);
	break;
      }
      case barq2d:{
	Barq2d->strains (lcid,i,0,0);
	break;
      }
      case barq3d:{
	Barq3d->strains (lcid,i,0,0);
	break;
      }
      case beam2d:{
	Beam2d->strains (i,lcid);
	break;
      }
      case beam3d:{
	Beam3d->nodal_displ (i,lcid);
	break;
      }
      case subsoilbeam:{
	Sbeam->strains (lcid, i, 0, 0);
	break;
      }
      case spring_1:
      case spring_2:
      case spring_3:
      case spring_4:
      case spring_5:
      case spring_6:
	{
	  Spring->strains (i,lcid);
	  break;
	}
	
      case planeelementlt:{
	Pelt->strains (lcid,i,0,0);
	break;
      }
      case planeelementqt:{
	Peqt->strains (lcid,i,0,0);
	break;
      }
      case planeelementrotlt:{
	Perlt->strains (lcid,i,0,0);
	break;
      }
	
      case planeelementlq:{
	Pelq->strains (lcid,i,0,0);
	break;
      }
      case planeelementqq:{
	Peqq->strains (lcid,i,0,0);
	break;
      }
      case planeelementrotlq:{
	Perlq->strains (lcid,i,0,0);
	break;
      }
	
      case planeelementsubqt:{
	Pesqt->strains (lcid,i,0,0);
	break;
      }
	
      case cctel:{
	Cct->strains (lcid,i,0,0);
	break;
      }
      case dktel:{
	Dkt->strains (lcid,i,0,0);
	break;
      }
      case dstel:{
	Dst->strains (lcid,i,0,0);
	break;
      }
      case q4plateel:{
	Q4pl->strains (lcid,i,0,0);
	break;
      }
      case subsoilplateq:{
	Splq->strains (lcid,i,0,0);
	break;
      }
	
	
      case axisymmlt:{
	Asymlt->strains (lcid,i,0,0);
	break;
      }
      case axisymmlq:{
	Asymlq->strains (lcid,i,0,0);
	break;
      }
      case axisymmqq:{
	Asymqq->strains (lcid,i,0,0);
	break;
      }
	
      case shelltrelem:{
	Shtr->strains (lcid,i);
	break;
      }
	
      case lineartet:{
	Ltet->strains (lcid,i,0,0);
	break;
      }
      case linearhex:{
	Lhex->strains (lcid,i,0,0);
	break;
      }
      case quadrhex:{
	Qhex->strains (lcid,i,0,0);
	break;
      }
	
      case particleelem:{
	break;
      }
	
      default:{
	fprintf (stderr,"\n unknown element type is required in function computestrains (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
    }
  }
  */
}



// ********************
// ********************
// ****  STRESSES  ****
// ********************
// ********************

/**
   function computes all stresses at integration points