/*****************************************************************************
 * $CAMITK_LICENCE_BEGIN$
 *
 * CamiTK - Computer Assisted Medical Intervention ToolKit
 * (c) 2001-2016 Univ. Grenoble Alpes, CNRS, TIMC-IMAG UMR 5525 (GMCAO)
 *
 * Visit http://camitk.imag.fr for more information
 *
 * This file is part of CamiTK.
 *
 * CamiTK is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License version 3
 * only, as published by the Free Software Foundation.
 *
 * CamiTK 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 version 3 for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * version 3 along with CamiTK.  If not, see <http://www.gnu.org/licenses/>.
 *
 * $CAMITK_LICENCE_END$
 ****************************************************************************/

#include <stdexcept>
#include <iostream>
#include <fstream>
#include <iomanip>


// Other includes
#include "PhysicalModel.h"
//#include "Object3D.h"
#include "Atom.h"
#include "Cell.h"
#include "CellProperties.h"
#include "MultiComponent.h"
#include "StructuralComponent.h"

#include "PhysicalModelVersion.h"

// pmlschema stuffs
#include <PhysicalModel.hxx>
#include <StructuralComponent.hxx>
#include <InformativeComponent.hxx>
#include <Atoms.hxx>
#include <MultiComponent.hxx>

// XercesC stuffs
#include <xercesc/util/PlatformUtils.hpp>

//--------------- Constructor/Destructor ------------------------------
PhysicalModel::PhysicalModel() {
    init();
}

PhysicalModel::PhysicalModel ( const char * fileName, PtrToSetProgressFunction pspf ) throw ( PMLAbortException ) {
    init();
    setProgressFunction = pspf;

    // load from the xml file
    xmlRead ( fileName );
}

// --------------- destructor ---------------
PhysicalModel::~PhysicalModel() {
    clear();
}

// --------------- init ---------------
void PhysicalModel::init() {
    properties = new Properties ( this );
    positionPtr = NULL;
    exclusiveComponents = NULL;
    informativeComponents = NULL;
    atoms = NULL;
    setProgressFunction = NULL;
    cellIndexOptimized = true; //always hopeful!
    isModifiedFlag = false;
}

// --------------- setProgress ---------------
void PhysicalModel::setProgress ( const float donePercentage ) {
    if ( setProgressFunction != NULL ) {
        setProgressFunction ( donePercentage );
    }
}

// --------------- clear ---------------
void PhysicalModel::clear() {
    if ( informativeComponents ) {
        delete informativeComponents;
    }

    informativeComponents = NULL;

    if ( exclusiveComponents ) {
        delete exclusiveComponents;
    }

    exclusiveComponents = NULL;

    // delete all atoms at the end (otherwise deletion of cells in informative
    // or exclusive components will try to tell already deleted (unaccessible) atoms
    // that some of their SC is to be removed from their list
    if ( atoms ) {
        delete atoms;
    }

    atoms = NULL;

    // clean position memory
    delete [] positionPtr;

    // reset all the unique indexes
    AtomProperties::resetUniqueIndex();

    CellProperties::resetUniqueIndex();

    atomMap.clear();

    cellMap.clear();

    delete properties;
    properties = NULL;
}

// --------------- getNumberOfCells   ---------------
unsigned int PhysicalModel::getNumberOfCells() const {
    unsigned int nrOfCells = 0;

    if ( exclusiveComponents ) {
        nrOfCells += exclusiveComponents->getNumberOfCells();
    }

    if ( informativeComponents ) {
        nrOfCells += informativeComponents->getNumberOfCells();
    }

    return nrOfCells;
}

// --------------- getPositionPointer ---------------
double * PhysicalModel::getPositionPointer() const {
    return positionPtr;
}

double * PhysicalModel::getPositionPointer ( const Atom *a ) const {
    unsigned int idInAtom = a->getIndexInAtoms();

    if ( idInAtom>=0 )
        // * 3 because the position are stored consequently, using a double for x, y and z (3 double then!)
    {
        return positionPtr + idInAtom*3;
    } else {
        // look for atom  #index in the atoms
        Atom *ai = NULL;
        unsigned int i = 0;

        while ( i<atoms->getNumberOfStructures() && ( ai != a ) )  {
            ai = dynamic_cast<Atom *> ( atoms->getStructure ( i ) );
            i++;
        }

        if ( ai == a )
            // the memory allocated to atom #index is (i-1) * 3 (because the position are stored consequently, using a double for x, y and z
        {
            return positionPtr + ( i-1 ) *3;
        } else {
            return NULL;
        }
    }
}

