Changeset 5489

Timestamp:

08/23/10 07:54:04 (15 years ago)

Author:

Mathieu Morlighem

Message:

As usual

Location:

issm/trunk/src/c/objects/Bamg

Files:

• GeometricalEdge.cpp (modified) (1 diff)
• Mesh.cpp (modified) (1 diff)
• Metric.cpp (modified) (1 diff)
• Metric.h (modified) (1 diff)

issm/trunk/src/c/objects/Bamg/GeometricalEdge.cpp

@@ -31 +31 @@
                 ISSMASSERT(theta>-1e-12 && theta<1+1e-12);
 
-                if (TgA()) 
-                 if (TgB()) // interpolation d'hermite
-                        { cb =  theta*theta*(3-2*theta);
-                         ca =  1-cb;     
-                         cta = (1-theta)*(1-theta)*theta;
-                         ctb = (theta-1)*theta*theta ;
+                if (TgA()){ 
+                        if (TgB()){ //Hermite interpolation
+                                cb =  theta*theta*(3-2*theta);
+                                ca =  1-cb;     
+                                cta = (1-theta)*(1-theta)*theta;
+                                ctb = (theta-1)*theta*theta ;
                         }
-                 else { // 1-t*t, t-t*t, t*t
-                         double t = theta;
-                         cb = t*t;
-                         ca = 1-cb;
-                         cta= t-cb;
-                         ctb=0;    
-                 }    
-                else
-                 if (TgB()){
-                         double t = 1-theta;
-                         ca = t*t;
-                         cb = 1-ca;
-                         ctb= -t+ca;
-                         cta=0;    
-                 }
-                 else {
-                         ca =(1-theta),cb = theta,cta=ctb=0; // lagrange P1
-                 }
-                return A*ca + B*cb + tg[0]* cta + tg[1] * ctb;
+                        else {
+                                double t = theta;
+                                cb = t*t;
+                                ca = 1-cb;
+                                cta= t-cb;
+                                ctb=0;    
+                        }
+                }
+                else{
+                        if (TgB()){
+                                double t = 1-theta;
+                                ca = t*t;
+                                cb = 1-ca;
+                                ctb= -t+ca;
+                                cta=0;    
+                        }
+                        else { // lagrange P1
+                                ca =(1-theta);
+                                cb = theta;
+                                cta=ctb=0;
+                        }
+                }
+                return A*ca + B*cb + tg[0]*cta + tg[1]*ctb;
           }
         /*}}}1*/

issm/trunk/src/c/objects/Bamg/Mesh.cpp

@@ -5047 +5047 @@
                                                                                         x += xstep;
                                                                                         B =  e->F(k ? x : 1-x);
-                                                                                        MBs= background ? BTh.MetricAt(B) : Metric(1-x, MA, x ,MB);
+                                                                                        MBs= background ? BTh.MetricAt(B) : Metric(1-x,MA,x,MB);
                                                                                         AB = A-B;
-                                                                                        lSubEdge[kk]= (ledge += (MAs(AB)+MBs(AB))/2);
+                                                                                        lSubEdge[kk]= (ledge+=(MAs(AB)+MBs(AB))/2);
                                                                                 }
                                                                         }
 
-                                                                        double lcurveb = lcurve+ ledge ;
-
-                                                                        //Now, create corresponding points
-                                                                        while (lcurve<=s && s <= lcurveb && nbv < nbvend){
+                                                                        double lcurveb = lcurve+ledge;
+
+                                                                        /*Now, create corresponding points*/
+                                                                        while (lcurve<=s && s<=lcurveb && nbv<nbvend){
 
                                                                                 double ss = s-lcurve;
 
                                                                                 /*Find the SubEdge containing ss using Dichotomy*/
-
                                                                                 int kk0=-1,kk1=NbSubEdge-1,kkk;
                                                                                 double ll0=0,ll1=ledge,llk;

issm/trunk/src/c/objects/Bamg/Metric.cpp

@@ -51 +51 @@
                 /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/EigenMetric)*/
 
+                /*Compute metric (linear combination of ma and mb)*/
                 Metric mab(a*ma.a11+b*mb.a11,a*ma.a21+b*mb.a21,a*ma.a22+b*mb.a22);
+
+                /*Get Eigen values and vectors*/
                 EigenMetric vab(mab);
-
                 R2 v1(vab.v.x,vab.v.y);
                 R2 v2(-v1.y,v1.x);
 
-
-                double h1 = a / ma(v1) + b / mb(v1);
-                double h2 = a / ma(v2) + b / mb(v2);
-                vab.lambda1 =  1 / (h1*h1);
-                vab.lambda2 =  1 / (h2*h2);
-                *this = vab;
+                /*Modify eigen values (a+b=1)*/
+                double h1 = a/ma(v1) + b/mb(v1);
+                double h2 = a/ma(v2) + b/mb(v2);
+                vab.lambda1 =  1/(h1*h1);
+                vab.lambda2 =  1/(h2*h2);
+                *this=vab;
         }
         /*}}}1*/

issm/trunk/src/c/objects/Bamg/Metric.h

@@ -46 +46 @@
                         Metric operator/(double c) const {double c2=1/(c*c);return  Metric(a11*c2,a21*c2,a22*c2);} 
                         operator D2xD2(){ return D2xD2(a11,a21,a21,a22);}
-                        double  operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};
+                        double  operator()(R2 x) const { return sqrt(x.x*x.x*a11+2*x.x*x.y*a21+x.y*x.y*a22);};        // length of x in metric sqrt(<Mx,x>)
                         double  operator()(R2 x,R2 y) const { return x.x*y.x*a11+(x.x*x.y+x.y*y.x)*a21+x.y*y.y*a22;};