Index: sm/trunk/src/c/Bamgx/Metric.cpp
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.cpp	(revision 2836)
+++ 	(revision )
@@ -1,302 +1,0 @@
-// -*- Mode : c++ -*-
-//
-// SUMMARY  :      
-// USAGE    :        
-// ORG      : 
-// AUTHOR   : Frederic Hecht
-// E-MAIL   : hecht@ann.jussieu.fr
-//
-
-/*
- 
- This file is part of Freefem++
- 
- Freefem++ 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 2.1 of the License, or
- (at your option) any later version.
- 
- Freefem++  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
- along with Freefem++; if not, write to the Free Software
- Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
- */
-
-#include <stdio.h>
-#include "Mesh2.h"
-
-namespace bamg {
-
-SaveMetricInterpole  LastMetricInterpole;
-
-void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) {
-  double a11=M1.a11,a21=M1.a21,a22=M1.a22;
-  double b11=M2.a11,b21=M2.a21,b22=M2.a22;
-  //  M1 v = l M2 v
-  // (M1 - l M2) v =0
-  // det (M1 - l M2) =0
-  // det (M1 - l M2) = a l^2 + b l + c;
-  // = (a11 - l * b11) * (a22 - l * b22) - (a21 - l * b21 ) ^2
-  //  const double eps = 1.e-5;
-  const double /*c11 = a11*a11,*/ c21= a21*a21;
-  const double /*d11 = b11*b11,*/ d21= b21*b21;
-  const double a=b11*b22 - d21;
-  const double b=-a11*b22-a22*b11+2*a21*b21;
-  const double c=-c21+a11*a22;
-  const double bb = b*b,ac= a*c;
-  const double delta = bb - 4 * ac;
-  //  const double kk=c11+c22+c21+d11+d21+d22;
-  // modif F Hecht feb 1998 
-  // cerr.precision(14);
-  //cerr  <<  bb << " " << ac << " " <<  bb <<  " " <<a << endl;
-  // cerr << a11 << " " << a21 << " " << a22 << endl;
-  //cerr << b11 << " " << b21 << " " << b22 << endl;
-  if (bb + Abs(ac) < 1.0e-20 || (delta< 1.0E-4 * bb ) )
-   {
-   // racine double;
-     // cerr << "double " << endl ;
-    if (Abs(a) < 1.e-30 )
-     l1 = l2 = 0;
-    else 
-     l1=l2=-b/(2*a); 
-    V= D2xD2(1,0,0,1);
-   }
-  else {
-    // cerr << "  -- " << a << endl ;
-    const double delta2 = sqrt(delta);
-    l1= (-b - delta2)/(2*a);
-    l2= (-b + delta2)/(2*a);
-    // M1 v = l M2 v
-    //  ( (M1 - I M2) x,y)  = (x,(M1 - I M2) y) \forall y
-    // so Ker((M1 - I M2)) = Im((M1 - I M2))^\perp
-      double v0 = a11-l1*b11, v1 = a21-l1*b21,v2 = a22 - l1*b22;
-      double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-      double vp1x,vp1y,vp2x,vp2y;
-
-      if(s1 < s0)
-	s0=sqrt(s0),vp1x=v1/s0,vp1y=-v0/s0;
-      else
-	s1=sqrt(s1),vp1x=v2/s1,vp1y=-v1/s1;
-
-      v0 = a11-l2*b11, v1 = a21-l2*b21,v2 = a22 - l2*b22;
-      s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
-      if(s1 < s0)
-	s0=sqrt(s0),vp2x=v1/s0,vp2y=-v0/s0;
-      else
-	s1=sqrt(s1),vp2x=v2/s1,vp2y=-v1/s1;
-		V=D2xD2(vp1x,vp2x,vp1y,vp2y);
-  }
-  return;
-
-}
-
-MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
-MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) {
-      D2xD2 M;
-      double l1,l2;
-      ReductionSimultanee(M1,M2,l1,l2,M);
-      R2 v0(M.x.x,M.y.x);
-      R2 v1(M.x.y,M.y.y);
-      D2xD2 M_1(M.inv());
-      D2xD2 D(Max(M1(v0,v0),M2(v0,v0)),0,0,Max(M1(v1,v1),M2(v1,v1)));
-      D2xD2 Mi(M_1.t()*D*M_1);
-      return MetricAnIso(Mi.x.x,0.5*(Mi.x.y+Mi.y.x),Mi.y.y);
-}
-
-MetricAnIso::MetricAnIso(const Real8  a[3],const  MetricAnIso m0, const  MetricAnIso m1,const  MetricAnIso m2 ){
-  MetricAnIso mab(a[0]*m0.a11 + a[1]*m1.a11 + a[2]*m2.a11,
-		  a[0]*m0.a21 + a[1]*m1.a21 + a[2]*m2.a21,
-		  a[0]*m0.a22 + a[1]*m1.a22 + a[2]*m2.a22);
-  
-  MatVVP2x2 vab(mab);
- 
-  R2 v1(vab.v.x,vab.v.y);
-  R2 v2(-v1.y,v1.x);
-  
-  Real8 h1 = a[0] / m0(v1) + a[1] / m1(v1) + a[2] / m2(v1);
-  Real8 h2 = a[0] / m0(v2) + a[1] / m1(v2) + a[2] / m2(v2);
-
-  vab.lambda1 =  1 / (h1*h1);
-  vab.lambda2 =  1 / (h2*h2);
-  *this = vab;
-}
-
- MetricAnIso::MetricAnIso( Real8  a,const  MetricAnIso ma, Real8  b,const  MetricAnIso mb) { 
-  MetricAnIso mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
-  MatVVP2x2 vab(mab);
-  
-  R2 v1(vab.v.x,vab.v.y);
-  R2 v2(-v1.y,v1.x);
-  
-
-  Real8 h1 = a / ma(v1) + b / mb(v1);
-  Real8 h2 = a / ma(v2) + b / mb(v2);
-  vab.lambda1 =  1 / (h1*h1);
-  vab.lambda2 =  1 / (h2*h2);
-  *this = vab;
-}
-
-
-
- MatVVP2x2::MatVVP2x2(const MetricAnIso M){
-  double a11=M.a11,a21=M.a21,a22=M.a22;
-  const double eps = 1.e-5;
-  double c11 = a11*a11, c22 = a22*a22, c21= a21*a21;
-  double b=-a11-a22,c=-c21+a11*a22;
-  double   delta = b*b - 4 * c ;
-  double n2=(c11+c22+c21);
-  if ( n2 < 1e-30) 
-     lambda1=lambda2=0,v.x=1,v.y=0;
-  else if (delta < eps*n2)
-    { 
-      lambda1=lambda2=-b/2, v.x=1,v.y=0;
-    }
-  else 
-    {  //    ---  construction  de 2 vecteur dans (Im ( A - D(i) Id) ortogonal 
-      delta = sqrt(delta);
-      lambda1 = (-b-delta)/2.0,lambda2 = (-b+delta)/2.0;
-      double v0 = a11-lambda1, v1 = a21,v2 = a22 - lambda1;
-      double s0 = v0*v0 + v1*v1, s1 = v1*v1 +v2*v2;
- 
-      if(s1 < s0)
-	s0=sqrt(s0),v.x=v1/s0,v.y=-v0/s0;
-      else
-	s1=sqrt(s1),v.x=v2/s1,v.y=-v1/s1;
-    };
-}
-
-
- int MetricAnIso::IntersectWith(const MetricAnIso M2) 
-{
-  //cerr << " - " << *this << M2 <<  endl;
-      int r=0;
-      MetricAnIso & M1 = *this;
-      D2xD2 M;
-      double l1,l2;
-      
-      ReductionSimultanee(*this,M2,l1,l2,M);
-      // cerr << M << endl;
-      R2 v1(M.x.x,M.y.x);
-      R2 v2(M.x.y,M.y.y);
-      double l11=M1(v1,v1);
-      double l12=M1(v2,v2);
-      double l21=M2(v1,v1);
-      double l22=M2(v2,v2);
-      if ( l11 < l21 )  r=1,l11=l21;
-      if ( l12 < l22 )  r=1,l12=l22; 
-      // cerr << r << endl;
-      if (r) { // change
-        D2xD2 M_1(M.inv());
-        D2xD2 D(l11,0,0,l12); 
-        D2xD2 Mi(M_1.t()*D*M_1);
-        a11=Mi.x.x;
-        a21=0.5*(Mi.x.y+Mi.y.x);
-        a22=Mi.y.y; }
-      return r;
-}
-
-Real8 LengthInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB)
-{
-  Real8 k=1./2.;
-  int level=0;
-  static int kkk=0;
-  static  Metric Ms1[32],Ms2[32];
-  static Real8 lMs1[32],lMs2[32];
-  static double K[32];
-  Real8 l=0,sss=0;
-  Ms1[level]=Ma;
-  Ms2[level]=Mb;
-  Real8 sa =  Ma(AB);
-  Real8 sb =  Mb(AB);
-  lMs1[level]=sa;
-  lMs2[level]=sb;
-  K[level]=k;
-  level++;
-  int i=0;
-  Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
-  Real8  sstop = 0.1; // Max(0.6,(sa+sb)/5000);
-  while (level) {
-    level--;
-    Metric M1=Ms1[level];
-    Metric M2=Ms2[level];
-    k=K[level];
-    Real8 s1=  lMs1[level];
-    Real8 s2=  lMs2[level];
-
-    Real8 s= (s1+s2)*k;
-//    if (level >20  && i < 2030-level)
-//    cout << "                  level " << level << " " << i << " " << s << " " << k <<endl;
-    if( s > sstop   && level < 30 && i < 500-level ) {
-      Metric Mi(0.5,M1,0.5,M2);
-      Real8 si = Mi(AB);
-      if( Abs((s1+s2)-(si+si)) > s1*0.001) 
-	{
-	  k=k/2;
-	  // we begin by the end to walk in the correct sens from a to b
-	       // due to the stack 
-	       Ms1[level]=Mi;
-	  Ms2[level]=M2;
-	  lMs1[level]=si;
-	  lMs2[level]=s2;
-	  K[level]=k;
-	  level++;
-	  Ms1[level]=M1;
-	  Ms2[level]=Mi;
-	  lMs1[level]=s1;
-	  lMs2[level]=si;
-	  K[level]=k;
-	  level++;
-	}
-      else
-	L[i]= l += s,S[i]=sss+=k,i++;
-    }
-    else 
-      L[i]= l += s,S[i]=sss+=k,i++;//cout << i << " l = " << l << " sss = " << sss << endl;
-  }
-  // warning for optimisation S is in [0:0.5] not in [0:1]
-  assert(i<512);
-  LastMetricInterpole.lab=l;
-  LastMetricInterpole.opt=i;
-  if (i>200 && kkk++<10)
-     cout << "Warning LengthInterpole: ( i = " << i << " l = " << l << " sss " << sss << " ) " << sstop <<endl;
-  return l;
-}
-
-Real8 abscisseInterpole(const MetricAnIso Ma,const  MetricAnIso Mb, R2 AB,Real8 s,int optim)
-{ 
-  if(!optim)  LengthInterpole(Ma,Mb,AB);
-  Real8 l  = s* LastMetricInterpole.lab,r;
-  int j=LastMetricInterpole.opt-1,i=0,k;
-  
-  Real8 * L= LastMetricInterpole.L, *S = LastMetricInterpole.S;
-  // warning for optimisation S is the abcisse in [0:0.5]
-  // and L is le lenght 
-  if(l<=L[0])
-    r=2*S[0]*l/L[0];
-  else if (l>=L[j])
-    r=1;
-  else 
-    {
-      while (j-i>1)
-	{
-	  k= (i+j)/2;
-	  if(l<=L[k])
-	    j=k;// l<=L[j] 
-	  else
-	    i=k; //  L[i]<l
-	};
-      //   cout  << i << " " << j  <<" " << L[i] << " " << L[j] << " " << S[i] << " " <<  S[j]  << " l=" << l << endl;
-      if (i==j)
-	r = 2*S[i];
-      else
-	r =  2*(S[i]*(L[j]-l)+ S[j]*(l-L[i]))/(L[j]-L[i]);
-    }
-  assert(r<=1 && r>=0);
-  return r ;
-    
-}
-}   // end of namespace bamg 
Index: /issm/trunk/src/c/Bamgx/Metric.h
===================================================================
--- /issm/trunk/src/c/Bamgx/Metric.h	(revision 2836)
+++ /issm/trunk/src/c/Bamgx/Metric.h	(revision 2837)
@@ -163,4 +163,5 @@
 
 void ReductionSimultanee( MetricAnIso M1,  MetricAnIso M2,double & l1,double & l2, D2xD2 & V) ;