double * PhysicalModel::getPositionPointer ( const unsigned int index ) const {
    // look for atom  #index in the atoms
    bool found = false;
    unsigned int i = 0;

    while ( i<atoms->getNumberOfStructures() && !found ) {
        Atom *a = dynamic_cast<Atom *> ( atoms->getStructure ( i ) );
        found = ( a->getIndex() == index );
        i++;
    }

    if ( found )
        // the memory allocated to atom #index is (i-1) * 3 (because the position are stored consequently, using a double for x, y and z
    {
        return positionPtr + ( i-1 ) *3;
    } else {
        return NULL;
    }
}

// --------------- optimizeIndexes ---------------
void PhysicalModel::optimizeIndexes ( MultiComponent * mc, unsigned int * index ) {
    Component *c;
    StructuralComponent *sc;
    Cell *cell;

    for ( unsigned int i = 0; i < mc->getNumberOfSubComponents(); i++ ) {
        c = mc->getSubComponent ( i );

        if ( c->isInstanceOf ( "MultiComponent" ) ) {
            optimizeIndexes ( ( MultiComponent * ) c, index );
        } else {
            if ( c->isInstanceOf ( "StructuralComponent" ) ) {
                sc = ( StructuralComponent * ) c;
                // check all cells

                for ( unsigned int j = 0; j < sc->getNumberOfStructures(); j++ ) {
                    if ( sc->getStructure ( j )->isInstanceOf ( "Cell" ) ) {
                        cell = ( Cell * ) sc->getStructure ( j );
                        // if this is the sc that make the cell print its data, change cell index

                        if ( cell->makePrintData ( sc ) ) {
                            cell->setIndex ( *index );
                            *index = ( *index ) + 1;
                        }
                    }
                }
            }
        }
    }
}

void PhysicalModel::optimizeIndexes() {
    // to optimize the indexes: do as if it was a print/read operation (same order
    // and change the cell index (as everyone is linked with ptrs, that should not
    // change anything else

    // first: the atoms
    if ( atoms ) {
        for ( unsigned int i = 0; i < atoms->getNumberOfStructures(); i++ ) {
            ( ( Atom * ) atoms->getStructure ( i ) )->setIndex ( i );
        }
    }

    // then the cells
    unsigned int newIndex =  0;

    if ( exclusiveComponents ) {
        optimizeIndexes ( exclusiveComponents, &newIndex );
    }

    if ( informativeComponents ) {
        optimizeIndexes ( informativeComponents, &newIndex );
    }

}

// --------------- xmlPrint   ---------------
void PhysicalModel::xmlPrint ( std::ostream &o, bool opt ) {

    // should we optimize the cell indexes ?
    if ( !cellIndexOptimized && opt ) {
        optimizeIndexes();
    }

    // print out the whole thing
    o << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << std::endl;

    o << "<!-- physical model (PML) is a generic representation for 3D physical model." << std::endl
      << "     PML supports not continous indexes and multiple non-exclusive labelling." << std::endl
      << "  --> " << std::endl;

    o << "<physicalModel";

    if ( getName() != "" ) {
        o << " name=\"" << getName().c_str() << "\"";
    }

    if ( atoms ) {
        o << " nrOfAtoms=\"" << atoms->getNumberOfStructures() << "\"" << std::endl;
    }

    if ( exclusiveComponents ) {
        o << " nrOfExclusiveComponents=\"" << exclusiveComponents->getNumberOfSubComponents() << "\"" << std::endl;
    }

    if ( informativeComponents ) {
        o << " nrOfInformativeComponents=\"" << informativeComponents->getNumberOfSubComponents() << "\"" << std::endl;
    }

    o << " nrOfCells=\"" << getNumberOfCells() << "\"" << std::endl;

    for ( unsigned int i = 0; i<properties->numberOfFields(); i++ ) {
        o << " " <<  properties->getField ( i ) <<"=\"" << properties->getString ( properties->getField ( i ) ) << "\"" << std::endl;
    }

    o << ">" << std::endl;

    o << "<!-- list of atoms: -->" << std::endl;

    o << "<atoms>" << std::endl;

    if ( atoms ) {
        atoms->xmlPrint ( o );
    }

    o << "</atoms>" << std::endl;

    o << "<!-- list of exclusive components : -->" << std::endl;

    o << "<exclusiveComponents>" << std::endl;

    if ( exclusiveComponents ) {
        exclusiveComponents->xmlPrint ( o );
    }

    o << "</exclusiveComponents>" << std::endl;

    if ( informativeComponents ) {
        o << "<!-- list of informative components : -->" << std::endl;
        o << "<informativeComponents>" << std::endl;
        informativeComponents->xmlPrint ( o );
        o << "</informativeComponents>" << std::endl;
    }

    o << "</physicalModel>" << std::endl;

    // serialized/marshalled -> saved
    isModifiedFlag = false;
}

