homogenisation_disp_atom_test.cpp

Go to the documentation of this file.

/* This file is part of MoFEM.
 * MoFEM 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 3 of the License, or (at your
 * option) any later version.
 *
 * MoFEM 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 MoFEM. If not, see <http://www.gnu.org/licenses/>. */
 
#include <MoFEM.hpp>
using namespace MoFEM;
 
#include <Projection10NodeCoordsOnField.hpp>
#include <petsctime.h>
 
#include <FEMethod_LowLevelStudent.hpp>
#include <FEMethod_UpLevelStudent.hpp>
 
#include <adolc/adolc.h>
#include <NonLinearElasticElement.hpp>
#include <Hooke.hpp>
 
#include <boost/shared_ptr.hpp>
#include <boost/numeric/ublas/vector_proxy.hpp>
#include <boost/iostreams/tee.hpp>
#include <boost/iostreams/stream.hpp>
#include <fstream>
#include <iostream>
 
#include <MethodForForceScaling.hpp>
#include <TimeForceScale.hpp>
#include <BCs_RVELagrange_Disp.hpp>
#include <VolumeCalculation.hpp>
 
#include <boost/ptr_container/ptr_map.hpp>
#include <PostProcOnRefMesh.hpp>
#include <PostProcStresses.hpp>
 
 
 
 
static char help[] = "...\n\n";
#define RND_EPS 1e-6
 
 
//Rounding
double roundn(double n)
{
    //break n into fractional part (fract) and integral part (intp)
  double fract, intp;
  fract = modf(n,&intp);
 
  //    //round up
  //    if (fract>=.5)
  //    {
  //        n*=10;
  //        ceil(n);
  //        n/=10;
  //    }
  //    //round down
  //    if (fract<=.5)
  //    {
  //        n*=10;
  //        floor(n);
  //        n/=10;
  //    }
  // case where n approximates zero, set n to "positive" zero
  if (abs(intp)==0)
  {
    if(abs(fract)<=RND_EPS)
    {
      n=0.000;
    }
  }
  return n;
}
 
 
int main(int argc, char *argv[]) {
 
  MoFEM::Core::Initialize(&argc,&argv,(char *)0,help);
 
  try {
 
  moab::Core mb_instance;
  moab::Interface& moab = mb_instance;
  int rank;
  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 
  //Reade parameters from line command
  PetscBool flg = PETSC_TRUE;
  char mesh_file_name[255];
  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
  if(flg != PETSC_TRUE) {
    SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
  }
  PetscInt order;
  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
  if(flg != PETSC_TRUE) {
    order = 1;
  }
 
  //Read mesh to MOAB
  const char *option;
  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
  if(pcomm == NULL) pcomm =  new ParallelComm(&moab,PETSC_COMM_WORLD);
 
  //We need that for code profiling
  PetscLogDouble t1,t2;
  PetscLogDouble v1,v2;
  ierr = PetscTime(&v1); CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
 
  //Create MoFEM (Joseph) database
  MoFEM::Core core(moab);
  MoFEM::Interface& m_field = core;
 
  //ref meshset ref level 0
  ierr = m_field.seed_ref_level_3D(0,0); CHKERRQ(ierr);
 
  // stl::bitset see for more details
  BitRefLevel bit_level0;
  bit_level0.set(0);
  EntityHandle meshset_level0;
  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
  ierr = m_field.get_entities_by_ref_level(bit_level0,BitRefLevel().set(),meshset_level0); CHKERRQ(ierr);
 
  /***/
  //Define problem
 
  //Fields
  int field_rank=3;
  ierr = m_field.add_field("DISPLACEMENT",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,field_rank,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
  ierr = m_field.add_field("LAGRANGE_MUL_DISP",H1,AINSWORTH_LEGENDRE_BASE,field_rank); CHKERRQ(ierr);
 
  //Add entitities to the fields
  //============================================================================================================
  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"DISPLACEMENT"); CHKERRQ(ierr);
  ierr = m_field.add_ents_to_field_by_type(0,MBTET,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
  Range SurfacesFaces;
  ierr = m_field.get_cubit_msId_entities_by_dimension(103,SIDESET,2,SurfacesFaces,true); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SideSet 103 = %d\n",SurfacesFaces.size()); CHKERRQ(ierr);
 
  EntityHandle BoundFacesMeshset;
  rval = moab.create_meshset(MESHSET_SET,BoundFacesMeshset); CHKERRQ_MOAB(rval);
    rval = moab.add_entities(BoundFacesMeshset,SurfacesFaces); CHKERRQ_MOAB(rval);
  ierr = m_field.seed_ref_level_MESHSET(BoundFacesMeshset,BitRefLevel().set()); CHKERRQ(ierr);
  ierr = m_field.add_ents_to_field_by_type(BoundFacesMeshset,MBTRI,"LAGRANGE_MUL_DISP"); CHKERRQ(ierr);
  //============================================================================================================
 
  //set approximaiton order
  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
  ierr = m_field.set_field_order(0,MBTET,"DISPLACEMENT",order); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBTRI,"DISPLACEMENT",order); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBEDGE,"DISPLACEMENT",order); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBVERTEX,"DISPLACEMENT",1); CHKERRQ(ierr);
 
  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
 
  ierr = m_field.set_field_order(0,MBTRI,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBEDGE,"LAGRANGE_MUL_DISP",order); CHKERRQ(ierr);
  ierr = m_field.set_field_order(0,MBVERTEX,"LAGRANGE_MUL_DISP",1); CHKERRQ(ierr);
 
 
  //Define FE
  //define eleatic element
  boost::shared_ptr<Hooke<adouble> > hooke_adouble(new Hooke<adouble>);
  boost::shared_ptr<Hooke<double> > hooke_double(new Hooke<double>);
  NonlinearElasticElement elastic(m_field,2);
  ierr = elastic.setBlocks(hooke_double,hooke_adouble); CHKERRQ(ierr);
  ierr = elastic.addElement("ELASTIC","DISPLACEMENT"); CHKERRQ(ierr);
  ierr = elastic.setOperators("DISPLACEMENT","MESH_NODE_POSITIONS",false,true); CHKERRQ(ierr);
 
 
  BCs_RVELagrange_Disp lagrangian_element(m_field);
  lagrangian_element.addLagrangiangElement("LAGRANGE_ELEM","DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS");
 
    //define problems
  ierr = m_field.add_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
 
  //set finite elements for problem
  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","ELASTIC"); CHKERRQ(ierr);
  ierr = m_field.modify_problem_add_finite_element("ELASTIC_MECHANICS","LAGRANGE_ELEM"); CHKERRQ(ierr);
 
 
  //set refinement level for problem
  ierr = m_field.modify_problem_ref_level_add_bit("ELASTIC_MECHANICS",bit_level0); CHKERRQ(ierr);
 
  /***/
  //Declare problem
 
  /****/
  //build database
 
  //build field
  ierr = m_field.build_fields(); CHKERRQ(ierr);
  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
 
  //build finite elemnts
  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
 
  //build adjacencies
  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
 
  //build problem
  ierr = m_field.build_problems(); CHKERRQ(ierr);
 
 
  /****/
  //mesh partitioning
 
  //partition
  ierr = m_field.partition_problem("ELASTIC_MECHANICS"); CHKERRQ(ierr);
  ierr = m_field.partition_finite_elements("ELASTIC_MECHANICS"); CHKERRQ(ierr);
  //what are ghost nodes, see Petsc Manual
  ierr = m_field.partition_ghost_dofs("ELASTIC_MECHANICS"); CHKERRQ(ierr);
 
  //print bcs
  ierr = m_field.print_cubit_displacement_set(); CHKERRQ(ierr);
  ierr = m_field.print_cubit_force_set(); CHKERRQ(ierr);
  //print block sets with materials
  ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
 
  //create matrices (here F, D and Aij are matrices for the full problem)
  vector<Vec> F(6);
  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",ROW,&F[0]); CHKERRQ(ierr);
  for(int ii = 1;ii<6;ii++) {
    ierr = VecDuplicate(F[0],&F[ii]); CHKERRQ(ierr);
  }
 
  Vec D;
  ierr = m_field.VecCreateGhost("ELASTIC_MECHANICS",COL,&D); CHKERRQ(ierr);
 
  Mat Aij;
  ierr = m_field.MatCreateMPIAIJWithArrays("ELASTIC_MECHANICS",&Aij); CHKERRQ(ierr);
 
  ierr = VecZeroEntries(D); CHKERRQ(ierr);
  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = m_field.set_global_ghost_vector(
         "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
         ); CHKERRQ(ierr);
 
  for(int ii = 0;ii<6;ii++) {
    ierr = VecZeroEntries(F[ii]); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F[ii],INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  }
  ierr = MatZeroEntries(Aij); CHKERRQ(ierr);
 
  //Assemble F and Aij
  elastic.getLoopFeLhs().snes_B = Aij;
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC",elastic.getLoopFeLhs()); CHKERRQ(ierr);
 
  lagrangian_element.setRVEBCsOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",Aij,F);
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCLhs()); CHKERRQ(ierr);
  cout<<"after lagrangian_element.getLoopFeRVEBCLhs "<<endl;
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCRhs()); CHKERRQ(ierr);
  cout<<"after lagrangian_element.getLoopFeRVEBCRhs "<<endl;
