Changeset 15688


Ignore:
Timestamp:
08/02/13 11:47:37 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: implemented FS condensed mini

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

Legend:

Unmodified
Added
Removed

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

    r15670 r15688  
    32683268                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
    32693269                        break;
    3270                 case P1P1Enum: case P1P1GLSEnum: case MINIcondensedEnum:
     3270                case P1P1Enum: case P1P1GLSEnum:
    32713271                        numnodes         = 12;
    32723272                        penta_node_ids   = xNew<int>(numnodes);
    … …  
    32853285                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+5];
    32863286                        break;
    3287                 case MINIEnum:
     3287                case MINIEnum: case MINIcondensedEnum:
    32883288                        numnodes         = 13;
    32893289                        penta_node_ids   = xNew<int>(numnodes);
    … …  
    57925792
    57935793        /*Fetch number of nodes and dof for this finite element*/
    5794         int vnumnodes = this->NumberofNodesVelocityFinal();
     5794        int vnumnodes = this->NumberofNodesVelocity();
    57955795        int pnumnodes = this->NumberofNodesPressure();
    57965796        int vnumdof   = vnumnodes*NDOF3;
    … …  
    74187418        Ke =new ElementMatrix(Ke1,Ke2);
    74197419
     7420        /*Condense if requested*/
     7421        if(this->element_type==MINIcondensedEnum){
     7422                int indices[3]={18,19,20};
     7423                Ke->StaticCondensation(3,&indices[0]);
     7424        }
     7425
    74207426        /*clean-up and return*/
    74217427        delete Ke1;
    … …  
    84418447        pe =new ElementVector(pe1,pe2,pe3);
    84428448
     8449        /*Condense if requested*/
     8450        if(this->element_type==MINIcondensedEnum){
     8451                int indices[3]={18,19,20};
     8452
     8453                this->element_type=MINIEnum;
     8454                ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
     8455                this->element_type=MINIcondensedEnum;
     8456
     8457                pe->StaticCondensation(Ke,3,&indices[0]);
     8458                delete Ke;
     8459        }
     8460
    84438461        /*clean-up and return*/
    84448462        delete pe1;
    … …  
    92509268
    92519269        /*Fetch number of nodes and dof for this finite element*/
    9252         int vnumnodes = this->NumberofNodesVelocityFinal();
     9270        int vnumnodes = this->NumberofNodesVelocity();
    92539271        int pnumnodes = this->NumberofNodesPressure();
    92549272        int vnumdof   = vnumnodes*NDOF3;
    … …  
    1009910117
    1010010118        /*Fetch number of nodes and dof for this finite element*/
    10101         int vnumnodes = this->NumberofNodesVelocityFinal();
     10119        int vnumnodes = this->NumberofNodesVelocity();
    1010210120        int pnumnodes = this->NumberofNodesPressure();
    1010310121        int vnumdof   = vnumnodes*NDOF3;

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

    r15656 r15688  
    19291929                case P1P1Enum:          return NUMNODESP1*2;
    19301930                case P1P1GLSEnum:       return NUMNODESP1*2;
    1931                 case MINIcondensedEnum: return NUMNODESP1*2;
     1931                case MINIcondensedEnum: return NUMNODESP1b+NUMNODESP1;
    19321932                case MINIEnum:          return NUMNODESP1b+NUMNODESP1;
    19331933                case TaylorHoodEnum:    return NUMNODESP2+NUMNODESP1;
    … …  
    19681968}
    19691969/*}}}*/
    1970 /*FUNCTION PentaRef::NumberofNodesVelocityFinal{{{*/
    1971 int PentaRef::NumberofNodesVelocityFinal(void){
    1972 
    1973         /*When static condensation is applied, the final number of nodes might be different*/
    1974 
    1975         switch(this->element_type){
    1976                 case P1P1Enum:          return NUMNODESP1;
    1977                 case P1P1GLSEnum:       return NUMNODESP1;
    1978                 case MINIcondensedEnum: return NUMNODESP1;
    1979                 case MINIEnum:          return NUMNODESP1b;
    1980                 case TaylorHoodEnum:    return NUMNODESP2;
    1981                 default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    1982         }
    1983 
    1984         return -1;
    1985 }
    1986 /*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/PentaRef.h

    r15654 r15688  
    7171                int  NumberofNodes(void);
    7272                int  NumberofNodesVelocity(void);
    73                 int  NumberofNodesVelocityFinal(void);
    7473                int  NumberofNodesPressure(void);
    7574};

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

    r15643 r15688  
    833833
    834834        /*intermediary: */
    835         int      i;
     835        int      i,M,N;
    836836        int      configuration_type;
    837837        Element *element = NULL;
    … …  
    846846
    847847        /*Display message*/
    848         if(VerboseModule()) _printf0_("   Generating matrices\n");
     848        if(VerboseModule()) _printf0_("   Generating matrices");
    849849
    850850        /*retrive parameters: */
    … …  
    877877        /*Allocate stiffness matrices and load vector*/
    878878        this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf);
     879
     880        /*Display size*/
     881        if(VerboseModule()){
     882                Kff->GetSize(&M,&N);
     883                _printf0_(" (Kff stiffness matrix size: "<<M<<" x "<<N<<")\n");
     884        }
    879885
    880886        /*Fill stiffness matrix from elements and loads */

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

    r15654 r15688  
    372372}
    373373/*}}}*/
     374/*FUNCTION ElementMatrix::StaticCondensation{{{*/
     375void ElementMatrix::StaticCondensation(int bsize,int* bindices){
     376        /*
     377         * | Kii  Kib | | Ui |    |Fi|
     378         * | Kbi  Kbb | | Ub |  = |Fb|
     379         *
     380         * Kii Ui + Kib Ub = Fi
     381         * Kbi Ui + Kbb Ub = Fb
     382         *
     383         * We want to remove Ub from the equation:
     384         *
     385         * Kii Ui + Kib inv(Kbb) (Fb - Kbi Ui) = Fi
     386         *
     387         * which gives:
     388         *
     389         * (Kii - Kib inv(Kbb) Kbi) Ui = Fi - Kib inv(Kbb) Fb
     390         */
     391
     392        /*Checks in debugging mode*/
     393        _assert_(this->nrows==this->ncols && bsize>0 && bsize<this->ncols && this->values);
     394
     395        /*Intermediaries*/
     396        int         counter,i,j,isize;
     397        IssmDouble *Kii         = NULL;
     398        IssmDouble *Kib         = NULL;
     399        IssmDouble *Kbi         = NULL;
     400        IssmDouble *Kbb         = NULL;
     401        IssmDouble *Kbbinv      = NULL;
     402        IssmDouble *Ktemp       = NULL;
     403        int        *iindices    = NULL;
     404        bool        flag;
     405
     406        /*Get new sizes and indices*/
     407        isize    = this->nrows - bsize;
     408        iindices = xNew<int>(isize);
     409        counter  = 0;
     410        for(i=0;i<this->nrows;i++){
     411                flag = true;
     412                for(j=0;j<bsize;j++){
     413                        if(i==bindices[j]){
     414                                flag = false;
     415                                break;
     416                        }
     417                }
     418                if(flag){
     419                        _assert_(counter<isize);
     420                        iindices[counter++] = i;
     421                }
     422        }
     423        _assert_(counter == isize);
     424
     425        /*Get submatrices*/
     426        Kii = xNew<IssmDouble>(isize*isize);
     427        Kib = xNew<IssmDouble>(isize*bsize);
     428        Kbi = xNew<IssmDouble>(bsize*isize);
     429        Kbb = xNew<IssmDouble>(bsize*bsize);
     430        for(i=0;i<isize;i++) for(j=0;j<isize;j++) Kii[i*isize+j] = this->values[iindices[i]*this->ncols + iindices[j]];
     431        for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = this->values[iindices[i]*this->ncols + bindices[j]];
     432        for(i=0;i<bsize;i++) for(j=0;j<isize;j++) Kbi[i*isize+j] = this->values[bindices[i]*this->ncols + iindices[j]];
     433        for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = this->values[bindices[i]*this->ncols + bindices[j]];
     434
     435        /*Invert Kbb*/
     436        Kbbinv = xNew<IssmDouble>(bsize*bsize);
     437        switch(bsize){
     438                case 1:
     439                        Kbbinv[0] = 1./Kbb[0];
     440                        break;
     441                case 2:
     442                        Matrix2x2Invert(Kbbinv,Kbb);
     443                        break;
     444                case 3:
     445                        Matrix3x3Invert(Kbbinv,Kbb);
     446                        break;
     447                default:
     448                        MatrixInverse(Kbbinv,bsize,bsize,NULL,0,NULL);
     449                        break;
     450        }
     451
     452        /*Calculate  Kib inv(Kbb) Kbi*/
     453        Ktemp = xNew<IssmDouble>(isize*isize);
     454        TripleMultiply(Kib,isize,bsize,0, Kbbinv,bsize,bsize,0, Kbi,bsize,isize,0, Ktemp,0);
     455
     456        /*New Ke*/
     457        for(i=0;i<isize*isize;i++) Ktemp[i] = Kii[i] - Ktemp[i];
     458
     459        /*Update matrix values*/
     460        for(i=0;i<this->nrows*this->ncols;i++) this->values[i]=0.;
     461        for(i=0;i<isize;i++){
     462                for(j=0;j<isize;j++){
     463                        this->values[iindices[i]*this->ncols + iindices[j]] = Ktemp[i*isize+j];
     464                }
     465        }
     466
     467        /*Clean up and return*/
     468        xDelete<IssmDouble>(Kii);
     469        xDelete<IssmDouble>(Kib);
     470        xDelete<IssmDouble>(Kbi);
     471        xDelete<IssmDouble>(Kbb);
     472        xDelete<IssmDouble>(Kbbinv);
     473        xDelete<IssmDouble>(Ktemp);
     474        xDelete<int>(iindices);
     475        return;
     476}
     477/*}}}*/
    374478/*FUNCTION ElementMatrix::Echo{{{*/
    375479void ElementMatrix::Echo(void){

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

    r15067 r15688  
    1010
    1111/*Headers:*/
    12 /*{{{*/
    1312#include "../../datastructures/datastructures.h"
    1413#include "../../toolkits/toolkits.h"
    … …  
    1716template <class doublematrix> class Matrix;
    1817class Parameters;
    19 /*}}}*/
    2018
    2119class ElementMatrix{
    … …  
    5149                int*     col_sglobaldoflist;
    5250
    53                 /*ElementMatrix constructors, destructors {{{*/
     51                /*ElementMatrix constructors, destructors*/
    5452                ElementMatrix();
    5553                ElementMatrix(ElementMatrix* Ke);
    … …  
    5856                ElementMatrix(Node** nodes,int numnodes,Parameters* parameters,int approximation=NoneApproximationEnum);
    5957                ~ElementMatrix();
    60                 /*}}}*/
    61                 /*ElementMatrix specific routines {{{*/
     58
     59                /*ElementMatrix specific routines*/
    6260                void AddToGlobal(Matrix<IssmDouble>* Kff, Matrix<IssmDouble>* Kfs);
    6361                void AddToGlobal(Matrix<IssmDouble>* Jff);
    6462                void Echo(void);
    6563                void CheckConsistency(void);
     64                void StaticCondensation(int numindices,int* indices);
    6665                void Transpose(void);
    6766                void Init(ElementMatrix* Ke);
    6867                void SetDiag(IssmDouble scalar);
    69                 /*}}}*/
    7068};
    7169#endif //#ifndef _ELEMENT_MATRIX_H_

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

    r15104 r15688  
    274274}
    275275/*}}}*/
     276/*FUNCTION ElementVector::StaticCondensation{{{*/
     277void ElementVector::StaticCondensation(ElementMatrix* Ke,int bsize,int* bindices){
     278        /*
     279         * | Kii  Kib | | Ui |    |Fi|
     280         * | Kbi  Kbb | | Ub |  = |Fb|
     281         *
     282         * Kii Ui + Kib Ub = Fi
     283         * Kbi Ui + Kbb Ub = Fb
     284         *
     285         * We want to remove Ub from the equation:
     286         *
     287         * Kii Ui + Kib inv(Kbb) (Fb - Kbi Ui) = Fi
     288         *
     289         * which gives:
     290         *
     291         * (Kii - Kib inv(Kbb) Kbi) Ui = Fi - Kib inv(Kbb) Fb
     292         */
     293
     294        /*Checks in debugging mode*/
     295        _assert_(bsize>0 && bsize<this->nrows && this->values && Ke);
     296        _assert_(Ke->nrows==Ke->ncols && this->nrows==Ke->nrows);
     297
     298        /*Intermediaries*/
     299        int         counter,i,j,isize;
     300        IssmDouble *Fb          = NULL;
     301        IssmDouble *Fi          = NULL;
     302        IssmDouble *Kib         = NULL;
     303        IssmDouble *Kbb         = NULL;
     304        IssmDouble *Kbbinv      = NULL;
     305        IssmDouble *Ftemp       = NULL;
     306        int        *iindices    = NULL;
     307        bool        flag;
     308
     309        /*Get new sizes and indices*/
     310        isize    = this->nrows - bsize;
     311        iindices = xNew<int>(isize);
     312        counter  = 0;
     313        for(i=0;i<this->nrows;i++){
     314                flag = true;
     315                for(j=0;j<bsize;j++){
     316                        if(i==bindices[j]){
     317                                flag = false;
     318                                break;
     319                        }
     320                }
     321                if(flag){
     322                        _assert_(counter<isize);
     323                        iindices[counter++] = i;
     324                }
     325        }
     326        _assert_(counter == isize);
     327
     328        /*Get submatrices*/
     329        Kib = xNew<IssmDouble>(isize*bsize);
     330        Kbb = xNew<IssmDouble>(bsize*bsize);
     331        Fb  = xNew<IssmDouble>(bsize);
     332        Fi  = xNew<IssmDouble>(isize);
     333        for(i=0;i<isize;i++) for(j=0;j<bsize;j++) Kib[i*bsize+j] = Ke->values[iindices[i]*Ke->ncols + bindices[j]];
     334        for(i=0;i<bsize;i++) for(j=0;j<bsize;j++) Kbb[i*bsize+j] = Ke->values[bindices[i]*Ke->ncols + bindices[j]];
     335        for(i=0;i<bsize;i++) Fb[i] = this->values[bindices[i]];
     336        for(i=0;i<isize;i++) Fi[i] = this->values[iindices[i]];
     337
     338        /*Invert Kbb*/
     339        Kbbinv = xNew<IssmDouble>(bsize*bsize);
     340        switch(bsize){
     341                case 1:
     342                        Kbbinv[0] = 1./Kbb[0];
     343                        break;
     344                case 2:
     345                        Matrix2x2Invert(Kbbinv,Kbb);
     346                        break;
     347                case 3:
     348                        Matrix3x3Invert(Kbbinv,Kbb);
     349                        break;
     350                default:
     351                        MatrixInverse(Kbbinv,bsize,bsize,NULL,0,NULL);
     352                        break;
     353        }
     354
     355        /*Calculate  Kib inv(Kbb) Fb*/
     356        Ftemp = xNew<IssmDouble>(isize);
     357        TripleMultiply(Kib,isize,bsize,0, Kbbinv,bsize,bsize,0, Fb,bsize,1,0, Ftemp,0);
     358
     359        /*New Pe*/
     360        for(i=0;i<isize;i++) Ftemp[i] = Fi[i] - Ftemp[i];
     361
     362        /*Update matrix values*/
     363        for(i=0;i<this->nrows;i++) this->values[i]=0.;
     364        for(i=0;i<isize;i++){
     365                this->values[iindices[i]] = Ftemp[i];
     366        }
     367
     368        /*Clean up and return*/
     369        xDelete<IssmDouble>(Kib);
     370        xDelete<IssmDouble>(Kbb);
     371        xDelete<IssmDouble>(Kbbinv);
     372        xDelete<IssmDouble>(Fb);
     373        xDelete<IssmDouble>(Fi);
     374        xDelete<IssmDouble>(Ftemp);
     375        xDelete<int>(iindices);
     376        return;
     377}
     378/*}}}*/

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

    r15067 r15688  
    1010
    1111/*Headers:*/
    12 /*{{{*/
    1312#include "../../datastructures/datastructures.h"
    1413#include "../../toolkits/toolkits.h"
    … …  
    1716template <class doubletype> class Vector;
    1817class Parameters;
    19 /*}}}*/
     18class ElementMatrix;
    2019
    2120class ElementVector{
    2221
    2322        public:
    24 
    25                 int      nrows;
    26                 IssmDouble*  values;
     23                int         nrows;
     24                IssmDouble* values;
    2725
    2826                //gset
    29                 int*     gglobaldoflist;
     27                int* gglobaldoflist;
    3028
    3129                //fset
    32                 int      fsize;
    33                 int*     flocaldoflist;
    34                 int*     fglobaldoflist;
     30                int  fsize;
     31                int* flocaldoflist;
     32                int* fglobaldoflist;
    3533
    36                 /*ElementVector constructors, destructors {{{*/
     34                /*ElementVector constructors, destructors*/
    3735                ElementVector();
    3836                ElementVector(ElementVector* pe1,ElementVector* pe2);
    … …  
    4038                ElementVector(Node** nodes,int numnodes,Parameters* parameters,int approximation=NoneApproximationEnum);
    4139                ~ElementVector();
    42                 /*}}}*/
    43                 /*ElementVector specific routines {{{*/
     40
     41                /*ElementVector specific routines*/
    4442                void AddToGlobal(Vector<IssmDouble>* pf);
    4543                void InsertIntoGlobal(Vector<IssmDouble>* pf);
    … …  
    4846                void Init(ElementVector* pe);
    4947                void SetValue(IssmDouble scalar);
    50                 /*}}}*/
     48                void StaticCondensation(ElementMatrix* Ke,int numindices,int* indices);
    5149};
    5250#endif //#ifndef _ELEMENT_VECTOR_H_

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp

    r15658 r15688  
    130130                case MINIcondensedEnum:
    131131                        _assert_(approximation==FSApproximationEnum);
    132                         /*P1 velocity (bubble statically condensed)*/
    133                         for(i=0;i<iomodel->numberofvertices;i++){
    134                                 if(iomodel->my_vertices[i]){
    135                                         nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
    136                                 }
    137                         }
    138                         /*P1 pressure*/
    139                         for(i=0;i<iomodel->numberofvertices;i++){
    140                                 if(iomodel->my_vertices[i]){
    141                                         nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum));
     132                        /*P1+ velocity (bubble statically condensed)*/
     133                        for(i=0;i<iomodel->numberofvertices;i++){
     134                                if(iomodel->my_vertices[i]){
     135                                        nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum));
     136                                }
     137                        }
     138                        for(i=0;i<iomodel->numberofelements;i++){
     139                                if(iomodel->my_elements[i]){
     140                                        node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,FSvelocityEnum);
     141                                        node->Deactivate();
     142                                        nodes->AddObject(node);
     143                                }
     144                        }
     145                        /*P1 pressure*/
     146                        for(i=0;i<iomodel->numberofvertices;i++){
     147                                if(iomodel->my_vertices[i]){
     148                                        nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,i,iomodel,analysis,FSpressureEnum));
    142149                                }
    143150                        }

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r15657 r15688  
    9898                                case 0 : finiteelement = P1Enum;       break;//P1P1
    9999                                case 1 : finiteelement = P1Enum;       break;//P1P1GSL
    100                                 case 2 : finiteelement = P1Enum;      break;//MINIcondensed
     100                                case 2 : finiteelement = P1bubbleEnum; break;//MINIcondensed
    101101                                case 3 : finiteelement = P1bubbleEnum; break;//MINI
    102102                                case 4 : finiteelement = P2Enum;       break;//TaylorHood (P2P1)

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp

    r15656 r15688  
    7070                        }
    7171                }
    72 
    7372                iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
    7473                                        MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum);

  • issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

    r15104 r15688  
    2121        int         idima,idimb,idimc,idimd;
    2222        IssmDouble *dtemp;
     23        _assert_(a && b && c && d);
    2324
    2425        /*  set up dimensions for triple product  */
    25 
    2626        if (!itrna){
    2727                idima=nrowa;
    … …  
    8383        int i,j,k,ipta,iptb,iptc;
    8484        int nrowc,ncolc,iinca,jinca,iincb,jincb,ntrma,ntrmb,nterm;
     85
     86        _assert_(a && b && c);
    8587
    8688        /*  set up dimensions and increments for matrix a  */
    … …  
    326328        /*Compute determinant*/
    327329        Matrix2x2Determinant(&det,A);
    328         if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller that machine epsilon");
     330        if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
    329331
    330332        /*Compute invert*/
    … …  
    351353        /*Compute determinant*/
    352354        Matrix3x3Determinant(&det,A);
    353         if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller that machine epsilon");
     355        if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
    354356
    355357        /*Compute invert*/
Note: See TracChangeset for help on using the changeset viewer.