mettools.c

Go to the documentation of this file.

/* =============================================================================
**  This file is part of the mmg software package for the tetrahedral
**  mesh modification.
**  Copyright (c) Bx INP/CNRS/Inria/UBordeaux/UPMC, 2004-
**
**  mmg is free software: you can redistribute it and/or modify it
**  under the terms of the GNU Lesser General Public License as published
**  by the Free Software Foundation, either version 3 of the License, or
**  (at your option) any later version.
**
**  mmg 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 Lesser General Public
**  License for more details.
**
**  You should have received a copy of the GNU Lesser General Public
**  License and of the GNU General Public License along with mmg (in
**  files COPYING.LESSER and COPYING). If not, see
**  <http://www.gnu.org/licenses/>. Please read their terms carefully and
**  use this copy of the mmg distribution only if you accept them.
** =============================================================================
*/
 
#include "mmgcommon_private.h"
#include "inlined_functions_private.h"
 
static inline
void MMG5_eigenvmat_buildsym(MMG5_pMesh mesh,int8_t dim,double m[],
                             double lambda[],double v[]) {
  int8_t i,j,k,ij;
 
  /* Storage of a matrix as a one-dimensional array: (dim+1)*dim/2 entries for a
   * symmetric matrix (loop on each symmetric entry only once). */
  ij = 0;
 
  /* Loop on matrix rows */
  for( i = 0; i < dim; i++ ) {
    /* Loop on the upper-triangular part of the matrix */
    for( j = i; j < dim; j++ ) {
      /* Initialize matrix entry */
      m[ij] = 0.0;
      /* Compute matrix entry as the recomposition of eigenvalues and
       * eigenvectors:
       *
       * M_{ij} = \sum_{k,l} V_{ik} Lambda_{kl} V_{jl} =
       *        = \sum_{k} lambda_{k} V_{ik} V_{jk}
       *
       * Eigenvectors are stored as rows in Mmg (not columns) so their indices
       * have to be exchanged when implementing the above formula. */
      for( k = 0; k < dim; k++ ) {
        m[ij] += lambda[k]*v[k*dim+i]*v[k*dim+j];
      }
      /* Go to next entry */
      ++ij;
    }
  }
  assert( ij == (dim+1)*dim/2 );
}
 
static inline
void MMG5_eigenvmat_buildnonsym(MMG5_pMesh mesh,int8_t dim,double m[],
                                double lambda[],double v[],double iv[]){
  int8_t i,j,k,ij;
 
  /* Storage of a matrix as a one-dimensional array: dim^2 entries for a
   * non-symmetric matrix. */
  ij = 0;
 
  /* Loop on matrix rows */
  for( i = 0; i < dim; i++ ) {
    /* Loop on matrix columns */
    for( j = 0; j < dim; j++ ) {
      /* Initialize matrix entry */
      m[ij] = 0.0;
      /* Compute matrix entry as the recomposition of eigenvalues and
       * eigenvectors:
       *
       * M_{ij} = \sum_{k,l} V_{ik} Lambda_{kl} V^{-1}_{lj} =
       *        = \sum_{k} lambda_{k} V_{ik} V^{-1}_{kj}
       *
       * Eigenvectors are stored as rows in Mmg (not columns) so their indices
       * have to be exchanged when implementing the above formula. */
      for( k = 0; k < dim; k++ ) {
        m[ij] += lambda[k]*v[k*dim+i]*iv[j*dim+k];
      }
      /* Go to next entry */
      ++ij;
    }
  }
  assert( ij == dim*dim );
}
 
int MMG5_eigenvmatsym2d(MMG5_pMesh mesh,double m[],double lambda[],double v[][2]) {
 
  /* Build matrix */
  MMG5_eigenvmat_buildsym(mesh,2,m,lambda,(double *)v);
 
  return 1;
}
 
int MMG5_eigenvmatsym3d(MMG5_pMesh mesh,double m[],double lambda[],double v[][3]) {
 
  /* Build matrix */
  MMG5_eigenvmat_buildsym(mesh,3,m,lambda,(double *)v);
 
  return 1;
}
 
int MMG5_eigenvmatnonsym2d(MMG5_pMesh mesh,double m[],double lambda[],double v[][2]) {
  double iv[2][2];
 
  /* Compute the inverse of the eigenvectors matrix */
  if( !MMG5_invmat22(v,iv) )
    return 0;
 
  /* Build matrix */
  MMG5_eigenvmat_buildnonsym(mesh,2,m,lambda,(double *)v,(double *)iv);
 
  return 1;
}
 
int MMG5_eigenvmatnonsym3d(MMG5_pMesh mesh,double m[],double lambda[],double v[][3]) {
  double iv[3][3];
 
  /* Compute the inverse of the eigenvectors matrix */
  if( !MMG5_invmat33(v,iv) )
    return 0;
 
  /* Build matrix */
  MMG5_eigenvmat_buildnonsym(mesh,3,m,lambda,(double *)v,(double *)iv);
 
  return 1;
}
 
