geometry/ioMesh/mel.cpp

/*
XLiFE++ is an extended library of finite elements written in C++
    Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 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 General Public License for more details.
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

/*!
  \file mel.cpp
  \author N. Kielbasiewicz
  \since 25 mar 2013
  \date 28 mar 2013

  \brief Implementation of xlifepp::Mesh class members related to input and output mel file format (melina)

  \verbatim
  * comment line
  ----------------------------------------------------------------------------
  TITRE  nbligne(s)
  ...
  ----------------------------------------------------------------------------
  [ FORMAT [DE] LECTURE
    [ [DES] COORDONNEES < ||C*50 ^ '*' > ]
    [ [DE LA] NUMEROTATION [GLOBALE] < ||C*50 ^ '*' > ]
    [ < SANS ^ AVEC > COMMENTAIRE ]
  ]
  ---------------------------- default values --------------------------------
      FORMAT DE LECTURE DES COORDONNEES '6E12.4'
      DE LA NUMEROTATION GLOBALE '18I4'
      AVEC COMMENTAIRE
  ----------------------------------------------------------------------------
  DESCRIPTION GLOBALE [DU] [MAILLAGE]
    [ [NOM] [DES] VARIABLES [D''] ESPACE : + ||C ]
  NOMBRE [D''] ELEMENTS : ||I ]
  --------------------------------- example ----------------------------------
      DESCRIPTION GLOBALE DU MAILLAGE
      NOMS DES VARIABLES D''ESPACE : 'X' 'Y' 'Z'
      NOMBRE D''ELEMENTS : 100
  ----------------------------------------------------------------------------
  + BLOC [DE]
    < < HEXAEDRES ^ PRISMES ^ TETRAEDRES ^ QUADRANGLES ^ TRIANGLES ^ SEGMENTS >
      DE LAGRANGE [ AUX POINTS DE GAUSS-LOBATTO ]
      < < P1 ^ Q1 ^ P2 ^ Q2 ^ P3 ^ Q3 ^ Q4 ... Q9 ^ Q10 >
      ^ < D'ORDRE ^ DE DEGRE > ||I
      >
      ^ [TYPE] [GEOMETRIQUE] ||A
    > : ||IE ELEMENTS
  --------------------------------- examples ---------------------------------
      BLOC DE TRIANGLES DE LAGRANGE P2 : 10 ELEMENTS         BLOC DE TYPE GEOMETRIQUE TR02 : 10 ELEMENTS
              QUADRANGLES DE LAGRANGE Q2 : 5 ELEMENTS   ou                            QU02 : 5 ELEMENTS
              TRIANGLES DE LAGRANGE P2 : 15 ELEMENTS                                  TR02 : 15 ELEMENTS
  ----------------------------------------------------------------------------
  X1 Y1 Z1 X2 Y2 Z2 X3 Y3 Z3 .... (coordinates of nodes of element 1)
  I1 I2 I3 ....                   (number of nodes of element 1)
  ...
  ---------------------------------- example ---------------------------------
      BLOC DE TRIANGLES DE LAGRANGE P1 : 3 ELEMENTS
      1.5000 0.0000 1.3858 0.5740 1.0000 0.0000
      5  6  1
      0.8660 0.5000 1.0000 0.0000 1.3858 0.5740
      2  1  6
      1.3858 0.5740 1.0607 1.0607 0.8660 0.5000
      6  7  2
  ----------------------------------------------------------------------------
  DOMAINE 'nomdomaine'
  + E[LEMENT] ||I [ / ||J ]
  + E[LEMENT] ||I1
    < F[ACE] ||I2 ^ A[RETE] ||I2 ^ P[OINT] ||I2 >
  --------------------------------- examples ---------------------------------
      DOMAINE 'OMEGA'
      E 1 E 3 E 5 E 7 E 8 E 9                                 equivalent to E 1 3 5 7 ELEMENT 8 9
      DOMAINE 'OMEGA' ELEMENTS 1 / 7                          equivalent to E 1 2 3 4 5 6 7
      DOMAINE 'GAMMA'
      ELEMENT 1 FACE 3 ELEMENT 2 FACE 2 ELEMENT 7 FACE 2      equivalent to  E 1 F 3 E 2 F 2 E 7 F 2
  ----------------------------------------------------------------------------
  FIN
  \endverbatim

  !!!! WARNING !!!!
    - do not write more than 80 characters per line ...
*/

