flsubdom.cpp

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

#include "flsubdom.h"
#include "gtopology.h"
#include "global.h"
#include "gmatrix.h"

flsubdom::flsubdom ()
{
  enrbm=0;  nlm=0;

  de = NULL;
  //  array for rigid body motions
  rbm = NULL;
  gg = NULL;
  igg = NULL;
  e = NULL;
  //  array of Lagrange multipliers
  w = NULL;
  //  array of total Lagrange multipliers
  tw = NULL;
  //  array of compliances
  compli=NULL;

  mu=NULL;
  muo=NULL;
  dmu=NULL;
  ddmu=NULL;

}

flsubdom::~flsubdom ()
{
  delete [] de;
  delete [] rbm;
  delete [] gg;
  delete [] igg;
  delete [] e;
  delete [] w;
  delete [] tw;
  delete [] compli;
  
  delete [] mu;
  delete [] muo;
  delete [] dmu;
  delete [] ddmu;
}

/**
   function initializes variables and arrays
   
   @param ndofs - number of DOFs
   @param esnrbm - estimated number of rigid body motions
   
   JK, 20.6.2006
*/
void flsubdom::initialization (long ndofs,long esnrbm,double *rhs)
{
  long i,j;
  double *ee,*gggg;

  //  number of displacements
  ndof = ndofs;
  //  estimated number of rigid body motions
  enrbm = esnrbm;
  
  //  allocation of array for rigid body motions
  if (rbm==NULL)
    rbm = new double [ndof*enrbm];
  //  allocation of array for list of dependent rows and columns
  if (de==NULL)
    de = new long [enrbm];
  
  
  //  computation of rigid body motions and factorization of the matrix
  Smat->kernel (rbm,nrbm,de,enrbm,Mp->limit,3);
  
  
  //  number of Lagrange multipliers
  nlm=Gtm->flsub.number_of_lagr_mult (Gtm->nln,Gtm->lgnodes,Gtm->gnodes);
  
  if (w==NULL)
    w = new double [nlm];
  
  if (compli==NULL){
    compli = new double [nlm];
    for (i=0;i<nlm;i++){
      compli[i]=5.0;
    }
  }
  
  if (Mp->tprob == nonlin_floating_subdomain){
    
    if (tw==NULL)
      tw = new double [nlm];
    
    /*
    if (mu==NULL){
      mu = new double [nlm];
      for (i=0;i<nlm;i++)
	mu[i]=0.0;
    }
    if (muo==NULL){
      muo = new double [nlm];
      for (i=0;i<nlm;i++)
	muo[i]=0.0;
    }
    if (dmu==NULL){
      dmu = new double [nlm];
      for (i=0;i<nlm;i++)
	dmu[i]=0.0;
    }
    if (ddmu==NULL){
      ddmu = new double [nlm];
      for (i=0;i<nlm;i++)
	ddmu[i]=0.0;
    }
    */
  }
  
  
  //  list of displacements and multipliers
  Gtm->flsub.displ_extract (Gtm->nln,Gtm->lgnodes,Gtm->gnodes);
  
  //  allocation of vector e
  if (e==NULL)
    e = new double [nrbm];
  nullv (e,nrbm);
  
  //  computation of vector e
  evector (rhs);
  
  
  //  allocation of array for matrix G
  if (gg==NULL)
    gg = new double [nlm*nrbm];
  nullv (gg,nlm*nrbm);
  
  //  assembling of G matrix
  g_matrix ();
  
  // *********************************************
  //  computation of inverse matrix to matrix GG
  // *********************************************
  
  //  allocation of array for inverse matrix to matrix GG
  if (igg==NULL)
    igg = new double [nrbm*nrbm];
  
  //  auxiliary arrays
  gggg = new double [nrbm*nrbm];
  ee = new double [nrbm*nrbm];
  
  //  unit (identity) matrix
  for (i=0;i<nrbm;i++){
    for (j=0;j<nrbm;j++){
      ee[i*nrbm+j]=0.0;
    }
    ee[i*nrbm+i]=1.0;
  }
  
  //  computation of GG
  mtxm (gg,gg,gggg,nlm,nrbm,nrbm);
  
  //  computation of inverse matrix to matrix GG
  gemp (gggg,igg,ee,nrbm,nrbm,Mp->zero,1);
  
  
  delete [] ee;
  delete [] gggg;
  
}


/**
   function assembles %vector e
   notation can be found in Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.
   
   @param f - right hand side (nodal forces)
   
   JK, 20.6.2006 (18.11.2005)
*/
void flsubdom::evector (double *f)
{
  long i;
  
  for (i=0;i<nrbm;i++){
    e[i]=0.0-ss(rbm+i*ndof,f,ndof);
  }
  
}

/**
   function assembles %matrix G
   columns are RBM after localization
   notation can be found in Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.
   
   JK, 20.6.2006 (18.11.2005)
*/
void flsubdom::g_matrix (void)
{
  long i,j;
  double *buff;
  
  buff = new double [nlm*nrbm];
  nullv (buff,nlm*nrbm);
  
  //  extraction of RBM from local vectors to buffer
  for (i=0;i<nrbm;i++){
    Gtm->flsub.local_coarse (buff+i*nlm,rbm+i*ndof);
  }
  
  //  assembling of G matrix
  for (i=0;i<nlm;i++){
    for (j=0;j<nrbm;j++){
      gg[i*nrbm+j]=buff[j*nlm+i]*(-1.0);
    }
  }
  
  delete [] buff;
}


