MElement.cpp

// Gmsh - Copyright (C) 1997-2012 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.

#include <stdlib.h>
#include <math.h>
#include "GmshConfig.h"
#include "GmshMessage.h"
#include "MElement.h"
#include "MPoint.h"
#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
#include "MElementCut.h"
#include "MSubElement.h"
#include "GEntity.h"
#include "StringUtils.h"
#include "Numeric.h"
#include "Context.h"

#define SQU(a)      ((a)*(a))

int MElement::_globalNum = 0;
double MElement::_isInsideTolerance = 1.e-6;

MElement::MElement(int num, int part) : _visible(1)
{
#pragma omp critical
  {
    if(num){
      _num = num;
      _globalNum = std::max(_globalNum, _num);
    }
    else{
      _num = ++_globalNum;
    }
    _partition = (short)part;
  }
}

void MElement::_getEdgeRep(MVertex *v0, MVertex *v1,
                           double *x, double *y, double *z, SVector3 *n,
                           int faceIndex)
{
  x[0] = v0->x(); y[0] = v0->y(); z[0] = v0->z();
  x[1] = v1->x(); y[1] = v1->y(); z[1] = v1->z();
  if(faceIndex >= 0){
    n[0] = n[1] = getFace(faceIndex).normal();
  }
  else{
    MEdge e(v0, v1);
    n[0] = n[1] = e.normal();
  }
}

void MElement::_getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2,
                           double *x, double *y, double *z, SVector3 *n)
{
  x[0] = v0->x(); x[1] = v1->x(); x[2] = v2->x();
  y[0] = v0->y(); y[1] = v1->y(); y[2] = v2->y();
  z[0] = v0->z(); z[1] = v1->z(); z[2] = v2->z();
  SVector3 t1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
  SVector3 t2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
  SVector3 normal = crossprod(t1, t2);
  normal.normalize();
  for(int i = 0; i < 3; i++) n[i] = normal;
}

char MElement::getVisibility() const
{
  if(CTX::instance()->hideUnselected && _visible < 2) return false;
  return _visible;
}

double MElement::minEdge()
{
  double m = 1.e25;
  for(int i = 0; i < getNumEdges(); i++){
    MEdge e = getEdge(i);
    m = std::min(m, e.getVertex(0)->distance(e.getVertex(1)));
  }
  return m;
}

double MElement::maxEdge()
{
  double m = 0.;
  for(int i = 0; i < getNumEdges(); i++){
    MEdge e = getEdge(i);
    m = std::max(m, e.getVertex(0)->distance(e.getVertex(1)));
  }
  return m;
}

double MElement::rhoShapeMeasure()
{
  double min = minEdge();
  double max = maxEdge();
  if(max)
    return min / max;
  else
    return 0.;
}

void MElement::scaledJacRange(double &jmin, double &jmax)
{
  jmin = jmax = 1.0;
#if defined(HAVE_MESH)
  extern double mesh_functional_distorsion(MElement*,double,double);
  if (getPolynomialOrder() == 1) return;
  const bezierBasis *jac = getJacobianFuncSpace()->bezier;
  fullVector<double>Ji(jac->points.size1());
  for (int i=0;i<jac->points.size1();i++){
    double u = jac->points(i,0);
    double v = jac->points(i,1);
    if (getType() == TYPE_QUA){
      u = -1 + 2*u;
      v = -1 + 2*v;
    }
    Ji(i) = mesh_functional_distorsion(this,u,v);
  }
  fullVector<double> Bi( jac->matrixLag2Bez.size1() );
  jac->matrixLag2Bez.mult(Ji,Bi);
  jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
  jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
#endif
}

void MElement::getNode(int num, double &u, double &v, double &w)
{
  // only for MElements that don't have a lookup table for this
  // (currently only 1st order elements have)
  double uvw[3];
  MVertex* ver = getVertex(num);
  double xyz[3] = {ver->x(), ver->y(), ver->z()};
  xyz2uvw(xyz, uvw);
  u = uvw[0];
  v = uvw[1];
  w = uvw[2];
}

