globmat.cpp

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

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

  switch (te){

  case bar2d:{
    Bar2d->mass_matrix (eid,mm);
    break;
  }
    /*
  case barq2d:{
    Barq2d->res_mass_matrix (eid,mm);
    break;
  }
    */
  case beam2d:{
    Beam2d->res_mass_matrix (eid,mm);
    break;
  }
  case beam3d:{
    Beam3d->res_mass_matrix (eid,mm);
    break;
  }

  
  case planeelementlt:{
    Pelt->res_mass_matrix (eid,mm);
    break;
  }
  case planeelementqt:{
    Peqt->mass_matrix (eid,mm);
    break;
  }
  case planeelementrotlt:{
    Perlt->res_mass_matrix (eid,mm);
    break;
  }
  case planeelementlq:{
    Pelq->res_mass_matrix (eid,mm);
    break;
  }
  case planeelementqq:{
    Peqq->res_mass_matrix (eid,mm);
    break;
  }
  case planeelementrotlq:{
    Perlq->res_mass_matrix (eid,mm);
    break;
  }

  case planeelementsubqt:{
    Pesqt->res_mass_matrix (eid,mm);
    break;
  }

  case cctel:{
    Cct->res_mass_matrix (eid,mm);
    break;
  }

  case axisymmlt:{
    Asymlt->mass_matrix (eid,mm);
    break;
  }
  case axisymmlq:{
    Asymlq->mass_matrix (eid,mm);
    break;
  }
  case axisymmqq:{
    Asymqq->mass_matrix (eid,mm);
    break;
  }

  case lineartet:{
    Ltet->mass_matrix (eid,mm);
    break;
  }
  case linearhex:{
    Lhex->mass_matrix (eid,mm);
    break;
  }
  case quadrhex:{
    Qhex->mass_matrix (eid,mm);
    break;
  }

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

void initstiffmat (long lcid,long eid,matrix &sm)
{
  elemtype te;
  
  te = Mt->give_elem_type (eid);

  switch (te){

  case beam2d:{
    Beam2d->initstr_matrix_expl (eid,0,0,sm);
    break;
  }

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

}


/**
   function assembles constraint matrix of one layer on one node
   
   @param nid - node id
   @param cid - contact id
   @param lcm - local constraint matrix
   
   7.12.2002
*/
void constr_matrix (long nid,long cid,matrix &lcm)
{
  long i,idcs;
  crsectype crst;
  double tl,tu;
  
  i=Gtm->lgnodes[nid].nodes[cid-1];
  crst = Mt->nodes[i].crst;
  idcs = Mt->nodes[i].idcs;
  tl = Mc->give_onethickness (crst,idcs);
  i=Gtm->lgnodes[nid].nodes[cid];
  crst = Mt->nodes[i].crst;
  idcs = Mt->nodes[i].idcs;
  tu = Mc->give_onethickness (crst,idcs);
  
  fillm (0.0,lcm);

  if (Mp->tlm==1){
    lcm[0][0] = -1.0;
    lcm[1][1] = -1.0;
    lcm[2][0] = tl/(-2.0);
    
    lcm[3][0] = 1.0;
    lcm[4][1] = 1.0;
    lcm[5][0] = tu/(-2.0);
  }
  
  if (Mp->tlm==2){
    
    /*
      uplne stara verze
      lcm[0][0] = -1.0;
      lcm[1][1] = -1.0;
      lcm[2][2] = -1.0;
      lcm[1][3] = tl/2.0;
      lcm[0][4] = tl/(-2.0);
      lcm[3][5] = -1.0;
      
      lcm[0][6]  = 1.0;
      lcm[1][7]  = 1.0;
      lcm[2][8]  = 1.0;
      lcm[1][9]  = tu/2.0;
      lcm[0][10] = tu/(-2.0);
      lcm[3][11] = 1.0;
    */
    
    
    lcm[0][0] = -1.0;
    lcm[1][1] = -1.0;
    lcm[2][2] = -1.0;
    lcm[3][1] = tl/2.0;
    lcm[4][0] = tl/(-2.0);
    lcm[5][3] = -1.0;
    
    lcm[6][0]  = 1.0;
    lcm[7][1]  = 1.0;
    lcm[8][2]  = 1.0;
    lcm[9][1]  = tu/2.0;
    lcm[10][0] = tu/(-2.0);
    lcm[11][3] = 1.0;
  }
  
}

void mefel_right_hand_side (long lcid,double *rhs)
{
  long ne;
  ivector cn;
  vector nfor;
  
  ne=Mt->ne;

  switch (Mp->tprob){
  case linear_statics:{
    nullv (rhs+lcid*Ndofm,Ndofm);
    Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,1.0);
    break;
  }
  case eigen_dynamics:{
    break;
  }
  case forced_dynamics:{
    nullv (rhs+lcid*Ndofm,Ndofm);
    //Mb->lc[lcid].assemble (lcid,av);
    Mb->dlc[lcid].assemble (lcid,rhs,Ndofm,Mp->time);
    break;
  }
  case linear_stability:{

    break;
  }
  case mat_nonlinear_statics:{
    nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
    nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
    Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,1.0);
    Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,1.0);
    break;
  }
  case geom_nonlinear_statics:{
    nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
    nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
    Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,1.0);
    Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,1.0);
    break;
  }
  case mech_timedependent_prob:{
    nullv (rhs+lcid*Ndofm, Ndofm);
    Mb->dlc[lcid].assemble (lcid,rhs+lcid*Ndofm,Ndofm,Mp->time);
    break;
  }
  case growing_mech_structure:{
    //  lcid must be equal to zero
    nullv (rhs, Ndofm);
    Mb->dlc[lcid].assemble (0,rhs,Ndofm,Mp->time);
    break;
  }
  case earth_pressure:{

    break;
  }
  case layered_linear_statics:{

    break;
  }
  case lin_floating_subdomain:{
    nullv (rhs+lcid*Ndofm,Ndofm);
    Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,1.0);
    break;
  }
  case nonlin_floating_subdomain:{
    nullv (rhs+(lcid*2+0)*Ndofm, Ndofm);
    nullv (rhs+(lcid*2+1)*Ndofm, Ndofm);
    Mb->lc[lcid*2+0].assemble (lcid*2+0,rhs+(lcid*2+0)*Ndofm,1.0);
    Mb->lc[lcid*2+1].assemble (lcid*2+1,rhs+(lcid*2+1)*Ndofm,1.0);
    break;
  }

  case var_stiff_method:{
    nullv (rhs+lcid*Ndofm,Ndofm);
    Mb->lc[lcid].assemble (lcid,rhs+lcid*Ndofm,1.0);
    break;
  }

  default:{
    fprintf (stderr,"\n\n unknown problem type is required in function mefel_right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
  
  
  if (Mp->eigstrains==1){

    //if (Mb->
    Mm->est=eigstrain;
    nodal_eigstrain_forces (lcid,rhs);
    
  /*
    for (i=0;i<ne;i++){
      te = Mt->give_elem_type (i);
      ndofe=Mt->give_ndofe (i);
      allocv (ndofe,nfor);
      allocv (ndofe,cn);
      
      switch (te){
	
      case barq2d:{
	break;
      }
      case planeelementlq:{
	Pelq->res_eigstrain_rhs (i,nfor);
	break;
      }
	
      case lineartet:{
	Ltet->res_eigstrain_rhs (i,nfor);
	break;
      }
      case linearhex:{
	Lhex->res_eigstrain_rhs (i,nfor);
	break;
      }
      case quadrhex:{
	Qhex->res_eigstrain_rhs (i,nfor);
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function mefel_right_hand_side (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
      Mt->give_code_numbers (i,cn.a);
      locglob (rhs,nfor.a,cn.a,ndofe);
      
      destrv (nfor);
      destrv (cn);
    }
  */

  }
  
}

/**
   function computes all required values
   this function is called e.g. before output print

   @param lcid - load case id

   1.4.2004, JK
*/
void compute_req_val (long lcid)
{
  
  /*
  //  strains computation
  if (Mp->straincomp==1){
    if (Mp->tprob != nonlinear_statics)
      Mm->computestrains (lcid);
    else
    {
      if (Mp->strainaver==1)
        Mm->compute_nodestrains (lcid);
    }
    Mm->stra.transformvalues (1);
  }
  if (Mp->strainaver==1)
    Mm->compute_nodestrains (lcid);
  
  //Mm->compute_averstrains ();
  
  //  stresses computation
  if (Mp->stresscomp==1){
    if (Mp->tprob != nonlinear_statics)
      Mm->computestresses (lcid);
    else
    {
      if (Mp->strainaver==1)
        Mm->compute_nodestresses (lcid);
    }
    Mm->stre.transformvalues(0);
  }

  //Mm->compute_averstrains ();

  if (Mp->otheraver)
    Mm->compute_nodeothers (lcid);
  //  reactions computation
  if (Mp->reactcomp==1){
    switch (Mp->tprob){
    case linear_statics:
    case nonlinear_statics:
    case earth_pressure:
      Mb->lc[lcid].compute_reactions (lcid);
      break;
    case forced_dynamics:
    case mech_timedependent_prob:
      break;
    default:{
      fprintf (stderr,"\n\n unknown type of problem is required in function compute_req_val (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
  }
  */


  //  
  Mm->totstrains ();

  //Mm->computestresses (lcid);
  if (Mp->straincomp == 1){
    if (Mp->strainstate==0)
      compute_ipstrains (lcid);
    if (Mp->strainpos == 2)
      compute_nodestrains (lcid);
  }
  
  if (Mp->stresscomp == 1){
    if (Mp->stressstate==0)
      compute_ipstresses (lcid);
    if (Mp->stresspos == 2)
      compute_nodestresses (lcid);
  }
  
  //if (Mp->othercomp == 1)
  //compute_nodeothers (lcid);
  
  if (Mp->reactcomp==1){
    //Mb->lc[lcid].compute_reactions (lcid);
  }

  
}





/**
   Function initializes material models with initial values.
   
   @param gv - array containing values at all nodes of the mesh
   @param nmq - type of non-mechanical quantity
   
   21.6.2004, JK
*/
void initmaterialmodels (void)
{
  long i,j,nip,ipp,nm,ido;
  elemtype te;
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      te = Mt->give_elem_type (i);
      nip = Mt->give_tnip (i);
      //  number of the first integration point
      ipp=Mt->elements[i].ipp[0][0];
      for(j=0; j<nip; j++){
	nm = Mm->ip[ipp+j].nm;
	Mm->initvalues(ipp+j, 0, 0);
	if (Mm->ip[ipp].hmt & 1){
	  // thermal dilatancy is present and ncompeqother 
	  switch (Mm->ip[ipp].tm[nm-1]){
	  case therisodilat:
	    break;
	  case therisodilattime:    
	    ido = Mm->givencompeqother(ipp+j, 0);
	    ido -= Mm->givencompeqother(ipp+j, nm-1);
	    Mm->initvalues(ipp+j, nm-1, ido);
	    break;
	  default:{
	    fprintf (stderr,"\n\n Unknown thermal material type is required in function");
	    fprintf (stderr,"\n givencompeqother (file %s, line %d).\n",__FILE__,__LINE__);
	  }
	  }
	}
      }
    }
  }
}


