source: issm/oecreview/Archive/14312-15392/ISSM-14778-14779.diff@ 15393

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

@@ -92 +92 @@
         xDelete<double>(seg_weights);
 }
 /*}}}*/
+/*FUNCTION GaussTria::GaussTria(int index,double r1,double r2,int order) {{{*/
+GaussTria::GaussTria(int index,IssmDouble r1,IssmDouble r2,bool mainlyfloating,int order){
+
+        /*
+         *  ^ 
+         *  |
+         * 1|\
+         *  |  \
+         *  |    \
+         *  |      \
+         *  |        \
+         *  |          \
+         *  |    +(x,y)  \
+         *  |              \
+         *  +---------------+-->
+         *  0               1
+         *
+         */
+        int         ig;
+        IssmDouble x,y;
+        IssmDouble xy_list[3][2];
+
+        if(mainlyfloating){
+                /*Get gauss points*/
+                GaussLegendreTria(&this->numgauss,&this->coords1,&this->coords2,&this->coords3,&this->weights,order);
+
+                xy_list[0][0]=0;  xy_list[0][1]=0; 
+                xy_list[1][0]=r1; xy_list[1][1]=0; 
+                xy_list[2][0]=0;  xy_list[2][1]=r2; 
+
+                for(ig=0;ig<this->numgauss;ig++){
+                        x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
+                        y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
+
+                        switch(index){
+                                case 0:
+                                        this->coords1[ig] = 1.-x-y;
+                                        this->coords2[ig] = x;
+                                        this->coords3[ig] = y;
+                                        break;
+                                case 1:
+                                        this->coords1[ig] = y;
+                                        this->coords2[ig] = 1.-x-y;
+                                        this->coords3[ig] = x;
+                                        break;
+                                case 2:
+                                        this->coords1[ig] = x;
+                                        this->coords2[ig] = y;
+                                        this->coords3[ig] = 1.-x-y;
+                                        break;
+                                default:
+                                        _error_("index "<<index<<" not supported yet");
+                        }
+                        this->weights[ig] = this->weights[ig]*r1*r2;
+                }
+        }
+        else{
+                /*Double number of gauss points*/
+                GaussTria *gauss1    = NULL;
+                GaussTria *gauss2    = NULL;
+                gauss1=new GaussTria(order);
+                gauss2=new GaussTria(order);
+
+                xy_list[0][0]=r1; xy_list[0][1]=0; 
+                xy_list[1][0]=0;  xy_list[1][1]=1.; 
+                xy_list[2][0]=0;  xy_list[2][1]=r2; 
+
+                        //gauss1->Echo();
+                for(ig=0;ig<gauss1->numgauss;ig++){
+                        x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
+                        y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
+
+                        switch(index){
+                                case 0:
+                                        gauss1->coords1[ig] = 1.-x-y;
+                                        gauss1->coords2[ig] = x;
+                                        gauss1->coords3[ig] = y;
+                                        break;
+                                case 1:
+                                        gauss1->coords1[ig] = y;
+                                        gauss1->coords2[ig] = 1.-x-y;
+                                        gauss1->coords3[ig] = x;
+                                        break;
+                                case 2:
+                                        gauss1->coords1[ig] = x;
+                                        gauss1->coords2[ig] = y;
+                                        gauss1->coords3[ig] = 1.-x-y;
+                                        break;
+                                default:
+                                        _error_("index "<<index<<" not supported yet");
+                        }
+                        gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
+                }
+                        //gauss1->Echo();
+                xy_list[0][0]=r1; xy_list[0][1]=0; 
+                xy_list[1][0]=1.; xy_list[1][1]=0; 
+                xy_list[2][0]=0;  xy_list[2][1]=1.; 
+
+                        //gauss2->Echo();
+                for(ig=0;ig<gauss2->numgauss;ig++){
+                        x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
+                        y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
+
+                        switch(index){
+                                case 0:
+                                        gauss2->coords1[ig] = 1.-x-y;
+                                        gauss2->coords2[ig] = x;
+                                        gauss2->coords3[ig] = y;
+                                        break;
+                                case 1:
+                                        gauss2->coords1[ig] = y;
+                                        gauss2->coords2[ig] = 1.-x-y;
+                                        gauss2->coords3[ig] = x;
+                                        break;
+                                case 2:
+                                        gauss2->coords1[ig] = x;
+                                        gauss2->coords2[ig] = y;
+                                        gauss2->coords3[ig] = 1.-x-y;
+                                        break;
+                                default:
+                                        _error_("index "<<index<<" not supported yet");
+                        }
+                        gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
+                }
+
+                this->numgauss = gauss1->numgauss + gauss2->numgauss;
+                this->coords1=xNew<IssmPDouble>(this->numgauss);
+                this->coords2=xNew<IssmPDouble>(this->numgauss);
+                this->coords3=xNew<IssmPDouble>(this->numgauss);
+                this->weights=xNew<IssmPDouble>(this->numgauss);
+
+                for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
+                        this->coords1[ig]=gauss1->coords1[ig];
+                        this->coords2[ig]=gauss1->coords2[ig];
+                        this->coords3[ig]=gauss1->coords3[ig];
+                        this->weights[ig]=gauss1->weights[ig];
+                }
+                for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
+                        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];
+                }
+
+                /*Delete gauss points*/
+                delete gauss1;
+                delete gauss2;
+        }
+
+        /*Initialize static fields as undefined*/
+        weight=UNDEF;
+        coord1=UNDEF;
+        coord2=UNDEF;
+        coord3=UNDEF;
+}
+/*}}}*/
 /*FUNCTION GaussTria::~GaussTria(){{{*/
 GaussTria::~GaussTria(){
         xDelete<IssmPDouble>(weights);

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

@@ -31 +31 @@
                 GaussTria();
                 GaussTria(int order);
                 GaussTria(int index1,int index2,int order);
+                GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
                 ~GaussTria();
 
                 /*Methods*/