void MElement::getShapeFunctions(double u, double v, double w, double s[], int o)
{
  const polynomialBasis* fs = getFunctionSpace(o);
  if(fs) fs->f(u, v, w, s);
  else Msg::Error("Function space not implemented for this type of element");
}

void MElement::getGradShapeFunctions(double u, double v, double w, double s[][3],
                                     int o)
{
  const polynomialBasis* fs = getFunctionSpace(o);
  if(fs) fs->df(u, v, w, s);
  else Msg::Error("Function space not implemented for this type of element");
}

void MElement::getHessShapeFunctions(double u, double v, double w, double s[][3][3],
                                     int o)
{
  const polynomialBasis* fs = getFunctionSpace(o);
  if(fs) fs->ddf(u, v, w, s);
  else Msg::Error("Function space not implemented for this type of element");
}

void MElement::getThirdDerivativeShapeFunctions(double u, double v, double w, double s[][3][3][3],
                                     int o){
  const polynomialBasis* fs = getFunctionSpace(o);
  if(fs) fs->dddf(u, v, w, s);
  else Msg::Error("Function space not implemented for this type of element");
};

SPoint3 MElement::barycenter_infty ()
{
  double xmin =  getVertex(0)->x();
  double xmax = xmin;
  double ymin =  getVertex(0)->y();
  double ymax = ymin;
  double zmin =  getVertex(0)->z();
  double zmax = zmin;
  int n = getNumVertices();
  for(int i = 0; i < n; i++) {
    MVertex *v = getVertex(i);
    xmin = std::min(xmin,v->x());
    xmax = std::max(xmax,v->x());
    ymin = std::min(ymin,v->y());
    ymax = std::max(ymax,v->y());
    zmin = std::min(zmin,v->z());
    zmax = std::max(zmax,v->z());
  }
  return SPoint3(0.5*(xmin+xmax),0.5*(ymin+ymax),0.5*(zmin+zmax));
}

SPoint3 MElement::barycenter()
{
  SPoint3 p(0., 0., 0.);
  int n = getNumVertices();
  for(int i = 0; i < n; i++) {
    MVertex *v = getVertex(i);
    p[0] += v->x();
    p[1] += v->y();
    p[2] += v->z();
  }
  p[0] /= (double)n;
  p[1] /= (double)n;
  p[2] /= (double)n;
  return p;
}

SPoint3 MElement::barycenterUVW()
{
  SPoint3 p(0., 0., 0.);
  int n = getNumVertices();
  for(int i = 0; i < n; i++) {
    double x, y, z;
    getNode(i, x, y, z);
    p[0] += x;
    p[1] += y;
    p[2] += z;
  }
  p[0] /= (double)n;
  p[1] /= (double)n;
  p[2] /= (double)n;
  return p;
}

double MElement::getVolume()
{
  int npts;
  IntPt *pts;
  getIntegrationPoints(getDim() * (getPolynomialOrder() - 1), &npts, &pts);
  double vol = 0.;
  for (int i = 0; i < npts; i++){
    vol += getJacobianDeterminant(pts[i].pt[0], pts[i].pt[1], pts[i].pt[2])
      * pts[i].weight;
  }
  return vol;
}

int MElement::getVolumeSign()
{
  double v = getVolume();
  if(v < 0.) return -1;
  else if(v > 0.) return 1;
  else return 0;
}

bool MElement::setVolumePositive()
{
  if(getDim() < 3) return true;
  int s = getVolumeSign();
  if(s < 0) revert();
  if(!s) return false;
  return true;
}

std::string MElement::getInfoString()
{
  char tmp[256];
  sprintf(tmp, "Element %d", getNum());
  return std::string(tmp);
}

