Index: ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp	(revision 14778)
+++ ../trunk-jpl/src/c/classes/gauss/GaussTria.cpp	(revision 14779)
@@ -92,6 +92,162 @@
 	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);
Index: ../trunk-jpl/src/c/classes/gauss/GaussTria.h
===================================================================
--- ../trunk-jpl/src/c/classes/gauss/GaussTria.h	(revision 14778)
+++ ../trunk-jpl/src/c/classes/gauss/GaussTria.h	(revision 14779)
@@ -31,6 +31,7 @@
 		GaussTria();
 		GaussTria(int order);
 		GaussTria(int index1,int index2,int order);
+		GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
 		~GaussTria();
 
 		/*Methods*/