globmat.cpp

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

#include "globmat.h"
#include "genfile.h"
#include "global.h"
#include "elemhead.h"
#include "loadcase.h"
#include "dloadcase.h"
#include "dloadpd.h"
#include "node.h"
#include "element.h"
#include "intpoints.h"
#include "problem.h"
#include "mechmat.h"
#include "creepb.h"
#include "creepbs.h"
#include "mechcrsec.h"
#include "elemswitch.h"


/**
   function assembles stiffness %matrix
   
   @param lcid - load case id
   
   JK
*/
void stiffness_matrix (long lcid)
{
  long i,j,nl,ndofe,ndofn,nmult,nid,*cn,*cna;
  matrix lm,lmt;
  
  if (Smat==NULL)  Smat = new gmatrix;
  Smat->ts=Mp->tstorsm;
  Smat->initiate (Gtm,Ndofm,Mespr);
  //gmat_initialize (Mp,Smat);
  //Mp->provizorium (Smat);
  Smat->setval (Mp->ssle);

  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      ndofe = Mt->give_ndofe (i);
      
      allocm (ndofe,ndofe,lm);
      //  computation of stiffness matrix of one element
      stiffmat (lcid,i,lm);
      
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);
      
      //  localization of matrix of one element to matrix of problem
      Smat->localize (lm,cn,i);
      
      destrm (lm);
      delete [] cn;
    }
  }
  /*
  DSS::SparseMatrixF sm((unsigned long)Ndofm,Smat->scr->a,(unsigned long*)Smat->scr->ci,(unsigned long*)Smat->scr->adr);
  FILE *out;
  out = fopen ("matice.sm","wb");
  sm.SaveMatrix(out);
  */

  if (Mp->tprob==layered_linear_statics){
    for (i=0;i<Mt->nln;i++){
      nl=Gtm->lgnodes[i].nl;
      ndofn=Gtm->lgnodes[i].ndofn;
      nmult=Gtm->lgnodes[i].nmult;
      //allocm (nmult,ndofn*2,lm);
      allocm (ndofn*2,nmult,lm);
      allocm (nmult,ndofn*2,lmt);
      cn = new long [ndofn*2];
      cna = new long [nmult];
      
      for (j=1;j<nl;j++){
	constr_matrix (i,j,lm);

	nid=Gtm->lgnodes[i].nodes[j-1];
	Mt->give_node_code_numbers (nid,cn);
	nid=Gtm->lgnodes[i].nodes[j];
	Mt->give_node_code_numbers (nid,cn+ndofn);
	Mt->give_mult_code_numbers (i,j,cna);
	
	tranm (lm,lmt);
	
	Smat->glocalize (lmt,cna,cn);
	Smat->glocalize (lm,cn,cna);
	

      }
      destrm (lm);
      destrm (lmt);
      delete [] cn;
      delete [] cna;
    }
  }
  
  Smat->prepmat (Gtm,0.0,Mespr);
  
  //Smat->printdiag(Out);
  //Smat->printmat(Out);
  
  
  
}

/**
   function assembles mass %matrix
   
   @param lcid - load case id
   
   JK
*/
void mass_matrix (long lcid)
{
  if (Mmat==NULL)  Mmat = new gmatrix;
  Mmat->ts=Mp->tstorsm;
  Mmat->initiate (Gtm,Ndofm,Mespr);
  Mmat->setval (Mp->ssle);
  
  
  long i,j,k,ndofn,ndofe,idcs,*cn;
  double m;
  matrix lm;
  crsectype cr;
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      ndofe = Mt->give_ndofe (i);
      
      allocm (ndofe,ndofe,lm);
      //  computation of mass matrix of one element
      massmat (lcid,i,lm);
      
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);
      
      //  localization of matrix of element to matrix of problem
      Mmat->localize (lm,cn,i);
      
      destrm (lm);
      delete [] cn;
    }
  }
  
  for (i=0;i<Mt->nn;i++){
    cr = Mt->nodes[i].crst;
    if (cr==nocrosssection){
      continue;
    }
    else{
      ndofn = Mt->give_ndofn (i);
      for (j=0;j<ndofn;j++){
	k = Mt->give_dof (i,j);
	if (k>0){
	  idcs = Mt->nodes[i].idcs;
	  m=Mc->give_weight (cr,idcs);
	  Mmat->add_entry (m,k-1,k-1);
	}
      }
    }
  }
  
  
  
  Mmat->prepmat (Gtm,0.0,Mespr);

}

