comparison.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           comparison.cpp
// author(s):      Ladislav Svoboda
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

// toto je kvuli WIN32_LEAN_AND_MEAN
#include "types.h"

#include "comparison.h"
#include "problem.h"
#include "mesh.h"


namespace mumech {


double Comparison :: cmp_mm_ansys_total_energy (const char *basename, const char *vtkGeomFile, int lc)
{
  fprintf (stdout, "cmp: \n");

  char name[1023];

  // MESH
  Mesh *mesh = new Mesh(0, NULL);
  sprintf (name, "%sANSYS/rslts_lc%d/list.IB.mesh.vtu", basename, lc);   // orig: IB
  mesh->read_geometry_file_vtk (name, -1, -1, false);
  
  // muMECH reseni
  Problem *pMM = new Problem;
  pMM->read_inclusions_plus_initialize_and_print (vtkGeomFile, NULL, DT_NUMERICAL);
  int pidMM = mesh->add_problem(pMM);
  mesh->compute_element_fields ('t', false, true,  true,  pidMM, lc, 1, PFCM_FULL);  // termitovo ERROR kdyz tam dam optimize tak se to sere
  mesh->compute_node_fields    ('t', true,  false, false, pidMM, lc, 1, PFCM_FULL);  // termitovo ERROR kdyz tam dam optimize tak se to sere
  
  // ANSYS reseni
  ProblemFEM *pAS = new ProblemFEM (pMM->give_dimension(), pMM->give_nLC());
  int pidAS = mesh->add_problem(pAS);
  sprintf (name, "%sANSYS/rslts_lc%d/list.IB.rslt.vtu", basename, lc);
  mesh->read_geometry_file_vtk (name, pidAS, lc, true);
  
  
  // ENERGY
  // energy error
  double Err_all = mesh->give_error(pidMM, pidAS, lc);
  fprintf (stdout, "Error: %.3lf%%\n", Err_all);
  

  // PRINT
  if (1) {
    sprintf (name, "%s%s", basename, "cmp_energy.meshMM.vtk");
    mesh->print_geometry_file_vtk (name, pidMM, lc, 1);
    
    sprintf (name, "%s%s", basename, "cmp_energy.meshAS.vtk");
    mesh->print_geometry_file_vtk (name, pidAS, lc, 1);
  }
  
  delete mesh;
  delete pMM;
  delete pAS;

  return Err_all;
}


void Comparison :: give_energy_MM ()
{
  fprintf (stdout, "cmp: \n");
  
  char name[1023];
  
  Mesh *mesh = NULL;
  
  int pidMM;
  double Eh, Et, Ep;

  double p1[3];  p1[2] = 0.0;
  double p2[3];  p2[2] = 0.0;
  long   nn[3];  nn[2] = 0;
  double x1, x2, y1, y2, dx, dy;
  x1 = y1 = -5.0;
  x2 = y2 =  5.0;
  dx = dy = 0.02;
  double delta = 1.0;
  
  // muMECH reseni
  Problem *pMM = new Problem;
  pMM->generate_equiv_1x2dI_2d(1.0, 0.2, NULL, 0.0, 0.0, 10, 0.3, 1.0, 0.5, 0.0, DT_NUMERICAL);
  
  int steps = 13;

  FILE *file = fopen ("energy_MM.txt", "wt");

  sprintf (name, "energy_MM.lc0.%2ld.%4.2lf.dat", steps, dx);  FILE *file0 = fopen (name, "wt");
  sprintf (name, "energy_MM.lc1.%2ld.%4.2lf.dat", steps, dx);  FILE *file1 = fopen (name, "wt");
  sprintf (name, "energy_MM.lc2.%2ld.%4.2lf.dat", steps, dx);  FILE *file2 = fopen (name, "wt");
  
  
  for (int step=0; step<steps; step++) {
    // regular mesh
    p1[0] = x1 - step*delta;       p1[1] = y1 - step*delta;
    p2[0] = x2 + step*delta;       p2[1] = y2 + step*delta;
    nn[0] = (p2[0] - p1[0]) / dx;  nn[1] = (p2[1] - p1[1]) / dy;
    
    fprintf (file, "BB %lf \n", p2[0]);
    
    // MESH
    mesh = new Mesh(0, NULL);
    mesh->generate_regular_mesh(p1, p2, nn);
    pidMM = mesh->add_problem(pMM);
    mesh->compute_element_fields ('t', false, true,  true,  pidMM, 0, 0, PFCM_FULL);  // termitovo ERROR kdyz tam dam optimize tak se to sere
    //mesh->compute_node_fields    ('t', true,  false, false, pidMM, lc, 1, PFCM_FULL);  // termitovo ERROR kdyz tam dam optimize tak se to sere
    
    
    for (int lc=0; lc<3; lc++) {
      fprintf (file, " lc = %ld\n", lc);

      Eh = mesh->give_homog_energy(pidMM, lc);
      fprintf (file, "  lc%ld Ehom = %lf\n", lc, Eh);
      
      Et = mesh->give_total_energy(pidMM, lc);
      fprintf (file, "  lc%ld Etot = %lf\n", lc, Et);
      
      fprintf (file, "  lc%ld Epert num = Etotal - Ehomog = %lf\n", lc, Et - Eh);
      
      if (lc == 0)  fprintf (file0, "%lf %lf\n", p2[0], Et - Eh);
      if (lc == 1)  fprintf (file1, "%lf %lf\n", p2[0], Et - Eh);
      if (lc == 2)  fprintf (file2, "%lf %lf\n", p2[0], Et - Eh);
      
      //Ep = pMM->compute_supplement_energy (lc);
      //fprintf (file, "Epert anal = %lf\n", Ep);
    }
    
    delete mesh;
  }
  
  
  fclose(file);
  fclose(file0);
  fclose(file1);
  fclose(file2);
  
  //if (1) {
  //  sprintf (name, "%s%s", basename, "cmp_energy.meshMM.vtk");
  //  mesh->print_geometry_file_vtk (name, pidMM, lc, 1);
  //  
  //  sprintf (name, "%s%s", basename, "cmp_energy.meshAS.vtk");
  //  mesh->print_geometry_file_vtk (name, pidAS, lc, 1);
  //}
  
  //delete mesh;
  delete pMM;
}



} // end of namespace mumech

/*end of file*/