/**
   function computes projection in the FETI method
   v_new = v_old - H . (H^T . H)^{-1} . H^T . v_old
   function overwrites %vector v by projected %vector
   notation can be found in Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.

   @param v - %vector
   
   JK, 20.6.2006 (20.11.2005)
*/
void flsubdom::feti_projection (double *v)
{
  long i;
  double *p,*q,*u;
  
  p = new double [nrbm];
  q = new double [nrbm];
  u = new double [nlm];
  
  //  p = G^T.v
  mtxv (gg,v,p,nlm,nrbm);
  
  
  /*
  fprintf (out,"\n\n\n kontrola H^T.v\n");
  for (i=0;i<nrbm;i++){
    fprintf (out," %lf",p[i]);
  }
  */

  
  //  q = (G^T.G)^{-1}.p
  mxv (igg,p,q,nrbm,nrbm);


  /*
  fprintf (out,"\n\n\n kontrola (H^T.H)^{-1}.H^T.v\n");
  for (i=0;i<nrbm;i++){
    fprintf (out," %lf",q[i]);
  }
  */


  //  w = G.q
  mxv (gg,q,u,nlm,nrbm);
  
  
  /*
  fprintf (out,"\n\n\n kontrola H.(H^T.H)^{-1}.H^T.v\n");
  for (i=0;i<ndofcp;i++){
    fprintf (out," %lf",w[i]);
  }
  */
  
  //  v_new = v_old - w
  for (i=0;i<nlm;i++){
    v[i]-=u[i];
  }
  
  delete [] u;  delete [] q;  delete [] p;
}



/**
   function performs preconditioned modified conjugate gradient method
   for more details see Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.

   @param rhs - right hand side (%vector of nodal forces)
   @param out - pointer to output stream
   
   JK, 20.6.2006 (20.11.2005)
*/
void flsubdom::pmcg (double *rhs,FILE *out)
{
  long i,j,nicg,anicg;
  double nom,denom,alpha,beta,zero,errcg,aerrcg;
  double *s,*ls,*r,*p,*pp,*igge;
  //  zaloha smerovych vektoru pro plnou ortogonalizaci
  
  //zero = Mp->limit;
  zero = 1.0e-20;
  
  //  maximum number of iterations
  nicg = Mp->ssle.ni;
  //  required error
  errcg = Mp->ssle.err;
  
  /*
  //  array for compliance coefficients
  //  compli = new double [nlm];
  j=0;
  for (i=0;i<nlm/2;i++){
    compli[j]=0.0;
    j++;
    compli[j]=0.0;
    j++;
  }
  */
  
  //  direction vector allocation
  s = new double [nlm];
  //zald = new double [nlm*50];
  //zalp = new double [nlm*50];
  //jkkj = new double [50];
  //  residuum vector allocation
  r = new double [nlm];
  //  auxiliary vector allocation
  p = new double [nlm];
  
  ls = new double [ndof];
  pp = new double [ndof];
  
  // ************************************
  //  initiation of conjugate gradients
  // ************************************

  //  computation of vector e
  evector (rhs);

  //  auxiliary array
  igge = new double [nrbm];
  //  (G^T.G)^{-1} e
  mxv (igg,e,igge,nrbm,nrbm);
  //  lambda_0 = G [(G^T.G)^{-1} e]
  mxv (gg,igge,w,nlm,nrbm);
  delete [] igge;
  
  //  E^T lamba_0
  nullv (ls,ndof);
  Gtm->flsub.coarse_local (w,ls);
  
  //  f - E^T lambda_0
  //  terms are changed in order to not rewrite nodal forces
  //  therefore multiplication by -1 is necessary
  subv (ls,rhs,ndof);
  cmulv (-1.0,ls,ndof);
  
  //  K^+ (f - E^T \lambda_0) = p
  Smat->ldl_feti (pp,ls,nrbm,de,zero);
  
  //  r = E p
  nullv (r,nlm);
  Gtm->flsub.local_coarse (r,pp);
  
  
  //  doplnek na diagonalu
  for (j=0;j<nlm;j++){
    r[j]-=compli[j]*w[j];
    //g[j]-=compli[j];
  }
  //  konec doplnku na diagonale
  
  //for (i=0;i<nlm;i++){
  //zalp[i]=g[i];
  //}
  
  //  P r
  feti_projection (r);
  
  //  initialization of direction vector
  for (i=0;i<nlm;i++){
    s[i]=r[i];
    //zald[i]=d[i];
  }
  
  //  nominator evaluation
  nom = ss (r,r,nlm);
  
  // *************************
  //  main iteration loop   **
  // *************************
  for (i=0;i<nicg;i++){
    
    // ******************************************
    //   auxiliary vector computation F.s = p  **
    // ******************************************