sequent.cpp

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

#include <string.h>
#include <limits.h>
#include <math.h>
#include "sequent.h"
#include "global.h"
#include "matrix.h"
#include "vector.h"
#include "gtopology.h"
#include "mathem.h"
#include "element.h"
#include "plelemlt.h"
#include "plelemqt.h"
#include "plelemrotlt.h"
#include "plelemlq.h"
#include "plelemqq.h"
#include "plelemrotlq.h"
#include "lintet.h"
#include "linhex.h"
#include "quadhex.h"




/**
   function creates list of adjacent integration points
   and their distances
   
   results are stored in object Mt (class mechtop)
   Mt->nadjip[i] - number of adjacent integration points to the i-th integration point
   Mt->adjip[i] - array of numbers of adjacent integration points to the i-th integration point
   Mt->dist[i] - array of distances of adjacent integration points
                 to the i-th integration point
*/
void adjacip (void)
{
  long i,j,k,l,ii,jj,nb,ipp,nip,totnip,nae,nnae,nte;
  long *adjelel,*adjelelback,*treatedel,*change;
  double rlim;
  vector coord(3),auxcoord(3);
  
  //  total number of integration points
  totnip=Mm->tnip;
  Mt->nadjip = new long [totnip];
  Mt->adjip = new long* [totnip];
  Mt->dist = new double* [totnip];
  
  for (i=0;i<Gtm->ne;i++){
    nb = Mt->give_nb (i);
    
    for (ii=0;ii<nb;ii++){
      for (jj=0;jj<nb;jj++){
	
	ipp=Mt->elements[i].ipp[ii][jj];
	nip = Mt->give_nip (i,ii,jj);
	
	for (j=0;j<nip;j++){
	  //  radius of nonlocal volume
	  rlim = Mm->nonlocradius (ipp,0);
	  
	  //  coordinates of integration points
	  ipcoord (i,ipp,ii,jj,coord);
	  
	  //  number of adjacent integration points determination
	  Mt->nadjip[ipp]=0;
	  //  number of adjacent elements to required element
	  nae=Gtm->nadjelel[i];
	  nte=nae;
	  //  list of adjacent elements to required element
	  adjelel = new long [nae];
	  treatedel = new long [nae];
	  for (k=0;k<nae;k++){
	    adjelel[k]=Gtm->adjelel[i][k];
	    treatedel[k]=adjelel[k];
	  }
	  
	  do{
	    in_dist (ipp,ii,jj,coord,adjelel,nae,rlim);
	    adjelelback = new long [nae];
	    newnadjelel (ipp,adjelelback,adjelel,nae,nnae);
	    if (nnae>0){
	      delete [] adjelel;
	      adjelel = new long [nnae];
	      newadjelel (adjelelback,adjelel,nae,nnae,treatedel,nte);
	      delete [] adjelelback;
	      change = new long [nte];
	      for (k=0;k<nte;k++){
		change[k]=treatedel[k];
	      }
	      delete [] treatedel;
	      treatedel = new long [nte+nnae];
	      for (k=0;k<nte;k++){
		treatedel[k]=change[k];
	      }
	      l=0;
	      for (k=nte;k<nte+nnae;k++){
		treatedel[k]=adjelel[l];  l++;
	      }
	      delete [] change;
	      nte+=nnae;
	      nae=nnae;
	    }
	    else{}
	  }
	  while (nnae!=0);
	  
	  delete [] adjelel;
	  delete [] treatedel;
	  
	  
	  //  adjacent integration points determination
	  Mt->adjip[ipp] = new long [Mt->nadjip[ipp]];
	  Mt->dist[ipp] = new double [Mt->nadjip[ipp]];
	  Mt->nadjip[ipp]=0;
	  nae=Gtm->nadjelel[i];
	  nte=nae;
	  adjelel = new long [nae];
	  treatedel = new long [nae];
	  for (k=0;k<nae;k++){
	    adjelel[k]=Gtm->adjelel[i][k];
	    treatedel[k]=adjelel[k];
	  }
	  
	  do{
	    dist (ipp,ii,jj,coord,adjelel,nae,rlim);
	    adjelelback = new long [nae];
	    newnadjelel (ipp,adjelelback,adjelel,nae,nnae);
	    if (nnae>0){
	      delete [] adjelel;
	      adjelel = new long [nnae];
	      newadjelel (adjelelback,adjelel,nae,nnae,treatedel,nte);
	      delete [] adjelelback;
	      change = new long [nte];
	      for (k=0;k<nte;k++){
		change[k]=treatedel[k];
	      }
	      delete [] treatedel;
	      treatedel = new long [nte+nnae];
	      for (k=0;k<nte;k++){
		treatedel[k]=change[k];
	      }
	      l=0;
	      for (k=nte;k<nte+nnae;k++){
		treatedel[k]=adjelel[l];  l++;
	      }
	      delete [] change;
	      nte+=nnae;
	      nae=nnae;
	    }
	  }
	  while (nnae!=0);
	  
	  delete [] adjelel;
	  delete [] treatedel;
	  
	  ipp++;
	}
      }
    }
  }
  
  if (Mespr==1){
    fprintf (Out,"\n\n\n LIST OF ADJACENT INTEGRATION POINTS \n");
    for (i=0;i<totnip;i++){
      fprintf (Out,"\n %4ld   %4ld  :",i,Mt->nadjip[i]);
      
      for (j=0;j<Mt->nadjip[i];j++){
	fprintf (Out,"     %ld %f",Mt->adjip[i][j],Mt->dist[i][j]);
      }
      
    }
  }
  
  
}

