MElement.cpp

    fwrite(verts, sizeof(int), n + 1, fp);
  }
  else{
    fprintf(fp, "%d", n);
    for(int i = 0; i < n; i++)
      fprintf(fp, " %d", getVertexVTK(i)->getIndex() - 1);
    fprintf(fp, "\n");
  }
}

void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
{
  int type = getTypeForUNV();
  if(!type) return;

  setVolumePositive();
  int n = getNumVertices();
  int physical_property = elementary;
  int material_property = abs(physical);
  int color = 7;
  fprintf(fp, "%10d%10d%10d%10d%10d%10d\n",
          num ? num : _num, type, physical_property, material_property, color, n);
  if(type == 21 || type == 24) // linear beam or parabolic beam
    fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);

  if(physical < 0) revert();

  for(int k = 0; k < n; k++) {
    fprintf(fp, "%10d", getVertexUNV(k)->getIndex());
    if(k % 8 == 7)
      fprintf(fp, "\n");
  }
  if(n - 1 % 8 != 7)
    fprintf(fp, "\n");

  if(physical < 0) revert();
}

void MElement::writeMESH(FILE *fp, int elementTagType, int elementary,
                         int physical)
{
  setVolumePositive();
  for(int i = 0; i < getNumVertices(); i++)
    fprintf(fp, " %d", getVertex(i)->getIndex());
  fprintf(fp, " %d\n", (elementTagType == 3) ? _partition :
          (elementTagType == 2) ? physical : elementary);
}

void MElement::writeIR3(FILE *fp, int elementTagType, int num, int elementary,
                        int physical)
{
  int numVert = getNumVertices();
  bool ok = setVolumePositive();
  if(getDim() == 3 && !ok) Msg::Error("Element %d has zero volume", num);
  fprintf(fp, "%d %d %d", num, (elementTagType == 3) ? _partition :
          (elementTagType == 2) ? physical : elementary, numVert);
  for(int i = 0; i < numVert; i++)
    fprintf(fp, " %d", getVertex(i)->getIndex());
  fprintf(fp, "\n");
}

void MElement::writeBDF(FILE *fp, int format, int elementTagType, int elementary,
                        int physical)
{
  const char *str = getStringForBDF();
  if(!str) return;

  setVolumePositive();
  int n = getNumVertices();
  const char *cont[4] = {"E", "F", "G", "H"};
  int ncont = 0;

  int tag =  (elementTagType == 3) ? _partition : (elementTagType == 2) ?
    physical : elementary;

  if(format == 0){ // free field format
    fprintf(fp, "%s,%d,%d", str, _num, tag);
    for(int i = 0; i < n; i++){
      fprintf(fp, ",%d", getVertexBDF(i)->getIndex());
      if(i != n - 1 && !((i + 3) % 8)){
        fprintf(fp, ",+%s%d\n+%s%d", cont[ncont], _num, cont[ncont], _num);
        ncont++;
      }
    }
    if(n == 2) // CBAR
      fprintf(fp, ",0.,0.,0.");
    fprintf(fp, "\n");
  }
  else{ // small or large field format
    fprintf(fp, "%-8s%-8d%-8d", str, _num, tag);
    for(int i = 0; i < n; i++){
      fprintf(fp, "%-8d", getVertexBDF(i)->getIndex());
      if(i != n - 1 && !((i + 3) % 8)){
        fprintf(fp, "+%s%-6d\n+%s%-6d", cont[ncont], _num, cont[ncont], _num);
        ncont++;
      }
    }
    if(n == 2) // CBAR
      fprintf(fp, "%-8s%-8s%-8s", "0.", "0.", "0.");
    fprintf(fp, "\n");
  }
}

