adaptivity.cpp

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

   Function computes error of FEM solution from difference between
   apriximate derivative fied(from FEM solution) and much more exact derivative fied(from array rderfull).
   Considering 'printflags' printing of any files and information is performed.
   
   created  3.6.2002, Ladislav Svoboda, [e-mail]
*/
void adaptivity::nonlin_error (void)
{
  long i;
  double e2     ,u2     ,volume   ,corr1,corr2;
  vector ei2(ne),ui2(ne),sizel(ne),refsizel(ne),epsil(ne);
  
  u2 = e2 = volume = 0.0;
  
  for (i=0;i<ne;i++){
    
    switch (Mt->give_elem_type (i)){
    case planeelementlt:{  Pelt->compute_error (i,volume,ei2[i],ui2[i],sizel[i],rderfull,otherflags);  corr1 = 1.00; corr2 = 1.11;   break; }
    case planeelementlq:{  Pelq->compute_error (i,volume,ei2[i],ui2[i],sizel[i],rderfull);             corr1 = 1.00; corr2 = 1.20;   break; }
    case planeelementqt:{  Peqt->compute_error (i,volume,ei2[i],ui2[i],sizel[i],rderfull);             corr1 = 1.07; corr2 = 1.55;   break; }
    case planeelementqq:{  Peqq->compute_error (i,volume,ei2[i],ui2[i],sizel[i],rderfull);             corr1 = 1.00; corr2 = 1.60;   break; }
    case lineartet:{       Ltet->compute_error (i,volume,ei2[i],ui2[i],sizel[i],rderfull);             corr1 = 1.00; corr2 = 1.00;   break; }
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function");
      fprintf (stderr," adaptivity::compute_error (%s, line %d)",__FILE__,__LINE__);
    }
    }
    
    if (tad == 1) ei2[i] *= corr1*corr1;
    if (tad == 2) ei2[i] *= corr2*corr2;
    
    e2 += ei2[i];
    u2 += ui2[i];
  }
  
  compute_refsizel_nonlin (refsizel.a,sizel.a,ei2.a,u2,e2,ord,ne);
  
  /*
  compute_epsil_nonlin (ei2.a,ui2.a,ne,epsil.a,e2,u2);
  
  //  double corr = 0.8;  //<0.7;1.0>  test = 0.8
  
  for (i=0;i<ne;i++){
    if (epsil[i]<0.05) epsil[i] = 0.05;
    refsizel[i] = sizel[i] / pow(epsil[i] , 1.0/(double)ord * corr);
  }
  */
  // ***********************************************
  if (answer){    
    char file[255];
    sprintf (file,"%srefsizel.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"refsizel",refsizel.a,'d');
  }
  // ***********************************************
  
  print (volume,e2,u2,ei2.a,ui2.a,refsizel.a);
}

/**
   
   
 */
