flsubdom.cpp

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

    
    //  E^T s
    nullv (ls,ndof);
    Gtm->flsub.coarse_local (s,ls);
    
    //  K^+ E^T s = p
    Smat->ldl_feti (pp,ls,nrbm,de,zero);
      
    //  E p
    nullv (p,nlm);
    Gtm->flsub.local_coarse (p,pp);
    
    //for (j=0;j<nlm;j++){
    //zalp[i*nlm+j]=p[j];
    //}
    

    //  doplnek na diagonalu
    for (j=0;j<nlm;j++){
      p[j]+=compli[j]*s[j];
    }
    //  konec doplnku na diagonalu


    //  denominator of alpha
    denom = ss (s,p,nlm);
    
    //fprintf (stdout,"\n denominator u alpha   %le",denom);
    
    fprintf (out,"\n\n krok %ld      kontrola citatele a jmenovatele pred alpha %e / %e",i,nom,denom);
    //fprintf (stderr,"\n\n kontrola citatele a jmenovatele pred alpha %lf / %lf",nom,denom);
    
    if (fabs(denom)<zero){
      fprintf (stderr,"\n\n zero denominator in modified conjugate gradient method (file %s, line %d).\n",__FILE__,__LINE__);
      break;
    }
    
    // *******************
    //   vypocet alpha  **
    // *******************
    alpha = nom/denom;
    
    // **************************************************************
    //  vypocet noveho gradientu g a nove aproximace neznamych x   **
    // **************************************************************
    for (j=0;j<nlm;j++){
      w[j]+=alpha*s[j];
      r[j]-=alpha*p[j];
    }
    
    feti_projection (r);
    
    
    denom = nom;
    if (fabs(denom)<zero){
      fprintf (stderr,"\n\n zero denominator of beta in function (file %s, line %d).\n",__FILE__,__LINE__);
      break;
    }
    
    nom = ss(r,r,nlm);
    
    fprintf (stdout,"\n CG iteration  %4ld        norm g   %e",i,sqrt(nom));
    fprintf (out,"\n CG iteration  %4ld        norm g   %e",i,sqrt(nom));
    
    fprintf (out,"\n kontrola citatele a jmenovatele pred betou  %e / %e",nom,denom);
    //fprintf (stderr,"\n kontrola citatele a jmenovatele pred betou  %lf / %lf",nom,denom);
    
    if (sqrt(nom)<errcg){
      break;
    }
    
    //  computation of beta coefficient
    beta = nom/denom;
    
    //  new direction vector
    for (j=0;j<nlm;j++){
    s[j]=beta*s[j]+r[j];
    }
    
    /*
    for (j=0;j<i;j++){
      cit=0.0;  jmen=0.0;
      for (k=0;k<nlm;k++){
	cit+=g[k]*zalp[j*nlm+k];
	jmen+=zald[j*nlm+k]*zalp[j*nlm+k];
      }
      jkkj[j]=cit/jmen;
    }
    for (j=0;j<nlm;j++){
      d[j]=0.0;
      for (k=0;k<i;k++){
	d[j]-=jkkj[k]*zald[k*nlm+j];
      }
      d[j]-=g[j];
      zald[i*nlm+j]=d[j];
    }
    */

    anicg=i;  aerrcg=nom;
    
    //fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
    //for (j=0;j<nlm;j++){
    //fprintf (out,"\n lagr. mult %4ld    %e",j,w[j]);
    //}
    
  }
  
  fprintf (out,"\n\n\n\n kontrola Lagrangeovych multiplikatoru \n");
  for (j=0;j<nlm;j++){
    fprintf (out,"\n lagr. mult %4ld    %e",j,w[j]);
  }
  
  long nlgn = Gtm->nln;
  long nod;
  fprintf (out,"\n\nprubeh napeti \n\n");
  double sum=0.0;
  for (i=0;i<nlgn;i++){
    nod = Gtm->lgnodes[i].nodes[0];
    fprintf (out,"%lf %le\n",Gtm->gnodes[nod].y,w[2*i+1]);
    sum+=w[2*i+1];
  }
  fprintf (out,"\n\n soucet multiplikatoru %lf\n",sum);
  fprintf (out,"\n\n");




  Mp->ssle.ani=anicg;
  Mp->ssle.aerr=aerrcg;
  
  delete [] s;
  delete [] p;
  delete [] r;
  
  delete [] ls;
  delete [] pp;
}



