meshKCL.cpp

#include "Mesh.hpp"
#include "meshKCL.hpp"
#include "Chrono.hpp"
#include<stdio.h>
#include<stdlib.h>
#include<iostream>
#include<fstream>
#include<limits>
#include <math.h>
#include<map>


//empty constructor
MeshKCL::MeshKCL()
:Mesh()
{}


MeshKCL::MeshKCL(Mesh & mesh)
:Mesh(mesh)
{}


MeshKCL::MeshKCL(std::string mfilename)
:Mesh()
{
    readMeshFromFile(mfilename);
}


MeshKCL::~MeshKCL()
{}


void MeshKCL::initFromFile(std::string mfilename)
{
    clearMesh();
    readMeshFromFile(mfilename);
}


void MeshKCL::readMeshFromFile(std::string meshFileName)
{
    std::string line;
    Chrono chrono;
    chrono.start();
    read_points(meshFileName);
    read_elems(meshFileName);
    finalizeMesh();
    chrono.stop();
    std::cout<<"mesh read in  "<<chrono<<std::endl;
}


// Reads the file of the nodes
void MeshKCL::read_points(const std::string & meshFileName)
{
    Chrono chrono;
    std::string pts_filename=meshFileName+".pts";
    std::ifstream meshFile(pts_filename.c_str());
    
    if(meshFile.fail())
    {
        std::cerr<<" FILE "<<pts_filename<<" NOT FOUND OR CORRUPTED"<<std::endl;
        exit(1);
    }
  
    chrono.start();
    std::cout<<"Reading nodes"<<std::endl;
    meshFile>>_num_Points;
    _points.resize(_num_Points);
    for(size_t iPt=0; iPt<_num_Points; iPt++)
    {
        meshFile>>(_points[iPt]).coord[0];
        meshFile>>(_points[iPt]).coord[1];
        meshFile>>(_points[iPt]).coord[2];
    }
    meshFile.close();
    chrono.stop();
    std::cout<<"nodes read in  "<<chrono<<std::endl;
}


// Reads the file of the elements
void MeshKCL::read_elems(const std::string & meshFileName)
{    
    std::string line;
    Chrono chrono;
    std::string elem_filename=meshFileName+".elem";
    std::ifstream meshFile(elem_filename.c_str());
    
    if(meshFile.fail())
    {
        std::cerr<<" FILE "<<elem_filename<<" NOT FOUND OR CORRUPTED"<<std::endl;
        exit(1);
    }
    size_t nbElem=0;
    chrono.start();
    std::cout<<"Reading elements"<<std::endl;
    meshFile>>nbElem;
    
    _num_Edges=0;
    _num_Trias=0;
    _num_Quads=0;
    _num_Tetras=0;
    _num_Hexas=0;
    size_t num_unk_elem=0;
    
    _edges.resize(nbElem);
    _triangles.resize(nbElem);
    _quadrilaterals.resize(nbElem);
    _tetrahedras.resize(nbElem);
    _hexahedras.resize(nbElem);

    for(size_t iElem=0; iElem<nbElem; iElem++)
    {
      std::string TypeOfElem;
      meshFile>>TypeOfElem;
      
      switch(elemCode(TypeOfElem))
      {
        // Edges
        case(3):
        {
            _num_Edges=_num_Edges+1;
            meshFile>>(_edges[_num_Edges-1]).vertices(0);
            meshFile>>(_edges[_num_Edges-1]).vertices(1);
            meshFile>>(_edges[_num_Edges-1]).label();
            break;
        }
        // Trianlges
        case(5):
        {
            _num_Trias=_num_Trias+1;
            meshFile>>(_triangles[_num_Trias-1]).vertices(0);
            meshFile>>(_triangles[_num_Trias-1]).vertices(1);
            meshFile>>(_triangles[_num_Trias-1]).vertices(2);
            meshFile>>(_triangles[_num_Trias-1]).label();        
            break;
        }
        // Quadrilaterals
        case(9):
        {
            _num_Quads=_num_Quads+1;
            meshFile>>(_quadrilaterals[_num_Quads-1]).vertices(0);
            meshFile>>(_quadrilaterals[_num_Quads-1]).vertices(1);
            meshFile>>(_quadrilaterals[_num_Quads-1]).vertices(2);
            meshFile>>(_quadrilaterals[_num_Quads-1]).vertices(3);
            meshFile>>(_quadrilaterals[_num_Quads-1]).label();
            break;
        }        
        // Tetrahedra
        case(10):
        {
            _num_Tetras=_num_Tetras+1;
            meshFile>>(_tetrahedras[_num_Tetras-1]).vertices(0);
            meshFile>>(_tetrahedras[_num_Tetras-1]).vertices(1);
            meshFile>>(_tetrahedras[_num_Tetras-1]).vertices(2);
            meshFile>>(_tetrahedras[_num_Tetras-1]).vertices(3);
            meshFile>>(_tetrahedras[_num_Tetras-1]).label();
            break;
        }
        // Hexahedra
        case(12):
        {
            _num_Hexas=_num_Hexas+1;        
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(0);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(1);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(2);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(3);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(4);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(5);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(6);
            meshFile>>(_hexahedras[_num_Hexas-1]).vertices(7);
            meshFile>>(_hexahedras[_num_Hexas-1]).label();
            break;
        }
        //unknown elements
        default:
        {
            num_unk_elem=num_unk_elem+1;
            meshFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n' );
            break;
        }
      }//end switch
    }//end for
    meshFile.close();
    std::cout<<num_unk_elem<<" unknown elements"<<std::endl;

    if(_num_Edges<nbElem)
    {
      size_t delta=nbElem-_num_Edges;
      _edges.erase(_edges.begin()+_num_Edges,_edges.begin()+_num_Edges+delta);
    }
    if(_num_Trias<nbElem)
    {
      size_t delta=nbElem-_num_Trias;
      _triangles.erase(_triangles.begin()+_num_Trias,_triangles.begin()+_num_Trias+delta);
    }
    if(_num_Quads<nbElem)
    {
      size_t delta=nbElem-_num_Quads;
      _quadrilaterals.erase(_quadrilaterals.begin()+_num_Quads,_quadrilaterals.begin()+_num_Quads+delta);
    }
    if(_num_Tetras<nbElem)
    {
      size_t delta=nbElem-_num_Tetras;
      _tetrahedras.erase(_tetrahedras.begin()+_num_Tetras,_tetrahedras.begin()+_num_Tetras+delta);
    }
    
    if(_num_Hexas<nbElem)
    {
      size_t delta=nbElem-_num_Hexas;
      _hexahedras.erase(_hexahedras.begin()+_num_Hexas,_hexahedras.begin()+_num_Hexas+delta);
    }
    
    chrono.stop();
    std::cout<<"elements read in  "<<chrono<<std::endl;
}


const unsigned char MeshKCL::elemCode(const std::string  & elemName)
{
    //codes are thos used in vtk
    static std::map<std::string, unsigned char> elem_to_code;
    elem_to_code["Cx"]=3; // edges
    elem_to_code["Tr"]=5; // Triangles
    elem_to_code["Qd"]=9; // Quadrilateral
    elem_to_code["Tt"]=10; // TetraHedra
    elem_to_code["Hx"]=12; // Hexahedra
    elem_to_code["Py"]=14; // Pyramids
    elem_to_code["Pr"]=13; // Prisms
    
    return(elem_to_code[elemName]);
}