void adaptivity::compute_refsizel_nonlin (double *refsizel,double *sizel,double *ei2,double u2,double e2,long ord,long ne)
{
  long i,ri,ci,ipp,nle;
  double em2,ff,limf,specacc;
  double *rsemin,*rsemax,*f,*epsil;
  char file[255];
  vector q,sig;
  matrix sigt(3,3);
  
  ri = ci = 0;
  specacc = Mp->adapt_accuracy;
  limf = 0.1;
  
  f = new double [ne];
  epsil = new double [ne];
  rsemin = new double [ne];
  rsemax = new double [ne];
  
  allocv (4,sig);  // only 2D planestress
  allocv (1,q);    // only 2D planestress
  
  if (sqrt(e2/(u2+e2)) > specacc)
    answer = 1;
  
  // epsil
  em2 = specacc*specacc*(u2+e2)/ne;
  
  for (i=0;i<ne;i++){
    epsil[i] = sqrt(ei2[i]/em2);
    if (epsil[i]<0.05) epsil[i] = 0.05;
    rsemax[i] = sizel[i]/pow(epsil[i],1.0/(double)ord);
  }
  
  for (i=0;i<ne;i++){
    epsil[i] = epsil[i]/0.3;
    if (epsil[i]<0.05) epsil[i] = 0.05;
    rsemin[i] = sizel[i]/pow(epsil[i],1.0/(double)ord);
  }
  
  // f
  vector_tensor (sig,sigt,Mm->ip[0].ssst,stress);
  ff = Mm->yieldfunction (0,0,sigt,q);
  
  for (i=0;i<ne;i++){
    ipp = Mt->elements[i].ipp[ri][ci];
    
    Mm->givestress (0,ipp,sig);
    sig[3] = 0.0;
    
    switch (Mm->ip[ipp].tm[0]){
    case jflow:{  q[0] = Mm->ip[ipp].eqother[Mm->ip[ipp].ncompstr+1];  break; }
    default:   {  fprintf (stderr,"\n\n unknown material type is required in function adaptivity::compute_refsizel_nonlin (%s, line %d).\n",__FILE__,__LINE__); }
    }
    
    vector_tensor (sig,sigt,Mm->ip[ipp].ssst,stress);
    f[i] = Mm->yieldfunction (ipp,0,sigt,q);
  }
  
  if (0){
  // ************************************************************
  // print ff.dat
  FILE *fx; sprintf (file,"%sfxlambda.dat",Mp->path);
  fx = fopen (file,"a"); if (fx==NULL) fprintf (stderr,"\n fxlambda file has not been opened.\n\n");
  fprintf (fx,"  %f\n",f[85]); fclose (fx);
  
  // print fo f.dx
  sprintf (file,"%sf.dx",Mp->path);
  print_valel (Gtm,Mp,file,"f",f,'d');
  
  // print sts  
  double *sts;
  sts = new double [ne];
  
  for (i=0;i<ne;i++){
    if (Mm->ip[i].eqother[4] != 0.0)    sts[i] = 1.0;
    else if (f[i] < limf*ff)          sts[i] = 3.0;
    else                              sts[i] = 2.0;
  }
  sprintf (file,"%ssts.dx",Mp->path);
  print_valel (Gtm,Mp,file,"sts",sts,'d');
  
  delete [] sts;
  // ************************************************************
  }
  
  nle = 0;
  for (i=0;i<ne;i++){
    
    if (Mm->ip[i].eqother[4] != 0.0)
      refsizel[i] = rsemin[i];  //20.0*rsemax[i];
    else if (f[i] < limf*ff)
      refsizel[i] = rsemax[i];
    else{
      refsizel[i] = rsemin[i] + (rsemax[i] - rsemin[i]) * f[i] / ff / limf ;
      nle++;
    }
  }
  
  if (nle)
    answer = 1;
  
  // print fo f.dx
  sprintf (file,"%srefsizel.dx",Mp->path);
  print_valel (Gtm,Mp,file,"refsizel",refsizel,'d');
  
  delete [] f;
  delete [] epsil;
  delete [] rsemin;
  delete [] rsemax;
}




/**
   
 */