// --------------- xmlRead   ---------------
void PhysicalModel::xmlRead ( const char * filename ) throw ( PMLAbortException ) {
    // clear all the current data
    clear();

    // Set the locale to C for using dot as decimal point dispite locale
    // Set utf8 for output to enforce using utf8 strings.
    char * statusOk = setlocale ( LC_CTYPE, "C.UTF-8" );
    if ( statusOk!=NULL ) {
        statusOk = setlocale ( LC_NUMERIC, "C.UTF-8" );
    }
    if ( statusOk!=NULL ) {
        statusOk = setlocale ( LC_TIME, "C.UTF-8" );
    }

    // french is: "fr_FR.UTF8"
    if ( !statusOk ) {
        // try without UTF-8
        statusOk = setlocale ( LC_CTYPE,"C" );
        if ( statusOk!=NULL ) {
            statusOk = setlocale ( LC_NUMERIC, "C" );
        }
        if ( statusOk!=NULL ) {
            statusOk = setlocale ( LC_TIME, "C" );
        }
        if ( statusOk==NULL ) {
            std::cerr << "Could not set the locale to C. This is mandatory to enforce using dot as decimal separator (platform independency)." << std::endl;
            std::cerr << "This can cause a lot of trouble for XML I/O... Beware of decimal dots..." << std::endl;
        }
    }

    // Initialize manually the XercesC++ runtime, required by using wildcards in XSD schema.
    // Implies that we pass the xml_shema::flags::dont_initialize to any serialization method of the pmlschema library
    // as we already have initialized the XercesC runtime manually.
    xercesc::XMLPlatformUtils::Initialize ();

    // Use XSD with the pmlschema library to read the loads described in the xml file.
    try {
        std::auto_ptr<physicalModel::PhysicalModel> root = physicalModel::physicalModel ( filename, xml_schema::flags::dont_initialize | xml_schema::flags::dont_validate );

        // get the basename
        std::string basename = filename;
        // remove the path
        unsigned lastSeparator = basename.find_last_of ( "/\\" ); // works for unix and windows
        if ( lastSeparator != std::string::npos ) {
            basename = basename.substr ( lastSeparator+1 );
        }
        // remove the extension
        lastSeparator = basename.find_last_of ( "." );
        if ( lastSeparator != std::string::npos ) {
            basename.erase ( lastSeparator );
        }

        // Parse the xml content
        this->parseTree ( root, basename );
    } catch ( const xml_schema::exception& e ) {
        std::ostringstream os;
        os << "Library-pml Error: In PhysicalModel::xmlRead(..): Failed to read the xml file, reason:" << std::endl;
        os << e << std::endl;
        std::cerr << os.str() << std::endl;
        throw PMLAbortException ( os.str () );
    }

    // Terminate the XercesC++ runtime manually
    xercesc::XMLPlatformUtils::Terminate ();
}


