t3dinterface.C

Go to the documentation of this file.

/*
 *
 *                 #####    #####   ######  ######  ###   ###
 *               ##   ##  ##   ##  ##      ##      ## ### ##
 *              ##   ##  ##   ##  ####    ####    ##  #  ##
 *             ##   ##  ##   ##  ##      ##      ##     ##
 *            ##   ##  ##   ##  ##      ##      ##     ##
 *            #####    #####   ##      ######  ##     ##
 *
 *
 *             OOFEM : Object Oriented Finite Element Code
 *
 *               Copyright (C) 1993 - 2013   Borek Patzak
 *
 *
 *
 *       Czech Technical University, Faculty of Civil Engineering,
 *   Department of Structural Mechanics, 166 29 Prague, Czech Republic
 *
 *  This library 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.1 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 library; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "t3dinterface.h"
#include "errorestimator.h"
#include "domain.h"
#include "node.h"
#include "element.h"
#include "engngm.h"
#include "remeshingcrit.h"
#include "oofemtxtdatareader.h"
#include "dynamicinputrecord.h"

#include "crosssection.h"
#include "classfactory.h"

#include "nonlocalbarrier.h"
#include "initialcondition.h"
#include "classfactory.h"
//#include "loadtimefunction.h"
#include "function.h"
#include "outputmanager.h"

#include <cmath>

namespace oofem {
REGISTER_Mesher(T3DInterface, MPT_T3D);


MesherInterface :: returnCode
T3DInterface :: createMesh(TimeStep *tStep, int domainNumber, int domainSerNum, Domain **dNew)
{
    * dNew = NULL;
    if ( this->createInput(this->domain, tStep) ) {
        return MI_NEEDS_EXTERNAL_ACTION;
    } else {
        return MI_FAILED;
    }
}

int
T3DInterface :: createInput(Domain *d, TimeStep *tStep)
{
    int nnodes = d->giveNumberOfDofManagers(), nelem = d->giveNumberOfElements();
    double density;
    FILE *outputStrem;
    Node *inode;
    Element *ielem;
    int edges, trias, quads, tetras, pyrams, wedges, hexas;
    IntArray edgeIdArray, triaIdArray, quadIdArray, tetraIdArray, pyramIdArray, wedgeIdArray, hexaIdArray;
    bool tri_tetra = false;
    char fileName [ 32 ];

    edges = trias = quads = tetras = pyrams = wedges = hexas = 0;
    for ( int i = 1; i <= nelem; i++ ) {
        ielem = d->giveElement(i);
        switch ( ielem->giveGeometryType() ) {
        case EGT_point:
            break;
        case EGT_line_1:
        case EGT_line_2:
            edges++;
            break;
        case EGT_triangle_1:
        case EGT_triangle_2:
            trias++;
            break;
        case EGT_quad_1:
        case EGT_quad_2:
            quads++;
            break;
        case EGT_tetra_1:
        case EGT_tetra_2:
            tetras++;
            break;
        case EGT_wedge_1:
        case EGT_wedge_2:
            wedges++;
            break;
        case EGT_hexa_1:
        case EGT_hexa_2:
            hexas++;
            break;
        default:
            OOFEM_ERROR( "unknown element type (%s)",
                         __Element_Geometry_TypeToString( ielem->giveGeometryType() ) );
        }
    }

    edgeIdArray.resize(edges);
    triaIdArray.resize(trias);
    quadIdArray.resize(quads);
    tetraIdArray.resize(tetras);
    pyramIdArray.resize(pyrams);
    wedgeIdArray.resize(wedges);
    hexaIdArray.resize(hexas);

    edges = trias = quads = tetras = pyrams = wedges = hexas = 0;
    for ( int i = 1; i <= nelem; i++ ) {
        ielem = d->giveElement(i);
        switch ( ielem->giveGeometryType() ) {
        case EGT_point:
            break;
        case EGT_line_1:
        case EGT_line_2:
            edgeIdArray.at(++edges) = i;
            break;
        case EGT_triangle_1:
        case EGT_triangle_2:
            triaIdArray.at(++trias) = i;
            break;
        case EGT_quad_1:
        case EGT_quad_2:
            quadIdArray.at(++quads) = i;
            break;
        case EGT_tetra_1:
        case EGT_tetra_2:
            tetraIdArray.at(++tetras) = i;
            break;
        case EGT_wedge_1:
        case EGT_wedge_2:
            wedgeIdArray.at(++wedges) = i;
            break;
        case EGT_hexa_1:
        case EGT_hexa_2:
            hexaIdArray.at(++hexas) = i;
            break;
        default:
            break;
        }
    }

    if ( quads + hexas + pyrams + wedges == 0 ) {
        tri_tetra = true;
    }

    if ( d->giveEngngModel()->isParallel() ) {
        sprintf( fileName, "%s.%d", BMF_FILENAME, d->giveEngngModel()->giveRank() );
    } else {
        sprintf(fileName, "%s", BMF_FILENAME);
    }

    outputStrem = fopen(fileName, "w");
    if ( tri_tetra == true ) {
        fprintf(outputStrem, "3 1 -1\n");
        fprintf(outputStrem, "%d %d %d %d\n", nnodes, edges, trias, tetras);
    } else {
        fprintf(outputStrem, "7 1 -1\n");
        fprintf(outputStrem, "%d %d %d %d %d %d %d %d\n", nnodes, edges, trias, quads, tetras, pyrams, wedges, hexas);
    }

    // loop over nodes
    for ( int i = 1; i <= nnodes; i++ ) {
        density = d->giveErrorEstimator()->giveRemeshingCrit()->giveRequiredDofManDensity(i, tStep);
        inode = d->giveNode(i);
        fprintf(outputStrem, "%d %e %e %e  %e\n", i, inode->giveCoordinate(1), inode->giveCoordinate(2), inode->giveCoordinate(3), density);
    }

    // loop separately for each type of element
    // since T3ds support only linear elements in bg mesh, only corner nodes are stored
    // alternatively quadratic elements may be splitted to linear (not supported now)

    if ( edges != 0 ) {
        for ( int i = 1; i <= edges; i++ ) {
            ielem = d->giveElement( edgeIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 2; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( trias != 0 ) {
        for ( int i = 1; i <= trias; i++ ) {
            ielem = d->giveElement( triaIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 3; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( quads != 0 ) {
        for ( int i = 1; i <= quads; i++ ) {
            ielem = d->giveElement( quadIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 4; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( tetras != 0 ) {
        for ( int i = 1; i <= tetras; i++ ) {
            ielem = d->giveElement( tetraIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 4; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( pyrams != 0 ) {
        for ( int i = 1; i <= pyrams; i++ ) {
            ielem = d->giveElement( pyramIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 5; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( wedges != 0 ) {
        for ( int i = 1; i <= wedges; i++ ) {
            ielem = d->giveElement( wedgeIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 6; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    if ( hexas != 0 ) {
        for ( int i = 1; i <= hexas; i++ ) {
            ielem = d->giveElement( hexaIdArray.at(i) );
            fprintf(outputStrem, "%d", i);
            for ( int j = 1; j <= 8; j++ ) {
                fprintf( outputStrem, " %d", ielem->giveNode(j)->giveNumber() );
            }

            fprintf(outputStrem, "\n");
        }
    }

    fclose(outputStrem);

    OOFEM_LOG_INFO("t3d.bmf file created\n");
    return 1;
}





int
T3DInterface :: t3d_2_OOFEM(const char *t3dOutFile, Domain **dNew)
{
  
  std::ifstream inputStream;
  inputStream.open( t3dOutFile );
  if ( !inputStream.is_open() ) {
    OOFEM_ERROR("OOFEMTXTDataReader::OOFEMTXTDataReader: Can't open T3D input stream (%s)", t3dOutFile);
    return 0;
  }
   // create new domain
  ( * dNew ) = new Domain( 2, this->domain->giveSerialNumber() + 1, this->domain->giveEngngModel() );
  ( * dNew )->setDomainType( this->domain->giveDomainType() );

  std::string line;

  //read first line from t3d out file - 4 numbers
  // 2 degree of interpolation
  // other are not important so far 
  std::getline(inputStream, line);
  //convert const char to char in order to use strtok
  char *currentLine = new char[line.size() + 1];
  std::strcpy ( currentLine, line.c_str() );
  // tokenizing line
  char *token = std::strtok(currentLine, " ");
  // set counter to 0
  int i = 0;
  while (token != NULL) {
    if(i == 0)
  {}//int n1 = atoi(token);
    else if (i == 1)
  {}//int interp = atoi(token);
    else if(i == 2)
  {}//int n3 = atoi(token);
    else if(i == 3)
  {}//int n4 = atoi(token);
    else
      break;
    token = std::strtok(NULL, " ");
    i++;
  }  


  int nnodes,ntriangles,ntetras; // nedges,
    
  /*read second line from t3d out file
    4 numbers
    1 - number of nodes
    2 - number of edges
    3 - number of triangles
    4 - number of tetras
  */
  std::getline(inputStream, line);
  //convert const char to char in order to use strtok
  currentLine = new char[line.size() + 1];
  std::strcpy ( currentLine, line.c_str() );
  // tokenizing line
  token = std::strtok(currentLine, " ");
  // set counter to 0
  i = 0;
  while (token != NULL) {
    if(i == 0)
      nnodes = atoi(token);
    else if (i == 1)
  {}//nedges = atoi(token);
    else if(i == 2)
      ntriangles = atoi(token);
    else if(i == 3)
      ntetras = atoi(token);
    else
      break;
    token = std::strtok(NULL, " ");
    i++;
  }  
  // create new domain
  (*dNew)->resizeDofManagers(nnodes);


  //one empty line
  std::getline(inputStream, line);  
  // create nodes
  Node *node;
  // read dofs
  const IntArray dofIDArrayPtr = domain->giveDefaultNodeDofIDArry();
  int ndofs = dofIDArrayPtr.giveSize();
  // loop over number of nodes, read coordinates and create new nodes
  for ( int inode = 1; inode <= nnodes; inode++ ) {
      FloatArray coords(3);
      std::getline(inputStream, line);  
      //convert const char to char in order to use strtok
      currentLine = new char[line.size() + 1];
      std::strcpy ( currentLine, line.c_str() );
      // tokenizing line
      token = std::strtok(currentLine, " ");
      // set counter to 0
      i = 0;
      while (token != NULL) {
  if(i == 0)
      {}//int nodeNum = atoi(token);
  else if (i == 1)
    coords.at(1) = atof(token);
  else if(i == 2)
    coords.at(2) = atof(token);
  else if(i == 3)
    coords.at(3) = atof(token);
  else
    break;
  token = std::strtok(NULL, " ");
  i++;
      }  
      // newly created node
      node = new Node(inode, * dNew);
      //create new node with default DOFs
      node->setNumberOfDofs(ndofs);
      node->setCoordinates( coords );
      // how to set boundary condition ??
      //      node->setBoundaryFlag( mesh->giveNode(inode)->isBoundary() );
      // ( * dNew )->setDofManager(inode, node);

  }//end loop over nodes

  //read empty line 
  std::getline(inputStream, line);  



  Element *parentElementPtr, *elem;
  parentElementPtr = domain->giveElement(1); //??km??
  // loop over triangles, read dofman numbers and create new elements
  for (int itriangle = 1; itriangle <= ntriangles; itriangle++) {
    int elemNumber;
    IntArray dofManagers(3);
    std::getline(inputStream, line); 
    //convert const char to char in order to use strtok
    currentLine = new char[line.size() + 1];
    std::strcpy ( currentLine, line.c_str() );
    // tokenizing line
    token = std::strtok(currentLine, " ");
    // set counter to 0
    i = 0;
    while (token != NULL) {
      if(i == 0)
  elemNumber = atoi(token);
      else if (i >= 1 && i<=4)
  dofManagers.at(i) = atof(token);
      else
  break;
      token = std::strtok(NULL, " ");
      i++;
      } 
    elem = classFactory.createElement(parentElementPtr->giveClassName(), elemNumber, * dNew);
    elem->setDofManagers( dofManagers );
    elem->setMaterial( parentElementPtr->giveMaterial()->giveNumber() );
    elem->setCrossSection( parentElementPtr->giveCrossSection()->giveNumber() );
    //
    // ...
    //

  }//end loop over triangles 

  