static double _computeDeterminantAndRegularize(MElement *ele, double jac[3][3])
{
  double dJ = 0;

  switch (ele->getDim()) {

  case 0:
    {
      dJ = 1.0;
      jac[0][0] = jac[1][1] = jac[2][2] = 1.0;
      jac[0][1] = jac[1][0] = jac[2][0] = 0.0;
      jac[0][2] = jac[1][2] = jac[2][1] = 0.0;
      break;
    }
  case 1:
    {
      dJ = sqrt(SQU(jac[0][0]) + SQU(jac[0][1]) + SQU(jac[0][2]));

      // regularize matrix
      double a[3], b[3], c[3];
      a[0] = jac[0][0];
      a[1] = jac[0][1];
      a[2] = jac[0][2];
      if((fabs(a[0]) >= fabs(a[1]) && fabs(a[0]) >= fabs(a[2])) ||
         (fabs(a[1]) >= fabs(a[0]) && fabs(a[1]) >= fabs(a[2]))) {
        b[0] = a[1]; b[1] = -a[0]; b[2] = 0.;
      }
      else {
        b[0] = 0.; b[1] = a[2]; b[2] = -a[1];
      }
      norme(b);
      prodve(a, b, c);
      norme(c);
      jac[1][0] = b[0]; jac[1][1] = b[1]; jac[1][2] = b[2];
      jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
      break;
    }
  case 2:
    {
      dJ = sqrt(SQU(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0]) +
                SQU(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2]) +
                SQU(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1]));

      // regularize matrix
      double a[3], b[3], c[3];
      a[0] = jac[0][0];
      a[1] = jac[0][1];
      a[2] = jac[0][2];
      b[0] = jac[1][0];
      b[1] = jac[1][1];
      b[2] = jac[1][2];
      prodve(a, b, c);
      norme(c);
      jac[2][0] = c[0]; jac[2][1] = c[1]; jac[2][2] = c[2];
      break;
    }
  case 3:
    {
      dJ = (jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
	    jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
	    jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
      break;
    }
  }
  return dJ;
}

double MElement::getJacobian(double u, double v, double w, double jac[3][3])
{
  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
  jac[2][0] = jac[2][1] = jac[2][2] = 0.;

  double gsf[1256][3];
  getGradShapeFunctions(u, v, w, gsf);
  for (int i = 0; i < getNumShapeFunctions(); i++) {
    const MVertex *ver = getShapeFunctionNode(i);
    double* gg = gsf[i];
    for (int j = 0; j < getDim(); j++) {
      jac[j][0] += ver->x() * gg[j];
      jac[j][1] += ver->y() * gg[j];
      jac[j][2] += ver->z() * gg[j];
    }
    //    printf("GSF (%d,%g %g) = %g %g \n",i,u,v,gg[0],gg[1]);
  }
  return _computeDeterminantAndRegularize(this, jac);
}

double MElement::getJacobian(const fullMatrix<double> &gsf, double jac[3][3])
{
  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
  jac[2][0] = jac[2][1] = jac[2][2] = 0.;

  for (int i = 0; i < getNumShapeFunctions(); i++) {
    const MVertex *v = getShapeFunctionNode(i);
    for (int j = 0; j < gsf.size2(); j++) {
      jac[j][0] += v->x() * gsf(i, j);
      jac[j][1] += v->y() * gsf(i, j);
      jac[j][2] += v->z() * gsf(i, j);
    }
  }
  return _computeDeterminantAndRegularize(this, jac);
}

double MElement::getJacobian(const std::vector<SVector3> &gsf, double jac[3][3])
{
  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
  jac[2][0] = jac[2][1] = jac[2][2] = 0.;

  for (int i = 0; i < getNumShapeFunctions(); i++) {
    const MVertex *v = getShapeFunctionNode(i);
    for (int j = 0; j < 3; j++) {
      double mult = gsf[i][j];
      jac[j][0] += v->x() * mult;
      jac[j][1] += v->y() * mult;
      jac[j][2] += v->z() * mult;
    }
  }
  return _computeDeterminantAndRegularize(this, jac);
}

double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][3])
{
  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
  jac[2][0] = jac[2][1] = jac[2][2] = 0.;

  double gsf[1256][3];
  getGradShapeFunctions(u, v, w, gsf, 1);
  for(int i = 0; i < getNumPrimaryShapeFunctions(); i++) {
    const MVertex *v = getShapeFunctionNode(i);
    double* gg = gsf[i];
    for (int j = 0; j < 3; j++) {
      jac[j][0] += v->x() * gg[j];
      jac[j][1] += v->y() * gg[j];
      jac[j][2] += v->z() * gg[j];
    }
  }

  return _computeDeterminantAndRegularize(this, jac);
}

