double approximationError(simpleFunction<double> &f, MElement *element)
{
  std::vector<double> VALS(element->getNumVertices());

  for(std::size_t i = 0; i < element->getNumVertices(); i++) {
    MVertex *v = element->getVertex(i);
    VALS[i] = f(v->x(), v->y(), v->z());
  }

  int npts;
  IntPt *pts;
  element->getIntegrationPoints(2 * element->getPolynomialOrder() + 2, &npts,
                                &pts);
  double errSqr = 0.0;
  for(int k = 0; k < npts; k++) {
    const double u = pts[k].pt[0];
    const double v = pts[k].pt[1];
    const double w = pts[k].pt[2];
    SPoint3 p;
    element->pnt(u, v, w, p);
    const double Jac = element->getJacobianDeterminant(u, v, w);
    const double C = element->interpolate(&VALS[0], u, v, w);
    const double F = f(p.x(), p.y(), p.z());
    errSqr += pts[k].weight * Jac * std::pow(C - F, 2);
  }
  return std::sqrt(errSqr);
}

// Bezier
static Vertex InterpolateBezier(Curve *Curve, double u, int derivee)
{
  int NbCurves = (List_Nbr(Curve->Control_Points) - 1) / 3;
  int iCurve = (int)floor(u * (double)NbCurves);
  if(iCurve >= NbCurves) iCurve = NbCurves - 1; // u = 1
  if(iCurve <= 0) iCurve = 0;
  double t1 = (double)(iCurve) / (double)(NbCurves);
  double t2 = (double)(iCurve+1) / (double)(NbCurves);
  double t = (u - t1) / (t2 - t1);
  Vertex *v[4];
  for(int i = 0; i < 4; i++) {
    List_Read(Curve->Control_Points, iCurve * 3 + i , &v[i]);
  }

  if(Curve->geometry){
    SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, t1, t2, Curve->geometry,derivee);
    SPoint3 pt = Curve->geometry->point(pp);
    Vertex V;
    V.Pos.X = pt.x();
    V.Pos.Y = pt.y();
    V.Pos.Z = pt.z();
    return V;
  }
  else
    return InterpolateCubicSpline(v, t, Curve->mat, derivee, t1, t2);
}

bool gmshFace::containsPoint(const SPoint3 &pt) const
{
  if(s->Typ == MSH_SURF_PLAN){
    // OK to use the normal from the mean plane here: we compensate
    // for the (possibly wrong) orientation at the end
    double n[3] = {meanPlane.a, meanPlane.b, meanPlane.c};
    norme(n);
    double angle = 0.;
    double v[3] = {pt.x(), pt.y(), pt.z()};
    for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
      Curve *c;
      List_Read(s->Generatrices, i, &c);
      int N = (c->Typ == MSH_SEGM_LINE) ? 1 : 10;
      for(int j = 0; j < N; j++) {
        double u1 = (double)j / (double)N;
        double u2 = (double)(j + 1) / (double)N;
        Vertex p1 = InterpolateCurve(c, u1, 0);
        Vertex p2 = InterpolateCurve(c, u2, 0);
        double v1[3] = {p1.Pos.X, p1.Pos.Y, p1.Pos.Z};
        double v2[3] = {p2.Pos.X, p2.Pos.Y, p2.Pos.Z};
        angle += angle_plan(v, v1, v2, n);
      }
    }
    // we're inside if angle equals 2 * pi
    if(fabs(angle) > 2 * M_PI - 0.5 && fabs(angle) < 2 * M_PI + 0.5)
      return true;
    return false;
  }

  return false;
}

double Filler::improvement(GEntity* ge,MElementOctree* octree,SPoint3 point,double h1,SVector3 direction){
  double x,y,z;
  double average;
  double h2;
  double coeffA,coeffB;

  x = point.x() + h1*direction.x();
  y = point.y() + h1*direction.y();
  z = point.z() + h1*direction.z();
  
  if(inside_domain(octree,x,y,z)){
    h2 = get_size(x,y,z);
  }
  else h2 = h1;

  coeffA = 1.0;
  coeffB = 0.16;
  
  if(h2>h1){
    average = coeffA*h1 + (1.0-coeffA)*h2;
  }
  else{
    average = coeffB*h1 + (1.0-coeffB)*h2;
  }
	
  return average;
}