int MMG5_eigenvmat_check(MMG5_pMesh mesh,int8_t dim,int8_t symmat,double m[],
                         double mnew[]) {
 
  /* Compute eigendecomposition, recompose matrix from eigendecomposition */
  if( dim == 2 ) {
    double lambda[2],v[2][2];
 
    if( !MMG5_eigenv2d(symmat,m,lambda,v) )
      return 0;
 
    if( symmat ) {
      if( !MMG5_eigenvmatsym2d(mesh,mnew,lambda,v) )
        return 0;
    } else {
      if( !MMG5_eigenvmatnonsym2d(mesh,mnew,lambda,v) )
        return 0;
    }
 
  } else if( dim == 3 ) {
    double lambda[3],v[3][3];
 
    if( !MMG5_eigenv3d(symmat,m,lambda,v) )
      return 0;
 
    if( symmat ) {
      if( !MMG5_eigenvmatsym3d(mesh,mnew,lambda,v) )
        return 0;
    } else {
      if( !MMG5_eigenvmatnonsym3d(mesh,mnew,lambda,v) )
        return 0;
    }
  }
 
  return 1;
}
 
void MMG5_sort_eigenv( int8_t dim,double *lambda,double *vp,
                       double *swap,int8_t *perm ) {
  MMG5_nsort(dim,lambda,perm);
  MMG5_nperm(dim,0,1,lambda,swap,perm);
  for( int8_t i = 0; i < dim; i++ )
    MMG5_nperm(dim,i,dim,vp,swap,perm);
}
 