void MElement::pnt(double u, double v, double w, SPoint3 &p)
{
  double x = 0., y = 0., z = 0.;
  double sf[1256];
  getShapeFunctions(u, v, w, sf);
  for (int j = 0; j < getNumShapeFunctions(); j++) {
    const MVertex *v = getShapeFunctionNode(j);
    x += sf[j] * v->x();
    y += sf[j] * v->y();
    z += sf[j] * v->z();
  }
  p = SPoint3(x, y, z);
}

void MElement::pnt(const std::vector<double> &sf, SPoint3 &p)
{
  double x = 0., y = 0., z = 0.;
  for (int j = 0; j < getNumShapeFunctions(); j++) {
    const MVertex *v = getShapeFunctionNode(j);
    x += sf[j] * v->x();
    y += sf[j] * v->y();
    z += sf[j] * v->z();
  }
  p = SPoint3(x, y, z);
}

void MElement::primaryPnt(double u, double v, double w, SPoint3 &p)
{
  double x = 0., y = 0., z = 0.;
  double sf[1256];
  getShapeFunctions(u, v, w, sf, 1);
  for (int j = 0; j < getNumPrimaryShapeFunctions(); j++) {
    const MVertex *v = getShapeFunctionNode(j);
    x += sf[j] * v->x();
    y += sf[j] * v->y();
    z += sf[j] * v->z();
  }
  p = SPoint3(x,y,z);
}

void MElement::xyz2uvw(double xyz[3], double uvw[3])
{
  // general Newton routine for the nonlinear case (more efficient
  // routines are implemented for simplices, where the basis functions
  // are linear)
  uvw[0] = uvw[1] = uvw[2] = 0.;
  int iter = 1, maxiter = 20;
  double error = 1., tol = 1.e-6;

  while (error > tol && iter < maxiter){
    double jac[3][3];
    if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
    double xn = 0., yn = 0., zn = 0.;
    double sf[1256];
    getShapeFunctions(uvw[0], uvw[1], uvw[2], sf);
    for (int i = 0; i < getNumShapeFunctions(); i++) {
      MVertex *v = getShapeFunctionNode(i);
      xn += v->x() * sf[i];
      yn += v->y() * sf[i];
      zn += v->z() * sf[i];
    }
    double inv[3][3];
    inv3x3(jac, inv);
    double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
      inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
    double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) +
      inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn);
    double wn = uvw[2] + inv[0][2] * (xyz[0] - xn) +
      inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn);
    error = sqrt(SQU(un - uvw[0]) + SQU(vn - uvw[1]) + SQU(wn - uvw[2]));
    uvw[0] = un;
    uvw[1] = vn;
    uvw[2] = wn;
    iter++ ;
  }
}

void MElement::movePointFromParentSpaceToElementSpace(double &u, double &v, double &w)
{
  if(!getParent()) return;
  SPoint3 p;
  getParent()->pnt(u, v, w, p);
  double xyz[3] = {p.x(), p.y(), p.z()};
  double uvwE[3];
  xyz2uvw(xyz, uvwE);
  u = uvwE[0]; v = uvwE[1]; w = uvwE[2];
}

void MElement::movePointFromElementSpaceToParentSpace(double &u, double &v, double &w)
{
  if(!getParent()) return;
  SPoint3 p;
  pnt(u, v, w, p);
  double xyz[3] = {p.x(), p.y(), p.z()};
  double uvwP[3];
  getParent()->xyz2uvw(xyz, uvwP);
  u = uvwP[0]; v = uvwP[1]; w = uvwP[2];
}

double MElement::interpolate(double val[], double u, double v, double w, int stride,
                             int order)
{
  double sum = 0;
  int j = 0;
  double sf[1256];
  getShapeFunctions(u, v, w, sf, order);
  for(int i = 0; i < getNumShapeFunctions(); i++){
    sum += val[j] * sf[i];
    j += stride;
  }
  return sum;
}