void MElement::writeDIFF(FILE *fp, int num, bool binary, int physical_property)
{
  const char *str = getStringForDIFF();
  if(!str) return;

  setVolumePositive();

  int n = getNumVertices();
  if(binary){
    // TODO
  }
  else{
    fprintf(fp, "%d %s %d ", num, str, physical_property);
    for(int i = 0; i < n; i++)
      fprintf(fp, " %d", getVertexDIFF(i)->getIndex());
    fprintf(fp, "\n");
  }
}

void MElement::writeINP(FILE *fp, int num)
{
  setVolumePositive();
  fprintf(fp, "%d", num);
  for(int i = 0; i < getNumVertices(); i++)
    fprintf(fp, ", %d", getVertexINP(i)->getIndex());
  fprintf(fp, "\n");
}

int MElement::getInfoMSH(const int typeMSH, const char **const name)
{
  switch(typeMSH){
  case MSH_PNT     : if(name) *name = "Point";            return 1;
  case MSH_LIN_2   : if(name) *name = "Line 2";           return 2;
  case MSH_LIN_3   : if(name) *name = "Line 3";           return 2 + 1;
  case MSH_LIN_4   : if(name) *name = "Line 4";           return 2 + 2;
  case MSH_LIN_5   : if(name) *name = "Line 5";           return 2 + 3;
  case MSH_LIN_6   : if(name) *name = "Line 6";           return 2 + 4;
  case MSH_LIN_7   : if(name) *name = "Line 7";           return 2 + 5;
  case MSH_LIN_8   : if(name) *name = "Line 8";           return 2 + 6;
  case MSH_LIN_9   : if(name) *name = "Line 9";           return 2 + 7;
  case MSH_LIN_10  : if(name) *name = "Line 10";          return 2 + 8;
  case MSH_LIN_11  : if(name) *name = "Line 11";          return 2 + 9;
  case MSH_LIN_B   : if(name) *name = "Line Border";      return 2;
  case MSH_LIN_C   : if(name) *name = "Line Child";       return 2;
  case MSH_TRI_3   : if(name) *name = "Triangle 3";       return 3;
  case MSH_TRI_6   : if(name) *name = "Triangle 6";       return 3 + 3;
  case MSH_TRI_9   : if(name) *name = "Triangle 9";       return 3 + 6;
  case MSH_TRI_10  : if(name) *name = "Triangle 10";      return 3 + 6 + 1;
  case MSH_TRI_12  : if(name) *name = "Triangle 12";      return 3 + 9;
  case MSH_TRI_15  : if(name) *name = "Triangle 15";      return 3 + 9 + 3;
  case MSH_TRI_15I : if(name) *name = "Triangle 15I";     return 3 + 12;
  case MSH_TRI_21  : if(name) *name = "Triangle 21";      return 3 + 12 + 6;
  case MSH_TRI_28  : if(name) *name = "Triangle 28";      return 3 + 15 + 10;
  case MSH_TRI_36  : if(name) *name = "Triangle 36";      return 3 + 18 + 15;
  case MSH_TRI_45  : if(name) *name = "Triangle 45";      return 3 + 21 + 21;
  case MSH_TRI_55  : if(name) *name = "Triangle 55";      return 3 + 24 + 28;
  case MSH_TRI_66  : if(name) *name = "Triangle 66";      return 3 + 27 + 36;
  case MSH_TRI_18  : if(name) *name = "Triangle 18";      return 3 + 15;
  case MSH_TRI_21I : if(name) *name = "Triangle 21I";     return 3 + 18;
  case MSH_TRI_24  : if(name) *name = "Triangle 24";      return 3 + 21;
  case MSH_TRI_27  : if(name) *name = "Triangle 27";      return 3 + 24;
  case MSH_TRI_30  : if(name) *name = "Triangle 30";      return 3 + 27;
  case MSH_TRI_B   : if(name) *name = "Triangle Border";  return 3;
  case MSH_QUA_4   : if(name) *name = "Quadrilateral 4";  return 4;
  case MSH_QUA_8   : if(name) *name = "Quadrilateral 8";  return 4 + 4;
  case MSH_QUA_9   : if(name) *name = "Quadrilateral 9";  return 9;
  case MSH_QUA_16  : if(name) *name = "Quadrilateral 16"; return 16;
  case MSH_QUA_25  : if(name) *name = "Quadrilateral 25"; return 25;
  case MSH_QUA_36  : if(name) *name = "Quadrilateral 36"; return 36;
  case MSH_QUA_49  : if(name) *name = "Quadrilateral 49"; return 49;
  case MSH_QUA_64  : if(name) *name = "Quadrilateral 64"; return 64;
  case MSH_QUA_81  : if(name) *name = "Quadrilateral 81"; return 81;
  case MSH_QUA_100 : if(name) *name = "Quadrilateral 100";return 100;
  case MSH_QUA_121 : if(name) *name = "Quadrilateral 121";return 121;
  case MSH_QUA_12  : if(name) *name = "Quadrilateral 12"; return 12;
  case MSH_QUA_16I : if(name) *name = "Quadrilateral 16I";return 16;
  case MSH_QUA_20  : if(name) *name = "Quadrilateral 20"; return 20;
  case MSH_POLYG_  : if(name) *name = "Polygon";          return 0;
  case MSH_POLYG_B : if(name) *name = "Polygon Border";   return 0;
  case MSH_TET_4   : if(name) *name = "Tetrahedron 4";    return 4;
  case MSH_TET_10  : if(name) *name = "Tetrahedron 10";   return 4 + 6;
  case MSH_TET_20  : if(name) *name = "Tetrahedron 20";   return 4 + 12 + 4;
  case MSH_TET_34  : if(name) *name = "Tetrahedron 34";   return 4 + 18 + 12 + 0;
  case MSH_TET_35  : if(name) *name = "Tetrahedron 35";   return 4 + 18 + 12 + 1;
  case MSH_TET_52  : if(name) *name = "Tetrahedron 52";   return 4 + 24 + 24 + 0;
  case MSH_TET_56  : if(name) *name = "Tetrahedron 56";   return 4 + 24 + 24 + 4;
  case MSH_TET_84  : if(name) *name = "Tetrahedron 84";   return (7*8*9)/6;
  case MSH_TET_120 : if(name) *name = "Tetrahedron 120";  return (8*9*10)/6;
  case MSH_TET_165 : if(name) *name = "Tetrahedron 165";  return (9*10*11)/6;
  case MSH_TET_220 : if(name) *name = "Tetrahedron 220";  return (10*11*12)/6;
  case MSH_TET_286 : if(name) *name = "Tetrahedron 286";  return (11*12*13)/6;
  case MSH_HEX_8   : if(name) *name = "Hexahedron 8";     return 8;
  case MSH_HEX_20  : if(name) *name = "Hexahedron 20";    return 8 + 12;
  case MSH_HEX_27  : if(name) *name = "Hexahedron 27";    return 8 + 12 + 6 + 1;
  case MSH_HEX_64  : if(name) *name = "Hexahedron 64";    return 64;
  case MSH_HEX_125 : if(name) *name = "Hexahedron 125";   return 125;
  case MSH_HEX_216 : if(name) *name = "Hexahedron 216";   return 216;
  case MSH_HEX_343 : if(name) *name = "Hexahedron 343";   return 343;
  case MSH_HEX_512 : if(name) *name = "Hexahedron 512";   return 512;
  case MSH_HEX_729 : if(name) *name = "Hexahedron 729";   return 729;
  case MSH_HEX_1000: if(name) *name = "Hexahedron 1000";  return 1000;
  case MSH_HEX_56  : if(name) *name = "Hexahedron 56";    return 56;
  case MSH_HEX_98  : if(name) *name = "Hexahedron 98";    return 98;
  case MSH_HEX_152 : if(name) *name = "Hexahedron 152";   return 152;
  case MSH_HEX_222 : if(name) *name = "Hexahedron 222";   return 222;
  case MSH_HEX_296 : if(name) *name = "Hexahedron 296";   return 296;
  case MSH_HEX_386 : if(name) *name = "Hexahedron 386";   return 386;
  case MSH_HEX_488 : if(name) *name = "Hexahedron 488";   return 488;
  case MSH_PRI_6   : if(name) *name = "Prism 6";          return 6;
  case MSH_PRI_15  : if(name) *name = "Prism 15";         return 6 + 9;
  case MSH_PRI_18  : if(name) *name = "Prism 18";         return 6 + 9 + 3;
  case MSH_PYR_5   : if(name) *name = "Pyramid 5";        return 5;
  case MSH_PYR_13  : if(name) *name = "Pyramid 13";       return 5 + 8;
  case MSH_PYR_14  : if(name) *name = "Pyramid 14";       return 5 + 8 + 1;
  case MSH_POLYH_  : if(name) *name = "Polyhedron";       return 0;
  case MSH_PNT_SUB : if(name) *name = "Point Xfem";       return 1;
  case MSH_LIN_SUB : if(name) *name = "Line Xfem";        return 2;
  case MSH_TRI_SUB : if(name) *name = "Triangle Xfem";    return 3;
  case MSH_TET_SUB : if(name) *name = "Tetrahedron Xfem"; return 4;
  default:
    Msg::Error("Unknown type of element %d", typeMSH);
    if(name) *name = "Unknown";
    return 0;
  }
}