void local_global_displ_transf (long lcid)
{
  long i,j,k,ndofn;
  long *cn;
  vector l,g;
  matrix tm(3,3);
  
  for (i=0;i<Mt->nn;i++){
    if (Mt->nodes[i].transf>0){

      //fprintf (Outm,"\n uzel %ld   %ld",i,Mt->nodes[i].transf);

      ndofn = Mt->give_ndofn (i);
      allocv (ndofn,l);
      allocv (ndofn,g);
      
      noddispl (lcid,l.a,i);
      
      tm[0][0]=Mt->nodes[i].e1[0];
      tm[1][0]=Mt->nodes[i].e1[1];
      tm[2][0]=Mt->nodes[i].e1[2];
      
      tm[0][1]=Mt->nodes[i].e2[0];
      tm[1][1]=Mt->nodes[i].e2[1];
      tm[2][1]=Mt->nodes[i].e2[2];
      
      tm[0][2]=Mt->nodes[i].e3[0];
      tm[1][2]=Mt->nodes[i].e3[1];
      tm[2][2]=Mt->nodes[i].e3[2];
      
      mxv(tm,l,g);
      
      cn = new long [ndofn];
      Mt->give_node_code_numbers (i,cn);
      
      for (j=0;j<ndofn;j++){
	k=cn[j]-1;
	if (k>-1)
	  Lsrs->lhs[k]=g[j];
      }
      
      delete [] cn;
      destrv (l);
      destrv (g);
    }
  }

}