void MElement::interpolateGrad(double val[], double u, double v, double w, double f[],
                               int stride, double invjac[3][3], int order)
{
  double dfdu[3] = {0., 0., 0.};
  int j = 0;
  double gsf[1256][3];
  getGradShapeFunctions(u, v, w, gsf, order);
  for(int i = 0; i < getNumShapeFunctions(); i++){
    dfdu[0] += val[j] * gsf[i][0];
    dfdu[1] += val[j] * gsf[i][1];
    dfdu[2] += val[j] * gsf[i][2];
    j += stride;
  }
  if(invjac){
    matvec(invjac, dfdu, f);
  }
  else{
    double jac[3][3], inv[3][3];
    getJacobian(u, v, w, jac);
    inv3x3(jac, inv);
    matvec(inv, dfdu, f);
  }
}

void MElement::interpolateCurl(double val[], double u, double v, double w, double f[],
                               int stride, int order)
{
  double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
  getJacobian(u, v, w, jac);
  inv3x3(jac, inv);
  interpolateGrad(&val[0], u, v, w, fx, stride, inv, order);
  interpolateGrad(&val[1], u, v, w, fy, stride, inv, order);
  interpolateGrad(&val[2], u, v, w, fz, stride, inv, order);
  f[0] = fz[1] - fy[2];
  f[1] = -(fz[0] - fx[2]);
  f[2] = fy[0] - fx[1];
}

double MElement::interpolateDiv(double val[], double u, double v, double w,
                                int stride, int order)
{
  double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
  getJacobian(u, v, w, jac);
  inv3x3(jac, inv);
  interpolateGrad(&val[0], u, v, w, fx, stride, inv, order);
  interpolateGrad(&val[1], u, v, w, fy, stride, inv, order);
  interpolateGrad(&val[2], u, v, w, fz, stride, inv, order);
  return fx[0] + fy[1] + fz[2];
}

double MElement::integrate(double val[], int pOrder, int stride, int order)
{
  int npts; IntPt *gp;
  getIntegrationPoints(pOrder, &npts, &gp);
  double sum = 0;
  for (int i = 0; i < npts; i++){
    double u = gp[i].pt[0];
    double v = gp[i].pt[1];
    double w = gp[i].pt[2];
    double weight = gp[i].weight;
    double detuvw = getJacobianDeterminant(u, v, w);
    sum += interpolate(val, u, v, w, stride, order)*weight*detuvw;
  }
  return sum;
}

static int getTriangleType (int order)
{
  switch(order) {
  case 0 : return MSH_TRI_1;
  case 1 : return MSH_TRI_3;
  case 2 : return MSH_TRI_6;
  case 3 : return MSH_TRI_10;
  case 4 : return MSH_TRI_15;
  case 5 : return MSH_TRI_21;
  case 6 : return MSH_TRI_28;
  case 7 : return MSH_TRI_36;
  case 8 : return MSH_TRI_45;
  case 9 : return MSH_TRI_55;
  case 10 : return MSH_TRI_66;
  default : Msg::Error("triangle order %i unknown", order); return 0;
  }
}

static int getQuadType (int order)
{
  switch(order) {
  case 0 : return MSH_QUA_1;
  case 1 : return MSH_QUA_4;
  case 2 : return MSH_QUA_9;
  case 3 : return MSH_QUA_16;
  case 4 : return MSH_QUA_25;
  case 5 : return MSH_QUA_36;
  case 6 : return MSH_QUA_49;
  case 7 : return MSH_QUA_64;
  case 8 : return MSH_QUA_81;
  case 9 : return MSH_QUA_100;
  case 10 : return MSH_QUA_121;
  default : Msg::Error("quad order %i unknown", order); return 0;
  }
}

static int getLineType (int order)
{
  switch(order) {
  case 0 : return MSH_LIN_1;
  case 1 : return MSH_LIN_2;
  case 2 : return MSH_LIN_3;
  case 3 : return MSH_LIN_4;
  case 4 : return MSH_LIN_5;
  case 5 : return MSH_LIN_6;
  case 6 : return MSH_LIN_7;
  case 7 : return MSH_LIN_8;
  case 8 : return MSH_LIN_9;
  case 9 : return MSH_LIN_10;
  case 10 : return MSH_LIN_11;
  default : Msg::Error("line order %i unknown", order); return 0;
  }
}

