src/Elasticity/main.cc

Go to the documentation of this file.

// ------------------------------------------------------------
// main_enw.cc
//
// ------------------------------------------------------------

#include "entw.h"
#include "printlas.h"
#include "eingabe.h"      


// new
#define ddx(U) DX_FE(U)*Grid_stretch
#define ddy(U) DY_FE(U)*Grid_stretch
#define ddz(U) DZ_FE(U)
               
#define u_exakt zz


double Grid_stretch, finest_mesh_size;
double mesh_z;
double stretch_x,stretch_y,stretch_z;
double eps;

exit_type Solver(// parameter for grid generation
                 Grid_gen_parameters ggen_parameter,
                 int test_processors,
                 double meshsize,
                 double grid_stretch, // stretch grid in z-direction
                 All_Domains* domain,
                 // Parameter for Poisson
                 double kappa_x, double kappa_y, double kappa_z,
                 bool (*Neumann_boundary)(double x, double y, double z),
                 bool (*Third_boundary)(double x, double y, double z),
                 double beta_third,
                 double (*T_bulk)(double x,double y,double z),
                 double (*temp_Dirichlet)(double x,double y,double z),
                 double (*heat_source)(double x,double y,double z), 
                           // heat beschrieben durch Funktion
                           // dies wird verwendet falls input_object==NULL
                 Input_data_object* input_object,   // heat beschrieben
                                                    // falls nicht verwendet,
                                                    // muss Zeiger NULL sein
                 bool normalize_heat_integral,  
                           // falls normalize_heat_integral==true dann
                           // ist das Integral ueber heat folgender Wert:
                 double heat_integral, 
                 // Parameter for elasticity
                 double alpha,
                 double mu, double E,
                 double T_0,
                 // Parameter for linear-nonlinear
                 bool nonlinear_kappa,
                 bool nonlinear_alpha,
                 double (*Material_A)(double x,double y,double z),
                 double (*nonlinear_kappa_A)(double x),
                 double (*nonlinear_kappa_B)(double x),
                 double (*nonlinear_alpha_A)(double x),
                 double (*nonlinear_alpha_B)(double x),
                 // Parameter solver
                 int iteration_Poi,
                 int iteration_elas,
                 int solver_typ_elas,
                 double eps,
                 double bound_poiss,
                 double bound_elas, 
                 // Print Infos?
                 bool infos,
                 // pointer on ThreadComm-Class
                 ThreadComm* comm, 
                 // for Parallelization
                 MPI_Comm mpi_comm);

// Parameter for linear-nonlinear
double Material_A(double x,double y,double z) {
  // Rueckgabe nur 1.0 und 0.0
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;

  if(z > (0.2))
    return 1.0;
  return 0.0;
}

// Kappa_x = nonlinear_kappa * kappa_x
// Kappa_y = nonlinear_kappa * kappa_y
// Kappa_z = nonlinear_kappa * kappa_z
double nonlinear_kappa_A(double T) {
  return 33.132 / (T - 42.871) / 0.0103;
}

double nonlinear_kappa_B(double T) {
  return 33.132 / (T - 42.871) / 0.0103;
}
  
double nonlinear_alpha_A(double T) {
  return 0.0000069;
}

double nonlinear_alpha_B(double T) {
  return 0.0000069;
}

/*
// Markierung des Randes
bool Neumann_boundary(double x, double y, double z) {
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;
  double r;
  r = sqrt(x*x+y*y);

  if(z > 0.8-eps) return true;
  else return false;
}

bool Third_boundary(double x, double y, double z) {
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;
  double r;
  r = sqrt(x*x+y*y);
  return false;
}

// Funktionen fuer den Rand
double T_bulk(double x, double y, double z) {
  double r;
  r = sqrt(z*z+y*y);

  return 0.0;
}

double temp_Dirichlet(double x, double y, double z) {
  double r;
  r = sqrt(x*x+y*y);

  return 100.0;
}
*/

// Markierung des Randes
bool Neumann_boundary(double x, double y, double z) {
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;
  double r;
  r = sqrt(x*x+y*y);

// Test ob dann nicht NaN

  return false;

  if(x < 0.0) 
    return true;
  else 
    return false;


  /*
  if(z < 3.0-eps) return true;
  // if(z < 2.5) return true;
  else return false;
  */
}

bool Third_boundary(double x, double y, double z) {
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;
  double r;
  r = sqrt(x*x+y*y);
  return false;
}

// Funktionen fuer den Rand
double T_bulk(double x, double y, double z) {
  double r;
  r = sqrt(z*z+y*y);

  return 0.0;
}

double temp_Dirichlet(double x, double y, double z) {
  double r;
  r = sqrt(x*x+y*y);

  return 0.0;
}


// Funktion fuer Waerme
double Heat_source(double x, double y, double z) {
  double radius2;

  return 1.0;
 
  /*
  // fuer Streckung des Gitters! wichtig!
  z = z * Grid_stretch;

  z = z - 1.5;

  //  return 0.5-x;

  radius2 = z*z+y*y;

  if(radius2 <= 0.2)
    return 1.0;
  else
    return 0.0;
  */
}

