adaptivity.cpp

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

    }
    case planeelementlt:{
      fprintf (f_test," %20.10f "  ,Mm->ip[ipp].strain[0]);
      fprintf (f_test," %20.10f "  ,Mm->ip[ipp].strain[1]);
      fprintf (f_test," %20.10f \n",Mm->ip[ipp].strain[2]);
      break;
    }
    case planeelementqt:
    case planeelementsubqt:{
      for (j=0;j<3;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+3].strain[2]);
      }
      break;
    }
    case lineartet:{
      for (j=0;j<6;j++)
	fprintf (f_test," %20.10f ",Mm->ip[ipp].strain[j]);
      
      fprintf (f_test,"\n");
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function");
      fprintf (stderr," adaptivity::print_test  (%s, line %d)",__FILE__,__LINE__);
    }
    }
    
  }
  
  if (otherflags & 8){
    fprintf (f_test,"\n");
    
    for (i=0;i<ne;i++){
      ipp = Mt->elements[i].ipp[0][0];
      
      switch ((long)Mt->elements[i].te){
      case planeelementlt:{
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp].stress[0]);
	fprintf (f_test," %20.10f "  ,Mm->ip[ipp].stress[1]);
	fprintf (f_test," %20.10f \n",Mm->ip[ipp].stress[2]);
	break;
      }
      default:{
	fprintf (stderr,"\n\n unknown element type is required in function");
	fprintf (stderr," adaptivity::print_test  (%s, line %d)",__FILE__,__LINE__);
      }
      }
    }
  }
  
  
  // ***************************************************************************  
  fprintf (f_test,"\n\n RDERFULL \n\n");
  
  for (i=0;i<ncomp*nn;i++){
    fprintf (f_test," %20.10f ",rderfull[i]);
    if ((i+1)%7==0)  fprintf (f_test,"\n");
  }
  
  if (otherflags & 8){
    fprintf (f_test,"\n");
    
    for (i=0;i<ncomp*nn;i++){
      fprintf (f_test," %20.10f ",rderfull[i + ncomp*nn]);
      if ((i+1)%7==0)  fprintf (f_test,"\n");
    }
  }
  
  // ***************************************************************************  
  fprintf (f_test,"\n\n REFSIZEL \n\n");
  
  for (i=0;i<Mt->ne;i++){
    fprintf (f_test," %20.10f ",refsizel[i]);
    if ((i+1)%7==0)  fprintf (f_test,"\n");
  }
  
  // ***************************************************************************  
  if (Mm->ip[0].eqother!=NULL){
    fprintf (f_test,"\n\n OTHER \n\n");
    
    for (i=0;i<Mm->tnip;i++){
      for (j=0;j<6;j++)
	fprintf (f_test," %20.10f "  ,Mm->ip[i].eqother[j]);
      
      fprintf (f_test,"\n");
    }
  }
  
  fclose (f_test);
}


//**********************
//***  FINE PROBLEM  ***
//**********************