void adaptivity::compute_epsil_nonlin ( double *ei2, double *ui2,long ne,double *epsil,double &e2,double &u2)
{
  long i,ri,ci,ipp,nelel,nclosel,nresel;
  double f0,limf,em2,um2,sael,sapl;
  vector q,sig,beta;
  char file[255];
  matrix sigt(3,3);
  strastrestate ssst;
  
  e2 = u2 = nelel = nclosel = nresel = ri = ci = 0;
  limf = 0.25;
  sael = Mp->adapt_accuracy;   // elastic accuracy
  sapl = sael / plra;           // plastic accuracy
  
  // predpoklad jednoho typu ssst na domene
  ssst = Mm->ip[0].ssst;
  
  allocv (ne,beta);
  // zde to predelat aby se to alocovalo automaticky, ted je tam nastvrdo 2D planestress
  allocv (4,sig);
  allocv (1,q);
  
  vector_tensor (sig,sigt,ssst,stress);  // zde predpoklad ze sig a g jsou vynulovane !!
  f0 = Mm->yieldfunction (ipp,0,sigt,q);   // !!!!!!! proc neni ipp vynulovane ??????
  
  for (i=0;i<ne;i++){
    ipp = Mt->elements[i].ipp[ri][ci];
    
    if ( Mm->j2f[Mm->ip[ipp].idm[0]].give_consparam (ipp,0) ) {  beta[i] = 0.0;  continue;  }
    
    Mm->givestress (0,ipp,sig);
    sig[3] = 0.0;
    
    switch (Mm->ip[ipp].tm[0]){
    case jflow:{  q[0] = Mm->ip[ipp].eqother[Mm->ip[ipp].ncompstr+1];  break; }
    default:   {  fprintf (stderr,"\n\n unknown material type is required in function adaptivity::compute_epsil_nonlin (%s, line %d).\n",__FILE__,__LINE__); }
    }
    
    vector_tensor (sig,sigt,ssst,stress);
    beta[i] = Mm->yieldfunction(ipp,0,sigt,q) / f0;
    
    // ***********************************************
    if (i == 85){
      FILE *fx; sprintf (file,"%sfxlambda.dat",Mp->path);
      fx = fopen (file,"a"); if (fx==NULL) fprintf (stderr,"\n fxlambda file has not been opened.\n\n");
      fprintf (fx,"  %f\n",beta[i]*f0); fclose (fx);
    }
    // ***********************************************
    
    e2 += ei2[i];
    u2 += ui2[i];
    nelel++;
  }
  u2 += e2;
  
  um2 = u2/nelel;
  em2 = sael*sael*um2;
  
  for (i=0;i<ne;i++){
    //if      (beta[i] == 0.0)  epsil[i] = 1.0;
    if      (beta[i] == 0.0)  epsil[i] = sqrt(ei2[i] / em2 );
    else if (beta[i] > limf)  epsil[i] = sqrt(ei2[i] / em2 );
    else {
      epsil[i] = sqrt(ei2[i] / em2 );
      epsil[i] = sqrt(ei2[i] / ui2[i] ) / (sapl+(sael-sapl)*beta[i]/limf);
      epsil[i] = sqrt(ei2[i] / um2 ) / (sapl+(sael-sapl)*beta[i]/limf);
      nclosel++;
      if (epsil[i]>1.0) nresel++;
      //if (epsil[i]>2.0)   answer = 1;
    }
  }
  
  fprintf (stdout,"\n nclosel = %ld",nclosel);
  fprintf (stdout,"\n nresel = %ld",nresel);
  
  if (nclosel){
    if (nresel > clos*nclosel)
      answer = 1;
  }
  else
    if ( sael < sqrt(e2/u2) )
      answer = 1;
  
  // ***********************************************
  //if (answer){
  if (answer){
    
    sprintf (file,"%sbeta.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"beta",beta.a,'d');
    
    for (i=0;i<ne;i++){
      if      (beta[i] == 0.0)    beta[i] = 1.0;
      else if (beta[i] < limf)    beta[i] = 2.0;
      else                        beta[i] = 3.0;
    }
    sprintf (file,"%ssts.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"sts",beta.a,'d');
    
    for (i=0;i<ne;i++)  beta[i] = sqrt (ui2[i]);
    sprintf (file,"%sui.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"ui",beta.a,'d');
  
    for (i=0;i<ne;i++)  beta[i] = sqrt (ei2[i]);
    sprintf (file,"%sei.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"ei",beta.a,'d');
    
    for (i=0;i<ne;i++)  beta[i] = 100 * sqrt (ei2[i] / (ui2[i]+ei2[i]));
    sprintf (file,"%serror.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"error",beta.a,'d');

    sprintf (file,"%sepsil.dx.%s",Mp->path,ni);
    print_valel (Gtm,Mp,file,"epsil",epsil,'d');
    
  }
  
  // ***********************************************  
}


/**
 */
void adaptivity::print (double volume,double e2,double u2, double *ei2, double *ui2, double *refsizel)
{
  long i;
  double e_pere;
  vector valel(ne);
  char name[255];
  
  e_pere = 100.0*sqrt(e2/u2);
  
  // *** output to file.out ***
  fprintf (Out,"\n\n\n ADAPTIVITY\n\n e2 = %f\n u2 = %f\n e_pere = %f",e2,u2,e_pere);
  
  // *** standard output ***
  if (printflags & 1)
    fprintf (stdout,"\n e2celk = %20.10f\n u2celk = %20.10f\n e_pere = %f %%\n volume = %f",e2,u2,e_pere,volume);
  
  // *** print of background file ***
  if (printflags & 2  &&  answer){
    sprintf (name,"%smesh.bgm%s",Mp->path,ni);
    print_valel (Gtm,Mp,name,"chyba",refsizel,'e');
  }
  
  // *** print of contoures ***
  if (printflags & 4){
    for (i=0;i<ne;i++){
      valel[i] = 100.0*sqrt(ei2[i]/(ui2[i]+ei2[i]));
      //valel[i] = pow(sqrt(ei2[i]/em2),1.0/(double)ord);
    }
    //sprintf (name,"%schyba.el%s",Mp->path,ni);
    //print_valel (Gtm,Mp,name,"chyba",valel.a,'e');
    sprintf (name,"%schyba.dx%s",Mp->path,ni);
    print_valel (Gtm,Mp,name,"chyba",valel.a,'d');
  }
  
  // *** fine problem ***
  if (printflags & 8)
    print_coarsemesh (ei2);
  
  // *** print of test file ***
  if (printflags & 16)
    print_test (refsizel);
  
  // *** print of dx file ***
  if (printflags & 32  &&  answer){
    sprintf (name,"%s%s.ada.dx%s",Mp->path,Mp->filename,ni);
    print_default_dx (Gtm,Mp,Mt,Mm,0,name);
  }

  // *** print of rderfull in file.dx ***
  if (printflags & 64){
    sprintf (name,"%srderfull.1.dx%s",Mp->path,ni);
    print_valnod (Gtm,Mp,Mt,name,"rderfull",rderfull,'d');
    

    if (otherflags & 8){
      sprintf (name,"%srderfull.2.dx%s",Mp->path,ni);
      print_valnod (Gtm,Mp,Mt,name,"rderfull",rderfull + ncomp*nn,'d');
    }
  }
}


/**

   created  3.6.2002, Ladislav Svoboda, [e-mail]
*/
void adaptivity::print_test ( double *refsizel)
{
  long i,j,ipp;
  char f_name[255];
  FILE *f_test;
  
  sprintf (f_name,"%s%s.test",Mp->path,Mp->filename);
  f_test = fopen (f_name,"a");            if (f_test==NULL) fprintf (stderr,"\n .test file has not been opened.\n\n");
  
  
  // ***************************************************************************  
  fprintf (f_test,"\n\n LEFT HAND SIDE \n\n");
  
  for (i=0;i<Mb->nlc*Ndofm;i++){
    fprintf (f_test," %20.10f ",Lsrs->lhs[i]);
    if ((i+1)%7==0)  fprintf (f_test,"\n");
  }
  
  
  // ***************************************************************************  
  fprintf (f_test,"\n\n STRAINS AND STRESSES ON INTEGRATION POINTS \n\n");
  
  for (i=0;i<ne;i++){
    ipp = Mt->elements[i].ipp[0][0];
    
    switch ((long)Mt->elements[i].te){
    case planeelementlq:{
      for (j=0;j<4;j++){
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp+j].strain[0]);
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp+j].strain[1]);
	fprintf (f_test," %20.10f \n",Mm->ip[ipp+4].strain[2]);
      }
      break;
    }
    case planeelementqq:{
      for (j=0;j<4;j++){
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp+j  ].strain[0]);
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp+j  ].strain[1]);
	fprintf (f_test," %20.10f \n",Mm->ip[ipp+j+4].strain[2]);
	}
      break;