Ignore:
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:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp

    r5827 r5831  
    4545}
    4646/*}}}*/
    47 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){{{1*/
    48 ElementMatrix::ElementMatrix(int gsize,bool square,int* in_gglobaldoflist){
     47/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){{{1*/
     48ElementMatrix::ElementMatrix(ElementMatrix* Ke1, ElementMatrix* Ke2){
     49
     50        /*intermediaries*/
     51        int i,j,counter;
     52        int gsize,fsize,ssize;
     53        int* P=NULL;
     54        bool found;
     55
     56        /*If one of the two matrix is NULL, we copy the other one*/
     57        if(!Ke1 && !Ke2){
     58                ISSMERROR("Two input element matrices are NULL");
     59        }
     60        else if(!Ke1){
     61                this->Init(Ke2);
     62                return;
     63        }
     64        else if(!Ke2){
     65                this->Init(Ke1);
     66                return;
     67        }
     68
     69        /*General Case: Ke1 and Ke2 are not empty*/
     70        if(!Ke1->square || !Ke2->square) ISSMERROR("merging 2 non square matrices not implemented yet");
     71
     72        /*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/
     73        P=(int*)xmalloc(Ke2->nrows*sizeof(int));
     74
     75        /*1: Get the new numbering of Ke2 and get size of the new matrix*/
     76        gsize=Ke1->nrows;
     77        for(i=0;i<Ke2->nrows;i++){
     78                found=false;
     79                for(j=0;j<Ke1->nrows;j++){
     80                        if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){
     81                                found=true; P[i]=j; break;
     82                        }
     83                }
     84                if(!found){
     85                        P[i]=gsize; gsize++;
     86                }
     87        }
     88
     89        /*2: Initialize static fields*/
     90        this->nrows=gsize;
     91        this->ncols=gsize;
     92        this->square=true;
     93        this->kff=Ke1->kff;
     94
     95        /*Gset and values*/
     96        this->gglobaldoflist=(int*)xmalloc(this->nrows*sizeof(int));
     97        this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
     98        for(i=0;i<Ke1->nrows;i++){
     99                for(j=0;j<Ke1->ncols;j++){
     100                        this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j];
     101                }
     102                this->gglobaldoflist[i]=Ke1->gglobaldoflist[i];
     103        }
     104        for(i=0;i<Ke2->nrows;i++){
     105                for(j=0;j<Ke2->ncols;j++){
     106                        this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j];
     107                }
     108                this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i];
     109        }
     110
     111        /*Fset*/
     112        fsize=Ke1->row_fsize;
     113        for(i=0;i<Ke2->row_fsize;i++){
     114                if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++;
     115        }
     116        printf("fsize = %i\n",fsize);
     117        ISSMERROR("STOP");
     118        if(fsize){
     119                this->row_flocaldoflist =(int*)xmalloc(fsize*sizeof(int));
     120                this->row_fglobaldoflist=(int*)xmalloc(fsize*sizeof(int));
     121                for(i=0;i<Ke1->row_fsize;i++){
     122                        this->row_flocaldoflist[i]=Ke1->row_flocaldoflist[i];
     123                }
     124                counter=Ke1->row_fsize;
     125                for(i=0;i<Ke2->row_fsize;i++){
     126                        if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){
     127                                this->row_flocaldoflist[counter]=P[Ke2->row_flocaldoflist[i]];
     128                        }
     129                }
     130        }
     131        else{
     132                this->row_flocaldoflist=NULL;
     133                this->row_fglobaldoflist=NULL;
     134        }
     135
     136        ssize=Ke1->row_ssize+Ke2->row_ssize;
     137        this->row_ssize=ssize;
     138        if(ssize){
     139                ISSMERROR("STOP");
     140                this->row_slocaldoflist=(int*)xmalloc(ssize*sizeof(int));
     141                this->row_sglobaldoflist=(int*)xmalloc(ssize*sizeof(int));
     142                //memcpy(this->row_slocaldoflist,slocaldoflist,ssize*sizeof(int));
     143                //memcpy(this->row_sglobaldoflist,sglobaldoflist,ssize*sizeof(int));
     144        }
     145        else{
     146                this->row_slocaldoflist=NULL;
     147                this->row_sglobaldoflist=NULL;
     148        }
     149
     150        /*don't do cols, we can pick them up from the rows: */
     151        this->col_fsize=0;
     152        this->col_flocaldoflist=NULL;
     153        this->col_fglobaldoflist=NULL;
     154        this->col_ssize=0;
     155        this->col_slocaldoflist=NULL;
     156        this->col_sglobaldoflist=NULL;
     157
     158        /*clean-up*/
     159        xfree((void**)&P);
     160}
     161/*}}}*/
     162/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){{{1*/
     163ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize){
    49164
    50165        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
     
    83198}
    84199/*}}}*/
    85 /*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
    86 ElementMatrix::ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
     200/*FUNCTION ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){{{1*/
     201ElementMatrix::ElementMatrix(bool square,int* in_gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize){
    87202
    88203        if(square=false)ISSMERROR(" calling square constructor with false square flag!");
     204        ISSMASSERT(gsize);
    89205
    90206        this->nrows=gsize;
     
    95211        /*fill values with 0: */
    96212        this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double));
    97        
     213
     214        /*g dofs*/
     215        this->gglobaldoflist=(int*)xmalloc(gsize*sizeof(int));
     216        memcpy(this->gglobaldoflist,in_gglobaldoflist,gsize*sizeof(int));
     217
    98218        /*row dofs: */
    99219        this->row_fsize=fsize;
     
    121241        }
    122242
    123         /*don't do cols, we can't pick them up from the rows: */
     243        /*don't do cols, we can pick them up from the rows: */
    124244        this->col_fsize=0;
    125245        this->col_flocaldoflist=NULL;
     
    128248        this->col_slocaldoflist=NULL;
    129249        this->col_sglobaldoflist=NULL;
    130 
    131         /*don't do g-set: */
    132         this->gglobaldoflist=NULL;
    133250
    134251}
     
    151268
    152269/*ElementMatrix specific routines: */
    153 /*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/
     270/*FUNCTION ElementMatrix::AddValues{{{1*/
    154271void ElementMatrix::AddValues(double* Ke_gg){
    155272
     
    163280}
    164281/*}}}*/
    165 /*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/
     282/*FUNCTION ElementMatrix::AddToGlobal{{{1*/
    166283void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){
    167284
     
    215332}
    216333/*}}}*/
    217 /*FUNCTION ElementMatrix::Echo(void){{{1*/
     334/*FUNCTION ElementMatrix::Echo{{{1*/
    218335void ElementMatrix::Echo(void){
    219336
     
    262379}
    263380/*}}}*/
     381/*FUNCTION ElementMatrix::Init{{{1*/
     382void ElementMatrix::Init(ElementMatrix* Ke){
     383
     384        ISSMASSERT(Ke);
     385
     386        this->nrows =Ke->nrows;
     387        this->ncols =Ke->ncols;
     388        this->square=Ke->square;
     389        this->kff   =Ke->kff;
     390        this->values=(double*)xmalloc(this->nrows*this->ncols*sizeof(double));
     391        memcpy(this->values,Ke->values,this->nrows*this->ncols*sizeof(double));
     392
     393        this->row_fsize=Ke->row_fsize;
     394        if(this->row_fsize){
     395                this->row_flocaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
     396                memcpy(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize*sizeof(int));
     397                this->row_fglobaldoflist=(int*)xmalloc(this->row_fsize*sizeof(int));
     398                memcpy(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize*sizeof(int));
     399        }
     400        else{
     401                this->row_flocaldoflist=NULL;
     402                this->row_fglobaldoflist=NULL;
     403        }
     404
     405        this->row_ssize=Ke->row_ssize;
     406        if(this->row_ssize){
     407                this->row_slocaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
     408                memcpy(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize*sizeof(int));
     409                this->row_sglobaldoflist=(int*)xmalloc(this->row_ssize*sizeof(int));
     410                memcpy(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize*sizeof(int));
     411        }
     412        else{
     413                this->row_slocaldoflist=NULL;
     414                this->row_sglobaldoflist=NULL;
     415        }
     416
     417        this->col_fsize=Ke->col_fsize;
     418        if(this->col_fsize){
     419                this->col_flocaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
     420                memcpy(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize*sizeof(int));
     421                this->col_fglobaldoflist=(int*)xmalloc(this->col_fsize*sizeof(int));
     422                memcpy(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize*sizeof(int));
     423        }
     424        else{
     425                this->col_flocaldoflist=NULL;
     426                this->col_fglobaldoflist=NULL;
     427        }
     428
     429        this->col_ssize=Ke->col_ssize;
     430        if(this->col_ssize){
     431                this->col_slocaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
     432                memcpy(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize*sizeof(int));
     433                this->col_sglobaldoflist=(int*)xmalloc(this->col_ssize*sizeof(int));
     434                memcpy(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize*sizeof(int));
     435        }
     436        else{
     437                this->col_slocaldoflist=NULL;
     438                this->col_sglobaldoflist=NULL;
     439        }
     440}
     441/*}}}*/
  • issm/trunk/src/c/objects/Numerics/ElementMatrix.h

    r5827 r5831  
    2121                int      nrows;
    2222                int      ncols;
    23                 double*  values;
    2423                bool     square;
    2524                bool     kff;
     25                double*  values;
    2626
    2727                //gset
     
    5050                /*ElementMatrix constructors, destructors {{{1*/
    5151                ElementMatrix();
    52                 ElementMatrix(int gsize,bool square,int* gglobaldoflist);
    53                 ElementMatrix(int gsize,bool square,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
     52                ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2);
     53                ElementMatrix(bool square,int* gglobaldoflist,int gsize);
     54                ElementMatrix(bool square,int* gglobaldoflist,int gsize,int* flocaldoflist,int* fglobaldoflist,int fsize,int* slocaldoflist,int* sglobaldoflist,int ssize);
    5455                ~ElementMatrix();
    5556                /*}}}*/
     
    5859                void AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs);
    5960                void Echo(void);
     61                void Init(ElementMatrix* Ke);
    6062                /*}}}*/
    6163};
Note: See TracChangeset for help on using the changeset viewer.