double MElement::integrateCirc(double val[], int edge, int pOrder, int order)
{
  if(edge > getNumEdges() - 1){
    Msg::Error("No edge %d for this element", edge);
    return 0;
  }

  std::vector<MVertex*> v;
  getEdgeVertices(edge, v);
  MElementFactory f;
  int type = getLineType(getPolynomialOrder());
  MElement* ee = f.create(type, v);

  double intv[3];
  for(int i = 0; i < 3; i++){
    intv[i] = ee->integrate(&val[i], pOrder, 3, order);
  }
  delete ee;

  double t[3] = {v[1]->x() - v[0]->x(), v[1]->y() - v[0]->y(), v[1]->z() - v[0]->z()};
  norme(t);
  double result;
  prosca(t, intv, &result);
  return result;
}

double MElement::integrateFlux(double val[], int face, int pOrder, int order)
{
  if(face > getNumFaces() - 1){
    Msg::Error("No face %d for this element", face);
    return 0;
  }
  std::vector<MVertex*> v;
  getFaceVertices(face, v);
  MElementFactory f;
  int type = 0;
  switch(getType()) {
  case TYPE_TRI : type = getTriangleType(getPolynomialOrder()); break;
  case TYPE_TET : type = getTriangleType(getPolynomialOrder()); break;
  case TYPE_QUA : type = getQuadType(getPolynomialOrder());     break;
  case TYPE_HEX : type = getQuadType(getPolynomialOrder());     break;
  case TYPE_PYR :
    if(face < 4) type = getTriangleType(getPolynomialOrder());
    else type = getQuadType(getPolynomialOrder());
    break;
  case TYPE_PRI :
    if(face < 2) type = getTriangleType(getPolynomialOrder());
    else type = getQuadType(getPolynomialOrder());
    break;
  default: type = 0; break;
  }
  MElement* fe = f.create(type, v);

  double intv[3];
  for(int i = 0; i < 3; i++){
    intv[i] = fe->integrate(&val[i], pOrder, 3, order);
  }
  delete fe;

  double n[3];
  normal3points(v[0]->x(), v[0]->y(), v[0]->z(),
                v[1]->x(), v[1]->y(), v[1]->z(),
                v[2]->x(), v[2]->y(), v[2]->z(), n);
  double result;
  prosca(n, intv, &result);
  return result;
}

void MElement::writeMSH(FILE *fp, bool binary, int elementary,
                        std::vector<short> *ghosts)
{
  int num = getNum();
  int type = getTypeForMSH();
  if(!type) return;

  // if necessary, change the ordering of the vertices to get positive volume
  setVolumePositive();

  std::vector<int> info;
  info.push_back(0);
  info.push_back(elementary);
  if(getParent())
    info.push_back(getParent()->getNum());
  if(getPartition()){
    if(ghosts){
      info.push_back(1 + ghosts->size());
      info.push_back(getPartition());
      info.insert(info.end(), ghosts->begin(), ghosts->end());
    }
    else{
      info.push_back(1);
      info.push_back(getPartition());
    }
  }
  info[0] = info.size() - 1;
  std::vector<int> verts;
  getVerticesIdForMSH(verts);
  info.insert(info.end(), verts.begin(), verts.end());

  if(!binary){
    fprintf(fp, "%d %d", num, type);
    for(unsigned int i = 0; i < info.size(); i++)
      fprintf(fp, " %d", info[i]);
    fprintf(fp, "\n");
  }
  else{
    fwrite(&num, sizeof(int), 1, fp);
    fwrite(&type, sizeof(int), 1, fp);
    fwrite(&info[0], sizeof(int), info.size(), fp);
  }
}