/**
   function assembles initial stiffness %matrix

   it is also called geometric stiffness %matrix
   (in R.W. Clough and J. Penzien: Dynamics of structures)
   
   @param lcid - load case id
   
   JK
*/
void initial_stiffness_matrix (long lcid)
{
  if (Ismat==NULL)  Ismat = new gmatrix;
  Ismat->ts=Mp->tstorsm;
  Ismat->initiate (Gtm,Ndofm,Mespr);
  Ismat->setval (Mp->ssle);
  
  
  long i,ndofe,*cn;
  matrix lm;
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){

      ndofe = Mt->give_ndofe (i);
      
      allocm (ndofe,ndofe,lm);
      
      //  computation of initial stiffness matrix of one element
      initstiffmat (lcid,i,lm);
      
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);
      
      //  localization of matrix of one element to matrix of problem
      Ismat->localize (lm,cn,i);
      
      destrm (lm);
      delete [] cn;
    }
  }
  
  Ismat->prepmat (Gtm,0.0,Mespr);
  
}

/**
   function assembles damping %matrix
   
   JK, 10.8.2005
*/
void damping_matrix ()
{
  if (Dmat==NULL)  Dmat = new gmatrix;
  Dmat->ts=Mp->tstorsm;
  Dmat->initiate (Gtm,Ndofm,Mespr);
  Dmat->setval (Mp->ssle);
  
  switch (Mp->damp){
  case nodamp:{
    break;
  }
  case massdamp:{
    Dmat->addgm (Mp->dmass,*Mmat);
    break;
  }
  case stiffdamp:{
    Dmat->addgm (Mp->dstiff,*Smat);
    break;
  }
  case rayleighdamp:{
    Dmat->addgm (Mp->dmass,*Mmat);
    Dmat->addgm (Mp->dstiff,*Smat);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of damping is required in function probdesc::print (file %s, line %d),\n",__FILE__,__LINE__);
  }
  }
  

  
  Dmat->prepmat (Gtm,0.0,Mespr);
  
}


/**
   function computes %vector of internal forces
   
   @param lcid - load case id
   @param intfor - internal forces
   
   28.7.2001
*/
void internal_forces (long lcid,double *intfor)
{
  //fprintf (Out,"\n\n new internal forces  (globmat.cpp, line 188)\n");

  if (Mp->matmodel==local){
    loc_internal_forces (lcid,intfor);
  }
  if (Mp->matmodel==nonlocal){
    nonloc_internal_forces (lcid,intfor);
  }
  
  //  contribution to internal forces by Lagrange multipliers
  //  it is used in problems with floating subdomains solved by the FETI method
  //  JK, 22.6.2006
  if (Mp->tprob == lin_floating_subdomain || Mp->tprob == nonlin_floating_subdomain){
    lagrmultcontr_intforces (lcid,intfor);
  }
  
  //  indicator of strain computation
  //  it means, strains at integration points have been computed
  Mp->strainstate=1;
  //  indicator of stress computation
  //  it means, stresses at integration points have been computed
  Mp->stressstate=1;
  //  indicator of computation of other array
  //  it means, stresses at integration points have been computed
  Mp->otherstate=1;
}

/**
   function computes vector of internal forces
   local material models are used
   
   @param lcid - load case id
   @param intfor - internal forces
   
   28.7.2001
*/
void loc_internal_forces (long lcid,double *intfor)
{
  long i,j,ne,ndofe;
  ivector cn;
  vector ifor;
  
  //  function computes strains at integration points
  compute_ipstrains (lcid);
    
  Mt->nodedisplacements ();
  
  ne=Mt->ne;
  
  nullv (intfor,Ndofm);
  
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      
      ndofe=Mt->give_ndofe (i);
      allocv (ndofe,ifor);
      allocv (ndofe,cn);
      
      elem_internal_forces (i,lcid,ifor);
      
      
      /*      if (Mp->phase == 2)//debug!!!
	      if (Mp->time ==  10798800.0){//debug!!!
	      fprintf (Out,"\n element = %ld\n",i);
	      for (j=0;j<ndofe;j++){
	      fprintf (Out,"%e\n",ifor[j]);
	      }
	      }
      */

      Mt->give_code_numbers (i,cn.a);

      /*
      long p = -1;
      for (long j = 0; j < cn.n; j++)
	{
	  if ((cn[j] > 0) && (cn[j] < 33))
	    p = 1;
	}
      */
      locglob (intfor,ifor.a,cn.a,ndofe);
      destrv (ifor);
      destrv (cn);
    }
  }
  
}

