loadcase.cpp

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

#include <math.h>
#include "loadcase.h"
#include "global.h"
#include "globmat.h"
#include "elemswitch.h"
#include "gtopology.h"
#include "matrix.h"
#include "vector.h"
#include "loadn.h"
#include "loadel.h"
#include "element.h"
#include "node.h"
#include "elemhead.h"


loadcase::loadcase (void)
{
  nln = 0;  nle = 0;  npd = 0;  npt = 0;  tempchang=0;
  lon = NULL;  loe = NULL;  pd = NULL;  pt = NULL;
}

loadcase::~loadcase (void)
{
  delete [] lon;  delete [] loe;  delete [] pd;  delete [] pt;
}

/**
   function reads load case characteristics
   
   24.7.2001
*/
void loadcase::read (XFILE *in)
{
  long i;
  
  //  loaded nodes
  xfscanf (in,"%ld",&nln);
  lon = new loadn [nln];
  for (i=0;i<nln;i++){
    lon[i].read (in);
  }
  
  //  loaded elements
  xfscanf (in,"%ld",&nle);
  loe = new loadel [nle];
  for (i=0;i<nle;i++){
    loe[i].read (in);
  }
  
  //  prescribed displacements
  xfscanf (in,"%ld",&npd);
  pd = new double [npd];
  for (i=0;i<npd;i++){
    xfscanf (in,"%lf",pd+i);
  }
  
  //  prescribed temperature changes
  xfscanf (in,"%ld",&tempchang);
  if (tempchang>0)
    Mp->temperature=1;
  
  if (tempchang==1 || tempchang==2){
    npt = Mt->nn;
    pt = new double [npt];
    for (i=0;i<npt;i++){
      xfscanf (in,"%lf",pt+i);
    }
    
  }
  
}

/**
   function assembles vector of right hand side
   
   @param lcid - load case id
   @param rhs - pointer to the right hand side
   @param scale - scale factor, it is used for temperature scaling
                  it is useless for other loads
   24.7.2001
*/
void loadcase::assemble (long lcid,double *rhs,double scale)
{
  long i,j,k,l,ndofn,ndofe,*cn;
  vector r,f;
  matrix sm;

  //nullv (rhs,Ndofm);

  //  contributions from nodes
  for (i=0;i<nln;i++){
    ndofn=Mt->give_ndofn (lon[i].nid);
    for (j=0;j<ndofn;j++){
      k=Mt->give_dof (lon[i].nid,j);
      if (k<=0)  continue;
      if (k>0)  rhs[k-1]+=lon[i].f[j];
    }
  }
    
  //TKr:
  for (i=0;i<Ndofm;i++){
    rhs[i]*=scale;
  }

  //  contributions from elements
  for (i=0;i<nle;i++){
    l=loe[i].eid;
    if (Gtm->leso[l]==1){
      ndofe=Mt->give_ndofe (loe[i].eid);
      cn = new long [ndofe];
      Mt->give_code_numbers (loe[i].eid,cn);
      
      for (j=0;j<ndofe;j++){
	k=cn[j]-1;
	if (k<0)  continue;
	else  rhs[k]+=loe[i].nf[j];
      }
      delete [] cn;
    }
  }
  
  //  contributions from prescribed displacements
  if ((npd>0) && (Mp->tprob != mat_nonlinear_statics)){
    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 (lcid,i,sm);
	
	allocv (ndofe,r);
	allocv (ndofe,f);
	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;
	destrv (f);  destrv (r);
      }
    }
  }
  
  //JK:
  //for (i=0;i<Ndofm;i++){
  //rhs[i]*=scale;
  //}
  
  //  contributions from temperature changes
  if (tempchang==1){
    Mm->est=tempstrain;
    //  approximation of nodal temperatures to integration points
    intpointval (pt,temperature,scale);
    
    Mm->temprstrains (lcid);

    nodal_eigstrain_forces (lcid,rhs);
  }
  
  if (tempchang==2){
    Mm->est=tempstrain;
    //  approximation of nodal temperatures
    intpointval (pt,temperature,scale);
    
    Mm->temprstrains (lcid);
    
    //nodal_eigstrain_forces (lcid,rhs);
  }

  if (tempchang==3){
    Mm->est=tempstrain;
    
    Mm->temprstrains (lcid);
    
    //nodal_eigstrain_forces (lcid,rhs);
  }
}

void loadcase::compute_reactions (long lcid)
{
  long i,j,k,l,ii,nne,ndofn,nn,ne,ndofe,*cn;
  double *r;
  elemtype te;
  ivector nod;
  vector f,tf;
  matrix sm;
  
  nn=Mt->nn;
  for (i=0;i<nn;i++){
    if (Mt->nodes[i].react==1){
      ndofn=Mt->give_ndofn (i);
      for (k=0;k<ndofn;k++){
	Mt->nodes[i].r[k]=0.0;
      }
    }
  }
  
  
  ne=Mt->ne;
  for (i=0;i<ne;i++){
    if (Mt->elements[i].react==1){
      ndofe=Mt->give_ndofe (i);
      allocm (ndofe,ndofe,sm);
      stiffmat (lcid,i,sm);
      
      allocv (ndofe,f);
      r = new double [ndofe];
      cn = new long [ndofe];
      Mt->give_code_numbers (i,cn);

      //eldispl (lcid,i,r,cn,ndofe);
      //mxv (sm.a,r,f.a,ndofe,ndofe);
      elem_internal_forces (i,lcid,f);
      
      
      if (Mt->elements[i].presctemp==1){
	allocv (ndofe,tf);
	elem_eigstrain_forces (lcid,i,tf);
	subv (f,tf,f);
	destrv (tf);
      }

      nne=Mt->give_nne (i);
      allocv (nne,nod);
      Mt->give_elemnodes (i,nod);
      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);
      destrv (f);
      delete [] r;
      delete [] cn;
      destrv (nod);
    }
  }
  
  
  for (i=0;i<nle;i++){
    l=loe[i].eid;
    if (Mt->elements[l].react==1){
      ndofe=Mt->give_ndofe (l);
      cn = new long [ndofe];
      Mt->give_code_numbers (l,cn);
      nne=Mt->give_nne (l);
      allocv (nne,nod);
      Mt->give_elemnodes (l,nod);
      
      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]-=loe[i].nf[ii];
	    ii++;
	  }
	}
	else  ii+=ndofn;
      }
      
      destrv (nod);
      delete [] cn;
    }
  }
  
  
  
  if (Mp->eigstrains==1){
    
    for (i=0;i<ne;i++){
      if (Mt->elements[i].react==1){
	te = Mt->give_elem_type (i);
	ndofe=Mt->give_ndofe (i);
	allocv (ndofe,f);
	
	nne=Mt->give_nne (i);
	allocv (nne,nod);
	Mt->give_elemnodes (i,nod);
	

	elem_eigstrain_forces (lcid,i,f);
	
	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;
	}
	
	destrv (f);
      }
    }
  }
  
  
}