Ignore:
Timestamp:
06/15/10 11:13:28 (15 years ago)
Author:
Eric.Larour
Message:

Starting cleanup of control_core.
Introduced gradient_core and adjoint_core solutions. Trying to get the gradient and the adjoint
inside the elements, instead of as node vectors.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4046 r4047  
    41254125/*}}}*/
    41264126/*FUNCTION Tria::Gradj {{{1*/
    4127 void  Tria::Gradj(Vec grad_g,int control_type){
     4127void  Tria::Gradj(int control_type){
    41284128
    41294129        /*inputs: */
    … …  
    41374137
    41384138        if (control_type==DragCoefficientEnum){
    4139                 GradjDrag( grad_g);
     4139                GradjDrag();
    41404140        }
    41414141        else if (control_type==RheologyBEnum){
    4142                 GradjB( grad_g);
     4142                GradjB();
    41434143        }
    41444144        else ISSMERROR("%s%i","control type not supported yet: ",control_type);
     4145}
     4146/*}}}*/
     4147/*FUNCTION Tria::GradjDrag {{{1*/
     4148void  Tria::GradjDrag(void){
     4149
     4150
     4151        int i;
     4152
     4153        /* node data: */
     4154        const int    numgrids=3;
     4155        const int    NDOF2=2;
     4156        const int    numdof=NDOF2*numgrids;
     4157        double       xyz_list[numgrids][3];
     4158        int          doflist1[numgrids];
     4159        double       dh1dh3[NDOF2][numgrids];
     4160
     4161        /* grid data: */
     4162        double adjx_list[numgrids];
     4163        double adjy_list[numgrids];
     4164
     4165        /* gaussian points: */
     4166        int     num_gauss,ig;
     4167        double* first_gauss_area_coord  =  NULL;
     4168        double* second_gauss_area_coord =  NULL;
     4169        double* third_gauss_area_coord  =  NULL;
     4170        double* gauss_weights           =  NULL;
     4171        double  gauss_weight;
     4172        double  gauss_l1l2l3[3];
     4173
     4174        /* parameters: */
     4175        double  dk[NDOF2];
     4176        double  vx,vy;
     4177        double  lambda,mu;
     4178        double  bed,thickness,Neff;
     4179        double  alpha_complement;
     4180        int     drag_type;
     4181        double  drag;
     4182        Friction* friction=NULL;
     4183
     4184        /*element vector at the gaussian points: */
     4185        double  grade_g[numgrids]={0.0};
     4186        double  grade_g_gaussian[numgrids];
     4187
     4188        /* Jacobian: */
     4189        double Jdet;
     4190
     4191        /*nodal functions: */
     4192        double l1l2l3[3];
     4193
     4194        /* strain rate: */
     4195        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     4196
     4197        /*inputs: */
     4198        bool shelf;
     4199
     4200        /*parameters: */
     4201        double  cm_noisedmp;
     4202        double  cm_mindmp_slope;
     4203        double  cm_mindmp_value;
     4204        double  cm_maxdmp_value;
     4205        double  cm_maxdmp_slope;
     4206
     4207        int analysis_type;
     4208
     4209        /*retrive parameters: */
     4210        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4211
     4212        /*retrieve inputs :*/
     4213        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     4214
     4215        /*retrieve some parameters: */
     4216        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
     4217        this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
     4218        this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
     4219        this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
     4220        this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
     4221
     4222
     4223        /*Get out if shelf*/
     4224        if(shelf)return;
     4225
     4226        /* Get node coordinates and dof list: */
     4227        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4228        GetDofList1(&doflist1[0]);
     4229
     4230        /*Build frictoin element, needed later: */
     4231        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4232        friction=new Friction("2d",inputs,matpar,analysis_type);
     4233
     4234        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     4235        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     4236
     4237        /* Start  looping on the number of gaussian points: */
     4238        for (ig=0; ig<num_gauss; ig++){
     4239                /*Pick up the gaussian point: */
     4240                gauss_weight=*(gauss_weights+ig);
     4241                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     4242                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     4243                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4244
     4245                /*Build alpha_complement_list: */
     4246                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
     4247                else alpha_complement=0;
     4248       
     4249                /*Recover alpha_complement and k: */
     4250                inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
     4251
     4252                /*recover lambda and mu: */
     4253                inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
     4254                inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
     4255                       
     4256                /*recover vx and vy: */
     4257                inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
     4258                inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
     4259
     4260                /* Get Jacobian determinant: */
     4261                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4262               
     4263                /* Get nodal functions value at gaussian point:*/
     4264                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4265
     4266                /*Get nodal functions derivatives*/
     4267                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
     4268
     4269                /*Get k derivative: dk/dx */
     4270                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
     4271
     4272                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     4273                for (i=0;i<numgrids;i++){
     4274
     4275                        //standard term dJ/dki
     4276                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
     4277
     4278                        //noise dampening d/dki(1/2*(dk/dx)^2)
     4279                        grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     4280                       
     4281                        //min dampening
     4282                        if(drag<cm_mindmp_value){
     4283                                grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4284                        }
     4285
     4286                        //max dampening
     4287                        if(drag>cm_maxdmp_value){
     4288                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4289                        }
     4290                }
     4291               
     4292                /*Add gradje_g_gaussian vector to gradje_g: */
     4293                for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
     4294        }
     4295
     4296
     4297        /*Add grade_g to the inputs of this element: */
     4298        this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
     4299
     4300        cleanup_and_return:
     4301        xfree((void**)&first_gauss_area_coord);
     4302        xfree((void**)&second_gauss_area_coord);
     4303        xfree((void**)&third_gauss_area_coord);
     4304        xfree((void**)&gauss_weights);
     4305        delete friction;
     4306
     4307}
     4308/*}}}*/
     4309/*FUNCTION Tria::GradjDragStokes {{{1*/
     4310void  Tria::GradjDragStokes(Vec grad_g){
     4311
     4312        int i;
     4313
     4314        /* node data: */
     4315        const int    numgrids=3;
     4316        const int    NDOF2=2;
     4317        double       xyz_list[numgrids][3];
     4318        int          doflist1[numgrids];
     4319        double       dh1dh3[NDOF2][numgrids];
     4320
     4321        /* grid data: */
     4322        double drag;
     4323        double alpha_complement;
     4324        Friction* friction=NULL;
     4325
     4326        /* gaussian points: */
     4327        int     num_gauss,ig;
     4328        double* first_gauss_area_coord  =  NULL;
     4329        double* second_gauss_area_coord =  NULL;
     4330        double* third_gauss_area_coord  =  NULL;
     4331        double* gauss_weights           =  NULL;
     4332        double  gauss_weight;
     4333        double  gauss_l1l2l3[3];
     4334        double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
     4335
     4336        /* parameters: */
     4337        double  vx,vy,vz;
     4338        double  lambda,mu,xi;
     4339        double  bed,thickness,Neff;
     4340        double  surface_normal[3];
     4341        double  bed_normal[3];
     4342        double  dk[NDOF2];
     4343
     4344        /*element vector at the gaussian points: */
     4345        double  grade_g[numgrids]={0.0};
     4346        double  grade_g_gaussian[numgrids];
     4347
     4348        /* Jacobian: */
     4349        double Jdet;
     4350
     4351        /*nodal functions: */
     4352        double l1l2l3[3];
     4353
     4354        /* strain rate: */
     4355        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     4356
     4357        /*inputs: */
     4358        bool shelf;
     4359        int  drag_type;
     4360
     4361        /*parameters: */
     4362        double  cm_noisedmp;
     4363        double  cm_mindmp_slope;
     4364        double  cm_mindmp_value;
     4365        double  cm_maxdmp_value;
     4366        double  cm_maxdmp_slope;
     4367
     4368        int analysis_type;
     4369
     4370        /*retrive parameters: */
     4371        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     4372
     4373        /*retrieve inputs :*/
     4374        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     4375        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4376
     4377        /*retrieve some parameters: */
     4378        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
     4379        this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
     4380        this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
     4381        this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
     4382        this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
     4383
     4384        /*Get out if shelf*/
     4385        if(shelf)return;
     4386
     4387        /* Get node coordinates and dof list: */
     4388        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4389        GetDofList1(&doflist1[0]);
     4390
     4391        /*Build frictoin element, needed later: */
     4392        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4393        friction=new Friction("2d",inputs,matpar,analysis_type);
     4394
     4395        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     4396        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     4397
     4398        /* Start  looping on the number of gaussian points: */
     4399        for (ig=0; ig<num_gauss; ig++){
     4400                /*Pick up the gaussian point: */
     4401                gauss_weight=*(gauss_weights+ig);
     4402                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     4403                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     4404                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4405
     4406                /*Recover alpha_complement and drag: */
     4407                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
     4408                else alpha_complement=0;
     4409                inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
     4410
     4411                /*recover lambda mu and xi: */
     4412                inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum);
     4413                inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum);
     4414                inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum);
     4415
     4416                /*recover vx vy and vz: */
     4417                inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum);
     4418                inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum);
     4419                inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum);
     4420
     4421                /*Get normal vector to the bed */
     4422                SurfaceNormal(&surface_normal[0],xyz_list);
     4423
     4424                bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
     4425                bed_normal[1]=-surface_normal[1];
     4426                bed_normal[2]=-surface_normal[2];
     4427
     4428                /* Get Jacobian determinant: */
     4429                GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4430
     4431                /* Get nodal functions value at gaussian point:*/
     4432                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4433
     4434                /*Get nodal functions derivatives*/
     4435                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
     4436
     4437                /*Get k derivative: dk/dx */
     4438                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
     4439
     4440                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     4441                for (i=0;i<numgrids;i++){
     4442                        //standard gradient dJ/dki
     4443                        grade_g_gaussian[i]=(
     4444                                                -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
     4445                                                -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
     4446                                                -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
     4447                                                )*Jdet*gauss_weight*l1l2l3[i];
     4448
     4449                        //Add regularization term
     4450                        grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     4451
     4452                        //min dampening
     4453                        if(drag<cm_mindmp_value){
     4454                                grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4455                        }
     4456
     4457                        //max dampening
     4458                        if(drag>cm_maxdmp_value){
     4459                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4460                        }
     4461                }
     4462
     4463                /*Add gradje_g_gaussian vector to gradje_g: */
     4464                for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
     4465        }
     4466
     4467        /*Add grade_g to the inputs of this element: */
     4468        this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
     4469
     4470        cleanup_and_return:
     4471        xfree((void**)&first_gauss_area_coord);
     4472        xfree((void**)&second_gauss_area_coord);
     4473        xfree((void**)&third_gauss_area_coord);
     4474        xfree((void**)&gauss_weights);
     4475        delete friction;
     4476
    41454477}
    41464478/*}}}*/
    … …  
    42814613        }
    42824614
    4283         /*Add grade_g to global vector grad_g: */
    4284         VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    4285 
    4286 cleanup_and_return:
    4287         xfree((void**)&first_gauss_area_coord);
    4288         xfree((void**)&second_gauss_area_coord);
    4289         xfree((void**)&third_gauss_area_coord);
    4290         xfree((void**)&gauss_weights);
    4291 }
    4292 /*}}}*/
    4293 /*FUNCTION Tria::GradjDrag {{{1*/
    4294 void  Tria::GradjDrag(Vec grad_g){
    4295 
    4296 
    4297         int i;
    4298 
    4299         /* node data: */
    4300         const int    numgrids=3;
    4301         const int    NDOF2=2;
    4302         const int    numdof=NDOF2*numgrids;
    4303         double       xyz_list[numgrids][3];
    4304         int          doflist1[numgrids];
    4305         double       dh1dh3[NDOF2][numgrids];
    4306 
    4307         /* grid data: */
    4308         double adjx_list[numgrids];
    4309         double adjy_list[numgrids];
    4310 
    4311         /* gaussian points: */
    4312         int     num_gauss,ig;
    4313         double* first_gauss_area_coord  =  NULL;
    4314         double* second_gauss_area_coord =  NULL;
    4315         double* third_gauss_area_coord  =  NULL;
    4316         double* gauss_weights           =  NULL;
    4317         double  gauss_weight;
    4318         double  gauss_l1l2l3[3];
    4319 
    4320         /* parameters: */
    4321         double  dk[NDOF2];
    4322         double  vx,vy;
    4323         double  lambda,mu;
    4324         double  bed,thickness,Neff;
    4325         double  alpha_complement;
    4326         int     drag_type;
    4327         double  drag;
    4328         Friction* friction=NULL;
    4329 
    4330         /*element vector at the gaussian points: */
    4331         double  grade_g[numgrids]={0.0};
    4332         double  grade_g_gaussian[numgrids];
    4333 
    4334         /* Jacobian: */
    4335         double Jdet;
    4336 
    4337         /*nodal functions: */
    4338         double l1l2l3[3];
    4339 
    4340         /* strain rate: */
    4341         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    4342 
    4343         /*inputs: */
    4344         bool shelf;
    4345 
    4346         /*parameters: */
    4347         double  cm_noisedmp;
    4348         double  cm_mindmp_slope;
    4349         double  cm_mindmp_value;
    4350         double  cm_maxdmp_value;
    4351         double  cm_maxdmp_slope;
    4352 
    4353         int analysis_type;
    4354 
    4355         /*retrive parameters: */
    4356         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4357 
    4358         /*retrieve inputs :*/
    4359         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    4360 
    4361         /*retrieve some parameters: */
    4362         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    4363         this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
    4364         this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
    4365         this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
    4366         this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
    4367 
    4368 
    4369         /*Get out if shelf*/
    4370         if(shelf)return;
    4371 
    4372         /* Get node coordinates and dof list: */
    4373         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    4374         GetDofList1(&doflist1[0]);
    4375 
    4376         /*Build frictoin element, needed later: */
    4377         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4378         friction=new Friction("2d",inputs,matpar,analysis_type);
    4379 
    4380         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4381         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    4382 
    4383         /* Start  looping on the number of gaussian points: */
    4384         for (ig=0; ig<num_gauss; ig++){
    4385                 /*Pick up the gaussian point: */
    4386                 gauss_weight=*(gauss_weights+ig);
    4387                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4388                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4389                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    4390 
    4391                 /*Build alpha_complement_list: */
    4392                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
    4393                 else alpha_complement=0;
    4394        
    4395                 /*Recover alpha_complement and k: */
    4396                 inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
    4397 
    4398                 /*recover lambda and mu: */
    4399                 inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
    4400                 inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
    4401                        
    4402                 /*recover vx and vy: */
    4403                 inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
    4404                 inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
    4405 
    4406                 /* Get Jacobian determinant: */
    4407                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    4408                
    4409                 /* Get nodal functions value at gaussian point:*/
    4410                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    4411 
    4412                 /*Get nodal functions derivatives*/
    4413                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    4414 
    4415                 /*Get k derivative: dk/dx */
    4416                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    4417 
    4418                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    4419                 for (i=0;i<numgrids;i++){
    4420 
    4421                         //standard term dJ/dki
    4422                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
    4423 
    4424                         //noise dampening d/dki(1/2*(dk/dx)^2)
    4425                         grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    4426                        
    4427                         //min dampening
    4428                         if(drag<cm_mindmp_value){
    4429                                 grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4430                         }
    4431 
    4432                         //max dampening
    4433                         if(drag>cm_maxdmp_value){
    4434                                 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4435                         }
    4436                 }
    4437                
    4438                 /*Add gradje_g_gaussian vector to gradje_g: */
    4439                 for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
    4440         }
    4441 
    4442         /*Add grade_g to global vector grad_g: */
    4443         VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    4444        
     4615        /*Add grade_g to the inputs of this element: */
     4616        this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
     4617
    44454618        cleanup_and_return:
    44464619        xfree((void**)&first_gauss_area_coord);
    … …  
    44484621        xfree((void**)&third_gauss_area_coord);
    44494622        xfree((void**)&gauss_weights);
    4450         delete friction;
    4451 
    4452 }
    4453 /*}}}*/
    4454 /*FUNCTION Tria::GradjDragStokes {{{1*/
    4455 void  Tria::GradjDragStokes(Vec grad_g){
    4456 
    4457         int i;
    4458 
    4459         /* node data: */
    4460         const int    numgrids=3;
    4461         const int    NDOF2=2;
    4462         double       xyz_list[numgrids][3];
    4463         int          doflist1[numgrids];
    4464         double       dh1dh3[NDOF2][numgrids];
    4465 
    4466         /* grid data: */
    4467         double drag;
    4468         double alpha_complement;
    4469         Friction* friction=NULL;
    4470 
    4471         /* gaussian points: */
    4472         int     num_gauss,ig;
    4473         double* first_gauss_area_coord  =  NULL;
    4474         double* second_gauss_area_coord =  NULL;
    4475         double* third_gauss_area_coord  =  NULL;
    4476         double* gauss_weights           =  NULL;
    4477         double  gauss_weight;
    4478         double  gauss_l1l2l3[3];
    4479         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    4480 
    4481         /* parameters: */
    4482         double  vx,vy,vz;
    4483         double  lambda,mu,xi;
    4484         double  bed,thickness,Neff;
    4485         double  surface_normal[3];
    4486         double  bed_normal[3];
    4487         double  dk[NDOF2];
    4488 
    4489         /*element vector at the gaussian points: */
    4490         double  grade_g[numgrids]={0.0};
    4491         double  grade_g_gaussian[numgrids];
    4492 
    4493         /* Jacobian: */
    4494         double Jdet;
    4495 
    4496         /*nodal functions: */
    4497         double l1l2l3[3];
    4498 
    4499         /* strain rate: */
    4500         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    4501 
    4502         /*inputs: */
    4503         bool shelf;
    4504         int  drag_type;
    4505 
    4506         /*parameters: */
    4507         double  cm_noisedmp;
    4508         double  cm_mindmp_slope;
    4509         double  cm_mindmp_value;
    4510         double  cm_maxdmp_value;
    4511         double  cm_maxdmp_slope;
    4512 
    4513         int analysis_type;
    4514 
    4515         /*retrive parameters: */
    4516         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4517 
    4518         /*retrieve inputs :*/
    4519         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    4520         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4521 
    4522         /*retrieve some parameters: */
    4523         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    4524         this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
    4525         this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
    4526         this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
    4527         this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
    4528 
    4529         /*Get out if shelf*/
    4530         if(shelf)return;
    4531 
    4532         /* Get node coordinates and dof list: */
    4533         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    4534         GetDofList1(&doflist1[0]);
    4535 
    4536         /*Build frictoin element, needed later: */
    4537         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4538         friction=new Friction("2d",inputs,matpar,analysis_type);
    4539 
    4540         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4541         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    4542 
    4543         /* Start  looping on the number of gaussian points: */
    4544         for (ig=0; ig<num_gauss; ig++){
    4545                 /*Pick up the gaussian point: */
    4546                 gauss_weight=*(gauss_weights+ig);
    4547                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4548                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4549                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    4550 
    4551                 /*Recover alpha_complement and drag: */
    4552                 if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
    4553                 else alpha_complement=0;
    4554                 inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
    4555 
    4556                 /*recover lambda mu and xi: */
    4557                 inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum);
    4558                 inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum);
    4559                 inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum);
    4560 
    4561                 /*recover vx vy and vz: */
    4562                 inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum);
    4563                 inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum);
    4564                 inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum);
    4565 
    4566                 /*Get normal vector to the bed */
    4567                 SurfaceNormal(&surface_normal[0],xyz_list);
    4568 
    4569                 bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    4570                 bed_normal[1]=-surface_normal[1];
    4571                 bed_normal[2]=-surface_normal[2];
    4572 
    4573                 /* Get Jacobian determinant: */
    4574                 GetJacobianDeterminant3d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    4575 
    4576                 /* Get nodal functions value at gaussian point:*/
    4577                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    4578 
    4579                 /*Get nodal functions derivatives*/
    4580                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    4581 
    4582                 /*Get k derivative: dk/dx */
    4583                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    4584 
    4585                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    4586                 for (i=0;i<numgrids;i++){
    4587                         //standard gradient dJ/dki
    4588                         grade_g_gaussian[i]=(
    4589                                                 -lambda*(2*drag*alpha_complement*(vx - vz*bed_normal[0]*bed_normal[2]))
    4590                                                 -mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
    4591                                                 -xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
    4592                                                 )*Jdet*gauss_weight*l1l2l3[i];
    4593 
    4594                         //Add regularization term
    4595                         grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    4596 
    4597                         //min dampening
    4598                         if(drag<cm_mindmp_value){
    4599                                 grade_g_gaussian[i]+= cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4600                         }
    4601 
    4602                         //max dampening
    4603                         if(drag>cm_maxdmp_value){
    4604                                 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4605                         }
    4606                 }
    4607 
    4608                 /*Add gradje_g_gaussian vector to gradje_g: */
    4609                 for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
    4610         }
    4611 
    4612         /*Add grade_g to global vector grad_g: */
    4613         VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    4614 
    4615         cleanup_and_return:
    4616         xfree((void**)&first_gauss_area_coord);
    4617         xfree((void**)&second_gauss_area_coord);
    4618         xfree((void**)&third_gauss_area_coord);
    4619         xfree((void**)&gauss_weights);
    4620         delete friction;
    4621 
    46224623}
    46234624/*}}}*/
    … …  
    54375438}
    54385439/*}}}*/
     5440/*FUNCTION Tria::ScaleInput(int enum_type,double scale_factor){{{1*/
     5441void  Tria::ScaleInput(int enum_type,double scale_factor){
     5442
     5443        Input* input=NULL;
     5444
     5445        /*Make a copy of the original input: */
     5446        input=(Input*)this->inputs->GetInput(enum_type);
     5447
     5448        /*Scale: */
     5449        input->Scale(scale_factor);
     5450}
     5451/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.