// loop over tetras, read dofman numbers and create new elements
  for (int itetra = 1; itetra <= ntetras; itetra++) {
    int elemNumber;
    IntArray dofManagers(4);
    std::getline(inputStream, line); 
    //convert const char to char in order to use strtok
    currentLine = new char[line.size() + 1];
    std::strcpy ( currentLine, line.c_str() );
    // tokenizing line
    token = std::strtok(currentLine, " ");
       // set counter to 0
    i = 0;
    while (token != NULL) {
      if(i == 0)
  elemNumber = atoi(token);
      else if (i >= 1 && i<=4)
  dofManagers.at(i) = atof(token);
      else
  break;
      token = std::strtok(NULL, " ");
      i++;
    } 
    //elem = classFactory.createElement(parentElementPtr->giveClassName(), elemNumber, * dNew);
    elem = classFactory.createElement("LTRSpace", elemNumber, * dNew);
    elem->setDofManagers( dofManagers );
    elem->setMaterial( parentElementPtr->giveMaterial()->giveNumber() );
    elem->setCrossSection( parentElementPtr->giveCrossSection()->giveNumber() );
    (*dNew)->setElement(elemNumber, elem);
  }//end loop over tetras
  
  
  
  CrossSection *crossSection;
  Material *mat;
  NonlocalBarrier *barrier;
  GeneralBoundaryCondition *bc;
  InitialCondition *ic;
  //LoadTimeFunction *ltf;
  Function *ltf;
  std::string name;
  
  // copy of crossections from old domain
  int ncrosssect = domain->giveNumberOfCrossSectionModels();
  ( * dNew )->resizeCrossSectionModels(ncrosssect);
  for ( int i = 1; i <= ncrosssect; i++ ) {
    DynamicInputRecord ir;
    domain->giveCrossSection(i)->giveInputRecord(ir);
    ir.giveRecordKeywordField(name);
    
    crossSection = classFactory.createCrossSection(name.c_str(), i, * dNew);
    crossSection->initializeFrom(&ir);
    ( * dNew )->setCrossSection(i, crossSection);
  }
  
  // copy of materials  from old domain
  int nmat = domain->giveNumberOfMaterialModels();
  ( * dNew )->resizeMaterials(nmat);
  for ( int i = 1; i <= nmat; i++ ) {
    DynamicInputRecord ir;
      domain->giveMaterial(i)->giveInputRecord(ir);
      ir.giveRecordKeywordField(name);
      
      mat = classFactory.createMaterial(name.c_str(), i, * dNew);
      mat->initializeFrom(&ir);
      ( * dNew )->setMaterial(i, mat);
  }
  
  // copy of crossections from old domain
  int nbarriers = domain->giveNumberOfNonlocalBarriers();
  ( * dNew )->resizeNonlocalBarriers(nbarriers);
  for ( int i = 1; i <= nbarriers; i++ ) {
      DynamicInputRecord ir;
      domain->giveNonlocalBarrier(i)->giveInputRecord(ir);
      ir.giveRecordKeywordField(name);
      
      barrier = classFactory.createNonlocalBarrier(name.c_str(), i, * dNew);
      barrier->initializeFrom(&ir);
      ( * dNew )->setNonlocalBarrier(i, barrier);
  }
  
    // copy of boundary conditions from old domain
  int nbc = domain->giveNumberOfBoundaryConditions();
  ( * dNew )->resizeBoundaryConditions(nbc);
  for ( int i = 1; i <= nbc; i++ ) {
    DynamicInputRecord ir;
    domain->giveBc(i)->giveInputRecord(ir);
    ir.giveRecordKeywordField(name);
    
    bc = classFactory.createBoundaryCondition(name.c_str(), i, * dNew);
    bc->initializeFrom(&ir);
    ( * dNew )->setBoundaryCondition(i, bc);
  }
  
  // copy of initial conditions from old domain
  int nic = domain->giveNumberOfInitialConditions();
  ( * dNew )->resizeInitialConditions(nic);
  for ( int i = 1; i <= nic; i++ ) {
    DynamicInputRecord ir;
    domain->giveIc(i)->giveInputRecord(ir);
    ir.giveRecordKeywordField(name);
    
    ic = new InitialCondition(i, *dNew);
    ic->initializeFrom(&ir);
    ( * dNew )->setInitialCondition(i, ic);
  }
  
  // copy of load time functions from old domain
  //  int nltf = domain->giveNumberOfLoadTimeFunctions();
    int nltf = domain->giveNumberOfFunctions();
  // ( * dNew )->resizeLoadTimeFunctions(nltf);
   ( * dNew )->resizeFunctions(nltf);
  for ( int i = 1; i <= nltf; i++ ) {
    DynamicInputRecord ir;
    //domain->giveLoadTimeFunction(i)->giveInputRecord(ir);
    domain->giveFunction(i)->giveInputRecord(ir);
    ir.giveRecordKeywordField(name);
    
    //ltf = classFactory.createLoadTimeFunction(name.c_str(), i, * dNew);
    ltf = classFactory.createFunction(name.c_str(), i, * dNew);
    ltf->initializeFrom(&ir);
    //( * dNew )->setLoadTimeFunction(i, ltf);
    ( * dNew )->setFunction(i, ltf);
  }
  
  // copy output manager settings from old domain
  ( * dNew )->giveOutputManager()->beCopyOf( domain->giveOutputManager() );  
 return 1;  
  
}

