lhsrhs.cpp

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

#include <string.h>
#include "lhsrhs.h"
#include "global.h"
#include "mechprint.h"
#include "gtopology.h"
#include "dloadcase.h"



lhsrhs::lhsrhs ()
{
  ndof=Ndofm;  nlc=Mb->nlc;
  lhs=NULL;  tdlhs=NULL;  stdlhs=NULL;  rhs=NULL;
  w=NULL;
}

lhsrhs::~lhsrhs ()
{
  delete [] lhs;
  delete [] tdlhs;
  delete [] stdlhs;
  delete [] rhs;
  delete [] w;
}

/**
   function allocates arrays for left and right hand sides of problem
   
   JK, 29.6.2001
*/
void lhsrhs::alloc ()
{
  switch (Mp->tprob){
  case linear_statics:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*ndof];
    memset (rhs,0,(nlc*ndof)*sizeof(double));
    break;
  }
  case eigen_dynamics:{
    lhs = new double [Mp->eigsol.nev*ndof];
    memset (lhs,0,(Mp->eigsol.nev*ndof)*sizeof(double));
    
    switch (Mp->eigsol.teigsol){
    case inv_iteration:{
      w = new double [1];
      break;
    }
    case subspace_it_jacobi:{
      w = new double [Mp->eigsol.nev];
      break;
    }
    case subspace_it_gsortho:{
      w = new double [Mp->eigsol.nev];
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown eigenvalue solver is required in");
      fprintf (stderr,"\n function alloc (file %s, line %d).\n",__FILE__,__LINE__);
    }
    }

    break;
  }
  case forced_dynamics:{
    if (Mp->tforvib==modal_analysis){
      lhs = new double [ndof+Mp->eigsol.nev*ndof];
      memset (lhs,0,(ndof+Mp->eigsol.nev*ndof)*sizeof(double));

      switch (Mp->eigsol.teigsol){
      case inv_iteration:{
	w = new double [1];
	break;
      }
      case subspace_it_jacobi:{
	w = new double [Mp->eigsol.nev];
	break;
      }
      case subspace_it_gsortho:{
	w = new double [Mp->eigsol.nev];
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown eigenvalue solver is required in");
	fprintf (stderr,"\n function alloc (file %s, line %d).\n",__FILE__,__LINE__);
      }
      }
      
    }
    else{
      lhs = new double [nlc*ndof];
      memset (lhs,0,(nlc*ndof)*sizeof(double));
      tdlhs = new double [nlc*ndof];
      memset (tdlhs,0,(nlc*ndof)*sizeof(double));
      stdlhs = new double [nlc*ndof];
      memset (stdlhs,0,(nlc*ndof)*sizeof(double));
      rhs = new double [nlc*2*ndof];
      memset (rhs,0,(nlc*2*ndof)*sizeof(double));
    }
    break;
  }
  case mat_nonlinear_statics:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*2*ndof];
    memset (rhs,0,(nlc*2*ndof)*sizeof(double));
    break;
  }
  case geom_nonlinear_statics:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*2*ndof];
    memset (rhs,0,(nlc*2*ndof)*sizeof(double));
    break;
  }
  case mech_timedependent_prob:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [2*nlc*ndof];
    memset (rhs,0,(2*nlc*ndof)*sizeof(double));
    break;
  }
  case growing_mech_structure:{
    lhs = new double [ndof];
    memset (lhs,0,(ndof)*sizeof(double));
    rhs = new double [ndof];
    memset (rhs,0,(ndof)*sizeof(double));
    break;
  }
  case earth_pressure:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*2*ndof];
    memset (rhs,0,(nlc*2*ndof)*sizeof(double));
    break;
  }
  case layered_linear_statics:{    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*ndof];
    memset (rhs,0,(nlc*ndof)*sizeof(double));
    break;
  }
  case lin_floating_subdomain:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*ndof];
    memset (rhs,0,(nlc*ndof)*sizeof(double));
    break;
  }
  case nonlin_floating_subdomain:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [2*nlc*ndof];
    memset (rhs,0,(2*nlc*ndof)*sizeof(double));
    break;
  }
  case var_stiff_method:{
    lhs = new double [nlc*ndof];
    memset (lhs,0,(nlc*ndof)*sizeof(double));
    rhs = new double [nlc*ndof];
    memset (rhs,0,(nlc*ndof)*sizeof(double));
    break;
  }
    
  default:{
    fprintf (stderr,"\n\n unknown problem type is required in function lhsrhs::alloc (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}

/**
   function returns pointer to %vector of unknowns (lhs %vector)
   this %vector contains nodal values
   
   @param i - load case id
   
   JK
*/
double *lhsrhs::give_lhs (long i)
{
  return (lhs+i*ndof);
}


/**
   function returns pointer to %vector of first time derivatives of unknowns (tdlhs %vector)
   this %vector contains first time derivatives of nodal values
   
   @param i - load case id
   
   JK
*/
double *lhsrhs::give_tdlhs (long i)
{
  return (tdlhs+i*ndof);
}

/**
   function returns pointer to %vector of second time derivatives of unknowns (stdlhs %vector)
   this %vector contains second time derivatives of nodal values
   
   @param i - load case id
   
   JK
*/
double *lhsrhs::give_stdlhs (long i)
{
  return (stdlhs+i*ndof);
}

/**
   function returns pointer to %vector of right hand side (rhs %vector)
   this %vector contains prescribed nodal values
   
   @param i - load case id
   
   JK
*/
double *lhsrhs::give_rhs (long i)
{
  return (rhs+i*ndof);
}


void lhsrhs::output (FILE *out,long lcid)
{
  print_displacements (out,lcid);
  if (Mp->straincomp==1)  print_strains (out,lcid);
  if (Mp->stresscomp==1)  print_stresses (out,lcid);
  if (Mp->reactcomp==1)   print_reactions (out,lcid);
}

/**
   function reads initial conditions
   
   @param in - input data stream 

   17.1.2002
*/
void lhsrhs::initcond (XFILE *in)
{
  long i,j,k,l,ndofn;
  double v;
  
  if (Mp->tprob == forced_dynamics && Mb->dlc[0].tdl!=responsespectrum){
    for (i=0;i<Mt->nn;i++){
      ndofn=Mt->give_ndofn (i);
      xfscanf (in,"%ld",&l);
      
      //  nodal values
      for (j=0;j<ndofn;j++){
	k=Gtm->give_dof (l-1,j) - 1;
	xfscanf (in,"%lf",&v);
	if (k>-1)  lhs[k]=v;
      }
      
      //  time derivatives of nodal values
      for (j=0;j<ndofn;j++){
	k=Gtm->give_dof (l-1,j) - 1;
	xfscanf (in,"%lf",&v);
	if (k>-1)  tdlhs[k]=v;
      }
      
    }
  }
}

/**
   function cleans %vector of unknowns
   
   JK
*/
void lhsrhs::clean_lhs ()
{
  long i;
  for (i=0;i<ndof;i++){
    lhs[i]=0.0;
  }
}