// ------------------ parse tree ------------------
bool PhysicalModel::parseTree ( std::auto_ptr<physicalModel::PhysicalModel> root , std::string defaultName ) {

    // Pml file have a name which can be an empty string
    std::string name;
    if ( root->name().present() ) {
        name = root->name().get();
    } else {
        name = defaultName;
    }

    this->properties = new Properties ( this, name );

    // add additionnal attributes as properties
    physicalModel::PhysicalModel::any_attribute_set unknownAttrs = root->any_attribute();
    this->properties->xmlToFields ( unknownAttrs );

    // Parse the atoms
    parseAtoms ( root->atoms() );

    // allocate the big memory bunch
    positionPtr = new double[3*atoms->getNumberOfStructures()];
    // assign all position memory
    double * currentPositionPtr = positionPtr;

    for ( unsigned int i=0; i<atoms->getNumberOfStructures(); i++ ) {
        Atom *a = dynamic_cast<Atom *> ( atoms->getStructure ( i ) );
        a->getProperties()->setPositionPointer ( currentPositionPtr );
        // next position => jump 3 double memory space
        currentPositionPtr = currentPositionPtr + 3;
    }

    // Parse the exclusive components
    physicalModel::ExclusiveComponent ecs = root->exclusiveComponents();
    physicalModel::MultiComponent xmlExclusiveMC = ecs.multiComponent();
    // Get the top level multicomponent node (with name + unkown attributes)
    MultiComponent* exclusiveMC = new MultiComponent ( this );
    // recursively build the tree of exclusive components
    parseComponents ( xmlExclusiveMC, exclusiveMC, true );
    // Add it to the structure
    this->setExclusiveComponents ( exclusiveMC );

    // Parse the informative components
    if ( root->informativeComponents().present() ) {
        physicalModel::InformativeComponent ics = root->informativeComponents().get();
        physicalModel::MultiComponent xmlInformativeMC = ics.multiComponent();
        // Get the top level multicomponent node (with name + unkown attributes)
        MultiComponent* informativeMC = new MultiComponent ( this );
        // recursively build the tree of informative components
        parseComponents ( xmlInformativeMC, informativeMC, false );

        // Add it to the structure
        this->setInformativeComponents ( informativeMC );
    }

    return true;
}


// ------------------ parse atoms ------------------
bool PhysicalModel::parseAtoms ( physicalModel::PhysicalModel::atoms_type atomsRoot ) {

    // Parse the content of the structuralComponent
    physicalModel::Atoms:: structuralComponent_type xmlSC = atomsRoot.structuralComponent();
    StructuralComponent* sc = new StructuralComponent ( this, xmlSC );

    // Parse the number of structures
    if ( xmlSC.nrOfStructures().present() ) {
        sc->plannedNumberOfStructures ( xmlSC.nrOfStructures().get().value() );
    }

    // Parse the atoms
    unsigned int atomOrderNumber = 0; // index in the atoms SC
    physicalModel::StructuralComponent::atom_sequence &atoms = xmlSC.atom();
    for ( physicalModel::StructuralComponent::atom_iterator atomIt ( atoms.begin() ); atomIt != atoms.end(); atomIt++ ) {
        physicalModel::Atom& currentAtom = *atomIt;
        Atom *newAtom = new Atom ( this, currentAtom, atomOrderNumber );
        atomOrderNumber++;
        sc->addStructure ( newAtom );

    }
    // Identify the parsed structural component as the list of atoms
    this->setAtoms ( sc );

    return true;
}