int main(int argc, char** argv) {
  cout.precision(12);
  cout.setf(ios::fixed,ios::floatfield);
  ifstream PARAMETER;
  int test_processors;

  int my_rank;        /* Rank of process */

  //  Index3D(1,1,1).son(ENTd).coordinate().Print();


  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

  PARAMETER.open("para.dat",ios :: in);

  // stretch domain
  stretch_x =6.0;
  stretch_y =3.0;
  stretch_z =3.0;

  // grid mesh size
  PARAMETER >> finest_mesh_size;

  // test processors or not
  PARAMETER >> test_processors;

  // stretch grid in z-direction
  PARAMETER >> Grid_stretch;

  /*
  mesh_z = finest_mesh_size * Grid_stretch;
  if(my_rank==0) {
    fout << "\n Finest meshsize: " << finest_mesh_size << endl;    
    fout << "Finest meshsize along z axis: " <<mesh_z<< endl;    
  }
  */

  // Parameter for elasticity
  double alpha;
  double mu, E;
  double T_0;

  mu    = 0.25;
  E     = 30000.0;
  alpha = 0.0000069;
  T_0   = 0.0;

  // Parameter for Poisson
  double heat_integral, kappa_x, kappa_y, kappa_z;
  double beta_third;
 
  PARAMETER >> heat_integral;

  kappa_x = 0.0103;
  kappa_y = 0.0103;
  kappa_z = 0.0103;
  beta_third = 0.2;

  // Print Infos?
  bool infos;
  infos =true;
  
  // parameters for solver
  int iteration_Poi;
  int iteration_elas;
  int solver_typ_elas;
  double bound_poiss, bound_elas;
  double eps;

  PARAMETER >> solver_typ_elas;
  PARAMETER >> iteration_Poi;
  //  iteration_Poi  = 7;
  PARAMETER >> iteration_elas;
  //  iteration_elas = 1;
  bound_poiss = 0.000000001;
  bound_elas  = 0.000000001;

  eps = 0.0001;

  // Parameter for linear-nonlinear
  bool nonlinear_kappa;
  bool nonlinear_alpha;

  PARAMETER >> nonlinear_kappa;
  nonlinear_alpha = true;
  //  nonlinear_alpha = false;

  //Exits
  exit_type exit_hs;
  exit_hs = no_ex;


  // describtion of domain
  D3vector Vm(-0.5, -0.5, 0.0);
  // D3vector Vm(-0.0, -0.0, 0.0);
  D3vector Vstretch(stretch_x,stretch_y,stretch_z/Grid_stretch);

  Cylinder cylinder;
  /*
  //  Square cube;
  //  Skew_square cube(-0.5, 0.5);

  double alpha_d, alpha_t;
  alpha_d = -45;
  alpha_t = 45;

  double alpha_d_norm=atan(tan(alpha_d/180.*M_PI)/stretch_z*stretch_x);
  double alpha_t_norm=atan(tan(alpha_t/180.*M_PI)/stretch_z*stretch_x);
  Skew_square cube(alpha_d_norm, alpha_t_norm);
  */

  double alpha_d, alpha_t, alpha_l, alpha_r;
  /*
  alpha_d = -30;
  alpha_t = 30;
  
  alpha_l = -5;
  alpha_r = 7;
  */
  alpha_d = -60;
  alpha_t = 60;
  
  alpha_l = -5.0;
  alpha_r = 7.0;
  
  double alpha_d_norm=atan(tan(alpha_d/180.*M_PI)/stretch_x*stretch_z);
  double alpha_t_norm=atan(tan(alpha_t/180.*M_PI)/stretch_x*stretch_z);

  double alpha_l_norm=atan(tan(alpha_l/180.*M_PI)/stretch_y*stretch_z);
  double alpha_r_norm=atan(tan(alpha_r/180.*M_PI)/stretch_y*stretch_z);
  

  Skew_squareDTLR cube(alpha_d_norm, alpha_t_norm, alpha_l_norm, alpha_r_norm);


  Six_corner_pyramid Sechseck(D3vector(0.0,1.0,1.0), 
                              D3vector(1.0,0.0,1.0),
                              D3vector(2.0,1.0,1.0),
                              D3vector(2.0,2.0,1.0), 
                              D3vector(1.0,3.0,1.0),
                              D3vector(0.0,2.0,1.0),
                              
                              D3vector(-0.3,1.0,0.0), 
                              D3vector(1.0,-0.3,0.0),
                              D3vector(2.3,1.0,0.0),
                              D3vector(2.3,2.0,0.0), 
                              D3vector(1.0,3.3,0.0),
                              D3vector(-0.3,2.0,0.0));


  // Gebiet 1
  //  Domain domain((cube + Vm) * Vstretch);
  // Gebiet 2
  //  Domain domain(Sechseck * Vstretch);
  // Gebiet 3
  /*
  Square square;
  D3vector ecke(-5.0,-5.0,-5.0);
  double   Laenge = 10.0;
  Domain domain((square*Laenge+ecke));
  */

  Ball   ball;
  Square square;

  D3vector V1(-0.5,0.5,0.0);
  D3vector V2(0.5,0.0,0.5);
  D3vector V3(1.0,0.5,0.5);

  //  definition of the domain:
  Domain domain(ball || ((square*V3) + V1)  || ((square*V3) + V2));


  // Trace pro Sachen:
  Input_data_object* P_input_object;

  /*
  PARAMETER.open("file.txt",ios :: in);
  Reader_trace_pro reader(PARAMETER);
  PARAMETER.close();

  Input_data_object input_object(reader.ort,
                                 reader.value,
                                 reader.anzahl_punkte,
                                 reader.mesh_size);

  P_input_object = NULL;
  P_input_object = &input_object;
  */

  P_input_object = NULL;

  ThreadComm commObj;
  ThreadComm* comm;
  comm = &commObj;

  Grid_gen_parameters ggen_parameter;
  // ggen_parameter.Set_grid_point(D3vector(0.0,1.5,0.2/Grid_stretch)); 

  double off_par;
  double stre_par;
  int n_genau;

  PARAMETER >> off_par;
  PARAMETER >> stre_par; 
  PARAMETER >> n_genau;
  PARAMETER.close();


  ggen_parameter.Set_r_parallel(n_genau);
  ggen_parameter.Set_offset_square(off_par);
  ggen_parameter.Set_stretch_square(stre_par);

  exit_hs = Solver(// parameter for grid generation
                   ggen_parameter,
                   test_processors,
                   finest_mesh_size,
                   Grid_stretch, // stretch grid in z-direction
                   &domain,
                   // Parameter for Poisson
                   kappa_x, kappa_y, kappa_z,
                   Neumann_boundary,
                   Third_boundary,
                   beta_third,
                   T_bulk,
                   temp_Dirichlet,
                   Heat_source, // heat beschrieben 
                                // durch Funktion
                                // dies wird verwendet 
                                // falls input_object==NULL
                   P_input_object,        // heat beschrieben
                                          // falls nicht verwendet,
                                          // muss Zeiger NULL sein
                   true, // ( normalize_heat_integral )
                           // falls  normalize_heat_integral==true dann
                           // ist das Integral ueber heat folgender Wert:
                   heat_integral, 
                   // Parameter for elasticity
                   alpha,
                   mu, E,
                   T_0,
                   // Parameter for linear-nonlinear
                   nonlinear_kappa,
                   nonlinear_alpha,
                   Material_A,
                   nonlinear_kappa_A,
                   nonlinear_kappa_B,
                   nonlinear_alpha_A,
                   nonlinear_alpha_B,
                   // parameter solver
                   iteration_Poi,
                   iteration_elas,
                   solver_typ_elas,
                   eps,
                   bound_poiss,
                   bound_elas, 
                   // Print Infos?
                   infos,
                   // pointer on ThreadComm-Class
                   comm,MPI_COMM_WORLD);


  if (exit_hs == hard_ex) { 
    fout << "Calculation killed. " << endl;
  }    


  if(my_rank==0) cout << " \n End: " << endl;

  MPI_Finalize();
}