//
//  //Matrix View
//  MatView(Aij,PETSC_VIEWER_DRAW_WORLD);//PETSC_VIEWER_STDOUT_WORLD);
//  std::string wait;
//  std::cin >> wait;
 
  ierr = MatAssemblyBegin(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
  ierr = MatAssemblyEnd(Aij,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
 
  for(int ii = 0;ii<6;ii++) {
    ierr = VecGhostUpdateBegin(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(F[ii],ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(F[ii]); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(F[ii]); CHKERRQ(ierr);
  }
 
//  ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
//  ierr = MatView(Aij,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
 
 
  // Volume calculation
  //==============================================================================================================================
  VolumeElementForcesAndSourcesCore vol_elem(m_field);
  Vec volume_vec;
  int volume_vec_ghost[] = { 0 };
  ierr = VecCreateGhost(
                        PETSC_COMM_WORLD,(!m_field.get_comm_rank())?1:0,1,1,volume_vec_ghost,&volume_vec
                        );  CHKERRQ(ierr);
  ierr = VecZeroEntries(volume_vec); CHKERRQ(ierr);
  vol_elem.getOpPtrVector().push_back(new VolumeCalculation("DISPLACEMENT",volume_vec));
 
  ierr = m_field.loop_finite_elements("ELASTIC_MECHANICS","ELASTIC", vol_elem);  CHKERRQ(ierr);
 
  ierr = VecAssemblyBegin(volume_vec); CHKERRQ(ierr);
  ierr = VecAssemblyEnd(volume_vec); CHKERRQ(ierr);
  double rve_volume;
  ierr = VecSum(volume_vec,&rve_volume);  CHKERRQ(ierr);
 
  ierr = VecView(volume_vec,PETSC_VIEWER_STDOUT_WORLD); CHKERRQ(ierr);
  ierr = PetscPrintf(PETSC_COMM_WORLD,"RVE Volume %3.2g\n",rve_volume); CHKERRQ(ierr);
  ierr = VecDestroy(&volume_vec);
 
  //==============================================================================================================================
  //Post processing
  //==============================================================================================================================
  struct MyPostProc: public PostProcVolumeOnRefinedMesh {
    bool doPreProcess;
    bool doPostProcess;
    MyPostProc(MoFEM::Interface &m_field):
    PostProcVolumeOnRefinedMesh(m_field),
    doPreProcess(true),
    doPostProcess(true)
    {}
    void setDoPreProcess() { doPreProcess = true; }
    void unSetDoPreProcess() { doPreProcess = false; }
    void setDoPostProcess() { doPostProcess = true; }
    void unSetDoPostProcess() { doPostProcess = false; }
 
    
 
    PetscErrorCode preProcess() {
      PetscFunctionBegin;
      if(doPreProcess) {
        ierr = PostProcVolumeOnRefinedMesh::preProcess(); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0);
    }
    PetscErrorCode postProcess() {
      PetscFunctionBegin;
      if(doPostProcess) {
        ierr = PostProcVolumeOnRefinedMesh::postProcess(); CHKERRQ(ierr);
      }
      PetscFunctionReturn(0);
    }
  };
 
  MyPostProc post_proc(m_field);
 
  ierr = post_proc.generateReferenceElementMesh(); CHKERRQ(ierr);
  ierr = post_proc.addFieldValuesPostProc("DISPLACEMENT"); CHKERRQ(ierr);
  ierr = post_proc.addFieldValuesGradientPostProc("DISPLACEMENT"); CHKERRQ(ierr);
  ierr = post_proc.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
  ierr = post_proc.addFieldValuesGradientPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
 
  for(
      map<int,NonlinearElasticElement::BlockData>::iterator sit = elastic.setOfBlocks.begin();
      sit != elastic.setOfBlocks.end(); sit++
      ) {
    post_proc.getOpPtrVector().push_back(
   new PostPorcStress(
    post_proc.postProcMesh,
    post_proc.mapGaussPts,
    "DISPLACEMENT",
    sit->second,
    post_proc.commonData,
    false
    )
   );
  }
  ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
  ierr = m_field.set_global_ghost_vector(
   "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
   ); CHKERRQ(ierr);
 
  //==============================================================================================================================
 
  //Solver
  KSP solver;
  ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
  ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
  ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
  ierr = KSPSetUp(solver); CHKERRQ(ierr);
 
  MatrixDouble Dmat;
  Dmat.resize(6,6);
  Dmat.clear();
 
  //calculate honmogenised stress
  //create a vector for 6 components of homogenized stress
  Vec stress_homo;
  int stress_homo_ghost[] = { 0,1,2,3,4,5,6 };
  ierr = VecCreateGhost(
  PETSC_COMM_WORLD,(!m_field.get_comm_rank())?6:0,6,6,stress_homo_ghost,&stress_homo
  );  CHKERRQ(ierr);
 
  lagrangian_element.setRVEBCsHomoStressOperators("DISPLACEMENT","LAGRANGE_MUL_DISP","MESH_NODE_POSITIONS",stress_homo);
 
  PetscScalar *avec;
  ierr = VecGetArray(stress_homo,&avec); CHKERRQ(ierr);
  for(int ii = 0;ii<6;ii++) {
    ierr = VecZeroEntries(D); CHKERRQ(ierr);
    ierr = KSPSolve(solver,F[ii],D); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(D,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = m_field.set_global_ghost_vector(
     "ELASTIC_MECHANICS",ROW,D,INSERT_VALUES,SCATTER_REVERSE
     ); CHKERRQ(ierr);
 
    //post-processing the resutls
    post_proc.setDoPreProcess();
    post_proc.unSetDoPostProcess();
    ierr = m_field.loop_finite_elements(
            "ELASTIC_MECHANICS","ELASTIC",post_proc
            ); CHKERRQ(ierr);
    post_proc.unSetDoPreProcess();
    post_proc.setDoPostProcess();
    {
      ostringstream sss;
      sss << "mode_" << "Disp" << "_" << ii << ".h5m";
      rval = post_proc.postProcMesh.write_file(
             sss.str().c_str(),"MOAB","PARALLEL=WRITE_PART"
             ); CHKERRQ_MOAB(rval);
    }
 
    ierr = VecZeroEntries(stress_homo); CHKERRQ(ierr);
    ierr = m_field.loop_finite_elements(
    "ELASTIC_MECHANICS","LAGRANGE_ELEM",lagrangian_element.getLoopFeRVEBCStress()
    ); CHKERRQ(ierr);
    ierr = PetscOptionsGetReal(
     PETSC_NULL,"-my_rve_volume",&rve_volume,PETSC_NULL
     ); CHKERRQ(ierr);
    ierr = VecAssemblyBegin(stress_homo); CHKERRQ(ierr);
    ierr = VecAssemblyEnd(stress_homo); CHKERRQ(ierr);
    ierr = VecScale(stress_homo,1.0/rve_volume); CHKERRQ(ierr);
    ierr = VecGhostUpdateBegin(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    ierr = VecGhostUpdateEnd(stress_homo,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
    for(int jj=0; jj<6; jj++) {
      Dmat(jj,ii) = avec[jj];
    }
  }
 
  ierr = VecRestoreArray(stress_homo,&avec); CHKERRQ(ierr);
  PetscPrintf(
              PETSC_COMM_WORLD,"\nHomogenised Stiffens Matrix = \n\n"
              );
 
  for(int ii=0; ii<6; ii++) {
    PetscPrintf(
    PETSC_COMM_WORLD,
    "stress %d\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\t\t%8.5e\n",
    ii,Dmat(ii,0),Dmat(ii,1),Dmat(ii,2),Dmat(ii,3),Dmat(ii,4),Dmat(ii,5)
    );
  }
  //==============================================================================================================================
 
 
 
  //Open mesh_file_name.txt for writing
  ofstream myfile;
  myfile.open ((string(mesh_file_name)+".txt").c_str());
 
  //Output displacements
  myfile << "<<<< Homonenised stress >>>>>" << endl;
 
  if(pcomm->rank()==0){
    for(int ii=0; ii<6; ii++){
      myfile << boost::format("%.3lf") % roundn(Dmat(ii,0)) << endl;
      avec++;
    }
  }
  //Close mesh_file_name.txt
  myfile.close();
 
 
 
 
 
  //detroy matrices
  for(int ii = 0;ii<6;ii++) {
    ierr = VecDestroy(&F[ii]); CHKERRQ(ierr);
  }
  ierr = VecDestroy(&D); CHKERRQ(ierr);
  ierr = MatDestroy(&Aij); CHKERRQ(ierr);
  ierr = KSPDestroy(&solver); CHKERRQ(ierr);
  ierr = VecDestroy(&stress_homo); CHKERRQ(ierr);
 
  ierr = PetscTime(&v2);CHKERRQ(ierr);
  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
 
  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f CPU Time = %f\n",pcomm->rank(),v2-v1,t2-t1);
 
  }
  CATCH_ERRORS;
 
 
  MoFEM::Core::Finalize();
 
}