Changeset 8654


Ignore:
Timestamp:
06/17/11 09:01:02 (14 years ago)
Author:
Mathieu Morlighem
Message:

moved some Penta functions from Tria to Penta

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

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

    r8653 r8654  
    38743874}
    38753875/*}}}*/
     3876/*FUNCTION Penta::GetBasalElement{{{1*/
     3877Penta* Penta::GetBasalElement(void){
     3878
     3879        /*Output*/
     3880        Penta* penta=NULL;
     3881
     3882        /*Go through all elements till the bed is reached*/
     3883        penta=this;
     3884        for(;;){
     3885                /*Stop if we have reached the surface, else, take lower penta*/
     3886                if (penta->IsOnBed()) break;
     3887
     3888                /* get lower Penta*/
     3889                penta=penta->GetLowerElement();
     3890                _assert_(penta->Id()!=this->id);
     3891        }
     3892
     3893        /*return output*/
     3894        return penta;
     3895}
     3896/*}}}*/
     3897/*FUNCTION Penta::GetDofList {{{1*/
     3898void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
     3899
     3900        int  i,j,count=0;
     3901        int  numberofdofs=0;
     3902        int* doflist=NULL;
     3903
     3904        /*First, figure out size of doflist: */
     3905        for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     3906
     3907        /*Allocate: */
     3908        doflist=(int*)xmalloc(numberofdofs*sizeof(int));
     3909
     3910        /*Populate: */
     3911        count=0;
     3912        for(i=0;i<6;i++){
     3913                nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
     3914                count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
     3915        }
     3916
     3917        /*Assign output pointers:*/
     3918        *pdoflist=doflist;
     3919}
     3920/*}}}*/
     3921/*FUNCTION Penta::GetDofList1 {{{1*/
     3922void  Penta::GetDofList1(int* doflist){
     3923
     3924        int i;
     3925        for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1();
     3926
     3927}
     3928/*}}}*/
     3929/*FUNCTION Penta::GetElementType {{{1*/
     3930int Penta::GetElementType(){
     3931
     3932        /*return PentaRef field*/
     3933        return this->element_type;
     3934}
     3935/*}}}*/
     3936/*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/
     3937int* Penta::GetHorizontalNeighboorSids(){
     3938
     3939        /*return PentaRef field*/
     3940        return &this->horizontalneighborsids[0];
     3941
     3942}
     3943/*}}}*/
     3944/*FUNCTION Penta::GetLowerElement{{{1*/
     3945Penta* Penta::GetLowerElement(void){
     3946
     3947        Penta* upper_penta=NULL;
     3948
     3949        upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
     3950
     3951        return upper_penta;
     3952}
     3953/*}}}*/
     3954/*FUNCTION Penta::GetNodeIndex {{{1*/
     3955int Penta::GetNodeIndex(Node* node){
     3956
     3957        _assert_(nodes);
     3958        for(int i=0;i<NUMVERTICES;i++){
     3959                if(node==nodes[i])
     3960                 return i;
     3961        }
     3962        _error_("Node provided not found among element nodes");
     3963
     3964}
     3965/*}}}*/
     3966/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
     3967void Penta::GetParameterListOnVertices(double* pvalue,int enumtype){
     3968
     3969        /*Intermediaries*/
     3970        double     value[NUMVERTICES];
     3971        GaussPenta *gauss              = NULL;
     3972
     3973        /*Recover input*/
     3974        Input* input=inputs->GetInput(enumtype);
     3975        if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));
     3976
     3977        /*Checks in debugging mode*/
     3978        _assert_(pvalue);
     3979
     3980        /* Start looping on the number of vertices: */
     3981        gauss=new GaussPenta();
     3982        for (int iv=0;iv<NUMVERTICES;iv++){
     3983                gauss->GaussVertex(iv);
     3984                input->GetParameterValue(&pvalue[iv],gauss);
     3985        }
     3986
     3987        /*clean-up*/
     3988        delete gauss;
     3989}
     3990/*}}}*/
     3991/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
     3992void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
     3993
     3994        /*Intermediaries*/
     3995        double     value[NUMVERTICES];
     3996        GaussPenta *gauss              = NULL;
     3997
     3998        /*Recover input*/
     3999        Input* input=inputs->GetInput(enumtype);
     4000
     4001        /*Checks in debugging mode*/
     4002        _assert_(pvalue);
     4003
     4004        /* Start looping on the number of vertices: */
     4005        if (input){
     4006                gauss=new GaussPenta();
     4007                for (int iv=0;iv<NUMVERTICES;iv++){
     4008                        gauss->GaussVertex(iv);
     4009                        input->GetParameterValue(&pvalue[iv],gauss);
     4010                }
     4011        }
     4012        else{
     4013                for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
     4014        }
     4015
     4016        /*clean-up*/
     4017        delete gauss;
     4018}
     4019/*}}}*/
     4020/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
     4021void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){
     4022
     4023        Input* input=inputs->GetInput(enumtype);
     4024        if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));
     4025
     4026        GaussPenta* gauss=new GaussPenta();
     4027        gauss->GaussVertex(this->GetNodeIndex(node));
     4028
     4029        input->GetParameterValue(pvalue,gauss);
     4030        delete gauss;
     4031}
     4032/*}}}*/
     4033/*FUNCTION Penta::GetPhi {{{1*/
     4034void Penta::GetPhi(double* phi, double*  epsilon, double viscosity){
     4035        /*Compute deformational heating from epsilon and viscosity */
     4036
     4037        double epsilon_matrix[3][3];
     4038        double epsilon_eff;
     4039        double epsilon_sqr[3][3];
     4040
     4041        /* Build epsilon matrix */
     4042        epsilon_matrix[0][0]=*(epsilon+0);
     4043        epsilon_matrix[1][0]=*(epsilon+3);
     4044        epsilon_matrix[2][0]=*(epsilon+4);
     4045        epsilon_matrix[0][1]=*(epsilon+3);
     4046        epsilon_matrix[1][1]=*(epsilon+1);
     4047        epsilon_matrix[2][1]=*(epsilon+5);
     4048        epsilon_matrix[0][2]=*(epsilon+4);
     4049        epsilon_matrix[1][2]=*(epsilon+5);
     4050        epsilon_matrix[2][2]=*(epsilon+2);
     4051
     4052        /* Effective value of epsilon_matrix */
     4053        epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2);
     4054        epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2);
     4055        epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2);
     4056        epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2);
     4057        epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2);
     4058        epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2);
     4059        epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2);
     4060        epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2);
     4061        epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);
     4062        epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5);
     4063
     4064        /*Phi = Tr(sigma * eps)
     4065         *    = Tr(sigma'* eps)
     4066         *    = 2 * eps_eff * sigma'_eff
     4067         *    = 4 * mu * eps_eff ^2*/
     4068        *phi=4*pow(epsilon_eff,2.0)*viscosity;
     4069}
     4070/*}}}*/
     4071/*FUNCTION Penta::GetSidList{{{1*/
     4072void  Penta::GetSidList(int* sidlist){
     4073
     4074        int i;
     4075        for(i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
     4076
     4077}
     4078/*}}}*/
     4079/*FUNCTION Penta::GetSolutionFromInputs{{{1*/
     4080void  Penta::GetSolutionFromInputs(Vec solution){
     4081
     4082        int analysis_type;
     4083
     4084        /*retrive parameters: */
     4085        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4086
     4087        /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
     4088        if (analysis_type==DiagnosticHorizAnalysisEnum){
     4089                int approximation;
     4090                inputs->GetParameterValue(&approximation,ApproximationEnum);
     4091                if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
     4092                        GetSolutionFromInputsDiagnosticStokes(solution);
     4093                }
     4094                else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
     4095                        GetSolutionFromInputsDiagnosticHoriz(solution);
     4096                }
     4097                else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
     4098                        return; //the elements around will create the solution
     4099                }
     4100        }
     4101        else if(analysis_type==DiagnosticHutterAnalysisEnum){
     4102                GetSolutionFromInputsDiagnosticHutter(solution);
     4103        }
     4104        else if(analysis_type==DiagnosticVertAnalysisEnum){
     4105                GetSolutionFromInputsDiagnosticVert(solution);
     4106        }
     4107        else if(analysis_type==ThermalAnalysisEnum){
     4108                GetSolutionFromInputsThermal(solution);
     4109        }
     4110        else if(analysis_type==EnthalpyAnalysisEnum){
     4111                GetSolutionFromInputsEnthalpy(solution);
     4112        }
     4113        else{
     4114                _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
     4115        }
     4116}
     4117/*}}}*/
     4118/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
     4119void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
     4120
     4121        const int    numdof=NDOF2*NUMVERTICES;
     4122
     4123        int          i;
     4124        int          approximation;
     4125        int*         doflist=NULL;
     4126        double       vx,vy;
     4127        double       values[numdof];
     4128        GaussPenta*  gauss;
     4129
     4130        /*Get approximation enum and dof list: */
     4131        inputs->GetParameterValue(&approximation,ApproximationEnum);
     4132        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     4133        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     4134
     4135        /*If the element is a coupling, do nothing: every node is also on an other elements
     4136         * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
     4137        GetDofList(&doflist,approximation,GsetEnum);
     4138
     4139        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     4140        /*P1 element only for now*/
     4141        gauss=new GaussPenta();
     4142        for(i=0;i<NUMVERTICES;i++){
     4143
     4144                /*Recover vx and vy*/
     4145                gauss->GaussVertex(i);
     4146                vx_input->GetParameterValue(&vx,gauss);
     4147                vy_input->GetParameterValue(&vy,gauss);
     4148                values[i*NDOF2+0]=vx;
     4149                values[i*NDOF2+1]=vy;
     4150        }
     4151
     4152        /*Add value to global vector*/
     4153        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4154
     4155        /*Free ressources:*/
     4156        delete gauss;
     4157        xfree((void**)&doflist);
     4158}
     4159/*}}}*/
     4160/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/
     4161void  Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){
     4162
     4163        const int    numdof=NDOF2*NUMVERTICES;
     4164
     4165        int          i;
     4166        int*         doflist=NULL;
     4167        double       vx,vy;
     4168        double       values[numdof];
     4169        GaussPenta*  gauss=NULL;
     4170
     4171        /*Get dof list: */
     4172        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4173        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     4174        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     4175
     4176        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     4177        /*P1 element only for now*/
     4178        gauss=new GaussPenta();
     4179        for(i=0;i<NUMVERTICES;i++){
     4180                /*Recover vx and vy*/
     4181                gauss->GaussVertex(i);
     4182                vx_input->GetParameterValue(&vx,gauss);
     4183                vy_input->GetParameterValue(&vy,gauss);
     4184                values[i*NDOF2+0]=vx;
     4185                values[i*NDOF2+1]=vy;
     4186        }
     4187
     4188        /*Add value to global vector*/
     4189        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4190
     4191        /*Free ressources:*/
     4192        delete gauss;
     4193        xfree((void**)&doflist);
     4194}
     4195/*}}}*/
     4196/*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/
     4197void  Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){
     4198
     4199        const int    numdof=NDOF1*NUMVERTICES;
     4200
     4201        int          i;
     4202        int*         doflist=NULL;
     4203        double       vz;
     4204        double       values[numdof];
     4205        GaussPenta*  gauss=NULL;
     4206
     4207        /*Get dof list: */
     4208        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4209        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     4210
     4211        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     4212        /*P1 element only for now*/
     4213        gauss=new GaussPenta();
     4214        for(i=0;i<NUMVERTICES;i++){
     4215                /*Recover vz */
     4216                gauss->GaussVertex(i);
     4217                vz_input->GetParameterValue(&vz,gauss);
     4218                values[i]=vz;
     4219        }
     4220
     4221        /*Add value to global vector*/
     4222        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4223
     4224        /*Free ressources:*/
     4225        delete gauss;
     4226        xfree((void**)&doflist);
     4227}
     4228/*}}}*/
     4229/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/
     4230void  Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){
     4231
     4232        const int    numdof=NDOF4*NUMVERTICES;
     4233
     4234        int          i;
     4235        int*         doflist=NULL;
     4236        double       vx,vy,vz,p;
     4237        double       stokesreconditioning;
     4238        double       values[numdof];
     4239        GaussPenta   *gauss;
     4240
     4241        /*Get dof list: */
     4242        GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
     4243        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     4244        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     4245        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
     4246        Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input);
     4247
     4248        /*Recondition pressure: */
     4249        this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
     4250
     4251        /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
     4252        /*P1 element only for now*/
     4253        gauss=new GaussPenta();
     4254        for(i=0;i<NUMVERTICES;i++){
     4255                gauss->GaussVertex(i);
     4256                vx_input->GetParameterValue(&vx,gauss);
     4257                vy_input->GetParameterValue(&vy,gauss);
     4258                vz_input->GetParameterValue(&vz,gauss);
     4259                p_input ->GetParameterValue(&p ,gauss);
     4260                values[i*NDOF4+0]=vx;
     4261                values[i*NDOF4+1]=vy;
     4262                values[i*NDOF4+2]=vz;
     4263                values[i*NDOF4+3]=p/stokesreconditioning;
     4264        }
     4265
     4266        /*Add value to global vector*/
     4267        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4268
     4269        /*Free ressources:*/
     4270        delete gauss;
     4271        xfree((void**)&doflist);
     4272}
     4273/*}}}*/
     4274/*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/
     4275void  Penta::GetSolutionFromInputsThermal(Vec solution){
     4276
     4277        const int    numdof=NDOF1*NUMVERTICES;
     4278
     4279        int          i;
     4280        int*         doflist=NULL;
     4281        double       values[numdof];
     4282        double       temp;
     4283        GaussPenta   *gauss=NULL;
     4284
     4285        /*Get dof list: */
     4286        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4287        Input* t_input=inputs->GetInput(TemperatureEnum); _assert_(t_input);
     4288
     4289        gauss=new GaussPenta();
     4290        for(i=0;i<NUMVERTICES;i++){
     4291                /*Recover temperature*/
     4292                gauss->GaussVertex(i);
     4293                t_input->GetParameterValue(&temp,gauss);
     4294                values[i]=temp;
     4295        }
     4296
     4297        /*Add value to global vector*/
     4298        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4299
     4300        /*Free ressources:*/
     4301        delete gauss;
     4302        xfree((void**)&doflist);
     4303}
     4304/*}}}*/
     4305/*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/
     4306void  Penta::GetSolutionFromInputsEnthalpy(Vec solution){
     4307
     4308        const int    numdof=NDOF1*NUMVERTICES;
     4309
     4310        int          i;
     4311        int*         doflist=NULL;
     4312        double       values[numdof];
     4313        double       enthalpy;
     4314        GaussPenta   *gauss=NULL;
     4315
     4316        /*Get dof list: */
     4317        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4318        Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input);
     4319
     4320        gauss=new GaussPenta();
     4321        for(i=0;i<NUMVERTICES;i++){
     4322                /*Recover temperature*/
     4323                gauss->GaussVertex(i);
     4324                h_input->GetParameterValue(&enthalpy,gauss);
     4325                values[i]=enthalpy;
     4326        }
     4327
     4328        /*Add value to global vector*/
     4329        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     4330
     4331        /*Free ressources:*/
     4332        delete gauss;
     4333        xfree((void**)&doflist);
     4334}
     4335/*}}}*/
     4336/*FUNCTION Penta::GetStabilizationParameter {{{1*/
     4337double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
     4338        /*Compute stabilization parameter*/
     4339
     4340        double normu;
     4341        double tau_parameter;
     4342
     4343        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
     4344        if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
     4345                tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
     4346        }
     4347        else tau_parameter=diameter/(2*normu);
     4348
     4349        return tau_parameter;
     4350}
     4351/*}}}*/
     4352/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
     4353void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
     4354        /*Compute the 3d Blatter/PattynStrain Rate (5 components):
     4355         *
     4356         * epsilon=[exx eyy exy exz eyz]
     4357         *
     4358         * with exz=1/2 du/dz
     4359         *      eyz=1/2 dv/dz
     4360         *
     4361         * the contribution of vz is neglected
     4362         */
     4363
     4364        int i;
     4365        double epsilonvx[5];
     4366        double epsilonvy[5];
     4367
     4368        /*Check that both inputs have been found*/
     4369        if (!vx_input || !vy_input){
     4370                _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
     4371        }
     4372
     4373        /*Get strain rate assuming that epsilon has been allocated*/
     4374        vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
     4375        vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
     4376
     4377        /*Sum all contributions*/
     4378        for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
     4379}
     4380/*}}}*/
     4381/*FUNCTION Penta::GetStrainRate3d{{{1*/
     4382void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
     4383        /*Compute the 3d Strain Rate (6 components):
     4384         *
     4385         * epsilon=[exx eyy ezz exy exz eyz]
     4386         */
     4387
     4388        int i;
     4389        double epsilonvx[6];
     4390        double epsilonvy[6];
     4391        double epsilonvz[6];
     4392
     4393        /*Check that both inputs have been found*/
     4394        if (!vx_input || !vy_input || !vz_input){
     4395                _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input);
     4396        }
     4397
     4398        /*Get strain rate assuming that epsilon has been allocated*/
     4399        vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
     4400        vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
     4401        vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
     4402
     4403        /*Sum all contributions*/
     4404        for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
     4405}
     4406/*}}}*/
     4407/*FUNCTION Penta::GetUpperElement{{{1*/
     4408Penta* Penta::GetUpperElement(void){
     4409
     4410        Penta* upper_penta=NULL;
     4411
     4412        upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
     4413
     4414        return upper_penta;
     4415}
     4416/*}}}*/
     4417/*FUNCTION Penta::GetVectorFromInputs{{{1*/
     4418void  Penta::GetVectorFromInputs(Vec vector,int input_enum){
     4419
     4420        int doflist1[NUMVERTICES];
     4421
     4422        /*Get out if this is not an element input*/
     4423        if (!IsInput(input_enum)) return;
     4424
     4425        /*Prepare index list*/
     4426        this->GetDofList1(&doflist1[0]);
     4427
     4428        /*Get input (either in element or material)*/
     4429        Input* input=inputs->GetInput(input_enum);
     4430        if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
     4431
     4432        /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
     4433        input->GetVectorFromInputs(vector,&doflist1[0]);
     4434}
     4435/*}}}*/
     4436/*FUNCTION Penta::GetZcoord {{{1*/
     4437double Penta::GetZcoord(GaussPenta* gauss){
     4438
     4439        int    i;
     4440        double z;
     4441        double xyz_list[NUMVERTICES][3];
     4442        double z_list[NUMVERTICES];
     4443
     4444        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4445        for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
     4446        PentaRef::GetParameterValue(&z,z_list,gauss);
     4447
     4448        return z;
     4449}
     4450/*}}}*/
    38764451/*FUNCTION Penta::Gradj {{{1*/
    38774452void  Penta::Gradj(Vec gradient,int control_type){
    … …  
    38884463
    38894464                case DragCoefficientEnum:
    3890 
    38914465                        inputs->GetParameterValue(&approximation,ApproximationEnum);
    3892                         if(IsOnShelf() || !IsOnBed()) return;
    3893 
    38944466                        switch(approximation){
    38954467                                case MacAyealApproximationEnum:
    3896                                         tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3897                                         tria->GradjDragMacAyeal(gradient);
    3898                                         delete tria->matice; delete tria;
     4468                                        GradjDragMacAyeal(gradient);
    38994469                                        break;
    39004470                                case PattynApproximationEnum:
    3901                                         tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3902                                         tria->GradjDragMacAyeal(gradient);
    3903                                         delete tria->matice; delete tria;
    3904                                         //GradjDragPattyn(gradient);
     4471                                        GradjDragPattyn(gradient);
    39054472                                        break;
    39064473                                case StokesApproximationEnum:
    3907                                         tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
    3908                                         tria->GradjDragStokes(gradient);
    3909                                         delete tria->matice; delete tria;
    3910                                         //GradjDragStokes(gradient);
     4474                                        GradjDragStokes(gradient);
    39114475                                        break;
    39124476                                case NoneApproximationEnum:
    … …  
    39214485                        inputs->GetParameterValue(&approximation,ApproximationEnum);
    39224486                        switch(approximation){
    3923 
    39244487                                case MacAyealApproximationEnum:
    3925 
    3926                                         /*This element should be collapsed into a tria element at its base*/
    3927                                         if (!IsOnBed()) return;
    3928 
    3929                                         /*Depth Average B*/
    3930                                         this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    3931 
    3932                                         tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    3933                                         tria->GradjBMacAyeal(gradient);
    3934                                         delete tria->matice; delete tria;
    3935 
    3936                                         /*delete B average*/
    3937                                         this->matice->inputs->DeleteInput(RheologyBbarEnum);
     4488                                        GradjBbarMacAyeal(gradient);
    39384489                                        break;
    3939                                 case PattynApproximationEnum: case StokesApproximationEnum:
    3940 
    3941                                         /*Gradient is computed on bed only (Bbar)*/
    3942                                         if (!IsOnBed()) return;
    3943 
    3944                                         /*Depth Average B*/
    3945                                         this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
    3946 
    3947                                         /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
    3948                                         tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
    3949                                         tria->GradjBMacAyeal(gradient);
    3950                                         delete tria->matice; delete tria;
    3951 
    3952                                         /*delete B average*/
    3953                                         this->matice->inputs->DeleteInput(RheologyBbarEnum);
     4490                                case PattynApproximationEnum:
     4491                                        GradjBbarPattyn(gradient);
    39544492                                        break;
    3955 
     4493                                case StokesApproximationEnum:
     4494                                        GradjBbarStokes(gradient);
     4495                                        break;
    39564496                                case NoneApproximationEnum:
    39574497                                        /*Gradient is 0*/
    … …  
    39994539}
    40004540/*}}}*/
    4001 /*FUNCTION Penta::GetBasalElement{{{1*/
    4002 Penta* Penta::GetBasalElement(void){
    4003 
    4004         /*Output*/
    4005         Penta* penta=NULL;
    4006 
    4007         /*Go through all elements till the bed is reached*/
    4008         penta=this;
    4009         for(;;){
    4010                 /*Stop if we have reached the surface, else, take lower penta*/
    4011                 if (penta->IsOnBed()) break;
    4012 
    4013                 /* get lower Penta*/
    4014                 penta=penta->GetLowerElement();
    4015                 _assert_(penta->Id()!=this->id);
    4016         }
    4017 
    4018         /*return output*/
    4019         return penta;
    4020 }
    4021 /*}}}*/
    4022 /*FUNCTION Penta::GetDofList {{{1*/
    4023 void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
    4024 
    4025         int  i,j,count=0;
    4026         int  numberofdofs=0;
    4027         int* doflist=NULL;
    4028 
    4029         /*First, figure out size of doflist: */
    4030         for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    4031 
    4032         /*Allocate: */
    4033         doflist=(int*)xmalloc(numberofdofs*sizeof(int));
    4034 
    4035         /*Populate: */
    4036         count=0;
    4037         for(i=0;i<6;i++){
    4038                 nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
    4039                 count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
    4040         }
    4041 
    4042         /*Assign output pointers:*/
    4043         *pdoflist=doflist;
    4044 }
    4045 /*}}}*/
    4046 /*FUNCTION Penta::GetDofList1 {{{1*/
    4047 void  Penta::GetDofList1(int* doflist){
    4048 
    4049         int i;
    4050         for(i=0;i<6;i++) doflist[i]=nodes[i]->GetDofList1();
    4051 
    4052 }
    4053 /*}}}*/
    4054 /*FUNCTION Penta::GetElementType {{{1*/
    4055 int Penta::GetElementType(){
    4056 
    4057         /*return PentaRef field*/
    4058         return this->element_type;
    4059 }
    4060 /*}}}*/
    4061 /*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/
    4062 int* Penta::GetHorizontalNeighboorSids(){
    4063 
    4064         /*return PentaRef field*/
    4065         return &this->horizontalneighborsids[0];
    4066 
    4067 }
    4068 /*}}}*/
    4069 /*FUNCTION Penta::GetLowerElement{{{1*/
    4070 Penta* Penta::GetLowerElement(void){
    4071 
    4072         Penta* upper_penta=NULL;
    4073 
    4074         upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
    4075 
    4076         return upper_penta;
    4077 }
    4078 /*}}}*/
    4079 /*FUNCTION Penta::GetNodeIndex {{{1*/
    4080 int Penta::GetNodeIndex(Node* node){
    4081 
    4082         _assert_(nodes);
    4083         for(int i=0;i<NUMVERTICES;i++){
    4084                 if(node==nodes[i])
    4085                  return i;
    4086         }
    4087         _error_("Node provided not found among element nodes");
    4088 
    4089 }
    4090 /*}}}*/
    4091 /*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
    4092 void Penta::GetParameterListOnVertices(double* pvalue,int enumtype){
    4093 
    4094         /*Intermediaries*/
    4095         double     value[NUMVERTICES];
    4096         GaussPenta *gauss              = NULL;
    4097 
    4098         /*Recover input*/
    4099         Input* input=inputs->GetInput(enumtype);
    4100         if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));
    4101 
    4102         /*Checks in debugging mode*/
    4103         _assert_(pvalue);
    4104 
    4105         /* Start looping on the number of vertices: */
    4106         gauss=new GaussPenta();
    4107         for (int iv=0;iv<NUMVERTICES;iv++){
    4108                 gauss->GaussVertex(iv);
    4109                 input->GetParameterValue(&pvalue[iv],gauss);
    4110         }
    4111 
    4112         /*clean-up*/
    4113         delete gauss;
    4114 }
    4115 /*}}}*/
    4116 /*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
    4117 void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
    4118 
    4119         /*Intermediaries*/
    4120         double     value[NUMVERTICES];
    4121         GaussPenta *gauss              = NULL;
    4122 
    4123         /*Recover input*/
    4124         Input* input=inputs->GetInput(enumtype);
    4125 
    4126         /*Checks in debugging mode*/
    4127         _assert_(pvalue);
    4128 
    4129         /* Start looping on the number of vertices: */
    4130         if (input){
    4131                 gauss=new GaussPenta();
    4132                 for (int iv=0;iv<NUMVERTICES;iv++){
    4133                         gauss->GaussVertex(iv);
    4134                         input->GetParameterValue(&pvalue[iv],gauss);
    4135                 }
    4136         }
    4137         else{
    4138                 for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
    4139         }
    4140 
    4141         /*clean-up*/
    4142         delete gauss;
    4143 }
    4144 /*}}}*/
    4145 /*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
    4146 void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){
    4147 
    4148         Input* input=inputs->GetInput(enumtype);
    4149         if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));
    4150 
    4151         GaussPenta* gauss=new GaussPenta();
    4152         gauss->GaussVertex(this->GetNodeIndex(node));
    4153 
    4154         input->GetParameterValue(pvalue,gauss);
    4155         delete gauss;
    4156 }
    4157 /*}}}*/
    4158 /*FUNCTION Penta::GetPhi {{{1*/
    4159 void Penta::GetPhi(double* phi, double*  epsilon, double viscosity){
    4160         /*Compute deformational heating from epsilon and viscosity */
    4161 
    4162         double epsilon_matrix[3][3];
    4163         double epsilon_eff;
    4164         double epsilon_sqr[3][3];
    4165 
    4166         /* Build epsilon matrix */
    4167         epsilon_matrix[0][0]=*(epsilon+0);
    4168         epsilon_matrix[1][0]=*(epsilon+3);
    4169         epsilon_matrix[2][0]=*(epsilon+4);
    4170         epsilon_matrix[0][1]=*(epsilon+3);
    4171         epsilon_matrix[1][1]=*(epsilon+1);
    4172         epsilon_matrix[2][1]=*(epsilon+5);
    4173         epsilon_matrix[0][2]=*(epsilon+4);
    4174         epsilon_matrix[1][2]=*(epsilon+5);
    4175         epsilon_matrix[2][2]=*(epsilon+2);
    4176 
    4177         /* Effective value of epsilon_matrix */
    4178         epsilon_sqr[0][0]=pow(epsilon_matrix[0][0],2);
    4179         epsilon_sqr[1][0]=pow(epsilon_matrix[1][0],2);
    4180         epsilon_sqr[2][0]=pow(epsilon_matrix[2][0],2);
    4181         epsilon_sqr[0][1]=pow(epsilon_matrix[0][1],2);
    4182         epsilon_sqr[1][1]=pow(epsilon_matrix[1][1],2);
    4183         epsilon_sqr[2][1]=pow(epsilon_matrix[2][1],2);
    4184         epsilon_sqr[0][2]=pow(epsilon_matrix[0][2],2);
    4185         epsilon_sqr[1][2]=pow(epsilon_matrix[1][2],2);
    4186         epsilon_sqr[2][2]=pow(epsilon_matrix[2][2],2);
    4187         epsilon_eff=1/pow(2,0.5)*pow((epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]),0.5);
    4188 
    4189         /*Phi = Tr(sigma * eps)
    4190          *    = Tr(sigma'* eps)
    4191          *    = 2 * eps_eff * sigma'_eff
    4192          *    = 4 * mu * eps_eff ^2*/
    4193         *phi=4*pow(epsilon_eff,2.0)*viscosity;
    4194 }
    4195 /*}}}*/
    4196 /*FUNCTION Penta::GetSidList{{{1*/
    4197 void  Penta::GetSidList(int* sidlist){
    4198 
    4199         int i;
    4200         for(i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
    4201 
    4202 }
    4203 /*}}}*/
    4204 /*FUNCTION Penta::GetSolutionFromInputs{{{1*/
    4205 void  Penta::GetSolutionFromInputs(Vec solution){
    4206 
    4207         int analysis_type;
    4208 
    4209         /*retrive parameters: */
     4541/*FUNCTION Penta::GradjDragMacAyeal {{{1*/
     4542void  Penta::GradjDragMacAyeal(Vec gradient){
     4543
     4544        /*Gradient is 0 if on shelf or not on bed*/
     4545        if(IsOnShelf() || !IsOnBed()) return;
     4546
     4547        /*Spawn tria*/
     4548        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     4549        tria->GradjDragMacAyeal(gradient);
     4550        delete tria->matice; delete tria;
     4551
     4552} /*}}}*/
     4553/*FUNCTION Penta::GradjDragPattyn {{{1*/
     4554void  Penta::GradjDragPattyn(Vec gradient){
     4555
     4556        int        i,j,ig;
     4557        int        drag_type,analysis_type;
     4558        int        doflist1[NUMVERTICES];
     4559        double     vx,vy,lambda,mu,alpha_complement,Jdet;
     4560        double     bed,thickness,Neff,drag;
     4561        double     xyz_list[NUMVERTICES][3];
     4562        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4563        double     dk[NDOF2];
     4564        double     grade_g[NUMVERTICES]={0.0};
     4565        double     grade_g_gaussian[NUMVERTICES];
     4566        double     basis[6];
     4567        Friction*  friction=NULL;
     4568        GaussPenta  *gauss=NULL;
     4569
     4570        /*Gradient is 0 if on shelf or not on bed*/
     4571        if(IsOnShelf() || !IsOnBed()) return;
     4572
     4573        /*Retrieve all inputs and parameters*/
    42104574        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4211 
    4212         /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
    4213         if (analysis_type==DiagnosticHorizAnalysisEnum){
    4214                 int approximation;
    4215                 inputs->GetParameterValue(&approximation,ApproximationEnum);
    4216                 if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
    4217                         GetSolutionFromInputsDiagnosticStokes(solution);
    4218                 }
    4219                 else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
    4220                         GetSolutionFromInputsDiagnosticHoriz(solution);
    4221                 }
    4222                 else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
    4223                         return; //the elements around will create the solution
    4224                 }
    4225         }
    4226         else if(analysis_type==DiagnosticHutterAnalysisEnum){
    4227                 GetSolutionFromInputsDiagnosticHutter(solution);
    4228         }
    4229         else if(analysis_type==DiagnosticVertAnalysisEnum){
    4230                 GetSolutionFromInputsDiagnosticVert(solution);
    4231         }
    4232         else if(analysis_type==ThermalAnalysisEnum){
    4233                 GetSolutionFromInputsThermal(solution);
    4234         }
    4235         else if(analysis_type==EnthalpyAnalysisEnum){
    4236                 GetSolutionFromInputsEnthalpy(solution);
    4237         }
    4238         else{
    4239                 _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
    4240         }
    4241 }
    4242 /*}}}*/
    4243 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    4244 void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
    4245 
    4246         const int    numdof=NDOF2*NUMVERTICES;
    4247 
    4248         int          i;
    4249         int          approximation;
    4250         int*         doflist=NULL;
    4251         double       vx,vy;
    4252         double       values[numdof];
    4253         GaussPenta*  gauss;
    4254 
    4255         /*Get approximation enum and dof list: */
    4256         inputs->GetParameterValue(&approximation,ApproximationEnum);
    4257         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    4258         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    4259 
    4260         /*If the element is a coupling, do nothing: every node is also on an other elements
    4261          * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
    4262         GetDofList(&doflist,approximation,GsetEnum);
    4263 
    4264         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4265         /*P1 element only for now*/
    4266         gauss=new GaussPenta();
    4267         for(i=0;i<NUMVERTICES;i++){
    4268 
    4269                 /*Recover vx and vy*/
    4270                 gauss->GaussVertex(i);
     4575        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4576        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4577        GetDofList1(&doflist1[0]);
     4578        Input* adjointx_input=inputs->GetInput(AdjointxEnum);               _assert_(adjointx_input);
     4579        Input* adjointy_input=inputs->GetInput(AdjointyEnum);               _assert_(adjointy_input);
     4580        Input* vx_input=inputs->GetInput(VxEnum);                           _assert_(vx_input);
     4581        Input* vy_input=inputs->GetInput(VyEnum);                           _assert_(vy_input);
     4582        Input* dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); _assert_(dragcoefficient_input);
     4583
     4584        /*Build frictoin element, needed later: */
     4585        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4586        friction=new Friction("2d",inputs,matpar,analysis_type);
     4587
     4588        /* Start  looping on the number of gaussian points: */
     4589        gauss=new GaussPenta(0,1,2,4);
     4590        for (ig=gauss->begin();ig<gauss->end();ig++){
     4591
     4592                gauss->GaussPoint(ig);
     4593
     4594                GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
     4595                GetNodalFunctionsP1(basis, gauss);
     4596
     4597                /*Build alpha_complement_list: */
     4598                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     4599                else alpha_complement=0;
     4600
     4601                dragcoefficient_input->GetParameterValue(&drag, gauss);
     4602                adjointx_input->GetParameterValue(&lambda, gauss);
     4603                adjointy_input->GetParameterValue(&mu, gauss);
    42714604                vx_input->GetParameterValue(&vx,gauss);
    42724605                vy_input->GetParameterValue(&vy,gauss);
    4273                 values[i*NDOF2+0]=vx;
    4274                 values[i*NDOF2+1]=vy;
    4275         }
    4276 
    4277         /*Add value to global vector*/
    4278         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4279 
    4280         /*Free ressources:*/
     4606                dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     4607
     4608                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     4609                for (i=0;i<NUMVERTICES;i++){
     4610                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i]; /*basis are 0 for the 3 upper nodes*/
     4611                }
     4612
     4613                /*Add gradje_g_gaussian vector to gradje_g: */
     4614                for(i=0;i<NUMVERTICES;i++){
     4615                        _assert_(!isnan(grade_g[i]));
     4616                        grade_g[i]+=grade_g_gaussian[i];
     4617                }
     4618        }
     4619        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     4620
     4621        /*Clean up and return*/
    42814622        delete gauss;
    4282         xfree((void**)&doflist);
    4283 }
    4284 /*}}}*/
    4285 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{1*/
    4286 void  Penta::GetSolutionFromInputsDiagnosticHutter(Vec solution){
    4287 
    4288         const int    numdof=NDOF2*NUMVERTICES;
    4289 
    4290         int          i;
    4291         int*         doflist=NULL;
    4292         double       vx,vy;
    4293         double       values[numdof];
    4294         GaussPenta*  gauss=NULL;
    4295 
    4296         /*Get dof list: */
    4297         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4298         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    4299         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    4300 
    4301         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4302         /*P1 element only for now*/
    4303         gauss=new GaussPenta();
    4304         for(i=0;i<NUMVERTICES;i++){
    4305                 /*Recover vx and vy*/
    4306                 gauss->GaussVertex(i);
    4307                 vx_input->GetParameterValue(&vx,gauss);
    4308                 vy_input->GetParameterValue(&vy,gauss);
    4309                 values[i*NDOF2+0]=vx;
    4310                 values[i*NDOF2+1]=vy;
    4311         }
    4312 
    4313         /*Add value to global vector*/
    4314         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4315 
    4316         /*Free ressources:*/
     4623        delete friction;
     4624}
     4625/*}}}*/
     4626/*FUNCTION Penta::GradjDragStokes {{{1*/
     4627void  Penta::GradjDragStokes(Vec gradient){
     4628
     4629        int        i,j,ig;
     4630        int        drag_type,analysis_type;
     4631        int        doflist1[NUMVERTICES];
     4632        double     bed,thickness,Neff;
     4633        double     lambda,mu,xi,Jdet,vx,vy,vz;
     4634        double     alpha_complement,drag;
     4635        double     surface_normal[3],bed_normal[3];
     4636        double     xyz_list[NUMVERTICES][3];
     4637        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4638        double     dk[NDOF2];
     4639        double     basis[6];
     4640        double     grade_g[NUMVERTICES]={0.0};
     4641        double     grade_g_gaussian[NUMVERTICES];
     4642        Friction*  friction=NULL;
     4643        GaussPenta* gauss=NULL;
     4644
     4645        /*Gradient is 0 if on shelf or not on bed*/
     4646        if(IsOnShelf() || !IsOnBed()) return;
     4647
     4648        /*Retrieve all inputs and parameters*/
     4649        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4650        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4651        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4652        GetDofList1(&doflist1[0]);
     4653        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4654        Input* drag_input    =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input);
     4655        Input* vx_input      =inputs->GetInput(VxEnum);              _assert_(vx_input);
     4656        Input* vy_input      =inputs->GetInput(VyEnum);              _assert_(vy_input);
     4657        Input* vz_input      =inputs->GetInput(VzEnum);              _assert_(vz_input);
     4658        Input* adjointx_input=inputs->GetInput(AdjointxEnum);        _assert_(adjointx_input);
     4659        Input* adjointy_input=inputs->GetInput(AdjointyEnum);        _assert_(adjointy_input);
     4660        Input* adjointz_input=inputs->GetInput(AdjointzEnum);        _assert_(adjointz_input);
     4661
     4662        /*Build frictoin element, needed later: */
     4663        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4664        friction=new Friction("2d",inputs,matpar,analysis_type);
     4665
     4666        /* Start  looping on the number of gaussian points: */
     4667        gauss=new GaussPenta(0,1,2,4);
     4668        for(ig=gauss->begin();ig<gauss->end();ig++){
     4669
     4670                gauss->GaussPoint(ig);
     4671
     4672                /*Recover alpha_complement and drag: */
     4673                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
     4674                else alpha_complement=0;
     4675                drag_input->GetParameterValue(&drag,gauss);
     4676
     4677                /*recover lambda mu and xi: */
     4678                adjointx_input->GetParameterValue(&lambda,gauss);
     4679                adjointy_input->GetParameterValue(&mu    ,gauss);
     4680                adjointz_input->GetParameterValue(&xi    ,gauss);
     4681
     4682                /*recover vx vy and vz: */
     4683                vx_input->GetParameterValue(&vx, gauss);
     4684                vy_input->GetParameterValue(&vy, gauss);
     4685                vz_input->GetParameterValue(&vz, gauss);
     4686
     4687                /*Get normal vector to the bed */
     4688                SurfaceNormal(&surface_normal[0],xyz_list_tria);
     4689
     4690                bed_normal[0]=-surface_normal[0]; //Function is for upper surface, so the normal to the bed is the opposite of the result
     4691                bed_normal[1]=-surface_normal[1];
     4692                bed_normal[2]=-surface_normal[2];
     4693
     4694                /* Get Jacobian determinant: */
     4695                GetTriaJacobianDeterminant(&Jdet,&xyz_list_tria[0][0],gauss);
     4696                GetNodalFunctionsP1(basis, gauss);
     4697
     4698                /*Get k derivative: dk/dx */
     4699                drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
     4700
     4701                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     4702                for (i=0;i<NUMVERTICES;i++){
     4703                        //standard gradient dJ/dki
     4704                        grade_g_gaussian[i]=(
     4705                                                -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
     4706                                                -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
     4707                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
     4708                                                )*Jdet*gauss->weight*basis[i];
     4709                }
     4710
     4711                /*Add gradje_g_gaussian vector to gradje_g: */
     4712                for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i];
     4713        }
     4714
     4715        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
     4716
     4717        delete friction;
    43174718        delete gauss;
    4318         xfree((void**)&doflist);
    4319 }
    4320 /*}}}*/
    4321 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticVert{{{1*/
    4322 void  Penta::GetSolutionFromInputsDiagnosticVert(Vec solution){
    4323 
    4324         const int    numdof=NDOF1*NUMVERTICES;
    4325 
    4326         int          i;
    4327         int*         doflist=NULL;
    4328         double       vz;
    4329         double       values[numdof];
    4330         GaussPenta*  gauss=NULL;
    4331 
    4332         /*Get dof list: */
    4333         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4334         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    4335 
    4336         /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    4337         /*P1 element only for now*/
    4338         gauss=new GaussPenta();
    4339         for(i=0;i<NUMVERTICES;i++){
    4340                 /*Recover vz */
    4341                 gauss->GaussVertex(i);
    4342                 vz_input->GetParameterValue(&vz,gauss);
    4343                 values[i]=vz;
    4344         }
    4345 
    4346         /*Add value to global vector*/
    4347         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4348 
    4349         /*Free ressources:*/
    4350         delete gauss;
    4351         xfree((void**)&doflist);
    4352 }
    4353 /*}}}*/
    4354 /*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{1*/
    4355 void  Penta::GetSolutionFromInputsDiagnosticStokes(Vec solution){
    4356 
    4357         const int    numdof=NDOF4*NUMVERTICES;
    4358 
    4359         int          i;
    4360         int*         doflist=NULL;
    4361         double       vx,vy,vz,p;
    4362         double       stokesreconditioning;
    4363         double       values[numdof];
    4364         GaussPenta   *gauss;
    4365 
    4366         /*Get dof list: */
    4367         GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
    4368         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    4369         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    4370         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    4371         Input* p_input =inputs->GetInput(PressureEnum); _assert_(p_input);
    4372 
    4373         /*Recondition pressure: */
    4374         this->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
    4375 
    4376         /*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
    4377         /*P1 element only for now*/
    4378         gauss=new GaussPenta();
    4379         for(i=0;i<NUMVERTICES;i++){
    4380                 gauss->GaussVertex(i);
    4381                 vx_input->GetParameterValue(&vx,gauss);
    4382                 vy_input->GetParameterValue(&vy,gauss);
    4383                 vz_input->GetParameterValue(&vz,gauss);
    4384                 p_input ->GetParameterValue(&p ,gauss);
    4385                 values[i*NDOF4+0]=vx;
    4386                 values[i*NDOF4+1]=vy;
    4387                 values[i*NDOF4+2]=vz;
    4388                 values[i*NDOF4+3]=p/stokesreconditioning;
    4389         }
    4390 
    4391         /*Add value to global vector*/
    4392         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4393 
    4394         /*Free ressources:*/
    4395         delete gauss;
    4396         xfree((void**)&doflist);
    4397 }
    4398 /*}}}*/
    4399 /*FUNCTION Penta::GetSolutionFromInputsThermal{{{1*/
    4400 void  Penta::GetSolutionFromInputsThermal(Vec solution){
    4401 
    4402         const int    numdof=NDOF1*NUMVERTICES;
    4403 
    4404         int          i;
    4405         int*         doflist=NULL;
    4406         double       values[numdof];
    4407         double       temp;
    4408         GaussPenta   *gauss=NULL;
    4409 
    4410         /*Get dof list: */
    4411         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4412         Input* t_input=inputs->GetInput(TemperatureEnum); _assert_(t_input);
    4413 
    4414         gauss=new GaussPenta();
    4415         for(i=0;i<NUMVERTICES;i++){
    4416                 /*Recover temperature*/
    4417                 gauss->GaussVertex(i);
    4418                 t_input->GetParameterValue(&temp,gauss);
    4419                 values[i]=temp;
    4420         }
    4421 
    4422         /*Add value to global vector*/
    4423         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4424 
    4425         /*Free ressources:*/
    4426         delete gauss;
    4427         xfree((void**)&doflist);
    4428 }
    4429 /*}}}*/
    4430 /*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/
    4431 void  Penta::GetSolutionFromInputsEnthalpy(Vec solution){
    4432 
    4433         const int    numdof=NDOF1*NUMVERTICES;
    4434 
    4435         int          i;
    4436         int*         doflist=NULL;
    4437         double       values[numdof];
    4438         double       enthalpy;
    4439         GaussPenta   *gauss=NULL;
    4440 
    4441         /*Get dof list: */
    4442         GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
    4443         Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input);
    4444 
    4445         gauss=new GaussPenta();
    4446         for(i=0;i<NUMVERTICES;i++){
    4447                 /*Recover temperature*/
    4448                 gauss->GaussVertex(i);
    4449                 h_input->GetParameterValue(&enthalpy,gauss);
    4450                 values[i]=enthalpy;
    4451         }
    4452 
    4453         /*Add value to global vector*/
    4454         VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
    4455 
    4456         /*Free ressources:*/
    4457         delete gauss;
    4458         xfree((void**)&doflist);
    4459 }
    4460 /*}}}*/
    4461 /*FUNCTION Penta::GetStabilizationParameter {{{1*/
    4462 double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
    4463         /*Compute stabilization parameter*/
    4464 
    4465         double normu;
    4466         double tau_parameter;
    4467 
    4468         normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
    4469         if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
    4470                 tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
    4471         }
    4472         else tau_parameter=diameter/(2*normu);
    4473 
    4474         return tau_parameter;
    4475 }
    4476 /*}}}*/
    4477 /*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/
    4478 void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
    4479         /*Compute the 3d Blatter/PattynStrain Rate (5 components):
    4480          *
    4481          * epsilon=[exx eyy exy exz eyz]
    4482          *
    4483          * with exz=1/2 du/dz
    4484          *      eyz=1/2 dv/dz
    4485          *
    4486          * the contribution of vz is neglected
    4487          */
    4488 
    4489         int i;
    4490         double epsilonvx[5];
    4491         double epsilonvy[5];
    4492 
    4493         /*Check that both inputs have been found*/
    4494         if (!vx_input || !vy_input){
    4495                 _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
    4496         }
    4497 
    4498         /*Get strain rate assuming that epsilon has been allocated*/
    4499         vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
    4500         vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
    4501 
    4502         /*Sum all contributions*/
    4503         for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
    4504 }
    4505 /*}}}*/
    4506 /*FUNCTION Penta::GetStrainRate3d{{{1*/
    4507 void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
    4508         /*Compute the 3d Strain Rate (6 components):
    4509          *
    4510          * epsilon=[exx eyy ezz exy exz eyz]
    4511          */
    4512 
    4513         int i;
    4514         double epsilonvx[6];
    4515         double epsilonvy[6];
    4516         double epsilonvz[6];
    4517 
    4518         /*Check that both inputs have been found*/
    4519         if (!vx_input || !vy_input || !vz_input){
    4520                 _error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input);
    4521         }
    4522 
    4523         /*Get strain rate assuming that epsilon has been allocated*/
    4524         vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);
    4525         vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);
    4526         vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);
    4527 
    4528         /*Sum all contributions*/
    4529         for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];
    4530 }
    4531 /*}}}*/
    4532 /*FUNCTION Penta::GetUpperElement{{{1*/
    4533 Penta* Penta::GetUpperElement(void){
    4534 
    4535         Penta* upper_penta=NULL;
    4536 
    4537         upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above
    4538 
    4539         return upper_penta;
    4540 }
    4541 /*}}}*/
    4542 /*FUNCTION Penta::GetVectorFromInputs{{{1*/
    4543 void  Penta::GetVectorFromInputs(Vec vector,int input_enum){
    4544 
    4545         int doflist1[NUMVERTICES];
    4546 
    4547         /*Get out if this is not an element input*/
    4548         if (!IsInput(input_enum)) return;
    4549 
    4550         /*Prepare index list*/
    4551         this->GetDofList1(&doflist1[0]);
    4552 
    4553         /*Get input (either in element or material)*/
    4554         Input* input=inputs->GetInput(input_enum);
    4555         if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
    4556 
    4557         /*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
    4558         input->GetVectorFromInputs(vector,&doflist1[0]);
    4559 }
    4560 /*}}}*/
    4561 /*FUNCTION Penta::GetZcoord {{{1*/
    4562 double Penta::GetZcoord(GaussPenta* gauss){
    4563 
    4564         int    i;
    4565         double z;
    4566         double xyz_list[NUMVERTICES][3];
    4567         double z_list[NUMVERTICES];
    4568 
    4569         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4570         for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
    4571         PentaRef::GetParameterValue(&z,z_list,gauss);
    4572 
    4573         return z;
    4574 }
    4575 /*}}}*/
     4719}
     4720/*}}}*/
     4721/*FUNCTION Penta::GradjBbarMacAyeal {{{1*/
     4722void  Penta::GradjBbarMacAyeal(Vec gradient){
     4723
     4724        /*This element should be collapsed into a tria element at its base*/
     4725        if (!IsOnBed()) return;
     4726
     4727        /*Depth Average B*/
     4728        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     4729
     4730        /*Collapse element to the base*/
     4731        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria (lower face).
     4732        tria->GradjBMacAyeal(gradient);
     4733        delete tria->matice; delete tria;
     4734
     4735        /*delete Average B*/
     4736        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     4737
     4738} /*}}}*/
     4739/*FUNCTION Penta::GradjBbarPattyn {{{1*/
     4740void  Penta::GradjBbarPattyn(Vec gradient){
     4741
     4742        /*Gradient is computed on bed only (Bbar)*/
     4743        if (!IsOnBed()) return;
     4744
     4745        /*Depth Average B*/
     4746        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     4747
     4748        /*Collapse element to the base*/
     4749        Tria* tria=(Tria*)SpawnTria(0,1,2);
     4750        tria->GradjBMacAyeal(gradient);    //We use MacAyeal as an estimate for now
     4751        delete tria->matice; delete tria;
     4752
     4753        /*delete Average B*/
     4754        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     4755} /*}}}*/
     4756/*FUNCTION Penta::GradjBbarStokes {{{1*/
     4757void  Penta::GradjBbarStokes(Vec gradient){
     4758
     4759        /*Gradient is computed on bed only (Bbar)*/
     4760        if (!IsOnBed()) return;
     4761
     4762        /*Depth Average B*/
     4763        this->InputDepthAverageAtBase(RheologyBEnum,RheologyBbarEnum,MaterialsEnum);
     4764
     4765        /*Collapse element to the base*/
     4766        Tria* tria=(Tria*)SpawnTria(0,1,2);
     4767        tria->GradjBMacAyeal(gradient);    //We use MacAyeal as an estimate for now
     4768        delete tria->matice; delete tria;
     4769
     4770        /*delete Average B*/
     4771        this->matice->inputs->DeleteInput(RheologyBbarEnum);
     4772} /*}}}*/
    45764773/*FUNCTION Penta::Sid {{{1*/
    45774774int    Penta::Sid(){

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

    r8647 r8654  
    8888                void   GetVectorFromInputs(Vec vector,int NameEnum);
    8989                void   Gradj(Vec gradient,int control_type);
     90                void   GradjDragMacAyeal(Vec gradient);
     91                void   GradjDragPattyn(Vec gradient);
     92                void   GradjDragStokes(Vec gradient);
     93                void   GradjBbarMacAyeal(Vec gradient);
     94                void   GradjBbarPattyn(Vec gradient);
     95                void   GradjBbarStokes(Vec gradient);
    9096                int    Sid();
    9197                void   InputControlUpdate(double scalar,bool save_parameter);

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

    r8653 r8654  
    10441044        *Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
    10451045        if(*Jdet<0) _error_("negative jacobian determinant!");
    1046 
    10471046}
    10481047/*}}}*/

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

    r8653 r8654  
    29932993        double     xyz_list[NUMVERTICES][3];
    29942994        double     basis[3],epsilon[3];
    2995         double     dbasis[NDOF2][NUMVERTICES];
    29962995        double     grad[NUMVERTICES]={0.0};
    29972996        GaussTria *gauss = NULL;
    … …  
    30273026                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    30283027                GetNodalFunctions(basis,gauss);
    3029                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    30303028
    30313029                /*standard gradient dJ/dki*/
    … …  
    30503048        double     bed,thickness,Neff,drag;
    30513049        double     xyz_list[NUMVERTICES][3];
    3052         double     dbasis[NDOF2][NUMVERTICES];
    30533050        double     dk[NDOF2];
    30543051        double     grade_g[NUMVERTICES]={0.0};
    … …  
    30853082                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    30863083                GetNodalFunctions(basis, gauss);
    3087                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    30883084
    30893085                /*Build alpha_complement_list: */
    … …  
    31093105                }
    31103106        }
    3111 
    31123107        VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    31133108
    … …  
    31153110        delete gauss;
    31163111        delete friction;
    3117 }
    3118 /*}}}*/
    3119 /*FUNCTION Tria::GradjDragStokes {{{1*/
    3120 void  Tria::GradjDragStokes(Vec gradient){
    3121 
    3122         int        i,ig;
    3123         int        drag_type,analysis_type;
    3124         int        doflist1[NUMVERTICES];
    3125         double     bed,thickness,Neff;
    3126         double     lambda,mu,xi,Jdet,vx,vy,vz;
    3127         double     alpha_complement,drag;
    3128         double     surface_normal[3],bed_normal[3];
    3129         double     xyz_list[NUMVERTICES][3];
    3130         double     dbasis[NDOF2][NUMVERTICES];
    3131         double     dk[NDOF2];
    3132         double     basis[3];
    3133         double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    3134         double     grade_g[NUMVERTICES]={0.0};
    3135         double     grade_g_gaussian[NUMVERTICES];
    3136         Friction*  friction=NULL;
    3137         GaussTria* gauss=NULL;
    3138 
    3139         /*retrive parameters: */
    3140         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    3141 
    3142         /*retrieve inputs :*/
    3143         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    3144         Input* drag_input    =inputs->GetInput(DragCoefficientEnum); _assert_(drag_input);
    3145         Input* vx_input      =inputs->GetInput(VxEnum);              _assert_(vx_input);
    3146         Input* vy_input      =inputs->GetInput(VyEnum);              _assert_(vy_input);
    3147         Input* vz_input      =inputs->GetInput(VzEnum);              _assert_(vz_input);
    3148         Input* adjointx_input=inputs->GetInput(AdjointxEnum);        _assert_(adjointx_input);
    3149         Input* adjointy_input=inputs->GetInput(AdjointyEnum);        _assert_(adjointy_input);
    3150         Input* adjointz_input=inputs->GetInput(AdjointzEnum);        _assert_(adjointz_input);
    3151 
    3152         /*Get out if shelf*/
    3153         if(IsOnShelf())return;
    3154 
    3155         /* Get node coordinates and dof list: */
    3156         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3157         GetDofList1(&doflist1[0]);
    3158 
    3159         /*Build frictoin element, needed later: */
    3160         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    3161         friction=new Friction("2d",inputs,matpar,analysis_type);
    3162 
    3163         /* Start  looping on the number of gaussian points: */
    3164         gauss=new GaussTria(4);
    3165         for(ig=gauss->begin();ig<gauss->end();ig++){
    3166 
    3167                 gauss->GaussPoint(ig);
    3168 
    3169                 /*Recover alpha_complement and drag: */
    3170                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum);
    3171                 else alpha_complement=0;
    3172                 drag_input->GetParameterValue(&drag,gauss);
    3173 
    3174                 /*recover lambda mu and xi: */
    3175                 adjointx_input->GetParameterValue(&lambda,gauss);
    3176                 adjointy_input->GetParameterValue(&mu    ,gauss);
    3177                 adjointz_input->GetParameterValue(&xi    ,gauss);
    3178 
    3179                 /*recover vx vy and vz: */
    3180                 vx_input->GetParameterValue(&vx, gauss);
    3181                 vy_input->GetParameterValue(&vy, gauss);
    3182                 vz_input->GetParameterValue(&vz, gauss);
    3183 
    3184                 /*Get normal vector to the bed */
    3185                 SurfaceNormal(&surface_normal[0],xyz_list);
    3186 
    3187                 bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    3188                 bed_normal[1]=-surface_normal[1];
    3189                 bed_normal[2]=-surface_normal[2];
    3190 
    3191                 /* Get Jacobian determinant: */
    3192                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss);
    3193 
    3194                 /* Get nodal functions value at gaussian point:*/
    3195                 GetNodalFunctions(basis, gauss);
    3196 
    3197                 /*Get nodal functions derivatives*/
    3198                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    3199 
    3200                 /*Get k derivative: dk/dx */
    3201                 drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss);
    3202 
    3203                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    3204                 for (i=0;i<NUMVERTICES;i++){
    3205                         //standard gradient dJ/dki
    3206                         grade_g_gaussian[i]=(
    3207                                                 -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
    3208                                                 -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
    3209                                                 -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    3210                                                 )*Jdet*gauss->weight*basis[i];
    3211                 }
    3212 
    3213                 /*Add gradje_g_gaussian vector to gradje_g: */
    3214                 for( i=0; i<NUMVERTICES; i++)grade_g[i]+=grade_g_gaussian[i];
    3215         }
    3216 
    3217         VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
    3218 
    3219         delete friction;
    3220         delete gauss;
    32213112}
    32223113/*}}}*/
    … …  
    32903181        double     thickness,Jdet;
    32913182        double     basis[3];
    3292         double     dbasis[NDOF2][NUMVERTICES];
    32933183        double     Dlambda[2],dp[2];
    32943184        double     xyz_list[NUMVERTICES][3];
    … …  
    33123202                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    33133203                GetNodalFunctions(basis, gauss);
    3314                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    33153204               
    33163205                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    … …  
    33353224        double     thickness,Jdet;
    33363225        double     basis[3];
    3337         double     dbasis[NDOF2][NUMVERTICES];
    33383226        double     Dlambda[2],dp[2];
    33393227        double     xyz_list[NUMVERTICES][3];
    … …  
    33573245                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    33583246                GetNodalFunctions(basis, gauss);
    3359                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    33603247
    33613248                adjoint_input->GetParameterDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
    … …  
    53715258        double     xyz_list[NUMVERTICES][3];
    53725259        GaussTria *gauss = NULL;
    5373         double     dbasis[NDOF2][NUMVERTICES];
    53745260        double     dH[2];
    53755261
    … …  
    53915277                /* Get Jacobian determinant: */
    53925278                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
    5393                 GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    53945279
    53955280                /*Get parameters at gauss point*/

  • issm/trunk/src/c/objects/Loads/Friction.cpp

    r8224 r8654  
    185185}
    186186/*}}}*/
    187 /*FUNCTION Friction::GetAlphaComplement {{{1*/
     187/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum) {{{1*/
    188188void Friction::GetAlphaComplement(double* palpha_complement, GaussTria* gauss,int vxenum,int vyenum){
    189189
    … …  
    239239        if(vmag==0 && (s-1)<0) _error_("velocity is 0 ans (s-1)=%g<0, alpha_complement is Inf",s-1);
    240240
    241         alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    242         _assert_(!isnan(alpha_complement));
     241        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!isnan(alpha_complement));
    243242
    244243        /*Assign output pointers:*/
    245244        *palpha_complement=alpha_complement;
    246 
     245}
     246/*}}}*/
     247/*FUNCTION Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum) {{{1*/
     248void Friction::GetAlphaComplement(double* palpha_complement, GaussPenta* gauss,int vxenum,int vyenum){
     249
     250        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
     251         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
     252         * alpha_complement= Neff ^r * vel ^s*/
     253
     254        /*diverse: */
     255        int     i;
     256        double  Neff;
     257        double  r,s;
     258        double  vx;
     259        double  vy;
     260        double  vmag;
     261        double  drag_p,drag_q;
     262        double  drag_coefficient;
     263        double  bed,thickness;
     264        double  gravity,rho_ice,rho_water;
     265        double  alpha_complement;
     266
     267        /*Recover parameters: */
     268        inputs->GetParameterValue(&drag_p,DragPEnum);
     269        inputs->GetParameterValue(&drag_q,DragQEnum);
     270        this->GetParameterValue(&thickness, gauss,ThicknessEnum);
     271        this->GetParameterValue(&bed, gauss,BedEnum);
     272        this->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     273
     274        /*Get material parameters: */
     275        gravity=matpar->GetG();
     276        rho_ice=matpar->GetRhoIce();
     277        rho_water=matpar->GetRhoWater();
     278
     279
     280        //compute r and q coefficients: */
     281        r=drag_q/drag_p;
     282        s=1./drag_p;
     283
     284        //From bed and thickness, compute effective pressure when drag is viscous:
     285        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     286
     287        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     288          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     289          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     290          replace it by Neff=0 (ie, equival it to an ice shelf)*/
     291        if (Neff<0)Neff=0;
     292
     293        //We need the velocity magnitude to evaluate the basal stress:
     294        this->GetParameterValue(&vx, gauss,vxenum);
     295        this->GetParameterValue(&vy, gauss,vyenum);
     296        vmag=sqrt(pow(vx,2)+pow(vy,2));
     297
     298        /*Checks that s-1>0 if v=0*/
     299        if(vmag==0 && (s-1)<0) _error_("velocity is 0 ans (s-1)=%g<0, alpha_complement is Inf",s-1);
     300
     301        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));            _assert_(!isnan(alpha_complement));
     302
     303        /*Assign output pointers:*/
     304        *palpha_complement=alpha_complement;
    247305}
    248306/*}}}*/

  • issm/trunk/src/c/objects/Loads/Friction.h

    r5749 r8654  
    3030                void  GetAlpha2(double* palpha2, GaussPenta* gauss,int vxenum,int vyenum,int vzenum);
    3131                void  GetAlphaComplement(double* alpha_complement, GaussTria* gauss,int vxenum,int vyenum);
     32                void  GetAlphaComplement(double* alpha_complement, GaussPenta* gauss,int vxenum,int vyenum);
    3233                void  GetParameterValue(double* pvalue,GaussTria* gauss,int enum_type);
    3334                void  GetParameterValue(double* pvalue,GaussPenta* gauss,int enum_type);
Note: See TracChangeset for help on using the changeset viewer.