adaptivity.cpp

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

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

#include "adaptivity.h"
#include "global.h"
#include "elemhead.h"
#include "alias.h"
#include "vector.h"
#include "matrix.h"
#include "gtopology.h"
#include "gadaptivity.h"
#include "least_square.h"
#include "z2_smoothing.h"
#include "element.h"
#include "intpoints.h"
#include "loadcase.h"
#include "mechprint.h"
#include "j2flow.h"   // docasne

#include "mechmat.h"



#include "vecttens.h"


/*
predpoklad - na elementu je konstantni matice tuhosti d
funguje pro domenu s jednim materialem

printflags:
1 - print stdout
2 - print bgm
4 - print contoures
8 - print cmeshch.*
16 - print test
32 - print dx
64 - print rderfull
*/


adaptivity::adaptivity ()
{
  dim = Mt->give_dimension (0);
  ord = Mt->give_degree (0);
  ncomp = Mt->give_ncomp (0);
  nn = Mt->nn;
  ne = Mt->ne;
  
  tad = Mp->adaptivity;
  printflags = Mp->adaptflag;
  corr = Mp->corr;
  limf = Mp->limf;
  clos = Mp->clos;
  plra = Mp->plra;
  
  otherflags = -1;
  answer = 0;
  
  ni = new char[15];
  body = new long[16];
}

adaptivity::~adaptivity ()
{
  delete ni;
  delete [] body;
}


/**
   This is main function of class adaptivity, it runs in three steps:
   1. preparing strain
   2. computing of refined derivative (strain or stress, determined by otherflag) by spr_smoothing or z2_smoothing
   3. computing of error of FEM solution
   (case 3 = computing of accurate error)
   
   created  3.3.2002, Ladislav Svoboda, [e-mail]
*/
long adaptivity::run (long of,int width,long nincr)
{
  long hod,min;
  double sec;
  
  if (Mespr==1){
    //    fprintf (stdout,"\n\n *********************************");
    fprintf (stdout,  "\n\n ***   Adaptivity is started   ***");
    
    sec = clock();
  }
  
  otherflags = of;
  answer = 0;
  
  //if (otherflags & 8)  rderfull = new double [ncomp*nn*2];
  //else                 rderfull = new double [ncomp*nn];
  rderfull = new double [ncomp*nn];
  
  prepare_ni (width,nincr);
  prepare_strain ();
  
  //if (!(otherflags & 8)){
  switch (tad){
  case 1:{
    spr ();
    compute_error ();
    break;
  }
  case 2:{
    z2_smoothing z2_est(ncomp);
    z2_est.run (rderfull);
    compute_error ();
    break;
  }
  case 3:{
    spr ();
    compute_error_fine ();
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown type of error smoothing is required in function");
    fprintf (stderr," adaptivity::run (%s, line %d)",__FILE__,__LINE__);
  }
  }
  //}
  //else{
  //  spr_nonlinear ();
  //  compute_error ();
  //}
  
  
  delete [] rderfull;
  
  if (Mespr==1){
    sec = (clock() - sec) / (double)CLOCKS_PER_SEC;
    hod = (long)sec/3600;  sec -= hod*3600;
    min = (long)sec/60;    sec -= min*60;
    fprintf(stdout,"\n Consumed time by Ada   %2ld:%02ld:%05.2f",hod,min,sec);
    
    fprintf (stdout,"\n ***   Adaptivity is finished  ***");
    //    fprintf (stdout,"\n *********************************");
  }
  return (answer);
}


/**
   
   
 */
void adaptivity::prepare_ni (int width,long nincr)
{
  
  switch (Mp->tprob){
  case linear_statics:   {
    char *fix = strchr(Mp->suffix+1,'.');
    if (fix==NULL)
      ni[0] = '\0';
    else
      strcpy (ni,fix);
    
    break;
  }
  case mat_nonlinear_statics:{
    sprintf (ni,".%0*ld",width,nincr);
    
    break;
  }
  default:{
    fprintf (stderr,"\n\n unknown problem type is required in function prepare_ni (file %s, line %d).\n",__FILE__,__LINE__);
  }
  }
  
}


/**
   Function calls 'mipstrains' for computing strains in main int. points
   and sets up 'body' - indicator of integration points.
   body[i]=0 - unusual int. points , body[i]=1 - usual int. points
   
   created  3.4.2002, Ladislav Svoboda, [e-mail]
*/
void adaptivity::prepare_strain (void)
{
  if ( (!(otherflags & 8) && !(otherflags & 4)) || ((otherflags & 8) && !(otherflags & 2)) )
    //Mm->mipstrains (lcid);
  
  
  if (Pelt!=NULL)
    if ( Pelt->nb==1 && Pelt->nip[0][0]==1 && Pelt->ncomp[0]==3 )
      body[0] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Pelt has not usual intpoints\n");  body[0] = 0; }
  
  if (Peqt!=NULL)
    if ( Peqt->nb==2 && Peqt->nip[0][0]==3 && Peqt->nip[1][1]==3 && Peqt->nip[1][0]==0 && Peqt->nip[0][1]==0 && Peqt->ncomp[0]==2 && Peqt->ncomp[1]==1)
      body[2] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Peqt has not usual intpoints\n");  body[2] = 0; }
  
  if (Pelq!=NULL)
    if ( Pelq->nb==2 && Pelq->nip[0][0]==4 && Pelq->nip[1][1]==1 && Pelq->nip[1][0]==0 && Pelq->nip[0][1]==0 && Pelq->ncomp[0]==2 && Pelq->ncomp[1]==1)
      body[5] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Pelq has not usual intpoints\n");  body[5] = 0; }
  
  if (Peqq!=NULL)
    if ( Peqq->nb==2 && Peqq->nip[0][0]==4 && Peqq->nip[1][1]==4 && Peqq->nip[1][0]==0 && Peqq->nip[0][1]==0 && Peqq->ncomp[0]==2 && Peqq->ncomp[1]==1)
      body[7] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Peqq has not usual intpoints\n");  body[7] = 0; }
  
  if (Cct!=NULL)
    if ( Cct->nb==2 && Cct->nip[0][0]==1 && Cct->nip[1][1]==1 && Cct->nip[1][0]==0 && Cct->nip[0][1]==0 && Cct->ncomp[0]==3 && Cct->ncomp[1]==2)
      body[10] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Cct has not usual intpoints\n");   body[10] = 0; }
  
  if (Ltet!=NULL)
    if ( Ltet->nb==1 && Ltet->nip[0][0]==1 && Ltet->ncomp[0]==6 )
      body[15] = 1;
    else{
      fprintf (stderr,"\n\n Warning: Ltet has not usual intpoints\n");  body[15] = 0; }
}


