Changeset 16074 for issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

Timestamp:

09/04/13 15:39:30 (12 years ago)

Author:

seroussi

Message:

CHG: changd gausspenta for subelement2 migration

File:

• issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

@@ -236 +236 @@
         xDelete<double>(seg_vert_coords);
         xDelete<double>(seg_vert_weights);
+}
+/*}}}*/
+/*FUNCTION GaussPenta::GaussPenta(int index,double r1,double r2,int order) {{{*/
+GaussPenta::GaussPenta(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;
+                        this->coords4=xNew<IssmDouble>(numgauss);
+                        for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
+                }
+        }
+        else{
+                /*Double number of gauss points*/
+                GaussPenta *gauss1    = NULL;
+                GaussPenta *gauss2    = NULL;
+                gauss1=new GaussPenta(0,1,2,order);
+                gauss2=new GaussPenta(0,1,2,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<IssmDouble>(this->numgauss);
+                this->coords2=xNew<IssmDouble>(this->numgauss);
+                this->coords3=xNew<IssmDouble>(this->numgauss);
+                this->coords4=xNew<IssmDouble>(this->numgauss);
+                this->weights=xNew<IssmDouble>(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->coords4[ig]=gauss1->coords4[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->coords4[gauss1->numgauss+ig]=gauss2->coords4[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;
+        coord4=UNDEF;
 }
 /*}}}*/