Changeset 5831 for issm/trunk/src/c/objects/Numerics

Timestamp:

09/15/10 16:30:49 (15 years ago)

Author:

Mathieu Morlighem

Message:

Added merging of 2 matrices, to be completed (almost there)

Location:

issm/trunk/src/c/objects/Numerics

Files:

• ElementMatrix.cpp (modified) (9 diffs)
• ElementMatrix.h (modified) (3 diffs)

issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp

@@ -45 +45 @@
 }
 /*}}}*/
-/*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 +198 @@
 }
 /*}}}*/
-/*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 +211 @@
         /*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 +241 @@
         }
 
-        /*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 +248 @@
         this->col_slocaldoflist=NULL;
         this->col_sglobaldoflist=NULL;
-
-        /*don't do g-set: */
-        this->gglobaldoflist=NULL;
 
 }
@@ -151 +268 @@
 
 /*ElementMatrix specific routines: */
-/*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/
+/*FUNCTION ElementMatrix::AddValues{{{1*/
 void ElementMatrix::AddValues(double* Ke_gg){
 
@@ -163 +280 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/
+/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
 void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){
 
@@ -215 +332 @@
 }
 /*}}}*/
-/*FUNCTION ElementMatrix::Echo(void){{{1*/
+/*FUNCTION ElementMatrix::Echo{{{1*/
 void ElementMatrix::Echo(void){
 
@@ -262 +379 @@
 }
 /*}}}*/
+/*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;
+        }
+}
+/*}}}*/

issm/trunk/src/c/objects/Numerics/ElementMatrix.h

@@ -21 +21 @@
                 int      nrows;
                 int      ncols;
-                double*  values;
                 bool     square;
                 bool     kff;
+                double*  values;
 
                 //gset
@@ -50 +50 @@
                 /*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 +59 @@
                 void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs);
                 void Echo(void);
+                void Init(ElementMatrix* Ke);
                 /*}}}*/
 };