// ------------------ parse Components ------------------
bool PhysicalModel::parseComponents ( physicalModel::MultiComponent xmlFatherMC, Component* father, bool isExclusive ) {

    // Parse the name and unkown properties attribute of the main multiComponent
    MultiComponent* fatherMC = ( MultiComponent* ) father;
    fatherMC->setName ( xmlFatherMC.name().get() );
    fatherMC->getProperties()->xmlToFields ( xmlFatherMC.any_attribute() );

    // Recursively consider multi component children
    physicalModel::MultiComponent::multiComponent_sequence xmlChildren = xmlFatherMC.multiComponent();
    for ( physicalModel::MultiComponent::multiComponent_iterator MCIt = xmlChildren.begin(); MCIt != xmlChildren.end(); MCIt++ ) {
        physicalModel::MultiComponent xml_child = ( *MCIt );
        MultiComponent* child = new MultiComponent ( this );
        fatherMC->addSubComponent ( child );
        this->parseComponents ( xml_child, child, isExclusive );
    }

    // Consider the structural component children
    physicalModel::MultiComponent::structuralComponent_sequence xmlAllSC = xmlFatherMC.structuralComponent();
    for ( physicalModel::MultiComponent::structuralComponent_iterator SCIt = xmlAllSC.begin(); SCIt != xmlAllSC.end(); SCIt++ ) {
        physicalModel::StructuralComponent xmlSC = ( *SCIt );
        StructuralComponent* sc = new StructuralComponent ( this, xmlSC );
        sc->setExclusive ( isExclusive );
        fatherMC->addSubComponent ( sc );

        // StructuralComponent's cells
        physicalModel::StructuralComponent::cell_sequence xmlAllCells = xmlSC.cell();
        for ( physicalModel::StructuralComponent::cell_iterator  cellIt = xmlAllCells.begin();
                cellIt != xmlAllCells.end();
                cellIt++ ) {
            physicalModel::Cell xmlCell = * ( cellIt );
            // create the cell according to its geometric type
            Cell* cell = new Cell ( this, xmlCell, sc );
            cell->setExclusive ( isExclusive );
            if ( this->cellIndexOptimized ) {
                std::GlobalIndexStructurePair pair ( cell->getIndex(), cell );
                this->addGlobalIndexCellPair ( pair );
            }
            sc->addStructure ( cell, false );
        }

        // StructuralComponent's atoms references
        physicalModel::StructuralComponent::atomRef_sequence xmlAllAtomRefs = xmlSC.atomRef();
        for ( physicalModel::StructuralComponent::atomRef_iterator   atomRefIt = xmlAllAtomRefs.begin();
                atomRefIt != xmlAllAtomRefs.end();
                atomRefIt++ ) {
            physicalModel::AtomRef atomRef = * ( atomRefIt );
            // Get the atom corresponding to the reference
            Atom* a = this->getAtom ( atomRef.index() );
            if ( a ) {
                sc->addStructure ( a );
            } else {
                std::cerr << "PhysicalModel::parseComponents: cannot find atom of ref: " << atomRef.index() << std::endl;
            }
        }
    }

    return true;
}



// ------------------ getComponentByName ------------------
Component * PhysicalModel::getComponentByName ( const std::string n ) {
    //-- look for the component in exclusive component first
    Component * foundC;
    foundC = exclusiveComponents->getComponentByName ( n );

    //-- then look in the informative components
    if ( !foundC && informativeComponents ) {
        foundC = informativeComponents->getComponentByName ( n );
    }
    
    //-- at last, just in case, look if this is not the name of the atoms SC
    if ( !foundC && getAtoms()->getName() == n) {
        foundC = getAtoms();
    }
    
    return foundC;
}

// ----------------------- setAtoms ------------------
void PhysicalModel::setAtoms ( StructuralComponent *sc, bool deleteOld ) {
    Atom *a;

    if ( sc->composedBy() == StructuralComponent::ATOMS ) {
        if ( atoms && deleteOld ) {
            delete atoms;
        }

        atoms = sc;

        // register all the atoms in the map, and tell the atoms about its new status
        for ( unsigned int i = 0; i < sc->getNumberOfStructures(); i++ ) {
            a = ( Atom * ) sc->getStructure ( i );
            a->getProperties()->setPhysicalModel ( this );
            addGlobalIndexAtomPair ( std::GlobalIndexStructurePair ( a->getIndex(), a ) );
        }
    }
}

// ----------------------- addAtom ------------------
bool PhysicalModel::addAtom ( Atom *newA ) {
    // register the atom in the map if possible
    if ( atoms && addGlobalIndexAtomPair ( std::GlobalIndexStructurePair ( newA->getIndex(), newA ) ) ) {
        // add the atom in the atom structural component
        atoms->addStructure ( newA );
        return true;
    } else {
        return false;    // atom does not have a unique index
    }
}