// used by HTS element to create export mesh to vtk

int
T3DInterface :: createInput(Element *e, char *t3dInFile)
{ 
  FILE *inputStrem;
  inputStrem = fopen(t3dInFile, "w");
 
  int nnodes = e->giveNumberOfDofManagers();
  // loop over nodes and write vertexes into t3d input file
  for ( int i = 1; i <= nnodes; i++ ) {
    Node *inode;
    inode = e->giveNode(i);
    fprintf(inputStrem, "vertex %d xyz %e %e %e\n", i, inode->giveCoordinate(1), inode->giveCoordinate(2), inode->giveCoordinate(3));
  }
  for (int i = 1; i<=nnodes; i++) {
    if( i< nnodes)
      fprintf(inputStrem, "curve %d order 2 vertex %d %d\n", i, i, i+1);
    else
      fprintf(inputStrem, "curve %d order 2 vertex %d %d\n", i, i, 1);
  }
  fprintf(inputStrem, "patch 1 normal 0 0 1 boundary curve");
  for (int i = 1; i<=nnodes; i++) {
    fprintf(inputStrem, " %d", i);
  }
  fclose(inputStrem);
  
  return 1;

}







int
T3DInterface :: createVTKExportMesh(const char *t3dOutFile,std::vector<FloatArray> &nodeCoords, std::vector<IntArray> &cellNodes, IntArray &cellTypes )
{

   std::ifstream inputStream;
  inputStream.open( t3dOutFile );
  if ( !inputStream.is_open() ) {
    OOFEM_ERROR("OOFEMTXTDataReader::OOFEMTXTDataReader: Can't open T3D input stream (%s)", t3dOutFile);
    return 0;
  }
  std::string line;
  //read and skip first line from t3d out file - 4 numbers
  std::getline(inputStream, line);
  /*read second line from t3d out file
    4 numbers
    1 - number of nodes
    2 - number of edges
    3 - number of triangles
    4 - number of tetras
  */
  int nnodes,ntriangles;// nedges, ntetras
  std::getline(inputStream, line);
  //convert const char to char in order to use strtok
  char *currentLine = new char[line.size() + 1];
  std::strcpy ( currentLine, line.c_str() );
  // tokenizing line
  char *token = std::strtok(currentLine, " ");
  // set counter to 0
  int i = 0;
  while (token != NULL) {
    if(i == 0)
      nnodes = atoi(token);
    else if (i == 1)
  {}//nedges = atoi(token);
    else if(i == 2)
      ntriangles = atoi(token);
    else if(i == 3)
  {}//ntetras = atoi(token);
    else
      break;
    token = std::strtok(NULL, " ");
    i++;
  }  
  delete [] currentLine;
  //read one empty line
  std::getline(inputStream, line);  


 // loop over number of nodes, read coordinates and put them into nodeCoords Array
  for ( int inode = 1; inode <= nnodes; inode++ ) {
      // int nodeNum;
      FloatArray coords(3);
      std::getline(inputStream, line);  
      //convert const char to char in order to use strtok
      currentLine = new char[line.size() + 1];
      std::strcpy ( currentLine, line.c_str() );
      // tokenizing line
      token = std::strtok(currentLine, " ");
      // set counter to 0
      i = 0;
      while (token != NULL) {
  if(i == 0)
      {}// nodeNum = atoi(token);
  else if (i == 1)
    coords.at(1) = atof(token);
  else if(i == 2)
    coords.at(2) = atof(token);
  else if(i == 3)
    coords.at(3) = atof(token);
  else
    break;
  token = std::strtok(NULL, " ");
  i++;
      }  

      nodeCoords.push_back(coords);  
      delete [] currentLine;
  }//end loop over nodes



  //one empty line
  std::getline(inputStream, line);  
  cellTypes.resize(ntriangles);
  // loop over triangles, and fill nodeCoords, cellNodes nad cellTypes arrays
  for (int itriangle = 1; itriangle <= ntriangles; itriangle++) {
    // int elemNumber;
    IntArray dofManagers(3);
    std::getline(inputStream, line); 
    //convert const char to char in order to use strtok
    currentLine = new char[line.size() + 1];
    std::strcpy ( currentLine, line.c_str() );
    // tokenizing line
    token = std::strtok(currentLine, " ");
    // set counter to 0
    i = 0;
    while (token != NULL) {
      if(i == 0)
    {}//elemNumber = atoi(token);
      else if (i >= 1 && i<=3)
  dofManagers.at(i) = atof(token) - 1;
      else
  break;
      token = std::strtok(NULL, " ");
      i++;
    } 
    cellNodes.push_back(dofManagers);
    // in vtk number 5 corresponds to linear triangle
    cellTypes.at(itriangle) = 5;
    delete [] currentLine;
  }//end loop over triangles 
  
  return 1;
}