#include "../Mesh.hpp"

namespace xlifepp {
using std::endl;


// save mesh and related domains to files
void Mesh::saveToMel(const string_t& fname, bool withDomains) const {
  string_t f = fname + ".mel";
  std::ofstream out(f.c_str());
  out.precision(fullPrec);
  melExport(out);
  out.close();
}

/*!
  Write the string str on the output stream out, limiting the length of all
  output lines to 80 characters
  Nota: the length of str is logically assumed to be < 80 characters. This is not
        checked since this is always the case when this function is called by the
        function melExport ; moreover, this would not produce any error: the string
        str would just be written as is on a new line.

  ncar is the number of characters previously written on the current line.
  If the sum of the length of str and ncar does not exceed 80 characters,
  then str is written at the end of the current line ; otherwise str is
  written on a new line, which then become the (new) current line.
  The function returns the updated number of characters written on the current line.
*/
number_t writeLigne(std::ostream& out, const std::string& str, number_t ncar) {
  number_t ns=str.size();
  if (ns+ncar>=80) {
    out << endl;
    ncar=0;
  }
  out << str;
  ncar+=ns;
  return ncar;
}

void Mesh::melExport(std::ostream& out) const {
  using std::map;
  using std::vector;
  using std::set;
  using std::ostringstream;

  map<ShapeType, string_t> melNames;
  melNames[_point]       = "PO";
  melNames[_segment]     = "SE";
  melNames[_triangle]    = "TR";
  melNames[_quadrangle]  = "QU";
  melNames[_tetrahedron] = "TE";
  melNames[_hexahedron]  = "HE";
  melNames[_prism]       = "PR";
  melNames[_pyramid]     = "PY";

  vector<string_t> tl;
  ostringstream ssti;
  ssti << "Generated by xlife++ from mesh " << name();
  string_t ti=ssti.str();

  while (ti.size()>0) {
    tl.push_back(ti.substr(0,80));
    ti.erase(0,80);
  }
  number_t ntl=tl.size();
  out << "TITRE " << ntl << endl;
  for (number_t i=0;i<ntl;i++) {
    out << tl[i] << endl;
  }

  out << "FORMAT DE LECTURE DES COORDONNEES '*'" << endl;
  out << "       DE LA NUMEROTATION GLOBALE '*'" << endl;
  out << "       SANS COMMENTAIRE" << endl;
  out << "DESCRIPTION GLOBALE DU MAILLAGE" << endl;
  string_t coor;
  switch (meshDim()) {
    case 3 : coor = " 'Z'";
    case 2 : coor = " 'Y'" + coor;
    case 1 : coor =  "'X'" + coor;
  }
  out << "NOM DES VARIABLES D''ESPACE : " << coor << endl;

  number_t ne=0;
  std::multimap<ShapeType,number_t> liElt;
  vector<ShapeType> liType;
  set<number_t> liNo;// will contain the node numbers (unique occurrence of each)

  for (vector<GeomDomain*>::const_iterator it_dom=domains().begin();it_dom!=domains().end();it_dom++) {
    if ((*it_dom)->dim()==spaceDim() && (*it_dom)->domType()==_meshDomain) {
      MeshDomain* meshDom_p=(*it_dom)->meshDomain();
      number_t nbEltsInDom = meshDom_p->numberOfElements();
      ne += nbEltsInDom;
      for (number_t i=0; i < nbEltsInDom; i++) {
        ShapeType t = meshDom_p->geomElements[i]->shapeType();
        if (find(liType.begin(),liType.end(),t)==liType.end()) {
          liType.push_back(t);
        }
        GeomElement* ge_p = meshDom_p->geomElements[i];
        liElt.insert(std::make_pair(t,ge_p->number()));
        for (vector<number_t>::const_iterator it=ge_p->meshElement()->nodeNumbers.begin();
             it!=ge_p->meshElement()->nodeNumbers.end(); it++) { liNo.insert(*it); }
      }
    }
  }
  // noNo will contain the correspondence between the node numbers and the numbers to be used
  // in the file, which should be a continuous sequence of integers starting from 1
  map<number_t,number_t> noNo;
  number_t newNo=1;
  for (set<number_t>::const_iterator it=liNo.begin(); it != liNo.end(); it++) {
    noNo[*it] = newNo++;
  }
  liNo.clear();

  // write element of mesh dimension by type
  out << "NOMBRE D''ELEMENTS : " << ne << endl;
  // element type list
  string_t com="BLOC DE TYPE GEOMETRIQUE ";
  for (size_t i=0; i<liType.size(); i++) {
    ShapeType t = liType[i];
    number_t nt = liElt.count(t);
    out << com << melNames[t];
    if (order()<=9) { out << "0";}
    out << order() << " : " << nt << " ELEMENTS" << endl;
  }

  vector<number_t> renumbering(nbOfElements()+1);// first element unused
  // coordinates and indices of all elements of dimension dim, sorted by type
  number_t elNo=0; // for the reindexing of elements (continuous sequence of integers starting from 1)

  typedef std::multimap<ShapeType,number_t>::iterator it_m;
  for (number_t i=0; i<liType.size(); i++) {
    std::pair<it_m,it_m> pit=liElt.equal_range(liType[i]);

    for (it_m jt=pit.first; jt!=pit.second; ++jt) {
      number_t elNum = jt->second;
      renumbering[elNum] = ++elNo;
      const GeomElement& geomElt = element(elNum);
      if (geomElt.hasMeshElement()) {
        for (number_t j=0; j<geomElt.numberOfNodes(); j++) {
          number_t ncar=0;
          for (number_t k=0; k<geomElt.elementDim(); k++) {
            ostringstream ss;
            ss << (*geomElt.meshElement()->nodes[j])[k] << " ";
            ncar = writeLigne(out,ss.str(),ncar);
          }
        }
        out << endl;
        number_t ncar=0;
        for (number_t j=0; j<geomElt.numberOfNodes(); j++) {
          ostringstream ss;
          ss << noNo[geomElt.meshElement()->nodeNumbers[j]] << " ";
          ncar = writeLigne(out,ss.str(),ncar);
        }
      }
      out << endl;
    }
  }

  // list of domains (multi dimensional domain are not processed !)
  for (vector<GeomDomain*>::const_iterator it_dom=domains().begin(); it_dom!=domains().end(); it_dom++) {
    out << "DOMAINE " << "\'" << (*it_dom)->name() << "\'" << endl;
    if ((*it_dom)->domType()==_meshDomain) {
      MeshDomain* meshDom_p=(*it_dom)->meshDomain();
      number_t nbEltsInDom = meshDom_p->numberOfElements();

      if ((*it_dom)->dim()==meshDim() && (*it_dom)->domType()==_meshDomain ) { // mesh.dim domain
        out << "ELEMENT ";
        number_t ncar=8, lastEl;
        for (number_t i=0; i<nbEltsInDom; i++) {
          ostringstream ss ;
          lastEl = renumbering[meshDom_p->geomElements[i]->number()];
          ss << lastEl << " ";
          ncar = writeLigne(out,ss.str(),ncar);
        }
        if (nbEltsInDom == 1) { out << "/ " << lastEl; }// This allows correct reading of the file, which fails without that.
        out << endl;
      }

      dimen_t domdim = meshDom_p->dim();
      if (domdim == meshDim()-1) { // mesh.dim-1 domain
        string_t bord;
        switch (domdim) {
          case 2 : bord="F"; break;
          case 1 : bord="A"; break;
          case 0 : bord="P"; break;
          default:
            error("domain_dim_not_handled",domdim);
        }

        number_t ncar=0;
        for (number_t i=0; i<nbEltsInDom; i++) {
          vector<GeoNumPair>::iterator it_gnp=meshDom_p->geomElements[i]->parentSides().begin();
          ostringstream ss;
          ss << "E " << it_gnp->first->number() << " " << bord << " " << it_gnp->second << " ";
          ncar = writeLigne(out,ss.str(),ncar);
        }
        out << endl;
      }
    }
  }

  // end of mel file
  out << "FIN" << endl;
}

} // end of namespace xlifepp