void rhs_graphic (double *rhs,long lcid,long n)
{
  long i,j,k,m,ndofe,nne,ndofn;
  long *cn;
  ivector nodes;
  vector r,f;
  matrix sm;
  
  nullv (rhs,n);
  
  for (i=0;i<Mt->ne;i++){
    if (Gtm->leso[i]==1){
      
      ndofe=Mt->give_ndofe (i);
      nne=Mt->give_nne (i);
      allocv (nne,nodes);
      Mt->give_elemnodes (i,nodes);
      
      allocm (ndofe,ndofe,sm);
      stiffmat (lcid,i,sm);
      
      allocv (ndofe,r);
      allocv (ndofe,f);
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);
      for (j=0;j<ndofe;j++){
	r[j]=Mt->elements[i].initdispl[j];
      }
      
      m=0;
      for (j=0;j<nne;j++){
	ndofn=Mt->give_ndofn (nodes[j]);
	for (k=0;k<ndofn;k++){
	  r[m]+=Mt->nodes[nodes[j]].nvgraph[k];
	  m++;
	}
      }
      
      mxv (sm,r,f);
      cmulv (-1.0,f);
      locglob (rhs,f.a,cn,ndofe);
      
      delete [] cn;
      destrm (sm);
      destrv (f);  destrv (r);
      destrv (nodes);
    }
  }
}