int
T3DInterface :: createQCInterpolationMesh(const char *t3dOutFile,std::vector<FloatArray> &nodeCoords, std::vector<IntArray> &cellNodes, IntArray &cellTypes )
//
// create interpolation mest from t3dOutFile
// node coordinates and element nodes are saved in "matrix" vector<FloatArray/IntArray> 
//
{

   std::ifstream inputStream;
  inputStream.open( t3dOutFile );
  if ( !inputStream.is_open() ) {
    OOFEM_ERROR("OOFEMTXTDataReader::OOFEMTXTDataReader: Can't open T3D input stream (%s)", t3dOutFile);
    return 0;
  }
  std::string line;
  //read and skip first line from t3d out file - 4 numbers
  std::getline(inputStream, line);
  /*read second line from t3d out file
    4 numbers
    1 - number of nodes
    2 - number of edges
    3 - number of triangles
    4 - number of tetras
  */
  int nnodes,ntriangles,ntetras;//nedges,
  std::getline(inputStream, line);
  //convert const char to char in order to use strtok
  char *currentLine = new char[line.size() + 1];
  std::strcpy ( currentLine, line.c_str() );
  // tokenizing line
  char *token = std::strtok(currentLine, " ");
  // set counter to 0
  int i = 0;
  while (token != NULL) {
    if(i == 0)
      nnodes = atoi(token);
    else if (i == 1)
  {}//nedges = atoi(token);
    else if(i == 2)
      ntriangles = atoi(token);
    else if(i == 3)
      ntetras = atoi(token);
    else
      break;
    token = std::strtok(NULL, " ");
    i++;
  }  
  delete [] currentLine;
  //read one empty line
  std::getline(inputStream, line);  

  // check the number of interpolation element
  if (ntriangles!=0 && ntetras!=0 ) {
    OOFEM_ERROR( "T3DInterface: 2D and 3D interpolation elements are not supported together");  
  } else if (ntriangles!=0) {
    cellTypes.resize(ntriangles);
  } else if (ntetras!=0) {
    cellTypes.resize(ntetras);
  } else { 
  OOFEM_ERROR( "T3DInterface: No interpolation element found in %s", t3dOutFile);  
  } 
  
  // loop over number of nodes, read coordinates and put them into nodeCoords Array
  for ( int inode = 1; inode <= nnodes; inode++ ) {
      FloatArray coords(3);
      // int nodeNum;
      std::getline(inputStream, line);  
      //convert const char to char in order to use strtok
      currentLine = new char[line.size() + 1];
      std::strcpy ( currentLine, line.c_str() );
      // tokenizing line
      token = std::strtok(currentLine, " ");
      // set counter to 0
      i = 0;
      while (token != NULL) {
  if(i == 0)
      {}// nodeNum = atoi(token);
  else if (i == 1)
    coords.at(1) = atof(token);
  else if(i == 2)
    coords.at(2) = atof(token);
  else if(i == 3)
    coords.at(3) = atof(token);
  else
    break;
  token = std::strtok(NULL, " ");
  i++;
      }  

      nodeCoords.push_back(coords);  
      delete [] currentLine;
  }//end loop over nodes

  //one empty line
  std::getline(inputStream, line);  

// loop over triangles, and fill nodeCoords, cellNodes nad cellTypes arrays
  for (int itriangle = 1; itriangle <= ntriangles; itriangle++) {
    IntArray dofManagers(3);
    // int elemNumber;
    std::getline(inputStream, line); 
    //convert const char to char in order to use strtok
    currentLine = new char[line.size() + 1];
    std::strcpy ( currentLine, line.c_str() );
    // tokenizing line
    token = std::strtok(currentLine, " ");
    // set counter to 0
    i = 0;
    while (token != NULL) {
      if(i == 0)
    {}//elemNumber = atoi(token);
      else if (i >= 1 && i<=3)
  dofManagers.at(i) = atof(token);
      else
  break;
      token = std::strtok(NULL, " ");
      i++;
    } 
    cellNodes.push_back(dofManagers);
    // in vtk number 5 corresponds to linear triangle
    cellTypes.at(itriangle) = 5;
    delete [] currentLine;
  }//end loop over triangles 

// loop over tetras, read dofman numbers and fill nodeCoords, cellNodes nad cellTypes arrays
  for (int itetra = 1; itetra <= ntetras; itetra++) {
    IntArray dofManagers(4);
    // int elemNumber;
    std::getline(inputStream, line); 
    //convert const char to char in order to use strtok
    currentLine = new char[line.size() + 1];
    std::strcpy ( currentLine, line.c_str() );
    // tokenizing line
    token = std::strtok(currentLine, " ");
       // set counter to 0
    i = 0;
    while (token != NULL) {
      if(i == 0)
    {}//elemNumber = atoi(token);
      else if (i >= 1 && i<=4)
  dofManagers.at(i) = atof(token);
      else
  break;
      token = std::strtok(NULL, " ");
      i++;
    } 
    cellNodes.push_back(dofManagers);
    // in vtk number 5 corresponds to linear triangle  
    cellTypes.at(itetra) = 5;    // ??km?? TODO: number for tetras in vtk = ???
    delete [] currentLine;
  }//end loop over tetras
 
  
  return 1;
}


  
} // end namespace oofem