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

moved some Penta functions from Tria to Penta

File:
1 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(){
Note: See TracChangeset for help on using the changeset viewer.