globmat.cpp

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

  case linearhex:{
    Lhex->res_internal_forces (lcid,i,ifor);
    break;
  }
  case quadrhex:{
    Qhex->res_internal_forces (lcid,i,ifor);
    break;
  }
    
    
  case particleelem:{
    Pelem->res_internal_forces (i,ifor);
    break;
  }
    
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function elem_internal_forces (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
}


/**
   function adds contributions to internal forces caused by
   Lagrange multipliers
   function is used in problems with floating subdomains
   solved by the FETI method
   
   @param lcid - load case id
   @param intfor - array of internal forces
   
   JK, 22.6.2006
*/
void lagrmultcontr_intforces (long lcid,double *intfor)
{
  long i;
  double *av,*av1,*av2,*cv;
  
  av = new double [Ndofm];
  av1 = new double [Ndofm];
  av2 = new double [Ndofm];
  cv = new double [Fsd->nlm];
  
  //  E^T lamba
  nullv (av,Ndofm);
  nullv (av1,Ndofm);
  nullv (av2,Ndofm);
  //Gtm->flsub.coarse_local (Fsd->tw,av);
  Gtm->flsub.coarse_local (Fsd->muo,av1);
  Gtm->flsub.coarse_local (Fsd->ddmu,av2);
  
  for (i=0;i<Ndofm;i++){
    av[i]=av1[i]+av2[i];
  }

  //Gtm->flsub.local_coarse (cv,Lsrs->lhs);
  
  /*
  long i;
  for (i=0;i<Ndofm;i++){
    fprintf (Out,"\n lagr mult %le",av[i]);
  }
  for (i=0;i<Fsd->nlm;i++){
    fprintf (Out,"\n coarse vector %le",cv[i]);
    

    if (cv[i]>10.0){
      if (av[i]<0.0){
	av[i]=-10.0;
	Fsd->tw[i]=-10.0;
      }
      if (av[i]>0.0){
	av[i]=10.0;
	Fsd->tw[i]=10.0;
      }
      Fsd->compli[i]=2.0;
    }
    
  }
  */

  addv (intfor,av,Ndofm);
  
  delete [] av;
  delete [] cv;
}


/**
   function computes vector of nodal forces
   caused by eigenstrains
   
   @param lcid - load case id
   @param nodfor - %vector of nodal forces
   
   28.7.2001, JK
*/
void nodal_eigstrain_forces (long lcid,double *nodfor)
{
  long i,ne,ndofe;
  elemtype te;
  ivector cn;
  vector nfor;
  
  ne=Mt->ne;
  
  //  nullvr (nodfor,Ndofm);
  
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      
      te = Mt->give_elem_type (i);
      ndofe=Mt->give_ndofe (i);
      allocv (ndofe,nfor);
      allocv (ndofe,cn);
      
      elem_eigstrain_forces (lcid,i,nfor);
      
      Mt->give_code_numbers (i,cn.a);
      locglob (nodfor,nfor.a,cn.a,ndofe);
      destrv (nfor);
      destrv (cn);
    }
  }
  
}

