Changeset 4838


Ignore:
Timestamp:
07/27/10 15:20:20 (15 years ago)
Author:
seroussi
Message:

added GetInput for Penta

Location:
issm/trunk/src/c/objects/Elements
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/objects/Elements/Penta.cpp

    r4835 r4838  
    30243024        int          numberofdofspernode;
    30253025
    3026 
    30273026        /* 3d gaussian points: */
    30283027        int     num_gauss,ig;
    … …  
    30973096        bool shelf;
    30983097
     3098        /*inputs: */
     3099        Input* vx_input=NULL;
     3100        Input* vy_input=NULL;
     3101        Input* vxold_input=NULL;
     3102        Input* vyold_input=NULL;
     3103
    30993104        /*retrieve inputs :*/
    31003105        inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
    … …  
    31123117          bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build
    31133118          the stiffness matrix. */
    3114 
    31153119
    31163120        if ((collapse==1) && (onbed==0)){
    … …  
    31363140                GetDofList(&doflist[0],&numberofdofspernode);
    31373141
    3138 
    31393142                /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    31403143                  get tria gaussian points as well as segment gaussian points. For tria gaussian
    … …  
    31483151                GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    31493152
     3153                /*Retrieve all inputs we will be needing: */
     3154                vx_input=inputs->GetInput(VxEnum);
     3155                vy_input=inputs->GetInput(VyEnum);
     3156                vxold_input=inputs->GetInput(VxOldEnum);
     3157                vyold_input=inputs->GetInput(VyOldEnum);
     3158
    31503159                /* Start  looping on the number of gaussian points: */
    31513160                for (ig1=0; ig1<num_area_gauss; ig1++){
    … …  
    31653174
    31663175                                /*Get strain rate from velocity: */
    3167                                 inputs->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum);
    3168                                 inputs->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum);
     3176                                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input);
     3177                                this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input);
    31693178
    31703179                                /*Get viscosity: */
    … …  
    33443353        bool isstokes;
    33453354
     3355        /*inputs: */
     3356        Input* vx_input=NULL;
     3357        Input* vy_input=NULL;
     3358        Input* vz_input=NULL;
     3359
    33463360        /*retrive parameters: */
    33473361        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    … …  
    33793393        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    33803394        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     3395
     3396        /*Retrieve all inputs we will be needing: */
     3397        vx_input=inputs->GetInput(VxEnum);
     3398        vy_input=inputs->GetInput(VyEnum);
     3399        vz_input=inputs->GetInput(VzEnum);
    33813400
    33823401        /* Start  looping on the number of gaussian points: */
    … …  
    33933412
    33943413                        /*Compute strain rate: */
    3395                         inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
     3414                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    33963415
    33973416                        /*Get viscosity: */
    … …  
    34623481
    34633482                        /*Compute strain rate: */
    3464                         inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
     3483                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    34653484
    34663485                        /*Get viscosity at last iteration: */
    … …  
    40944113        bool collapse;
    40954114        bool onbed;
     4115        Input* surface_input=NULL;
     4116        Input* thickness_input=NULL;
    40964117
    40974118        /*retrieve inputs :*/
    … …  
    41354156                GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    41364157
     4158                /*Retrieve all inputs we will be needing: */
     4159                thickness_input=inputs->GetInput(ThicknessEnum);
     4160                surface_input=inputs->GetInput(SurfaceEnum);
     4161
    41374162                /* Start  looping on the number of gaussian points: */
    41384163                for (ig1=0; ig1<num_area_gauss; ig1++){
    … …  
    41504175
    41514176                                /*Compute thickness at gaussian point: */
    4152                                 inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum);
     4177                                thickness_input->GetParameterValue(&thickness, gauss_coord);
    41534178
    41544179                                /*Compute slope at gaussian point: */
    4155                                 inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum);
     4180                                surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord);
    41564181
    41574182                                /* Get Jacobian determinant: */
    … …  
    42934318        bool shelf;
    42944319        bool isstokes;
     4320        Input* vx_input=NULL;
     4321        Input* vy_input=NULL;
     4322        Input* vz_input=NULL;
     4323        Input* bed_input=NULL;
    42954324
    42964325        /*retrieve inputs :*/
    … …  
    43264355        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    43274356        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     4357
     4358        /*Retrieve all inputs we will be needing: */
     4359        vx_input=inputs->GetInput(VxEnum);
     4360        vy_input=inputs->GetInput(VyEnum);
     4361        vz_input=inputs->GetInput(VzEnum);
     4362        bed_input=inputs->GetInput(BedEnum);
    43284363
    43294364        /* Start  looping on the number of gaussian points: */
    … …  
    43404375
    43414376                        /*Compute strain rate and viscosity: */
    4342                         inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
     4377                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    43434378                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    43444379
    … …  
    44214456
    44224457                        /* Get bed at gaussian point */
    4423                         inputs->GetParameterValue(&bed, gauss_coord,BedEnum);
     4458                        bed_input->GetParameterValue(&bed, gauss_coord);
    44244459
    44254460                        /*Get L matrix : */
    … …  
    45144549        bool onwater;
    45154550        bool onbed;
     4551        Input* vx_input=NULL;
     4552        Input* vy_input=NULL;
    45164553
    45174554        /*retrieve inputs :*/
    … …  
    45414578        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    45424579
     4580        /*Retrieve all inputs we will be needing: */
     4581        vx_input=inputs->GetInput(VxEnum);
     4582        vy_input=inputs->GetInput(VyEnum);
     4583
    45434584        /* Start  looping on the number of gaussian points: */
    45444585        for (ig1=0; ig1<num_area_gauss; ig1++){
    … …  
    45574598                        /*Get velocity derivative, with respect to x and y: */
    45584599
    4559                         inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum);
    4560                         inputs->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord,VyEnum);
     4600                        vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord);
     4601                        vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord);
    45614602                        dudx=du[0];
    45624603                        dvdy=dv[1];
    … …  
    47254766        bool onbed;
    47264767        bool shelf;
     4768        Input* vx_input=NULL;
     4769        Input* vy_input=NULL;
     4770        Input* vz_input=NULL;
    47274771
    47284772        /*retrieve inputs :*/
    … …  
    47574801        /*Get gaussian points: */
    47584802        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
     4803
     4804        /*Retrieve all inputs we will be needing: */
     4805        vx_input=inputs->GetInput(VxEnum);
     4806        vy_input=inputs->GetInput(VyEnum);
     4807        vz_input=inputs->GetInput(VzEnum);
    47594808
    47604809        /* Start  looping on the number of gaussian points: */
    … …  
    47714820
    47724821                        /*Compute strain rate and viscosity: */
    4773                         inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);
     4822                        this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input);
    47744823                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    47754824
    … …  
    60146063}
    60156064/*}}}*/
     6065/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
     6066void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){
     6067        /*Compute the 3d Blatter/PattynStrain Rate (5 components):
     6068         *
     6069         * epsilon=[exx eyy exy exz eyz]
     6070         *
     6071         * with exz=1/2 du/dz
     6072         *      eyz=1/2 dv/dz
     6073         *
     6074         * the contribution of vz is neglected
     6075         */
     6076
     6077        int i;
     6078
     6079        double epsilonvx[5];
     6080        double epsilonvy[5];
     6081
     6082        /*Check that both inputs have been found*/
     6083        if (!vx_input || !vy_input){
     6084                ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
     6085        }
     6086
     6087        /*Get strain rate assuming that epsilon has been allocated*/
     6088        vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
     6089        vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     6090
     6091        /*Sum all contributions*/
     6092        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     6093
     6094}
     6095/*}}}*/
     6096/*FUNCTION Penta::GetStrainRate3d{{{1*/
     6097void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
     6098        /*Compute the 3d Strain Rate (6 components):
     6099         *
     6100         * epsilon=[exx eyy ezz exy exz eyz]
     6101         */
     6102
     6103        int i;
     6104
     6105        double epsilonvx[6];
     6106        double epsilonvy[6];
     6107        double epsilonvz[6];
     6108
     6109        /*Check that both inputs have been found*/
     6110        if (!vx_input || !vy_input || !vz_input){
     6111                ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input);
     6112        }
     6113
     6114        /*Get strain rate assuming that epsilon has been allocated*/
     6115        vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
     6116        vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
     6117        vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
     6118
     6119        /*Sum all contributions*/
     6120        for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
     6121
     6122}
     6123/*}}}*/
    60166124/*FUNCTION Penta::GradjB {{{1*/
    60176125void  Penta::GradjB(Vec gradient){

  • issm/trunk/src/c/objects/Elements/Penta.h

    r4835 r4838  
    168168                void      GetSolutionFromInputsDiagnosticVert(Vec solutiong);
    169169                void      GetSolutionFromInputsThermal(Vec solutiong);
    170                 void      GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
    171                 void      GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);
     170                void    GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
     171                void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
    172172                Penta*  GetUpperElement(void);
    173173                void      InputExtrude(int enum_type,bool only_if_collapsed);
Note: See TracChangeset for help on using the changeset viewer.