Index: issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp
===================================================================
--- issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp	(revision 5827)
+++ issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp	(revision 5831)
@@ -45,6 +45,121 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){{{1*/
-ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){
+/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){{{1*/
+ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){
+
+	/*intermediaries*/
+	int i,j,counter;
+	int gsize,fsize,ssize;
+	int* P=NULL;
+	bool found;
+
+	/*If one of the two matrix is NULL, we copy the other one*/
+	if(!Ke1 && !Ke2){
+		ISSMERROR("Two input element matrices are NULL");
+	}
+	else if(!Ke1){
+		this->Init(Ke2);
+		return;
+	}
+	else if(!Ke2){
+		this->Init(Ke1);
+		return;
+	}
+
+	/*General Case: Ke1 and Ke2 are not empty*/
+	if(!Ke1->square || !Ke2->square) ISSMERROR("merging 2 non square matrices not implemented yet");
+
+	/*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/
+	P=(int*)xmalloc(Ke2->nrows*sizeof(int));
+
+	/*1: Get the new numbering of Ke2 and get size of the new matrix*/
+	gsize=Ke1->nrows;
+	for(i=0;i<Ke2->nrows;i++){
+		found=false;
+		for(j=0;j<Ke1->nrows;j++){
+			if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){
+				found=true; P[i]=j; break;
+			}
+		}
+		if(!found){
+			P[i]=gsize; gsize++;
+		}
+	}
+
+	/*2: Initialize static fields*/
+	this->nrows=gsize;
+	this->ncols=gsize;
+	this->square=true;
+	this->kff=Ke1->kff;
+
+	/*Gset and values*/
+	this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int));
+	this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
+	for(i=0;i<Ke1->nrows;i++){
+		for(j=0;j<Ke1->ncols;j++){
+			this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j];
+		}
+		this->gglobaldoflist[i]=Ke1->gglobaldoflist[i];
+	}
+	for(i=0;i<Ke2->nrows;i++){
+		for(j=0;j<Ke2->ncols;j++){
+			this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j];
+		}
+		this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i];
+	}
+
+	/*Fset*/
+	fsize=Ke1->row_fsize;
+	for(i=0;i<Ke2->row_fsize;i++){
+		if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++;
+	}
+	printf("fsize = %i\n",fsize);
+	ISSMERROR("STOP");
+	if(fsize){
+		this->row_flocaldoflist =(int*)xmalloc(fsize*sizeof(int));
+		this->row_fglobaldoflist=(int*)xmalloc(fsize*sizeof(int));
+		for(i=0;i<Ke1->row_fsize;i++){
+			this->row_flocaldoflist[i]=Ke1->row_flocaldoflist[i];
+		}
+		counter=Ke1->row_fsize;
+		for(i=0;i<Ke2->row_fsize;i++){
+			if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){
+				this->row_flocaldoflist[counter]=P[Ke2->row_flocaldoflist[i]];
+			}
+		}
+	}
+	else{
+		this->row_flocaldoflist=NULL;
+		this->row_fglobaldoflist=NULL;
+	}
+
+	ssize=Ke1->row_ssize+Ke2->row_ssize;
+	this->row_ssize=ssize;
+	if(ssize){
+		ISSMERROR("STOP");
+		this->row_slocaldoflist=(int*)xmalloc(ssize*sizeof(int));
+		this->row_sglobaldoflist=(int*)xmalloc(ssize*sizeof(int));
+		//memcpy(this->row_slocaldoflist,slocaldoflist,ssize*sizeof(int));
+		//memcpy(this->row_sglobaldoflist,sglobaldoflist,ssize*sizeof(int));
+	}
+	else{
+		this->row_slocaldoflist=NULL;
+		this->row_sglobaldoflist=NULL;
+	}
+
+	/*don't do cols, we can pick them up from the rows: */
+	this->col_fsize=0;
+	this->col_flocaldoflist=NULL;
+	this->col_fglobaldoflist=NULL;
+	this->col_ssize=0;
+	this->col_slocaldoflist=NULL;
+	this->col_sglobaldoflist=NULL;
+
+	/*clean-up*/
+	xfree((void**)&P);
+}
+/*}}}*/
+/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){{{1*/
+ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){
 
 	if(square=false)ISSMERROR(" calling square constructor with false square flag!");
@@ -83,8 +198,9 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
-ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
+/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
+ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
 
 	if(square=false)ISSMERROR(" calling square constructor with false square flag!");
+	ISSMASSERT(gsize);
 
 	this->nrows=gsize;
@@ -95,5 +211,9 @@
 	/*fill values with 0: */
 	this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
-	
+
+	/*g dofs*/
+	this->gglobaldoflist=(int*)xmalloc(gsize*sizeof(int));
+	memcpy(this->gglobaldoflist,in_gglobaldoflist,gsize*sizeof(int));
+
 	/*row dofs: */
 	this->row_fsize=fsize;
@@ -121,5 +241,5 @@
 	}
 
-	/*don't do cols, we can't pick them up from the rows: */
+	/*don't do cols, we can pick them up from the rows: */
 	this->col_fsize=0;
 	this->col_flocaldoflist=NULL;
@@ -128,7 +248,4 @@
 	this->col_slocaldoflist=NULL;
 	this->col_sglobaldoflist=NULL;
-
-	/*don't do g-set: */
-	this->gglobaldoflist=NULL;
 
 }
@@ -151,5 +268,5 @@
 
 /*ElementMatrix specific routines: */
-/*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/
+/*FUNCTION ElementMatrix::AddValues{{{1*/
 void ElementMatrix::AddValues(double* Ke_gg){
 
@@ -163,5 +280,5 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/
+/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
 void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){
 
@@ -215,5 +332,5 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::Echo(void){{{1*/
+/*FUNCTION ElementMatrix::Echo{{{1*/
 void ElementMatrix::Echo(void){
 
@@ -262,2 +379,63 @@
 }
 /*}}}*/
+/*FUNCTION ElementMatrix::Init{{{1*/
+void ElementMatrix::Init(ElementMatrix* Ke){
+
+	ISSMASSERT(Ke);
+
+	this->nrows =Ke->nrows;
+	this->ncols =Ke->ncols;
+	this->square=Ke->square;
+	this->kff   =Ke->kff;
+	this->values=(double*)xmalloc(this->nrows*this->ncols*sizeof(double));
+	memcpy(this->values,Ke->values,this->nrows*this->ncols*sizeof(double));
+
+	this->row_fsize=Ke->row_fsize;
+	if(this->row_fsize){
+		this->row_flocaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
+		memcpy(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize*sizeof(int));
+		this->row_fglobaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
+		memcpy(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize*sizeof(int));
+	}
+	else{
+		this->row_flocaldoflist=NULL;
+		this->row_fglobaldoflist=NULL;
+	}
+
+	this->row_ssize=Ke->row_ssize;
+	if(this->row_ssize){
+		this->row_slocaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
+		memcpy(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize*sizeof(int));
+		this->row_sglobaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
+		memcpy(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize*sizeof(int));
+	}
+	else{
+		this->row_slocaldoflist=NULL;
+		this->row_sglobaldoflist=NULL;
+	}
+
+	this->col_fsize=Ke->col_fsize;
+	if(this->col_fsize){
+		this->col_flocaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
+		memcpy(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize*sizeof(int));
+		this->col_fglobaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
+		memcpy(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize*sizeof(int));
+	}
+	else{
+		this->col_flocaldoflist=NULL;
+		this->col_fglobaldoflist=NULL;
+	}
+
+	this->col_ssize=Ke->col_ssize;
+	if(this->col_ssize){
+		this->col_slocaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
+		memcpy(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize*sizeof(int));
+		this->col_sglobaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
+		memcpy(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize*sizeof(int));
+	}
+	else{
+		this->col_slocaldoflist=NULL;
+		this->col_sglobaldoflist=NULL;
+	}
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Numerics/ElementMatrix.h
===================================================================
--- issm/trunk/src/c/objects/Numerics/ElementMatrix.h	(revision 5827)
+++ issm/trunk/src/c/objects/Numerics/ElementMatrix.h	(revision 5831)
@@ -21,7 +21,7 @@
 		int      nrows;
 		int      ncols;
-		double*  values;
 		bool     square;
 		bool     kff;
+		double*  values;
 
 		//gset
@@ -50,6 +50,7 @@
 		/*ElementMatrix constructors, destructors {{{1*/
 		ElementMatrix();
-		ElementMatrix(int gsize,bool square,int* gglobaldoflist);
-		ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
+		ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2);
+		ElementMatrix(bool square,int* gglobaldoflist,int gsize);
+		ElementMatrix(bool square,int* gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
 		~ElementMatrix();
 		/*}}}*/
@@ -58,4 +59,5 @@
 		void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs);
 		void Echo(void);
+		void Init(ElementMatrix* Ke);
 		/*}}}*/
 };