void MElement::getVerticesIdForMSH(std::vector<int> &verts)
{
  int n = getNumVerticesForMSH();
  verts.resize(n);
  for(int i = 0; i < n; i++)
    verts[i] = getVertex(i)->getIndex();
}

MElement *MElement::copy(std::map<int, MVertex*> &vertexMap,
                         std::map<MElement*, MElement*> &newParents,
                         std::map<MElement*, MElement*> &newDomains)
{
  if(newDomains.count(this))
    return newDomains.find(this)->second;
  std::vector<MVertex*> vmv;
  int eType = getTypeForMSH();
  MElement *eParent = getParent();
  if(getNumChildren() == 0) {
    for(int i = 0; i < getNumVertices(); i++) {
      MVertex *v = getVertex(i);
      int numV = v->getNum(); //Index();
      if(vertexMap.count(numV))
        vmv.push_back(vertexMap[numV]);
      else {
        MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV);
        vmv.push_back(mv);
        vertexMap[numV] = mv;
      }
    }
  }
  else {
    for(int i = 0; i < getNumChildren(); i++) {
      for(int j = 0; j < getChild(i)->getNumVertices(); j++) {
        MVertex *v = getChild(i)->getVertex(j);
        int numV = v->getNum(); //Index();
        if(vertexMap.count(numV))
          vmv.push_back(vertexMap[numV]);
        else {
          MVertex *mv = new MVertex(v->x(), v->y(), v->z(), 0, numV);
          vmv.push_back(mv);
          vertexMap[numV] = mv;
        }
      }
    }
  }

  MElement *parent=0;
  if(eParent && !getDomain(0) && !getDomain(1)) {
    std::map<MElement*, MElement*>::iterator it = newParents.find(eParent);
    MElement *newParent;
    if(it == newParents.end()) {
      newParent = eParent->copy(vertexMap, newParents, newDomains);
      newParents[eParent] = newParent;
    }
    else
      newParent = it->second;
    parent = newParent;
  }

  MElementFactory factory;
  MElement *newEl = factory.create(eType, vmv, getNum(), _partition, ownsParent(), parent);

  for(int i = 0; i < 2; i++) {
    MElement *dom = getDomain(i);
    if(!dom) continue;
    std::map<MElement*, MElement*>::iterator it = newDomains.find(dom);
    MElement *newDom;
    if(it == newDomains.end()) {
      newDom = dom->copy(vertexMap, newParents, newDomains);
      newDomains[dom] = newDom;
    }
    else
      newDom = newDomains.find(dom)->second;
    newEl->setDomain(newDom, i);
  }
  return newEl;
}

MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
                                  int num, int part, bool owner, MElement *parent,
                                  MElement *d1, MElement *d2)
{
  switch (type) {
  case MSH_PNT:     return new MPoint(v, num, part);
  case MSH_LIN_2:   return new MLine(v, num, part);
  case MSH_LIN_3:   return new MLine3(v, num, part);
  case MSH_LIN_4:   return new MLineN(v, num, part);
  case MSH_LIN_5:   return new MLineN(v, num, part);
  case MSH_LIN_6:   return new MLineN(v, num, part);
  case MSH_LIN_7:   return new MLineN(v, num, part);
  case MSH_LIN_8:   return new MLineN(v, num, part);
  case MSH_LIN_9:   return new MLineN(v, num, part);
  case MSH_LIN_10:  return new MLineN(v, num, part);
  case MSH_LIN_11:  return new MLineN(v, num, part);
  case MSH_LIN_B:   return new MLineBorder(v, num, part, d1, d2);
  case MSH_LIN_C:   return new MLineChild(v, num, part, owner, parent);
  case MSH_TRI_3:   return new MTriangle(v, num, part);
  case MSH_TRI_6:   return new MTriangle6(v, num, part);
  case MSH_TRI_9:   return new MTriangleN(v, 3, num, part);
  case MSH_TRI_10:  return new MTriangleN(v, 3, num, part);
  case MSH_TRI_12:  return new MTriangleN(v, 4, num, part);
  case MSH_TRI_15:  return new MTriangleN(v, 4, num, part);
  case MSH_TRI_15I: return new MTriangleN(v, 5, num, part);
  case MSH_TRI_21:  return new MTriangleN(v, 5, num, part);
  case MSH_TRI_28:  return new MTriangleN(v, 6, num, part);
  case MSH_TRI_36:  return new MTriangleN(v, 7, num, part);
  case MSH_TRI_45:  return new MTriangleN(v, 8, num, part);
  case MSH_TRI_55:  return new MTriangleN(v, 9, num, part);
  case MSH_TRI_66:  return new MTriangleN(v,10, num, part);
  case MSH_TRI_B:   return new MTriangleBorder(v, num, part, d1, d2);
  case MSH_QUA_4:   return new MQuadrangle(v, num, part);
  case MSH_QUA_8:   return new MQuadrangle8(v, num, part);
  case MSH_QUA_9:   return new MQuadrangle9(v, num, part);
  case MSH_QUA_12:  return new MQuadrangleN(v, 3, num, part);
  case MSH_QUA_16:  return new MQuadrangleN(v, 3, num, part);
  case MSH_QUA_25:  return new MQuadrangleN(v, 4, num, part);
  case MSH_QUA_36:  return new MQuadrangleN(v, 5, num, part);
  case MSH_QUA_49:  return new MQuadrangleN(v, 6, num, part);
  case MSH_QUA_64:  return new MQuadrangleN(v, 7, num, part);
  case MSH_QUA_81:  return new MQuadrangleN(v, 8, num, part);
  case MSH_QUA_100: return new MQuadrangleN(v, 9, num, part);
  case MSH_QUA_121: return new MQuadrangleN(v, 10, num, part);
  case MSH_POLYG_:  return new MPolygon(v, num, part, owner, parent);
  case MSH_POLYG_B: return new MPolygonBorder(v, num, part, d1, d2);
  case MSH_TET_4:   return new MTetrahedron(v, num, part);
  case MSH_TET_10:  return new MTetrahedron10(v, num, part);
  case MSH_HEX_8:   return new MHexahedron(v, num, part);
  case MSH_HEX_20:  return new MHexahedron20(v, num, part);
  case MSH_HEX_27:  return new MHexahedron27(v, num, part);
  case MSH_PRI_6:   return new MPrism(v, num, part);
  case MSH_PRI_15:  return new MPrism15(v, num, part);
  case MSH_PRI_18:  return new MPrism18(v, num, part);
  case MSH_PYR_5:   return new MPyramid(v, num, part);
  case MSH_PYR_13:  return new MPyramid13(v, num, part);
  case MSH_PYR_14:  return new MPyramid14(v, num, part);
  case MSH_TET_20:  return new MTetrahedronN(v, 3, num, part);
  case MSH_TET_34:  return new MTetrahedronN(v, 3, num, part);
  case MSH_TET_35:  return new MTetrahedronN(v, 4, num, part);
  case MSH_TET_52:  return new MTetrahedronN(v, 5, num, part);
  case MSH_TET_56:  return new MTetrahedronN(v, 5, num, part);
  case MSH_TET_84:  return new MTetrahedronN(v, 6, num, part);
  case MSH_TET_120: return new MTetrahedronN(v, 7, num, part);
  case MSH_TET_165: return new MTetrahedronN(v, 8, num, part);
  case MSH_TET_220: return new MTetrahedronN(v, 9, num, part);
  case MSH_TET_286: return new MTetrahedronN(v, 10, num, part);
  case MSH_POLYH_:  return new MPolyhedron(v, num, part, owner, parent);
  case MSH_HEX_56:  return new MHexahedronN(v, 3, num, part);
  case MSH_HEX_64:  return new MHexahedronN(v, 3, num, part);
  case MSH_HEX_125: return new MHexahedronN(v, 4, num, part);
  case MSH_HEX_216: return new MHexahedronN(v, 5, num, part);
  case MSH_HEX_343: return new MHexahedronN(v, 6, num, part);
  case MSH_HEX_512: return new MHexahedronN(v, 7, num, part);
  case MSH_HEX_729: return new MHexahedronN(v, 8, num, part);
  case MSH_HEX_1000:return new MHexahedronN(v, 9, num, part);
  case MSH_PNT_SUB: return new MSubPoint(v, num, part, owner, parent);
  case MSH_LIN_SUB: return new MSubLine(v, num, part, owner, parent);
  case MSH_TRI_SUB: return new MSubTriangle(v, num, part, owner, parent);
  case MSH_TET_SUB: return new MSubTetrahedron(v, num, part, owner, parent);
  default:          return 0;
  }
}