/**
   Function fills up array rderfull by derivatives(strain or stresses - it depends on otherflags) in nodes.
   It fills up array 'spder' by strain or stresses in sampling points
   and calls spr_smoothing to compute refined derivatives in nodes.
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
*/
void adaptivity::spr (void)
{
  long i;
  double **spder;
  least_square spr(dim,ncomp,Gtm,0);
  
  spder = new double* [ne];
  
  for (i=0;i<ne;i++){
    switch (Mt->elements[i].te){
    case planeelementlt:   { Pelt->elchar  (i,spder[i],otherflags); break;}
    case planeelementqt:   { Peqt->elchar  (i,spder[i]); break;}
    case planeelementlq:   { Pelq->elchar  (i,spder[i]); break;}
    case planeelementqq:   { Peqq->elchar  (i,spder[i]); break;}
    case lineartet:        { Ltet->elchar  (i,spder[i]); break;}
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function spr (%s, line %d)",__FILE__,__LINE__);
    }
    }
  }
  
  spr.spr_default (Out,spder,rderfull);
  
  for (i=0;i<ne;i++)  delete [] spder[i];
  delete [] spder;
}


/**
   Function fills up array rderfull by strain(in first half) and stress(in second half) in nodes.
   1. it fills up array 'spder' by strain in sampling points and calls spr_smoothing to compute refined strain in nodes.
   1. it fills up array 'spder' by stress in sampling points and calls spr_smoothing to compute refined stress in nodes.
   
   created  3.5.2002, Ladislav Svoboda, [e-mail]
*/
void adaptivity::spr_nonlinear (void)
{
  long i,of;
  double **spder;
  least_square spr(dim,ncomp,Gtm,0);
  
  spder = new double* [ne];
  
  of = ((otherflags & 4) ? 6 : 2);
  
  for (i=0;i<ne;i++){
    switch (Mt->elements[i].te){
      
    case planeelementlt:   { Pelt->elchar  (i,spder[i],of); break;}
      
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function spr_nonlinear (%s, line %d)",__FILE__,__LINE__);
    }
    }
  }
  
  spr.spr_default (Out,spder,rderfull + ncomp*nn);
  
  for (i=0;i<ne;i++)  delete [] spder[i];
  
  
  for (i=0;i<ne;i++){
    switch (Mt->elements[i].te){
      
    case planeelementlt:   { Pelt->elchar  (i,spder[i],0); break;}
      
    default:{
      fprintf (stderr,"\n\n unknown element type is required in function spr_nonlinear (%s, line %d)",__FILE__,__LINE__);
    }
    }
  }
  
  spr.spr_default (Out,spder,rderfull);
  
  for (i=0;i<ne;i++)  delete [] spder[i];
  delete [] spder;
}


/**
 */
void adaptivity::compute_error (void)
{
  if (otherflags & 1)
    nonlin_error ();
  else
    lin_error ();
}


/**
   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::lin_error (void)
{
  long i;
  double e2     ,u2     ,volume   ,corr1,corr2;
  vector ei2(ne),ui2(ne),sizel(ne),refsizel(ne);
  
  e2 = u2 = 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];
  }
  u2 += e2;
  
  compute_refsizel_lin (sizel.a,ei2.a,e2,u2,ord,ne,refsizel.a);
  
  // docasne
  answer = 1;
  print (volume,e2,u2,ei2.a,ui2.a,refsizel.a);
}


/**
   Function computes new=required characteristic size of elements.
   
   @param sizel - array of old characteristic sizes of element
   @param
   
   created  3.6.2002, Ladislav Svoboda, [e-mail]
 */
void adaptivity::compute_refsizel_lin ( double *sizel, double *ei2,double e2,double u2,long ord,long ne,double *refsizel)
{
  long i;
  double em2,specacc;
  vector epsil;
  
  allocv (ne,epsil);
  
  specacc = Mp->adapt_accuracy;
  em2 = specacc*specacc * u2/ne;
  
  //corr = 0.8;  //<0.7;1.0>
  
  for (i=0;i<ne;i++){
    epsil[i] = sqrt(ei2[i]/em2);
    if (epsil[i]<0.05) epsil[i] = 0.05;
    refsizel[i] = sizel[i] / pow(epsil[i] , 1.0/(double)ord * corr);
  }
  
  // ***********************************************
  //char file[255];
  //sprintf (file,"%sepsil.dx%s",Mp->path,ni);
  //print_valel (Gtm,Mp,file,"epsil",epsil.a,'d');
  // ***********************************************
}


/**