void MElement::writeMSH2(FILE *fp, double version, bool binary, int num,
                         int elementary, int physical, int parentNum,
                         int dom1Num, int dom2Num, std::vector<short> *ghosts)
{
  int type = getTypeForMSH();

  if(!type) return;

  // if necessary, change the ordering of the vertices to get positive
  // volume
  setVolumePositive();
  int n = getNumVerticesForMSH();
  int par = (parentNum) ? 1 : 0;
  int dom = (dom1Num) ? 2 : 0;
  bool poly = (type == MSH_POLYG_ || type == MSH_POLYH_ || type == MSH_POLYG_B);

  // if polygon loop over children (triangles and tets)
  if(CTX::instance()->mesh.saveTri){
    if(poly){
      for (int i = 0; i < getNumChildren() ; i++){
         MElement *t = getChild(i);
         t->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
      }
      return;
    }
    if(type == MSH_TRI_B){
      MTriangle *t = new MTriangle(getVertex(0), getVertex(1), getVertex(2));
      t->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
      delete t;
      return;
    }
    if(type == MSH_LIN_B || type == MSH_LIN_C){
      MLine *l = new MLine(getVertex(0), getVertex(1));
      l->writeMSH2(fp, version, binary, num++, elementary, physical, 0, 0, 0, ghosts);
      delete l;
      return;
    }
  }

  if(!binary){
    fprintf(fp, "%d %d", num ? num : _num, type);
    if(version < 2.0)
      fprintf(fp, " %d %d %d", abs(physical), elementary, n);
    else if (version < 2.2)
      fprintf(fp, " %d %d %d", abs(physical), elementary, _partition);
    else if(!_partition && !par && !dom)
      fprintf(fp, " %d %d %d", 2 + par + dom, abs(physical), elementary);
    else if(!ghosts)
      fprintf(fp, " %d %d %d 1 %d", 4 + par + dom, abs(physical), elementary, _partition);
    else{
      int numGhosts = ghosts->size();
      fprintf(fp, " %d %d %d %d %d", 4 + numGhosts + par + dom, abs(physical),
              elementary, 1 + numGhosts, _partition);
      for(unsigned int i = 0; i < ghosts->size(); i++)
        fprintf(fp, " %d", -(*ghosts)[i]);
    }
    if(version >= 2.0 && par)
      fprintf(fp, " %d", parentNum);
    if(version >= 2.0 && dom)
      fprintf(fp, " %d %d", dom1Num, dom2Num);
    if(version >= 2.0 && poly)
      fprintf(fp, " %d", n);
  }
  else{
    int numTags, numGhosts = 0;
    if(!_partition) numTags = 2;
    else if(!ghosts) numTags = 4;
    else{
      numGhosts = ghosts->size();
      numTags = 4 + numGhosts;
    }
    numTags += par;
    // we write elements in blobs of single elements; this will lead
    // to suboptimal reads, but it's much simpler when the number of
    // tags change from element to element (third-party codes can
    // still write MSH file optimized for reading speed, by grouping
    // elements with the same number of tags in blobs)
    int blob[60] = {type, 1, numTags, num ? num : _num, abs(physical), elementary,
                    1 + numGhosts, _partition};
    if(ghosts)
      for(int i = 0; i < numGhosts; i++) blob[8 + i] = -(*ghosts)[i];
    if(par) blob[8 + numGhosts] = parentNum;
    if(poly) Msg::Error("Unable to write polygons/polyhedra in binary files.");
    fwrite(blob, sizeof(int), 4 + numTags, fp);
  }

  if(physical < 0) revert();

  std::vector<int> verts;
  getVerticesIdForMSH(verts);

  if(!binary){
    for(int i = 0; i < n; i++)
      fprintf(fp, " %d", verts[i]);
    fprintf(fp, "\n");
  }
  else{
    fwrite(&verts[0], sizeof(int), n, fp);
  }

  if(physical < 0) revert();
}