void adaptivity::compute_error_fine (void)
{
  long i,j,auxl,te,necm,nncm,nnelcm,nesure;
  long *nadjelnod,**adjelnod;
  long *loccfmesh,**etnodes;
  double **xyrr,thetaglob,e2cm,e2cmfine;
  vector ei2(ne),ei2cm,ei2cmfine,theta,pom;
  char auxfile[1023];
  FILE *cmc,*edg,*cor;


  if (Pelt==NULL)  Pelt = new planeelemlt;
  if (Peqt==NULL)  Peqt = new planeelemqt;
  if (Pesqt==NULL)  Pesqt = new planeelemsubqt;
  if (Pelq==NULL)  Pelq = new planeelemlq;
  if (Peqq==NULL)  Peqq = new planeelemqq;

  // reading of edges.in
  sprintf (auxfile,"%sedges.in%s",Mp->path,ni);
  edg=fopen(auxfile,"r"); if (edg==NULL){ fprintf (stderr,"\n edg file has not been specified.\n Try it again!\n"); return ;}

  fscanf (edg,"%ld",&nesure);
  if (nesure!=ne) fprintf (stdout,"\n\n\n nesure!=ne \n\n\n ");
  loccfmesh = new long [ne];

  for (i=0;i<ne;i++){
    fscanf (edg,"%ld",&auxl);
    fscanf (edg,"%ld",loccfmesh+auxl-1);
  }
  fclose(edg);

  // reading of cmchar.in
  sprintf (auxfile,"%sproblem/cmeshchr%s",Mp->path,ni);
  cmc=fopen(auxfile,"r"); if (cmc==NULL){fprintf (stderr,"\n cmc file has not been specified.\n Try it again!\n"); return ;}

  fscanf (cmc,"%ld",&nncm);
  nadjelnod = new long [nncm];
  adjelnod = new long* [nncm];
  xyrr = new double* [nncm];       for (i=0;i<nncm;i++)    xyrr[i] = new double [7];

  for (i=0;i<nncm;i++){
    for (j=0;j<7;j++)
      fscanf (cmc,"%lf",xyrr[i]+j);
    fscanf (cmc,"%ld",nadjelnod+i);
    adjelnod[i] = new long [nadjelnod[i]];
    for (j=0;j<nadjelnod[i];j++)
      fscanf (cmc,"%ld",adjelnod[i]+j);
  }

  fscanf (cmc,"%ld",&necm);
  etnodes = new long* [necm];
  allocv (necm,ei2cm);

  for (i=0;i<necm;i++){
    fscanf (cmc,"%ld %ld",&te,&nnelcm);

    etnodes[i] = new long [nnelcm+2];
    etnodes[i][0] = te;
    etnodes[i][1] = nnelcm;
    for (j=2;j<nnelcm+2;j++)
      fscanf (cmc,"%ld",etnodes[i]+j);
    fscanf (cmc,"%lf",ei2cm.a+i);
  }
  fclose(cmc);

  // computing
  for (i=0;i<ne;i++){
    te = Mt->give_elem_type (i);
    
    switch (te){
    case planeelementqt:{
      Peqt->compute_error_fine (i,etnodes[loccfmesh[i]-1],xyrr,rderfull,ei2[i]);
      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function");
      fprintf (stderr," adaptivity::compute_error_fine (%s, line %d)",__FILE__,__LINE__);
    }
    }
  }   //fprintf_vector (stdout,ei2,"ei2:",10);
  
  
  allocv (necm,ei2cmfine);
  allocv (necm,theta);
  
  for (i=0;i<ne;i++)
    ei2cmfine[loccfmesh[i]-1] += ei2[i];
  
  e2cm = e2cmfine = 0.0;
  for (i=0;i<necm;i++){
    e2cm     += ei2cm    [i];              fprintf (stdout,"\n [%ld] %f",i,ei2cm[i]);
    e2cmfine += ei2cmfine[i];              fprintf (stdout,"   %f",ei2cmfine[i]);
  }                                        //fprintf (stdout,"\n %f",e2cm); fprintf (stdout,"\n %f",e2cmfine);
  fprintf (stdout,"\n e2cm = %f",e2cm);
  fprintf (stdout,"\n e2cmfine = %f",e2cmfine);
  
  // compute theta
  thetaglob = sqrt(e2cm/e2cmfine);
  fprintf (stdout,"\n thetaglob =  %f",thetaglob);
  fprintf (stdout,"\n correction =  %f",1/thetaglob);
  
  for (i=0;i<necm;i++)
    theta[i] = sqrt(ei2cm[i]/ei2cmfine[i]);
  
  for (i=0;i<necm;i++)
    if (theta[i] > 3.0) theta[i] = 1.5;
  
  //fprintf_vector (stdout,theta,"theta",10);
  sprintf (auxfile,"%stheta%s",Mp->path,ni);
  print_thetacon (auxfile,nadjelnod,adjelnod,nncm,necm,xyrr,etnodes,theta.a);
  
  //  print theta in interval (-**;0.65> A <1.5;**)
  fprintf(stdout,"\n ");
  for (j=0,i=0;i<necm;i++)
    if (theta[i]>1.5 ){
      fprintf (stdout,"[%3ld] = %5.3f  ",i,theta[i]);   j++;
      if ((j)%6==0) fprintf(stdout,"\n ");
    }
  fprintf (stdout,"\n j =  %ld",j);
  
  //  print correc_file
  sprintf (auxfile,"%scorrec.txt",Mp->path);
  cor = fopen (auxfile,"a"); if (cor==NULL){fprintf(stderr,"\n cor file has not been specified.\n Try it again!\n"); return ;}
  fprintf (cor,"%.3f\n",1/thetaglob);
  fclose (cor);
  
  
  for (i=0;i<necm;i++){ delete [] etnodes[i];     }
  for (i=0;i<nncm;i++){ delete [] adjelnod[i];
                        delete [] xyrr[i];        }
  delete [] etnodes;
  delete [] adjelnod;
  delete [] xyrr;
  delete [] loccfmesh;
  delete [] nadjelnod;

  /*
  for (i=0;i<necm;i++){
    etnodes [i][0] = 1;
    etnodes [i][1] = etnodes [i][2]+1 ;
    etnodes [i][2] = etnodes [i][3]+1 ;
    etnodes [i][3] = etnodes [i][4]+1 ;
  }
  // print deltasigxcon
  for (i=0;i<nncm;i++){
    for (j=0;j<nn;j++)
      if (Gt->gnodes[j].x == xyrr[i][0] && Gt->gnodes[j].y == xyrr[i][1])
	if ((rderfull[j] - xyrr[i][4]) < 0.0)
	  xyrr[i][3] = -(rderfull[j] - xyrr[i][4]);
	else
	  xyrr[i][3] =   rderfull[0*nn+j] - xyrr[i][4] ;

    xyrr[i][2] = 0.0;
    xyrr[i][3] *= 1000.0;
  }
  print_confile ("fine_problem/deltasigx.",nncm,necm,necm,xyrr,etnodes,ni,dim);
 */
}


