Changeset 18179


Ignore:
Timestamp:
06/24/14 16:47:09 (11 years ago)
Author:
seroussi
Message:

NEW: starting to implement pressure augmentation

Location:
issm/trunk-jpl/src
Files:
5 added
19 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r18135 r18179  
    354354                                        ./analyses/EnumToAnalysis.cpp\
    355355                                        ./analyses/Analysis.h\
     356                                        ./solutionsequences/solutionsequence_la.cpp\
     357                                        ./solutionsequences/solutionsequence_la_theta.cpp\
    356358                                        ./solutionsequences/solutionsequence_linear.cpp\
    357359                                        ./solutionsequences/solutionsequence_nonlinear.cpp\
    358360                                        ./solutionsequences/solutionsequence_newton.cpp\
    359                                         ./solutionsequences/solutionsequence_la_theta.cpp\
    360361                                        ./solutionsequences/convergence.cpp\
    361362                                        ./classes/Options/Options.h\
    … …  
    538539issm_sources += ./analyses/StressbalanceAnalysis.cpp
    539540endif
     541if UZAWAPRESSURE
     542issm_sources += ./analyses/UzawaPressureAnalysis.cpp
     543endif
    540544if STRESSBALANCESIA
    541545issm_sources += ./analyses/StressbalanceSIAAnalysis.cpp

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

    r17933 r18179  
    9595                case ThermalAnalysisEnum : return new ThermalAnalysis();
    9696                #endif
     97                #ifdef _HAVE_UZAWAPRESSURE_
     98                case UzawaPressureAnalysisEnum : return new UzawaPressureAnalysis();
     99                #endif
    97100                #ifdef _HAVE_GIA_
    98101                case GiaAnalysisEnum : return new GiaAnalysis();

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

    r18173 r18179  
    867867                if (fe_FS==XTaylorHoodEnum)
    868868                 solutionsequence_la_theta(femmodel);
     869                else if (fe_FS==LATaylorHoodEnum)
     870                 solutionsequence_la(femmodel);
    869871                else if(newton>0)
    870872                 solutionsequence_newton(femmodel);
    … …  
    28522854        if(fe_FS==XTaylorHoodEnum)
    28532855         Ke1=CreateKMatrixFSViscousXTH(element);
     2856        else if(fe_FS==LATaylorHoodEnum)
     2857         Ke1=CreateKMatrixFSViscousLATH(element);
    28542858        else
    28552859         Ke1=CreateKMatrixFSViscous(element);
    … …  
    28632867        delete Ke2;
    28642868        delete Ke3;
     2869        return Ke;
     2870}/*}}}*/
     2871ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousLATH(Element* element){/*{{{*/
     2872
     2873        /*Intermediaries*/
     2874        int         i,dim,epssize;
     2875        IssmDouble  r,Jdet,viscosity,DU;
     2876        IssmDouble *xyz_list = NULL;
     2877
     2878        /*Get problem dimension*/
     2879        element->FindParam(&dim,DomainDimensionEnum);
     2880        element->FindParam(&r,AugmentedLagrangianREnum);
     2881        if(dim==2) epssize = 3;
     2882        else       epssize = 6;
     2883
     2884        /*Fetch number of nodes and dof for this finite element*/
     2885        int vnumnodes = element->NumberofNodesVelocity();
     2886        int pnumnodes;
     2887        if(dim==2) pnumnodes=3;
     2888        else pnumnodes=6;
     2889        //int pnumnodes = element->NumberofNodes(P1Enum);
     2890        int numdof    = vnumnodes*dim;
     2891        int pnumdof   = pnumnodes;
     2892
     2893        /*Prepare coordinate system list*/
     2894        int* cs_list = xNew<int>(vnumnodes);
     2895        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     2896        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     2897
     2898        /*Initialize Element matrix and vectors*/
     2899        ElementMatrix* Ke       = element->NewElementMatrix(FSvelocityEnum);
     2900        IssmDouble*    B        = xNew<IssmDouble>(epssize*numdof);
     2901        IssmDouble*    Bprime   = xNew<IssmDouble>(epssize*numdof);
     2902        IssmDouble*    BtBUzawa = xNew<IssmDouble>(numdof*pnumdof);
     2903        IssmDouble*    BU       = xNew<IssmDouble>(pnumdof);
     2904        IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
     2905        IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
     2906
     2907        /*Retrieve all inputs and parameters*/
     2908        element->GetVerticesCoordinates(&xyz_list);
     2909        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     2910        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     2911        Input* vz_input;
     2912        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     2913
     2914        /* Start  looping on the number of gaussian points: */
     2915        Gauss* gauss = element->NewGauss(5);
     2916        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2917                gauss->GaussPoint(ig);
     2918
     2919                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2920                this->GetBFS(B,element,dim,xyz_list,gauss);
     2921                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
     2922
     2923                element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     2924                for(i=0;i<epssize;i++)   D[i*epssize+i] = 2*viscosity*gauss->weight*Jdet;
     2925
     2926                TripleMultiply(B,epssize,numdof,1,
     2927                                        D,epssize,epssize,0,
     2928                                        Bprime,epssize,numdof,0,
     2929                                        &Ke->values[0],1);
     2930
     2931                this->GetBFSUzawa(BU,element,dim,xyz_list,gauss);
     2932                this->GetBFSprimeUzawa(BprimeU,element,dim,xyz_list,gauss);
     2933
     2934                DU = gauss->weight*Jdet;
     2935
     2936                TripleMultiply(BU,1,pnumdof,1,
     2937                                        &DU,1,1,0,
     2938                                        BprimeU,1,numdof,0,
     2939                                        BtBUzawa,1);
     2940
     2941        }
     2942
     2943        MatrixMultiply(BtBUzawa,pnumdof,numdof,1,
     2944                                BtBUzawa,pnumdof,numdof,0,
     2945                                &Ke->values[0],1);
     2946
     2947        /*Transform Coordinate System*/
     2948        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     2949
     2950        /*Clean up and return*/
     2951        delete gauss;
     2952        xDelete<IssmDouble>(xyz_list);
     2953        xDelete<IssmDouble>(D);
     2954        xDelete<IssmDouble>(Bprime);
     2955        xDelete<IssmDouble>(B);
     2956        xDelete<IssmDouble>(BprimeU);
     2957        xDelete<IssmDouble>(BU);
     2958        xDelete<IssmDouble>(BtBUzawa);
     2959        xDelete<int>(cs_list);
    28652960        return Ke;
    28662961}/*}}}*/
    … …  
    33053400                ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
    33063401                ElementVector* pe4=CreatePVectorFSViscousXTH(element);
     3402                pe = new ElementVector(petemp,pe4);
     3403                delete pe1;
     3404                delete pe2;
     3405                delete pe3;
     3406                delete petemp;
     3407                delete pe4;
     3408        }
     3409        else if(fe_FS==LATaylorHoodEnum){
     3410                ElementVector* pe1=CreatePVectorFSViscous(element);
     3411                ElementVector* pe2=CreatePVectorFSShelf(element);
     3412                ElementVector* pe3=CreatePVectorFSFront(element);
     3413                ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
     3414                ElementVector* pe4=CreatePVectorFSViscousLATH(element);
    33073415                pe = new ElementVector(petemp,pe4);
    33083416                delete pe1;
    … …  
    35653673        xDelete<IssmDouble>(vdbasis);
    35663674        xDelete<IssmDouble>(tbasis);
     3675        return pe;
     3676}/*}}}*/
     3677ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLATH(Element* element){/*{{{*/
     3678
     3679        int         i,dim;
     3680        IssmDouble  Jdet,r,pressure;
     3681        IssmDouble *xyz_list = NULL;
     3682        Gauss*      gauss    = NULL;
     3683
     3684        /*Get problem dimension*/
     3685        element->FindParam(&dim,DomainDimensionEnum);
     3686
     3687        /*Fetch number of nodes and dof for this finite element*/
     3688        int numnodes = element->GetNumberOfNodes();
     3689
     3690        /*Prepare coordinate system list*/
     3691        int* cs_list = xNew<int>(numnodes);
     3692        if(dim==2) for(i=0;i<numnodes;i++) cs_list[i] = XYEnum;
     3693        else       for(i=0;i<numnodes;i++) cs_list[i] = XYZEnum;
     3694
     3695        /*Initialize vectors*/
     3696        ElementVector* pe      = element->NewElementVector(FSvelocityEnum);
     3697        IssmDouble*    dbasis  = xNew<IssmDouble>(3*numnodes);
     3698
     3699        /*Retrieve all inputs and parameters*/
     3700        element->FindParam(&r,AugmentedLagrangianREnum);
     3701        element->GetVerticesCoordinates(&xyz_list);
     3702
     3703        /*Get d and tau*/
     3704        Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
     3705
     3706        gauss=element->NewGauss(5);
     3707        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3708                gauss->GaussPoint(ig);
     3709                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3710               
     3711                pressure_input->GetInputValue(&pressure, gauss);
     3712                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     3713
     3714                for(i=0;i<numnodes;i++){
     3715                        pe->values[i*dim+0] += r*pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
     3716                        pe->values[i*dim+1] += r*pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
     3717                        if(dim==3){
     3718                                pe->values[i*dim+2]+=r*pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
     3719                        }
     3720                }
     3721        }
     3722
     3723        /*Transform coordinate system*/
     3724        element->TransformLoadVectorCoord(pe,cs_list);
     3725
     3726        /*Clean up and return*/
     3727        delete gauss;
     3728        xDelete<int>(cs_list);
     3729        xDelete<IssmDouble>(xyz_list);
     3730        xDelete<IssmDouble>(dbasis);
    35673731        return pe;
    35683732}/*}}}*/
    … …  
    39554119        xDelete<IssmDouble>(pbasis);
    39564120}/*}}}*/
     4121void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4122        /*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi.
     4123         */
     4124
     4125        /*Fetch number of nodes for this finite element*/
     4126        int pnumnodes;
     4127        if(dim==2) pnumnodes=3;
     4128        else pnumnodes=6;
     4129        //int pnumnodes = element->NumberofNodes(P1Enum);
     4130
     4131        /*Get nodal functions derivatives*/
     4132        IssmDouble* basis =xNew<IssmDouble>(pnumnodes);
     4133        element->NodalFunctionsP1(basis,gauss);
     4134
     4135        /*Build B: */
     4136        for(int i=0;i<pnumnodes;i++){
     4137                B[i] = basis[i];
     4138        }
     4139
     4140        /*Clean up*/
     4141        xDelete<IssmDouble>(basis);
     4142}/*}}}*/
     4143void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
     4144        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2.
     4145         *      For node i, Bi' can be expressed in the actual coordinate system
     4146         *      by:
     4147         *                      Bvi' = [  dphi/dx   dphi/dy ]
     4148         *
     4149         *      In 3d
     4150         *         Bvi=[ dh/dx   dh/dy    dh/dz  ]
     4151         *      where phi is the finiteelement function for node i.
     4152         */
     4153
     4154        /*Fetch number of nodes for this finite element*/
     4155        int vnumnodes = element->NumberofNodesVelocity();
     4156
     4157        /*Get nodal functions derivatives*/
     4158        IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
     4159        element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
     4160
     4161        /*Build B_prime: */
     4162        if(dim==2){
     4163                for(int i=0;i<vnumnodes;i++){
     4164                        Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
     4165                        Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
     4166                }
     4167        }
     4168        else{
     4169                for(int i=0;i<vnumnodes;i++){
     4170                        Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
     4171                        Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
     4172                        Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
     4173                }
     4174        }
     4175
     4176        /*Clean up*/
     4177        xDelete<IssmDouble>(vdbasis);
     4178}/*}}}*/
    39574179void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
    39584180        /* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.

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

    r18057 r18179  
    6868                ElementMatrix* CreateJacobianMatrixFS(Element* element);
    6969                ElementMatrix* CreateKMatrixFS(Element* element);
     70                ElementMatrix* CreateKMatrixFSViscousLATH(Element* element);
    7071                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
    7172                ElementMatrix* CreateKMatrixFSViscous(Element* element);
    … …  
    7475                ElementVector* CreatePVectorFS(Element* element);
    7576                ElementVector* CreatePVectorFSViscous(Element* element);
     77                ElementVector* CreatePVectorFSViscousLATH(Element* element);
    7678                ElementVector* CreatePVectorFSViscousXTH(Element* element);
    7779                ElementVector* CreatePVectorFSShelf(Element* element);
    … …  
    8082                void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    8183                void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     84                void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
     85                void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    8286                void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
    8387                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);

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

    r17933 r18179  
    3737#include "./StressbalanceSIAAnalysis.h"
    3838#include "./StressbalanceVerticalAnalysis.h"
     39#include "./UzawaPressureAnalysis.h"
    3940#include "./L2ProjectionBaseAnalysis.h"
    4041#include "./L2ProjectionEPLAnalysis.h"

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

    r18173 r18179  
    945945
    946946        switch(finiteelement){
     947                case NoneEnum:              return 0;
    947948                case P0Enum:                return NUMNODESP0;
    948949                case P1Enum:                return NUMNODESP1;
    … …  
    962963                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    963964                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
     965                case LATaylorHoodEnum:      return NUMNODESP2;
    964966                case OneLayerP4zEnum:       return NUMNODESP2xP4+NUMNODESP1;
    965967                case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
    … …  
    978980                case MINIEnum:           return P1bubbleEnum;
    979981                case TaylorHoodEnum:     return P2Enum;
     982                case LATaylorHoodEnum:   return P2Enum;
    980983                case OneLayerP4zEnum:    return P2xP4Enum;
    981984                case CrouzeixRaviartEnum:return P2bubbleEnum;
    … …  
    994997                case MINIEnum:           return P1Enum;
    995998                case TaylorHoodEnum:     return P1Enum;
     999                case LATaylorHoodEnum:   return P1Enum;
    9961000                case OneLayerP4zEnum:    return P1Enum;
    9971001                case CrouzeixRaviartEnum:return P1DGEnum;

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

    r18163 r18179  
    863863                        tetra_node_ids[13]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[4*index+3];
    864864                        break;
     865                case LATaylorHoodEnum:
     866                        numnodes        = 10;
     867                        tetra_node_ids  = xNew<int>(numnodes);
     868                        tetra_node_ids[0]=iomodel->nodecounter+iomodel->elements[4*index+0];
     869                        tetra_node_ids[1]=iomodel->nodecounter+iomodel->elements[4*index+1];
     870                        tetra_node_ids[2]=iomodel->nodecounter+iomodel->elements[4*index+2];
     871                        tetra_node_ids[3]=iomodel->nodecounter+iomodel->elements[4*index+3];
     872                        tetra_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+0]+1;
     873                        tetra_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+1]+1;
     874                        tetra_node_ids[6]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+2]+1;
     875                        tetra_node_ids[7]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+3]+1;
     876                        tetra_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+4]+1;
     877                        tetra_node_ids[9]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[6*index+5]+1;
     878                        break;
    865879                default:
    866880                        _error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");

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

    r18143 r18179  
    354354                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    355355                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
     356                case LATaylorHoodEnum:      return NUMNODESP2;
    356357                case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
    357358                default: _error_("Element type "<<EnumToStringx(finiteelement)<<" not supported yet");
    … …  
    369370                case MINIEnum:          return P1bubbleEnum;
    370371                case TaylorHoodEnum:    return P2Enum;
     372                case LATaylorHoodEnum:  return P2Enum;
    371373                case XTaylorHoodEnum:   return P2Enum;
    372374                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");
    … …  
    384386                case MINIEnum:          return P1Enum;
    385387                case TaylorHoodEnum:    return P1Enum;
     388                case LATaylorHoodEnum:  return NoneEnum;
    386389                case XTaylorHoodEnum:   return P1Enum;
    387390                default:       _error_("Element type "<<EnumToStringx(fe_stokes)<<" not supported yet");

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

    r18173 r18179  
    19751975                        tria_node_ids[8]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2];
    19761976                        break;
     1977                case LATaylorHoodEnum:
     1978                        numnodes        = 6;
     1979                        tria_node_ids   = xNew<int>(numnodes);
     1980                        tria_node_ids[0]=iomodel->nodecounter+iomodel->elements[3*index+0];
     1981                        tria_node_ids[1]=iomodel->nodecounter+iomodel->elements[3*index+1];
     1982                        tria_node_ids[2]=iomodel->nodecounter+iomodel->elements[3*index+2];
     1983                        tria_node_ids[3]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
     1984                        tria_node_ids[4]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
     1985                        tria_node_ids[5]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
     1986                        break;
    19771987                case CrouzeixRaviartEnum:
    19781988                        numnodes        = 10;

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

    r18173 r18179  
    406406
    407407        switch(finiteelement){
     408                case NoneEnum:              return 0;
    408409                case P0Enum:                return NUMNODESP0;
    409410                case P1Enum:                return NUMNODESP1;
    … …  
    419420                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    420421                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
     422                case LATaylorHoodEnum:      return NUMNODESP2;
    421423                case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
    422424                case CrouzeixRaviartEnum:   return NUMNODESP2b+NUMNODESP1;
    … …  
    435437                case MINIEnum:           return P1bubbleEnum;
    436438                case TaylorHoodEnum:     return P2Enum;
     439                case LATaylorHoodEnum:   return P2Enum;
    437440                case XTaylorHoodEnum:    return P2Enum;
    438441                case CrouzeixRaviartEnum:return P2bubbleEnum;
    … …  
    451454                case MINIEnum:            return P1Enum;
    452455                case TaylorHoodEnum:      return P1Enum;
     456                case LATaylorHoodEnum:    return NoneEnum;
    453457                case XTaylorHoodEnum:     return P1Enum;
    454458                case CrouzeixRaviartEnum: return P1DGEnum;

  • issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp

    r17831 r18179  
    2626
    2727                case StressbalanceSolutionEnum:
    28                         numanalyses=6;
     28                        numanalyses=7;
    2929                        analyses=xNew<int>(numanalyses);
    3030                        analyses[0]=StressbalanceAnalysisEnum;
    … …  
    3434                        analyses[4]=ExtrudeFromBaseAnalysisEnum;
    3535                        analyses[5]=DepthAverageAnalysisEnum;
     36                        analyses[6]=UzawaPressureAnalysisEnum;
    3637                        break;
    3738

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

    r18173 r18179  
    378378                        }
    379379                        break;
     380                case LATaylorHoodEnum:
     381                        _assert_(approximation==FSApproximationEnum);
     382                        /*P2 velocity*/
     383                        EdgesPartitioning(&my_edges,iomodel);
     384                        for(i=0;i<iomodel->numberofvertices;i++){
     385                                if(iomodel->my_vertices[i]){
     386                                        nodes->AddObject(new Node(id0+i+1,i,lid++,i,iomodel,analysis,FSvelocityEnum));
     387                                }
     388                        }
     389                        for(i=0;i<iomodel->numberofedges;i++){
     390                                if(my_edges[i]){
     391                                        nodes->AddObject(new Node(id0+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,analysis,FSvelocityEnum));
     392                                }
     393                        }
     394                        id0 = id0+iomodel->numberofvertices+iomodel->numberofedges;
     395              if(iomodel->meshelementtype==PentaEnum){
     396                                FacesPartitioning(&my_faces,iomodel);
     397                                for(i=0;i<iomodel->numberoffaces;i++){
     398                                        if(iomodel->faces[i*iomodel->facescols+2]==2){/*Vertical quads*/
     399                                                if(my_faces[i]){
     400                                                        node = new Node(id0+i+1,iomodel->numberofvertices+iomodel->numberofedges+i,lid++,0,iomodel,analysis,FSvelocityEnum);
     401                                                        nodes->AddObject(node);
     402                                                }
     403                                        }
     404                                        else if(iomodel->faces[i*iomodel->facescols+2]==1){/*Triangular base/top*/
     405                                                /*Nothing*/
     406                                        }
     407                                        else{
     408                                                _error_("not supported");
     409                                        }
     410                                }
     411                        }
     412
     413                        /*No pressure*/
     414                        break;
    380415                case OneLayerP4zEnum:
    381416                        _assert_(approximation==FSApproximationEnum);

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

    r18173 r18179  
    390390        ThermalSolutionEnum,
    391391        TransientSolutionEnum,
     392        UzawaPressureAnalysisEnum,
    392393        GiaSolutionEnum,
    393394        GiaAnalysisEnum,
    … …  
    613614        MINIcondensedEnum,
    614615        TaylorHoodEnum,
     616        LATaylorHoodEnum,
    615617        XTaylorHoodEnum,
    616618        OneLayerP4zEnum,

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

    r18173 r18179  
    393393                case ThermalSolutionEnum : return "ThermalSolution";
    394394                case TransientSolutionEnum : return "TransientSolution";
     395                case UzawaPressureAnalysisEnum : return "UzawaPressureAnalysis";
    395396                case GiaSolutionEnum : return "GiaSolution";
    396397                case GiaAnalysisEnum : return "GiaAnalysis";
    … …  
    602603                case MINIcondensedEnum : return "MINIcondensed";
    603604                case TaylorHoodEnum : return "TaylorHood";
     605                case LATaylorHoodEnum : return "LATaylorHood";
    604606                case XTaylorHoodEnum : return "XTaylorHood";
    605607                case OneLayerP4zEnum : return "OneLayerP4z";

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

    r18173 r18179  
    402402              else if (strcmp(name,"ThermalSolution")==0) return ThermalSolutionEnum;
    403403              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
     404              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
    404405              else if (strcmp(name,"GiaSolution")==0) return GiaSolutionEnum;
    405406              else if (strcmp(name,"GiaAnalysis")==0) return GiaAnalysisEnum;
    … …  
    505506              else if (strcmp(name,"FractionIncrement")==0) return FractionIncrementEnum;
    506507              else if (strcmp(name,"Friction")==0) return FrictionEnum;
    507               else if (strcmp(name,"Internal")==0) return InternalEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
     511              if (strcmp(name,"Internal")==0) return InternalEnum;
     512              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    512513              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    513514              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    … …  
    614615              else if (strcmp(name,"MINIcondensed")==0) return MINIcondensedEnum;
    615616              else if (strcmp(name,"TaylorHood")==0) return TaylorHoodEnum;
     617              else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum;
    616618              else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    617619              else if (strcmp(name,"OneLayerP4z")==0) return OneLayerP4zEnum;
    … …  
    627629              else if (strcmp(name,"Time")==0) return TimeEnum;
    628630              else if (strcmp(name,"WaterColumnOld")==0) return WaterColumnOldEnum;
    629               else if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
    630               else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
     634              if (strcmp(name,"Outputdefinition")==0) return OutputdefinitionEnum;
     635              else if (strcmp(name,"OutputdefinitionList")==0) return OutputdefinitionListEnum;
     636              else if (strcmp(name,"Massfluxatgate")==0) return MassfluxatgateEnum;
    635637              else if (strcmp(name,"MassfluxatgateName")==0) return MassfluxatgateNameEnum;
    636638              else if (strcmp(name,"MassfluxatgateSegments")==0) return MassfluxatgateSegmentsEnum;

  • issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h

    r17525 r18179  
    1818void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
    1919void solutionsequence_linear(FemModel* femmodel);
     20void solutionsequence_la(FemModel* femmodel);
    2021void solutionsequence_la_theta(FemModel* femmodel);
    2122void solutionsequence_adjoint_linear(FemModel* femmodel);

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

    r18173 r18179  
    134134                                md = checkfield(md,'fieldname','flowequation.fe_SSA','values',{'P1','P1bubble','P1bubblecondensed','P2','P2bubble'});
    135135                                md = checkfield(md,'fieldname','flowequation.fe_HO' ,'values',{'P1','P1bubble','P1bubblecondensed','P1xP2','P2xP1','P2','P2bubble','P1xP3','P2xP4'});
    136                                 md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'});
     136                                md = checkfield(md,'fieldname','flowequation.fe_FS' ,'values',{'P1P1','P1P1GLS','MINIcondensed','MINI','TaylorHood','LaTaylorHood','XTaylorHood','OneLayerP4z','CrouzeixRaviart'});
    137137                                md = checkfield(md,'fieldname','flowequation.XTH_r','numel',[1],'>',0.);
    138138                                md = checkfield(md,'fieldname','flowequation.XTH_theta','numel',[1],'>=',0.,'<',0.5);

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

    r18173 r18179  
    4747                string="%s\n%s"%(string,fielddisplay(self,'fe_SSA',"Finite Element for SSA: 'P1', 'P1bubble' 'P1bubblecondensed' 'P2'"))
    4848                string="%s\n%s"%(string,fielddisplay(self,'fe_HO' ,"Finite Element for HO:  'P1' 'P1bubble' 'P1bubblecondensed' 'P1xP2' 'P2xP1' 'P2'"))
    49                 string="%s\n%s"%(string,fielddisplay(self,'fe_FS' ,"Finite Element for FS:  'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'XTaylorHood'"))
     49                string="%s\n%s"%(string,fielddisplay(self,'fe_FS' ,"Finite Element for FS:  'P1P1' (debugging only) 'P1P1GLS' 'MINIcondensed' 'MINI' 'TaylorHood' 'LATaylorHood' 'XTaylorHood'"))
    5050                string="%s\n%s"%(string,fielddisplay(self,'vertex_equation',"flow equation for each vertex"))
    5151                string="%s\n%s"%(string,fielddisplay(self,'element_equation',"flow equation for each element"))

  • issm/trunk-jpl/src/m/enum/EnumDefinitions.py

    r18170 r18179  
    385385def ThermalSolutionEnum(): return StringToEnum("ThermalSolution")[0]
    386386def TransientSolutionEnum(): return StringToEnum("TransientSolution")[0]
     387def UzawaPressureAnalysisEnum(): return StringToEnum("UzawaPressureAnalysis")[0]
    387388def GiaSolutionEnum(): return StringToEnum("GiaSolution")[0]
    388389def GiaAnalysisEnum(): return StringToEnum("GiaAnalysis")[0]
    … …  
    594595def MINIcondensedEnum(): return StringToEnum("MINIcondensed")[0]
    595596def TaylorHoodEnum(): return StringToEnum("TaylorHood")[0]
     597def LATaylorHoodEnum(): return StringToEnum("LATaylorHood")[0]
    596598def XTaylorHoodEnum(): return StringToEnum("XTaylorHood")[0]
    597599def OneLayerP4zEnum(): return StringToEnum("OneLayerP4z")[0]
     600def CrouzeixRaviartEnum(): return StringToEnum("CrouzeixRaviart")[0]
    598601def SaveResultsEnum(): return StringToEnum("SaveResults")[0]
    599602def BoolExternalResultEnum(): return StringToEnum("BoolExternalResult")[0]
Note: See TracChangeset for help on using the changeset viewer.