GPoint OCCEdge::closestPoint(const SPoint3 &qp, double &param) const
{
  if(curve.IsNull()){
    Msg::Error("OCC curve is null in closestPoint");
    return GPoint(0, 0);
  }

  gp_Pnt pnt(qp.x(), qp.y(), qp.z());
  GeomAPI_ProjectPointOnCurve proj(pnt, curve, s0, s1);

  if(!proj.NbPoints()){
    Msg::Error("OCC ProjectPointOnCurve failed");
    return GPoint(0, 0);
  }

  param = proj.LowerDistanceParameter();

  if(param < s0 || param > s1){
    Msg::Error("Point projection is out of edge bounds");
    return GPoint(0, 0);
  }

  pnt = proj.NearestPoint();
  return GPoint(pnt.X(), pnt.Y(), pnt.Z(), this, param);
}

SPoint2 fourierFace::parFromPoint(const SPoint3 &p, bool onSurface) const
{
  double u, v, x, y, z;
  x = p.x(); y = p.y(); z = p.z();
  face->Inverse(x,y,z,u,v);

  return SPoint2(u, v);
}

void elasticitySolver::computeEffectiveStiffness(std::vector<double> stiff)
{
  double st[6] = {0., 0., 0., 0., 0., 0.};
  double volTot = 0.;
  for(std::size_t i = 0; i < elasticFields.size(); ++i) {
    double E = elasticFields[i]._e;
    double nu = elasticFields[i]._nu;
    SolverField<SVector3> Field(pAssembler, LagSpace);
    for(groupOfElements::elementContainer::const_iterator it =
          elasticFields[i].g->begin();
        it != elasticFields[i].g->end(); ++it) {
      MElement *e = *it;
      double vol = e->getVolume() * e->getVolumeSign();
      int nbVertex = e->getNumVertices();
      std::vector<SVector3> val(nbVertex);

      double valx[256];
      double valy[256];
      double valz[256];
      for(int k = 0; k < nbVertex; k++) {
        MVertex *v = e->getVertex(k);
        MPoint p(v);
        Field.f(&p, 0, 0, 0, val[k]);
        valx[k] = val[k](0);
        valy[k] = val[k](1);
        valz[k] = val[k](2);
      }

      double gradux[3];
      double graduy[3];
      double graduz[3];
      SPoint3 center = e->barycenterUVW();
      double u = center.x(), v = center.y(), w = center.z();
      e->interpolateGrad(valx, u, v, w, gradux);
      e->interpolateGrad(valy, u, v, w, graduy);
      e->interpolateGrad(valz, u, v, w, graduz);

      double eps[6] = {gradux[0],
                       graduy[1],
                       graduz[2],
                       0.5 * (gradux[1] + graduy[0]),
                       0.5 * (gradux[2] + graduz[0]),
                       0.5 * (graduy[2] + graduz[1])};

      double A = E / (1. + nu);
      double B = A * (nu / (1. - 2 * nu));
      double trace = eps[0] + eps[1] + eps[2];
      st[0] += (A * eps[0] + B * trace) * vol;
      st[1] += (A * eps[1] + B * trace) * vol;
      st[2] += (A * eps[2] + B * trace) * vol;
      st[3] += (A * eps[3]) * vol;
      st[4] += (A * eps[4]) * vol;
      st[5] += (A * eps[5]) * vol;
      volTot += vol;
    }
  }
  for(int i = 0; i < 6; i++) stiff[i] = st[i] / volTot;
}

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

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

 bool operator()(const SPoint3 &b1, const SPoint3 &b2)
 {
   double p1[3] = {b1.x(), b1.y(), b1.z()};
   double p2[3] = {b2.x(), b2.y(), b2.z()};
   double c[3] = {v.x(), v.y(), v.z()};
   double a1 = myangle(c, p1);
   double a2 = myangle(c, p2);
   return a1 < a2;
 }