void adaptivity::print_thetacon ( char *file, long *nadjelnod,long **adjelnod,long nncm,long necm,double **xyrr,long **etnodes, double *theta)
{
  long i,j,te,nmidp;
  long **auxnod;
  double **auxxyzv;
  vector valnod(nncm);

  auxxyzv = new double* [nncm+necm];
  auxnod = new long* [necm];

  // computing of valnod
  for (i=0;i<nncm;i++){
    for (j=0;j<nadjelnod[i];j++){
      valnod[i] += theta[adjelnod[i][j]]; 
    }
    valnod[i] /= (double)nadjelnod[i]; 
  }

  // assemble aux*
  for (i=0;i<nncm;i++){
    auxxyzv[i] = new double [4];
    auxxyzv[i][0] = xyrr[i][0];
    auxxyzv[i][1] = xyrr[i][1];
    auxxyzv[i][2] = 0.0;
    auxxyzv[i][3] = valnod[i];
  }
  for (i=0;i<necm;i++){
    if (etnodes[i][0] == 23 || etnodes[i][0] == 24){
      auxxyzv[i+nncm] = new double [4];
      auxxyzv[i+nncm][0] = (xyrr[etnodes[i][0+2]][0] + xyrr[etnodes[i][1+2]][0] + xyrr[etnodes[i][2+2]][0] + xyrr[etnodes[i][3+2]][0]) / 4;
      auxxyzv[i+nncm][1] = (xyrr[etnodes[i][0+2]][1] + xyrr[etnodes[i][1+2]][1] + xyrr[etnodes[i][2+2]][1] + xyrr[etnodes[i][3+2]][1]) / 4;
      auxxyzv[i+nncm][2] = 0.0;
      auxxyzv[i+nncm][3] = theta[i];
    }
    else
      auxxyzv[i+nncm] = NULL;
  }

  nmidp = nncm;
  for (i=0;i<necm;i++){
    te = etnodes[i][0];
    nmidp++;

    switch (te){
    case planeelementlt:{
      auxnod[i] = new long [4];
      auxnod[i][0] = 1;   auxnod[i][ 1] = etnodes[i][0+2]+1;   auxnod[i][ 2] = etnodes[i][1+2]+1;   auxnod[i][ 3] = etnodes[i][2+2]+1;

      break;
    }
    case planeelementqt:{
      auxnod[i] = new long [13];
      auxnod[i][0] = 4;   auxnod[i][ 1] = etnodes[i][0+2]+1;   auxnod[i][ 2] = etnodes[i][3+2]+1;   auxnod[i][ 3] = etnodes[i][5+2]+1;
                          auxnod[i][ 4] = etnodes[i][3+2]+1;   auxnod[i][ 5] = etnodes[i][1+2]+1;   auxnod[i][ 6] = etnodes[i][4+2]+1;
                          auxnod[i][ 7] = etnodes[i][4+2]+1;   auxnod[i][ 8] = etnodes[i][2+2]+1;   auxnod[i][ 9] = etnodes[i][5+2]+1;
                          auxnod[i][10] = etnodes[i][3+2]+1;   auxnod[i][11] = etnodes[i][4+2]+1;   auxnod[i][12] = etnodes[i][5+2]+1;

      break;
    }
    case planeelementlq:{
      auxnod[i] = new long [13];
      auxnod[i][0] = 4;   auxnod[i][ 1] = etnodes[i][0+2]+1;   auxnod[i][ 2] = etnodes[i][1+2]+1;   auxnod[i][ 3] = nmidp;
                          auxnod[i][ 4] = etnodes[i][1+2]+1;   auxnod[i][ 5] = etnodes[i][2+2]+1;   auxnod[i][ 6] = nmidp;
                          auxnod[i][ 7] = etnodes[i][2+2]+1;   auxnod[i][ 8] = etnodes[i][3+2]+1;   auxnod[i][ 9] = nmidp;
                          auxnod[i][10] = etnodes[i][3+2]+1;   auxnod[i][11] = etnodes[i][0+2]+1;   auxnod[i][12] = nmidp;

      break;
    }
    case planeelementqq:{
      auxnod[i] = new long [25];
      auxnod[i][0] = 8;   auxnod[i][ 1] = etnodes[i][0+2]+1;   auxnod[i][ 2] = etnodes[i][4+2]+1;   auxnod[i][ 3] = nmidp;
                          auxnod[i][ 4] = etnodes[i][4+2]+1;   auxnod[i][ 5] = etnodes[i][1+2]+1;   auxnod[i][ 6] = nmidp;
                          auxnod[i][ 7] = etnodes[i][1+2]+1;   auxnod[i][ 8] = etnodes[i][5+2]+1;   auxnod[i][ 9] = nmidp;
                          auxnod[i][10] = etnodes[i][5+2]+1;   auxnod[i][11] = etnodes[i][2+2]+1;   auxnod[i][12] = nmidp;
                          auxnod[i][13] = etnodes[i][2+2]+1;   auxnod[i][14] = etnodes[i][6+2]+1;   auxnod[i][15] = nmidp;
                          auxnod[i][16] = etnodes[i][6+2]+1;   auxnod[i][17] = etnodes[i][3+2]+1;   auxnod[i][18] = nmidp;
                          auxnod[i][19] = etnodes[i][3+2]+1;   auxnod[i][20] = etnodes[i][7+2]+1;   auxnod[i][21] = nmidp;
                          auxnod[i][22] = etnodes[i][7+2]+1;   auxnod[i][23] = etnodes[i][0+2]+1;   auxnod[i][24] = nmidp;

      break;
    }
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function");
      fprintf (stderr," adaptivity::print_thetacon (%s, line %d)",__FILE__,__LINE__);
    }
    }
  }

  print_confile (file,nncm,necm,auxxyzv,auxnod,dim);

  for (i=0;i<nncm+necm;i++)
    delete [] auxxyzv[i];
  for (i=0;i<necm;i++)
    delete [] auxnod[i];
  delete [] auxxyzv;  delete [] auxnod;
}