int MMG5_test_eigenvmatsym2d(MMG5_pMesh mesh,double *mex,double lambdaex[],
                             double vpex[][2]) {
  double mnum[3],lambdanum[2],vpnum[2][2]; /* Numerical quantities */
  double swap[2],maxerr,err;
  int8_t perm[2] = {0,1}; /* eigenvalues permutation array */
 
 
  MMG5_eigenvmat_buildsym(mesh,2,mnum,lambdaex,(double *)vpex);
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(3,(double *)mex,(double *)mnum);
  if( maxerr > 10.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenv2d(1,mex,lambdanum,vpnum) )
    return 0;
 
  /* Naively sort eigenpairs in increasing order */
  MMG5_sort_eigenv(2,lambdanum,(double *)vpnum,swap,perm);
 
  /* Check eigenvalues error in norm inf */
  maxerr = MMG5_test_mat_error(2,(double *)lambdaex,(double *)lambdanum);
  if( maxerr > 10*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvalues: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  /* Check eigenvectors error through scalar product */
  maxerr = 0.;
  for( int8_t i = 0; i < 2; i++ ) {
    err = 0.;
    for( int8_t j = 0; j < 2; j++ )
      err += vpex[i][j] * vpnum[i][j];
    err = 1.-fabs(err);
    maxerr = MG_MAX(maxerr,err);
  }
  if( maxerr > MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvectors: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenvmat_check(mesh,2,1,mex,mnum) )
    return 0;
  maxerr = MMG5_test_mat_error(3,(double *)mex,(double *)mnum);
  if( maxerr > 10*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigendecomposition and recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  return 1;
}
 
int MMG5_test_eigenvmatnonsym2d(MMG5_pMesh mesh,double *mex,double lambdaex[],
                                double vpex[][2],double ivpex[][2]) {
  double mnum[4],lambdanum[2],vpnum[2][2]; /* Numerical quantities */
  double swap[2],maxerr,err;
  int8_t perm[2] = {0,1}; /* eigenvalues permutation array */
 
 
  MMG5_eigenvmat_buildnonsym(mesh,2,mnum,lambdaex,(double *)vpex,(double *)ivpex);
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(4,(double *)mex,(double *)mnum);
  if( maxerr > 100.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenv2d(0,mex,lambdanum,vpnum) )
    return 0;
 
  /* Naively sort eigenpairs in increasing order */
  MMG5_sort_eigenv(2,lambdanum,(double *)vpnum,swap,perm);
 
  /* Check eigenvalues error in norm inf */
  maxerr = MMG5_test_mat_error(2,(double *)lambdaex,(double *)lambdanum);
  if( maxerr > 10*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvalues: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  /* Check eigenvectors error through scalar product */
  maxerr = 0.;
  for( int8_t i = 0; i < 2; i++ ) {
    err = 0.;
    for( int8_t j = 0; j < 2; j++ )
      err += vpex[i][j] * vpnum[i][j];
    err = 1.-fabs(err);
    maxerr = MG_MAX(maxerr,err);
  }
  if( maxerr > MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvectors: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenvmat_check(mesh,2,0,mex,mnum) )
    return 0;
  maxerr = MMG5_test_mat_error(4,(double *)mex,(double *)mnum);
  if( maxerr > 100*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigendecomposition and recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  return 1;
}
 
int MMG5_test_eigenvmatsym3d(MMG5_pMesh mesh,double *mex,double lambdaex[],
                             double vpex[][3]) {
  double mnum[6],lambdanum[3],vpnum[3][3]; /* Numerical quantities */
  double swap[3],maxerr,err;
  int8_t perm[3] = {0,1,2}; /* eigenvalues permutation array */
 
 
  MMG5_eigenvmat_buildsym(mesh,3,mnum,lambdaex,(double *)vpex);
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(6,(double *)mex,(double *)mnum);
  if( maxerr > 100.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenv3d(1,mex,lambdanum,vpnum) )
    return 0;
 
  /* Naively sort eigenpairs in increasing order */
  MMG5_sort_eigenv(3,lambdanum,(double *)vpnum,swap,perm);
 
  /* Check eigenvalues error in norm inf */
  maxerr = MMG5_test_mat_error(3,(double *)lambdaex,(double *)lambdanum);
  if( maxerr > 10*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvalues: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  /* Check eigenvectors orthogonality (checking numerical eigenvectors against
   * the exact ones does not make sense in case of multiple eigenvalues) */
  maxerr = 0.;
  for( int8_t i = 0; i < 2; i++ ) {
    for( int8_t j = i+1; j < 3; j++ ) {
      err = 0.;
      for( int8_t k = 0; k < 3; k++ )
        err += vpnum[i][k] * vpnum[j][k];
      err = fabs(err);
      maxerr = MG_MAX(maxerr,err);
    }
  }
  if( maxerr > MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvectors orthogonality: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenvmat_check(mesh,3,1,mex,mnum) )
    return 0;
  maxerr = MMG5_test_mat_error(6,(double *)mex,(double *)mnum);
  if( maxerr > 100*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigendecomposition and recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  return 1;
}
 
int MMG5_test_eigenvmatnonsym3d(MMG5_pMesh mesh,double *mex,double lambdaex[],
                                double vpex[][3],double ivpex[][3]) {
  double mnum[9],lambdanum[3],vpnum[3][3]; /* Numerical quantities */
  double swap[3],maxerr;
  int8_t perm[3] = {0,1,2}; /* eigenvalues permutation array */
 
 
  MMG5_eigenvmat_buildnonsym(mesh,3,mnum,lambdaex,(double *)vpex,(double *)ivpex);
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(9,(double *)mex,(double *)mnum);
  if( maxerr > 1000.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
 
  if( !MMG5_eigenv3d(0,mex,lambdanum,vpnum) )
    return 0;
 
  /* Naively sort eigenpairs in increasing order */
  MMG5_sort_eigenv(3,lambdanum,(double *)vpnum,swap,perm);
 
  /* Check eigenvalues error in norm inf */
  maxerr = MMG5_test_mat_error(3,(double *)lambdaex,(double *)lambdanum);
  if( maxerr > 1000*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigenvalues: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  /* skip eigenvectors check */
 
 
  if( !MMG5_eigenvmat_check(mesh,3,0,mex,mnum) )
    return 0;
  maxerr = MMG5_test_mat_error(9,(double *)mex,(double *)mnum);
  if( maxerr > 1000*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error matrix eigendecomposition and recomposition: in function %s, max error %e\n",
      __func__,maxerr);
    return 0;
  }
 
  return 1;
}
 
inline int
MMG5_buildridmetfic(MMG5_pMesh mesh,double t[3],double n[3],double dtan,
                     double dv,double dn,double m[6]) {
  double u[3],r[3][3];
 
  u[0] = n[1]*t[2] - n[2]*t[1];
  u[1] = n[2]*t[0] - n[0]*t[2];
  u[2] = n[0]*t[1] - n[1]*t[0];
 
  /* If u = n1 ^ t, matrix of the desired metric in (t,u,n1) = diag(dtan,dv,dn)*/
  r[0][0] = t[0];  r[0][1] = u[0];  r[0][2] = n[0];
  r[1][0] = t[1];  r[1][1] = u[1];  r[1][2] = n[1];
  r[2][0] = t[2];  r[2][1] = u[2];  r[2][2] = n[2];
 
  m[0] = dtan*r[0][0]*r[0][0] + dv*r[0][1]*r[0][1] + dn*r[0][2]*r[0][2];
  m[1] = dtan*r[0][0]*r[1][0] + dv*r[0][1]*r[1][1] + dn*r[0][2]*r[1][2];
  m[2] = dtan*r[0][0]*r[2][0] + dv*r[0][1]*r[2][1] + dn*r[0][2]*r[2][2];
  m[3] = dtan*r[1][0]*r[1][0] + dv*r[1][1]*r[1][1] + dn*r[1][2]*r[1][2];
  m[4] = dtan*r[1][0]*r[2][0] + dv*r[1][1]*r[2][1] + dn*r[1][2]*r[2][2];
  m[5] = dtan*r[2][0]*r[2][0] + dv*r[2][1]*r[2][1] + dn*r[2][2]*r[2][2];
 
  return 1;
}
 
int MMG5_intmetsavedir(MMG5_pMesh mesh, double *m,double *n,double *mr) {
  int    i;
  double lambda[2],vp[2][2],siz,isqhmin;
 
  isqhmin = 1.0 / (mesh->info.hmin * mesh->info.hmin);
  MMG5_eigensym(m,lambda,vp);
 
  for (i=0; i<2; i++) {
    siz = n[0]*vp[i][0]*vp[i][0] + 2.0*n[1]*vp[i][0]*vp[i][1]
      + n[2]*vp[i][1]*vp[i][1];
    lambda[i] = MG_MAX(lambda[i],siz);
    lambda[i] = MG_MIN(lambda[i],isqhmin);
  }
  mr[0] = lambda[0]*vp[0][0]*vp[0][0] + lambda[1]*vp[1][0]*vp[1][0];
  mr[1] = lambda[0]*vp[0][0]*vp[0][1] + lambda[1]*vp[1][0]*vp[1][1];
  mr[2] = lambda[0]*vp[0][1]*vp[0][1] + lambda[1]*vp[1][1]*vp[1][1];
 
  return 1;
}
 
int MMG5_buildridmet(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np0,
                     double ux,double uy,double uz,double mr[6],
                     double r[3][3] ) {
  MMG5_pPoint  p0;
  MMG5_pxPoint go;
  double       ps1,ps2,*n1,*n2,*t,*m,dv,dn,u[3];
 
  p0 = &mesh->point[np0];
  if ( !(MG_GEO & p0->tag) )  return 0;
  m = &met->m[6*np0];
  t = &p0->n[0];
 
  /* Check that ridge tangent is not null */
  assert ( t[0]*t[0] + t[1]*t[1] + t[2]*t[2] > 0. && "Null tangent");
 
  go = &mesh->xpoint[p0->xp];
 
  /* Decide between the two possible configurations */
  n1 = &go->n1[0];
  n2 = &go->n2[0];
 
  ps1 = ux*n1[0] + uy*n1[1] + uz*n1[2];
  ps2 = ux*n2[0] + uy*n2[1] + uz*n2[2];
 
  if ( fabs(ps2)<fabs(ps1) ) {
    n1 = &go->n2[0];
    dv = m[2];
    dn = m[4];
  }
  else{
    dv = m[1];
    dn = m[3];
  }
 
  /* Check that choosed normal is not null */
  /* Remark: a null second normal along ridge point in mmg3d may be the
   * consequence of the computation of the same normal with opposite sign for n1
   * and n2 at a previously inserted ridge point. This append when the ridge
   * delimits a closed angle an we choose the wrong point normal in bezierCP. */
  assert ( n1[0]*n1[0] + n1[1]*n1[1] + n1[2]*n1[2] > 0. && "Null normal");
 
  u[0] = n1[1]*t[2] - n1[2]*t[1];
  u[1] = n1[2]*t[0] - n1[0]*t[2];
  u[2] = n1[0]*t[1] - n1[1]*t[0];
 
  /* If u = n1 ^ t, matrix of the desired metric in (t,u,n1) =
   * diag(m[0],dv,dn). Now, compute the metric in the canonical basis. */
  r[0][0] = t[0];  r[0][1] = u[0];  r[0][2] = n1[0];
  r[1][0] = t[1];  r[1][1] = u[1];  r[1][2] = n1[1];
  r[2][0] = t[2];  r[2][1] = u[2];  r[2][2] = n1[2];
 
  mr[0] = m[0]*r[0][0]*r[0][0] + dv*r[0][1]*r[0][1] + dn*r[0][2]*r[0][2];
  mr[1] = m[0]*r[0][0]*r[1][0] + dv*r[0][1]*r[1][1] + dn*r[0][2]*r[1][2];
  mr[2] = m[0]*r[0][0]*r[2][0] + dv*r[0][1]*r[2][1] + dn*r[0][2]*r[2][2];
  mr[3] = m[0]*r[1][0]*r[1][0] + dv*r[1][1]*r[1][1] + dn*r[1][2]*r[1][2];
  mr[4] = m[0]*r[1][0]*r[2][0] + dv*r[1][1]*r[2][1] + dn*r[1][2]*r[2][2];
  mr[5] = m[0]*r[2][0]*r[2][0] + dv*r[2][1]*r[2][1] + dn*r[2][2]*r[2][2];
  return 1;
}
 
int MMG5_buildridmetnor(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np0,double nt[3],
                        double mr[6],double r[3][3] ) {
  MMG5_pPoint  p0;
  MMG5_pxPoint go;
  double       ps1,ps2,*n1,*n2,*t,*m,dv,dn,u[3];
  int          ier = 0;
 
  p0 = &mesh->point[np0];
  if ( !(MG_GEO & p0->tag) )  return 0;
  m = &met->m[6*np0];
  t = &p0->n[0];
  go = &mesh->xpoint[p0->xp];
 
  /* Decide between the two possible configurations */
  n1 = &go->n1[0];
  n2 = &go->n2[0];
 
  ps1 = nt[0]*n1[0] + nt[1]*n1[1] + nt[2]*n1[2];
  ps2 = nt[0]*n2[0] + nt[1]*n2[1] + nt[2]*n2[2];
 
  if ( fabs(ps2) > fabs(ps1) ) {
    n1 = &go->n2[0];
    dv = m[2];
    dn = m[4];
    ier = 2;
  }
  else{
    dv = m[1];
    dn = m[3];
    ier = 1;
  }
 
  u[0] = n1[1]*t[2] - n1[2]*t[1];
  u[1] = n1[2]*t[0] - n1[0]*t[2];
  u[2] = n1[0]*t[1] - n1[1]*t[0];
 
  /* If u = n1 ^ t, matrix of the desired metric in (t,u,n1) = diag(m[0],dv,0).
     Now, compute the metric in the canonical basis.*/
  r[0][0] = t[0];  r[0][1] = u[0];  r[0][2] = n1[0];
  r[1][0] = t[1];  r[1][1] = u[1];  r[1][2] = n1[1];
  r[2][0] = t[2];  r[2][1] = u[2];  r[2][2] = n1[2];
 
  mr[0] = m[0]*r[0][0]*r[0][0] + dv*r[0][1]*r[0][1] + dn*r[0][2]*r[0][2];
  mr[1] = m[0]*r[0][0]*r[1][0] + dv*r[0][1]*r[1][1] + dn*r[0][2]*r[1][2];
  mr[2] = m[0]*r[0][0]*r[2][0] + dv*r[0][1]*r[2][1] + dn*r[0][2]*r[2][2];
  mr[3] = m[0]*r[1][0]*r[1][0] + dv*r[1][1]*r[1][1] + dn*r[1][2]*r[1][2];
  mr[4] = m[0]*r[1][0]*r[2][0] + dv*r[1][1]*r[2][1] + dn*r[1][2]*r[2][2];
  mr[5] = m[0]*r[2][0]*r[2][0] + dv*r[2][1]*r[2][1] + dn*r[2][2]*r[2][2];
 
  return ier;
}
 
int MMG5_intersecmet22(MMG5_pMesh mesh, double *m,double *n,double *mr) {
  double        det,imn[4],lambda[2],vp[2][2],dm[2],dn[2],d0,d1,ip[4];
  double        isqhmin,isqhmax;
  static int8_t mmgWarn0 = 0;
  int           order;
 
  isqhmin  = 1.0 / (mesh->info.hmin*mesh->info.hmin);
  isqhmax  = 1.0 / (mesh->info.hmax*mesh->info.hmax);
 
  /* Compute imn = M^{-1}N */
  det = m[0]*m[2] - m[1]*m[1];
  if ( fabs(det) < MMG5_EPS*MMG5_EPS ) {
    if ( !mmgWarn0 ) {
      fprintf(stderr,"\n  ## Warning: %s: null metric det : %E \n",__func__,det);
      mmgWarn0 = 1;
    }
    return 0;
  }
  det = 1.0 / det;
 
  imn[0] = det * ( m[2]*n[0] - m[1]*n[1]);
  imn[1] = det * ( m[2]*n[1] - m[1]*n[2]);
  imn[2] = det * (-m[1]*n[0] + m[0]*n[1]);
  imn[3] = det * (-m[1]*n[1] + m[0]*n[2]);
 
  /* Find eigenvalues of imn */
  order = MMG5_eigenv2d(0,imn,lambda,vp);
 
  if ( !order ) {
    if ( !mmgWarn0 ) {
      mmgWarn0 = 1;
      fprintf(stderr,"\n  ## Warning: %s: at least 1 failing"
              " simultaneous reduction.\n",__func__);
    }
    return 0;
  }
 
  /* First case : matrices m and n are homothetic : n = lambda0*m */
  if ( order == 2 ) {
    /* Diagonalize m and truncate eigenvalues : trimn, det, etc... are reused */
    if ( fabs(m[1]) < MMG5_EPS ) {
      dm[0]    = m[0];
      dm[1]    = m[2];
      vp[0][0] = 1;
      vp[0][1] = 0;
      vp[1][0] = 0;
      vp[1][1] = 1;
    }
    else {
      MMG5_eigensym(m,dm,vp);
    }
    /* Eigenvalues of the resulting matrix*/
    dn[0] = MG_MAX(dm[0],lambda[0]*dm[0]);
    dn[0] = MG_MIN(isqhmin,MG_MAX(isqhmax,dn[0]));
    dn[1] = MG_MAX(dm[1],lambda[0]*dm[1]);
    dn[1] = MG_MIN(isqhmin,MG_MAX(isqhmax,dn[1]));
 
    /* Intersected metric = P diag(d0,d1){^t}P, P = (vp0, vp1) stored in columns */
    mr[0] = dn[0]*vp[0][0]*vp[0][0] + dn[1]*vp[1][0]*vp[1][0];
    mr[1] = dn[0]*vp[0][0]*vp[0][1] + dn[1]*vp[1][0]*vp[1][1];
    mr[2] = dn[0]*vp[0][1]*vp[0][1] + dn[1]*vp[1][1]*vp[1][1];
 
    return 1;
  }
 
  /* Second case : both eigenvalues of imn are distinct ; theory says qf associated to m and n
     are diagonalizable in basis (vp0, vp1) - the coreduction basis */
  else if( order == 1 ) {
 
    /* Compute diagonal values in simultaneous reduction basis */
    dm[0] = m[0]*vp[0][0]*vp[0][0] + 2.0*m[1]*vp[0][0]*vp[0][1] + m[2]*vp[0][1]*vp[0][1];
    dm[1] = m[0]*vp[1][0]*vp[1][0] + 2.0*m[1]*vp[1][0]*vp[1][1] + m[2]*vp[1][1]*vp[1][1];
    dn[0] = n[0]*vp[0][0]*vp[0][0] + 2.0*n[1]*vp[0][0]*vp[0][1] + n[2]*vp[0][1]*vp[0][1];
    dn[1] = n[0]*vp[1][0]*vp[1][0] + 2.0*n[1]*vp[1][0]*vp[1][1] + n[2]*vp[1][1]*vp[1][1];
 
    /* Diagonal values of the intersected metric */
    d0 = MG_MAX(dm[0],dn[0]);
    d0 = MG_MIN(isqhmin,MG_MAX(d0,isqhmax));
 
    d1 = MG_MAX(dm[1],dn[1]);
    d1 = MG_MIN(isqhmin,MG_MAX(d1,isqhmax));
 
    /* Intersected metric = tP^-1 diag(d0,d1)P^-1, P = (vp0, vp1) stored in columns */
    det = vp[0][0]*vp[1][1] - vp[0][1]*vp[1][0];
    if ( fabs(det) < MMG5_EPS )  return 0;
    det = 1.0 / det;
 
    ip[0] =  vp[1][1]*det;
    ip[1] = -vp[1][0]*det;
    ip[2] = -vp[0][1]*det;
    ip[3] =  vp[0][0]*det;
 
    mr[0] = d0*ip[0]*ip[0] + d1*ip[2]*ip[2];
    mr[1] = d0*ip[0]*ip[1] + d1*ip[2]*ip[3];
    mr[2] = d0*ip[1]*ip[1] + d1*ip[3]*ip[3];
  }
  return 1;
}
 
int MMG5_intersecmet33(MMG5_pMesh mesh, double *m,double *n,double *mr) {
  double  vp[3][3],dm[3],dn[3],d[3],ivp[3][3];
  double  isqhmin,isqhmax;
  int8_t  i;
 
  isqhmin  = 1.0 / (mesh->info.hmin*mesh->info.hmin);
  isqhmax  = 1.0 / (mesh->info.hmax*mesh->info.hmax);
 
  /* Simultaneous reduction */
  if( !MMG5_simred3d(mesh,m,n,dm,dn,vp) )
    return 0;
 
  /* Diagonal values of the intersected metric */
  for( i = 0; i < 3; i++ ) {
    d[i] = MG_MAX(dm[i],dn[i]);
    d[i] = MG_MIN(isqhmin,MG_MAX(d[i],isqhmax));
  }
 
  /* Intersected metric = tP^-1 diag(d0,d1,d2)P^-1, P = (vp0, vp1,vp2) stored in
   * columns */
 
  /* Compute the inverse of the eigenvectors matrix */
  if( !MMG5_invmat33(vp,ivp) )
    return 0;
 
  /* Recompose matrix */
  MMG5_simredmat(3,mr,d,(double *)ivp);
 
  return 1;
}
 
int MMG5_test_intersecmet22(MMG5_pMesh mesh) {
  double m[3] = { 508., -504,  502.}; /* Test matrix 1 */
  double n[3] = {4020.,-2020.,1020.}; /* Test matrix 2 */
  double intex[3] = {4500.,-2500.,1500.}; /* Exact intersection */
  double intnum[3],maxerr; /* Numerical quantities */
 
  if( !MMG5_intersecmet22(mesh,m,n,intnum) )
    return 0;
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(3,intex,intnum);
  if( maxerr > 1000.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error metric intersection: in function %s, line %d, max error %e\n",
      __func__,__LINE__,maxerr);
    return 0;
  }
 
  for( int8_t it = 0; it < 20; it++ ) {
 
    if( !MMG5_intersecmet22(mesh,n,intnum,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(3,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet22(mesh,intnum,n,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(3,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet22(mesh,m,intnum,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(3,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet22(mesh,intnum,m,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(3,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
  }
 
  return 1;
}
 
int MMG5_test_intersecmet33(MMG5_pMesh mesh) {
  double m[6] = {111./2.,-109./2.,  89./2.,111./2.,-91./2.,111./2.}; /* Test matrix 1 */
  double n[6] = {409./2.,-393./2.,-407./2.,409./2.,391./2.,409./2.}; /* Test matrix 2 */
  double intex[6] = {254.,-246.,-154.,254.,146.,254.}; /* Exact intersection */
  double intnum[6],maxerr; /* Numerical quantities */
 
  if( !MMG5_intersecmet33(mesh,m,n,intnum) )
    return 0;
 
  /* Check error in norm inf */
  maxerr = MMG5_test_mat_error(6,intex,intnum);
  if( maxerr > 1000.*MMG5_EPSOK ) {
    fprintf(stderr,"  ## Error metric intersection: in function %s, line %d, max error %e\n",
      __func__,__LINE__,maxerr);
    return 0;
  }
 
  for( int8_t it = 0; it < 20; it++ ) {
 
    if( !MMG5_intersecmet33(mesh,n,intnum,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(6,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet33(mesh,intnum,n,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(6,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet33(mesh,m,intnum,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(6,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
 
    if( !MMG5_intersecmet33(mesh,intnum,m,intnum) )
      return 0;
 
    /* Check error in norm inf */
    maxerr = MMG5_test_mat_error(6,intex,intnum);
    if( maxerr > 1.e5*MMG5_EPSOK ) {
      fprintf(stderr,"  ## Error metric re-intersection: in function %s, line %d, iteration %d, max error %e\n",
        __func__,__LINE__,it,maxerr);
      return 0;
    }
 
  }
 
  return 1;
}
 
int MMG5_mmgIntextmet(MMG5_pMesh mesh,MMG5_pSol met,MMG5_int np,double me[6],
                       double n[3]) {
  MMG5_pPoint         p0;
  MMG5_pxPoint        go;
  double              hu,isqhmin,isqhmax,dd,alpha1,alpha2,alpha3,u[3];
  double              lambda[3],vp[3][3];
  double              *m,*n1,*n2,*t,r[3][3],mrot[6],mr[3],mtan[3],metan[3];
  int                 order;
  int8_t              i;
  static int8_t       mmgWarn=0, mmgWarn1=0, mmgWarn2=0;
 
  isqhmin = 1.0 / (mesh->info.hmin*mesh->info.hmin);
  isqhmax = 1.0 / (mesh->info.hmax*mesh->info.hmax);
 
  p0 = &mesh->point[np];
  m  = &met->m[6*np];
 
  if ( MG_SIN(p0->tag) || (p0->tag & MG_NOM) ) {
 
    order = MMG5_eigenv3d(1,me,lambda,vp);
    if ( !order ) {
      if ( !mmgWarn ) {
        fprintf(stderr,"\n  ## Warning: %s: Unable to diagonalize at least"
                " 1 metric.\n",__func__);
        mmgWarn = 1;
      }
      return 0;
    }
 
    hu = 0.0;
    for( i = 0; i < 3; i++ ) {
      hu = MG_MAX(hu,lambda[i]);
    }
 
    hu = MG_MIN(isqhmin,hu);
    hu = MG_MAX(isqhmax,hu);
 
    if( hu > m[0] )
      m[0] = m[3] = m[5] = hu;
 
  }
  else if ( p0->tag & MG_GEO ) {
    /* Size prescribed by metric me in direction t */
    t = n;
    hu = me[0]*t[0]*t[0] + me[3]*t[1]*t[1] + me[5]*t[2]*t[2] \
      + 2.0*me[1]*t[0]*t[1] + 2.0*me[2]*t[0]*t[2] + 2.0*me[4]*t[1]*t[2];
 
    hu = MG_MIN(isqhmin,hu);
    hu = MG_MAX(isqhmax,hu);
    m[0] = MG_MAX(m[0],hu);
 
    assert ( p0->xp );
    go = &mesh->xpoint[p0->xp];
    n1 = &go->n1[0];
    n2 = &go->n2[0];
 
    u[0] = n1[1]*t[2] - n1[2]*t[1];
    u[1] = n1[2]*t[0] - n1[0]*t[2];
    u[2] = n1[0]*t[1] - n1[1]*t[0];
    dd = u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
    if ( dd > MMG5_EPSD ) {
      dd = 1.0 / sqrt(dd);
 
      u[0] *= dd;
      u[1] *= dd;
      u[2] *= dd;
 
      hu = me[0]*u[0]*u[0] + me[3]*u[1]*u[1] + me[5]*u[2]*u[2]          \
        + 2.0*me[1]*u[0]*u[1] + 2.0*me[2]*u[0]*u[2] + 2.0*me[4]*u[1]*u[2];
 
      hu = MG_MIN(isqhmin,hu);
      hu = MG_MAX(isqhmax,hu);
      m[1] = MG_MAX(m[1],hu);
    }
    u[0] = n2[1]*t[2] - n2[2]*t[1];
    u[1] = n2[2]*t[0] - n2[0]*t[2];
    u[2] = n2[0]*t[1] - n2[1]*t[0];
    dd = u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
    if ( dd > MMG5_EPSD ) {
      dd = 1.0 / sqrt(dd);
 
      u[0] *= dd;
      u[1] *= dd;
      u[2] *= dd;
 
      hu =     me[0]*u[0]*u[0] +     me[3]*u[1]*u[1] +     me[5]*u[2]*u[2] \
        + 2.0*me[1]*u[0]*u[1] + 2.0*me[2]*u[0]*u[2] + 2.0*me[4]*u[1]*u[2];
 
      hu = MG_MIN(isqhmin,hu);
      hu = MG_MAX(isqhmax,hu);
      m[2] = MG_MAX(m[2],hu);
    }
 
    hu = me[0]*n1[0]*n1[0] + me[3]*n1[1]*n1[1] + me[5]*n1[2]*n1[2]
      + 2.0*me[1]*n1[0]*n1[1] + 2.0*me[2]*n1[0]*n1[2] + 2.0*me[4]*n1[1]*n1[2];
 
    hu = MG_MIN(isqhmin,hu);
    hu = MG_MAX(isqhmax,hu);
    m[3] = hu;
 
    hu = me[0]*n2[0]*n2[0] + me[3]*n2[1]*n2[1] + me[5]*n2[2]*n2[2]
      + 2.0*me[1]*n2[0]*n2[1] + 2.0*me[2]*n2[0]*n2[2] + 2.0*me[4]*n2[1]*n2[2];
 
    hu = MG_MIN(isqhmin,hu);
    hu = MG_MAX(isqhmax,hu);
    m[4] = hu;
  }
  else {
    MMG5_rotmatrix(n,r);
 
    MMG5_rmtr(r,m,mrot);
    mtan[0] = mrot[0];
    mtan[1] = mrot[1];
    mtan[2] = mrot[3];
 
 
    MMG5_rmtr(r,me,mrot);
    metan[0] = mrot[0];
    metan[1] = mrot[1];
    metan[2] = mrot[3];
 
    if ( !MMG5_intersecmet22(mesh,mtan,metan,mr) ) {
      if ( !mmgWarn1 ) {
        fprintf(stderr,"\n  ## Warning: %s: impossible metric intersection:"
                " surfacic metric skipped.\n",__func__);
        mmgWarn1 = 1;
      }
      m[0] = me[0];
      m[1] = me[1];
      m[2] = me[2];
      m[3] = me[3];
      m[4] = me[4];
      m[5] = me[5];
 
      return 0;
    }
 
    mtan[0]  = mr[0]*r[0][0] + mr[1]*r[1][0];
    mtan[1]  = mr[0]*r[0][1] + mr[1]*r[1][1];
    mtan[2]  = mr[0]*r[0][2] + mr[1]*r[1][2];
    metan[0] = mr[1]*r[0][0] + mr[2]*r[1][0];
    metan[1] = mr[1]*r[0][1] + mr[2]*r[1][1];
    metan[2] = mr[1]*r[0][2] + mr[2]*r[1][2];
 
    alpha1 = r[2][0]*mrot[5];
    alpha2 = r[2][1]*mrot[5];
    alpha3 = r[2][2]*mrot[5];
 
    m[0] = r[0][0] * mtan[0] + r[1][0] * metan[0] + r[2][0]*alpha1;
    m[1] = r[0][0] * mtan[1] + r[1][0] * metan[1] + r[2][0]*alpha2;
    m[2] = r[0][0] * mtan[2] + r[1][0] * metan[2] + r[2][0]*alpha3;
    m[3] = r[0][1] * mtan[1] + r[1][1] * metan[1] + r[2][1]*alpha2;
    m[4] = r[0][1] * mtan[2] + r[1][1] * metan[2] + r[2][1]*alpha3;
    m[5] = r[0][2] * mtan[2] + r[1][2] * metan[2] + r[2][2]*alpha3;
 
    order = MMG5_eigenv3d(1,m,lambda,vp);
    if ( !order ) {
      if ( !mmgWarn ) {
        fprintf(stderr,"\n  ## Warning: %s: Unable to diagonalize at least"
                " 1 metric.\n",__func__);
        mmgWarn = 1;
      }
      return 0;
    }
 
    for (i=0; i<3; i++) {
      if( lambda[i]<=0) {
        if ( !mmgWarn2 ) {
          fprintf(stderr,"\n  ## Warning: %s: at least 1 wrong metric "
                  "(eigenvalues : %e %e %e): surfacic metric skipped.\n",
                  __func__,lambda[0],lambda[1],lambda[2]);
          mmgWarn2 = 1;
        }
        m[0] = me[0];
        m[1] = me[1];
        m[2] = me[2];
        m[3] = me[3];
        m[4] = me[4];
        m[5] = me[5];
        return 0;
      }
      lambda[i]=MG_MIN(isqhmin,lambda[i]);
      lambda[i]=MG_MAX(isqhmax,lambda[i]);
    }
 
    m[0] = vp[0][0]*vp[0][0]*lambda[0] + vp[1][0]*vp[1][0]*lambda[1]
      + vp[2][0]*vp[2][0]*lambda[2];
    m[1] = vp[0][0]*vp[0][1]*lambda[0] + vp[1][0]*vp[1][1]*lambda[1]
      + vp[2][0]*vp[2][1]*lambda[2];
    m[2] = vp[0][0]*vp[0][2]*lambda[0] + vp[1][0]*vp[1][2]*lambda[1]
      + vp[2][0]*vp[2][2]*lambda[2];
    m[3] = vp[0][1]*vp[0][1]*lambda[0] + vp[1][1]*vp[1][1]*lambda[1]
      + vp[2][1]*vp[2][1]*lambda[2];
    m[4] = vp[0][1]*vp[0][2]*lambda[0] + vp[1][1]*vp[1][2]*lambda[1]
      + vp[2][1]*vp[2][2]*lambda[2];
    m[5] = vp[0][2]*vp[0][2]*lambda[0] + vp[1][2]*vp[1][2]*lambda[1]
      + vp[2][2]*vp[2][2]*lambda[2];
  }
 
  return 1;
}
 
int MMG5_paratmet(double c0[3],double n0[3],double m[6],double c1[3],double n1[3],double mt[6]) {
  double  r[3][3],mrot[6],mtan[3],lambda[2],vp[2][2],u[3],ps,ll;
 
  /* Take the induced metric tensor in the tangent plane by change of basis : R * M * {^t}R*/
  if ( !MMG5_rotmatrix(n0,r) )  return 0;
  MMG5_rmtr(r,m,mrot);
  mtan[0] = mrot[0];
  mtan[1] = mrot[1];
  mtan[2] = mrot[3];
 
  /* Take eigenvectors of metric tensor in tangent plane */
  MMG5_eigensym(mtan,lambda,vp);
 
  /* Eigenvector in canonical basis = {t}R*vp[0] */
  u[0] = r[0][0]*vp[0][0] + r[1][0]*vp[0][1];
  u[1] = r[0][1]*vp[0][0] + r[1][1]*vp[0][1];
  u[2] = r[0][2]*vp[0][0] + r[1][2]*vp[0][1];
 
  /* Projection in the tangent plane of c1 */
  ps = u[0]*n1[0] + u[1]*n1[1] + u[2]*n1[2];
  u[0] -= ps*n1[0];
  u[1] -= ps*n1[1];
  u[2] -= ps*n1[2];
  ll = u[0]*u[0] + u[1]*u[1] + u[2]*u[2];
  if ( ll < MMG5_EPSD )  return 0;
  ll = 1.0 / sqrt(ll);
  u[0] *= ll;
  u[1] *= ll;
  u[2] *= ll;
 
  /* And the transported metric is diag(lambda[0], lambda[1], mrot[5]) in basis
   * (u,n1^u,n1) */
  r[0][0] = u[0];
  r[1][0] = u[1];
  r[2][0] = u[2];
 
  r[0][1] = n1[1]*u[2] - n1[2]*u[1];
  r[1][1] = n1[2]*u[0] - n1[0]*u[2];
  r[2][1] = n1[0]*u[1] - n1[1]*u[0];
 
  ll = r[0][1]*r[0][1] + r[1][1]*r[1][1] + r[2][1]*r[2][1];
  if ( ll < MMG5_EPSD )  return 0;
  ll = 1.0 / sqrt(ll);
  r[0][1] *= ll;
  r[1][1] *= ll;
  r[2][1] *= ll;
 
  r[0][2] = n1[0];
  r[1][2] = n1[1];
  r[2][2] = n1[2];
 
  /*mt = R * diag(lambda[0], lambda[1], mrot[5])*{^t}R */
  mt[0] = lambda[0]*r[0][0]*r[0][0] + lambda[1]*r[0][1]*r[0][1]
    + mrot[5]*r[0][2]*r[0][2];
 
  mt[1] = lambda[0]*r[0][0]*r[1][0]
    + lambda[1]*r[0][1]*r[1][1] + mrot[5]*r[0][2]*r[1][2];
 
  mt[2] = lambda[0]*r[0][0]*r[2][0]
    + lambda[1]*r[0][1]*r[2][1] + mrot[5]*r[0][2]*r[2][2];
 
  mt[3] = lambda[0]*r[1][0]*r[1][0] + lambda[1]*r[1][1]*r[1][1]
    + mrot[5]*r[1][2]*r[1][2];
 
  mt[4] = lambda[0]*r[2][0]*r[1][0]
    + lambda[1]*r[2][1]*r[1][1] + mrot[5]*r[2][2]*r[1][2];
 
  mt[5] = lambda[0]*r[2][0]*r[2][0] + lambda[1]*r[2][1]*r[2][1]
    + mrot[5]*r[2][2]*r[2][2];
 
  return 1;
}