+MetricAnIso Intersection(const MetricAnIso M1,const MetricAnIso M2) ;
 
 inline MetricAnIso::MetricAnIso(const MatVVP2x2 M)  
Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 2836)
+++ /issm/trunk/src/c/Makefile.am	(revision 2837)
@@ -323,9 +323,10 @@
 					./Bamgx/Mesh2.h \
 					./Bamgx/meshtype.h \
+					./Bamgx/objects/MatVVP2x2.cpp \
+					./Bamgx/objects/MetricAnIso.cpp \
 					./Bamgx/objects/GeometricalEdge.cpp \
 					./Bamgx/objects/Triangles.cpp	\
 					./Bamgx/objects/Triangle.cpp	\
 					./Bamgx/objects/Geometry.cpp	\
-					./Bamgx/Metric.cpp \
 					./Bamgx/Metric.h \
 					./Bamgx/QuadTree.cpp \
@@ -663,9 +664,10 @@
 					./Bamgx/Mesh2.h \
 					./Bamgx/meshtype.h \
+					./Bamgx/objects/MatVVP2x2.cpp \
+					./Bamgx/objects/MetricAnIso.cpp \
 					./Bamgx/objects/GeometricalEdge.cpp \
 					./Bamgx/objects/Triangles.cpp	\
 					./Bamgx/objects/Triangle.cpp	\
 					./Bamgx/objects/Geometry.cpp	\
-					./Bamgx/Metric.cpp \
 					./Bamgx/Metric.h \
 					./Bamgx/QuadTree.cpp \