Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 25436) @@ -1451,6 +1451,5 @@ /* Start looping on the number of gaussian points: */ - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); + while(gauss->next()){ friction->GetAlpha2(&alpha2,gauss); @@ -1532,6 +1531,5 @@ /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(xyz_list,xyz_list_boundary,3); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); + while(gauss->next()){ this->GetBSSAFriction(B,element,dim,xyz_list,gauss); @@ -1599,6 +1597,5 @@ /* Start looping on the number of gaussian points: */ Gauss* gauss = element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - gauss->GaussPoint(ig); + while(gauss->next()){ element->JacobianDeterminant(&Jdet,xyz_list,gauss); @@ -1713,7 +1710,5 @@ /* Start looping on the number of gaussian points: */ Gauss* gauss=element->NewGauss(2); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); + while(gauss->next()){ element->JacobianDeterminant(&Jdet,xyz_list,gauss); @@ -1783,7 +1778,5 @@ /*Start looping on Gaussian points*/ Gauss* gauss=element->NewGauss(xyz_list,xyz_list_front,3); - for(int ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); + while(gauss->next()){ thickness_input->GetInputValue(&thickness,gauss); sealevel_input->GetInputValue(&sealevel,gauss); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25436) @@ -230,4 +230,9 @@ break; case P1Enum: + for(int i=0;i(values[i])) values[i] = 0.; + if(xIsNan(values_min[i])) values_min[i] = 0.; + if(xIsNan(values_max[i])) values_max[i] = 0.; + } control_input->SetControl(interpolation_enum,NUMVERTICES,&vertexlids[0],values,values_min,values_max); break; Index: /issm/trunk-jpl/src/c/classes/gauss/Gauss.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/Gauss.h (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/Gauss.h (revision 25436) @@ -11,8 +11,9 @@ IssmDouble weight; - virtual ~Gauss(){}; + virtual ~Gauss(){}; virtual int begin(void)=0; virtual void Echo(void)=0; virtual int end(void)=0; + virtual bool next(void)=0; virtual int Enum(void)=0; virtual void GaussNode(int finitelement,int iv)=0; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp (revision 25436) @@ -16,4 +16,5 @@ GaussPenta::GaussPenta(){/*{{{*/ + ig = -1; numgauss=-1; @@ -34,13 +35,12 @@ /*Intermediaries*/ - int ighoriz,igvert; - int numgauss_horiz; - int numgauss_vert; - IssmDouble *coords1_horiz = NULL; - IssmDouble *coords2_horiz = NULL; - IssmDouble *coords3_horiz = NULL; + int numgauss_horiz; + int numgauss_vert; + IssmDouble *coords1_horiz = NULL; + IssmDouble *coords2_horiz = NULL; + IssmDouble *coords3_horiz = NULL; IssmDouble *weights_horiz = NULL; - double *coords_vert = NULL; - double *weights_vert = NULL; + double *coords_vert = NULL; + double *weights_vert = NULL; /*Get gauss points*/ @@ -50,4 +50,5 @@ /*Allocate GaussPenta fields*/ + ig = -1; numgauss=numgauss_horiz*numgauss_vert; coords1=xNew(numgauss); @@ -58,6 +59,6 @@ /*Combine Horizontal and vertical points*/ - for (ighoriz=0; ighorizig = -1; + } /*}}}*/ @@ -182,4 +186,5 @@ /*Allocate GaussPenta fields*/ + this->ig = -1; numgauss=order_horiz*order_vert; coords1=xNew(numgauss); @@ -236,20 +241,4 @@ 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]; @@ -263,31 +252,31 @@ xy_list[2][0]=0; xy_list[2][1]=r2; - for(ig=0;ignumgauss;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]; + for(int ii=0;iinumgauss;ii++){ + x = this->coords1[ii]*xy_list[0][0] + this->coords2[ii]*xy_list[1][0] + this->coords3[ii]*xy_list[2][0]; + y = this->coords1[ii]*xy_list[0][1] + this->coords2[ii]*xy_list[1][1] + this->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - this->coords1[ig] = 1.-x-y; - this->coords2[ig] = x; - this->coords3[ig] = y; + this->coords1[ii] = 1.-x-y; + this->coords2[ii] = x; + this->coords3[ii] = y; break; case 1: - this->coords1[ig] = y; - this->coords2[ig] = 1.-x-y; - this->coords3[ig] = x; + this->coords1[ii] = y; + this->coords2[ii] = 1.-x-y; + this->coords3[ii] = x; break; case 2: - this->coords1[ig] = x; - this->coords2[ig] = y; - this->coords3[ig] = 1.-x-y; + this->coords1[ii] = x; + this->coords2[ii] = y; + this->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = this->weights[ig]*r1*r2; + this->weights[ii] = this->weights[ii]*r1*r2; } this->coords4=xNew(numgauss); - for(ig=0;igcoords4[ig]=-1.0; + for(int ii=0;iicoords4[ii]=-1.0; } else{ @@ -303,28 +292,28 @@ //gauss1->Echo(); - for(ig=0;ignumgauss;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]; + for(int ii=0;iinumgauss;ii++){ + x = gauss1->coords1[ii]*xy_list[0][0] + gauss1->coords2[ii]*xy_list[1][0] + gauss1->coords3[ii]*xy_list[2][0]; + y = gauss1->coords1[ii]*xy_list[0][1] + gauss1->coords2[ii]*xy_list[1][1] + gauss1->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - gauss1->coords1[ig] = 1.-x-y; - gauss1->coords2[ig] = x; - gauss1->coords3[ig] = y; + gauss1->coords1[ii] = 1.-x-y; + gauss1->coords2[ii] = x; + gauss1->coords3[ii] = y; break; case 1: - gauss1->coords1[ig] = y; - gauss1->coords2[ig] = 1.-x-y; - gauss1->coords3[ig] = x; + gauss1->coords1[ii] = y; + gauss1->coords2[ii] = 1.-x-y; + gauss1->coords3[ii] = x; break; case 2: - gauss1->coords1[ig] = x; - gauss1->coords2[ig] = y; - gauss1->coords3[ig] = 1.-x-y; + gauss1->coords1[ii] = x; + gauss1->coords2[ii] = y; + gauss1->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = gauss1->weights[ig]*r1*(1-r2); + gauss1->weights[ii] = gauss1->weights[ii]*r1*(1-r2); } //gauss1->Echo(); @@ -334,28 +323,28 @@ //gauss2->Echo(); - for(ig=0;ignumgauss;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]; + for(int ii=0;iinumgauss;ii++){ + x = gauss2->coords1[ii]*xy_list[0][0] + gauss2->coords2[ii]*xy_list[1][0] + gauss2->coords3[ii]*xy_list[2][0]; + y = gauss2->coords1[ii]*xy_list[0][1] + gauss2->coords2[ii]*xy_list[1][1] + gauss2->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - gauss2->coords1[ig] = 1.-x-y; - gauss2->coords2[ig] = x; - gauss2->coords3[ig] = y; + gauss2->coords1[ii] = 1.-x-y; + gauss2->coords2[ii] = x; + gauss2->coords3[ii] = y; break; case 1: - gauss2->coords1[ig] = y; - gauss2->coords2[ig] = 1.-x-y; - gauss2->coords3[ig] = x; + gauss2->coords1[ii] = y; + gauss2->coords2[ii] = 1.-x-y; + gauss2->coords3[ii] = x; break; case 2: - gauss2->coords1[ig] = x; - gauss2->coords2[ig] = y; - gauss2->coords3[ig] = 1.-x-y; + gauss2->coords1[ii] = x; + gauss2->coords2[ii] = y; + gauss2->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = gauss2->weights[ig]*(1-r1); + gauss2->weights[ii] = gauss2->weights[ii]*(1-r1); } @@ -367,17 +356,17 @@ this->weights=xNew(this->numgauss); - for(ig=0;ignumgauss;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;ignumgauss;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]; + for(int ii=0;iinumgauss;ii++){ // Add the first triangle gauss points + this->coords1[ii]=gauss1->coords1[ii]; + this->coords2[ii]=gauss1->coords2[ii]; + this->coords3[ii]=gauss1->coords3[ii]; + this->coords4[ii]=gauss1->coords4[ii]; + this->weights[ii]=gauss1->weights[ii]; + } + for(int ii=0;iinumgauss;ii++){ // Add the second triangle gauss points + this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii]; + this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii]; + this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii]; + this->coords4[gauss1->numgauss+ii]=gauss2->coords4[ii]; + this->weights[gauss1->numgauss+ii]=gauss2->weights[ii]; } @@ -393,4 +382,5 @@ coord3=UNDEF; coord4=UNDEF; + ig = -1; } /*}}}*/ @@ -408,4 +398,5 @@ /*Allocate GaussPenta fields*/ + ig = -1; numgauss=order_horiz*order_vert; coords1=xNew(numgauss); @@ -443,4 +434,5 @@ /*Allocate GaussPenta fields*/ + ig = -1; numgauss=order_horiz; coords1=xNew(numgauss); @@ -581,4 +573,22 @@ } /*}}}*/ +bool GaussPenta::next(void){/*{{{*/ + + /*Increment Gauss point*/ + this->ig++; + + /*Have we reached the end?*/ + if(this->ig==this->numgauss) return false; + + /*If not let's go to the next point*/ + _assert_(this->ig>=0 && this->ig< numgauss); + weight=weights[ig]; + coord1=coords1[ig]; + coord2=coords2[ig]; + coord3=coords3[ig]; + coord4=coords4[ig]; + + return true; +}/*}}}*/ void GaussPenta::GaussNode(int finiteelement,int iv){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h (revision 25436) @@ -14,5 +14,6 @@ private: - int numgauss; + int numgauss; /*Total number of gauss points*/ + int ig; /*Current gauss point index*/ IssmDouble* weights; IssmDouble* coords1; @@ -41,4 +42,5 @@ /*Methods*/ + bool next(void); int begin(void); void Echo(void); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.cpp (revision 25436) @@ -14,4 +14,5 @@ GaussSeg::GaussSeg(){/*{{{*/ + ig=-1; numgauss=-1; @@ -29,4 +30,5 @@ /*Get gauss points*/ + this->ig = -1; this->numgauss = order; GaussLegendreLinear(&pcoords1,&pweights,order); @@ -54,4 +56,5 @@ /*Get gauss points*/ this->numgauss = 1; + this->ig = -1; this->coords1=xNew(numgauss); this->weights=xNew(numgauss); @@ -117,4 +120,19 @@ } /*}}}*/ +bool GaussSeg::next(void){/*{{{*/ + + /*Increment Gauss point*/ + this->ig++; + + /*Have we reached the end?*/ + if(this->ig==this->numgauss) return false; + + /*If not let's go to the next point*/ + _assert_(this->ig>=0 && this->ig< numgauss); + weight=weights[ig]; + coord1=coords1[ig]; + + return true; +}/*}}}*/ int GaussSeg::Enum(void){/*{{{*/ return GaussSegEnum; Index: /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussSeg.h (revision 25436) @@ -13,6 +13,7 @@ private: - int numgauss; - IssmDouble* weights; + int numgauss; /*Total number of gauss points*/ + int ig; /*Current gauss point index*/ + IssmDouble* weights; /*List of weights*/ IssmDouble* coords1; @@ -29,4 +30,5 @@ /*Methods*/ + bool next(void); int begin(void); void Echo(void); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.cpp (revision 25436) @@ -15,4 +15,5 @@ GaussTetra::GaussTetra(){/*{{{*/ + ig = -1; numgauss=-1; @@ -54,4 +55,5 @@ coord2=UNDEF; coord3=UNDEF; + ig = -1; } /*}}}*/ @@ -82,4 +84,5 @@ _error_(index1 <<" "<ig = -1; } /*}}}*/ @@ -184,4 +187,22 @@ } /*}}}*/ +bool GaussTetra::next(void){/*{{{*/ + + /*Increment Gauss point*/ + this->ig++; + + /*Have we reached the end?*/ + if(this->ig==this->numgauss) return false; + + /*If not let's go to the next point*/ + _assert_(this->ig>=0 && this->ig< numgauss); + weight=weights[ig]; + coord1=coords1[ig]; + coord2=coords2[ig]; + coord3=coords3[ig]; + coord4=coords4[ig]; + + return true; +}/*}}}*/ void GaussTetra::GaussNode(int finiteelement,int iv){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.h (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTetra.h (revision 25436) @@ -13,5 +13,6 @@ private: - int numgauss; + int numgauss; /*Total number of gauss points*/ + int ig; /*Current gauss point index*/ IssmDouble* weights; IssmDouble* coords1; @@ -35,4 +36,5 @@ /*Methods*/ + bool next(void); int begin(void); void Echo(void); Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp (revision 25436) @@ -9,4 +9,5 @@ GaussTria::GaussTria(){/*{{{*/ + ig = -1; numgauss=-1; @@ -32,4 +33,5 @@ coord2=UNDEF; coord3=UNDEF; + ig = -1; } @@ -87,4 +89,5 @@ /*Initialize static fields as undefined*/ + ig = -1; weight=UNDEF; coord1=UNDEF; @@ -122,4 +125,5 @@ /*Initialize static fields as undefined*/ + ig = -1; weight=UNDEF; coord1=UNDEF; @@ -149,5 +153,5 @@ * */ - int ig; + int ii; IssmDouble x,y; IssmDouble xy_list[3][2]; @@ -161,28 +165,28 @@ xy_list[2][0]=0; xy_list[2][1]=r2; - for(ig=0;ignumgauss;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]; + for(ii=0;iinumgauss;ii++){ + x = this->coords1[ii]*xy_list[0][0] + this->coords2[ii]*xy_list[1][0] + this->coords3[ii]*xy_list[2][0]; + y = this->coords1[ii]*xy_list[0][1] + this->coords2[ii]*xy_list[1][1] + this->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - this->coords1[ig] = 1.-x-y; - this->coords2[ig] = x; - this->coords3[ig] = y; + this->coords1[ii] = 1.-x-y; + this->coords2[ii] = x; + this->coords3[ii] = y; break; case 1: - this->coords1[ig] = y; - this->coords2[ig] = 1.-x-y; - this->coords3[ig] = x; + this->coords1[ii] = y; + this->coords2[ii] = 1.-x-y; + this->coords3[ii] = x; break; case 2: - this->coords1[ig] = x; - this->coords2[ig] = y; - this->coords3[ig] = 1.-x-y; + this->coords1[ii] = x; + this->coords2[ii] = y; + this->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = this->weights[ig]*r1*r2; + this->weights[ii] = this->weights[ii]*r1*r2; } } @@ -199,28 +203,28 @@ //gauss1->Echo(); - for(ig=0;ignumgauss;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]; + for(ii=0;iinumgauss;ii++){ + x = gauss1->coords1[ii]*xy_list[0][0] + gauss1->coords2[ii]*xy_list[1][0] + gauss1->coords3[ii]*xy_list[2][0]; + y = gauss1->coords1[ii]*xy_list[0][1] + gauss1->coords2[ii]*xy_list[1][1] + gauss1->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - gauss1->coords1[ig] = 1.-x-y; - gauss1->coords2[ig] = x; - gauss1->coords3[ig] = y; + gauss1->coords1[ii] = 1.-x-y; + gauss1->coords2[ii] = x; + gauss1->coords3[ii] = y; break; case 1: - gauss1->coords1[ig] = y; - gauss1->coords2[ig] = 1.-x-y; - gauss1->coords3[ig] = x; + gauss1->coords1[ii] = y; + gauss1->coords2[ii] = 1.-x-y; + gauss1->coords3[ii] = x; break; case 2: - gauss1->coords1[ig] = x; - gauss1->coords2[ig] = y; - gauss1->coords3[ig] = 1.-x-y; + gauss1->coords1[ii] = x; + gauss1->coords2[ii] = y; + gauss1->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = gauss1->weights[ig]*r1*(1-r2); + gauss1->weights[ii] = gauss1->weights[ii]*r1*(1-r2); } //gauss1->Echo(); @@ -230,28 +234,28 @@ //gauss2->Echo(); - for(ig=0;ignumgauss;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]; + for(ii=0;iinumgauss;ii++){ + x = gauss2->coords1[ii]*xy_list[0][0] + gauss2->coords2[ii]*xy_list[1][0] + gauss2->coords3[ii]*xy_list[2][0]; + y = gauss2->coords1[ii]*xy_list[0][1] + gauss2->coords2[ii]*xy_list[1][1] + gauss2->coords3[ii]*xy_list[2][1]; switch(index){ case 0: - gauss2->coords1[ig] = 1.-x-y; - gauss2->coords2[ig] = x; - gauss2->coords3[ig] = y; + gauss2->coords1[ii] = 1.-x-y; + gauss2->coords2[ii] = x; + gauss2->coords3[ii] = y; break; case 1: - gauss2->coords1[ig] = y; - gauss2->coords2[ig] = 1.-x-y; - gauss2->coords3[ig] = x; + gauss2->coords1[ii] = y; + gauss2->coords2[ii] = 1.-x-y; + gauss2->coords3[ii] = x; break; case 2: - gauss2->coords1[ig] = x; - gauss2->coords2[ig] = y; - gauss2->coords3[ig] = 1.-x-y; + gauss2->coords1[ii] = x; + gauss2->coords2[ii] = y; + gauss2->coords3[ii] = 1.-x-y; break; default: _error_("index "<weights[ig] = gauss2->weights[ig]*(1-r1); + gauss2->weights[ii] = gauss2->weights[ii]*(1-r1); } @@ -262,15 +266,15 @@ this->weights=xNew(this->numgauss); - for(ig=0;ignumgauss;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(ii=0;iinumgauss;ii++){ // Add the first triangle gauss points + this->coords1[ii]=gauss1->coords1[ii]; + this->coords2[ii]=gauss1->coords2[ii]; + this->coords3[ii]=gauss1->coords3[ii]; + this->weights[ii]=gauss1->weights[ii]; } - for(ig=0;ignumgauss;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]; + for(ii=0;iinumgauss;ii++){ // Add the second triangle gauss points + this->coords1[gauss1->numgauss+ii]=gauss2->coords1[ii]; + this->coords2[gauss1->numgauss+ii]=gauss2->coords2[ii]; + this->coords3[gauss1->numgauss+ii]=gauss2->coords3[ii]; + this->weights[gauss1->numgauss+ii]=gauss2->weights[ii]; } @@ -281,4 +285,5 @@ /*Initialize static fields as undefined*/ + ig = -1; weight=UNDEF; coord1=UNDEF; @@ -490,4 +495,21 @@ } /*}}}*/ +bool GaussTria::next(void){/*{{{*/ + + /*Increment Gauss point*/ + this->ig++; + + /*Have we reached the end?*/ + if(this->ig==this->numgauss) return false; + + /*If not let's go to the next point*/ + _assert_(this->ig>=0 && this->ig< numgauss); + weight=weights[ig]; + coord1=coords1[ig]; + coord2=coords2[ig]; + coord3=coords3[ig]; + + return true; +}/*}}}*/ void GaussTria::GaussNode(int finiteelement,int iv){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 25435) +++ /issm/trunk-jpl/src/c/classes/gauss/GaussTria.h (revision 25436) @@ -13,5 +13,6 @@ private: - int numgauss; + int numgauss; /*Total number of gauss points*/ + int ig; /*Current gauss point index*/ IssmDouble* weights; IssmDouble* coords1; @@ -36,4 +37,5 @@ /*Methods*/ + bool next(void); int begin(void); void Echo(void);