MElement *MElementFactory::create(int num, int type, const std::vector<int> &tags,
                                  std::vector<MVertex*> &v,
                                  std::map<int, MElement*> &elementCache,
                                  std::vector<short> &ghosts)
{
  MElement *parent = 0;
  int part = 0;
  if(tags.size() > 2 && (type == MSH_PNT_SUB || type == MSH_LIN_SUB ||
                         type == MSH_TRI_SUB || type == MSH_TET_SUB)){
    parent = elementCache[tags[1]];
    if(tags.size() > 3 && tags[2]){ // num partitions
      part = tags[3];
      for(int i = 0; i < tags[2] - 1; i++)
        ghosts.push_back(tags[4 + i]);
    }
  }
  else if(tags.size() > 1){
    if(tags[1]){ // num partitions
      part = tags[2];
      for(int i = 0; i < tags[1] - 1; i++)
        ghosts.push_back(tags[3 + i]);
    }
  }
  create(type, v, num, part, false, parent);
}

void MElement::xyzTouvw(fullMatrix<double> *xu)
{
  double _xyz[3] = {(*xu)(0,0),(*xu)(0,1),(*xu)(0,2)}, _uvw[3];
  xyz2uvw(_xyz, _uvw);
  (*xu)(1,0) = _uvw[0];
  (*xu)(1,1) = _uvw[1];
  (*xu)(1,2) = _uvw[2];
}