// ----------------------- addGlobalIndexAtomPair ------------------
bool PhysicalModel::addGlobalIndexAtomPair ( std::GlobalIndexStructurePair p ) {
    std::GlobalIndexStructureMapIterator mapIt;
    // check if the atom's index is unique
    mapIt = atomMap.find ( p.first );

    // if the index was found, one can not add the atom

    if ( mapIt != atomMap.end() ) {
        return false;
    }

    // if the atom is present in the map then replace the pair <atomIndex, Atom*>
    mapIt = atomMap.begin();

    while ( mapIt != atomMap.end() && mapIt->second != p.second ) {
        mapIt++;
    }

    // if found then remove the pair
    if ( mapIt != atomMap.end() ) {
        atomMap.erase ( mapIt );
    }

    // insert or re-insert (and return true if insertion was ok)
    return atomMap.insert ( p ).second;
}

// ----------------------- addGlobalIndexCellPair ------------------
bool PhysicalModel::addGlobalIndexCellPair ( std::GlobalIndexStructurePair p ) {
    std::GlobalIndexStructureMapIterator mapIt;
    // check if the cell index is unique
    mapIt = cellMap.find ( p.first );

    // if the index was found, one can not add the cell

    if ( mapIt != cellMap.end() ) {
        return false;
    }

    // if the cell is present in the map then replace the pair <cellIndex, Cell*>
    mapIt = cellMap.begin();

    while ( mapIt != cellMap.end() && mapIt->second != p.second ) {
        mapIt++;
    }

    // if found then remove the pair
    if ( mapIt != cellMap.end() ) {
        cellMap.erase ( mapIt );
    }

    // insert or re-insert
    bool insertionOk = cellMap.insert ( p ).second;

    // is that optimized?
    cellIndexOptimized = cellIndexOptimized && ( ( Cell * ) p.second )->getIndex() == optimizedCellList.size();

    if ( cellIndexOptimized ) {
        optimizedCellList.push_back ( ( Cell * ) p.second );
    }

    // insert or re-insert (and return true if insertion was ok)
    return insertionOk;
}

// ----------------------- setAtomPosition ------------------
void PhysicalModel::setAtomPosition ( Atom *atom, const double pos[3] ) {
    atom->setPosition ( pos );
}

// ----------------------- setExclusiveComponents ------------------
void PhysicalModel::setExclusiveComponents ( MultiComponent *mc ) {
    if ( exclusiveComponents ) {
        delete exclusiveComponents;
    }

    exclusiveComponents = mc;
    mc->setPhysicalModel ( this );
}

// ----------------------- setInformativeComponents ------------------
void PhysicalModel::setInformativeComponents ( MultiComponent *mc ) {
    if ( informativeComponents ) {
        delete informativeComponents;
    }

    informativeComponents = mc;
    mc->setPhysicalModel ( this );
}

// ----------------------- getNumberOfExclusiveComponents ------------------
unsigned int PhysicalModel::getNumberOfExclusiveComponents() const {
    if ( !exclusiveComponents ) {
        return 0;
    } else {
        return exclusiveComponents->getNumberOfSubComponents();
    }
}

// ----------------------- getNumberOfInformativeComponents ------------------
unsigned int PhysicalModel::getNumberOfInformativeComponents() const {
    if ( !informativeComponents ) {
        return 0;
    } else {
        return informativeComponents->getNumberOfSubComponents();
    }
}

// ----------------------- getNumberOfAtoms ------------------
unsigned int PhysicalModel::getNumberOfAtoms() const {
    if ( !atoms ) {
        return 0;
    } else {
        return atoms->getNumberOfStructures();
    }
}

// ----------------------- getExclusiveComponent ------------------
Component * PhysicalModel::getExclusiveComponent ( const unsigned int i ) const {
    if ( !exclusiveComponents ) {
        return 0;
    } else {
        return exclusiveComponents->getSubComponent ( i );
    }
}

// ----------------------- getInformativeComponent ------------------
Component * PhysicalModel::getInformativeComponent ( const unsigned int i ) const {
    if ( !informativeComponents ) {
        return 0;
    } else {
        return informativeComponents->getSubComponent ( i );
    }
}


