mechbclc.cpp

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

#include <stdlib.h>
#include <string.h>

#include "mechbclc.h"
#include "global.h"
#include "alias.h"
#include "loadcase.h"
#include "dloadcase.h"
#include "inicd.h"
#include "vector.h"
#include "matrix.h"
#include "elemhead.h"
#include "gtopology.h"
#include "element.h"
#include "elemswitch.h"



mechbclc::mechbclc ()
{
  nlc=0; nico = 0L;
  lc=NULL;
  dlc = NULL;
  ico = NULL;
  ncsum = 0;
  sumcomp = NULL;
}

mechbclc::~mechbclc ()
{
  delete [] lc;
  delete [] dlc;
  delete [] ico;
  delete [] sumcomp;
}

/**
   function reads data about boundary conditions and load cases
   
   @param in - input stream
   
   JK
*/
void mechbclc::read (XFILE *in)
{
  long i,j,k;
  long ncomp,nne,nip,ipid;
  strastrestate strastre;
  elemtype te;
  double *eigstr;
  vector nodval,ipval;
  
  switch (Mp->tprob){
  case linear_statics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc<0){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
    }
    break;
  }
  case mat_nonlinear_statics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [2*nlc];
    for (i=0;i<2*nlc;i++){
      lc[i].read (in);
    }
    readinic(in);
    break;
  }
  case geom_nonlinear_statics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [2*nlc];
    for (i=0;i<2*nlc;i++){
      lc[i].read (in);
    }
    readinic(in);
    break;
  }
  case eigen_dynamics:{
    if (Mp->straincomp==1 || Mp->stresscomp==1)
      nlc=Mp->eigsol.neigv;
    break;
  }
  case forced_dynamics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc<0){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    dlc = new dloadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
      dlc[i].read (in);
    }
    break;
  }
  case mech_timedependent_prob:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    dlc = new dloadcase [nlc];
    for (i=0;i<nlc;i++){
      dlc[i].read (in);
    }
    break;
  }
  case growing_mech_structure:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    dlc = new dloadcase [nlc];
    for (i=0;i<nlc;i++){
      dlc[i].read (in);
    }
    break;
  }
  case earth_pressure:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (Mp->tepsol < epplast_sol)
    {
      lc = new loadcase [nlc];
      dlc = new dloadcase [nlc];
      for (i=0;i<nlc;i++){
        lc[i].read (in);
        dlc[i].read (in);
      }
      readinic(in);
    }
    else
    {
      lc = new loadcase [2*nlc];
      for (i=0;i<2*nlc;i++){
        lc[i].read (in);
      }
      readinic(in);
    }
    break;
  }
  case layered_linear_statics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc<0){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
    }
    break;
  }
    
  case lin_floating_subdomain:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc<0){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
    }
    break;
  }

  case nonlin_floating_subdomain:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [2*nlc];
    for (i=0;i<2*nlc;i++){
      lc[i].read (in);
    }
    break;
  }

  case var_stiff_method:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc<0){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
    }
    break;
  }
    
    
  case nonlinear_dynamics:{
    xfscanf (in,"%k%ld","number_of_load_cases",&nlc);
    if (nlc!=1){
      fprintf (stderr,"\n\n wrong number of load cases in function mechbclc::read (file %s, line %d).\n",__FILE__,__LINE__);
      abort ();
    }
    lc = new loadcase [nlc];
    dlc = new dloadcase [nlc];
    for (i=0;i<nlc;i++){
      lc[i].read (in);
      dlc[i].read (in);
    }
    break;
  }
  default:{
    fprintf(stderr, "\n\n unknown problem type in function mechbclc::read (file %s, line %d)",__FILE__,__LINE__);
  }
  }
  
  // ********************
  //  eigstrain reading
  // ********************
  xfscanf (in,"%k%ld","eigenstrains",&(Mp->eigstrains));
  if (Mp->eigstrains==0 || Mp->eigstrains==1){}
  else{
    fprintf (stderr,"\n\n wrong definition of eigenstrains in function read (file %s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }
  if (Mespr == 1){
    if (Mp->eigstrains==0)
      fprintf(stdout, "\n eigenstrains are not defined");
    if (Mp->eigstrains==1)
      fprintf(stdout, "\n eigenstrains are defined");
  }
  
  if (Mp->eigstrains==1){
    //  type of strain/stress state
    xfscanf (in,"%d",(int *)&strastre);
    
    switch (strastre){
    case planestress:{
      eigstr = new double [3];
      xfscanf (in,"%lf %lf %lf",eigstr+0,eigstr+1,eigstr+2);
      ncomp=3;
      break;
    }
    case spacestress:{
      eigstr = new double [6];
      xfscanf (in,"%lf %lf %lf %lf %lf %lf",eigstr+0,eigstr+1,eigstr+2,eigstr+3,eigstr+4,eigstr+5);
      ncomp=6;
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown type of strain-stress state in function read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }
    

    if (Mm->eigstrains==NULL){
      Mm->eigstrains = new double* [Mm->tnip];
      for (i=0;i<Mm->tnip;i++){
	Mm->eigstrains[i] = new double [ncomp];
      }
    }
    
    for (i=0;i<Mt->ne;i++){
      if (Gtm->leso[i]==1){
	
	te = Mt->give_elem_type (i);
	//  number of nodes on element
	nne = Mt->give_nne (i);
	//  number of integration points on element
	nip = Mt->give_tnip (i);
	
	allocv (nne,nodval);
	allocv (nip,ipval);
	
	for (j=0;j<ncomp;j++){
	  //  loop over number of eigenstrain components
	  
	  //  selection of appropriate eigstrain component
	  for (k=0;k<nne;k++){
	    nodval[k]=eigstr[j];
	  }
	  
	  //  interpolation of nodal values to integration points
	  elem_intpointval (i,nodval,ipval);
	  
	  //  number of the first integration point
	  ipid=Mt->elements[i].ipp[0][0];
	  
	  for (k=0;k<nip;k++){
	    Mm->eigstrains[ipid][j]=ipval[k];
	    
	    //assign eigenstrains only to empty porosity (26) or water-filled porosity (phase 28), Vit Smilauer, 10.1.2006
	    //check the material phase, in idm[0] - the phase is one less
	    //if(Mp->homog==2 && (Mt->elements[i].idm[0]!=(26-1) && Mt->elements[i].idm[0]!=(28-1))){
	    if(Mp->homog==2 && Mt->elements[i].idm[0]!=(28-1)){
	      Mm->eigstrains[ipid][j]=0.;//zero
	    }
	    ipid++;
	  }
	}
	
	destrv (nodval);
	destrv (ipval);
      }
    }
  }
  
}

/**
  This method reads section with data about initial conditions in the nodes
  and if it is necessarry it allocates memory for mechmat ic array.
  
  @param in is pointer to the structure with opened file.

  @retval 0 - on succes
  @retval 1 - error reading number of nodes with initial condition
  @retval 2 - error reading node number
  @retval 3 - error reading type or values of given condition
*/
long mechbclc::readinic (XFILE *in)
{
  long n = xfscanf(in, "%ld", &nico);
  if (n != 1)
    return 0;
  if ((nico < 0) && (nico > Mt->nn))
  {
    fprintf (stderr,"\n\n Error reading number of nodes with inital condition\n in function mechbclc::readinic (file %s, line %d).\n",__FILE__,__LINE__);
    return 1;
  }
  if (nico == 0)
    return 0;
  ico = new inicd[Mt->nn];
  for (long i = 0L; i < nico; i++)
  {
    if (xfscanf(in, "%ld", &n) != 1)
    {
      fprintf (stderr,"\n\n Error reading node number\n in function mechbclc::readinic (file %s, line %d).\n",__FILE__,__LINE__);
      return 2;
    }
    if ((n < 0) || (n > Mt->nn))
    {
      fprintf (stderr,"\n\n Error reading node number\n in function mechbclc::readinic (file %s, line %d).\n",__FILE__,__LINE__);
      return 2;
    }
    if (ico[n-1].read(in))
    {
      fprintf (stderr,"\n\n Error reading values of initial condition of node %ld\n in function mechbclc::readinic (file %s, line %d).\n",n,__FILE__,__LINE__);
      return 3;
    }
    if ((ico[n-1].type & inicond) && (Mm->ic == NULL))
    {
      Mm->ic = new double *[Mm->tnip];
      memset(Mm->ic, 0, sizeof(*Mm->ic)*Mm->tnip);
    }
  }

  return 0;
}

/**
   function computes initial values at integration points
   from initial nodal values
   
   TKo
*/
void mechbclc::inicipval(void)
{
  long i, j, k, l, nne, nval, eid;
  elemtype et;
  ivector enod;
  matrix  nodval;
  inictype *ictn;   // type of initial condition in the nodes


  for (i = 0; i < Mt->ne; i++){
    if (Gtm->leso[i]==1){
      
      nne = Mt->give_nne(i);
      allocv(nne, enod);
      Mt->give_elemnodes (i, enod);
      nval = ico[enod[0]].nval;
      for (j = 0; j < nne; j++)
	{
	  if (nval < ico[enod[0]].nval)
	    nval = ico[enod[0]].nval;
	}
      allocm (nne, nval, nodval);
      fillm(0.0, nodval);
      ictn = new inictype[nne];
      for (j = 0; j < nne; j++)
	{
	  ictn[j] = ico[enod[j]].type;
	  if (ictn[j] & inidisp)  // skip initial nodal displacements
	    k = Gtm->give_ndofn(enod[j]);
	  else
	    k = 0;
	  for (k = 0; k < nval; k++)
	    {
	      if (ictn[j] & inidisp)  // skip initial nodal displacements
		{
		  l = k + Gtm->give_ndofn(enod[j]);
		  if (k < Gtm->give_ndofn(enod[j]))
		    continue;
		}
	      else
		l = k;
	      nodval[j][k] = ico[enod[j]].val[l];
	    }
	}
      eid = i;
      et = Mt->give_elem_type(eid);
      switch (et)
	{
	case bar2d:{              Bar2d->inicipval(eid, 0, 0, nodval, ictn);   break; }
	case barq2d:{             Barq2d->inicipval(eid, 0, 0, nodval, ictn);  break; }
	  //      case beam2d:{             Beam2d->inicipval(eid, 0, 0, nodval, ictn);  break; }
	  //      case beam3d:{             Beam3d->inicipval(eid, 0, 0, nodval, ictn);  break; }
	  //      case subsoilbeam:{        Sbeam->inicipval(eid, 0, 0, nodval, ictn);   break; }
	case spring_1:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case spring_2:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case spring_3:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case spring_4:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case spring_5:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case spring_6:{           Spring->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case planeelementlt:{     Pelt->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case planeelementqt:{     Peqt->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case planeelementrotlt:{  Perlt->inicipval(eid, 0, 0, nodval, ictn);   break; }
	case planeelementlq:{     Pelq->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case planeelementqq:{     Peqq->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case planeelementrotlq:{  Perlq->inicipval(eid, 0, 0, nodval, ictn);   break; }
	case planeelementsubqt:{  Pesqt->inicipval(eid, 0, 0, nodval, ictn);   break; }
	case cctel:{              Cct->inicipval(eid, 0, 0, nodval, ictn);     break; }
	case dktel:{              Dkt->inicipval(eid, 0, 0, nodval, ictn);     break; }
	case dstel:{              Dst->inicipval(eid, 0, 0, nodval, ictn);     break; }
	case q4plateel:{          Q4pl->inicipval(eid, 0, 0, nodval, ictn);    break; }
	  //      case subsoilplatetr:{     Spltr->inicipval(eid, 0, 0, nodval, ictn);   break; }
	  //      case subsoilplateq:{      Splq->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case axisymmlq:{          Asymlq->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case axisymmlt:{          Asymlt->inicipval(eid, 0, 0, nodval, ictn);  break; }
	case shelltrelem:{        Shtr->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case shellqelem:{         Shq->inicipval(eid, 0, 0, nodval, ictn);     break; }
	case lineartet:{          Ltet->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case linearhex:{          Lhex->inicipval(eid, 0, 0, nodval, ictn);    break; }
	case quadrhex:{           Qhex->inicipval(eid, 0, 0, nodval, ictn);    break; }
	default:
	  {
	    fprintf (stderr,"\n\n unknown element type is required in function");
	    fprintf (stderr,"\n mechbclc::inicipval (file %s, line %d).\n",__FILE__,__LINE__);
	  }
	}
      destrv(enod);
      destrm(nodval);
      delete [] ictn;
    }
  }
}

/**
   function allocates array containing sums of components in particular directions of the nodal force vector
   
   6.10.2003
*/
void mechbclc::alloc_sumcomp ()
{
  ncsum = Mt->give_ndofn (0);
  sumcomp = new double [ncsum];
}

/**
   function computes sums of components in particular directions of the nodal force vector
   
   @param rhs - %vector of right hand side (nodal forces)
   
   6.10.2003, JK
*/
void mechbclc::comp_sum (double *rhs)
{
  long i,j,ndofn,cn;
  
  nullv (sumcomp,ncsum);
  
  for (i=0;i<Mt->nn;i++){
    ndofn=Gtm->gnodes[i].ndofn;
    for (j=0;j<ndofn;j++){
      cn=Mt->give_dof (i,j);
      if (cn>0){
        sumcomp[j]+=rhs[cn-1];
      }
    }
  }
  
}

/**
   function returns sums of components in particular directions of the nodal force vector
   
   6.10.2003, JK
*/
void mechbclc::give_comp_sum (double *sum)
{
  long i;
  for (i=0;i<ncsum;i++){
    sum[i]=sumcomp[i];
  }
}