void MElement::writePOS(FILE *fp, bool printElementary, bool printElementNumber,
                        bool printGamma, bool printEta, bool printRho,
                        bool printDisto, double scalingFactor, int elementary)
{
  const char *str = getStringForPOS();
  if(!str) return;

  setVolumePositive();
  int n = getNumVertices();
  fprintf(fp, "%s(", str);
  for(int i = 0; i < n; i++){
    if(i) fprintf(fp, ",");
    fprintf(fp, "%g,%g,%g", getVertex(i)->x() * scalingFactor,
            getVertex(i)->y() * scalingFactor, getVertex(i)->z() * scalingFactor);
  }
  fprintf(fp, "){");
  bool first = true;
  if(printElementary){
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%d", elementary);
    }
  }
  if(printElementNumber){
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%d", getNum());
    }
  }
  if(printGamma){
    double gamma = gammaShapeMeasure();
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%g", gamma);
    }
  }
  if(printEta){
    double eta = etaShapeMeasure();
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%g", eta);
    }
  }
  if(printRho){
    double rho = rhoShapeMeasure();
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%g", rho);
    }
  }
  if(printDisto){
    double disto = distoShapeMeasure();
    for(int i = 0; i < n; i++){
      if(first) first = false; else fprintf(fp, ",");
      fprintf(fp, "%g", disto);
    }
  }
  fprintf(fp, "};\n");
}

void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
{
  if(getType() != TYPE_TRI && getType() != TYPE_QUA) return;
  int qid[3] = {0, 2, 3};
  SVector3 n = getFace(0).normal();
  if(!binary){
    fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
    fprintf(fp, "  outer loop\n");
    for(int j = 0; j < 3; j++)
      fprintf(fp, "    vertex %g %g %g\n",
              getVertex(j)->x() * scalingFactor,
              getVertex(j)->y() * scalingFactor,
              getVertex(j)->z() * scalingFactor);
    fprintf(fp, "  endloop\n");
    fprintf(fp, "endfacet\n");
    if(getNumVertices() == 4){
      fprintf(fp, "facet normal %g %g %g\n", n[0], n[1], n[2]);
      fprintf(fp, "  outer loop\n");
      for(int j = 0; j < 3; j++)
        fprintf(fp, "    vertex %g %g %g\n",
                getVertex(qid[j])->x() * scalingFactor,
                getVertex(qid[j])->y() * scalingFactor,
                getVertex(qid[j])->z() * scalingFactor);
      fprintf(fp, "  endloop\n");
      fprintf(fp, "endfacet\n");
    }
  }
  else{
    char data[50];
    float *coords = (float*)data;
    coords[0] = (float)n[0];
    coords[1] = (float)n[1];
    coords[2] = (float)n[2];
    for(int j = 0; j < 3; j++){
      coords[3 + 3 * j] = (float)(getVertex(j)->x() * scalingFactor);
      coords[3 + 3 * j + 1] = (float)(getVertex(j)->y() * scalingFactor);
      coords[3 + 3 * j + 2] = (float)(getVertex(j)->z() * scalingFactor);
    }
    data[48] = data[49] = 0;
    fwrite(data, sizeof(char), 50, fp);
    if(getNumVertices() == 4){
      for(int j = 0; j < 3; j++){
        coords[3 + 3 * j] = (float)(getVertex(qid[j])->x() * scalingFactor);
        coords[3 + 3 * j + 1] = (float)(getVertex(qid[j])->y() * scalingFactor);
        coords[3 + 3 * j + 2] = (float)(getVertex(qid[j])->z() * scalingFactor);
      }
      fwrite(data, sizeof(char), 50, fp);
    }
  }
}
void MElement::writePLY2(FILE *fp)
{
  setVolumePositive();
  fprintf(fp, "3 ");
  for(int i = 0; i < getNumVertices(); i++)
    fprintf(fp, " %d", getVertex(i)->getIndex() - 1);
  fprintf(fp, "\n");
}

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

void MElement::writeVTK(FILE *fp, bool binary, bool bigEndian)
{
  if(!getTypeForVTK()) return;

  setVolumePositive();

  int n = getNumVertices();
  if(binary){
    int verts[60];
    verts[0] = n;
    for(int i = 0; i < n; i++)
      verts[i + 1] = getVertexVTK(i)->getIndex() - 1;
    // VTK always expects big endian binary data
    if(!bigEndian) SwapBytes((char*)verts, sizeof(int), n + 1);