// ----------------------- exportAnsysMesh ------------------
void PhysicalModel::exportAnsysMesh ( std::string filename ) {
    //--- Writing nodes
    std::ofstream nodeFile;
    nodeFile.open ( ( filename + ".node" ).c_str() );

    if ( !nodeFile.is_open() ) {
        std::cerr << "Error in PhysicalModel::exportAnsysMesh : unable to create .node output file" << std::endl;
        return;
    }

    for ( unsigned int i = 0; i < getNumberOfAtoms(); i++ ) {
        double pos[3];
        //WARNING getAtom(i) do not work if indexes do not follow each others
        Atom* at= ( Atom* ) ( atoms->getStructure ( i ) );
        // WARNING : indexes are in base 1 !!!!
        unsigned int ansysIndex = at->getIndex() + 1;

        // coordinates of this node
        at->getPosition ( pos );

        nodeFile << std::setw ( 8 ) << ansysIndex << " ";
        for ( unsigned int j=0; j<3; j++ ) {
            nodeFile << std::setprecision ( 8 ) << std::setw ( 3 ) << std::fixed << std::scientific << pos[j] << "     ";
        }
        nodeFile << std::endl;
    }

    nodeFile.close();

    //--- Writing elements : exlusive cells
    std::ofstream elemFile;
    elemFile.open ( ( filename + ".elem" ).c_str() );

    if ( !elemFile.is_open() ) {
        std::cerr << "Error in PhysicalModel::exportAnsysMesh : unable to create .elem output file" << std::endl;
        return;
    }

    int MAT, TYPE;
    Component *elements = this->getComponentByName ( "Elements" );
    for ( unsigned int i = 0; i < elements->getNumberOfCells(); i++ ) {

        // get the cell
        Cell * cell = elements->getCell ( i );
        unsigned int ansysIndex = cell->getIndex() + 1;

        switch ( cell->getType() ) {

        case StructureProperties::HEXAHEDRON:

            // I,J,K,L,M,N,O,P,MAT,TYPE,REAL,SECNUM,ESYS,IEL
            //
            // Format hex:
            //   I,J,K,L,M,N,O,P                = indices des noeuds
            //   MAT,TYPE,REAL,SECNUM et ESYS   = attributes numbers
            //   SECNUM                         = beam section number
            //   IEL                            = element number

            MAT  = 1;
            TYPE = 1;

            for ( unsigned int k = 0; k < cell->getNumberOfStructures(); k++ ) {
                elemFile << " " << std::setw ( 5 ) << cell->getStructure ( k )->getIndex() + 1 ;
            }

            elemFile << " " << std::setw ( 5 ) << MAT << " " << std::setw ( 5 ) << TYPE << "     1     1     0 " << std::setw ( 5 ) << ansysIndex << std::endl;
            break;

        case StructureProperties::WEDGE:

            // I,J,K,L,M,N,O,P,MAT,TYPE,REAL,SECNUM,ESYS,IEL
            //
            // Format prism: on repete les noeuds 3 et 7:
            //   I,J,K,K,M,N,O,O                = indices des noeuds
            //   MAT,TYPE,REAL,SECNUM et ESYS   = attributes numbers
            //   SECNUM                         = beam section number
            //   IEL                            = element number
            MAT  = 1;
            TYPE = 1;

            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 0 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 1 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 2 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 2 )->getIndex() + 1;

            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 4 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 5 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 5 )->getIndex() + 1;

            elemFile << " " << std::setw ( 5 ) << MAT << " " << std::setw ( 5 ) << TYPE << "     1     1     0 " << std::setw ( 5 ) << ansysIndex << std::endl;

            break;


        case StructureProperties::TETRAHEDRON:

            MAT  = 1;
            TYPE = 1;

            elemFile << " " << std::setw ( 5 ) <<  cell->getStructure ( 0 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 1 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 2 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;
            elemFile << "     0     0     0     0" ;

            elemFile << " " << std::setw ( 5 ) << MAT << " " << std::setw ( 5 ) << TYPE << "     1     1     0 " << std::setw ( 5 ) << ansysIndex << std::endl;
            break;



        case StructureProperties::QUAD:
            // I,J,K,L,M,N,O,P,MAT,TYPE,REAL,SECNUM,ESYS,IEL
            //
            // Format quad:
            //   I,J,K,K,L,L,L,L                = indices des noeuds
            //   MAT,TYPE,REAL,SECNUM et ESYS   = attributes numbers
            //   SECNUM                         = beam section number
            //   IEL                            = element number

            MAT  = 1;
            TYPE = 1;

            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 0 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 1 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 2 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 2 )->getIndex() + 1;

            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;
            elemFile << " " << std::setw ( 5 ) << cell->getStructure ( 3 )->getIndex() + 1;

            elemFile << " " << std::setw ( 5 ) << MAT << " " << std::setw ( 5 ) << TYPE << "     1     1     0 " << std::setw ( 5 ) << ansysIndex << std::endl;
            break;


        default:
            std::cerr << "PhysicalModel::exportAnsysMesh : unknown type for cell " << cell->getIndex() + 1 << ", neither HEXAHEDRON, WEDGE, THETRAHEDRON nor QUAD. Cant' export in Patran format." << std::endl;
            continue;
        }

    }

    elemFile.close();
}