SPoint2 gmshFace::parFromPoint(const SPoint3 &qp, bool onSurface) const
{
  if(s->Typ == MSH_SURF_PLAN){
    double x, y, z, VX[3], VY[3];
    getMeanPlaneData(VX, VY, x, y, z);
    double u, v, vec[3] = {qp.x() - x, qp.y() - y, qp.z() - z};
    prosca(vec, VX, &u);
    prosca(vec, VY, &v);
    return SPoint2(u, v);
  }
  else{
    return GFace::parFromPoint(qp, onSurface);
  }
}

void Filler::print_node(Node* node,std::ofstream& file){
  double x,y,z;
  double x1,y1,z1;
  double x2,y2,z2;
  double x3,y3,z3;
  double x4,y4,z4;
  double x5,y5,z5;
  double x6,y6,z6;
  double h;
  Metric m;
  SPoint3 point;

  point = node->get_point();
  x = point.x();
  y = point.y();
  z = point.z();
  h = node->get_size();
  m = node->get_metric();

  x1 = x + k1*h*m.get_m11();
  y1 = y + k1*h*m.get_m21();
  z1 = z + k1*h*m.get_m31();

  x2 = x - k1*h*m.get_m11();
  y2 = y - k1*h*m.get_m21();
  z2 = z - k1*h*m.get_m31();

  x3 = x + k1*h*m.get_m12();
  y3 = y + k1*h*m.get_m22();
  z3 = z + k1*h*m.get_m32();

  x4 = x - k1*h*m.get_m12();
  y4 = y - k1*h*m.get_m22();
  z4 = z - k1*h*m.get_m32();

  x5 = x + k1*h*m.get_m13();
  y5 = y + k1*h*m.get_m23();
  z5 = z + k1*h*m.get_m33();

  x6 = x - k1*h*m.get_m13();
  y6 = y - k1*h*m.get_m23();
  z6 = z - k1*h*m.get_m33();

  print_segment(SPoint3(x,y,z),SPoint3(x1,y1,z1),file);
  print_segment(SPoint3(x,y,z),SPoint3(x2,y2,z2),file);
  print_segment(SPoint3(x,y,z),SPoint3(x3,y3,z3),file);
  print_segment(SPoint3(x,y,z),SPoint3(x4,y4,z4),file);
  print_segment(SPoint3(x,y,z),SPoint3(x5,y5,z5),file);
  print_segment(SPoint3(x,y,z),SPoint3(x6,y6,z6),file);
}

static void drawTangents(drawContext *ctx, std::vector<T*> &elements)
{
  glColor4ubv((GLubyte *) & CTX::instance()->color.mesh.tangents);
  for(unsigned int i = 0; i < elements.size(); i++){
    MElement *ele = elements[i];
    if(!isElementVisible(ele)) continue;
    SVector3 t = ele->getEdge(0).tangent();
    for(int j = 0; j < 3; j++)
      t[j] *= CTX::instance()->mesh.tangents * ctx->pixel_equiv_x / ctx->s[j];
    SPoint3 pc = ele->barycenter();
    ctx->drawVector(CTX::instance()->vectorType, 0, pc.x(), pc.y(), pc.z(),
                    t[0], t[1], t[2], CTX::instance()->mesh.light);
  }
}

