Changeset 18284


Ignore:
Timestamp:
07/22/14 17:04:28 (11 years ago)
Author:
seroussi
Message:

NEW: finished augmentation

Location:
issm/trunk-jpl/src/c/analyses
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r18277 r18284  
    245245        if(fe_FS==LATaylorHoodEnum){
    246246                iomodel->FetchDataToInput(elements,PressureEnum);
     247                InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
    247248        }
    248249
    … …  
    28962897        /*Fetch number of nodes and dof for this finite element*/
    28972898        int vnumnodes = element->NumberofNodesVelocity();
    2898         int pnumnodes;
    2899         if(dim==2) pnumnodes=3;
    2900         else pnumnodes=6;
    2901         int lnumnodes;
    2902         if(dim==2) lnumnodes=2;
    2903         else lnumnodes=3;
    2904         //int pnumnodes = element->NumberofNodes(P1Enum);
     2899        int pnumnodes = element->GetNumberOfNodes(P1Enum);
     2900        int lnumnodes = element->GetNumberOfNodes(P2Enum);
    29052901        int numdof    = vnumnodes*dim;
    29062902        int pnumdof   = pnumnodes;
    … …  
    29202916        IssmDouble*    BprimeU  = xNew<IssmDouble>(numdof);
    29212917        IssmDouble*    D        = xNewZeroInit<IssmDouble>(epssize*epssize);
    2922         IssmDouble*    CtCUzawa = xNewZeroInit<IssmDouble>(numdof*pnumdof);
    2923         IssmDouble*    C        = xNew<IssmDouble>(lnumdof);
    2924         IssmDouble*    Cprime   = xNew<IssmDouble>(numdof);
    29252918
    29262919        /*Retrieve all inputs and parameters*/
    29272920        element->GetVerticesCoordinates(&xyz_list);
    2928         Input* vx_input=element->GetInput(VxEnum);     _assert_(vx_input);
    2929         Input* vy_input=element->GetInput(VyEnum);     _assert_(vy_input);
    2930         Input* vz_input;
    2931         if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_input);}
     2921        Input* vx_input = element->GetInput(VxEnum);     _assert_(vx_input);
     2922        Input* vy_input = element->GetInput(VyEnum);     _assert_(vy_input);
     2923        Input* vz_input = NULL;
     2924        if(dim==3){vz_input = element->GetInput(VzEnum); _assert_(vz_input);}
    29322925
    29332926        /* Start  looping on the number of gaussian points: */
    … …  
    29572950                                        BprimeU,1,numdof,0,
    29582951                                        BtBUzawa,1);
    2959 
    29602952        }
    29612953
    … …  
    29652957                element->NormalBase(&normal[0],xyz_list_base);
    29662958
    2967                 int         lsize;
    2968                 IssmDouble* Dlambda = xNewZeroInit<IssmDouble>(dim*dim);
    2969                 IssmDouble* C       = xNewZeroInit<IssmDouble>(dim*lnumdof);
    2970                 IssmDouble* Cprime  = xNewZeroInit<IssmDouble>(dim*numdof);
     2959                IssmDouble* Dlambda  = xNewZeroInit<IssmDouble>(dim*dim);
     2960                IssmDouble* C        = xNewZeroInit<IssmDouble>(dim*lnumdof);
     2961                IssmDouble* Cprime   = xNewZeroInit<IssmDouble>(dim*numdof);
     2962                IssmDouble* CtCUzawa = xNewZeroInit<IssmDouble>(numdof*lnumdof);
    29712963
    29722964                delete gauss;
    2973                 Gauss* gauss=element->NewGaussBase(5);
     2965                gauss = element->NewGaussBase(5);
    29742966                for(int ig=gauss->begin();ig<gauss->end();ig++){
    29752967                        gauss->GaussPoint(ig);
    … …  
    29782970                        this->GetCFS(C,element,dim,xyz_list,gauss);
    29792971                        this->GetCFSprime(Cprime,element,dim,xyz_list,gauss);
    2980                         for(i=0;i<dim;i++)   Dlambda[i*epssize+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
     2972                        for(i=0;i<dim;i++) Dlambda[i*dim+i] = gauss->weight*Jdet*sqrt(normal[i]*normal[i])*sqrt(rl);
    29812973                        TripleMultiply(C,dim,lnumdof,1,
    29822974                                                Dlambda,dim,dim,0,
    … …  
    29842976                                                CtCUzawa,1);
    29852977                }
     2978
     2979                /*The sigma naugmentation should not be transformed*/
     2980                MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
     2981                                        CtCUzawa,lnumdof,numdof,0,
     2982                                        &Ke->values[0],1);
     2983
    29862984                /*Delete base part*/
    2987                 xDelete<IssmDouble>(xyz_list_base);
    29882985                xDelete<IssmDouble>(Dlambda);
    29892986                xDelete<IssmDouble>(C);
    29902987                xDelete<IssmDouble>(Cprime);
     2988                xDelete<IssmDouble>(CtCUzawa);
     2989                xDelete<IssmDouble>(xyz_list_base);
    29912990        }
    29922991
    … …  
    29972996        MatrixMultiply(BtBUzawa,pnumdof,numdof,1,
    29982997                                BtBUzawa,pnumdof,numdof,0,
    2999                                 &Ke->values[0],1);
    3000 
    3001         /*The sigma naugmentation should not be transformed*/
    3002         MatrixMultiply(CtCUzawa,lnumdof,numdof,1,
    3003                                 CtCUzawa,lnumdof,numdof,0,
    30042998                                &Ke->values[0],1);
    30052999
    … …  
    30133007        xDelete<IssmDouble>(BU);
    30143008        xDelete<IssmDouble>(BtBUzawa);
    3015         xDelete<IssmDouble>(Cprime);
    3016         xDelete<IssmDouble>(C);
    3017         xDelete<IssmDouble>(CtCUzawa);
    30183009        xDelete<int>(cs_list);
    30193010        return Ke;
    … …  
    32973288}/*}}}*/
    32983289ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
     3290
     3291        ElementVector* pe = NULL;
     3292
     3293        ElementVector* pe1=CreatePVectorFSViscous(element);
     3294        ElementVector* pe2=CreatePVectorFSFriction(element);
     3295        ElementVector* pe3=CreatePVectorFSStress(element);
     3296        pe =new ElementVector(pe1,pe2,pe3);
     3297        delete pe1;
     3298        delete pe2;
     3299        delete pe3;
     3300
     3301        /*clean-up and return*/
     3302        return pe;
     3303}/*}}}*/
     3304ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
    32993305
    33003306        int         i,dim,fe_FS;
    … …  
    33693375                return pe3;
    33703376        }
     3377        return pe;
     3378}/*}}}*/
     3379ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
     3380
     3381        if(!element->IsOnBase()) return NULL;
     3382
     3383        /*Intermediaries*/
     3384        int         dim;
     3385        IssmDouble  alpha2,Jdet;
     3386        IssmDouble  bed_normal[3];
     3387        IssmDouble *xyz_list_base = NULL;
     3388        Gauss*      gauss         = NULL;
     3389
     3390        /*Get problem dimension*/
     3391        element->FindParam(&dim,DomainDimensionEnum);
     3392
     3393        /*Fetch number of nodes and dof for this finite element*/
     3394        int vnumnodes = element->NumberofNodesVelocity();
     3395
     3396        /*Initialize Element matrix and vectors*/
     3397        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     3398        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3399
     3400        /*Retrieve all inputs and parameters*/
     3401        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3402        Input*  alpha2_input=element->GetInput(FrictionCoefficientEnum); _assert_(alpha2_input);
     3403
     3404        /* Start  looping on the number of gaussian points: */
     3405        gauss=element->NewGaussBase(3);
     3406        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3407                gauss->GaussPoint(ig);
     3408
     3409                alpha2_input->GetInputValue(&alpha2, gauss);
     3410                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3411                element->NodalFunctionsVelocity(vbasis,gauss);
     3412                element->NormalBase(&bed_normal[0],xyz_list_base);
     3413
     3414                for(int i=0;i<vnumnodes;i++){
     3415                        pe->values[i*dim+0] += - alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[1];
     3416                        pe->values[i*dim+1] += alpha2*gauss->weight*Jdet*vbasis[i]*bed_normal[0];
     3417                        if(dim==3){
     3418                                pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
     3419                        }
     3420                }
     3421
     3422        }
     3423
     3424        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3425
     3426        /*Clean up and return*/
     3427        delete gauss;
     3428        xDelete<IssmDouble>(xyz_list_base);
     3429        xDelete<IssmDouble>(vbasis);
     3430        return pe;
     3431}/*}}}*/
     3432ElementVector* StressbalanceAnalysis::CreatePVectorFSStress(Element* element){/*{{{*/
     3433
     3434        if(!element->IsOnBase()) return NULL;
     3435
     3436        /*Intermediaries*/
     3437        int         dim;
     3438        IssmDouble  sigmann,sigmant,Jdet,bedslope,beta;
     3439        IssmDouble *xyz_list_base = NULL;
     3440        Gauss*      gauss         = NULL;
     3441
     3442        /*Get problem dimension*/
     3443        element->FindParam(&dim,DomainDimensionEnum);
     3444
     3445        /*Fetch number of nodes and dof for this finite element*/
     3446        int vnumnodes = element->NumberofNodesVelocity();
     3447
     3448        /*Initialize Element matrix and vectors*/
     3449        ElementVector* pe = element->NewElementVector(FSvelocityEnum);
     3450        IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
     3451
     3452        /*Retrieve all inputs and parameters*/
     3453        element->GetVerticesCoordinatesBase(&xyz_list_base);
     3454        Input*  sigmann_input=element->GetInput(VzEnum); _assert_(sigmann_input);
     3455        Input*  sigmant_input=element->GetInput(TemperatureEnum); _assert_(sigmant_input);
     3456        Input*  bedslope_input=element->GetInput(BedSlopeXEnum);     _assert_(bedslope_input);
     3457
     3458        /* Start  looping on the number of gaussian points: */
     3459        gauss=element->NewGaussBase(3);
     3460        for(int ig=gauss->begin();ig<gauss->end();ig++){
     3461                gauss->GaussPoint(ig);
     3462
     3463                sigmann_input->GetInputValue(&sigmann, gauss);
     3464                sigmant_input->GetInputValue(&sigmant, gauss);
     3465                bedslope_input->GetInputValue(&bedslope, gauss);
     3466                element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
     3467                element->NodalFunctionsVelocity(vbasis,gauss);
     3468
     3469                beta=sqrt(1+bedslope*bedslope);
     3470                for(int i=0;i<vnumnodes;i++){
     3471                        pe->values[i*dim+0] += - (1./beta)*(-bedslope*sigmann + sigmant)*gauss->weight*Jdet*vbasis[i];
     3472                        pe->values[i*dim+1] += - (1./beta)*(sigmann + bedslope*sigmant)*gauss->weight*Jdet*vbasis[i];
     3473                        if(dim==3){
     3474                                //pe->values[i*dim+2]+= alpha2*gauss->weight*Jdet*vbasis[i];
     3475                                _error_("3d not supported yet");
     3476                        }
     3477                }
     3478
     3479        }
     3480
     3481        /*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
     3482
     3483        /*Clean up and return*/
     3484        delete gauss;
     3485        xDelete<IssmDouble>(xyz_list_base);
     3486        xDelete<IssmDouble>(vbasis);
    33713487        return pe;
    33723488}/*}}}*/
    … …  
    35023618        return pe;
    35033619}/*}}}*/
    3504 #endif
    35053620ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
    35063621
    … …  
    35723687        return pe;
    35733688}/*}}}*/
     3689#endif
    35743690ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousXTH(Element* element){/*{{{*/
    35753691
    … …  
    37723888        element->GetVerticesCoordinates(&xyz_list);
    37733889
    3774         /*Get d and tau*/
     3890        /*Get pressure and sigmann*/
    37753891        Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input);
    37763892        Input* sigmann_input =element->GetInput(SigmaNNEnum);  _assert_(sigmann_input);
    … …  
    37873903                        pe->values[i*dim+0] += pressure*gauss->weight*Jdet*dbasis[0*numnodes+i];
    37883904                        pe->values[i*dim+1] += pressure*gauss->weight*Jdet*dbasis[1*numnodes+i];
    3789                         if(dim==3){
    3790                                 pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
    3791                         }
     3905                        if(dim==3) pe->values[i*dim+2]+= pressure*gauss->weight*Jdet*dbasis[2*numnodes+i];
    37923906                }
    37933907        }
    37943908
    37953909        if(element->IsOnBase()){
    3796 
    37973910                IssmDouble   sigmann;
    37983911                IssmDouble*  vbasis = xNew<IssmDouble>(numnodes);
    … …  
    38023915
    38033916                delete gauss;
    3804                 Gauss* gauss=element->NewGaussBase(5);
     3917                gauss=element->NewGaussBase(5);
    38053918                for(int ig=gauss->begin();ig<gauss->end();ig++){
    38063919                        gauss->GaussPoint(ig);
    … …  
    38113924
    38123925                        for(i=0;i<numnodes;i++){
    3813                                 pe->values[i*dim+0] += - sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
    3814                                 pe->values[i*dim+1] += - sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
    3815                                 if(dim==3){
    3816                                         pe->values[i*dim+2]+= - sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
    3817                                 }
     3926                                pe->values[i*dim+0] += + sigmann*bed_normal[0]*gauss->weight*Jdet*vbasis[i];
     3927                                pe->values[i*dim+1] += + sigmann*bed_normal[1]*gauss->weight*Jdet*vbasis[i];
     3928                                if(dim==3) pe->values[i*dim+2] += + sigmann*bed_normal[2]*gauss->weight*Jdet*vbasis[i];
    38183929                        }
    38193930                }
    … …  
    38233934
    38243935        /*Transform coordinate system*/
    3825         //element->TransformLoadVectorCoord(pe,cs_list); Do not transform pressure augmentation
     3936        //element->TransformLoadVectorCoord(pe,cs_list); Do not transform augmentation
    38263937
    38273938        /*Clean up and return*/
    … …  
    44714582
    44724583        /*Fetch number of nodes for this finite element*/
    4473         int lnumnodes;
    4474         if(dim==2) lnumnodes=3;
    4475         else lnumnodes=6;
    4476         //int pnumnodes = element->NumberofNodes(P1Enum);
     4584        int lnumnodes = element->GetNumberOfNodes(P2Enum);
    44774585
    44784586        /*Get nodal functions derivatives*/
    44794587        IssmDouble* basis =xNew<IssmDouble>(lnumnodes);
    4480         element->NodalFunctions(basis,gauss);
     4588        element->NodalFunctionsP2(basis,gauss);
    44814589
    44824590        /*Build B: */
    44834591        for(int i=0;i<lnumnodes;i++){
    4484                 C[i*lnumnodes+0] = basis[i];
    4485                 C[i*lnumnodes+1] = basis[i];
     4592                C[lnumnodes*0+i] = basis[i];
     4593                C[lnumnodes*1+i] = basis[i];
    44864594        }
    44874595
    … …  
    45024610
    45034611        /*Fetch number of nodes for this finite element*/
    4504         int vnumnodes = element->GetNumberOfNodes();
     4612        int vnumnodes = element->NumberofNodesVelocity();
    45054613        int vnumdof   = vnumnodes*dim;
    45064614

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

    r18227 r18284  
    7979                ElementVector* CreatePVectorFSShelf(Element* element);
    8080                ElementVector* CreatePVectorFSFront(Element* element);
     81                ElementVector* CreatePVectorFSFriction(Element* element);
     82                ElementVector* CreatePVectorFSStress(Element* element);
    8183                void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
    8284                void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
Note: See TracChangeset for help on using the changeset viewer.