/**
   function computes vector of internal forces
   local material models are used

   @param lcid - load case id
   @param intfor - internal forces

   28.7.2001
*/
void nonloc_internal_forces (long lcid,double *intfor)
{
  long i,ne,te,ndofe;
  ivector cn;
  vector ifor;

  //  function computes strains at integration points
  compute_ipstrains (lcid);

  ne=Mt->ne;
  nullv (intfor,Ndofm);
  Mp->phase=1;

  for (i = 0; i < ne; i++){
    if (Gtm->leso[i]==1){
      te = Mt->give_elem_type (i);
      switch (te)
	{
	case planeelementlt :
	  Pelt->compute_nlstress (lcid,i,0,0);
	  break;
	case lineartet:
	  Ltet->local_values (lcid, i, 0, 0);
	  break;
	case linearhex:
	  Lhex->compute_nlstress (lcid,i,0,0);
	  break;
	default:
	  fprintf (stderr,"\n\n unknown element type is required in function nonloc_internal_forces (file %s, line %d).\n",
		   __FILE__,__LINE__);
	}
    }
  }
  
  for (i = 0; i < Mm->tnip; i++)
    Mm->nonlocaverage(i,0);
  Mp->phase=2;
  
  for (i = 0; i < ne; i++){
    if (Gtm->leso[i]==1){
      
      te = Mt->give_elem_type (i);
      ndofe=Mt->give_ndofe (i);
      allocv (ndofe,ifor);
      allocv (ndofe,cn);
      switch (te)
	{
	case planeelementlt :
	  Pelt->res_nonloc_internal_forces (lcid,i,ifor);
	  break;
	case lineartet:
	  Ltet->nonloc_internal_forces (lcid, i, 0, 0, ifor);
	  break;
	  
	case linearhex:
	  Lhex->nonloc_internal_forces (lcid, i, 0, 0, ifor);
	  break;
	  
	default:
	  fprintf (stderr,"\n\n unknown element type is required in function nonloc_internal_forces (file %s, line %d).\n",
		   __FILE__,__LINE__);
	}
      Mt->give_code_numbers (i,cn.a);
      locglob (intfor,ifor.a,cn.a,ndofe);
      destrv (ifor);
      destrv (cn);
    }
  }
}

/**
   function computes contributions to internal forces from one finite element
   
   @param i - element id
   @param lcid - load case id
   @param ifor - %vector of internal forces on one element
   
   JK, 3.11.2006
*/
void elem_internal_forces (long i,long lcid,vector &ifor)
{
  elemtype te;
  te = Mt->give_elem_type (i);
  
  switch (te){
    
  case bar2d:{
    Bar2d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case bar3d:{
    Bar3d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case barq2d:{
    Barq2d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case barq3d:{
    Barq3d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case beam2d:{
    Beam2d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case beam3d:{
    Beam3d->res_internal_forces (lcid,i,ifor);
    break;
  }
  case subsoilbeam:{
    Sbeam->res_internal_forces (lcid,i,ifor);
    break;
  }
  case beam2dsp:{
    Spbeam2d->res_internal_forces (lcid,i,ifor);
    break;
  }
    
    
  case spring_1:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
  case spring_2:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
  case spring_3:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
  case spring_4:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
  case spring_5:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
  case spring_6:{
    Spring->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case planeelementlt:{
    Pelt->res_internal_forces (lcid,i,ifor);
    break;
  }
  case planeelementqt:{
    Peqt->res_internal_forces (lcid,i,ifor);
    break;
  }
  case planeelementrotlt:{
    Perlt->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case planeelementlq:{
    Pelq->res_internal_forces (lcid,i,ifor);
    break;
  }
  case planeelementqq:{
    Peqq->res_internal_forces (lcid,i,ifor);
    break;
  }
  case planeelementrotlq:{
    Perlq->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case planeelementsubqt:{
    Pesqt->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case planequadcontact:{
    Pqcon->res_internal_forces (lcid,i,ifor);
    break;
  }
    
    
  case cctel:{
    Cct->internal_forces (lcid,i,ifor);
    break;
  }
    
  case q4plateel:{
    Q4pl->res_internal_forces (lcid,i,ifor);
    break;
  }
  case subsoilplateq:{
    Splq->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case axisymmlt:{
    Asymlt->res_internal_forces (lcid,i,ifor);
    break;
  }
  case axisymmlq:{
    Asymlq->res_internal_forces (lcid,i,ifor);
    break;
  }
  case axisymmqq:{
    Asymqq->res_internal_forces (lcid,i,ifor);
    break;
  }
    
  case lineartet:{
    Ltet->res_internal_forces (lcid,i,ifor);
    break;
  }