// Output 
void Print_Datas(Variable& ux, Variable& uy, Variable& uz,
                 Variable& Temp,
                 Variable& rx, Variable& ry, Variable& rz,
                 Variable& One,
                 Variable& zz,           
                 double T_0,
                 double F,
                 double mu,
                 double alpha,
                 Cell_Variable& Mat_A,
                 ThreadComm* comm) {
  double spmax, spmin;
  int my_rank;

  my_rank = ux.Give_grid()->Give_my_rank();

  // Output for AVS
  ofstream DATEI;

  // Output of stress 
  zz  = (((ddx(ux) * (1.0-mu) + (ddy(uy) + ddz(uz)) * mu)
          / One) - ((1.0-mu) * alpha + mu * (alpha + alpha)) * (Temp - T_0))
    * F  | all_points;

  /*
  if(my_rank==0) DATEI.open("stress.dat",ios :: out);
  Print_UCD_moved(zz,rx,ry,rz,&DATEI);
  if(my_rank==0) DATEI.close();
  */

  spmax = Maximum(zz);
  spmin = Minimum(zz);
  if(my_rank==0) fout <<" Maximum sxx: "<< spmax <<" Minimum sxx: "
                      << spmin << endl;

  /*
  if(my_rank==0) DATEI.open("temp.dat",ios :: out);
  Print_UCD_moved(Temp,rx,ry,rz,&DATEI);
  if(my_rank==0) DATEI.close();
  */

  /*
  if(my_rank==0) DATEI.open("mat_A.dat",ios :: out);
  Print_UCD_moved(Mat_A,rx,ry,rz,&DATEI);
  if(my_rank==0) DATEI.close();
  */
}

---