Changeset 17525


Ignore:
Timestamp:
03/24/14 11:39:40 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: working on XTaylorhood

Location:
issm/trunk-jpl/src/c
Files:
1 added
7 edited

Legend:

Unmodified
Added
Removed

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

    r17517 r17525  
    355355                                        ./solutionsequences/solutionsequence_nonlinear.cpp\
    356356                                        ./solutionsequences/solutionsequence_newton.cpp\
     357                                        ./solutionsequences/solutionsequence_la_theta.cpp\
    357358                                        ./solutionsequences/convergence.cpp\
    358359                                        ./classes/Options/Options.h\

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

    r17516 r17525  
    407407                                case MINIEnum          : finiteelement = P1bubbleEnum; break;
    408408                                case TaylorHoodEnum    : finiteelement = P2Enum;       break;
     409                                case XTaylorHoodEnum   : finiteelement = P2Enum;       break;
    409410                                case OneLayerP4zEnum   : finiteelement = P2xP4Enum;    break;
    410                                 default: _error_("finite element "<<finiteelement<<" not supported");
     411                                default: _error_("finite element "<<EnumToStringx(finiteelement)<<" not supported");
    411412                        }
    412413                }
    … …  
    826827        bool isSSA,isL1L2,isHO,isFS;
    827828        bool conserve_loads = true;
    828         int  newton,meshtype;
     829        int  newton,meshtype,fe_FS;
    829830
    830831        /* recover parameters:*/
    … …  
    833834        femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
    834835        femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
     836        femmodel->parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
    835837        femmodel->parameters->FindParam(&meshtype,MeshTypeEnum);
    836838        femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum);
    837839
    838         if((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor
     840        if(isFS){
     841                if(VerboseSolution()) _printf0_("   computing velocities\n");
     842
     843                femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
     844                if (fe_FS==XTaylorHoodEnum)
     845                 solutionsequence_la_theta(femmodel);
     846                else if(newton>0)
     847                 solutionsequence_newton(femmodel);
     848                else
     849                 solutionsequence_nonlinear(femmodel,conserve_loads);
     850        }
     851        else if(isSSA || isHO || isL1L2){
    839852                if(VerboseSolution()) _printf0_("   computing velocities\n");
    840853
    … …  
    851864                        extrudefrombase_core(femmodel);
    852865                }
     866        }
     867        else{
     868                _error_("not supported");
    853869        }
    854870
    … …  
    28572873ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/
    28582874
     2875        /*Get type of algorithm*/
     2876        int fe_FS;
     2877        element->FindParam(&fe_FS,FlowequationFeFSEnum);
     2878
    28592879        /*compute all stiffness matrices for this element*/
    2860         ElementMatrix* Ke1=CreateKMatrixFSViscous(element);
     2880        ElementMatrix* Ke1=NULL;
     2881        if(fe_FS==XTaylorHoodEnum)
     2882         Ke1=CreateKMatrixFSViscousXTH(element);
     2883        else
     2884         Ke1=CreateKMatrixFSViscous(element);
     2885
    28612886        ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
    28622887        ElementMatrix* Ke3=CreateKMatrixFSShelf(element);
    … …  
    28672892        delete Ke2;
    28682893        delete Ke3;
     2894        return Ke;
     2895}/*}}}*/
     2896ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscousXTH(Element* element){/*{{{*/
     2897
     2898        /*Intermediaries*/
     2899        int         i,meshtype,dim,epssize;
     2900        IssmDouble  r,FSreconditioning,Jdet;
     2901        IssmDouble *xyz_list = NULL;
     2902
     2903        /*FIXME this should change*/
     2904        r = 1.;
     2905
     2906        /*Get problem dimension*/
     2907        element->FindParam(&meshtype,MeshTypeEnum);
     2908        switch(meshtype){
     2909                case Mesh2DverticalEnum: dim = 2; epssize = 3; break;
     2910                case Mesh3DEnum:         dim = 3; epssize = 6; break;
     2911                case Mesh3DtetrasEnum:   dim = 3; epssize = 6; break;
     2912                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     2913        }
     2914
     2915        /*Fetch number of nodes and dof for this finite element*/
     2916        int vnumnodes = element->NumberofNodesVelocity();
     2917        int pnumnodes = element->NumberofNodesPressure();
     2918        int numdof    = vnumnodes*dim + pnumnodes;
     2919        int bsize     = epssize + 2;
     2920
     2921        /*Prepare coordinate system list*/
     2922        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     2923        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     2924        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     2925        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     2926
     2927        /*Initialize Element matrix and vectors*/
     2928        ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
     2929        IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
     2930        IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
     2931        IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
     2932
     2933        /*Retrieve all inputs and parameters*/
     2934        element->GetVerticesCoordinates(&xyz_list);
     2935        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     2936        Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
     2937        Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
     2938        Input* vz_input;
     2939        if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     2940
     2941        /* Start  looping on the number of gaussian points: */
     2942        Gauss* gauss = element->NewGauss(5);
     2943        for(int ig=gauss->begin();ig<gauss->end();ig++){
     2944                gauss->GaussPoint(ig);
     2945
     2946                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     2947                this->GetBFS(B,element,dim,xyz_list,gauss);
     2948                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
     2949
     2950                for(i=0;i<epssize;i++)     D[i*bsize+i] = + r*gauss->weight*Jdet;
     2951                for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
     2952
     2953                TripleMultiply(B,bsize,numdof,1,
     2954                                        D,bsize,bsize,0,
     2955                                        Bprime,bsize,numdof,0,
     2956                                        &Ke->values[0],1);
     2957        }
     2958
     2959        /*Transform Coordinate System*/
     2960        element->TransformStiffnessMatrixCoord(Ke,cs_list);
     2961
     2962        /*Clean up and return*/
     2963        delete gauss;
     2964        xDelete<IssmDouble>(xyz_list);
     2965        xDelete<IssmDouble>(D);
     2966        xDelete<IssmDouble>(Bprime);
     2967        xDelete<IssmDouble>(B);
     2968        xDelete<int>(cs_list);
    28692969        return Ke;
    28702970}/*}}}*/
    … …  
    31793279ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
    31803280
    3181         /*compute all load vectors for this element*/
    3182         ElementVector* pe1=CreatePVectorFSViscous(element);
    3183         ElementVector* pe2=CreatePVectorFSShelf(element);
    3184         ElementVector* pe3=CreatePVectorFSFront(element);
    3185         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3281        ElementVector* pe = NULL;
     3282
     3283        int fe_FS;
     3284        element->FindParam(&fe_FS,FlowequationFeFSEnum);
     3285
     3286        if(fe_FS==XTaylorHoodEnum){
     3287                ElementVector* pe1=CreatePVectorFSViscous(element);
     3288                ElementVector* pe2=CreatePVectorFSShelf(element);
     3289                ElementVector* pe3=CreatePVectorFSFront(element);
     3290                ElementVector* petemp =new ElementVector(pe1,pe2,pe3);
     3291                ElementVector* pe4=CreatePVectorFSViscousXTH(element);
     3292                ElementVector* pe = new ElementVector(petemp,pe4);
     3293                delete pe1;
     3294                delete pe2;
     3295                delete pe3;
     3296                delete petemp;
     3297                delete pe4;
     3298        }
     3299        else{
     3300                ElementVector* pe1=CreatePVectorFSViscous(element);
     3301                ElementVector* pe2=CreatePVectorFSShelf(element);
     3302                ElementVector* pe3=CreatePVectorFSFront(element);
     3303                pe =new ElementVector(pe1,pe2,pe3);
     3304                delete pe1;
     3305                delete pe2;
     3306                delete pe3;
     3307        }
    31863308
    31873309        /*clean-up and return*/
    3188         delete pe1;
    3189         delete pe2;
    3190         delete pe3;
    3191 
    31923310        return pe;
    31933311}/*}}}*/
    31943312#endif
    31953313ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
     3314
     3315        int         i,meshtype,dim;
     3316        IssmDouble  Jdet,forcex,forcey,forcez;
     3317        IssmDouble *xyz_list = NULL;
     3318
     3319        /*Get problem dimension*/
     3320        element->FindParam(&meshtype,MeshTypeEnum);
     3321        switch(meshtype){
     3322                case Mesh2DverticalEnum: dim = 2; break;
     3323                case Mesh3DEnum:         dim = 3; break;
     3324                case Mesh3DtetrasEnum:   dim = 3; break;
     3325                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     3326        }
     3327
     3328        /*Fetch number of nodes and dof for this finite element*/
     3329        int vnumnodes = element->NumberofNodesVelocity();
     3330        int pnumnodes = element->NumberofNodesPressure();
     3331
     3332        /*Prepare coordinate system list*/
     3333        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     3334        if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     3335        else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     3336        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     3337
     3338        /*Initialize vectors*/
     3339        ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
     3340        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3341
     3342        /*Retrieve all inputs and parameters*/
     3343        element->GetVerticesCoordinates(&xyz_list);
     3344        IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     3345        IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
     3346        Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
     3347        Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
     3348        Input*      loadingforcez_input=NULL;
     3349        if(dim==3){
     3350                loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
     3351        }
     3352
     3353        /* Start  looping on the number of gaussian points: */
     3354        Gauss* gauss=element->NewGauss(5);
     3355        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3356                gauss->GaussPoint(ig);
     3357
     3358                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     3359                element->NodalFunctionsVelocity(vbasis,gauss);
     3360
     3361                loadingforcex_input->GetInputValue(&forcex,gauss);
     3362                loadingforcey_input->GetInputValue(&forcey,gauss);
     3363                if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
     3364
     3365                for(i=0;i<vnumnodes;i++){
     3366                        pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
     3367                        pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
     3368                        if(dim==3){
     3369                                pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
     3370                                pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     3371                        }
     3372                        else{
     3373                                pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
     3374                        }
     3375                }
     3376        }
     3377
     3378        /*Transform coordinate system*/
     3379        element->TransformLoadVectorCoord(pe,cs_list);
     3380
     3381        /*Clean up and return*/
     3382        delete gauss;
     3383        xDelete<int>(cs_list);
     3384        xDelete<IssmDouble>(vbasis);
     3385        xDelete<IssmDouble>(xyz_list);
     3386        return pe;
     3387}/*}}}*/
     3388ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
     3389        _error_("STOP");
    31963390
    31973391        int         i,meshtype,dim;

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

    r17411 r17525  
    6767                ElementMatrix* CreateJacobianMatrixFS(Element* element);
    6868                ElementMatrix* CreateKMatrixFS(Element* element);
     69                ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
    6970                ElementMatrix* CreateKMatrixFSViscous(Element* element);
    7071                ElementMatrix* CreateKMatrixFSFriction(Element* element);
    … …  
    7273                ElementVector* CreatePVectorFS(Element* element);
    7374                ElementVector* CreatePVectorFSViscous(Element* element);
     75                ElementVector* CreatePVectorFSViscousXTH(Element* element);
    7476                ElementVector* CreatePVectorFSShelf(Element* element);
    7577                ElementVector* CreatePVectorFSFront(Element* element);

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

    r17516 r17525  
    15961596void  Tria::ResetFSBasalBoundaryCondition(void){
    15971597
    1598         int numnodes = this->GetNumberOfNodes();
     1598        int numnodes = this->NumberofNodesVelocity();
    15991599
    16001600        int          approximation;
    … …  
    19191919                        break;
    19201920                case TaylorHoodEnum:
     1921                case XTaylorHoodEnum:
    19211922                        numnodes        = 9;
    19221923                        tria_node_ids   = xNew<int>(numnodes);

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

    r17242 r17525  
    421421                        return;
    422422                case TaylorHoodEnum:
     423                case XTaylorHoodEnum:
    423424                        this->element_type = P2Enum;
    424425                        this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    … …  
    588589                case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
    589590                case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
     591                case XTaylorHoodEnum:       return NUMNODESP2+NUMNODESP1;
    590592                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    591593        }
    … …  
    603605                case MINIEnum:          return NUMNODESP1;
    604606                case TaylorHoodEnum:    return NUMNODESP1;
     607                case XTaylorHoodEnum:   return NUMNODESP1;
    605608                default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    606609        }
    … …  
    618621                case MINIEnum:          return NUMNODESP1b;
    619622                case TaylorHoodEnum:    return NUMNODESP2;
     623                case XTaylorHoodEnum:   return NUMNODESP2;
    620624                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    621625        }
    … …  
    633637                case MINIEnum:          return P1bubbleEnum;
    634638                case TaylorHoodEnum:    return P2Enum;
     639                case XTaylorHoodEnum:   return P2Enum;
    635640                default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
    636641        }

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

    r17419 r17525  
    294294                        break;
    295295                case TaylorHoodEnum:
     296                case XTaylorHoodEnum:
    296297                        _assert_(approximation==FSApproximationEnum);
    297298                        /*P2 velocity*/

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

    r17516 r17525  
    1818void solutionsequence_FScoupling_nonlinear(FemModel* femmodel,bool conserve_loads);
    1919void solutionsequence_linear(FemModel* femmodel);
     20void solutionsequence_la_theta(FemModel* femmodel);
    2021void solutionsequence_adjoint_linear(FemModel* femmodel);
    2122
Note: See TracChangeset for help on using the changeset viewer.