dloadcase.cpp

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

#include <math.h>
#include "dloadcase.h"
#include "loadcase.h"
#include "global.h"
#include "gtopology.h"
#include "matrix.h"
#include "vector.h"
#include "globmat.h"
#include "dloadn.h"
#include "dloadpd.h"
#include "element.h"
#include "node.h"



dloadcase::dloadcase (void)
{
  nslc=0;
  nln = 0;  nle = 0;  npd = 0;
  lon = NULL;  pd = NULL;
  slc = NULL;  gf=NULL;
}

dloadcase::~dloadcase (void)
{
  delete [] lon;  delete [] pd;
  delete [] slc;  delete [] gf;
}

/**
   function reads load case characteristics
   
   @param in - input stream
   
   24.7.2001
*/
void dloadcase::read (XFILE *in)
{
  long i;
  
  //  type of dynamic load
  xfscanf (in,"%k%d","type_of_dyn_load",(int*)&tdl);
  
  switch (tdl){
  case forcedload:{

    //  number of subload cases
    xfscanf (in,"%ld",&nslc);
    
    slc = new loadcase [nslc];
    gf = new gfunct [nslc];
    
    for (i=0;i<nslc;i++){
      slc[i].read (in);
      gf[i].read (in);
    }
    
    break;
  }
  case seismicload:{
    stool.read (in);
    break;
  }
  case responsespectrum:{
    stool.read (in);
    break;
  }
  case forceload_dyn:{
    
    //  loaded nodes
    xfscanf (in,"%ld",&nln);
    lon = new dloadn [nln];
    for (i=0;i<nln;i++){
      lon[i].read (in);
    }
    
    
    /*  
    //  loaded elements
    fscanf (in,"%ld",&nle);
    loe = new dloadel [nle];
    for (i=0;i<nle;i++){
    loe[i].read (in);
    }*/
    
    
    
    //  prescribed displacements
    xfscanf (in,"%ld",&npd);
    pd = new dloadpd [npd];
    for (i=0;i<npd;i++){
      pd[i].read (in);
    }
    
    break;
  }
    
  default:{
    fprintf (stderr,"\n\n unknown type of load case is required in function read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }

}



/**
   function assembles vector of right hand side
   
   @param lcid - load case id
   @param rhs - right hand side vector
   @param t - actual time

   26.10.2001
*/
void dloadcase::assemble (long lcid,double *rhs,long n,double t)
{
  long i,j,k,ndofn,ndofe,*cn;
  vector r,f;
  matrix sm;
  
  double a,*aux;
  
  // nullv (rhs,Ndofm);
  
  switch (tdl){
    /*
  case forcedload:{
    //  auxiliary array
    aux = new double [n];
    tottempr = NULL;      
    for (i=0;i<nslc;i++){
      nullv (aux,n);
      btemp = 0;
      if (slc[i].tempchang == 1)
      // temperature load needs separate time coeficient multiplication
      {
        btemp = 1;
        slc[i].tempchang = 0;
      }
      //  assembling of vector of contributions
      slc[i].assemble (0,aux);
      //  scale factor
      a=gf[i].getval(t);
      
      if (btemp==1)
      // temperature load is present in the current subloadcase
      {
        //  approximation of nodal temperatures
        intpointval (slc[i].pt, temperature);
        Mm->temprstrains (lcid);
        nodal_eigstrain_forces (lcid, aux);
        slc[i].tempchang = 1;
      }
      //  updating of rhs array
      for (j=0;j<n;j++){
	rhs[j]+=a*aux[j];
      }
      if ((tottempr == NULL) && (Mp->temperature))
      {
        tottempr = new double[Mm->tnip];
        nullv(tottempr, Mm->tnip);
      }
      if (Mp->temperature)
      {
        for(j=0; j<Mm->tnip; j++)
          tottempr[j] += a*Mm->tempr[j];
      }
    }
    if (Mp->temperature)
    {
      for(j=0; j<Mm->tnip; j++)
        Mm->tempr[j] = tottempr[j];
    }

    delete [] tottempr;
    delete [] aux;

    break;
  }
  */

  case forcedload:{
    aux = new double [n];
    
    for (i=0;i<nslc;i++){
      nullv (aux,n);
      //  scale factor
      a=gf[i].getval(t);
      
      //  assembling of vector of contributions
      slc[i].assemble (0,aux,a);
      //  updating of rhs array
      for (j=0;j<n;j++){
	rhs[j]+=aux[j];
      }
    }

    delete [] aux;

    break;
  }
  case seismicload:{
    
    //  auxiliary array
    aux = new double [n];
    
    stool.assemble (rhs,t);

    Mmat->gmxv (rhs,aux);
    
    copyv (aux,rhs,n);
    
    delete [] aux;
    
    break;
  }
  case responsespectrum:{
    //stool.assemble (rhs,t);
    break;
  }
  case forceload_dyn:
    {
      
      //  contributions from nodes
      for (i=0;i<nln;i++){
	ndofn=Mt->give_ndofn (lon[i].idn);
	for (j=0;j<ndofn;j++){
	  k=Mt->give_dof (lon[i].idn,j);
	  if (k<=0)  continue;
	  if (k>0)  rhs[k-1]+=lon[i].getval(t,j);
	}
      }
      
      
      /*
      //  contributions from elements
      for (i=0;i<nle;i++){
      ndofe=Mt->elements[loe[i].ide].ndofe;
      cn = new long [ndofe];
      Mt->give_code_numbers (loe[i].ide,cn);
      
      for (j=0;j<ndofe;j++){
      k=cn[j]-1;
      if (k<0)  continue;
      else  rhs[k]+=loe[i].fe[j];
      }
      }
      */
      
      
      //  contributions from prescribed displacements
      if (npd>0){
	long ne = Mt->ne;
	for (i=0;i<ne;i++){
	  if (Mt->elements[i].prescdispl==1){
	    ndofe=Mt->give_ndofe (i);
	    allocm (ndofe,ndofe,sm);
	    stiffmat (i,lcid,sm);
	    
	    r.a = new double [ndofe];  r.n=ndofe;
	    f.a = new double [ndofe];  f.n=ndofe;
	    cn = new long [ndofe];
	    Mt->give_code_numbers (i,cn);
	    eldispl (lcid,i,r.a,cn,ndofe);
	    mxv (sm,r,f);
	    cmulv (-1.0,f);
	    locglob (rhs,f.a,cn,ndofe);
	    destrm (sm);
	    delete [] cn;
	    //delete [] f.x;  delete [] r.x;
	  }
	}
      }
      
      break;
    }
    
  default:{
    fprintf (stderr,"\n\n unknown type of load case is required in function assemble (file %s, line %d).\n",
	     __FILE__,__LINE__);
  }
  }
  
}

/**
   function assembles vector of right hand side with respect of left hand side
   
   @param rhs - right hand side vector
   @param lhs - left hand side

   26.10.2001
*/
void dloadcase::assemble (double *rhs, double *lhs)
{
  long i,j,k,ndofn;
  vector r,f;
  matrix sm;

  //nullv (rhs,Ndofm);

  //  contributions from nodes
  for (i=0;i<nln;i++){
    ndofn=Mt->give_ndofn (lon[i].idn);
    for (j=0;j<ndofn;j++){
      k=Mt->give_dof (lon[i].idn,j);
      if (k<=0)  continue;
      if (k>0)  rhs[k-1]+=lon[i].getval(lhs[k-1], j);
    }
  }
}


/**
   function evaluates reactions
   
   @param lcid - load case id
   
*/
void dloadcase::compute_reactions (long lcid)
{
  long i,j,k,nne,ndofn,ne,te,ndofe,*cn;
  double *r,*f;
  ivector *nod;
  matrix sm;

  ne=Mt->ne;
  for (i=0;i<ne;i++){
    if (Mt->elements[i].react==1){
      te = Mt->give_elem_type (i);
      ndofe=Mt->give_ndofe (i);
      allocm (ndofe,ndofe,sm);

      stiffmat (i,lcid,sm);
      
      r = new double [ndofe];
      f = new double [ndofe];
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);
      eldispl (lcid,i,r,cn,ndofe);
      mxv (sm.a,r,f,ndofe,ndofe);
      nne=Mt->give_nne (i);
      
      nod = new ivector (nne);
      Mt->give_elemnodes (i,*nod);
      long ii=0;
      for (j=0;j<nne;j++){
	ndofn=Mt->give_ndofn (nod->a[j]);
	if (Mt->nodes[nod->a[j]].react==1){
	  for (k=0;k<ndofn;k++){
	    Mt->nodes[nod->a[j]].r[k]+=f[ii];
	    ii++;
	  }
	}
	else  ii+=ndofn;
      }
      
      destrm (sm);
      delete [] r;  delete [] f;  delete [] cn;  delete [] nod;
    }
  }
  
}