/**
   function calculates number of non-eliminated adjacent elements
   function is used in function adjacip
   
   @param ipp - integration point pointer
   @param adjelelback - auxiliary array of numbers of adjacent elements
   @param adjelel - array of numbers of adjacent elements
   @param nae - number of adjacent elements
   @param nnae - new number of adjacent elements
*/
void newnadjelel (long ipp,long *adjelelback,long *adjelel,long &nae,long &nnae)
{
  long i,nee;
  
  //  number of eliminated elements
  nee=0;
  for (i=0;i<nae;i++){
    if (adjelel[i]<0)  nee++;
  }
  
  if (nae!=nee){
    nnae=0;
    for (i=0;i<nae;i++){
      adjelelback[i]=adjelel[i];
      if (adjelel[i]>-1){
	nnae+=Gtm->nadjelel[adjelel[i]];
      }
    }
  }
  else{
    nnae=0;
  }
}

/**
   function creates list of adjacent elements to required element
   function is used in function adjacip
   
   @param adjelelback - auxiliary array of adjacent elements
   @param adjelel - array of adjacent elements
   @param nae - number of adjacent elements
   @param nnae - new number of adjacent elements
   @param treatedel - array of treated elements
   @param nte - number of treated elements
   
*/
void newadjelel (long *adjelelback,long *adjelel,long &nae,long &nnae,long *treatedel,long &nte)
{
  long i,j,k,l,prev,min;
  
  j=0;
  for (i=0;i<nae;i++){
    if (adjelelback[i]>-1){
      for (k=0;k<Gtm->nadjelel[adjelelback[i]];k++){
	adjelel[j]=Gtm->adjelel[adjelelback[i]][k];
	j++;
      }
    }
  }
  
  //  sorting
  prev=-1;
  for (i=0;i<nnae;i++){
    min=LONG_MAX;
    for (j=i;j<nnae;j++){
      if (min>adjelel[j]){
	min=adjelel[j];  k=j;
      }
    }
    if (min==prev){
      nnae--;  i--;
      adjelel[k]=adjelel[nnae];
    }
    else{
      l=adjelel[i];
      adjelel[i]=min;
      adjelel[k]=l;
      prev=min;
    }
  }
  
  for (i=0;i<nnae;i++){
    k=adjelel[i];
    for (j=0;j<nte;j++){
      if (treatedel[j]==k){
	nnae--;
	adjelel[i]=adjelel[nnae];
	i--;
	break;
      }
    }
  }

}

/**
   function computes number of integration points which are closer
   to required integration point (ipp) than rlim
   
   function is used in adjacip
   
   @param ipp - integration point pointer
   @param ri,ci - row and column indices
   @param coord - coordinates of one of integration points
   @param adjelel - array of adjacent elements
   @param nae - number of adjacent elements
   @param rlim - limit distance
   
*/
void in_dist (long ipp,long ri,long ci,vector &coord,long *adjelel,long &nae,double rlim)
{
  long i,j,lj,uj,ii,eid,nip;
  double r;
  vector auxcoord(3);
  
  for (i=0;i<nae;i++){
    eid = adjelel[i];
    nip = Mt->give_nip (eid,ri,ci);
    lj = Mt->elements[eid].ipp[ri][ci];
    uj = lj+nip;
    
    ii=0;
    for (j=lj;j<uj;j++){
      ipcoord (eid,j,ri,ci,auxcoord);
      r = length (coord,auxcoord);
      
      if (r<rlim){
        Mt->nadjip[ipp]++;
	ii++;
      }
    }

    if (ii==0)  adjelel[i]=-1-adjelel[i];
    
  }
}

/**
   function computes distances of integration points which are closer
   to required integration point (ipp) than rlim
   
   function is used in adjacip
   
   @param ipp - integration point pointer
   @param ri,ci - row and column indices
   @param coord - coordinates of one of integration points
   @param adjelel - array of adjacent elements
   @param nae - number of adjacent elements
   @param rlim - limit distance
   
*/
void dist (long ipp,long ri,long ci,vector &coord,long *adjelel,long &nae,double rlim)
{
  long i,j,lj,uj,ii,eid,nip;
  double r;
  vector auxcoord(3);
  
  for (i=0;i<nae;i++){
    eid = adjelel[i];
    nip = Mt->give_nip (eid,ri,ci);
    lj = Mt->elements[eid].ipp[ri][ci];
    uj = lj+nip;
    
    ii=0;
    for (j=lj;j<uj;j++){
      ipcoord (eid,j,ri,ci,auxcoord);
      r = length (coord,auxcoord);

      if (r<rlim){
	Mt->adjip[ipp][Mt->nadjip[ipp]]=j;
	Mt->dist[ipp][Mt->nadjip[ipp]]=r;
	Mt->nadjip[ipp]++;
	ii++;
      }
    }
    
    if (ii==0)  adjelel[i]=-1-adjelel[i];
  }
}

/**
   function returns coordinates of integration point (ipp) on
   element (eid)
   
   function is used in adjacip
   
   @param eid - element id
   @param ipp - integration point pointer
   @param ri,ci - row and column indices
   @param ipcoord - array containing coordinates of integration point
   
*/
void ipcoord (long eid,long ipp,long ri,long ci,vector &ipcoord)
{
  switch (Mt->give_elem_type(eid)){
    
  case planeelementlt:{
    Pelt->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case planeelementqt:{
    Peqt->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case planeelementrotlt:{
    Perlt->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case planeelementlq:{
    Pelq->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case planeelementqq:{
    Peqq->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case planeelementrotlq:{
    Perlq->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }

  case lineartet:{
    Ltet->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case linearhex:{
    Lhex->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  case quadrhex:{
    Qhex->ipcoord (eid,ipp,ri,ci,ipcoord);
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown element type is required in function ipcoord (%s, line %d).\n",__FILE__,__LINE__);
  }
  }

}