/**
   function computes displacements from Lagrange multipliers
   function is used in FETI method
   for more details see Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.
   
   @param d - array containing displacements
   @param f - array containing load forces

   JK, 20.6.2006 (21.11.2005)
*/
void flsubdom::lagrmultdispl (double *d,double *f)
{
  long i,j;
  double *r,*pp,*gamma,*u;
  
  /*
  //  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++;
  }
  */

  pp = new double [ndof];
  //  residuum vector
  r = new double [nlm];
  
  //  E^T lambda
  nullv (pp,ndof);
  Gtm->flsub.coarse_local (w,pp);
  
  //  f - E^T lambda
  //  terms are changed in order to not rewrite nodal forces
  //  therefore multiplication by -1 is necessary
  subv (pp,f,Ndofm);
  cmulv (-1.0,pp,ndof);
  
  
  //  K^+ (f - E^T lambda) = d
  Smat->ldl_feti (d,pp,nrbm,de,Mp->zero);
  
  
  //  E d
  nullv (r,nlm);
  Gtm->flsub.local_coarse (r,d);
  
  
  gamma = new double [nrbm];
  u = new double [nrbm];
  
  for (i=0;i<nlm;i++){
    r[i]-=compli[i]*w[i];
  }
  
  //  G^T r
  mtxv (gg,r,u,nlm,nrbm);
  //  (G^T G)^{-1} G^T r
  mxv (igg,u,gamma,nrbm,nrbm);

  
  for (i=0;i<ndof;i++){
    for (j=0;j<nrbm;j++){
      d[i]+=rbm[j*ndof+i]*gamma[j];
    }
  }
  
  
  delete [] pp;
  delete [] r;
  delete [] gamma;
  delete [] u;
}


/**
   function computes displacements from Lagrange multipliers
   function is used in FETI method
   for more details see Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.
   
   @param d - array containing displacements
   @param f - array containing load forces

   JK, 20.6.2006 (21.11.2005)
*/
void flsubdom::nonlinlagrmultdispl (double *d,double *f)
{
  long i,j;
  double *r,*pp,*gamma,*u;
  
  /*
  //  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++;
  }
  */

  pp = new double [ndof];
  //  residuum vector
  r = new double [nlm];
  
  //  E^T lambda
  nullv (pp,ndof);
  Gtm->flsub.coarse_local (tw,pp);
  
  //  f - E^T lambda
  //  terms are changed in order to not rewrite nodal forces
  //  therefore multiplication by -1 is necessary
  subv (pp,f,Ndofm);
  cmulv (-1.0,pp,ndof);
  
  
  //  K^+ (f - E^T lambda) = d
  Smat->ldl_feti (d,pp,nrbm,de,Mp->zero);
  
  
  //  E d
  nullv (r,nlm);
  Gtm->flsub.local_coarse (r,d);
  
  
  gamma = new double [nrbm];
  u = new double [nrbm];
  
  for (i=0;i<nlm;i++){
    r[i]-=compli[i]*tw[i];
  }
  
  //  G^T r
  mtxv (gg,r,u,nlm,nrbm);
  //  (G^T G)^{-1} G^T r
  mxv (igg,u,gamma,nrbm,nrbm);

  
  for (i=0;i<ndof;i++){
    for (j=0;j<nrbm;j++){
      d[i]+=rbm[j*ndof+i]*gamma[j];
    }
  }
  
  
  delete [] pp;
  delete [] r;
  delete [] gamma;
  delete [] u;
}


/**
   function solves system of linear algebraic equations
   by modified conjugate gradient method
   Lagrange multipliers are computed first and then displacements are obtained
   for more details see Kruis: Domain Decomposition Methods for
   Distributed Computing. Saxe-Coburg, 2006.
   
   @param lhs - left hand side %vector (%vector of unknown nodal displacement)
   @param rhs - right hand side (%vector of nodal forces)
   
   JK, 20.6.2006
*/
void flsubdom::solve_lin_alg_system (double *lhs,double *rhs)
{
  //  preconditioned modified conjugate gradient method
  pmcg (rhs,Out);
  
  //  computation of displacements from Lagrange multipliers
  lagrmultdispl (lhs,rhs);
  
}

/**
   function adds increments of Lagrange multipliers to
   the %vector of total multipliers in nonlinear problems
   
   @param dlambda - increment from arclength method
   
   JK, 22.6.2006
*/
void flsubdom::add_mult (double dlambda)
{
  long i;
  /*
  for (i=0;i<nlm;i++){
    tw[i]+=dlambda*w[i];
  }
  */
  for (i=0;i<nlm;i++){
    w[i]*=dlambda;
    tw[i]+=w[i];
  }
}
/**
   function corrects Lagrange multipliers with respect to
   defined law
   
   JK, 24.6.2006
*/
void flsubdom::mult_correction ()
{
  long i;
  
  //cv = new double [nlm];
  //Gtm->flsub.local_coarse (cv,Lsrs->lhs);
  
  nch=0;
  for (i=0;i<nlm;i++){
    
    if (tw[i]>1.0 && tw[i]<-1.0){
      compli[i]=5.0;
    }


    /*
    if (mu[i]>100.0){
      //tw[i]=10.0;
      //compli[i]=20000.0;
      //compli[i]=2.0;
      compli[i]=(mu[i]-100.0)/10000.0;
      nch++;
    }
    */
    /*
    if (tw[i]>10.0 && cv[i]<10.0){
      //tw[i]=10.0;
      //compli[i]=20000.0;
      compli[i]=2.0;
      //compli[i]=tw[i]-10.0;
      nch++;
    }
    if (tw[i]>10.0 && cv[i]>10.0){
      //tw[i]=10.0;
      compli[i]=20000.0;
      //compli[i]=2.0;
      //compli[i]=tw[i]-10.0;
      nch++;
    }
    */
    /*
    if (tw[i]>100.0)
      tw[i]=100.0;
    */
  }
  
  //delete [] cv;
}

void flsubdom::factorize ()
{
  Smat->kernel (rbm,nrbm,de,enrbm,Mp->limit,3);
}