void printJacobians(GModel *m, const char *nm)
{
  const int n = 100;
  double D[n][n], X[n][n], Y[n][n], Z[n][n];

  FILE *f = Fopen(nm,"w");
  fprintf(f,"View \"\"{\n");
  for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
    for(unsigned int j = 0; j < (*it)->triangles.size(); j++){
      MTriangle *t = (*it)->triangles[j];
      for(int i = 0; i < n; i++){
        for(int k = 0; k < n - i; k++){
          SPoint3 pt;
          double u = (double)i / (n - 1);
          double v = (double)k / (n - 1);
          t->pnt(u, v, 0, pt);
          D[i][k] = 0.; //mesh_functional_distorsion_2D(t, u, v);
          //X[i][k] = u;
          //Y[i][k] = v;
          //Z[i][k] = 0.0;
          X[i][k] = pt.x();
          Y[i][k] = pt.y();
          Z[i][k] = pt.z();
        }
      }
      for(int i= 0; i < n -1; i++){
        for(int k = 0; k < n - i -1; k++){
          fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
                  X[i][k],Y[i][k],Z[i][k],
                  X[i+1][k],Y[i+1][k],Z[i+1][k],
                  X[i][k+1],Y[i][k+1],Z[i][k+1],
                  D[i][k],
                  D[i+1][k],
                  D[i][k+1]);
          if (i != n-2 && k != n - i -2)
            fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%22.15E,%22.15E,%22.15E};\n",
                    X[i+1][k],Y[i+1][k],Z[i+1][k],
                    X[i+1][k+1],Y[i+1][k+1],Z[i+1][k+1],
                    X[i][k+1],Y[i][k+1],Z[i][k+1],
                    D[i+1][k],
                    D[i+1][k+1],
                    D[i][k+1]);
        }
      }
    }
  }
  fprintf(f,"};\n");
  fclose(f);
}

static void drawElementLabels(drawContext *ctx, GEntity *e,
                              std::vector<T*> &elements, int forceColor=0,
                              unsigned int color=0)
{
  unsigned col = forceColor ? color : getColorByEntity(e);
  glColor4ubv((GLubyte *) & col);

  int labelStep = CTX::instance()->mesh.labelSampling;
  if(labelStep <= 0) labelStep = 1;

  for(unsigned int i = 0; i < elements.size(); i++){
    MElement *ele = elements[i];
    if(!isElementVisible(ele)) continue;
    if(i % labelStep == 0) {
      SPoint3 pc = ele->barycenter();
      char str[256];
      if(CTX::instance()->mesh.labelType == 4)
        sprintf(str, "(%g,%g,%g)", pc.x(), pc.y(), pc.z());
      else if(CTX::instance()->mesh.labelType == 3)
        sprintf(str, "%d", ele->getPartition());
      else if(CTX::instance()->mesh.labelType == 2){
        int np = e->physicals.size();
        int p = np ? e->physicals[np - 1] : 0;
        sprintf(str, "%d", p);
      }
      else if(CTX::instance()->mesh.labelType == 1)
        sprintf(str, "%d", e->tag());
      else
        sprintf(str, "%d", ele->getNum());
      glRasterPos3d(pc.x(), pc.y(), pc.z());
      ctx->drawString(str);
    }
  }
}

bool iSRuledSurfaceASphere(Surface *s, SPoint3 &center, double &radius)
{
  if(s->Typ != MSH_SURF_REGL && s->Typ != MSH_SURF_TRIC) return false;

  bool isSphere = true;
  Vertex *O = 0;
  Curve *C[4] = {0, 0, 0, 0};
  for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++)
    List_Read(s->Generatrices, i, &C[i]);

  if(List_Nbr(s->InSphereCenter)) {
    // it's on a sphere: get the center
    List_Read(s->InSphereCenter, 0, &O);
  }
  else{
    // try to be intelligent (hum)
    for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++) {
      if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV){
        isSphere = false;
      }
      else if(isSphere){
        if(!i){
          List_Read(C[i]->Control_Points, 1, &O);
          ((double *)center)[0] = O->Pos.X;
          ((double *)center)[1] = O->Pos.Y;
          ((double *)center)[2] = O->Pos.Z;
        }
        else{
          Vertex *tmp;
          List_Read(C[i]->Control_Points, 1, &tmp);
          if(compareVertex(&O, &tmp))
            isSphere = false;
        }
      }
    }
  }
  if (isSphere && C[0]){
    Vertex *p = C[0]->beg;
    radius = sqrt ((p->Pos.X - center.x())+
                   (p->Pos.Y - center.y())+
                   (p->Pos.Z - center.z()));
  }

  return isSphere;
}

