mechtop.cpp

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <limits.h>
#include <math.h>
#include "mechtop.h"
#include "global.h"
#include "elemhead.h"
#include "gtopology.h"
#include "difcalc.h"
#include "mathem.h"
#include "element.h"
#include "node.h"
#include "siftop_element_types.h"
#include "loadcase.h"
#include "intpoints.h"
#include "globmat.h"
#include "siftop.h"





mechtop::mechtop ()
{
  nn=0;  ncn=0;  ne=0;
  nodes=NULL;  elements=NULL;
  nadjip=NULL;  adjip=NULL;  dist=NULL;
  nodedispl=NULL;
}

mechtop::~mechtop ()
{
  long i;
  
  delete [] nodes;  
  delete [] elements;
  delete [] nadjip;  
  delete [] adjip;  
  delete [] dist;
  
  if (nodedispl != NULL){
    for (i=0;i<nn;i++){
      delete [] nodedispl[i];
    }
    delete [] nodedispl;
  }
}

/**
   function reads basic topology informations
   
   @param in - input stream
*/
void mechtop::read (XFILE *in)
{
  long i,j,k,l,ndofn,nne,dim,order;
  
  if (Mp->tprob == growing_mech_structure)
    Gtm->dofcontr=1;
  
  // ********************
  //  node data reading
  // ********************
  xfscanf (in,"%k%ld","number_of_nodes",&nn);
  Gtm->alloc_nodes (nn);
  if (Mespr==1)  fprintf (stdout,"\n number of nodes  %ld",nn);
  nodes = new node [nn];
  
  switch (Mp->tprob){
  case linear_statics:
  case eigen_dynamics:
  case forced_dynamics:
  case linear_stability:
  case mat_nonlinear_statics:
  case geom_nonlinear_statics:
  case mech_timedependent_prob:
  case growing_mech_structure:
  case var_stiff_method:
  case earth_pressure:{
    for (i=0;i<nn;i++){
      j=nn+1;
      xfscanf (in,"%ld",&j);
      if (j<1)
	fprintf (stderr,"\n\n node number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
      if (j>nn)
	fprintf (stderr,"\n\n node number in function read (file %s, line %d) is greater than number of nodes.\n",__FILE__,__LINE__);
      Gtm->gnodes[j-1].read (in);
      ndofn = Gtm->gnodes[j-1].ndofn;
      nodes[j-1].read (in,ndofn);
    }
    break;
  }
  case layered_linear_statics:{
    xfscanf (in,"%k%ld","number_of_layered_nodes",&nln);
    Gtm->alloc_lnodes (nln);
    for (i=0;i<nln;i++){
      xfscanf (in,"%ld",&l);
      xfscanf (in,"%ld",&(Gtm->lgnodes[l-1].nl));
      Gtm->lgnodes[l-1].nodes = new long [Gtm->lgnodes[l-1].nl];
      for (k=0;k<Gtm->lgnodes[l-1].nl;k++){
	j=nn+1;
	xfscanf (in,"%ld",&j);
	if (j<1)
	  fprintf (stderr,"\n\n node number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
	if (j>nn)
	  fprintf (stderr,"\n\n node number in function read (file %s, line %d) is greater than number of nodes.\n",__FILE__,__LINE__);
	Gtm->gnodes[j-1].read (in);
	ndofn = Gtm->gnodes[j-1].ndofn;
	nodes[j-1].read (in,ndofn);
	Gtm->lgnodes[l-1].nodes[k]=j-1;
      }
    }
    Gtm->unodelnode ();
    break;
  }
  case lin_floating_subdomain:
  case nonlin_floating_subdomain:{
    //  number of subdomains
    xfscanf (in,"%k%ld","number_of_floating_subdomains",&(Gtm->flsub.nsd));
    //  allocation of arrays
    Gtm->flsub.nnsd = new long [Gtm->flsub.nsd];
    Gtm->flsub.fnnsd = new long [Gtm->flsub.nsd+1];
    
    Gtm->flsub.fnnsd[0]=0;
    for (i=0;i<Gtm->flsub.nsd;i++){
      xfscanf (in,"%ld",&(Gtm->flsub.nnsd[i]));
      Gtm->flsub.fnnsd[i+1] = Gtm->flsub.fnnsd[i]+Gtm->flsub.nnsd[i];
    }
    Gtm->flsub.fnnsd[Gtm->flsub.nsd]--;

    //  nodes
    for (i=0;i<nn;i++){
      j=nn+1;
      xfscanf (in,"%ld",&j);
      if (j<1)
	fprintf (stderr,"\n\n node number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
      if (j>nn)
	fprintf (stderr,"\n\n node number in function read (file %s, line %d) is greater than number of nodes.\n",__FILE__,__LINE__);
      Gtm->gnodes[j-1].read (in);
      ndofn = Gtm->gnodes[j-1].ndofn;
      nodes[j-1].read (in,ndofn);
    }
    
    
    //  layered nodes
    xfscanf (in,"%ld",&nln);
    Gtm->alloc_lnodes (nln);
    for (i=0;i<nln;i++){
      xfscanf (in,"%ld",&l);
      xfscanf (in,"%ld",&(Gtm->lgnodes[l-1].nl));
      Gtm->lgnodes[l-1].nodes = new long [Gtm->lgnodes[l-1].nl];
      for (k=0;k<Gtm->lgnodes[l-1].nl;k++){
	j=nn+1;
	xfscanf (in,"%ld",&j);
	if (j<1)
	  fprintf (stderr,"\n\n node number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
	if (j>nn)
	  fprintf (stderr,"\n\n node number in function read (file %s, line %d) is greater than number of nodes.\n",__FILE__,__LINE__);
	Gtm->lgnodes[l-1].nodes[k]=j-1;
      }
    }
    //Gtm->unodelnode ();
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of problem is required in function read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }



  //  constrained nodes
  xfscanf (in,"%k%ld","number_of_constraints",&ncn);
  if (Mespr==1)  fprintf (stdout,"\n number of constrained nodes  %ld",ncn);
  for (i=0;i<ncn;i++){
    k=nn+1;
    xfscanf (in,"%ld",&k);
    if (k<1)
      fprintf (stderr,"\n\n number of constrained node in function mechtop::read is less than 1.\n");
    if (k>nn){
      fprintf (stderr,"\n\n number of constrained node in function mechtop::read");
      fprintf (stderr,"\n is greater than number of all nodes.\n");
    }
    Gtm->gnodes[k-1].constr (Gtm->dofcontr,in);
  }
  

  // ***********************
  // ***********************
  //  element data reading
  // ***********************
  // ***********************
  xfscanf (in,"%k%ld","number_of_elements",&ne);
  Gtm->alloc_elements (ne);
  if (Mespr==1)  fprintf (stdout,"\n number of elements  %ld",ne);
  elements = new element [ne];

  switch (Mp->tprob){
  case linear_statics:
  case eigen_dynamics:
  case forced_dynamics:
  case linear_stability:
  case mat_nonlinear_statics:
  case geom_nonlinear_statics:
  case mech_timedependent_prob:
  case growing_mech_structure:
  case earth_pressure:
  case var_stiff_method:
  case layered_linear_statics:{
    for (i=0;i<ne;i++){
      j=ne+1;
      xfscanf (in,"%ld",&j);
      if (j<1)
	fprintf (stderr,"\n\n element number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
      if (j>ne){
	fprintf (stderr,"\n\n element number in function read (file %s, line %d) is",__FILE__,__LINE__);
	fprintf (stderr,"\n greater than total number of elements.\n");
      }
      
      elements[j-1].read (in,j-1);
    }
    break;
  }
  case lin_floating_subdomain:
  case nonlin_floating_subdomain:{

    Gtm->flsub.nesd = new long [Gtm->flsub.nsd];
    Gtm->flsub.fensd = new long [Gtm->flsub.nsd+1];
    Gtm->flsub.fensd[0]=0;
    //  number of elements on particular subdomains
    for (i=0;i<Gtm->flsub.nsd;i++){
      xfscanf (in,"%ld",&(Gtm->flsub.nesd[i]));
      Gtm->flsub.fensd[i+1]=Gtm->flsub.fensd[i]+Gtm->flsub.nesd[i];
    }
    
    if (ne!=Gtm->flsub.fensd[Gtm->flsub.nsd]){
      fprintf (stderr,"\n\n the sum of numbers of elements on subdomains is not equal to the total number of elements");
      fprintf (stderr,"\n in function read (file %s, line %d).\n",__FILE__,__LINE__);
    }
    
    for (i=0;i<ne;i++){
      j=ne+1;
      xfscanf (in,"%ld",&j);
      if (j<1)
	fprintf (stderr,"\n\n element number in function read (file %s, line %d) is less than 1.\n",__FILE__,__LINE__);
      if (j>ne){
	fprintf (stderr,"\n\n element number in function read (file %s, line %d) is",__FILE__,__LINE__);
	fprintf (stderr,"\n greater than total number of elements.\n");
      }
      
      elements[j-1].read (in,j-1);
    }
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of problem is required in function read (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
    
    
  /*
  if (Mp->tprob==layered_linear_statics){
    fscanf (in,"%ld",&nle);
    Gtm->alloc_lelements (nle);
    for (i=0;i<nle;i++){
      fscanf (in,"%ld",&l);
      fscanf (in,"%ld",&(Gtm->lgelements[l-1].nl));
      Gtm->lgelements[l-1].elem = new long [Gtm->lgelements[l-1].nl];
      for (k=0;k<Gtm->lgelements[l-1].nl;k++){
	j=ne+1;
	fscanf (in,"%ld",&j);
	if (j<1)
	  fprintf (stderr,"\n\n element number in function mechtop::read is less than 1.\n");
	if (j>ne){
	  fprintf (stderr,"\n\n element number in function mechtop::read is");
	  fprintf (stderr,"\n greater than total number of elements.\n");
	}
	elements[j-1].read (in,j-1);
	Gtm->lgelements[l-1].elem[k]=j-1;
      }
      Gtm->lgelements[l-1].elemnmult ();
    }
  }
  else{
  */

  

  if (Mp->adaptivity > 0){
    //  initialization of auxinf variable
    for (i=0;i<ne;i++){
      nne = give_nne (i);
      order = give_degree (i);
      dim = give_dimension (i);
      
      Gtm->gelements[i].auxinf = nne*100+order*10+dim;
    }
  }
  
  
  //  allocation and initiation of arrays lnso and leso
  //  lnso - list of nodes switched on
  //  leso - list of elements switched on
  //  3.11.2006
  Gtm->lneso_init ();


  if (Mp->tprob == growing_mech_structure){
    Gtm->auxinf_init ();
  }
  
}

/**
   function investigates presence of prescribed displacements on elements
   
   25.7.2001
*/
void mechtop::elemprescdisp (void)
{
  long i,j,ndofe,*cn;
  
  for (i=0;i<ne;i++){
    if (Gtm->leso[i]==1){
      ndofe=give_ndofe (i);
      cn = new long [ndofe];
      give_code_numbers (i,cn);
      for (j=0;j<ndofe;j++){
	if (cn[j]<0){
	  elements[i].prescdispl=1;
	  break;
	}
      }
      delete [] cn;
    }
  }
}

/**
   function investigates presence of temperature changes on elements
   
   30.11.2002
*/
void mechtop::elempresctemp (long lcid)
{
  long i;
  ivector nod;

  if (Mp->temperature==1 || Mp->temperature==2){
    for (i=0;i<ne;i++){
      if (Gtm->leso[i]==1){
	elements[i].presctemp=1;
      }
    }
  }
  
  /*
  if (Mb->lc[lcid].tempchang==1){
    for (i=0;i<ne;i++){
      nne=give_nne (i);
      allocv (nne,nod);
      dtelem = new double [nne];
      give_elemnodes (i,nod);
      Mb->lc[lcid].tempchanges (dtelem,nod);
      for (j=0;j<nne;j++){
	if (dtelem[j]!=0.0){
	  elements[i].presctemp=1;
	  break;
	}
      }
      delete [] dtelem;
      destrv (nod);
    }
  }
*/
}

/**
   function returns element type
   
   @param eid - element id
*/
elemtype mechtop::give_elem_type (long eid)
{
  return (elements[eid].te);
}

/**
   function returns number of DOFs of the element
   
   @param eid - element id
*/
long mechtop::give_ndofe (long eid)
{
  long ndofe=0;
  elemtype te;
  
  te = give_elem_type (eid);
  
  switch (te){
  case bar2d:{               ndofe=Bar2d->ndofe;              break; }
  case bar3d:{               ndofe=Bar3d->ndofe;              break; }
  case barq2d:{              ndofe=Barq2d->ndofe;             break; }
  case barq3d:{              ndofe=Barq3d->ndofe;             break; }
  case beam2d:{              ndofe=Beam2d->ndofe;             break; }
  case beam3d:{              ndofe=Beam3d->ndofe;             break; }