/**
   function computes nodal forces caused by eigenstrains
   
   @param lcid - 
   @param eid - element id
   @param nfor - %vector of nodal values
   
   JK
*/
void elem_eigstrain_forces (long lcid,long eid,vector &nfor)
{
  switch (Mt->give_elem_type(eid)){
    
  case bar2d:{
    Bar2d->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case bar3d:{
    Bar3d->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case barq2d:{
    Barq2d->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case barq3d:{
    Barq3d->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case planeelementlt:{
    Pelt->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case planeelementqt:{
    Peqt->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case planeelementlq:{
    Pelq->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case planeelementqq:{
    Peqq->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case axisymmlt:{
    Asymlt->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case axisymmlq:{
    Asymlq->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
  case axisymmqq:{
    Asymqq->res_eigstrain_forces (lcid,eid,nfor);
    break;
  }
    
  case linearhex:{
    Lhex->res_eigstrain_forces (lcid,eid,nfor);
    break;
 }

  default:{
    fprintf (stderr,"\n\n unknown element type is required in function res_temp_forces (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
}




/**
   function extracts displacements on one element
   
   @param lcid - number of load case
   @param cn - array containing code numbers
   @param ndofe - number of DOFs on actual element
   
   9.7.2001
*/
void eldispl (long lcid,long eid,double *r,long *cn,long ndofe)
{
  long i,ii;
  elemtype te;
  
  switch (Mp->tprob){
  case forced_dynamics:
  case mech_timedependent_prob:{
    for (i=0;i<ndofe;i++){
      ii=cn[i];
      if (ii<0)   r[i]=Mb->dlc[lcid].pd[0-ii-1].getval(Mp->time);
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
    break;
  }

  case mat_nonlinear_statics:{
    for (i=0;i<ndofe;i++){
      ii=cn[i];
      if (ii<0)   r[i]=Mb->lc[lcid].pd[0-ii-1] * Mp->lambda + Mb->lc[lcid+1].pd[0-ii-1];
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
    break;
  }
    
  case growing_mech_structure:{
    //  lcid must be equal to zero
    for (i=0;i<ndofe;i++){
      ii=cn[i];
      
      //fprintf (Out,"eid = %ld,cn = %ld,e\n",eid,ii);
      
      if (ii<0)   r[i]=Mb->lc[0].pd[0-ii-1];
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[ii-1];
    }
    
    //jeste to promyslet??!!
    //pridano TKr 22.7.2006
    te = Mt->give_elem_type (eid);
    
    switch (te){
      
    case spring_1:
    case spring_2:
    case spring_3:
    case spring_4:
    case spring_5:
    case spring_6:{
      //  subtraction of initial displacements
      Mt->elements[eid].subtrinitdispl (r,ndofe);
      break;
    }
    default:{
      //  subtraction of initial displacements
      Mt->elements[eid].subtrinitdispl (r,ndofe);
    }
    }
    
    break;
  }
    
  default:{
    for (i=0;i<ndofe;i++){
      ii=cn[i];
      if (ii<0)   r[i]=Mb->lc[lcid].pd[0-ii-1];
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
  }
  }
}

/**
   function extracts displacements on one node
   
   @param lcid  - number of load case
   @param r     - allocated array for displacement
   @param nid   - node number (node id)

   9.7.2001
*/
void noddispl (long lcid,double *r, long nid)
{
  long i,ii;
  
  switch (Mp->tprob){
  case forced_dynamics:
  case growing_mech_structure://tady??!!
  case mech_timedependent_prob:{
    for (i=0;i<Mt->give_ndofn (nid);i++){
      ii=Mt->give_dof(nid,i);
      
      //fprintf (Out,"nid = %ld,cn = %ld,n\n",nid,ii-1);
      
      if (ii<0)   r[i]=Mb->dlc[lcid].pd[0-ii-1].getval(Mp->time);
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
    break;
  }
  case mat_nonlinear_statics:{
    for (i=0;i<Mt->give_ndofn (nid);i++){
      ii=Mt->give_dof(nid,i);
      if (ii<0)   r[i]=Mb->lc[lcid].pd[0-ii-1]*Mp->lambda + Mb->lc[lcid+1].pd[0-ii-1];
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
    break;
  }
  default:{
    for (i=0;i<Mt->give_ndofn (nid);i++){
      ii=Mt->give_dof(nid,i);
      if (ii<0)   r[i]=Mb->lc[lcid].pd[0-ii-1];
      if (ii==0)  r[i]=0.0;
      if (ii>0)   r[i]=Lsrs->lhs[lcid*Ndofm+ii-1];
    }
  }
  }
}

/**
   function extracts displacements on one node for given problem
   
   @param lcid - number of load case
   @param idn  - node number
   @param p    - pointr on problem
   @param r    - answer - allocated array for displacement
   
   11.10.2003
*/
void noddispl (long lcid, long idn, problem *p, double *r)
{
  long i,ii;
  
  switch (p->mp->tprob)
  {
  case forced_dynamics :
  case growing_mech_structure://tady??!!
  case mech_timedependent_prob :
    {
      for (i=0;i<p->gt->give_ndofn (idn);i++){
        ii = p->gt->give_dof(idn,i);
        if (ii<0)   r[i]=p->mb->dlc[lcid].pd[0-ii-1].getval(p->mp->time);
        if (ii==0)  r[i]=0.0;
        if (ii>0)   r[i]=p->lsrs->lhs[lcid*p->lsrs->ndof+ii-1];
      }
      break;
    }
    case mat_nonlinear_statics :
    {
      for (i=0;i<p->gt->give_ndofn (idn);i++){
        ii = p->gt->give_dof(idn,i);
        if (ii<0)   r[i]=p->mb->lc[lcid].pd[0-ii-1] * p->mp->lambda + p->mb->lc[lcid+1].pd[0-ii-1];
        if (ii==0)  r[i]=0.0;
        if (ii>0)   r[i]=p->lsrs->lhs[lcid*p->lsrs->ndof+ii-1];
      }
      break;
    }
    default :
    {
      for (i=0;i<p->gt->give_ndofn (idn);i++)
      {
         ii = p->gt->give_dof(idn,i);
        if (ii<0)   r[i]=p->mb->lc[lcid].pd[0-ii-1];
        if (ii==0)  r[i]=0.0;
        if (ii>0)   r[i]=p->lsrs->lhs[lcid*p->lsrs->ndof+ii-1];
      }
    }
  }
}


/**
   function computes stiffness matrix of required element
   
   @param lcid - load case id
   @param eid - element id
   @param sm - stiffness matrix of the element
*/
void stiffmat (long lcid,long eid,matrix &sm)
{
  elemtype te;
  
  te = Mt->give_elem_type (eid);

  switch (te){

  case bar2d:{
    Bar2d->res_stiffness_matrix (eid,sm);
    break;
  }
  case bar3d:{
    Bar3d->res_stiffness_matrix (eid,sm);
    break;
  }
  case barq2d:{
    Barq2d->res_stiffness_matrix (eid,sm);
    break;
  }
  case barq3d:{
    Barq3d->res_stiffness_matrix (eid,sm);
    break;
  }
  case beam2d:{
    Beam2d->res_stiffness_matrix (eid,sm);
    break;
  }
  case beam3d:{
    Beam3d->res_stiffness_matrix (eid,sm);
    break;
  }
  case beamg3d:{
    Beam3dg->res_stiffness_matrix (eid,sm);
    break;
  }
  case subsoilbeam:{
    Sbeam->res_stiffness_matrix (eid,sm);
    break;
  }
  case beam2dsp:{
    Spbeam2d->res_stiffness_matrix (eid,sm);
    break;
  }

  case spring_1:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }
  case spring_2:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }
  case spring_3:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }
  case spring_4:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }
  case spring_5:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }
  case spring_6:{
    Spring->res_stiffness_matrix (eid,sm);
    break;
  }

  case planeelementlt:{
    Pelt->res_stiffness_matrix (eid,sm);
    break;
  }
  case planeelementqt:{
    Peqt->res_stiffness_matrix (eid,sm);
    break;
  }
  case planeelementrotlt:{
    Perlt->res_stiffness_matrix (eid,sm);
    break;
  }

  case planeelementlq:{
    Pelq->res_stiffness_matrix (lcid,eid,sm);
    break;
  }
  case planeelementqq:{
    Peqq->res_stiffness_matrix (eid,sm);
    break;
  }
  case planeelementrotlq:{
    Perlq->res_stiffness_matrix (eid,sm);
    break;
  }
    
  case planeelementsubqt:{
    Pesqt->res_stiffness_matrix (eid,sm);
    break;
  }

  case planequadcontact:{
    Pqcon->res_stiffness_matrix (lcid,eid,sm);
    break;
  }
    
    
  case cctel:{
    Cct->res_stiffness_matrix (eid,sm);
    break;
  }
  case dktel:{
    Dkt->res_stiffness_matrix (eid,sm);
    break;
  }
  case dstel:{
    Dst->res_stiffness_matrix (eid,sm);
    break;
  }
  case q4plateel:{
    Q4pl->res_stiffness_matrix (eid,sm);
    break;
  }
  case subsoilplatetr:{
    Spltr->res_stiffness_matrix (eid,sm);
    break;
  }
  case subsoilplateq:{
    Splq->res_stiffness_matrix (eid,sm);
    break;
  }
    
  case axisymmlt:{
    Asymlt->res_stiffness_matrix (eid,sm);
    break;
  }
  case axisymmlq:{
    Asymlq->res_stiffness_matrix (eid,sm);
    break;
  }
  case axisymmqq:{
    Asymqq->res_stiffness_matrix (eid,sm);
    break;
  }

  case shelltrelem:{
    Shtr->res_stiffness_matrix (eid,sm);
    break;
  }
  case shellqelem:{
    Shq->res_stiffness_matrix (eid,sm);
    break;
  }

  case lineartet:{
    Ltet->res_stiffness_matrix (eid,sm);
    break;
  }
  case linearhex:{
    Lhex->res_stiffness_matrix (lcid,eid,sm);
    break;
  }
  case quadrhex:{
    Qhex->res_stiffness_matrix (eid,sm);
    break;
  }
  case linearwed:{
    Lwed->res_stiffness_matrix (eid,sm);
    break;
  }
  case quadrwed:{
    Qwed->res_stiffness_matrix (eid,sm);
    break;
  }
    
  case particleelem:{
    Pelem->res_stiffness_matrix (eid,sm);
    break;
  }
  default:{
    fprintf (stderr,"\n unknown element type is required in function");
    fprintf (stderr,"\n stiffmat (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}