static bool _kaboom(fullVector<double> &uvt,
                    fullVector<double> &res, void *_data)
{
  intersectCurveSurfaceData *data = (intersectCurveSurfaceData*)_data;
  SPoint3 s = data->s(uvt(0), uvt(1));
  SPoint3 c = data->c(uvt(2));
  res(0) = s.x() - c.x();
  res(1) = s.y() - c.y();
  res(2) = s.z() - c.z();
  return true;
}

template<class T1> void LoadTerm<T1>::get(MElement *ele, int npts, IntPt *GP, fullVector<double> &m) const
{
  if(ele->getParent()) ele = ele->getParent();
  int nbFF = LinearTerm<T1>::space1.getNumKeys(ele);
  double jac[3][3];
  m.resize(nbFF);
  m.scale(0.);
  for(int i = 0; i < npts; i++){
    const double u = GP[i].pt[0]; const double v = GP[i].pt[1]; const double w = GP[i].pt[2];
    const double weight = GP[i].weight; const double detJ = ele->getJacobian(u, v, w, jac);
    std::vector<typename TensorialTraits<T1>::ValType> Vals;
    LinearTerm<T1>::space1.f(ele, u, v, w, Vals);
    SPoint3 p;
    ele->pnt(u, v, w, p);
    typename TensorialTraits<T1>::ValType load = (*Load)(p.x(), p.y(), p.z());
    for(int j = 0; j < nbFF ; ++j)
    {
      m(j) += dot(Vals[j], load) * weight * detJ;
    }
  }
}

PView *thermicSolver::buildErrorEstimateView(const std::string &errorFileName,
                                             simpleFunction<double> *sol)
{
  std::cout << "build Error View" << std::endl;
  std::map<int, std::vector<double> > data;

  SolverField<double> solField(pAssembler, LagSpace);
  for(std::size_t i = 0; i < thermicFields.size(); ++i) {
    for(groupOfElements::elementContainer::const_iterator it =
          thermicFields[i].g->begin();
        it != thermicFields[i].g->end(); ++it) {
      MElement *e = *it;
      int npts;
      IntPt *GP;
      double jac[3][3];
      int integrationOrder = 2 * (e->getPolynomialOrder() + 5);
      e->getIntegrationPoints(integrationOrder, &npts, &GP);
      double val = 0.0;
      for(int j = 0; j < npts; j++) {
        double u = GP[j].pt[0];
        double v = GP[j].pt[1];
        double w = GP[j].pt[2];
        double weight = GP[j].weight;
        double detJ = fabs(e->getJacobian(u, v, w, jac));
        SPoint3 p;
        e->pnt(u, v, w, p);
        double FEMVALUE;
        solField.f(e, u, v, w, FEMVALUE);
        double diff = (*sol)(p.x(), p.y(), p.z()) - FEMVALUE;
        val += diff * diff * detJ * weight;
      }
      std::vector<double> vec;
      vec.push_back(sqrt(val));
      data[e->getNum()] = vec;
    }
  }

  PView *pv = new PView(errorFileName, "ElementData", pModel, data, 0.0, 1);
  return pv;
}

bool Filler::far_from_boundary(MElementOctree* octree,Node* node){
  double x,y,z;
  double h;
  SPoint3 point;
  MElement *e1,*e2,*e3,*e4,*e5,*e6;

  point = node->get_point();
  x = point.x();
  y = point.y();
  z = point.z();
  h = node->get_size();

  e1 = (MElement*)octree->find(x+k2*h,y,z,3,true);
  e2 = (MElement*)octree->find(x-k2*h,y,z,3,true);
  e3 = (MElement*)octree->find(x,y+k2*h,z,3,true);
  e4 = (MElement*)octree->find(x,y-k2*h,z,3,true);
  e5 = (MElement*)octree->find(x,y,z+k2*h,3,true);
  e6 = (MElement*)octree->find(x,y,z-k2*h,3,true);

  if(e1!=NULL && e2!=NULL && e3!=NULL && e4!=NULL && e5!=NULL && e6!=NULL) return 1;
  else return 0;
}