geometry/ioMesh/loadMedit.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 loadMedit.cpp
  \author N. Kielbasiewicz
  \since 9 may 2016
  \date 9 may 2016

  \brief Read a file containing a mesh in Medit format

  The useful function is loadMedit, which is a member function of the class Mesh.
*/

#include "ioutils.hpp"

namespace xlifepp
{

//----------------------------------------------------------------------------------------
// Comparison function of two Medit elements according to type and domain numbers:
//  a < b if its type number is smaller, then if its domain number is smaller.
struct MELT  // Gmsh element
{
  ShapeType elty; // type of the element
  number_t ndom,   // number of the physical entity
           readOrder; // only used to make the sort function give the same result on different systems
  std::vector<number_t> nums; // node numbers defining the element
  GeomElement* ge_p; // address of the corresponding GeomElement
  bool plain;        // true if this GeomElement is created as a plain element, false if it is a side element

  MELT(): ge_p(0), plain(true) {} // plain element by default
  void side() { plain = false; }
};

bool compareMELTs(const MELT &a, const MELT &b)
{
  if (a.elty < b.elty) { return true; }
  else
  {
    if (a.elty == b.elty)
    {
      if (a.ndom < b.ndom) { return true; }
      else
      {
        if (a.ndom == b.ndom)
        {
          if (a.readOrder < b.readOrder) { return true; }
          else { return false; }
        }
        else { return false; }
      }
    }
    else { return false; }
  }
}
//----------------------------------------------------------------------------------------
// Returns the name of the domain ndom: it is the name found in the data file if any,
// or it is a generated one whose form is "Omega"ndom.
string_t genDomName(number_t ndom) { return "Omega"+tostring(ndom); }

//----------------------------------------------------------------------------------------
// Returns the name of the sub-domain nsdom: it is the name found in the data file if any,
// or it is a generated one whose form is "dOmega"ndom_nsdom.
string_t genSDomName(number_t ndom, number_t nsdom) { return "d"+genDomName(ndom)+"_"+tostring(nsdom); }

void loadMeditElements(std::istream& data, number_t spaceDim, ShapeType type, number_t dim, number_t nbVertices, std::vector<MELT>& melts, std::map<number_t, std::set<number_t> >& doms, std::map<ShapeType, ELTDEF>& elMap)
{
  number_t nbElems;
  data >> nbElems;

  number_t index=melts.size();
  melts.resize(index+nbElems);
  for (number_t i=0; i < nbElems; ++i)
  {
    MELT tmp;
    tmp.readOrder=index;
    tmp.elty=type;
    tmp.nums.resize(nbVertices);
    for (number_t e=0; e < nbVertices; ++e) { data >> tmp.nums[e]; }
    data >> tmp.ndom;
    melts[index]=tmp;
    doms[dim].insert(tmp.ndom);
    // Find or create reference element
    Interpolation* interp_p = findInterpolation(Lagrange, standard, 1, H1);
    RefElement* re_p = findRefElement(type, interp_p);
    // Store data related to this element type for future use:
    elMap[type] = ELTDEF(nbVertices, dim, re_p);
    index++;
  }
}

//! loads mesh from input file filename according to medit mesh format
void Mesh::loadMedit(const string_t& filename, number_t nodesDim)
{
  trace_p->push("Mesh::loadMedit");
  // error("not_yet_implemented", "Mesh::loadMedit(String filename)");

  std::ifstream data(filename.c_str());
  if ( !data ) { error("file_failopen","Mesh::loadMedit",filename); }

  comment_ = string_t("Medit mesh format read from file ") + basenameWithExtension(filename);
  if ( theVerboseLevel > 1) { info("loadFile_info", "Medit", filename); }

  string_t strVal;
  number_t nVal = 0;

  data >> strVal;
  bool lineError = false;

  if (strVal != "MeshVersionFormatted") { lineError=true; }
  if (!lineError) { data >> nVal; }
  if (nVal > 2) { lineError=true; }
  if (lineError) { error("mel_first_line"); }
  data >> strVal;
  if (strVal != "Dimension") { error("mel_second_line"); }
  number_t spaceDim=0;
  if (!lineError) { data >> spaceDim; }
  if (nVal > 3) { lineError=true; }
  if (lineError) { error("mel_bad_dim"); }

  std::map<number_t, std::map<number_t, std::vector<GeomElement*> > > physicalDomains;
  std::vector<MELT> melts;
  std::map<number_t, std::set<number_t> > doms;
  // typedef std::map<ShapeType, ELTDEF> ELTDEFMAP;
  std::map<ShapeType, ELTDEF> elMap;         // stores data for each element type found in the Gmsh data file
  isMadeOfSimplices_=true;
  while (data >> strVal)
  {
    if (strVal == "Vertices")
    {
      number_t nbPoints;
      data >> nbPoints;

      // Definition of nodes
      nodes.resize(nbPoints,Point(std::vector<real_t>(spaceDim)));
      vertices_.resize(nbPoints);
      for (number_t i=0; i < nbPoints; ++i)
      {
        for(number_t j=0; j < spaceDim; ++j) { data >> nodes[i][j]; }
        data >> nVal; // reference of node NOT USED
        vertices_[i]=i+1;
      }
    }
    else if (strVal == "Triangles") { loadMeditElements(data, spaceDim, _triangle, 2, 3, melts, doms, elMap); }
    else if (strVal == "Quadrilaterals") { loadMeditElements(data, spaceDim, _quadrangle, 2, 4, melts, doms, elMap); isMadeOfSimplices_=false; }
    else if (strVal == "Tetrahedra") { loadMeditElements(data, spaceDim, _tetrahedron, 3, 4, melts, doms, elMap); }
    else if (strVal == "Hexaedra") { loadMeditElements(data, spaceDim, _hexahedron, 3, 8, melts, doms, elMap); isMadeOfSimplices_=false; }
    else if (strVal == "Edges") { loadMeditElements(data, spaceDim, _segment, 1, 2, melts, doms, elMap); }
    else if (strVal == "End") { break; }
    else { error("mel_block_not_handled", strVal); }
  }

  // Sort elements by type, then by domain number. The following actions rely on this.
  sort(melts.begin(), melts.end(), compareMELTs);

  number_t nb_elts=melts.size();
  number_t nbDoms=0;
  number_t maxDim=0;
  number_t elmDim=0;
  ShapeType elmType, prevelmType=_noShape; // initialized to a non existing element type
  number_t ndomx = 0;
  // std::map<GeomElement*, number_t> elements2melt;// std::vector<number_t> elements2melt(nb_elts,0);
  RefElement* re_p=0;
  for (std::map<number_t, std::set<number_t> >::const_iterator it=doms.begin(); it != doms.end(); ++it)
  {
    nbDoms+=it->second.size();
    if (it->first > maxDim) { maxDim = it->first; }
  }

  number_t no_elt=0;             // the number associated to a plain element, updated when a new GeomElement is created
  number_t no_sideElt=nb_elts+1; // the number associated to a side  element, updated when a new GeomElement is created
  // ==> the numbers of the plain elements are in [1, nb_elts]
  //     the numbers of the side elements start from nb_elts+1
  // This is a temporary situation, true only during the construction phase of the elements. At the end,
  // the nb_elts elements are renumbered so that their numbers are consecutive and lying in [1, nb_elts],
  // following the order of creation of the elements.
  std::vector<GeomElement*> allElts(nb_elts); // Used to temporary store adresses of all elements present in the file.
  // If they are created as side elements during one iteration, they are not created as plain elements during
  // next iteration, in the following "decreasing loop on elements dimension".
  elements_.clear();
  domains_.reserve(nbDoms);
  std::vector<GeomElement*>::const_iterator itb_Elts = allElts.begin(), // corresponds to no_elt(0)
      itb_EltsCurDim, ite_EltsCurDim = itb_Elts;
  // The following vector records the numbers of the subdomains created in order to
  // not also create them as plain domains:
  std::vector<number_t> sdNum;

  for (int_t curEltDim=maxDim;  curEltDim>=0; curEltDim--) // decreasing loop on elements dimension
  {
    //  b. Create elements of dimension curEltDim
    SIDELTMAP parentEl;
    itb_EltsCurDim = ite_EltsCurDim;
    number_t curEltDimM1 = curEltDim - 1;
    // startDom holds the rank (in allElts) of the first element created (of dimension curEltDim)
    // in a domain and the associated domain number.
    // -> a domain may be made of elements of different types
    std::vector<pairNN > startDom;
    prevelmType = _noShape;           // initialize to a non existing element type
    number_t prevndom = ndomx; // initialize to a non existing domain number

    for (number_t el_n=0; el_n<nb_elts; el_n++)
    {
      elmType = melts[el_n].elty;
      // Detect change of element type in the list:
      if (elmType != prevelmType)
      {
        elmDim = elMap.at(elmType).elmDim;
        re_p   = elMap.at(elmType).refElt_p;
        prevelmType = elmType;
      }
      // Create element only if the dimension is the current one.
      if (elmDim == static_cast<number_t>(curEltDim))
      {
        // Detect change of domain number and record it for future use
        if (melts[el_n].ndom != prevndom)
        {
          prevndom = melts[el_n].ndom;
          // no_elt is the first element of domain number prevndom
          startDom.push_back(std::make_pair(no_elt,prevndom));
        }

        if (melts[el_n].plain)  // Create new plain element.
        {
          allElts[no_elt] = melts[el_n].ge_p = new GeomElement(this, re_p, spaceDim, no_elt + 1);
          // Store the node numbers inside the GeomElement.
          // Apply a permutation to Gmsh numbering in order to get the corresponding one in XLiFE++.
          MeshElement* mshelt = allElts[no_elt]->meshElement();
          number_t *gnum_el = &melts[el_n].nums[0];
          for (number_t i=0; i<elMap.at(elmType).nbPts; i++) { mshelt->nodeNumbers[i] = gnum_el[i];}
          mshelt->setNodes(nodes); // update node pointers

          // Update vertex numbers: the vertexNumbers are the first of nodeNumbers
          // Mark them in the global list markvert, formerly initialized to 0 (vertex numbers are > 0)
          number_t nbVert = re_p->geomRefElement()->nbVertices();
          for (number_t i=0; i<nbVert; i++)
          {
            mshelt->vertexNumbers[i] = mshelt->nodeNumbers[i];
          }
        }
        else
        {
          // melts[el_n] has already been created as a side element in the previous loop iteration:
          // get the address of the corresponding (side) GeomElement
          allElts[no_elt] = melts[el_n].ge_p;
        }
        // Store sides of this new GeomElement in a map : side ----> (parent element, side number)
        storeElSides(allElts[no_elt], no_elt, parentEl);

        // elements2melt.insert(std::make_pair(allElts[no_elt], el_n));
        no_elt++;
      }
    }
    // Record position of fictive first element of fictive next domain:
    startDom.push_back(std::make_pair(no_elt,0));

    ite_EltsCurDim = std::vector<GeomElement*>::const_iterator(allElts.begin()+no_elt);

    // From now on, itb_EltsCurDim and ite_EltsCurDim define the range in allElts
    // of the elements of current dimension just created.

    //  c. Create domains made of elements of dimension curEltDim, except if they have
    //     already been created as subdomains.
    //     Store the elements of the created domains in std::vector<GeomElement*> elements_ of mesh.
    //     By construction, each domain is made of elements whose numbers are consecutive.
    std::vector<pairNN >::iterator it=startDom.begin();
    nVal = it->second;
    number_t first_no_elt(it->first);
    for (it++; it < startDom.end(); it++) {
      if (find(sdNum.begin(),sdNum.end(),nVal) == sdNum.end()) {
        // nVal has not been created as a subdomain ; create it as a plain domain.
        string_t domName = genDomName(nVal);
        std::vector<GeomElement*>::const_iterator begin_elt(allElts.begin()+first_no_elt),
                                               end_elt(allElts.begin()+it->first);
        number_t nelts = end_elt-begin_elt;
        std::ostringstream ss;   ss << "num. ref. " << nVal << ", made of " << nelts << " elements";
        MeshDomain* meshdom_p = (new GeomDomain(*this, domName, curEltDim, ss.str()))->meshDomain();
        meshdom_p->geomElements.assign(begin_elt, end_elt);
        domains_.push_back(meshdom_p);
        // Store plain elements addresses related to this domain in elements_ vector.
        for (std::vector<GeomElement*>::const_iterator it=begin_elt; it < end_elt; it++) {
          elements_.push_back(*it);
        }
        if (theVerboseLevel > 2) { info("gmsh_domain_info",domName,nVal,nelts); }
        // nbOfEltsInDim[curEltDim] += end_elt-begin_elt;// unused
      }
      first_no_elt = it->first;
      nVal = it->second;
    }

    //  d. Create domains made of elements of dimension curEltDimM1 if they are sub-domains of domains
    //     just created, i.e. if they are boundary or interface domains.
    //     Among the elements stored in melts, only those of dimension curEltDim have been taken into
    //     account. The remaining ones (of dimension curEltDimM1) may belong to sub-domains of the
    //     domains just created. If this occur, they should be part of the boundary of some of the
    //     elements of dimension curEltDim created before.
    //

    //  d.1. First, for each of the remaining elements R of dimension curEltDimM1, identify which element
    //       E of dimension curEltDim it belongs to (if any), and which boundary part F of E it corresponds to.
    //       Each of these elements R is thus defined as a couple (E,F) where E is the rank of an element
    //       in allElts, and F is a local number in XLiFE++ numbering convention. It is an edge in 2D,
    //       a face in 3D. Theses couples (E,F) are stored in a vector used later to create the domains.
    //       In the case of an interface domain, there exist 2 elements E for each R. The one considered
    //       is the first found, and due to the ordering of the elements, the elements E corresponding to
    //       an interface all belong to the same domain.

    number_t rkxel=0;
    std::vector<std::vector<GeomElement*> > vv_ge;
    std::vector<GeomElement*> v_ge;
    std::vector<std::vector<number_t> > vv_id, vv_ln;
    std::vector<number_t> v_id, v_ndom, v_ln, v_rkxel;
    // vv_ge and vv_ln contain, for each sub-domain, the list of elements E and the corresponding
    // side number F. They are used during the construction of subdomains (d.2).
    // v_ndom contains these sub-domain numbers. For each of them, v_rkxel contains the rank of
    // maximum value in allElts among all elements found, when searching which element a sub-domain
    // element belongs to. This information is used in conjunction with startDom to recover the
    // domain number corresponding to each sub-domain.
    // vv_id contains, for each sub-domain, the indices in melts corresponding to the elements that
    // are found to be side of the elements E stored in vv_ge, in the same order. These indices are
    // used to manage the creation of the plain and side elements.

    // vv_ndom2 that contains for each subdomain, the list of secondary ref num of elements. It is
    // used during the construction of cracked subdomains (d.2).
    prevelmType = _noShape;  // initialize to a non existing element type
    prevndom = ndomx; // initialize to a non existing domain number
    v_ndom.push_back(prevndom); // Initialization to make vv_ge, vv_ln, v_rkxel and v_ndom to have
                                // the same length at the end
    for (number_t el_n=0; el_n<nb_elts; el_n++) {
      elmType = melts[el_n].elty;
      // Detect change of element type in the list:
      if (elmType != prevelmType) {
        elmDim = elMap.at(elmType).elmDim;
        prevelmType = elmType;
      }
      // Select elements whose dimension is equal to curEltDimM1 and find which element of dimension
      // curEltDim it belongs to, among those created and stored in allElts at step b. above.
      if (elmDim == curEltDimM1) {
        // Detect change of domain number and record it for future use
        if (melts[el_n].ndom != prevndom) {
          prevndom = melts[el_n].ndom;
          v_ndom.push_back(prevndom);
          v_rkxel.push_back(rkxel); rkxel = 0;
          // push then clear previous vectors which are empty at the first time
          vv_id.push_back(v_id); v_id.clear();
          vv_ge.push_back(v_ge); v_ge.clear();
          vv_ln.push_back(v_ln); v_ln.clear();
          // vv_geInd.push_back(v_geInd); v_geInd.clear();
        }
        ITPARENTS itpar; // pair of iterators on the map parentEl
        if (foundParent(melts[el_n].nums, parentEl, itpar)) {
          // Note : For now, we use the first parent found, given by itpar.first.
          // If itpar.second equals itpar.first + 1, then there is only one parent. Otherwise,
          // there are at least 2 parents, but more than 2 parents is an erroneous situation !
          // This may happen in the case of an interface domain: this element is the side of 2
          // parent elements and in this case, the second parent is given by itpar.first + 1.
          number_t rkel = itpar.first->second.first;
          if (rkel > rkxel) {rkxel = rkel;}
          v_id.push_back(el_n);
          v_ge.push_back(allElts[rkel]);              // parent element
          v_ln.push_back(itpar.first->second.second); // local side number in parent element
          // v_geInd.push_back(rkel);
          // std::cout << rkel << " side " << itpar.first->second.second << ", " <<  nbPar(itpar) << " parents" << std::endl;
        }
        // else this element is not part of any sub-domain of previously created domains.
      }
    }
    // Store last information built:
    v_rkxel.push_back(rkxel);
    vv_id.push_back(v_id); v_id.clear();
    vv_ge.push_back(v_ge); v_ge.clear();
    vv_ln.push_back(v_ln); v_ln.clear();
    // vv_geInd.push_back(v_geInd); v_geInd.clear();

    //  d.2. Second, create sub-domains (boundary or interface domains) made of elements
    //       whose dimension is equal to curEltDimM1.
    for (number_t i=1; i<v_ndom.size(); ++i) {
      if (vv_ge[i].size() > 0) {  // this one is a sub-domain of one of the domains previously created
        number_t numDom=0;
        for (std::vector<pairNN >::iterator itsd=startDom.begin(); itsd<startDom.end(); itsd++) {
          if (v_rkxel[i] >= itsd->first) { numDom = itsd->second; }
        }
        string_t domName = genSDomName(numDom,v_ndom[i]);
        std::ostringstream ss;   ss << "num. ref. " << v_ndom[i] << ", subdomain of domain whose num. ref. is "
                               << numDom << ", made of " << vv_ge[i].size() << " elements";
        MeshDomain* meshdom_p = (new GeomDomain(*this, domName, curEltDimM1, ss.str()))->meshDomain();
        std::vector<GeomElement*>& v_ge = meshdom_p->geomElements;
        std::vector<GeomElement*>::const_iterator it_ge;
        std::vector<number_t>::const_iterator it_ln, it_id;

        for (it_id=vv_id[i].begin(), it_ln=vv_ln[i].begin(), it_ge=vv_ge[i].begin();
             it_ge < vv_ge[i].end();   it_id++, it_ge++, it_ln++) {
          v_ge.push_back( melts[*it_id].ge_p = new GeomElement(*it_ge,*it_ln, no_sideElt++) );
          melts[*it_id].side();
        }

        domains_.push_back(meshdom_p);
        if (theVerboseLevel > 2) { info("gmsh_subdomain_info","  ",domName,v_ndom[i],numDom,vv_ge[i].size()); }
        // record the creation of this subdomain:
        sdNum.push_back(v_ndom[i]);
      } // end of if vv_ge[i].size() > 0
    } // for number_t i=1;...
  } // for (number_t curEltDim=elementDimx...

//  lastIndex_ = no_sideElt - 1;
  // Assign numbers to the elements according to creation order, as is done in other constructors.
  number_t numEl = 0;
  for (std::vector<GeomElement*>::const_iterator it=allElts.begin(); it < allElts.end(); it++) {
    (*it)->number() = ++numEl;
  }
  lastIndex_ = numEl;

  order_=1;
  if (elMap.begin()->second.refElt_p != 0) { order_ = elMap.begin()->second.refElt_p->order(); }
  for (std::map<ShapeType, ELTDEF>::iterator it=elMap.begin(); it != elMap.end(); it++) {
    if (it->second.refElt_p != 0) {
      if (it->second.refElt_p->order() != order_) { // all elements in the mesh should have the same approximation order
        error("gmsh_multiple_orders",filename);
      }
    }
  }

  //compute measures and orientation of mesh elements
  buildGeomData();
  setShapeTypes();

  //sides and sidesOfSides lists are not built (see buildSides and buildSideOfSides)

  firstOrderMesh_p = this;

  trace_p->pop();
}

} // end of namespace xlifepp