Changeset 18474


Ignore:
Timestamp:
09/02/14 07:40:03 (11 years ago)
Author:
seroussi
Message:

NEW: starting to implement LACR

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

Legend:

Unmodified
Added
Removed

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

    r18293 r18474  
    112112        /*XTH LATH parameters*/
    113113        iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
    114         if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum){
     114        if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
    115115                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
    116116                parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
    … …  
    247247        /*LATH parameters*/
    248248        iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
    249         if(fe_FS==LATaylorHoodEnum){
     249        if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
    250250                iomodel->FetchDataToInput(elements,PressureEnum,0.);
    251251                InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
    … …  
    428428                        /*Deduce velocity interpolation from finite element*/
    429429                        switch(finiteelement){
    430                                 case P1P1Enum            : finiteelement = P1Enum;       break;
    431                                 case P1P1GLSEnum         : finiteelement = P1Enum;       break;
    432                                 case MINIcondensedEnum   : finiteelement = P1bubbleEnum; break;
    433                                 case MINIEnum            : finiteelement = P1bubbleEnum; break;
    434                                 case TaylorHoodEnum      : finiteelement = P2Enum;       break;
    435                                 case XTaylorHoodEnum     : finiteelement = P2Enum;       break;
    436                                 case LATaylorHoodEnum    : finiteelement = P2Enum;       break;
    437                                 case OneLayerP4zEnum     : finiteelement = P2xP4Enum;    break;
    438                                 case CrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
     430                                case P1P1Enum              : finiteelement = P1Enum;       break;
     431                                case P1P1GLSEnum           : finiteelement = P1Enum;       break;
     432                                case MINIcondensedEnum     : finiteelement = P1bubbleEnum; break;
     433                                case MINIEnum              : finiteelement = P1bubbleEnum; break;
     434                                case TaylorHoodEnum        : finiteelement = P2Enum;       break;
     435                                case XTaylorHoodEnum       : finiteelement = P2Enum;       break;
     436                                case LATaylorHoodEnum      : finiteelement = P2Enum;       break;
     437                                case LACrouzeixRaviartEnum : finiteelement = P2bubbleEnum; break;
     438                                case OneLayerP4zEnum       : finiteelement = P2xP4Enum;    break;
     439                                case CrouzeixRaviartEnum   : finiteelement = P2bubbleEnum; break;
    439440                                default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
    440441                        }
    … …  
    882883                if (fe_FS==XTaylorHoodEnum)
    883884                 solutionsequence_la_theta(femmodel);
    884                 else if (fe_FS==LATaylorHoodEnum)
     885                else if (fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
    885886                 solutionsequence_la(femmodel);
    886887                else if(newton>0)
    … …  
    29432944        else if(fe_FS==LATaylorHoodEnum)
    29442945         Ke1=CreateKMatrixFSViscousLATH(element);
     2946        else if(fe_FS==LACrouzeixRaviartEnum)
     2947         Ke1=CreateKMatrixFSViscousLACR(element);
    29452948        else
    29462949         Ke1=CreateKMatrixFSViscous(element);
    … …  
    30863089        return Ke;
    30873090}/*}}}*/
     3091ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLACR(Element* element){/*{{{*/
     3092
     3093        /*Intermediaries*/
     3094        int         i,dim,epssize;
     3095        IssmDouble  r,rl,Jdet,viscosity,DU,DUl;
     3096        IssmDouble      normal[3];
     3097        IssmDouble *xyz_list = NULL;
     3098        IssmDouble *xyz_list_base = NULL;
     3099
     3100        /*Get problem dimension*/
     3101        element->FindParam(&dim,DomainDimensionEnum);
     3102        element->FindParam(&r,AugmentedLagrangianREnum);
     3103        if(dim==2) epssize = 3;
     3104        else       epssize = 6;
     3105
     3106        /*Fetch number of nodes and dof for this finite element*/
     3107        int vnumnodes = element->NumberofNodesVelocity();
     3108        int pnumnodes = element->GetNumberOfNodes(P1DGEnum);
     3109        int lnumnodes = element->GetNumberOfNodes(P2Enum);
     3110        int numdof    = vnumnodes*dim;
     3111        int pnumdof   = pnumnodes;
     3112        int lnumdof   = lnumnodes;
     3113
     3114        /*Prepare coordinate system list*/
     3115        int* cs_list = xNew<int>(vnumnodes);
     3116        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3117        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3118
     3119        /*Initialize Element matrix and vectors*/
     3120        ElementMatrix* Ke       = element->NewElementMatrix(FSvelocityEnum);
     3121        IssmDouble*    B        = xNew<IssmDouble>(epssize*numdof);
     3122        IssmDouble*    Bprime   = xNew<IssmDouble>(epssize*numdof);
     3123        IssmDouble*    BtBUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
     3124        IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
     3125        IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
     3126        IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
     3127
     3128        /*Retrieve all inputs and parameters*/
     3129        element->GetVerticesCoordinates(&xyz_list);
     3130        Input* vx_input = element->GetInput(VxEnum);     _assert_(vx_input);
     3131        Input* vy_input = element->GetInput(VyEnum);     _assert_(vy_input);
     3132        Input* vz_input = NULL;
     3133        if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
     3134
     3135        /* Start  looping on the number of gaussian points: */
     3136        Gauss* gauss = element->NewGauss(5);
     3137        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3138                gauss->GaussPoint(ig);
     3139
     3140                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3141                this->GetBFSvel(B,element,dim,xyz_list,gauss);
     3142                this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss);
     3143
     3144                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3145                for(i=0;i<epssize;i++)   D[i*epssize+i] = 2*viscosity*gauss->weight*Jdet;
     3146
     3147                TripleMultiply(B,epssize,numdof,1,
     3148                                        D,epssize,epssize,0,
     3149                                        Bprime,epssize,numdof,0,
     3150                                        &Ke->values[0],1);
     3151
     3152                this->GetBFSUzawa(BU,element,dim,xyz_list,gauss);
     3153                this->GetBFSprimeUzawa(BprimeU,element,dim,xyz_list,gauss);
     3154
     3155                DU = gauss->weight*Jdet*sqrt(r);
     3156
     3157                TripleMultiply(BU,1,pnumdof,1,
     3158                                        &DU,1,1,0,
     3159                                        BprimeU,1,numdof,0,
     3160                                        BtBUzawa,1);
     3161        }
     3162
     3163        if(element->IsOnBase() && 0){
     3164                element->FindParam(&rl,AugmentedLagrangianRlambdaEnum);
     3165                element->GetVerticesCoordinatesBase(&xyz_list_base);
     3166                element->NormalBase(&normal[0],xyz_list_base);
     3167
     3168                IssmDouble* Dlambda  = xNewZeroInit<IssmDouble>(dim*dim);
     3169                IssmDouble* C        = xNewZeroInit<IssmDouble>(dim*lnumdof);
     3170                IssmDouble* Cprime   = xNewZeroInit<IssmDouble>(dim*numdof);
     3171                IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
     3172
     3173                delete gauss;
     3174                gauss = element->NewGaussBase(5);
     3175                for(int ig=gauss->begin();ig<gauss->end();ig++){
     3176                        gauss->GaussPoint(ig);
     3177
     3178                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3179                        this->GetCFS(C,element,dim,xyz_list,gauss);
     3180                        this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
     3181                        for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
     3182                        TripleMultiply(C,dim,lnumdof,1,
     3183                                                Dlambda,dim,dim,0,
     3184                                                Cprime,dim,numdof,0,
     3185                                                CtCUzawa,1);
     3186                }
     3187
     3188                /*The sigma naugmentation should not be transformed*/
     3189                MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
     3190                                        CtCUzawa,lnumdof,numdof,0,
     3191                                        &Ke->values[0],1);
     3192
     3193                /*Delete base part*/
     3194                xDelete<IssmDouble>(Dlambda);
     3195                xDelete<IssmDouble>(C);
     3196                xDelete<IssmDouble>(Cprime);
     3197                xDelete<IssmDouble>(CtCUzawa);
     3198                xDelete<IssmDouble>(xyz_list_base);
     3199        }
     3200
     3201        /*Transform Coordinate System*/
     3202        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     3203
     3204        /*The pressure augmentation should not be transformed*/
     3205        MatrixMultiply(BtBUzawa,pnumdof,numdof,1,
     3206                                BtBUzawa,pnumdof,numdof,0,
     3207                                &Ke->values[0],1);
     3208
     3209        /*Clean up and return*/
     3210        delete gauss;
     3211        xDelete<IssmDouble>(xyz_list);
     3212        xDelete<IssmDouble>(D);
     3213        xDelete<IssmDouble>(Bprime);
     3214        xDelete<IssmDouble>(B);
     3215        xDelete<IssmDouble>(BprimeU);
     3216        xDelete<IssmDouble>(BU);
     3217        xDelete<IssmDouble>(BtBUzawa);
     3218        xDelete<int>(cs_list);
     3219        return Ke;
     3220}/*}}}*/
    30883221ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousXTH(Element* element){/*{{{*/
    30893222
    … …  
    34463579        else if(fe_FS==LATaylorHoodEnum){
    34473580                ElementVector* pe2=CreatePVectorFSViscousLATH(element);
     3581                ElementVector* pe3 = new ElementVector(pe,pe2);
     3582                delete pe;
     3583                delete pe2;
     3584                return pe3;
     3585        }
     3586        else if(fe_FS==LACrouzeixRaviartEnum){
     3587                ElementVector* pe2=CreatePVectorFSViscousLACR(element);
    34483588                ElementVector* pe3 = new ElementVector(pe,pe2);
    34493589                delete pe;
    … …  
    36743814                ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
    36753815                ElementVector* pe4=CreatePVectorFSViscousLATH(element);
     3816                pe = new ElementVector(petemp,pe4);
     3817                delete pe1;
     3818                delete pe2;
     3819                delete pe3;
     3820                delete petemp;
     3821                delete pe4;
     3822        }
     3823        else if(fe_FS==LACrouzeixRaviartEnum){
     3824                ElementVector* pe1=CreatePVectorFSViscous(element);
     3825                ElementVector* pe2=CreatePVectorFSShelf(element);
     3826                ElementVector* pe3=CreatePVectorFSFront(element);
     3827                ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
     3828                ElementVector* pe4=CreatePVectorFSViscousLACR(element);
    36763829                pe = new ElementVector(petemp,pe4);
    36773830                delete pe1;
    … …  
    39374090}/*}}}*/
    39384091ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLATH(Element* element){/*{{{*/
     4092
     4093        int         i,dim;
     4094        IssmDouble  Jdet,r,pressure;
     4095        IssmDouble  bed_normal[3];
     4096        IssmDouble *xyz_list      = NULL;
     4097        IssmDouble *xyz_list_base = NULL;
     4098        Gauss*      gauss         = NULL;
     4099
     4100        /*Get problem dimension*/
     4101        element->FindParam(&dim,DomainDimensionEnum);
     4102
     4103        /*Fetch number of nodes and dof for this finite element*/
     4104        int numnodes = element->GetNumberOfNodes();
     4105
     4106        /*Prepare coordinate system list*/
     4107        int* cs_list = xNew<int>(numnodes);
     4108        if(dim==2) for(i=0;i<numnodes;i++) cs_list[i] = XYEnum;
     4109        else       for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
     4110
     4111        /*Initialize vectors*/
     4112        ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
     4113        IssmDouble*    dbasis  = xNew<IssmDouble>(3*numnodes);
     4114
     4115        /*Retrieve all inputs and parameters*/
     4116        element->FindParam(&r,AugmentedLagrangianREnum);
     4117        element->GetVerticesCoordinates(&xyz_list);
     4118
     4119        /*Get pressure and sigmann*/
     4120        Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
     4121        Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
     4122
     4123        gauss=element->NewGauss(5);
     4124        for(int ig=gauss->begin();ig<gauss->end();ig++){
     4125                gauss->GaussPoint(ig);
     4126                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     4127               
     4128                pressure_input->GetInputValue(&pressure, gauss);
     4129                element->NodalFunctionsDerivativesVelocity(dbasis,xyz_list,gauss);
     4130
     4131                for(i=0;i<numnodes;i++){
     4132                        pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
     4133                        pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
     4134                        if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
     4135                }
     4136        }
     4137
     4138        if(element->IsOnBase()){
     4139                IssmDouble   sigmann;
     4140                IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
     4141
     4142                element->GetVerticesCoordinatesBase(&xyz_list_base);
     4143                element->NormalBase(&bed_normal[0],xyz_list_base);
     4144
     4145                delete gauss;
     4146                gauss=element->NewGaussBase(5);
     4147                for(int ig=gauss->begin();ig<gauss->end();ig++){
     4148                        gauss->GaussPoint(ig);
     4149
     4150                        element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     4151                        element->NodalFunctionsVelocity(vbasis,gauss);
     4152                        sigmann_input->GetInputValue(&sigmann, gauss);
     4153
     4154                        for(i=0;i<numnodes;i++){
     4155                                pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
     4156                                pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
     4157                                if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
     4158                        }
     4159                }
     4160                xDelete<IssmDouble>(xyz_list_base);
     4161                xDelete<IssmDouble>(vbasis);
     4162        }
     4163
     4164        /*Transform coordinate system*/
     4165        //element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
     4166
     4167        /*Clean up and return*/
     4168        delete gauss;
     4169        xDelete<int>(cs_list);
     4170        xDelete<IssmDouble>(xyz_list);
     4171        xDelete<IssmDouble>(dbasis);
     4172        return pe;
     4173}/*}}}*/
     4174ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLACR(Element* element){/*{{{*/
    39394175
    39404176        int         i,dim;

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h

    r18293 r18474  
    7070                ElementMatrix* CreateKMatrixFS(Element* element);
    7171                ElementMatrix* CreateKMatrixFSViscousLATH(Element* element);
     72                ElementMatrix* CreateKMatrixFSViscousLACR(Element* element);
    7273                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
    7374                ElementMatrix* CreateKMatrixFSViscous(Element* element);
    … …  
    7778                ElementVector* CreatePVectorFSViscous(Element* element);
    7879                ElementVector* CreatePVectorFSViscousLATH(Element* element);
     80                ElementVector* CreatePVectorFSViscousLACR(Element* element);
    7981                ElementVector* CreatePVectorFSViscousXTH(Element* element);
    8082                ElementVector* CreatePVectorFSShelf(Element* element);

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

    r18290 r18474  
    1717
    1818        /*Update elements: */
    19         int finiteelement = P1Enum;
     19        int finiteelement;
    2020        int counter=0;
     21        int fe_FS;
     22
     23        iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
     24        if(fe_FS=LATaylorHoodEnum) finiteelement = P1Enum;
     25        else if(fe_FS=LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
     26        else _error_("solution not supported yet");
     27
    2128        for(int i=0;i<iomodel->numberofelements;i++){
    2229                if(iomodel->my_elements[i]){
    … …  
    3542void UzawaPressureAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
    3643
    37         int finiteelement = P1Enum;
     44        int finiteelement;
     45        int fe_FS;
     46
     47        iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
     48        if(fe_FS=LATaylorHoodEnum) finiteelement = P1Enum;
     49        else if(fe_FS=LACrouzeixRaviartEnum) finiteelement = P1DGEnum;
     50        else _error_("solution not supported yet");
     51
    3852        ::CreateNodes(nodes,iomodel,UzawaPressureAnalysisEnum,finiteelement);
    3953}/*}}}*/
    … …  
    176190        /*Fetch number of nodes and dof for this finite element*/
    177191        int numnodes      = element->GetNumberOfNodes();
    178         int numnodessigma = element->GetNumberOfNodes(P2Enum);
     192        int numnodessigma;
     193        if(element->element_type==P1Enum) numnodessigma=element->GetNumberOfNodes(P2Enum);
     194        else if(element->element_type==P1DGEnum) numnodessigma=element->GetNumberOfNodes(P2Enum);
     195        else _error_("finite element not supported yet");
    179196
    180197        element->FindParam(&dim,DomainDimensionEnum);

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

    r18450 r18474  
    24172417                        penta_node_ids[24]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+iomodel->numberofelements+6*index+6;
    24182418                        break;
     2419                case LACrouzeixRaviartEnum:
     2420                        numnodes         = 19;
     2421                        penta_node_ids   = xNew<int>(numnodes);
     2422                        penta_node_ids[ 0]=iomodel->nodecounter+iomodel->elements[6*index+0];
     2423                        penta_node_ids[ 1]=iomodel->nodecounter+iomodel->elements[6*index+1];
     2424                        penta_node_ids[ 2]=iomodel->nodecounter+iomodel->elements[6*index+2];
     2425                        penta_node_ids[ 3]=iomodel->nodecounter+iomodel->elements[6*index+3];
     2426                        penta_node_ids[ 4]=iomodel->nodecounter+iomodel->elements[6*index+4];
     2427                        penta_node_ids[ 5]=iomodel->nodecounter+iomodel->elements[6*index+5];
     2428                        penta_node_ids[ 6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+0]+1;
     2429                        penta_node_ids[ 7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+1]+1;
     2430                        penta_node_ids[ 8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+2]+1;
     2431                        penta_node_ids[ 9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+3]+1;
     2432                        penta_node_ids[10]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+4]+1;
     2433                        penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+5]+1;
     2434                        penta_node_ids[12]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+6]+1;
     2435                        penta_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+7]+1;
     2436                        penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1;
     2437                        penta_node_ids[15]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+2]+1;
     2438                        penta_node_ids[16]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+3]+1;
     2439                        penta_node_ids[17]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elementtofaceconnectivity[5*index+4]+1;
     2440                        penta_node_ids[18]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberoffaces+index+1;
     2441                        break;
    24192442                default:
    24202443                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");

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

    r18212 r18474  
    966966                case OneLayerP4zEnum:       return NUMNODESP2xP4+NUMNODESP1;
    967967                case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
     968                case LACrouzeixRaviartEnum: return NUMNODESP2b;
    968969                default:       _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    969970        }
    … …  
    975976
    976977        switch(fe_stokes){
    977                 case P1P1Enum:           return P1Enum;
    978                 case P1P1GLSEnum:        return P1Enum;
    979                 case MINIcondensedEnum:  return P1bubbleEnum;
    980                 case MINIEnum:           return P1bubbleEnum;
    981                 case TaylorHoodEnum:     return P2Enum;
    982                 case LATaylorHoodEnum:   return P2Enum;
    983                 case OneLayerP4zEnum:    return P2xP4Enum;
    984                 case CrouzeixRaviartEnum:return P2bubbleEnum;
     978                case P1P1Enum:              return P1Enum;
     979                case P1P1GLSEnum:           return P1Enum;
     980                case MINIcondensedEnum:     return P1bubbleEnum;
     981                case MINIEnum:              return P1bubbleEnum;
     982                case TaylorHoodEnum:        return P2Enum;
     983                case LATaylorHoodEnum:      return P2Enum;
     984                case OneLayerP4zEnum:       return P2xP4Enum;
     985                case CrouzeixRaviartEnum:   return P2bubbleEnum;
     986                case LACrouzeixRaviartEnum: return P2bubbleEnum;
    985987                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    986988        }
    … …  
    992994
    993995        switch(fe_stokes){
    994                 case P1P1Enum:           return P1Enum;
    995                 case P1P1GLSEnum:        return P1Enum;
    996                 case MINIcondensedEnum:  return P1Enum;
    997                 case MINIEnum:           return P1Enum;
    998                 case TaylorHoodEnum:     return P1Enum;
    999                 case LATaylorHoodEnum:   return NoneEnum;
    1000                 case OneLayerP4zEnum:    return P1Enum;
    1001                 case CrouzeixRaviartEnum:return P1DGEnum;
     996                case P1P1Enum:              return P1Enum;
     997                case P1P1GLSEnum:           return P1Enum;
     998                case MINIcondensedEnum:     return P1Enum;
     999                case MINIEnum:              return P1Enum;
     1000                case TaylorHoodEnum:        return P1Enum;
     1001                case LATaylorHoodEnum:      return NoneEnum;
     1002                case OneLayerP4zEnum:       return P1Enum;
     1003                case CrouzeixRaviartEnum:   return P1DGEnum;
     1004                case LACrouzeixRaviartEnum: return NoneEnum;
    10021005                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    10031006        }

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

    r18403 r18474  
    20512051                        tria_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+3;
    20522052                        break;
     2053                case LACrouzeixRaviartEnum:
     2054                        numnodes        = 7;
     2055                        tria_node_ids   = xNew<int>(numnodes);
     2056                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     2057                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     2058                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     2059                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
     2060                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
     2061                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
     2062                        tria_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+index+1;
     2063                        break;
    20532064                default:
    20542065                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");

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

    r18182 r18474  
    408408
    409409        switch(finiteelement){
    410                 case NoneEnum:              return 0;
    411                 case P0Enum:                return NUMNODESP0;
    412                 case P1Enum:                return NUMNODESP1;
    413                 case P1DGEnum:              return NUMNODESP1;
    414                 case P1bubbleEnum:          return NUMNODESP1b;
    415                 case P1bubblecondensedEnum: return NUMNODESP1b;
    416                 case P2Enum:                return NUMNODESP2;
    417                 case P2bubbleEnum:          return NUMNODESP2b;
    418                 case P2bubblecondensedEnum: return NUMNODESP2b;
    419                 case P1P1Enum:              return NUMNODESP1*2;
    420                 case P1P1GLSEnum:           return NUMNODESP1*2;
    421                 case MINIcondensedEnum:     return NUMNODESP1b+NUMNODESP1;
    422                 case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    423                 case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
    424                 case LATaylorHoodEnum:      return NUMNODESP2;
    425                 case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
    426                 case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
     410                case NoneEnum:                return 0;
     411                case P0Enum:                  return NUMNODESP0;
     412                case P1Enum:                  return NUMNODESP1;
     413                case P1DGEnum:                return NUMNODESP1;
     414                case P1bubbleEnum:            return NUMNODESP1b;
     415                case P1bubblecondensedEnum:   return NUMNODESP1b;
     416                case P2Enum:                  return NUMNODESP2;
     417                case P2bubbleEnum:            return NUMNODESP2b;
     418                case P2bubblecondensedEnum:   return NUMNODESP2b;
     419                case P1P1Enum:                return NUMNODESP1*2;
     420                case P1P1GLSEnum:             return NUMNODESP1*2;
     421                case MINIcondensedEnum:       return NUMNODESP1b+NUMNODESP1;
     422                case MINIEnum:                return NUMNODESP1b+NUMNODESP1;
     423                case TaylorHoodEnum:          return NUMNODESP2+NUMNODESP1;
     424                case LATaylorHoodEnum:        return NUMNODESP2;
     425                case XTaylorHoodEnum:         return NUMNODESP2+NUMNODESP1;
     426                case CrouzeixRaviartEnum:     return NUMNODESP2b+NUMNODESP1;
     427                case LACrouzeixRaviartEnum:   return NUMNODESP2b;
    427428                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    428429        }
    … …  
    434435
    435436        switch(fe_stokes){
    436                 case P1P1Enum:           return P1Enum;
    437                 case P1P1GLSEnum:        return P1Enum;
    438                 case MINIcondensedEnum:  return P1bubbleEnum;
    439                 case MINIEnum:           return P1bubbleEnum;
    440                 case TaylorHoodEnum:     return P2Enum;
    441                 case LATaylorHoodEnum:   return P2Enum;
    442                 case XTaylorHoodEnum:    return P2Enum;
    443                 case CrouzeixRaviartEnum:return P2bubbleEnum;
     437                case P1P1Enum:              return P1Enum;
     438                case P1P1GLSEnum:           return P1Enum;
     439                case MINIcondensedEnum:     return P1bubbleEnum;
     440                case MINIEnum:              return P1bubbleEnum;
     441                case TaylorHoodEnum:        return P2Enum;
     442                case LATaylorHoodEnum:      return P2Enum;
     443                case XTaylorHoodEnum:       return P2Enum;
     444                case CrouzeixRaviartEnum:   return P2bubbleEnum;
     445                case LACrouzeixRaviartEnum: return P2bubbleEnum;
    444446                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    445447        }
    … …  
    451453
    452454        switch(fe_stokes){
    453                 case P1P1Enum:            return P1Enum;
    454                 case P1P1GLSEnum:         return P1Enum;
    455                 case MINIcondensedEnum:   return P1Enum;
    456                 case MINIEnum:            return P1Enum;
    457                 case TaylorHoodEnum:      return P1Enum;
    458                 case LATaylorHoodEnum:    return NoneEnum;
    459                 case XTaylorHoodEnum:     return P1Enum;
    460                 case CrouzeixRaviartEnum: return P1DGEnum;
     455                case P1P1Enum:              return P1Enum;
     456                case P1P1GLSEnum:           return P1Enum;
     457                case MINIcondensedEnum:     return P1Enum;
     458                case MINIEnum:              return P1Enum;
     459                case TaylorHoodEnum:        return P1Enum;
     460                case LATaylorHoodEnum:      return NoneEnum;
     461                case XTaylorHoodEnum:       return P1Enum;
     462                case CrouzeixRaviartEnum:   return P1DGEnum;
     463                case LACrouzeixRaviartEnum: return NoneEnum;
    461464                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    462465        }

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

    r18179 r18474  
    533533                        }
    534534                        break;
     535                case LACrouzeixRaviartEnum:
     536                        _assert_(approximation==FSApproximationEnum);
     537                        /*P2b velocity*/
     538                        EdgesPartitioning(&my_edges,iomodel);
     539                        for(i=0;i<iomodel->numberofvertices;i++){
     540                                if(iomodel->my_vertices[i]){
     541                                        nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
     542                                }
     543                        }
     544                        for(i=0;i<iomodel->numberofedges;i++){
     545                                if(my_edges[i]){
     546                                        nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
     547                                }
     548                        }
     549                        id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
     550                        if(iomodel->meshelementtype==PentaEnum){
     551                                FacesPartitioning(&my_faces,iomodel);
     552                                for(i=0;i<iomodel->numberoffaces;i++){
     553                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
     554                                                if(my_faces[i]){
     555                                                        node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
     556                                                        nodes->AddObject(node);
     557                                                }
     558                                        }
     559                                        else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
     560                                                /*Nothing*/
     561                                        }
     562                                        else{
     563                                                _error_("not supported");
     564                                        }
     565                                }
     566                                id0 = id0+iomodel->numberoffaces;
     567                        }
     568                        for(i=0;i<iomodel->numberofelements;i++){
     569                                if(iomodel->my_elements[i]){
     570                                        nodes->AddObject(new Node(id0+i+1,id0-iomodel->nodecounter+i,lid++,0,iomodel,analysis,FSvelocityEnum));
     571                                }
     572                        }
     573
     574                        /*No pressure*/
     575                        break;
    535576
    536577                default:

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r18403 r18474  
    609609        OneLayerP4zEnum,
    610610        CrouzeixRaviartEnum,
     611        LACrouzeixRaviartEnum,
    611612        /*}}}*/
    612613        /*Results{{{*/

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r18403 r18474  
    598598                case OneLayerP4zEnum : return "OneLayerP4z";
    599599                case CrouzeixRaviartEnum : return "CrouzeixRaviart";
     600                case LACrouzeixRaviartEnum : return "LACrouzeixRaviart";
    600601                case SaveResultsEnum : return "SaveResults";
    601602                case BoolExternalResultEnum : return "BoolExternalResult";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r18403 r18474  
    610610              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
    611611              else if (strcmp(name,"CrouzeixRaviart")==0) return CrouzeixRaviartEnum;
     612              else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum;
    612613              else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum;
    613614              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    … …  
    628629              else if (strcmp(name,"MisfitModelEnum")==0) return MisfitModelEnumEnum;
    629630              else if (strcmp(name,"MisfitObservation")==0) return MisfitObservationEnum;
    630               else if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
     634              if (strcmp(name,"MisfitObservationEnum")==0) return MisfitObservationEnumEnum;
     635              else if (strcmp(name,"MisfitTimeinterpolation")==0) return MisfitTimeinterpolationEnum;
    635636              else if (strcmp(name,"MisfitWeights")==0) return MisfitWeightsEnum;
    636637              else if (strcmp(name,"MisfitWeightsEnum")==0) return MisfitWeightsEnumEnum;
Note: See TracChangeset for help on using the changeset viewer.