source: issm/oecreview/Archive/24684-25833/ISSM-25441-25442.diff

../trunk-jpl/src/c/classes/gauss/GaussTria.cpp

@@ -348 +348 @@
         coord3=UNDEF;
 }
 /*}}}*/
+GaussTria::GaussTria(IssmDouble r1,IssmDouble r2,int order){/*{{{*/
+
+        /*
+         *  ^ 
+         *  ------------------
+         * 1|\              |
+         *  |  \            |
+         *  |    \          |
+         *  |      \        |
+         *  |        \      |
+         *  |          \    |
+         *  |    +(x,y)  \  |
+         *  |              \|
+         *  +---------------+-->
+         *  0               1
+         *
+         */
+        int         ig;
+        IssmDouble x,y;
+        IssmDouble xy_list[3][2];
+
+        /*Double number of gauss points*/
+        GaussTria *gauss1    = NULL; //blue
+        GaussTria *gauss2    = NULL; //green
+        GaussTria *gauss3    = NULL; //red
+        gauss1=new GaussTria(order); 
+        gauss2=new GaussTria(order); 
+        gauss3=new GaussTria(order); 
+
+        this->numgauss = gauss1->numgauss + gauss2->numgauss + gauss3->numgauss;
+        this->coords1=xNew<IssmDouble>(this->numgauss);
+        this->coords2=xNew<IssmDouble>(this->numgauss);
+        this->coords3=xNew<IssmDouble>(this->numgauss);
+        this->weights=xNew<IssmDouble>(this->numgauss);
+
+        for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points (BLUE)
+                this->coords1[ig]=gauss1->coords1[ig];
+                this->coords2[ig]=gauss1->coords2[ig];
+                this->coords3[ig]=gauss1->coords3[ig];
+                this->weights[ig]=gauss1->weights[ig]*r1*r2;
+        }
+        for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points (GREEN)
+                this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
+                this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
+                this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
+                this->weights[gauss1->numgauss+ig]=gauss2->weights[ig]*r1*(1-r2);
+        }
+        for(ig=0;ig<gauss3->numgauss;ig++){ // Add the second triangle gauss points (RED)
+                this->coords1[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords1[ig];
+                this->coords2[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords2[ig];
+                this->coords3[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->coords3[ig];
+                this->weights[gauss1->numgauss+gauss2->numgauss+ig]=gauss3->weights[ig]*(1-r1);
+        }
+
+        /*Delete gauss points*/
+        delete gauss1;
+        delete gauss2;
+        delete gauss3;
+
+        /*Initialize static fields as undefined*/
+        ig    = -1;
+        weight=UNDEF;
+        coord1=UNDEF;
+        coord2=UNDEF;
+        coord3=UNDEF;
+}
+/*}}}*/
 GaussTria::~GaussTria(){/*{{{*/
         xDelete<IssmDouble>(weights);
         xDelete<IssmDouble>(coords3);

../trunk-jpl/src/c/classes/gauss/GaussTria.h

@@ -32 +32 @@
                 GaussTria(int index1,int index2,int order);
                 GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
                 GaussTria(int index,IssmDouble r1, IssmDouble r2,int order);
+                GaussTria(IssmDouble r1, IssmDouble r2,int order);
                 GaussTria(IssmDouble area_coordinates[2][3],int order);
                 ~GaussTria();
 

../trunk-jpl/src/c/classes/Elements/Tetra.h

@@ -113 +113 @@
       Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
       Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
       Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
+      Gauss*      NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
                 Gauss*      NewGaussBase(int order);
                 Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
                 Gauss*      NewGaussTop(int order);

../trunk-jpl/src/c/classes/Elements/Element.h

@@ -302 +302 @@
                 virtual Gauss*     NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0;
                 virtual Gauss*     NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order)=0;
                 virtual Gauss*     NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order)=0;
+      virtual Gauss*     NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order)=0;
                 virtual Gauss*     NewGaussBase(int order)=0;
                 virtual Gauss*     NewGaussLine(int vertex1,int vertex2,int order)=0;
                 virtual Gauss*     NewGaussTop(int order)=0;

../trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -3548 +3548 @@
         return new GaussTria(point1,fraction1,fraction2,order);
 }
 /*}}}*/
+Gauss*     Tria::NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order){/*{{{*/
+
+        return new GaussTria(fraction1,fraction2,order);
+}
+/*}}}*/
 Gauss*     Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
 
         IssmDouble  area_coordinates[2][3];

../trunk-jpl/src/c/classes/Elements/Tria.h

@@ -207 +207 @@
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
       Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
       Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order);
+      Gauss*         NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order);
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
                 Gauss*         NewGaussBase(int order);
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};

../trunk-jpl/src/c/classes/Elements/Penta.h

@@ -135 +135 @@
                 Gauss*         NewGaussBase(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz);
                 Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
                 Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
+                Gauss*         NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
                 Gauss*         NewGaussBase(int order);
                 Gauss*         NewGaussLine(int vertex1,int vertex2,int order);
                 Gauss*         NewGaussTop(int order);

../trunk-jpl/src/c/classes/Elements/Seg.h

@@ -107 +107 @@
       Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
       Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
       Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
+      Gauss*      NewGauss(IssmDouble fraction1,IssmDouble fraction2,int order){_error_("not implemented yet");};
                 Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
                 Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
                 Gauss*      NewGaussTop(int order){_error_("not implemented yet");};