Changeset 25436


Ignore:
Timestamp:
08/21/20 08:37:13 (5 years ago)
Author:
Mathieu Morlighem
Message:

NEW: using while loop instead of for, may speed up the code a bit

Location:
issm/trunk-jpl/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r25379 r25436  
    14511451
    14521452        /* Start  looping on the number of gaussian points: */
    1453         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1454                 gauss->GaussPoint(ig);
     1453        while(gauss->next()){
    14551454
    14561455                friction->GetAlpha2(&alpha2,gauss);
    … …  
    15321531        /* Start  looping on the number of gaussian points: */
    15331532        Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3);
    1534         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1535                 gauss->GaussPoint(ig);
     1533        while(gauss->next()){
    15361534
    15371535                this->GetBSSAFriction(B,element,dim,xyz_list,gauss);
    … …  
    15991597        /* Start  looping on the number of gaussian points: */
    16001598        Gauss* gauss = element->NewGauss(2);
    1601         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1602                 gauss->GaussPoint(ig);
     1599        while(gauss->next()){
    16031600
    16041601                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    17131710        /* Start  looping on the number of gaussian points: */
    17141711        Gauss* gauss=element->NewGauss(2);
    1715         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1716 
    1717                 gauss->GaussPoint(ig);
     1712        while(gauss->next()){
    17181713
    17191714                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    … …  
    17831778        /*Start looping on Gaussian points*/
    17841779        Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3);
    1785         for(int ig=gauss->begin();ig<gauss->end();ig++){
    1786 
    1787                 gauss->GaussPoint(ig);
     1780        while(gauss->next()){
    17881781                thickness_input->GetInputValue(&thickness,gauss);
    17891782                sealevel_input->GetInputValue(&sealevel,gauss);

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

    r25415 r25436  
    230230                        break;
    231231                case P1Enum:
     232                        for(int i=0;i<NUMVERTICES;i++){
     233                                if(xIsNan<IssmDouble>(values[i])) values[i] = 0.;
     234                                if(xIsNan<IssmDouble>(values_min[i])) values_min[i] = 0.;
     235                                if(xIsNan<IssmDouble>(values_max[i])) values_max[i] = 0.;
     236                        }
    232237                        control_input->SetControl(interpolation_enum,NUMVERTICES,&vertexlids[0],values,values_min,values_max);
    233238                        break;

  • issm/trunk-jpl/src/c/classes/gauss/Gauss.h

    r20827 r25436  
    1111                IssmDouble   weight;
    1212
    13                 virtual        ~Gauss(){};
     13                virtual      ~Gauss(){};
    1414                virtual int  begin(void)=0;
    1515                virtual void Echo(void)=0;
    1616                virtual int  end(void)=0;
     17                virtual bool next(void)=0;
    1718                virtual int  Enum(void)=0;
    1819                virtual void GaussNode(int finitelement,int iv)=0;

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

    r24119 r25436  
    1616GaussPenta::GaussPenta(){/*{{{*/
    1717
     18        ig = -1;
    1819        numgauss=-1;
    1920
    … …  
    3435
    3536        /*Intermediaries*/
    36         int     ighoriz,igvert;
    37         int     numgauss_horiz;
    38         int     numgauss_vert;
    39         IssmDouble *coords1_horiz = NULL;
    40         IssmDouble *coords2_horiz = NULL;
    41         IssmDouble *coords3_horiz = NULL;
     37        int         numgauss_horiz;
     38        int         numgauss_vert;
     39        IssmDouble *coords1_horiz  = NULL;
     40        IssmDouble *coords2_horiz  = NULL;
     41        IssmDouble *coords3_horiz  = NULL;
    4242        IssmDouble *weights_horiz  = NULL;
    43         double *coords_vert = NULL;
    44         double *weights_vert   = NULL;
     43        double     *coords_vert    = NULL;
     44        double     *weights_vert   = NULL;
    4545
    4646        /*Get gauss points*/
    … …  
    5050
    5151        /*Allocate GaussPenta fields*/
     52   ig      = -1;
    5253        numgauss=numgauss_horiz*numgauss_vert;
    5354        coords1=xNew<IssmDouble>(numgauss);
    … …  
    5859
    5960        /*Combine Horizontal and vertical points*/
    60         for (ighoriz=0; ighoriz<numgauss_horiz; ighoriz++){
    61                 for (igvert=0; igvert<numgauss_vert; igvert++){
     61        for(int ighoriz=0; ighoriz<numgauss_horiz; ighoriz++){
     62                for(int igvert=0; igvert<numgauss_vert; igvert++){
    6263                        coords1[numgauss_vert*ighoriz+igvert]=coords1_horiz[ighoriz];
    6364                        coords2[numgauss_vert*ighoriz+igvert]=coords2_horiz[ighoriz];
    … …  
    9293
    9394        /*Get Segment gauss points*/
     95   ig      = -1;
    9496        numgauss=order;
    9597        GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
    … …  
    166168        }
    167169
     170   this->ig = -1;
     171
    168172}
    169173/*}}}*/
    … …  
    182186
    183187        /*Allocate GaussPenta fields*/
     188   this->ig = -1;
    184189        numgauss=order_horiz*order_vert;
    185190        coords1=xNew<IssmDouble>(numgauss);
    … …  
    236241GaussPenta::GaussPenta(int index,IssmDouble r1,IssmDouble r2,bool mainlyfloating,int order){/*{{{*/
    237242
    238         /*
    239          *  ^
    240          *  |
    241          * 1|\
    242          *  |  \
    243          *  |    \
    244          *  |      \
    245          *  |        \
    246          *  |          \
    247          *  |    +(x,y)  \
    248          *  |              \
    249          *  +---------------+-->
    250          *  0               1
    251          *
    252          */
    253         int         ig;
    254243        IssmDouble x,y;
    255244        IssmDouble xy_list[3][2];
    … …  
    263252                xy_list[2][0]=0;  xy_list[2][1]=r2;
    264253
    265                 for(ig=0;ig<this->numgauss;ig++){
    266                         x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
    267                         y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
     254                for(int ii=0;ii<this->numgauss;ii++){
     255                        x = this->coords1[ii]*xy_list[0][0] + this->coords2[ii]*xy_list[1][0] + this->coords3[ii]*xy_list[2][0];
     256                        y = this->coords1[ii]*xy_list[0][1] + this->coords2[ii]*xy_list[1][1] + this->coords3[ii]*xy_list[2][1];
    268257
    269258                        switch(index){
    270259                                case 0:
    271                                         this->coords1[ig] = 1.-x-y;
    272                                         this->coords2[ig] = x;
    273                                         this->coords3[ig] = y;
     260                                        this->coords1[ii] = 1.-x-y;
     261                                        this->coords2[ii] = x;
     262                                        this->coords3[ii] = y;
    274263                                        break;
    275264                                case 1:
    276                                         this->coords1[ig] = y;
    277                                         this->coords2[ig] = 1.-x-y;
    278                                         this->coords3[ig] = x;
     265                                        this->coords1[ii] = y;
     266                                        this->coords2[ii] = 1.-x-y;
     267                                        this->coords3[ii] = x;
    279268                                        break;
    280269                                case 2:
    281                                         this->coords1[ig] = x;
    282                                         this->coords2[ig] = y;
    283                                         this->coords3[ig] = 1.-x-y;
     270                                        this->coords1[ii] = x;
     271                                        this->coords2[ii] = y;
     272                                        this->coords3[ii] = 1.-x-y;
    284273                                        break;
    285274                                default:
    286275                                        _error_("index "<<index<<" not supported yet");
    287276                        }
    288                         this->weights[ig] = this->weights[ig]*r1*r2;
     277                        this->weights[ii] = this->weights[ii]*r1*r2;
    289278                }
    290279                this->coords4=xNew<IssmDouble>(numgauss);
    291                 for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
     280                for(int ii=0;ii<numgauss;ii++) this->coords4[ii]=-1.0;
    292281        }
    293282        else{
    … …  
    303292
    304293                        //gauss1->Echo();
    305                 for(ig=0;ig<gauss1->numgauss;ig++){
    306                         x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
    307                         y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
     294                for(int ii=0;ii<gauss1->numgauss;ii++){
     295                        x = gauss1->coords1[ii]*xy_list[0][0] + gauss1->coords2[ii]*xy_list[1][0] + gauss1->coords3[ii]*xy_list[2][0];
     296                        y = gauss1->coords1[ii]*xy_list[0][1] + gauss1->coords2[ii]*xy_list[1][1] + gauss1->coords3[ii]*xy_list[2][1];
    308297
    309298                        switch(index){
    310299                                case 0:
    311                                         gauss1->coords1[ig] = 1.-x-y;
    312                                         gauss1->coords2[ig] = x;
    313                                         gauss1->coords3[ig] = y;
     300                                        gauss1->coords1[ii] = 1.-x-y;
     301                                        gauss1->coords2[ii] = x;
     302                                        gauss1->coords3[ii] = y;
    314303                                        break;
    315304                                case 1:
    316                                         gauss1->coords1[ig] = y;
    317                                         gauss1->coords2[ig] = 1.-x-y;
    318                                         gauss1->coords3[ig] = x;
     305                                        gauss1->coords1[ii] = y;
     306                                        gauss1->coords2[ii] = 1.-x-y;
     307                                        gauss1->coords3[ii] = x;
    319308                                        break;
    320309                                case 2:
    321                                         gauss1->coords1[ig] = x;
    322                                         gauss1->coords2[ig] = y;
    323                                         gauss1->coords3[ig] = 1.-x-y;
     310                                        gauss1->coords1[ii] = x;
     311                                        gauss1->coords2[ii] = y;
     312                                        gauss1->coords3[ii] = 1.-x-y;
    324313                                        break;
    325314                                default:
    326315                                        _error_("index "<<index<<" not supported yet");
    327316                        }
    328                         gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
     317                        gauss1->weights[ii] = gauss1->weights[ii]*r1*(1-r2);
    329318                }
    330319                        //gauss1->Echo();
    … …  
    334323
    335324                        //gauss2->Echo();
    336                 for(ig=0;ig<gauss2->numgauss;ig++){
    337                         x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
    338                         y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
     325                for(int ii=0;ii<gauss2->numgauss;ii++){
     326                        x = gauss2->coords1[ii]*xy_list[0][0] + gauss2->coords2[ii]*xy_list[1][0] + gauss2->coords3[ii]*xy_list[2][0];
     327                        y = gauss2->coords1[ii]*xy_list[0][1] + gauss2->coords2[ii]*xy_list[1][1] + gauss2->coords3[ii]*xy_list[2][1];
    339328
    340329                        switch(index){
    341330                                case 0:
    342                                         gauss2->coords1[ig] = 1.-x-y;
    343                                         gauss2->coords2[ig] = x;
    344                                         gauss2->coords3[ig] = y;
     331                                        gauss2->coords1[ii] = 1.-x-y;
     332                                        gauss2->coords2[ii] = x;
     333                                        gauss2->coords3[ii] = y;
    345334                                        break;
    346335                                case 1:
    347                                         gauss2->coords1[ig] = y;
    348                                         gauss2->coords2[ig] = 1.-x-y;
    349                                         gauss2->coords3[ig] = x;
     336                                        gauss2->coords1[ii] = y;
     337                                        gauss2->coords2[ii] = 1.-x-y;
     338                                        gauss2->coords3[ii] = x;
    350339                                        break;
    351340                                case 2:
    352                                         gauss2->coords1[ig] = x;
    353                                         gauss2->coords2[ig] = y;
    354                                         gauss2->coords3[ig] = 1.-x-y;
     341                                        gauss2->coords1[ii] = x;
     342                                        gauss2->coords2[ii] = y;
     343                                        gauss2->coords3[ii] = 1.-x-y;
    355344                                        break;
    356345                                default:
    357346                                        _error_("index "<<index<<" not supported yet");
    358347                        }
    359                         gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
     348                        gauss2->weights[ii] = gauss2->weights[ii]*(1-r1);
    360349                }
    361350
    … …  
    367356                this->weights=xNew<IssmDouble>(this->numgauss);
    368357
    369                 for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
    370                         this->coords1[ig]=gauss1->coords1[ig];
    371                         this->coords2[ig]=gauss1->coords2[ig];
    372                         this->coords3[ig]=gauss1->coords3[ig];
    373                         this->coords4[ig]=gauss1->coords4[ig];
    374                         this->weights[ig]=gauss1->weights[ig];
    375                 }
    376                 for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
    377                         this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
    378                         this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
    379                         this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
    380                         this->coords4[gauss1->numgauss+ig]=gauss2->coords4[ig];
    381                         this->weights[gauss1->numgauss+ig]=gauss2->weights[ig];
     358                for(int ii=0;ii<gauss1->numgauss;ii++){ // Add the first triangle gauss points
     359                        this->coords1[ii]=gauss1->coords1[ii];
     360                        this->coords2[ii]=gauss1->coords2[ii];
     361                        this->coords3[ii]=gauss1->coords3[ii];
     362                        this->coords4[ii]=gauss1->coords4[ii];
     363                        this->weights[ii]=gauss1->weights[ii];
     364                }
     365                for(int ii=0;ii<gauss2->numgauss;ii++){ // Add the second triangle gauss points
     366                        this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii];
     367                        this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii];
     368                        this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii];
     369                        this->coords4[gauss1->numgauss+ii]=gauss2->coords4[ii];
     370                        this->weights[gauss1->numgauss+ii]=gauss2->weights[ii];
    382371                }
    383372
    … …  
    393382        coord3=UNDEF;
    394383        coord4=UNDEF;
     384   ig    = -1;
    395385}
    396386/*}}}*/
    … …  
    408398
    409399        /*Allocate GaussPenta fields*/
     400   ig      = -1;
    410401        numgauss=order_horiz*order_vert;
    411402        coords1=xNew<IssmDouble>(numgauss);
    … …  
    443434
    444435        /*Allocate GaussPenta fields*/
     436   ig = -1;
    445437        numgauss=order_horiz;
    446438        coords1=xNew<IssmDouble>(numgauss);
    … …  
    581573}
    582574/*}}}*/
     575bool GaussPenta::next(void){/*{{{*/
     576
     577        /*Increment Gauss point*/
     578        this->ig++;
     579
     580        /*Have we reached the end?*/
     581        if(this->ig==this->numgauss) return false;
     582
     583        /*If not let's go to the next point*/
     584        _assert_(this->ig>=0 && this->ig< numgauss);
     585        weight=weights[ig];
     586        coord1=coords1[ig];
     587        coord2=coords2[ig];
     588        coord3=coords3[ig];
     589        coord4=coords4[ig];
     590
     591        return true;
     592}/*}}}*/
    583593void GaussPenta::GaussNode(int finiteelement,int iv){/*{{{*/
    584594

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r24119 r25436  
    1414
    1515        private:
    16                 int numgauss;
     16                int         numgauss;   /*Total number of gauss points*/
     17                int         ig;         /*Current gauss point index*/
    1718                IssmDouble* weights;
    1819                IssmDouble* coords1;
    … …  
    4142
    4243                /*Methods*/
     44                bool next(void);
    4345                int  begin(void);
    4446                void Echo(void);

  • issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp

    r23066 r25436  
    1414GaussSeg::GaussSeg(){/*{{{*/
    1515
     16        ig=-1;
    1617        numgauss=-1;
    1718
    … …  
    2930
    3031        /*Get gauss points*/
     32        this->ig       = -1;
    3133        this->numgauss = order;
    3234        GaussLegendreLinear(&pcoords1,&pweights,order);
    … …  
    5456        /*Get gauss points*/
    5557        this->numgauss = 1;
     58        this->ig       = -1;
    5659        this->coords1=xNew<IssmDouble>(numgauss);
    5760        this->weights=xNew<IssmDouble>(numgauss);
    … …  
    117120}
    118121/*}}}*/
     122bool GaussSeg::next(void){/*{{{*/
     123
     124        /*Increment Gauss point*/
     125        this->ig++;
     126
     127        /*Have we reached the end?*/
     128        if(this->ig==this->numgauss) return false;
     129
     130        /*If not let's go to the next point*/
     131         _assert_(this->ig>=0 && this->ig< numgauss);
     132         weight=weights[ig];
     133         coord1=coords1[ig];
     134
     135         return true;
     136}/*}}}*/
    119137int GaussSeg::Enum(void){/*{{{*/
    120138        return GaussSegEnum;

  • issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h

    r20827 r25436  
    1313
    1414        private:
    15                 int numgauss;
    16                 IssmDouble* weights;
     15                int         numgauss;   /*Total number of gauss points*/
     16                int         ig;         /*Current gauss point index*/
     17                IssmDouble* weights;    /*List of weights*/
    1718                IssmDouble* coords1;
    1819
    … …  
    2930
    3031                /*Methods*/
     32                bool next(void);
    3133                int  begin(void);
    3234                void Echo(void);

  • issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp

    r20827 r25436  
    1515GaussTetra::GaussTetra(){/*{{{*/
    1616
     17        ig     = -1;
    1718        numgauss=-1;
    1819
    … …  
    5455        coord2=UNDEF;
    5556        coord3=UNDEF;
     57        ig    = -1;
    5658}
    5759/*}}}*/
    … …  
    8284                _error_(index1 <<" "<<index2 <<" "<<index3 <<" Not supported yet");
    8385        }
     86        this->ig = -1;
    8487}
    8588/*}}}*/
    … …  
    184187}
    185188/*}}}*/
     189bool GaussTetra::next(void){/*{{{*/
     190
     191        /*Increment Gauss point*/
     192        this->ig++;
     193
     194        /*Have we reached the end?*/
     195        if(this->ig==this->numgauss) return false;
     196
     197        /*If not let's go to the next point*/
     198        _assert_(this->ig>=0 && this->ig< numgauss);
     199        weight=weights[ig];
     200        coord1=coords1[ig];
     201        coord2=coords2[ig];
     202        coord3=coords3[ig];
     203        coord4=coords4[ig];
     204
     205        return true;
     206}/*}}}*/
    186207void GaussTetra::GaussNode(int finiteelement,int iv){/*{{{*/
    187208

  • issm/trunk-jpl/src/c/classes/gauss/GaussTetra.h

    r20827 r25436  
    1313
    1414        private:
    15                 int numgauss;
     15                int         numgauss;   /*Total number of gauss points*/
     16                int         ig;         /*Current gauss point index*/
    1617                IssmDouble* weights;
    1718                IssmDouble* coords1;
    … …  
    3536
    3637                /*Methods*/
     38                bool next(void);
    3739                int  begin(void);
    3840                void Echo(void);

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

    r24089 r25436  
    99GaussTria::GaussTria(){/*{{{*/
    1010
     11        ig      = -1;
    1112        numgauss=-1;
    1213
    … …  
    3233        coord2=UNDEF;
    3334        coord3=UNDEF;
     35        ig = -1;
    3436
    3537}
    … …  
    8789
    8890        /*Initialize static fields as undefined*/
     91        ig    = -1;
    8992        weight=UNDEF;
    9093        coord1=UNDEF;
    … …  
    122125
    123126        /*Initialize static fields as undefined*/
     127        ig    = -1;
    124128        weight=UNDEF;
    125129        coord1=UNDEF;
    … …  
    149153         *
    150154         */
    151         int         ig;
     155        int         ii;
    152156        IssmDouble x,y;
    153157        IssmDouble xy_list[3][2];
    … …  
    161165                xy_list[2][0]=0;  xy_list[2][1]=r2;
    162166
    163                 for(ig=0;ig<this->numgauss;ig++){
    164                         x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
    165                         y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
     167                for(ii=0;ii<this->numgauss;ii++){
     168                        x = this->coords1[ii]*xy_list[0][0] + this->coords2[ii]*xy_list[1][0] + this->coords3[ii]*xy_list[2][0];
     169                        y = this->coords1[ii]*xy_list[0][1] + this->coords2[ii]*xy_list[1][1] + this->coords3[ii]*xy_list[2][1];
    166170
    167171                        switch(index){
    168172                                case 0:
    169                                         this->coords1[ig] = 1.-x-y;
    170                                         this->coords2[ig] = x;
    171                                         this->coords3[ig] = y;
     173                                        this->coords1[ii] = 1.-x-y;
     174                                        this->coords2[ii] = x;
     175                                        this->coords3[ii] = y;
    172176                                        break;
    173177                                case 1:
    174                                         this->coords1[ig] = y;
    175                                         this->coords2[ig] = 1.-x-y;
    176                                         this->coords3[ig] = x;
     178                                        this->coords1[ii] = y;
     179                                        this->coords2[ii] = 1.-x-y;
     180                                        this->coords3[ii] = x;
    177181                                        break;
    178182                                case 2:
    179                                         this->coords1[ig] = x;
    180                                         this->coords2[ig] = y;
    181                                         this->coords3[ig] = 1.-x-y;
     183                                        this->coords1[ii] = x;
     184                                        this->coords2[ii] = y;
     185                                        this->coords3[ii] = 1.-x-y;
    182186                                        break;
    183187                                default:
    184188                                        _error_("index "<<index<<" not supported yet");
    185189                        }
    186                         this->weights[ig] = this->weights[ig]*r1*r2;
     190                        this->weights[ii] = this->weights[ii]*r1*r2;
    187191                }
    188192        }
    … …  
    199203
    200204                        //gauss1->Echo();
    201                 for(ig=0;ig<gauss1->numgauss;ig++){
    202                         x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
    203                         y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
     205                for(ii=0;ii<gauss1->numgauss;ii++){
     206                        x = gauss1->coords1[ii]*xy_list[0][0] + gauss1->coords2[ii]*xy_list[1][0] + gauss1->coords3[ii]*xy_list[2][0];
     207                        y = gauss1->coords1[ii]*xy_list[0][1] + gauss1->coords2[ii]*xy_list[1][1] + gauss1->coords3[ii]*xy_list[2][1];
    204208
    205209                        switch(index){
    206210                                case 0:
    207                                         gauss1->coords1[ig] = 1.-x-y;
    208                                         gauss1->coords2[ig] = x;
    209                                         gauss1->coords3[ig] = y;
     211                                        gauss1->coords1[ii] = 1.-x-y;
     212                                        gauss1->coords2[ii] = x;
     213                                        gauss1->coords3[ii] = y;
    210214                                        break;
    211215                                case 1:
    212                                         gauss1->coords1[ig] = y;
    213                                         gauss1->coords2[ig] = 1.-x-y;
    214                                         gauss1->coords3[ig] = x;
     216                                        gauss1->coords1[ii] = y;
     217                                        gauss1->coords2[ii] = 1.-x-y;
     218                                        gauss1->coords3[ii] = x;
    215219                                        break;
    216220                                case 2:
    217                                         gauss1->coords1[ig] = x;
    218                                         gauss1->coords2[ig] = y;
    219                                         gauss1->coords3[ig] = 1.-x-y;
     221                                        gauss1->coords1[ii] = x;
     222                                        gauss1->coords2[ii] = y;
     223                                        gauss1->coords3[ii] = 1.-x-y;
    220224                                        break;
    221225                                default:
    222226                                        _error_("index "<<index<<" not supported yet");
    223227                        }
    224                         gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
     228                        gauss1->weights[ii] = gauss1->weights[ii]*r1*(1-r2);
    225229                }
    226230                        //gauss1->Echo();
    … …  
    230234
    231235                        //gauss2->Echo();
    232                 for(ig=0;ig<gauss2->numgauss;ig++){
    233                         x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
    234                         y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
     236                for(ii=0;ii<gauss2->numgauss;ii++){
     237                        x = gauss2->coords1[ii]*xy_list[0][0] + gauss2->coords2[ii]*xy_list[1][0] + gauss2->coords3[ii]*xy_list[2][0];
     238                        y = gauss2->coords1[ii]*xy_list[0][1] + gauss2->coords2[ii]*xy_list[1][1] + gauss2->coords3[ii]*xy_list[2][1];
    235239
    236240                        switch(index){
    237241                                case 0:
    238                                         gauss2->coords1[ig] = 1.-x-y;
    239                                         gauss2->coords2[ig] = x;
    240                                         gauss2->coords3[ig] = y;
     242                                        gauss2->coords1[ii] = 1.-x-y;
     243                                        gauss2->coords2[ii] = x;
     244                                        gauss2->coords3[ii] = y;
    241245                                        break;
    242246                                case 1:
    243                                         gauss2->coords1[ig] = y;
    244                                         gauss2->coords2[ig] = 1.-x-y;
    245                                         gauss2->coords3[ig] = x;
     247                                        gauss2->coords1[ii] = y;
     248                                        gauss2->coords2[ii] = 1.-x-y;
     249                                        gauss2->coords3[ii] = x;
    246250                                        break;
    247251                                case 2:
    248                                         gauss2->coords1[ig] = x;
    249                                         gauss2->coords2[ig] = y;
    250                                         gauss2->coords3[ig] = 1.-x-y;
     252                                        gauss2->coords1[ii] = x;
     253                                        gauss2->coords2[ii] = y;
     254                                        gauss2->coords3[ii] = 1.-x-y;
    251255                                        break;
    252256                                default:
    253257                                        _error_("index "<<index<<" not supported yet");
    254258                        }
    255                         gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
     259                        gauss2->weights[ii] = gauss2->weights[ii]*(1-r1);
    256260                }
    257261
    … …  
    262266                this->weights=xNew<IssmDouble>(this->numgauss);
    263267
    264                 for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
    265                         this->coords1[ig]=gauss1->coords1[ig];
    266                         this->coords2[ig]=gauss1->coords2[ig];
    267                         this->coords3[ig]=gauss1->coords3[ig];
    268                         this->weights[ig]=gauss1->weights[ig];
     268                for(ii=0;ii<gauss1->numgauss;ii++){ // Add the first triangle gauss points
     269                        this->coords1[ii]=gauss1->coords1[ii];
     270                        this->coords2[ii]=gauss1->coords2[ii];
     271                        this->coords3[ii]=gauss1->coords3[ii];
     272                        this->weights[ii]=gauss1->weights[ii];
    269273                }
    270                 for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
    271                         this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
    272                         this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
    273                         this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
    274                         this->weights[gauss1->numgauss+ig]=gauss2->weights[ig];
     274                for(ii=0;ii<gauss2->numgauss;ii++){ // Add the second triangle gauss points
     275                        this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii];
     276                        this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii];
     277                        this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii];
     278                        this->weights[gauss1->numgauss+ii]=gauss2->weights[ii];
    275279                }
    276280
    … …  
    281285
    282286        /*Initialize static fields as undefined*/
     287        ig    = -1;
    283288        weight=UNDEF;
    284289        coord1=UNDEF;
    … …  
    490495}
    491496/*}}}*/
     497bool GaussTria::next(void){/*{{{*/
     498
     499        /*Increment Gauss point*/
     500        this->ig++;
     501
     502        /*Have we reached the end?*/
     503        if(this->ig==this->numgauss) return false;
     504
     505        /*If not let's go to the next point*/
     506        _assert_(this->ig>=0 && this->ig< numgauss);
     507        weight=weights[ig];
     508        coord1=coords1[ig];
     509        coord2=coords2[ig];
     510        coord3=coords3[ig];
     511
     512        return true;
     513}/*}}}*/
    492514void GaussTria::GaussNode(int finiteelement,int iv){/*{{{*/
    493515

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

    r21892 r25436  
    1313
    1414        private:
    15                 int numgauss;
     15                int         numgauss;   /*Total number of gauss points*/
     16                int         ig;         /*Current gauss point index*/
    1617                IssmDouble* weights;
    1718                IssmDouble* coords1;
    … …  
    3637
    3738                /*Methods*/
     39                bool next(void);
    3840                int  begin(void);
    3941                void Echo(void);
Note: See TracChangeset for help on using the changeset viewer.