void adaptivity::print_coarsemesh ( double *ei2)
{
  long i,j,te,nne,*cn;
  double r[2];
  ivector nodes;
  char file[1023];
  FILE *cmchar;
  
  sprintf (file,"%scmeshchr%s",Mp->path,ni);
  cmchar = fopen (file,"w");                 if (cmchar==NULL) fprintf (stderr,"\n .cmchar file has not been specified.\n Try it again!\n");
  
  // ** nastaveno pro 2D , lcid = 0  , codenumbers na noudech , d na prvku konstantni!!!!!!!!!!!!!!!!!!
  
  fprintf (cmchar,"%ld\n\n",nn);
  
  for (i=0;i<nn;i++){
    cn = Gtm->gnodes[i].cn;
    for (j=0;j<2;j++){
      if        (cn[j]>0)   r[j]=Lsrs->lhs[cn[j]-1];
      else if   (cn[j]==0)  r[j]=0.0;
           else             r[j]=Mb->lc[0].pd[0-cn[j]-1];
    }
    fprintf (cmchar," %15.7e %15.7e %15.7e %15.7e",Gtm->gnodes[i].x,Gtm->gnodes[i].y,r[0],r[1]);
    fprintf (cmchar," %15.7e %15.7e %15.7e %6ld"  ,rderfull[i],rderfull[nn+i],rderfull[2*nn+i],Gtm->nadjelnod[i]);
    for (j=0;j<Gtm->nadjelnod[i];j++)
      fprintf (cmchar,"%6ld",Gtm->adjelnod[i][j]);
    fprintf (cmchar,"\n");
  }
  
  fprintf (cmchar,"\n\n%ld\n\n",ne);
  
  for (i=0;i<ne;i++){
    te = Mt->give_elem_type (i);
    nne = Mt->give_nne_inner ((elemtype)te);
    allocv (nne,nodes);
    
    Gtm->give_nodes (i,nodes);
    
    fprintf (cmchar,"  %2ld  %2ld",te,nne);
    for (j=0;j<nne;j++)
      fprintf (cmchar,"%6ld ",nodes.a[j]);
    
    fprintf (cmchar,"%15.7f\n",ei2[i]);
    
    destrv (nodes);
  }
  
  fclose (cmchar);
}