Changeset 26131


Ignore:
Timestamp:
03/23/21 11:52:52 (4 years ago)
Author:
Mathieu Morlighem
Message:

CHG: big simplification here: we now carry the vectors of fset and sset using the same size as gset, with -1 when they are not in the right set. Example: old_fset = [1 2 3], new_fset=[1 -1 2 3 -1 -1]; enjoy

Location:
issm/trunk-jpl/src/c
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp

    r26090 r26131  
    227227        /*Transpose and return Ke*/
    228228        Ke->Transpose();
    229         _assert_(Ke->nrows == numvertices && Ke->ncols == numvertices);
     229        _assert_(Ke->nrows == numvertices);
    230230
    231231        for(int i=0;i<numvertices;i++){
    232232                for(int j=0;j<numvertices;j++){
    233                         ge[i] += Ke->values[i*Ke->ncols + j] * lambda[j];
     233                        ge[i] += Ke->values[i*Ke->nrows + j] * lambda[j];
    234234                }
    235235                //ge[i]=-lambda[i];

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r26073 r26131  
    947947        /*First, figure out size of doflist and create it: */
    948948        int numberofdofs=0;
    949         for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     949        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
    950950
    951951        /*Allocate output*/
    … …  
    955955        int count=0;
    956956        for(int i=0;i<numnodes;i++){
    957                 nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
    958                 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     957                nodes[i]->GetDofList(&doflist[count],approximation_enum,setenum);
     958                count+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
    959959        }
    960960
    … …  
    970970        /*First, figure out size of doflist and create it: */
    971971        int numberofdofs=0;
    972         for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     972        for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
    973973
    974974        /*Allocate output*/
    … …  
    978978        int count=0;
    979979        for(int i=0;i<numnodes;i++){
     980                _assert_(count<numberofdofs);
    980981                nodes[i]->GetDofListLocal(doflist+count,approximation_enum,setenum);
    981                 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     982                count+=nodes[i]->GetNumberOfDofs(approximation_enum,GsetEnum);
    982983        }
    983984
    … …  
    43734374
    43744375        /*Copy current stiffness matrix*/
    4375         values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
    4376         for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];
     4376        values=xNew<IssmDouble>(Ke->nrows*Ke->nrows);
     4377        for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->nrows;j++) values[i*Ke->nrows+j]=Ke->values[i*Ke->nrows+j];
    43774378
    43784379        /*Get Coordinate Systems transform matrix*/
    … …  
    43814382        /*Transform matrix: R*Ke*R^T */
    43824383        TripleMultiply(transform,numdofs,numdofs,0,
    4383                                 values,Ke->nrows,Ke->ncols,0,
     4384                                values,Ke->nrows,Ke->nrows,0,
    43844385                                transform,numdofs,numdofs,1,
    43854386                                &Ke->values[0],0);
    … …  
    45494550
    45504551        /*Copy current stiffness matrix*/
    4551         values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
    4552         for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];
     4552        values=xNew<IssmDouble>(Ke->nrows*Ke->nrows);
     4553        for(int i=0;i<Ke->nrows*Ke->nrows;i++) values[i]=Ke->values[i];
    45534554
    45544555        /*Get Coordinate Systems transform matrix*/
    … …  
    45574558        /*Transform matrix: R^T*Ke*R */
    45584559        TripleMultiply(transform,numdofs,numdofs,1,
    4559                                 values,Ke->nrows,Ke->ncols,0,
     4560                                values,Ke->nrows,Ke->nrows,0,
    45604561                                transform,numdofs,numdofs,0,
    45614562                                &Ke->values[0],0);

  • issm/trunk-jpl/src/c/classes/Node.cpp

    r25610 r26131  
    6868                this->svalues        = xNew<IssmDouble>(this->gsize);
    6969                this->gdoflist       = xNew<int>(this->gsize);
    70                 this->gdoflist_local = xNewZeroInit<int>(this->gsize);
     70                this->gdoflist_local = xNew<int>(this->gsize);
    7171                this->fsize          = -1;
    7272                this->ssize          = -1;
    73                 this->fdoflist       = NULL;
    74                 this->sdoflist       = NULL;
    75                 this->fdoflist_local = NULL;
    76                 this->sdoflist_local = NULL;
     73                this->fdoflist       = xNew<int>(this->gsize);
     74                this->sdoflist       = xNew<int>(this->gsize);
     75                this->fdoflist_local = xNew<int>(this->gsize);
     76                this->sdoflist_local = xNew<int>(this->gsize);
    7777        }
    7878        else{
    … …  
    9696                this->svalues[i]  = 0.;
    9797                this->gdoflist[i] = -1;
     98                this->fdoflist[i] = -1;
     99                this->sdoflist[i] = -1;
     100                this->gdoflist_local[i] = -1;
     101                this->fdoflist_local[i] = -1;
     102                this->sdoflist_local[i] = -1;
    98103        }
    99104
    … …  
    224229        }
    225230        if(output->fsize>0){
    226                 output->fdoflist=xNew<int>(output->fsize);
    227                 output->fdoflist_local=xNew<int>(output->fsize);
     231                output->fdoflist=xNew<int>(output->gsize);
     232                output->fdoflist_local=xNew<int>(output->gsize);
    228233        }
    229234        if(output->ssize>0){
    230                 output->sdoflist=xNew<int>(output->ssize);
    231                 output->sdoflist_local=xNew<int>(output->ssize);
     235                output->sdoflist=xNew<int>(output->gsize);
     236                output->sdoflist_local=xNew<int>(output->gsize);
    232237        }
    233238
    … …  
    239244                memcpy(output->gdoflist,this->gdoflist,output->gsize*sizeof(int));
    240245                memcpy(output->gdoflist_local,this->gdoflist_local,output->gsize*sizeof(int));
    241         }
    242         if(output->fsize>0){
    243                 memcpy(output->fdoflist,this->fdoflist,output->fsize*sizeof(int));
    244                 memcpy(output->fdoflist_local,this->fdoflist_local,output->fsize*sizeof(int));
    245         }
    246         if(output->ssize>0){
    247                 memcpy(output->sdoflist,this->sdoflist,output->ssize*sizeof(int));
    248                 memcpy(output->sdoflist_local,this->sdoflist_local,output->ssize*sizeof(int));
     246                memcpy(output->fdoflist,this->fdoflist,output->gsize*sizeof(int));
     247                memcpy(output->fdoflist_local,this->fdoflist_local,output->gsize*sizeof(int));
     248                memcpy(output->sdoflist,this->sdoflist,output->gsize*sizeof(int));
     249                memcpy(output->sdoflist_local,this->sdoflist_local,output->gsize*sizeof(int));
    249250        }
    250251
    … …  
    277278        marshallhandle->call(this->doftype,gsize);
    278279        marshallhandle->call(this->gdoflist,gsize);
    279         marshallhandle->call(this->fdoflist,fsize);
    280         marshallhandle->call(this->sdoflist,ssize);
     280        marshallhandle->call(this->fdoflist,gsize);
     281        marshallhandle->call(this->sdoflist,gsize);
    281282        marshallhandle->call(this->gdoflist_local,gsize);
    282         marshallhandle->call(this->fdoflist_local,fsize);
    283         marshallhandle->call(this->sdoflist_local,ssize);
     283        marshallhandle->call(this->fdoflist_local,gsize);
     284        marshallhandle->call(this->sdoflist_local,gsize);
    284285} /*}}}*/
    285286
    … …  
    327328        _printf_("\n");
    328329
    329         _printf_("   f_doflist (" << this->fsize << "): |");
    330         for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist[i] << " |");
     330        _printf_("   f_doflist (" << this->gsize << "): |");
     331        for(i=0;i<this->gsize;i++) _printf_(" " << fdoflist[i] << " |");
    331332        _printf_("\n");
    332         _printf_("   f_doflist_local (" << this->fsize << "): |");
    333         for(i=0;i<this->fsize;i++) _printf_(" " << fdoflist_local[i] << " |");
     333        _printf_("   f_doflist_local (" << this->gsize << "): |");
     334        for(i=0;i<this->gsize;i++) _printf_(" " << fdoflist_local[i] << " |");
    334335        _printf_("\n");
    335336
    336         _printf_("   s_doflist (" << this->ssize << "): |");
    337         for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist[i] << " |");
     337        _printf_("   s_doflist (" << this->gsize << "): |");
     338        for(i=0;i<this->gsize;i++) _printf_(" " << sdoflist[i] << " |");
    338339        _printf_("\n");
    339         _printf_("   s_doflist_local (" << this->ssize << "): |");
    340         for(i=0;i<this->ssize;i++) _printf_(" " << sdoflist_local[i] << " |");
     340        _printf_("   s_doflist_local (" << this->gsize << "): |");
     341        for(i=0;i<this->gsize;i++) _printf_(" " << sdoflist_local[i] << " |");
    341342        _printf_("\n");
    342343
    … …  
    386387        }
    387388        else if(setenum==FsetEnum){
    388                 _assert_(dofindex>=0 && dofindex<fsize);
     389                _assert_(dofindex>=0 && dofindex<gsize);
    389390                return fdoflist[dofindex];
    390391        }
    391392        else if(setenum==SsetEnum){
    392                 _assert_(dofindex>=0 && dofindex<ssize);
     393                _assert_(dofindex>=0 && dofindex<gsize);
    393394                return sdoflist[dofindex];
    394395        }
    … …  
    397398} /*}}}*/
    398399void Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
     400        _assert_(!this->indexingupdate);
    399401        int i;
    400         int count=0;
    401         int count2=0;
     402
     403        int* doflistpointer = NULL;
     404        if(setenum==GsetEnum) doflistpointer = gdoflist;
     405        else if(setenum==FsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = fdoflist;
     406        else if(setenum==SsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = sdoflist;
     407        else _error_("not supported");
     408
     409        if(approximation_enum==NoneApproximationEnum){
     410                for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
     411        }
     412        else{
     413                if(doftype){
     414                        int count = 0;
     415                        for(i=0;i<this->gsize;i++){
     416                                if(doftype[i]==approximation_enum) outdoflist[count++]=doflistpointer[i];
     417                        }
     418                }
     419                else for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
     420        }
     421}/*}}}*/
     422void Node::GetDofListLocal(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
    402423
    403424        _assert_(!this->indexingupdate);
     425        int i;
     426
     427        int* doflistpointer = NULL;
     428        if(setenum==GsetEnum) doflistpointer = gdoflist_local;
     429        else if(setenum==FsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = fdoflist_local;
     430        else if(setenum==SsetEnum)for(i=0;i<this->gsize;i++) doflistpointer = sdoflist_local;
     431        else _error_("not supported");
    404432
    405433        if(approximation_enum==NoneApproximationEnum){
    406                 if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist[i];
    407                 if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
    408                 if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
     434                for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
    409435        }
    410436        else{
    411 
    412                 if(setenum==GsetEnum){
    413                         if(doftype){
    414                                 count=0;
    415                                 for(i=0;i<this->gsize;i++){
    416                                         if(doftype[i]==approximation_enum){
    417                                                 outdoflist[count]=gdoflist[i];
    418                                                 count++;
    419                                         }
    420                                 }
    421                                 _assert_(count); //at least one dof should be the approximation requested
    422                         }
    423                         else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist[i];
    424                 }
    425                 else if(setenum==FsetEnum){
    426                         if(doftype){
    427                                 count=0;
    428                                 count2=0;
    429                                 for(i=0;i<this->gsize;i++){
    430                                         if(f_set[i]){
    431                                                 if(doftype[i]==approximation_enum){
    432                                                         outdoflist[count]=fdoflist[count2];
    433                                                         count++;
    434                                                 }
    435                                                 count2++;
    436                                         }
    437                                 }
    438                         }
    439                         else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist[i];
    440                 }
    441                 else if(setenum==SsetEnum){
    442                         if(doftype){
    443                                 count=0;
    444                                 count2=0;
    445                                 for(i=0;i<this->gsize;i++){
    446                                         if(s_set[i]){
    447                                                 if(doftype[i]==approximation_enum){
    448                                                         outdoflist[count]=sdoflist[count2];
    449                                                         count++;
    450                                                 }
    451                                                 count2++;
    452                                         }
    453                                 }
    454                         }
    455                         else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist[i];
    456                 }
    457                 else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    458         }
    459 }
    460 /*}}}*/
    461 void Node::GetDofListLocal(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
    462         int i;
    463         int count=0;
    464         int count2=0;
    465 
    466         _assert_(!this->indexingupdate);
    467 
    468         if(approximation_enum==NoneApproximationEnum){
    469                 if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist_local[i];
    470                 if(setenum==FsetEnum)for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
    471                 if(setenum==SsetEnum)for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
    472         }
    473         else{
    474 
    475                 if(setenum==GsetEnum){
    476                         if(doftype){
    477                                 count=0;
    478                                 for(i=0;i<this->gsize;i++){
    479                                         if(doftype[i]==approximation_enum){
    480                                                 outdoflist[count]=gdoflist_local[i];
    481                                                 count++;
    482                                         }
    483                                 }
    484                                 _assert_(count); //at least one dof should be the approximation requested
    485                         }
    486                         else for(i=0;i<this->gsize;i++) outdoflist[i]=gdoflist_local[i];
    487                 }
    488                 else if(setenum==FsetEnum){
    489                         if(doftype){
    490                                 count=0;
    491                                 count2=0;
    492                                 for(i=0;i<this->gsize;i++){
    493                                         if(f_set[i]){
    494                                                 if(doftype[i]==approximation_enum){
    495                                                         outdoflist[count]=fdoflist_local[count2];
    496                                                         count++;
    497                                                 }
    498                                                 count2++;
    499                                         }
    500                                 }
    501                         }
    502                         else for(i=0;i<this->fsize;i++) outdoflist[i]=fdoflist_local[i];
    503                 }
    504                 else if(setenum==SsetEnum){
    505                         if(doftype){
    506                                 count=0;
    507                                 count2=0;
    508                                 for(i=0;i<this->gsize;i++){
    509                                         if(s_set[i]){
    510                                                 if(doftype[i]==approximation_enum){
    511                                                         outdoflist[count]=sdoflist_local[count2];
    512                                                         count++;
    513                                                 }
    514                                                 count2++;
    515                                         }
    516                                 }
    517                         }
    518                         else for(i=0;i<this->ssize;i++) outdoflist[i]=sdoflist_local[i];
    519                 }
    520                 else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    521         }
    522 }
    523 /*}}}*/
    524 void Node::GetLocalDofList(int* outdoflist,int approximation_enum,int setenum){/*{{{*/
    525         int i;
    526         int count=0;
    527         int count2=0;
    528 
    529         _assert_(!this->indexingupdate);
    530 
    531         if(approximation_enum==NoneApproximationEnum){
    532                 if(setenum==GsetEnum)for(i=0;i<this->gsize;i++) outdoflist[i]=i;
    533                 else if(setenum==FsetEnum){
    534                         count=0;
     437                if(doftype){
     438                        int count =0;
    535439                        for(i=0;i<this->gsize;i++){
    536                                 if(f_set[i]){
    537                                         outdoflist[count]=i;
    538                                         count++;
    539                                 }
    540                         }
    541                 }
    542                 else if(setenum==SsetEnum){
    543                         count=0;
    544                         for(i=0;i<this->gsize;i++){
    545                                 if(s_set[i]){
    546                                         outdoflist[count]=i;
    547                                         count++;
    548                                 }
    549                         }
    550                 }
    551                 else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    552         }
    553         else{
    554 
    555                 if(setenum==GsetEnum){
    556                         if(doftype){
    557                                 count=0;
    558                                 for(i=0;i<this->gsize;i++){
    559                                         if(doftype[i]==approximation_enum){
    560                                                 outdoflist[count]=count;
    561                                                 count++;
    562                                         }
    563                                 }
    564                                 _assert_(count);
    565                         }
    566                         else for(i=0;i<this->gsize;i++) outdoflist[i]=i;
    567                 }
    568                 else if(setenum==FsetEnum){
    569 
    570                         if(doftype){
    571                                 count=0;
    572                                 count2=0;
    573                                 for(i=0;i<this->gsize;i++){
    574                                         if(doftype[i]==approximation_enum){
    575                                                 if(f_set[i]){
    576                                                         outdoflist[count]=count2;
    577                                                         count++;
    578                                                 }
    579                                                 count2++;
    580                                         }
    581                                 }
    582                                 _assert_(count2);
    583                         }
    584                         else{
    585 
    586                                 count=0;
    587                                 for(i=0;i<this->gsize;i++){
    588                                         if(f_set[i]){
    589                                                 outdoflist[count]=i;
    590                                                 count++;
    591                                         }
    592                                 }
    593                         }
    594                 }
    595                 else if(setenum==SsetEnum){
    596                         if(doftype){
    597                                 count=0;
    598                                 count2=0;
    599                                 for(i=0;i<this->gsize;i++){
    600                                         if(doftype[i]==approximation_enum){
    601                                                 if(s_set[i]){
    602                                                         outdoflist[count]=count2;
    603                                                         count++;
    604                                                 }
    605                                                 count2++;
    606                                         }
    607                                 }
    608                                 _assert_(count2);
    609                         }
    610                         else{
    611                                 count=0;
    612                                 for(i=0;i<this->gsize;i++){
    613                                         if(s_set[i]){
    614                                                 outdoflist[count]=i;
    615                                                 count++;
    616                                         }
    617                                 }
    618                         }
    619                 }
    620                 else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
     440                                if(doftype[i]==approximation_enum) outdoflist[count++]=doflistpointer[i];
     441                        }
     442                }
     443                else for(i=0;i<this->gsize;i++) outdoflist[i]=doflistpointer[i];
    621444        }
    622445}
    … …  
    660483void Node::CreateNodalConstraints(Vector<IssmDouble>* ys){/*{{{*/
    661484
    662         int i;
    663         IssmDouble* values=NULL;
    664         int count;
    665 
    666         /*Recover values for s set and plug them in constraints vector: */
    667485        if(this->ssize){
    668                 values=xNew<IssmDouble>(this->ssize);
    669                 count=0;
    670                 for(i=0;i<this->gsize;i++){
    671                         if(this->s_set[i]){
    672                                 values[count]=this->svalues[i];
    673                                 _assert_(!xIsNan<IssmDouble>(values[count]));
    674                                 count++;
    675                         }
    676                 }
    677 
    678486                /*Add values into constraint vector: */
    679                 ys->SetValues(this->ssize,this->sdoflist,values,INS_VAL);
    680         }
    681 
    682         /*Free ressources:*/
    683         xDelete<IssmDouble>(values);
    684 
    685 }
    686 /*}}}*/
     487                ys->SetValues(this->gsize,this->sdoflist,this->svalues,INS_VAL);
     488        }
     489
     490}/*}}}*/
    687491void Node::Deactivate(void){/*{{{*/
    688492
    … …  
    705509        _assert_(this->active);
    706510
    707         if(this->f_set[dof] == 0){
     511        if(!this->f_set[dof]){
    708512                if(this->freeze) _error_("Cannot change dof of frozen node");
    709513                this->indexingupdate = true;
    710                 this->f_set[dof]=1;
    711                 this->s_set[dof]=0;
     514                this->f_set[dof]=true;
     515                this->s_set[dof]=false;
    712516        }
    713517}
    … …  
    719523        _assert_(dof<this->gsize);
    720524
    721         if(this->f_set[dof] == 1){
     525        if(this->f_set[dof]){
    722526                //if(this->freeze) _error_("Cannot change dof of frozen node");
    723527                this->indexingupdate = true;
    724                 this->f_set[dof]=0; //n splits into f (for which we solve) and s (single point constraints)
    725                 this->s_set[dof]=1;
     528                this->f_set[dof]=false; //n splits into f (for which we solve) and s (single point constraints)
     529                this->s_set[dof]=true;
    726530        }
    727531}
    … …  
    800604/*}}}*/
    801605bool Node::IsActive(void){/*{{{*/
    802 
    803606        return active;
    804 
    805 }
    806 /*}}}*/
     607}/*}}}*/
    807608int  Node::IsClone(){/*{{{*/
    808 
    809609        return clone;
    810 
    811 }
    812 /*}}}*/
     610}/*}}}*/
    813611void Node::ReindexingDone(void){/*{{{*/
    814 
    815612        this->indexingupdate = false;
    816 
    817 }
    818 /*}}}*/
     613}/*}}}*/
    819614void Node::RelaxConstraint(int dof){/*{{{*/
    820615
    … …  
    845640                        if(this->f_set[i]){
    846641                                _assert_(local_uf);
    847                                 _assert_(this->fdoflist[count]==indices_uf[this->fdoflist_local[count]]);
    848 
    849                                 values[count] =local_uf[this->fdoflist_local[count]];
     642                                _assert_(this->fdoflist[i]==indices_uf[this->fdoflist_local[i]]);
     643
     644                                values[count] =local_uf[this->fdoflist_local[i]];
    850645                                indices[count]=this->gdoflist[i];
    851646                                count++;
    … …  
    865660                        if(this->s_set[i]){
    866661                                _assert_(local_ys);
    867                                 _assert_(this->sdoflist[count]==indices_ys[this->sdoflist_local[count]]);
    868 
    869                                 values[count] =local_ys[this->sdoflist_local[count]];
     662                                _assert_(this->sdoflist[i]==indices_ys[this->sdoflist_local[i]]);
     663
     664                                values[count] =local_ys[this->sdoflist_local[i]];
    870665                                indices[count]=this->gdoflist[i];
    871666                                count++;
    … …  
    894689                                _assert_(local_ug);
    895690                                _assert_(this->gdoflist[i]==indices_ug[this->gdoflist_local[i]]);
     691                                _assert_(this->fdoflist[i]>=0);
    896692
    897693                                values[count] =local_ug[this->gdoflist_local[i]];
    898                                 indices[count]=this->fdoflist[count];
     694                                indices[count]=this->fdoflist[i];
    899695                                count++;
    900696                        }
    901697                }
     698                _assert_(count==this->fsize);
    902699                uf->SetValues(this->fsize,indices,values,INS_VAL);
    903700
    … …  
    917714                for(int i=0;i<this->gsize;i++) if(f_set[i])size++;
    918715                this->fsize=size;
    919                 xDelete<int>(this->fdoflist);
    920                 xDelete<int>(this->fdoflist_local);
    921 
    922                 if(this->fsize){
    923                         this->fdoflist=xNew<int>(size);
    924                         this->fdoflist_local=xNew<int>(size);
    925                 }
    926                 else{
    927                         this->fdoflist=NULL;
    928                         this->fdoflist_local=NULL;
    929                 }
    930716        }
    931717        else if(setenum==SsetEnum){
    932718                for(int i=0;i<this->gsize;i++) if(s_set[i])size++;
    933719                this->ssize=size;
    934                 xDelete<int>(this->sdoflist);
    935                 xDelete<int>(this->sdoflist_local);
    936 
    937                 if(this->ssize){
    938                         this->sdoflist=xNew<int>(size);
    939                         this->sdoflist_local=xNew<int>(size);
    940                 }
    941                 else{
    942                         this->sdoflist=NULL;
    943                         this->sdoflist_local=NULL;
    944                 }
    945720        }
    946721        /*TODO, we should never be here for GSET! add an assert and track down whether this happens*/
    … …  
    955730        if(setenum==GsetEnum){
    956731                _assert_(this->gsize==0 || this->gdoflist_local);
    957                 for(int i=0;i<this->gsize;i++) gdoflist_local[i]=dofcount+i;
    958                 dofcount+=this->gsize;
     732                for(int i=0;i<this->gsize;i++) gdoflist_local[i]=dofcount++;
    959733        }
    960734        else if(setenum==FsetEnum){
    961735                _assert_(this->fsize==0 || this->fdoflist_local);
    962                 for(int i=0;i<this->fsize;i++) fdoflist_local[i]=dofcount+i;
    963                 dofcount+=this->fsize;
     736                for(int i=0;i<this->gsize;i++){
     737                        if(this->f_set[i]) fdoflist_local[i]=dofcount++;
     738                        else               fdoflist_local[i]=-1;
     739                }
    964740        }
    965741        else if(setenum==SsetEnum){
    966742                _assert_(this->ssize==0 || this->sdoflist_local);
    967                 for(int i=0;i<this->ssize;i++) sdoflist_local[i]=dofcount+i;
    968                 dofcount+=this->ssize;
     743                for(int i=0;i<this->gsize;i++){
     744                        if(this->s_set[i]) sdoflist_local[i]=dofcount++;
     745                        else               sdoflist_local[i]=-1;
     746                }
    969747        }
    970748        else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    … …  
    972750        /*Assign output pointers: */
    973751        *pdofcount=dofcount;
    974 }
    975 /*}}}*/
     752}/*}}}*/
    976753void Node::DistributeGlobalDofsMasters(int dofcount,int setenum){/*{{{*/
    977754
    … …  
    986763        else if(setenum==FsetEnum){
    987764                _assert_(this->fsize==0 || this->fdoflist);
    988                 for(int i=0;i<this->fsize;i++) this->fdoflist[i]=this->fdoflist_local[i]+dofcount;
     765                for(int i=0;i<this->gsize;i++){
     766                        if(this->f_set[i]) fdoflist[i]=this->fdoflist_local[i]+dofcount;
     767                        else               fdoflist[i]=-1;
     768                }
    989769        }
    990770        else if(setenum==SsetEnum){
    991771                _assert_(this->ssize==0 || this->sdoflist);
    992                 for(int i=0;i<this->ssize;i++) this->sdoflist[i]=this->sdoflist_local[i]+dofcount;
     772                for(int i=0;i<this->gsize;i++){
     773                        if(this->s_set[i]) sdoflist[i]=this->sdoflist_local[i]+dofcount;
     774                        else               sdoflist[i]=-1;
     775                }
    993776        }
    994777        else _error_("set of enum type " << EnumToStringx(setenum) << " not supported yet!");
    … …  
    1005788                        break;
    1006789                case FsetEnum:
    1007                         for(int j=0;j<this->fsize;j++) truedofs[j]=fdoflist[j];
     790                        for(int j=0;j<this->gsize;j++) truedofs[j]=fdoflist[j];
    1008791                        break;
    1009792                case SsetEnum:
    1010                         for(int j=0;j<this->ssize;j++) truedofs[j]=sdoflist[j];
     793                        for(int j=0;j<this->gsize;j++) truedofs[j]=sdoflist[j];
    1011794                        break;
    1012795                default:
    … …  
    1027810                        break;
    1028811                case FsetEnum:
    1029                         for(int j=0;j<this->fsize;j++) fdoflist[j]=alltruedofs[j];
     812                        for(int j=0;j<this->gsize;j++) fdoflist[j]=alltruedofs[j];
    1030813                        break;
    1031814                case SsetEnum:
    1032                         for(int j=0;j<this->ssize;j++) sdoflist[j]=alltruedofs[j];
     815                        for(int j=0;j<this->gsize;j++) sdoflist[j]=alltruedofs[j];
    1033816                        break;
    1034817                default:
    … …  
    1041824int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation){/*{{{*/
    1042825
    1043         int  i,numdof,count;
    1044         int* ndof_list=NULL;
     826        /*output*/
    1045827        int *doflist = NULL;
    1046828
    … …  
    1048830
    1049831                /*Allocate:*/
    1050                 ndof_list=xNew<int>(numnodes);
     832                int* ndof_list=xNew<int>(numnodes);
    1051833
    1052834                /*First, figure out size of doflist: */
    1053                 numdof=0;
    1054                 for(i=0;i<numnodes;i++){
    1055                         ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
     835                int numdof=0;
     836                for(int i=0;i<numnodes;i++){
     837                        ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,GsetEnum);
    1056838                        numdof+=ndof_list[i];
    1057839                }
    … …  
    1062844
    1063845                        /*Populate: */
    1064                         count=0;
    1065                         for(i=0;i<numnodes;i++){
     846                        int count=0;
     847                        for(int i=0;i<numnodes;i++){
    1066848                                nodes[i]->GetDofList(&doflist[count],approximation,setenum);
    1067849                                count+=ndof_list[i];
    … …  
    1069851                }
    1070852                else doflist=NULL;
    1071         }
    1072         /*Free ressources:*/
    1073         xDelete<int>(ndof_list);
    1074 
    1075         return doflist;
    1076 }
    1077 /*}}}*/
    1078 int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation){ /*{{{*/
    1079 
    1080         int  i,j,count,numdof,numgdof;
    1081         int* ndof_list=NULL;
    1082         int* ngdof_list_cumulative=NULL;
    1083         int *doflist = NULL;
    1084 
    1085         if(numnodes){
    1086                 /*allocate: */
    1087                 ndof_list=xNew<int>(numnodes);
    1088                 ngdof_list_cumulative=xNew<int>(numnodes);
    1089 
    1090                 /*Get number of dofs per node, and total for this given set*/
    1091                 numdof=0;
    1092                 numgdof=0;
    1093                 for(i=0;i<numnodes;i++){
    1094 
    1095                         /*Cumulative list= number of dofs before node i*/
    1096                         ngdof_list_cumulative[i]=numgdof;
    1097 
    1098                         /*Number of dofs for node i for given set and for the g set*/
    1099                         ndof_list[i]=nodes[i]->GetNumberOfDofs(approximation,setenum);
    1100                         numgdof    +=nodes[i]->GetNumberOfDofs(approximation,GsetEnum);
    1101                         numdof     +=ndof_list[i];
    1102                 }
    1103 
    1104                 if(numdof){
    1105                         /*Allocate: */
    1106                         doflist=xNew<int>(numdof);
    1107 
    1108                         /*Populate: */
    1109                         count=0;
    1110                         for(i=0;i<numnodes;i++){
    1111                                 nodes[i]->GetLocalDofList(&doflist[count],approximation,setenum);
    1112                                 count+=ndof_list[i];
    1113                         }
    1114 
    1115                         /*We now have something like: [0 1 0 2 1 2]. Offset by gsize, to get something like: [0 1 2 4 6 7]:*/
    1116                         count=0;
    1117                         for(i=0;i<numnodes;i++){
    1118                                 for(j=0;j<ndof_list[i];j++){
    1119                                         doflist[count+j]+=ngdof_list_cumulative[i];
    1120                                 }
    1121                                 count+=ndof_list[i];
    1122                         }
    1123                 }
    1124                 else doflist=NULL;
    1125         }
    1126 
    1127         /*Free ressources:*/
    1128         xDelete<int>(ndof_list);
    1129         xDelete<int>(ngdof_list_cumulative);
    1130 
    1131         /*CLean-up and return*/
     853
     854                /*Free ressources:*/
     855                xDelete<int>(ndof_list);
     856        }
     857
    1132858        return doflist;
    1133859}

  • issm/trunk-jpl/src/c/classes/Node.h

    r25508 r26131  
    9494                void  GetDofList(int* poutdoflist,int approximation_enum,int setenum);
    9595                void  GetDofListLocal(int* poutdoflist,int approximation_enum,int setenum);
    96                 void  GetLocalDofList(int* poutdoflist,int approximation_enum,int setenum);
    9796                int   GetNumberOfDofs(int approximation_enum,int setenum);
    9897                void  HardDeactivate(void);
    … …  
    118117/*Methods inherent to Node: */
    119118int* GetGlobalDofList(Node** nodes,int numnodes,int setenum,int approximation);
    120 int* GetLocalDofList(Node** nodes,int numnodes,int setenum,int approximation);
    121119int  GetNumberOfDofs(Node** nodes,int numnodes,int setenum,int approximation);
    122120

  • issm/trunk-jpl/src/c/classes/Nodes.cpp

    r25539 r26131  
    187187         * object that is not a clone, tell them to show their dofs, so that later on, they can get picked
    188188         * up by their clones: */
    189         int  maxdofspernode = this->MaxNumDofs(setenum);
     189        int  maxdofspernode = this->MaxNumDofs(GsetEnum);
    190190        int* truedofs       = xNewZeroInit<int>(this->Size()*maxdofspernode); //only one alloc
    191191        for(int rank=0;rank<num_procs;rank++){

  • issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.cpp

    r25910 r26131  
    2121
    2222        this->nrows=0;
    23         this->ncols=0;
    2423        this->values=NULL;
    25         this->dofsymmetrical=false;
    26 
    27         this->row_fsize=0;
    28         this->row_flocaldoflist=NULL;
    29         this->row_fglobaldoflist=NULL;
    30         this->row_ssize=0;
    31         this->row_slocaldoflist=NULL;
    32         this->row_sglobaldoflist=NULL;
    33 
    34         this->col_fsize=0;
    35         this->col_flocaldoflist=NULL;
    36         this->col_fglobaldoflist=NULL;
    37         this->col_ssize=0;
    38         this->col_slocaldoflist=NULL;
    39         this->col_sglobaldoflist=NULL;
    40 
     24
     25        this->fsize=0;
     26        this->fglobaldoflist=NULL;
     27        this->ssize=0;
     28        this->sglobaldoflist=NULL;
    4129}
    4230/*}}}*/
    … …  
    5341        int i,j,counter;
    5442        int gsize,fsize,ssize;
    55         int* P=NULL;
    56         bool found;
    5743
    5844        /*If one of the two matrix is NULL, we copy the other one*/
    … …  
    7056
    7157        /*General Case: Ke1 and Ke2 are not empty*/
    72         if(!Ke1->dofsymmetrical || !Ke2->dofsymmetrical) _error_("merging 2 non dofsymmetrical matrices not implemented yet");
    7358
    7459        /*Initialize itransformation matrix Ke[P[i]] = Ke2[i]*/
    75         P=xNew<int>(Ke2->nrows);
     60        int* P=xNew<int>(Ke2->nrows);
    7661
    7762        /*1: Get the new numbering of Ke2 and get size of the new matrix*/
    7863        gsize=Ke1->nrows;
    7964        for(i=0;i<Ke2->nrows;i++){
    80                 found=false;
     65                bool found=false;
    8166                for(j=0;j<Ke1->nrows;j++){
    8267                        if(Ke2->gglobaldoflist[i]==Ke1->gglobaldoflist[j]){
    … …  
    9176        /*2: Initialize static fields*/
    9277        this->nrows=gsize;
    93         this->ncols=gsize;
    94         this->dofsymmetrical=true;
    9578
    9679        /*Gset and values*/
    9780        this->gglobaldoflist=xNew<int>(this->nrows);
    98         this->values=xNewZeroInit<IssmDouble>(this->nrows*this->ncols);
     81        this->fglobaldoflist=xNew<int>(this->nrows);
     82        this->sglobaldoflist=xNew<int>(this->nrows);
     83        this->values=xNewZeroInit<IssmDouble>(this->nrows*this->nrows);
    9984        for(i=0;i<Ke1->nrows;i++){
    100                 for(j=0;j<Ke1->ncols;j++){
    101                         this->values[i*this->ncols+j] += Ke1->values[i*Ke1->ncols+j];
     85                for(j=0;j<Ke1->nrows;j++){
     86                        this->values[i*this->nrows+j] += Ke1->values[i*Ke1->nrows+j];
    10287                }
    10388                this->gglobaldoflist[i]=Ke1->gglobaldoflist[i];
     89                this->fglobaldoflist[i]=Ke1->fglobaldoflist[i];
     90                this->sglobaldoflist[i]=Ke1->sglobaldoflist[i];
    10491        }
    10592        for(i=0;i<Ke2->nrows;i++){
    106                 for(j=0;j<Ke2->ncols;j++){
    107                         this->values[P[i]*this->ncols+P[j]] += Ke2->values[i*Ke2->ncols+j];
     93                for(j=0;j<Ke2->nrows;j++){
     94                        this->values[P[i]*this->nrows+P[j]] += Ke2->values[i*Ke2->nrows+j];
    10895                }
    10996                this->gglobaldoflist[P[i]]=Ke2->gglobaldoflist[i];
     97                this->fglobaldoflist[P[i]]=Ke2->fglobaldoflist[i];
     98                this->sglobaldoflist[P[i]]=Ke2->sglobaldoflist[i];
    11099        }
    111100
    112101        /*Fset*/
    113         fsize=Ke1->row_fsize;
    114         for(i=0;i<Ke2->row_fsize;i++){
    115                 if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows) fsize++;
    116         }
    117         this->row_fsize=fsize;
    118         if(fsize){
    119                 this->row_flocaldoflist =xNew<int>(fsize);
    120                 this->row_fglobaldoflist=xNew<int>(fsize);
    121                 for(i=0;i<Ke1->row_fsize;i++){
    122                         this->row_flocaldoflist[i] =Ke1->row_flocaldoflist[i];
    123                         this->row_fglobaldoflist[i]=Ke1->row_fglobaldoflist[i];
    124                 }
    125                 counter=Ke1->row_fsize;
    126                 for(i=0;i<Ke2->row_fsize;i++){
    127                         if(P[Ke2->row_flocaldoflist[i]] >= Ke1->nrows){
    128                                 this->row_flocaldoflist[counter] =P[Ke2->row_flocaldoflist[i]];
    129                                 this->row_fglobaldoflist[counter]=Ke2->row_fglobaldoflist[i];
    130                                 counter++;
    131                         }
    132                 }
    133         }
    134         else{
    135                 this->row_flocaldoflist=NULL;
    136                 this->row_fglobaldoflist=NULL;
    137         }
     102        this->fsize=0;
     103        for(i=0;i<this->nrows;i++) if(this->fglobaldoflist[i]>=0) this->fsize++;
    138104
    139105        /*Sset*/
    140         ssize=Ke1->row_ssize;
    141         for(i=0;i<Ke2->row_ssize;i++){
    142                 if(P[Ke2->row_slocaldoflist[i]] >= Ke1->nrows) ssize++;
    143         }
    144         this->row_ssize=ssize;
    145         if(ssize){
    146                 this->row_slocaldoflist =xNew<int>(ssize);
    147                 this->row_sglobaldoflist=xNew<int>(ssize);
    148                 for(i=0;i<Ke1->row_ssize;i++){
    149                         this->row_slocaldoflist[i] =Ke1->row_slocaldoflist[i];
    150                         this->row_sglobaldoflist[i]=Ke1->row_sglobaldoflist[i];
    151                 }
    152                 counter=Ke1->row_ssize;
    153                 for(i=0;i<Ke2->row_ssize;i++){
    154                         if(P[Ke2->row_slocaldoflist[i]] >= Ke1->nrows){
    155                                 this->row_slocaldoflist[counter] =P[Ke2->row_slocaldoflist[i]];
    156                                 this->row_sglobaldoflist[counter]=Ke2->row_sglobaldoflist[i];
    157                                 counter++;
    158                         }
    159                 }
    160         }
    161         else{
    162                 this->row_slocaldoflist=NULL;
    163                 this->row_sglobaldoflist=NULL;
    164         }
    165 
    166         /*don't do cols, we can pick them up from the rows: */
    167         this->col_fsize=0;
    168         this->col_flocaldoflist=NULL;
    169         this->col_fglobaldoflist=NULL;
    170         this->col_ssize=0;
    171         this->col_slocaldoflist=NULL;
    172         this->col_sglobaldoflist=NULL;
     106        this->ssize=0;
     107        for(i=0;i<this->nrows;i++) if(this->sglobaldoflist[i]>=0) this->ssize++;
    173108
    174109        /*clean-up*/
    … …  
    193128
    194129        /*get Matrix size and properties*/
    195         this->dofsymmetrical=true;
    196130        this->nrows=GetNumberOfDofs(nodes,numnodes,GsetEnum,approximation);
    197         this->ncols=this->nrows;
    198131
    199132        /*fill values with 0: */
    200         this->values=xNewZeroInit<IssmDouble>(this->nrows*this->ncols);
     133        this->values=xNewZeroInit<IssmDouble>(this->nrows*this->nrows);
    201134
    202135        /*g list*/
    … …  
    204137
    205138        /*get dof lists for f and s set: */
    206         this->row_fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
    207         this->row_flocaldoflist =GetLocalDofList( nodes,numnodes,FsetEnum,approximation);
    208         this->row_fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
    209         this->row_ssize=GetNumberOfDofs(nodes,numnodes,SsetEnum,approximation);
    210         this->row_slocaldoflist =GetLocalDofList( nodes,numnodes,SsetEnum,approximation);
    211         this->row_sglobaldoflist=GetGlobalDofList(nodes,numnodes,SsetEnum,approximation);
    212 
    213         /*Because this matrix is "dofsymmetrical" don't do cols, we can pick them up from the rows: */
    214         this->col_fsize=0;
    215         this->col_flocaldoflist=NULL;
    216         this->col_fglobaldoflist=NULL;
    217         this->col_ssize=0;
    218         this->col_slocaldoflist=NULL;
    219         this->col_sglobaldoflist=NULL;
     139        this->fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
     140        this->fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
     141        this->ssize=GetNumberOfDofs(nodes,numnodes,SsetEnum,approximation);
     142        this->sglobaldoflist=GetGlobalDofList(nodes,numnodes,SsetEnum,approximation);
    220143}
    221144/*}}}*/
    … …  
    224147        xDelete<IssmDouble>(this->values);
    225148        xDelete<int>(this->gglobaldoflist);
    226         xDelete<int>(this->row_flocaldoflist);
    227         xDelete<int>(this->row_fglobaldoflist);
    228         xDelete<int>(this->row_slocaldoflist);
    229         xDelete<int>(this->row_sglobaldoflist);
    230         xDelete<int>(this->col_flocaldoflist);
    231         xDelete<int>(this->col_fglobaldoflist);
    232         xDelete<int>(this->col_slocaldoflist);
    233         xDelete<int>(this->col_sglobaldoflist);
     149        xDelete<int>(this->fglobaldoflist);
     150        xDelete<int>(this->sglobaldoflist);
    234151}
    235152/*}}}*/
    … …  
    246163        this->CheckConsistency();
    247164
    248         if(this->dofsymmetrical){
    249                 /*only use row dofs to add values into global matrices: */
    250 
    251                 if(this->row_fsize){
    252                         /*first, retrieve values that are in the f-set from the g-set values matrix: */
    253                         localvalues=xNew<IssmDouble>(this->row_fsize);
    254                         for(int i=0;i<this->row_fsize;i++){
    255                                 localvalues[i] = this->values[this->ncols*this->row_flocaldoflist[i]+ this->row_flocaldoflist[i]];
    256                         }
    257 
    258                         /*add local values into global  matrix, using the fglobaldoflist: */
    259                         pf->SetValues(this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
    260 
    261                         /*Free ressources:*/
    262                         xDelete<IssmDouble>(localvalues);
    263                 }
    264         }
    265         else{
    266                 _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
     165        if(this->fsize){
     166                /*first, retrieve values that are in the f-set from the g-set values matrix: */
     167                localvalues=xNew<IssmDouble>(this->nrows);
     168                for(int i=0;i<this->nrows;i++){
     169                        localvalues[i] = this->values[this->nrows*i + i];
     170                }
     171
     172                /*add local values into global  matrix, using the fglobaldoflist: */
     173                pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
     174
     175                /*Free ressources:*/
     176                xDelete<IssmDouble>(localvalues);
    267177        }
    268178
    … …  
    283193        this->CheckConsistency();
    284194
    285         if(this->dofsymmetrical){
    286                 /*only use row dofs to add values into global matrices: */
    287 
    288                 if(this->row_fsize){
    289                         /*first, retrieve values that are in the f-set from the g-set values matrix: */
    290                         IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_fsize);
    291                         for(int i=0;i<this->row_fsize;i++){
    292                                 for(int j=0;j<this->row_fsize;j++){
    293                                         localvalues[this->row_fsize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]];
    294                                 }
    295                         }
    296 
    297                         /*add local values into global  matrix, using the fglobaldoflist: */
    298                         Kff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
    299 
    300                         /*Free ressources:*/
    301                         xDelete<IssmDouble>(localvalues);
    302                 }
    303 
    304                 if((this->row_ssize!=0) && (this->row_fsize!=0)){
    305                         /*first, retrieve values that are in the f and s-set from the g-set values matrix: */
    306                         IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_ssize);
    307                         for(int i=0;i<this->row_fsize;i++){
    308                                 for(int j=0;j<this->row_ssize;j++){
    309                                         localvalues[this->row_ssize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_slocaldoflist[j]];
    310                                 }
    311                         }
    312                         /*add local values into global  matrix, using the fglobaldoflist: */
    313                         Kfs->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_ssize,this->row_sglobaldoflist,localvalues,ADD_VAL);
    314 
    315                         /*Free ressources:*/
    316                         xDelete<IssmDouble>(localvalues);
    317                 }
    318         }
    319         else{
    320                 _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
    321         }
    322 
     195        /*only use row dofs to add values into global matrices: */
     196        if(this->fsize){
     197                /*add local values into global  matrix, using the fglobaldoflist: */
     198                Kff->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
     199        }
     200
     201        if((this->ssize!=0) && (this->fsize!=0)){
     202                Kfs->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->sglobaldoflist,this->values,ADD_VAL);
     203        }
    323204}
    324205/*}}}*/
    … …  
    331212        this->CheckConsistency();
    332213
    333         if(this->dofsymmetrical){
    334                 /*only use row dofs to add values into global matrices: */
    335 
    336                 if(this->row_fsize){
    337                         /*first, retrieve values that are in the f-set from the g-set values matrix: */
    338                         IssmDouble* localvalues=xNew<IssmDouble>(this->row_fsize*this->row_fsize);
    339                         for(int i=0;i<this->row_fsize;i++){
    340                                 for(int j=0;j<this->row_fsize;j++){
    341                                         localvalues[this->row_fsize*i+j]=this->values[this->ncols*this->row_flocaldoflist[i]+this->row_flocaldoflist[j]];
    342                                 }
    343                         }
    344                         /*add local values into global  matrix, using the fglobaldoflist: */
    345                         Jff->SetValues(this->row_fsize,this->row_fglobaldoflist,this->row_fsize,this->row_fglobaldoflist,localvalues,ADD_VAL);
    346 
    347                         /*Free ressources:*/
    348                         xDelete<IssmDouble>(localvalues);
    349                 }
    350 
    351         }
    352         else{
    353                 _error_("non dofsymmetrical matrix AddToGlobal routine not support yet!");
    354         }
    355 
     214        if(this->fsize){
     215                Jff->SetValues(this->nrows,this->fglobaldoflist,this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
     216        }
    356217}
    357218/*}}}*/
    … …  
    360221        #ifdef _ISSM_DEBUG_
    361222        for (int i=0;i<this->nrows;i++){
    362                 for(int j=0;j<this->ncols;j++){
    363                         if (xIsNan<IssmDouble>(this->values[i*this->ncols+j])) _error_("NaN found in Element Matrix");
    364                         if (xIsInf<IssmDouble>(this->values[i*this->ncols+j])) _error_("Inf found in Element Matrix");
    365                         if (fabs(this->values[i*this->ncols+j])>1.e+50) _error_("Element Matrix values exceeds 1.e+50");
     223                for(int j=0;j<this->nrows;j++){
     224                        if (xIsNan<IssmDouble>(this->values[i*this->nrows+j])) _error_("NaN found in Element Matrix");
     225                        if (xIsInf<IssmDouble>(this->values[i*this->nrows+j])) _error_("Inf found in Element Matrix");
     226                        if (fabs(this->values[i*this->nrows+j])>1.e+50) _error_("Element Matrix values exceeds 1.e+50");
    366227                }
    367228        }
    … …  
    374235        _printf_("Element Matrix echo:\n");
    375236        _printf_("   nrows: " << this->nrows << "\n");
    376         _printf_("   ncols: " << this->ncols << "\n");
    377         _printf_("   dofsymmetrical: " << (dofsymmetrical?"true":"false") << "\n");
    378237
    379238        _printf_("   values:\n");
    380239        for(i=0;i<nrows;i++){
    381240                _printf_(setw(4) << right << i << ": ");
    382                 for(j=0;j<ncols;j++) _printf_( " " << setw(11) << setprecision (5) << right << values[i*ncols+j]);
     241                for(j=0;j<nrows;j++) _printf_( " " << setw(11) << setprecision (5) << right << values[i*nrows+j]);
    383242                _printf_("\n");
    384243        }
    … …  
    387246        if(gglobaldoflist) for(i=0;i<nrows;i++) _printf_(" " << gglobaldoflist[i]); _printf_("\n");
    388247
    389         _printf_("   row_fsize: " << row_fsize << "\n");
    390         _printf_("   row_flocaldoflist  (" << row_flocaldoflist << "): ");
    391         if(row_flocaldoflist) for(i=0;i<row_fsize;i++) _printf_(" " << row_flocaldoflist[i] ); _printf_(" \n");
    392         _printf_("   row_fglobaldoflist  (" << row_fglobaldoflist << "): ");
    393         if(row_fglobaldoflist) for(i=0;i<row_fsize;i++)_printf_(" " << row_fglobaldoflist[i]); _printf_(" \n");
    394 
    395         _printf_("   row_ssize: " << row_ssize << "\n");
    396         _printf_("   row_slocaldoflist  (" << row_slocaldoflist << "): ");
    397         if(row_slocaldoflist) for(i=0;i<row_ssize;i++) _printf_(" " << row_slocaldoflist[i] ); _printf_(" \n");
    398         _printf_("   row_sglobaldoflist  (" << row_sglobaldoflist << "): ");
    399         if(row_sglobaldoflist) for(i=0;i<row_ssize;i++)_printf_(" " << row_sglobaldoflist[i]); _printf_(" \n");
    400 
    401         if(!dofsymmetrical){
    402                 _printf_("   col_fsize: " << col_fsize << "\n");
    403                 _printf_("   col_flocaldoflist  (" << col_flocaldoflist << "): ");
    404                 if(col_flocaldoflist) for(i=0;i<col_fsize;i++) _printf_(" " << col_flocaldoflist[i] ); _printf_(" \n");
    405                 _printf_("   col_fglobaldoflist  (" << col_fglobaldoflist << "): ");
    406                 if(col_fglobaldoflist) for(i=0;i<col_fsize;i++)_printf_(" " << col_fglobaldoflist[i]); _printf_(" \n");
    407 
    408                 _printf_("   col_ssize: " << col_ssize << "\n");
    409                 _printf_("   col_slocaldoflist  (" << col_slocaldoflist << "): ");
    410                 if(col_slocaldoflist) for(i=0;i<col_ssize;i++) _printf_(" " << col_slocaldoflist[i] ); _printf_(" \n");
    411                 _printf_("   col_sglobaldoflist  (" << col_sglobaldoflist << "): ");
    412                 if(col_sglobaldoflist) for(i=0;i<col_ssize;i++)_printf_(" " << col_sglobaldoflist[i]); _printf_(" \n");
    413         }
     248        _printf_("   fsize: " << fsize << "\n");
     249        _printf_("   fglobaldoflist  (" << fglobaldoflist << "): ");
     250        for(i=0;i<nrows;i++)_printf_(" " << fglobaldoflist[i]); _printf_(" \n");
     251
     252        _printf_("   ssize: " << ssize << "\n");
     253        _printf_("   sglobaldoflist  (" << sglobaldoflist << "): ");
     254        for(i=0;i<nrows;i++)_printf_(" " << sglobaldoflist[i]); _printf_(" \n");
    414255}
    415256/*}}}*/
    … …  
    420261
    421262        this->nrows =Ke->nrows;
    422         this->ncols =Ke->ncols;
    423         this->dofsymmetrical=Ke->dofsymmetrical;
    424 
    425         this->values=xNew<IssmDouble>(this->nrows*this->ncols);
    426         xMemCpy<IssmDouble>(this->values,Ke->values,this->nrows*this->ncols);
     263
     264        this->values=xNew<IssmDouble>(this->nrows*this->nrows);
     265        xMemCpy<IssmDouble>(this->values,Ke->values,this->nrows*this->nrows);
    427266
    428267        this->gglobaldoflist=xNew<int>(this->nrows);
    429268        xMemCpy<int>(this->gglobaldoflist,Ke->gglobaldoflist,this->nrows);
    430269
    431         this->row_fsize=Ke->row_fsize;
    432         if(this->row_fsize){
    433                 this->row_flocaldoflist=xNew<int>(this->row_fsize);
    434                 xMemCpy<int>(this->row_flocaldoflist,Ke->row_flocaldoflist,this->row_fsize);
    435                 this->row_fglobaldoflist=xNew<int>(this->row_fsize);
    436                 xMemCpy<int>(this->row_fglobaldoflist,Ke->row_fglobaldoflist,this->row_fsize);
    437         }
    438         else{
    439                 this->row_flocaldoflist=NULL;
    440                 this->row_fglobaldoflist=NULL;
    441         }
    442 
    443         this->row_ssize=Ke->row_ssize;
    444         if(this->row_ssize){
    445                 this->row_slocaldoflist=xNew<int>(this->row_ssize);
    446                 xMemCpy<int>(this->row_slocaldoflist,Ke->row_slocaldoflist,this->row_ssize);
    447                 this->row_sglobaldoflist=xNew<int>(this->row_ssize);
    448                 xMemCpy<int>(this->row_sglobaldoflist,Ke->row_sglobaldoflist,this->row_ssize);
    449         }
    450         else{
    451                 this->row_slocaldoflist=NULL;
    452                 this->row_sglobaldoflist=NULL;
    453         }
    454 
    455         this->col_fsize=Ke->col_fsize;
    456         if(this->col_fsize){
    457                 this->col_flocaldoflist=xNew<int>(this->col_fsize);
    458                 xMemCpy<int>(this->col_flocaldoflist,Ke->col_flocaldoflist,this->col_fsize);
    459                 this->col_fglobaldoflist=xNew<int>(this->col_fsize);
    460                 xMemCpy<int>(this->col_fglobaldoflist,Ke->col_fglobaldoflist,this->col_fsize);
    461         }
    462         else{
    463                 this->col_flocaldoflist=NULL;
    464                 this->col_fglobaldoflist=NULL;
    465         }
    466 
    467         this->col_ssize=Ke->col_ssize;
    468         if(this->col_ssize){
    469                 this->col_slocaldoflist=xNew<int>(this->col_ssize);
    470                 xMemCpy<int>(this->col_slocaldoflist,Ke->col_slocaldoflist,this->col_ssize);
    471                 this->col_sglobaldoflist=xNew<int>(this->col_ssize);
    472                 xMemCpy<int>(this->col_sglobaldoflist,Ke->col_sglobaldoflist,this->col_ssize);
    473         }
    474         else{
    475                 this->col_slocaldoflist=NULL;
    476                 this->col_sglobaldoflist=NULL;
    477         }
     270        this->fsize=Ke->fsize;
     271        this->fglobaldoflist=xNew<int>(this->nrows);
     272        xMemCpy<int>(this->fglobaldoflist,Ke->fglobaldoflist,this->nrows);
     273
     274        this->ssize=Ke->ssize;
     275        this->sglobaldoflist=xNew<int>(this->nrows);
     276        xMemCpy<int>(this->sglobaldoflist,Ke->sglobaldoflist,this->nrows);
    478277}
    479278/*}}}*/
    480279void ElementMatrix::Lump(void){/*{{{*/
    481280
    482         if(!dofsymmetrical) _error_("not supported yet");
    483 
    484281        for(int i=0;i<this->nrows;i++){
    485                 for(int j=0;j<this->ncols;j++){
     282                for(int j=0;j<this->nrows;j++){
    486283                        if(i!=j){
    487                                 this->values[i*this->ncols+i] += this->values[i*this->ncols+j];
    488                                 this->values[i*this->ncols+j]  = 0.;
     284                                this->values[i*this->nrows+i] += this->values[i*this->nrows+j];
     285                                this->values[i*this->nrows+j]  = 0.;
    489286                        }
    490287                }
    … …  
    499296        ElementMatrix* Ke_copy=new ElementMatrix(this);
    500297
    501         /*Update sizes*/
    502         this->nrows=Ke_copy->ncols;
    503         this->ncols=Ke_copy->nrows;
    504 
    505298        /*Transpose values*/
    506         for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[i*this->ncols+j]=Ke_copy->values[j*Ke_copy->ncols+i];
    507 
    508         /*Transpose indices*/
    509         if(!dofsymmetrical){
    510                 _error_("not supported yet");
    511         }
     299        for (int i=0;i<this->nrows;i++) for(int j=0;j<this->nrows;j++) this->values[i*this->nrows+j]=Ke_copy->values[j*Ke_copy->nrows+i];
    512300
    513301        /*Clean up and return*/
    … …  
    534322
    535323        /*Checks in debugging mode*/
    536         _assert_(this->nrows==this->ncols && bsize>0 && bsize<this->ncols && this->values);
     324        _assert_(bsize>0 && bsize<this->nrows && this->values);
    537325
    538326        /*Intermediaries*/
    … …  
    571359        Kbi = xNew<IssmDouble>(bsize*isize);
    572360        Kbb = xNew<IssmDouble>(bsize*bsize);
    573         for(i=0;i<isize;i++) for(j=0;j<isize;j++) Kii[i*isize+j] = this->values[iindices[i]*this->ncols + iindices[j]];
    574         for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = this->values[iindices[i]*this->ncols + bindices[j]];
    575         for(i=0;i<bsize;i++) for(j=0;j<isize;j++) Kbi[i*isize+j] = this->values[bindices[i]*this->ncols + iindices[j]];
    576         for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = this->values[bindices[i]*this->ncols + bindices[j]];
     361        for(i=0;i<isize;i++) for(j=0;j<isize;j++) Kii[i*isize+j] = this->values[iindices[i]*this->nrows + iindices[j]];
     362        for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = this->values[iindices[i]*this->nrows + bindices[j]];
     363        for(i=0;i<bsize;i++) for(j=0;j<isize;j++) Kbi[i*isize+j] = this->values[bindices[i]*this->nrows + iindices[j]];
     364        for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = this->values[bindices[i]*this->nrows + bindices[j]];
    577365
    578366        /*Invert Kbb*/
    … …  
    601389
    602390        /*Update matrix values*/
    603         for(i=0;i<this->nrows*this->ncols;i++) this->values[i]=0.;
     391        for(i=0;i<this->nrows*this->nrows;i++) this->values[i]=0.;
    604392        for(i=0;i<isize;i++){
    605393                for(j=0;j<isize;j++){
    606                         this->values[iindices[i]*this->ncols + iindices[j]] = Ktemp[i*isize+j];
     394                        this->values[iindices[i]*this->nrows + iindices[j]] = Ktemp[i*isize+j];
    607395                }
    608396        }

  • issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h

    r20827 r26131  
    2222
    2323                int      nrows;
    24                 int      ncols;
    25                 bool     dofsymmetrical;
     24                int      fsize;
     25                int      ssize;
    2626                IssmDouble*  values;
    2727
    28                 //gset
    2928                int*     gglobaldoflist;
    30 
    31                 /*row wise: */
    32                 //fset
    33                 int      row_fsize;
    34                 int*     row_flocaldoflist;
    35                 int*     row_fglobaldoflist;
    36                 //sset
    37                 int      row_ssize;
    38                 int*     row_slocaldoflist;
    39                 int*     row_sglobaldoflist;
    40 
    41                 /*column wise: */
    42                 //fset
    43                 int      col_fsize;
    44                 int*     col_flocaldoflist;
    45                 int*     col_fglobaldoflist;
    46                 //sset
    47                 int      col_ssize;
    48                 int*     col_slocaldoflist;
    49                 int*     col_sglobaldoflist;
     29                int*     fglobaldoflist;
     30                int*     sglobaldoflist;
    5031
    5132                /*ElementMatrix constructors, destructors*/

  • issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp

    r25910 r26131  
    2121
    2222        this->nrows=0;
     23        this->fsize=0;
    2324        this->values=NULL;
    24         this->fsize=0;
    25         this->flocaldoflist=NULL;
    2625        this->fglobaldoflist=NULL;
    2726
    … …  
    3130
    3231        /*intermediaries*/
    33         int i,j,counter;
     32        int i,j;
    3433        int gsize,fsize;
    35         int* P=NULL;
    36         bool found;
    3734
    3835        /*If one of the two matrix is NULL, we copy the other one*/
    … …  
    5047
    5148        /*Initialize itransformation matrix pe[P[i]] = pe2[i]*/
    52         P=xNew<int>(pe2->nrows);
     49        int* P=xNew<int>(pe2->nrows);
    5350
    5451        /*1: Get the new numbering of pe2 and get size of the new matrix*/
    5552        gsize=pe1->nrows;
    5653        for(i=0;i<pe2->nrows;i++){
    57                 found=false;
     54                bool found=false;
    5855                for(j=0;j<pe1->nrows;j++){
    5956                        if(pe2->gglobaldoflist[i]==pe1->gglobaldoflist[j]){
    … …  
    7168        /*Gset and values*/
    7269        this->gglobaldoflist=xNew<int>(this->nrows);
     70        this->fglobaldoflist=xNew<int>(this->nrows);
    7371        this->values=xNewZeroInit<IssmDouble>(this->nrows);
    7472        for(i=0;i<pe1->nrows;i++){
    7573                this->values[i] += pe1->values[i];
    7674                this->gglobaldoflist[i]=pe1->gglobaldoflist[i];
     75                this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
    7776        }
    7877        for(i=0;i<pe2->nrows;i++){
    7978                this->values[P[i]] += pe2->values[i];
    8079                this->gglobaldoflist[P[i]]=pe2->gglobaldoflist[i];
     80                this->fglobaldoflist[P[i]]=pe2->fglobaldoflist[i];
    8181        }
    8282
    8383        /*Fset*/
    84         fsize=pe1->fsize;
    85         for(i=0;i<pe2->fsize;i++){
    86                 if(P[pe2->flocaldoflist[i]] >= pe1->nrows) fsize++;
    87         }
    88         this->fsize=fsize;
    89         if(fsize){
    90                 this->flocaldoflist =xNew<int>(fsize);
    91                 this->fglobaldoflist=xNew<int>(fsize);
    92                 for(i=0;i<pe1->fsize;i++){
    93                         this->flocaldoflist[i] =pe1->flocaldoflist[i];
    94                         this->fglobaldoflist[i]=pe1->fglobaldoflist[i];
    95                 }
    96                 counter=pe1->fsize;
    97                 for(i=0;i<pe2->fsize;i++){
    98                         if(P[pe2->flocaldoflist[i]] >= pe1->nrows){
    99                                 this->flocaldoflist[counter] =P[pe2->flocaldoflist[i]];
    100                                 this->fglobaldoflist[counter]=pe2->fglobaldoflist[i];
    101                                 counter++;
    102                         }
    103                 }
    104         }
    105         else{
    106                 this->flocaldoflist=NULL;
    107                 this->fglobaldoflist=NULL;
    108         }
     84        this->fsize=0;
     85        for(i=0;i<this->nrows;i++) if(this->fglobaldoflist[i]>=0) this->fsize++;
    10986
    11087        /*clean-up*/
    … …  
    139116        /*Get fsize*/
    140117        this->fsize=GetNumberOfDofs(nodes,numnodes,FsetEnum,approximation);
    141         this->flocaldoflist =GetLocalDofList( nodes,numnodes,FsetEnum,approximation);
    142118        this->fglobaldoflist=GetGlobalDofList(nodes,numnodes,FsetEnum,approximation);
    143119}
    144120/*}}}*/
    145121ElementVector::~ElementVector(){/*{{{*/
    146 
    147122        xDelete<IssmDouble>(this->values);
    148123        xDelete<int>(this->gglobaldoflist);
    149         xDelete<int>(this->flocaldoflist);
    150124        xDelete<int>(this->fglobaldoflist);
    151125}
    … …  
    155129void ElementVector::AddToGlobal(Vector<IssmDouble>* pf){/*{{{*/
    156130
    157         int i;
    158         IssmDouble* localvalues=NULL;
    159 
    160131        /*In debugging mode, check consistency (no NaN, and values not too big)*/
    161132        this->CheckConsistency();
    162133
    163134        if(this->fsize){
    164                 /*first, retrieve values that are in the f-set from the g-set values vector: */
    165                 localvalues=xNew<IssmDouble>(this->fsize);
    166                 for(i=0;i<this->fsize;i++){
    167                         localvalues[i]=this->values[this->flocaldoflist[i]];
    168                 }
    169                 /*add local values into global  vector, using the fglobaldoflist: */
    170                 pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,ADD_VAL);
    171 
    172                 /*Free ressources:*/
    173                 xDelete<IssmDouble>(localvalues);
     135                pf->SetValues(this->nrows,this->fglobaldoflist,this->values,ADD_VAL);
    174136        }
    175137
    … …  
    193155        _printf_("Element Vector echo:\n");
    194156        _printf_("   nrows: " << nrows << "\n");
     157        _printf_("   fsize: " << fsize << "\n");
    195158        _printf_("   values:\n");
    196         for(i=0;i<nrows;i++){
    197                 _printf_(setw(4) << right << i << ": " << setw(10) << values[i] << "\n");
    198         }
     159        for(i=0;i<nrows;i++) _printf_(setw(4) << right << i << ": " << setw(10) << values[i] << "\n");
    199160
    200161        _printf_("   gglobaldoflist (" << gglobaldoflist << "): ");
    … …  
    202163        _printf_(" \n");
    203164
    204         _printf_("   fsize: " << fsize << "\n");
    205         _printf_("   flocaldoflist  (" << flocaldoflist << "): ");
    206         if(flocaldoflist) for(i=0;i<fsize;i++) _printf_(" " << flocaldoflist[i] );
    207         _printf_(" \n");
    208165        _printf_("   fglobaldoflist (" << fglobaldoflist << "): ");
    209166        if(fglobaldoflist) for(i=0;i<fsize;i++) _printf_(" " << fglobaldoflist[i] );
    … …  
    216173
    217174        this->nrows =pe->nrows;
     175        this->fsize =pe->fsize;
    218176
    219177        this->values=xNew<IssmDouble>(this->nrows);
    … …  
    223181        xMemCpy<int>(this->gglobaldoflist,pe->gglobaldoflist,this->nrows);
    224182
    225         this->fsize=pe->fsize;
     183        this->fglobaldoflist=xNew<int>(this->nrows);
     184        xMemCpy<int>(this->fglobaldoflist,pe->fglobaldoflist,this->nrows);
     185}
     186/*}}}*/
     187void ElementVector::InsertIntoGlobal(Vector<IssmDouble>* pf){/*{{{*/
     188
    226189        if(this->fsize){
    227                 this->flocaldoflist=xNew<int>(this->fsize);
    228                 xMemCpy<int>(this->flocaldoflist,pe->flocaldoflist,this->fsize);
    229                 this->fglobaldoflist=xNew<int>(this->fsize);
    230                 xMemCpy<int>(this->fglobaldoflist,pe->fglobaldoflist,this->fsize);
    231         }
    232         else{
    233                 this->flocaldoflist=NULL;
    234                 this->fglobaldoflist=NULL;
    235         }
    236 }
    237 /*}}}*/
    238 void ElementVector::InsertIntoGlobal(Vector<IssmDouble>* pf){/*{{{*/
    239 
    240         if(this->fsize){
    241                 /*first, retrieve values that are in the f-set from the g-set values vector: */
    242                 IssmDouble* localvalues=xNew<IssmDouble>(this->fsize);
    243                 for(int i=0;i<this->fsize;i++){
    244                         localvalues[i]=this->values[this->flocaldoflist[i]];
    245                 }
    246190                /*add local values into global  vector, using the fglobaldoflist: */
    247                 pf->SetValues(this->fsize,this->fglobaldoflist,localvalues,INS_VAL);
    248 
    249                 /*Free ressources:*/
    250                 xDelete<IssmDouble>(localvalues);
     191                pf->SetValues(this->nrows,this->fglobaldoflist,this->values,INS_VAL);
    251192        }
    252193
    … …  
    280221        /*Checks in debugging mode*/
    281222        _assert_(bsize>0 && bsize<this->nrows && this->values && Ke);
    282         _assert_(Ke->nrows==Ke->ncols && this->nrows==Ke->nrows);
     223        _assert_(this->nrows==Ke->nrows);
    283224
    284225        /*Intermediaries*/
    … …  
    317258        Fb  = xNew<IssmDouble>(bsize);
    318259        Fi  = xNew<IssmDouble>(isize);
    319         for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->ncols + bindices[j]];
    320         for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->ncols + bindices[j]];
     260        for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->nrows + bindices[j]];
     261        for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->nrows + bindices[j]];
    321262        for(i=0;i<bsize;i++) Fb[i] = this->values[bindices[i]];
    322263        for(i=0;i<isize;i++) Fi[i] = this->values[iindices[i]];

  • issm/trunk-jpl/src/c/classes/matrix/ElementVector.h

    r20827 r26131  
    2222        public:
    2323                int         nrows;
    24                 IssmDouble* values;
    25 
    26                 //gset
    27                 int* gglobaldoflist;
    28 
    29                 //fset
    30                 int  fsize;
    31                 int* flocaldoflist;
    32                 int* fglobaldoflist;
     24                int         fsize;
     25                IssmDouble *values;
     26                int        *gglobaldoflist;
     27                int        *fglobaldoflist;
    3328
    3429                /*ElementVector constructors, destructors*/

  • issm/trunk-jpl/src/c/toolkits/issm/Bucket.h

    r25508 r26131  
    3434                        this->Initialize();
    3535                } /*}}}*/
    36                 Bucket(int min,int* idxmin,int nin,int* idxnin,doubletype* valuesin,InsMode modein){ /*{{{*/
     36                Bucket(int m_in,int* idxm_in,int n_in,int* idxn_in,doubletype* values_in,InsMode mode_in){ /*{{{*/
    3737
    3838                        this->Initialize();
    3939
     40                        /*Check whether we have negative indices?*/
     41                        int m_numneg = 0;
     42                        int n_numneg = 0;
     43                        for(int i=0;i<m_in;i++) if(idxm_in[i]<0) m_numneg++;
     44                        for(int i=0;i<n_in;i++) if(idxn_in[i]<0) n_numneg++;
     45
    4046                        this->type=MATRIX_BUCKET;
    41                         this->m=min;
    42                         this->n=nin;
    43                         this->mode=modein;
     47                        this->m=m_in-m_numneg;
     48                        this->n=n_in-n_numneg;
     49                        this->mode=mode_in;
     50                        if(this->m*this->n){
     51                                this->idxm=xNew<int>(this->m);
     52                                this->idxn=xNew<int>(this->n);
     53                                this->values=xNew<doubletype>(this->n*this->m);
     54
     55                                if(m_numneg==0 && n_numneg==0){
     56                                        xMemCpy(this->values,values_in,this->n*this->m);
     57                                        xMemCpy(this->idxm,idxm_in,this->m);
     58                                        xMemCpy(this->idxn,idxn_in,this->n);
     59                                }
     60                                else{
     61                                        int count;
     62                                        count =0; for(int i=0;i<m_in;i++) if(idxm_in[i]>=0) this->idxm[count++] = idxm_in[i];
     63                                        count =0; for(int i=0;i<n_in;i++) if(idxn_in[i]>=0) this->idxn[count++] = idxn_in[i];
     64                                        count = 0;
     65                                        for(int i=0;i<m_in;i++){
     66                                                if(idxm_in[i]>=0){
     67                                                        for(int j=0;j<n_in;j++){
     68                                                                if(idxn_in[j]>=0){
     69                                                                        this->values[count++] = values_in[i*n_in+j];
     70                                                                }
     71                                                        }
     72                                                }
     73                                        }
     74                                }
     75                        }
     76                        else{
     77                                /*Set both as 0 to save time*/
     78                                this->m = 0;
     79                                this->n = 0;
     80                        }
     81                } /*}}}*/
     82                Bucket(int m_in,int* idxm_in,doubletype* values_in,InsMode mode_in){ /*{{{*/
     83                        this->Initialize();
     84
     85                        /*Check whether we have negative indices?*/
     86                        int numneg = 0;
     87                        for(int i=0;i<m_in;i++) if(idxm_in[i]<0) numneg++;
     88
     89                        this->type=VECTOR_BUCKET;
     90                        this->m=m_in-numneg;
     91                        this->mode=mode_in;
    4492                        if(this->m){
    45                                 this->idxm=xNew<int>(this->m);
    46                                 xMemCpy(this->idxm,idxmin,this->m);
    47                         }
    48                         if(this->n){
    49                                 this->idxn=xNew<int>(this->n);
    50                                 xMemCpy(this->idxn,idxnin,this->n);
    51                         }
    52                         if(this->m*this->n){
    53                                 this->values=xNew<doubletype>(this->n*this->m);
    54                                 xMemCpy(this->values,valuesin,this->n*this->m);
    55                         }
    56                 } /*}}}*/
    57                 Bucket(int min,int* idxmin,doubletype* valuesin,InsMode modein){ /*{{{*/
    58                         this->Initialize();
    59 
    60                         this->type=VECTOR_BUCKET;
    61                         this->m=min;
    62                         this->mode=modein;
    63                         if(this->m){
    64                                 this->idxm=xNew<int>(this->m);
    65                                 xMemCpy(this->idxm,idxmin,this->m);
    66 
     93                                this->idxm=xNew<int>(this->m);
    6794                                this->values=xNew<doubletype>(this->m);
    68                                 xMemCpy(this->values,valuesin,this->m);
     95
     96                                if(numneg==0){
     97                                        xMemCpy(this->idxm,idxm_in,this->m);
     98                                        xMemCpy(this->values,values_in,this->m);
     99                                }
     100                                else{
     101                                        int count = 0;
     102                                        for(int i=0;i<m_in;i++){
     103                                                if(idxm_in[i]>=0){
     104                                                        this->idxm[count]   = idxm_in[i];
     105                                                        this->values[count] = values_in[i];
     106                                                        count++;
     107                                                }
     108                                        }
     109                                }
    69110                        }
    70111                } /*}}}*/
    … …  
    88129                /*object virtual functions definitions:*/
    89130                void    Echo(){ /*{{{*/
    90                         _printf_("Bucket echo (cpu #: "<<IssmComm::GetRank()<<")\n");
     131                        _printf_("Bucket echo (cpu #"<<IssmComm::GetRank()<<")\n");
    91132                        _printf_("bucket type: " << type << "\n");
    92133                        _printf_("num rows: "<<this->m<<" num cols: "<<this->n << "\n");
    … …  
    95136                        int i,j;
    96137
    97                         _printf_("Bucket echo (cpu #: "<<IssmComm::GetRank()<<")\n");
     138                        _printf_("Bucket echo (cpu #"<<IssmComm::GetRank()<<")\n");
    98139                        _printf_("bucket type: " << type << "\n");
    99140                        _printf_("num rows: "<<this->m<<" num cols: "<<this->n << "\n");

  • issm/trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp

    r26110 r26131  
    3737        VecSetSizes(this->vector,m,M);
    3838        VecSetFromOptions(this->vector);
     39        VecSetOption(this->vector,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
    3940        VecSet(this->vector,0.0);
    4041

  • issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp

    r25710 r26131  
    111111        MatSetSizes(outmatrix,m,n,M,N);
    112112        MatSetFromOptions(outmatrix);
     113        MatSetOption(outmatrix,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
    113114
    114115        /*preallocation  according to type: */

  • issm/trunk-jpl/src/c/toolkits/petsc/patches/NewVec.cpp

    r17662 r26131  
    3535        VecSetSizes(vector,local_size,PETSC_DECIDE);
    3636        VecSetFromOptions(vector);
     37        VecSetOption(vector,VEC_IGNORE_NEGATIVE_INDICES,PETSC_TRUE);
    3738
    3839        return vector;
Note: See TracChangeset for help on using the changeset viewer.