Changeset 15694


Ignore:
Timestamp:
08/02/13 15:05:10 (12 years ago)
Author:
Mathieu Morlighem
Message:

NEW: extending P1bubble and P1bubblecondensed to SSA

Location:
issm/trunk-jpl/src
Files:
14 edited

Legend:

Unmodified
Added
Removed

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

    r15689 r15694  
    456456        }
    457457
    458         /*Add to global matrix*/
    459458        if(Ke){
     459                /*Condense if requested*/
     460                if(this->element_type==MINIcondensedEnum){
     461                        int indices[3]={18,19,20};
     462                        Ke->StaticCondensation(3,&indices[0]);
     463                }
     464
     465                /*Add to global matrix*/
    460466                Ke->AddToGlobal(Kff,Kfs);
    461467                delete Ke;
    … …  
    592598        }
    593599
    594         /*Add to global Vector*/
    595600        if(pe){
     601
     602                /*StaticCondensation if requested*/
     603                if(this->element_type==MINIcondensedEnum){
     604                        int indices[3]={18,19,20};
     605
     606                        this->element_type=MINIEnum;
     607                        ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
     608                        this->element_type=MINIcondensedEnum;
     609
     610                        pe->StaticCondensation(Ke,3,&indices[0]);
     611                        delete Ke;
     612                }
     613
     614                /*Add to global Vector*/
    596615                pe->AddToGlobal(pf);
    597616                delete pe;
    … …  
    11651184                GaussPenta* gauss=new GaussPenta();
    11661185                for(int iv=0;iv<this->NumberofNodes();iv++){
    1167                         gauss->GaussNode(numnodes,iv);
     1186                        gauss->GaussNode(this->element_type,iv);
    11681187                        input->GetInputValue(&pvalue[iv],gauss);
    11691188                }
    … …  
    11881207        GaussPenta* gauss=new GaussPenta();
    11891208        for (int iv=0;iv<this->NumberofNodes();iv++){
    1190                 gauss->GaussNode(numnodes,iv);
     1209                gauss->GaussNode(this->element_type,iv);
    11911210                input->GetInputValue(&pvalue[iv],gauss);
    11921211        }
    … …  
    74297448        Ke =new ElementMatrix(Ke1,Ke2);
    74307449
    7431         /*Condense if requested*/
    7432         if(this->element_type==MINIcondensedEnum){
    7433                 int indices[3]={18,19,20};
    7434                 Ke->StaticCondensation(3,&indices[0]);
    7435         }
    7436 
    74377450        /*clean-up and return*/
    74387451        delete Ke1;
    … …  
    84588471        pe =new ElementVector(pe1,pe2,pe3);
    84598472
    8460         /*Condense if requested*/
    8461         if(this->element_type==MINIcondensedEnum){
    8462                 int indices[3]={18,19,20};
    8463 
    8464                 this->element_type=MINIEnum;
    8465                 ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
    8466                 this->element_type=MINIcondensedEnum;
    8467 
    8468                 pe->StaticCondensation(Ke,3,&indices[0]);
    8469                 delete Ke;
    8470         }
    8471 
    84728473        /*clean-up and return*/
    84738474        delete pe1;
    … …  
    91829183        GaussPenta* gauss=new GaussPenta();
    91839184        for(int i=0;i<numnodes;i++){
    9184                 gauss->GaussNode(numnodes,i);
     9185                gauss->GaussNode(this->element_type,i);
    91859186
    91869187                /*Recover vx and vy*/
    … …  
    93019302        gauss = new GaussPenta();
    93029303        for(int i=0;i<vnumnodes;i++){
    9303                 gauss->GaussNode(vnumnodes,i);
     9304                gauss->GaussNode(this->VelocityInterpolation(),i);
    93049305
    93059306                vx_input->GetInputValue(&vx,gauss);
    … …  
    93119312        }
    93129313        for(int i=0;i<pnumnodes;i++){
    9313                 gauss->GaussNode(pnumnodes,i);
     9314                gauss->GaussNode(this->PressureInterpolation(),i);
    93149315
    93159316                p_input ->GetInputValue(&p ,gauss);

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

    r15689 r15694  
    19681968}
    19691969/*}}}*/
     1970/*FUNCTION PentaRef::VelocityInterpolation{{{*/
     1971int PentaRef::VelocityInterpolation(void){
     1972
     1973        switch(this->element_type){
     1974                case P1P1Enum:          return P1Enum;
     1975                case P1P1GLSEnum:       return P1Enum;
     1976                case MINIcondensedEnum: return P1bubbleEnum;
     1977                case MINIEnum:          return P1bubbleEnum;
     1978                case TaylorHoodEnum:    return P2Enum;
     1979                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1980        }
     1981
     1982        return -1;
     1983}
     1984/*}}}*/
     1985/*FUNCTION PentaRef::PressureInterpolation{{{*/
     1986int PentaRef::PressureInterpolation(void){
     1987
     1988        switch(this->element_type){
     1989                case P1P1Enum:          return P1Enum;
     1990                case P1P1GLSEnum:       return P1Enum;
     1991                case MINIcondensedEnum: return P1Enum;
     1992                case MINIEnum:          return P1Enum;
     1993                case TaylorHoodEnum:    return P1Enum;
     1994                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
     1995        }
     1996
     1997        return -1;
     1998}
     1999/*}}}*/

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

    r15688 r15694  
    7272                int  NumberofNodesVelocity(void);
    7373                int  NumberofNodesPressure(void);
     74                int  VelocityInterpolation(void);
     75                int  PressureInterpolation(void);
    7476};
    7577#endif

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

    r15689 r15694  
    231231        }
    232232
    233         /*Add to global matrix*/
    234233        if(Ke){
     234                /*Static condensation if requested*/
     235                if(this->element_type==P1bubblecondensedEnum){
     236                        int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     237                        int offset = 0;
     238                        for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     239                        int* indices=xNew<int>(size);
     240                        for(int i=0;i<size;i++) indices[i] = offset+i;
     241                        printarray(indices,1,size);
     242                        Ke->StaticCondensation(size,indices);
     243                        xDelete<int>(indices);
     244                }
     245
     246                /*Add to global matrix*/
    235247                Ke->AddToGlobal(Kff,Kfs);
    236248                delete Ke;
    … …  
    346358        /*Add to global Vector*/
    347359        if(pe){
     360                /*Static condensation if requested*/
     361                if(this->element_type==P1bubblecondensedEnum){
     362                        int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     363                        int offset = 0;
     364                        for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
     365                        int* indices=xNew<int>(size);
     366                        for(int i=0;i<size;i++) indices[i] = offset+i;
     367
     368                        this->element_type=P1bubbleEnum;
     369                        ElementMatrix* Ke = CreateKMatrixDiagnosticSSA();
     370                        this->element_type=P1bubblecondensedEnum;
     371
     372                        pe->StaticCondensation(Ke,size,indices);
     373                        xDelete<int>(indices);
     374                        delete Ke;
     375                }
     376
    348377                pe->AddToGlobal(pf);
    349378                delete pe;
    … …  
    11291158                GaussTria* gauss=new GaussTria();
    11301159                for(int iv=0;iv<this->NumberofNodes();iv++){
    1131                         gauss->GaussNode(numnodes,iv);
     1160                        gauss->GaussNode(this->element_type,iv);
    11321161                        input->GetInputValue(&pvalue[iv],gauss);
    11331162                }
    … …  
    11521181        GaussTria* gauss=new GaussTria();
    11531182        for (int iv=0;iv<this->NumberofNodes();iv++){
    1154                 gauss->GaussNode(numnodes,iv);
     1183                gauss->GaussNode(this->element_type,iv);
    11551184                input->GetInputValue(&pvalue[iv],gauss);
    11561185        }
    … …  
    33593388        gauss=new GaussTria();
    33603389        for(int i=0;i<numnodes;i++){
    3361                 gauss->GaussNode(numnodes,i);
     3390                gauss->GaussNode(this->element_type,i);
    33623391
    33633392                /*Recover vx and vy*/

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

    r15689 r15694  
    444444                        /*bubble*/
    445445                        basis[3]=27.*gauss->coord1*gauss->coord2*gauss->coord3;
     446                        return;
    446447                case P2Enum:
    447448                        /*Corner nodes*/
    … …  
    560561                        dbasis[NUMNODESP1b*1+2] = SQRT3/3.;
    561562                        /*Nodal function 4*/
    562                         dbasis[NUMNODESP1b*0+2] = 27.*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);;
    563                         dbasis[NUMNODESP1b*1+2] = 27.*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
     563                        dbasis[NUMNODESP1b*0+3] = 27.*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);;
     564                        dbasis[NUMNODESP1b*1+3] = 27.*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
    564565                        return;
    565566                case P2Enum:

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp

    r15654 r15694  
    88#include "../../shared/Exceptions/exceptions.h"
    99#include "../../shared/MemOps/MemOps.h"
     10#include "../../shared/Enum/Enum.h"
    1011#include "../../shared/Numerics/GaussPoints.h"
    1112#include "../../shared/Numerics/constants.h"
    … …  
    397398/*}}}*/
    398399/*FUNCTION GaussPenta::GaussNode{{{*/
    399 void GaussPenta::GaussNode(int numnodes,int iv){
     400void GaussPenta::GaussNode(int finiteelement,int iv){
    400401
    401402        /*in debugging mode: check that the default constructor has been called*/
    … …  
    403404
    404405        /*update static arrays*/
    405         switch(numnodes){
    406                 case 6: //P1 Lagrange element
     406        switch(finiteelement){
     407                case P1Enum: case P1DGEnum:
    407408                        switch(iv){
    408409                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
    … …  
    415416                        }
    416417                        break;
    417                 case 9: //P1xP2 Lagrange element
     418                case P1xP2Enum:
    418419                        switch(iv){
    419420                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
    … …  
    430431                        }
    431432                        break;
    432                 case 12: //P2xP1 Lagrange element
     433                case P2xP1Enum:
    433434                        switch(iv){
    434435                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
    … …  
    448449                        }
    449450                        break;
    450                 case 7: //P1+ Lagrange element
     451                case P1bubbleEnum:  case P1bubblecondensedEnum:
     452                        switch(iv){
     453                                case 0: coord1=1.;    coord2=0.;    coord3=0.;    coord4=-1.; break;
     454                                case 1: coord1=0.;    coord2=1.;    coord3=0.;    coord4=-1.; break;
     455                                case 2: coord1=0.;    coord2=0.;    coord3=1.;    coord4=-1.; break;
     456                                case 3: coord1=1.;    coord2=0.;    coord3=0.;    coord4=+1.; break;
     457                                case 4: coord1=0.;    coord2=1.;    coord3=0.;    coord4=+1.; break;
     458                                case 5: coord1=0.;    coord2=0.;    coord3=1.;    coord4=+1.; break;
     459                                case 6: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.;  break;
     460                                default: _error_("node index should be in [0 5]");
     461                        }
     462                        break;
     463                case P2Enum:
    451464                        switch(iv){
    452465                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
    … …  
    457470                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
    458471
    459                                 case 6: coord1=0.; coord2=0.; coord3=0.; coord4=0.;break;
    460                                 default: _error_("node index should be in [0 5]");
    461                         }
    462                         break;
    463                 case 15: //P2 Lagrange element
    464                         switch(iv){
    465                                 case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
    466                                 case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
    467                                 case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
    468                                 case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
    469                                 case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
    470                                 case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
    471 
    472472                                case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
    473473                                case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
    … …  
    483483                        }
    484484                        break;
    485                 default: _error_("Number of nodes "<<numnodes<<" not supported");
     485                default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
    486486        }
    487487

  • issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h

    r15625 r15694  
    4444                void GaussPoint(int ig);
    4545                void GaussVertex(int iv);
    46                 void GaussNode(int numnodes,int iv);
     46                void GaussNode(int finitelement,int iv);
    4747                void GaussFaceTria(int index1, int index2, int index3, int order);
    4848                void GaussCenter(void);

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.cpp

    r15538 r15694  
    442442/*}}}*/
    443443/*FUNCTION GaussTria::GaussNode{{{*/
    444 void GaussTria::GaussNode(int numnodes,int iv){
     444void GaussTria::GaussNode(int finiteelement,int iv){
    445445
    446446        /*in debugging mode: check that the default constructor has been called*/
    … …  
    448448
    449449        /*update static arrays*/
    450         switch(numnodes){
    451                 case 3: //P1 Lagrange element
     450        switch(finiteelement){
     451                case P1Enum: case P1DGEnum:
    452452                        switch(iv){
    453453                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
    … …  
    457457                        }
    458458                        break;
    459                 case 6: //P2 Lagrange element
     459                case P1bubbleEnum: case P1bubblecondensedEnum:
     460                        switch(iv){
     461                                case 0: coord1=1.;    coord2=0.;    coord3=0.;    break;
     462                                case 1: coord1=0.;    coord2=1.;    coord3=0.;    break;
     463                                case 2: coord1=0.;    coord2=0.;    coord3=1.;    break;
     464                                case 3: coord1=1./3.; coord2=1./3.; coord3=1./3.; break;
     465                                default: _error_("node index should be in [0 3]");
     466                        }
     467                        break;
     468                case P2Enum:
    460469                        switch(iv){
    461470                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
    … …  
    468477                        }
    469478                        break;
    470                 default: _error_("supported number of nodes are 3 and 6");
     479                default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
    471480        }
    472481

  • issm/trunk-jpl/src/c/classes/gauss/GaussTria.h

    r15517 r15694  
    4141                void GaussPoint(int ig);
    4242                void GaussVertex(int iv);
    43                 void GaussNode(int numnodes,int iv);
     43                void GaussNode(int finitelement,int iv);
    4444                void GaussCenter(void);
    4545                void GaussEdgeCenter(int index1,int index2);

  • issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp

    r15689 r15694  
    4545
    4646        switch(finite_element){
    47                 case P1Enum:
     47                case P1Enum: 
    4848                        /*Nothing else to do*/
    4949                        break;
    5050                case P1bubbleEnum:
    51                         if(iomodel->dim!=3) _error_("3d is the only supported dimension");
    52                         elementnbv = 6;
     51                        switch(iomodel->dim){
     52                                case 2: elementnbv = 3; break;
     53                                case 3: elementnbv = 6; break;
     54                                default: _error_("3d is the only supported dimension");
     55                        }
     56                        break;
     57                case P1bubblecondensedEnum:
     58                        /*Nothing else to do*/
    5359                        break;
    5460                case P1xP2Enum:

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

    r15688 r15694  
    7474                        iomodel->Constant(&temp,FlowequationFeSSAEnum);
    7575                        switch(temp){
    76                                 case 0 : finiteelement = P1Enum; break;
    77                                 case 1 : finiteelement = P2Enum; break;
     76                                case 0 : finiteelement = P1Enum;                break;
     77                                case 1 : finiteelement = P2Enum;                break;
     78                                case 2 : finiteelement = P1bubblecondensedEnum; break;
     79                                case 3 : finiteelement = P1bubbleEnum;          break;
    7880                                default: _error_("finite element "<<temp<<" not supported");
    7981                        }
    … …  
    8587                        iomodel->Constant(&temp,FlowequationFeHOEnum);
    8688                        switch(temp){
    87                                 case 0 : finiteelement = P1Enum;    break;
    88                                 case 1 : finiteelement = P1xP2Enum; break;
    89                                 case 2 : finiteelement = P2xP1Enum; break;
    90                                 case 3 : finiteelement = P2Enum;    break;
     89                                case 0 : finiteelement = P1Enum;                break;
     90                                case 1 : finiteelement = P1xP2Enum;             break;
     91                                case 2 : finiteelement = P2xP1Enum;             break;
     92                                case 3 : finiteelement = P2Enum;                break;
     93                                case 4 : finiteelement = P1bubblecondensedEnum; break;
     94                                case 5 : finiteelement = P1bubbleEnum;          break;
    9195                                default: _error_("finite element "<<temp<<" not supported");
    9296                        }

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

    r15688 r15694  
    3838                        iomodel->Constant(&temp,FlowequationFeSSAEnum);
    3939                        switch(temp){
    40                                 case 0 : finiteelement = P1Enum; break;
    41                                 case 1 : finiteelement = P2Enum; break;
     40                                case 0 : finiteelement = P1Enum;                break;
     41                                case 1 : finiteelement = P2Enum;                break;
     42                                case 2 : finiteelement = P1bubblecondensedEnum; break;
     43                                case 3 : finiteelement = P1bubbleEnum;          break;
    4244                                default: _error_("finite element "<<temp<<" not supported");
    4345                        }
    … …  
    5153                        iomodel->Constant(&temp,FlowequationFeHOEnum);
    5254                        switch(temp){
    53                                 case 0 : finiteelement = P1Enum;    break;
    54                                 case 1 : finiteelement = P1xP2Enum; break;
    55                                 case 2 : finiteelement = P2xP1Enum; break;
    56                                 case 3 : finiteelement = P2Enum;    break;
     55                                case 0 : finiteelement = P1Enum;                break;
     56                                case 1 : finiteelement = P1xP2Enum;             break;
     57                                case 2 : finiteelement = P2xP1Enum;             break;
     58                                case 3 : finiteelement = P2Enum;                break;
     59                                case 4 : finiteelement = P1bubblecondensedEnum; break;
     60                                case 5 : finiteelement = P1bubbleEnum;          break;
    5761                                default: _error_("finite element "<<temp<<" not supported");
    5862                        }

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

    r15656 r15694  
    4545                        iomodel->Constant(&temp,FlowequationFeSSAEnum);
    4646                        switch(temp){
    47                                 case 0 : finiteelement = P1Enum; break;
    48                                 case 1 : finiteelement = P2Enum; break;
     47                                case 0 : finiteelement = P1Enum;                break;
     48                                case 1 : finiteelement = P2Enum;                break;
     49                                case 2 : finiteelement = P1bubblecondensedEnum; break;
     50                                case 3 : finiteelement = P1bubbleEnum;          break;
    4951                                default: _error_("finite element "<<temp<<" not supported");
    5052                        }
    … …  
    5658                        iomodel->Constant(&temp,FlowequationFeHOEnum);
    5759                        switch(temp){
    58                                 case 0 : finiteelement = P1Enum;    break;
    59                                 case 1 : finiteelement = P1xP2Enum; break;
    60                                 case 2 : finiteelement = P2xP1Enum; break;
    61                                 case 3 : finiteelement = P2Enum;    break;
     60                                case 0 : finiteelement = P1Enum;                break;
     61                                case 1 : finiteelement = P1xP2Enum;             break;
     62                                case 2 : finiteelement = P2xP1Enum;             break;
     63                                case 3 : finiteelement = P2Enum;                break;
     64                                case 4 : finiteelement = P1bubblecondensedEnum; break;
     65                                case 5 : finiteelement = P1bubbleEnum;          break;
    6266                                default: _error_("finite element "<<temp<<" not supported");
    6367                        }

  • issm/trunk-jpl/src/m/classes/flowequation.m

    r15689 r15694  
    8181                                md = checkfield(md,'flowequation.isHO','numel',[1],'values',[0 1]);
    8282                                md = checkfield(md,'flowequation.isFS','numel',[1],'values',[0 1]);
    83                                 md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0 1]);
    84                                 md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0:3]);
     83                                md = checkfield(md,'flowequation.fe_SSA','numel',[1],'values',[0:3]);
     84                                md = checkfield(md,'flowequation.fe_HO','numel',[1],'values',[0:5]);
    8585                                md = checkfield(md,'flowequation.fe_FS','numel',[1],'values',[0:4]);
    8686                                md = checkfield(md,'flowequation.borderSSA','size',[md.mesh.numberofvertices 1],'values',[0 1]);
    … …  
    115115                        fielddisplay(obj,'isHO','is the Higher-Order (HO) approximation used ?');
    116116                        fielddisplay(obj,'isFS','are the Full-FS (FS) equations used ?');
    117                         fielddisplay(obj,'fe_SSA','Finite Element for SSA  0: Lagrange P1 (linear), 1: Lagrange P2 (quadratic)');
    118                         fielddisplay(obj,'fe_HO', 'Finite Element for HO   0: P1xP1, 1: P1xP2, 2: P2xP1, 3: P2xP2');
     117                        fielddisplay(obj,'fe_SSA','Finite Element for SSA  0: Lagrange P1 (linear), 1: Lagrange P2 (quadratic), 2: P1+ condensed, 3: P1+');
     118                        fielddisplay(obj,'fe_HO', 'Finite Element for HO   0: P1xP1, 1: P1xP2, 2: P2xP1, 3: P2xP2, 4: P1+ condensed, 5: P1+');
    119119                        fielddisplay(obj,'fe_FS', 'Finite Element for FS   0: P1P1 (debugging), 1: P1P1GSL (under dev), 2: MINI condensed, 3: MINI, 4: P2P1 (Taylor-Hood)');
    120120                        fielddisplay(obj,'vertex_equation','flow equation for each vertex');
Note: See TracChangeset for help on using the changeset viewer.