// ----------------------- exportPatran ------------------
void PhysicalModel::exportPatran ( std::string filename ) {
    std::ofstream outputFile;
    outputFile.open ( filename.c_str() );
    if ( !outputFile.is_open() ) {
        std::cerr << "Error in PhysicalModel::exportPatran : unable to create output file" << std::endl;
        return;
    }

    //--- patran header -> mostly useless info in our case...
    outputFile << "25       0       0       1       0       0       0       0       0\n";
    outputFile << "PATRAN File from: " << getName().c_str() << std::endl;
    outputFile << "26       0       0       1   " << this->getNumberOfAtoms() << "    " << this->getExclusiveComponent ( 0 )->getNumberOfCells() << "       3       4      -1\n";
    outputFile << "17-Mar-00   08:00:00         3.0\n";


    //--- Nodes (atoms)

    for ( unsigned int i = 0; i < this->getNumberOfAtoms(); i++ ) {
        double pos[3];

        // first line
        // WARNING : indexes are in base 1 !!!!
        outputFile << " 1" << std::setw ( 8 ) << getAtom ( i )->getIndex() + 1 << "       0       2       0       0       0       0       0\n";

        // coordinates of this node
        //  fscanf(inputFile, "%d %f %f %f", &j, &x, &y, &z);
        getAtom ( i )->getPosition ( pos );

        // second line : node coordinates
        for ( unsigned int j=0; j<3; j++ ) {
            outputFile << std::setprecision ( 8 ) << std::setw ( 16 ) << std::fixed << std::scientific << pos[j];
        }
        outputFile << " " << std::endl;

        // third line : ??
        outputFile << "1G       6       0       0  000000\n";
    }



    //--- Elements : exlusive cells
    for ( unsigned int i = 0; i < this->getExclusiveComponent ( 0 )->getNumberOfCells(); i++ ) {
        int typeElement;

        // get the cell
        Cell * cell = this->getExclusiveComponent ( 0 )->getCell ( i );

        switch ( cell->getType() ) {

        case StructureProperties::HEXAHEDRON:
            typeElement = 8;
            break;

        case StructureProperties::WEDGE:
            typeElement = 7;
            break;

        default:
            std::cerr << "PhysicalModel::exportPatran : unknown type for cell " << cell->getIndex() + 1 << ", neither HEXAHEDRON nor WEDGE. Cant' export in Patran format." << std::endl;
            continue;
        }

        // first element line
        outputFile << " 2" << std::setw ( 8 ) <<  cell->getIndex() + 1 << std::setw ( 8 ) << typeElement << "       2       0       0       0       0       0\n";

        // second element line
        outputFile << std::setw ( 8 ) << cell->getNumberOfStructures() << "       0       1       0 0.000000000E+00 0.000000000E+00 0.000000000E+00\n";

        // third element line : list of nodes
        for ( unsigned int k = 0; k < cell->getNumberOfStructures(); k++ ) {
            outputFile << std::setw ( 8 ) << cell->getStructure ( k )->getIndex() + 1;
        }

        outputFile << "\n";
    }

    //--- final line
    outputFile << "99       0       0       1       0       0       0       0       0\n" ;

    outputFile.close();
}