Changeset 3833


Ignore:
Timestamp:
05/18/10 21:22:39 (15 years ago)
Author:
Eric.Larour
Message:

Moved Friction2, now that it compiles, to Friction

Location:
issm/trunk/src/c
Files:
2 deleted
6 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r3832 r3833  
    5050                                        ./objects/Loads/Friction.h\
    5151                                        ./objects/Loads/Friction.cpp\
    52                                         ./objects/Loads/Friction2.h\
    53                                         ./objects/Loads/Friction2.cpp\
    5452                                        ./objects/DakotaPlugin.h\
    5553                                        ./objects/DakotaPlugin.cpp\
    … …  
    499497                                        ./objects/Loads/Friction.h\
    500498                                        ./objects/Loads/Friction.cpp\
    501                                         ./objects/Loads/Friction2.h\
    502                                         ./objects/Loads/Friction2.cpp\
    503499                                        ./objects/DakotaPlugin.h\
    504500                                        ./objects/DakotaPlugin.cpp\

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

    r3824 r3833  
    11801180        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    11811181        double  viscosity;
    1182         double  alpha2_list[numgrids2d];
    11831182        double  alpha2_gauss;
    1184 
    1185         double  vx_list[numgrids];
    1186         double  vy_list[numgrids];
    1187         double  vz_list[numgrids];
    1188         double  thickness_list[numgrids];
    1189         double  bed_list[numgrids];
    1190         double  dragcoefficient_list[numgrids];
    1191         double  drag_p,drag_q;
     1183        Friction* friction=NULL;
    11921184
    11931185        /*parameters: */
    … …  
    12811273        if((onbed==1) && (shelf==0)){
    12821274
    1283                 /*Build alpha2_list used by drag stiffness matrix*/
    1284                 Friction* friction=NewFriction();
    1285 
    1286                 /*Initialize all fields: */
    1287                 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    1288                 strcpy(friction->element_type,"2d");
    1289 
    1290                 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum);
    1291                 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum);
    1292                 inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum);
    1293                 inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum);
    1294                 inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum);
    1295                 inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum);
    1296                 inputs->GetParameterValue(&drag_p,DragPEnum);
    1297                 inputs->GetParameterValue(&drag_q,DragQEnum);
    1298 
    1299                 friction->gravity=matpar->GetG();
    1300                 friction->rho_ice=matpar->GetRhoIce();
    1301                 friction->rho_water=matpar->GetRhoWater();
    1302                 friction->K=&dragcoefficient_list[0];
    1303                 friction->bed=&bed_list[0];
    1304                 friction->thickness=&thickness_list[0];
    1305                 friction->vx=&vx_list[0];
    1306                 friction->vy=&vy_list[0];
    1307                 friction->vz=&vz_list[0];
    1308                 friction->p=drag_p;
    1309                 friction->q=drag_q;
    1310 
    1311                 /*Compute alpha2_list: */
    1312                 FrictionGetAlpha2(&alpha2_list[0],friction);
    1313 
    1314                 /*Erase friction object: */
    1315                 DeleteFriction(&friction);
     1275                /*build friction object, used later on: */
     1276                friction=new Friction("2d",inputs,matpar);
    13161277
    13171278                for(i=0;i<numgrids2d;i++){
    … …  
    13571318
    13581319                        /*Calculate DL on gauss point */
    1359                         tria->GetParameterValue(&alpha2_gauss,&alpha2_list[0],gauss_coord_tria);
     1320                        friction->GetAlpha2(&alpha2_gauss, gauss_coord_tria,VxEnum,VyEnum,VzEnum);
    13601321
    13611322                        DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d;
    … …  
    13841345                        }
    13851346                }
     1347       
     1348                /*Free ressources:*/
     1349                delete friction;
     1350
    13861351        } //if ( (onbed==1) && (shelf==0))
    13871352

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

    r3832 r3833  
    12991299
    13001300        /*friction: */
    1301         Friction2* friction2=NULL;
     1301        Friction* friction=NULL;
    13021302        double alpha2;
    13031303
    … …  
    13241324        /*build friction object, used later on: */
    13251325        if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
    1326         friction2=new Friction2("2d",inputs,matpar);
     1326        friction=new Friction("2d",inputs,matpar);
    13271327
    13281328        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    13391339
    13401340                /*Friction: */
    1341                 friction2->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum);
     1341                friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum);
    13421342
    13431343                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    … …  
    13801380        xfree((void**)&third_gauss_area_coord);
    13811381        xfree((void**)&gauss_weights);
    1382         delete friction2;
     1382        delete friction;
    13831383
    13841384}       
    … …  
    29022902        int    drag_type;
    29032903        double basalfriction;
    2904         Friction2* friction2=NULL;
     2904        Friction* friction=NULL;
    29052905        double alpha2,vx,vy;
    29062906        double geothermalflux_value;
    … …  
    29362936        inputs->GetParameterValue(&drag_type,DragTypeEnum);
    29372937        if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
    2938         friction2=new Friction2("3d",inputs,matpar);
     2938        friction=new Friction("3d",inputs,matpar);
    29392939       
    29402940        /* Ice/ocean heat exchange flux on ice shelf base */
    … …  
    29572957                inputs->GetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum);
    29582958       
    2959                 friction2->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum);
     2959                friction->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum);
    29602960                inputs->GetParameterValue(&vx, &gauss_coord[0],VxAverageEnum);
    29612961                inputs->GetParameterValue(&vy, &gauss_coord[0],VyAverageEnum);
    … …  
    29812981        xfree((void**)&third_gauss_area_coord);
    29822982        xfree((void**)&gauss_weights);
    2983         delete friction2;
     2983        delete friction;
    29842984
    29852985}
    … …  
    39433943
    39443944        /* grid data: */
    3945         double vx_list[numgrids];
    3946         double vy_list[numgrids];
    39473945        double adjx_list[numgrids];
    39483946        double adjy_list[numgrids];
    3949         double thickness_list[numgrids];
    3950         double bed_list[numgrids];
    3951         double dragcoefficient_list[numgrids];
    3952         double drag_p;
    3953         double drag_q;
    3954         int    drag_type;
    3955 
     3947
     3948        /* gaussian points: */
     3949        int     num_gauss,ig;
     3950        double* first_gauss_area_coord  =  NULL;
     3951        double* second_gauss_area_coord =  NULL;
     3952        double* third_gauss_area_coord  =  NULL;
     3953        double* gauss_weights           =  NULL;
     3954        double  gauss_weight;
     3955        double  gauss_l1l2l3[3];
     3956
     3957        /* parameters: */
     3958        double  dk[NDOF2];
     3959        double  vx,vy;
     3960        double  lambda,mu;
     3961        double  bed,thickness,Neff;
     3962        double  alpha_complement;
     3963        int     drag_type;
     3964        double  drag;
     3965        Friction* friction=NULL;
     3966
     3967        /*element vector at the gaussian points: */
     3968        double  grade_g[numgrids]={0.0};
     3969        double  grade_g_gaussian[numgrids];
     3970
     3971        /* Jacobian: */
     3972        double Jdet;
     3973
     3974        /*nodal functions: */
     3975        double l1l2l3[3];
     3976
     3977        /* strain rate: */
     3978        double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
     3979
     3980        /*inputs: */
     3981        bool shelf;
     3982
     3983        /*parameters: */
     3984        double  cm_noisedmp;
     3985        double  cm_mindmp_slope;
     3986        double  cm_mindmp_value;
     3987        double  cm_maxdmp_value;
     3988        double  cm_maxdmp_slope;
     3989
     3990        /*retrieve inputs :*/
     3991        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     3992
     3993        /*retrieve some parameters: */
     3994        this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
     3995        this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
     3996        this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
     3997        this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
     3998        this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
     3999
     4000
     4001        /*Get out if shelf*/
     4002        if(shelf)return;
     4003
     4004        /* Get node coordinates and dof list: */
     4005        GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
     4006        GetDofList1(&doflist1[0]);
     4007
     4008        /*Build frictoin element, needed later: */
     4009        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4010        friction=new Friction("2d",inputs,matpar);
     4011
     4012        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     4013        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
     4014
     4015        /* Start  looping on the number of gaussian points: */
     4016        for (ig=0; ig<num_gauss; ig++){
     4017                /*Pick up the gaussian point: */
     4018                gauss_weight=*(gauss_weights+ig);
     4019                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     4020                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     4021                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     4022
     4023                /*Build alpha_complement_list: */
     4024                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
     4025                else alpha_complement=0;
     4026       
     4027                /*Recover alpha_complement and k: */
     4028                inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
     4029
     4030                /*recover lambda and mu: */
     4031                inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
     4032                inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
     4033                       
     4034                /*recover vx and vy: */
     4035                inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
     4036                inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
     4037
     4038                /* Get Jacobian determinant: */
     4039                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     4040               
     4041                /* Get nodal functions value at gaussian point:*/
     4042                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     4043
     4044                /*Get nodal functions derivatives*/
     4045                GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
     4046
     4047                /*Get k derivative: dk/dx */
     4048                inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
     4049
     4050                /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
     4051                for (i=0;i<numgrids;i++){
     4052
     4053                        //standard term dJ/dki
     4054                        grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
     4055
     4056                        //noise dampening d/dki(1/2*(dk/dx)^2)
     4057                        grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
     4058                       
     4059                        //min dampening
     4060                        if(drag<cm_mindmp_value){
     4061                                grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4062                        }
     4063
     4064                        //max dampening
     4065                        if(drag>cm_maxdmp_value){
     4066                                grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
     4067                        }
     4068                }
     4069               
     4070                /*Add gradje_g_gaussian vector to gradje_g: */
     4071                for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
     4072        }
     4073
     4074        /*Add grade_g to global vector grad_g: */
     4075        VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
     4076       
     4077        cleanup_and_return:
     4078        xfree((void**)&first_gauss_area_coord);
     4079        xfree((void**)&second_gauss_area_coord);
     4080        xfree((void**)&third_gauss_area_coord);
     4081        xfree((void**)&gauss_weights);
     4082        delete friction;
     4083
     4084}
     4085/*}}}*/
     4086/*FUNCTION Tria::GradjDragStokes {{{1*/
     4087void  Tria::GradjDragStokes(Vec grad_g,int analysis_type,int sub_analysis_type){
     4088
     4089        int i;
     4090
     4091        /* node data: */
     4092        const int    numgrids=3;
     4093        const int    NDOF2=2;
     4094        double       xyz_list[numgrids][3];
     4095        int          doflist1[numgrids];
     4096        double       dh1dh3[NDOF2][numgrids];
     4097
     4098        /* grid data: */
    39564099        double drag;
     4100        double alpha_complement;
     4101        Friction* friction=NULL;
    39574102
    39584103        /* gaussian points: */
    … …  
    39674112
    39684113        /* parameters: */
     4114        double  vx,vy,vz;
     4115        double  lambda,mu,xi;
     4116        double  bed,thickness,Neff;
     4117        double  surface_normal[3];
     4118        double  bed_normal[3];
    39694119        double  dk[NDOF2];
    3970         double  vx,vy;
    3971         double  lambda,mu;
    3972         double  bed,thickness,Neff;
    3973        
    3974         /*drag: */
    3975         double alpha_complement_list[numgrids];
    3976         double alpha_complement;
    39774120
    39784121        /*element vector at the gaussian points: */
    … …  
    39914134        /*inputs: */
    39924135        bool shelf;
     4136        int  drag_type;
    39934137
    39944138        /*parameters: */
    … …  
    40014145        /*retrieve inputs :*/
    40024146        inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
     4147        inputs->GetParameterValue(&drag_type,DragTypeEnum);
    40034148
    40044149        /*retrieve some parameters: */
    … …  
    40094154        this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
    40104155
    4011 
    40124156        /*Get out if shelf*/
    40134157        if(shelf)return;
    … …  
    40174161        GetDofList1(&doflist1[0]);
    40184162
    4019         /*Recover inputs: */
    4020         inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
    4021         inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
    4022         inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
    4023         inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
    4024         inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
    4025         inputs->GetParameterValue(&drag_p,DragPEnum);
    4026         inputs->GetParameterValue(&drag_q,DragQEnum);
     4163        /*Build frictoin element, needed later: */
    40274164        inputs->GetParameterValue(&drag_type,DragTypeEnum);
     4165        friction=new Friction("2d",inputs,matpar);
    40284166
    40294167        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    40384176                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    40394177
    4040                 /*Build alpha_complement_list: */
    4041                 if (drag_type==2){
    4042                
    4043                         /*Allocate friction object: */
    4044                         Friction* friction=NewFriction();
    4045                        
    4046                         /*Initialize all fields: */
    4047                         friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    4048                         strcpy(friction->element_type,"2d");
    4049                        
    4050                        
    4051                         friction->gravity=matpar->GetG();
    4052                         friction->rho_ice=matpar->GetRhoIce();
    4053                         friction->rho_water=matpar->GetRhoWater();
    4054                         friction->K=&dragcoefficient_list[0];
    4055                         friction->bed=&bed_list[0];
    4056                         friction->thickness=&thickness_list[0];
    4057                         friction->vx=&vx_list[0];
    4058                         friction->vy=&vy_list[0];
    4059                         friction->p=drag_p;
    4060                         friction->q=drag_q;
    4061 
    4062                        
    4063                         if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
    4064 
    4065                         /*Compute alpha complement: */
    4066                         FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
    4067 
    4068                         /*Erase friction object: */
    4069                         DeleteFriction(&friction);
    4070                 }
    4071                 else{
    4072                         alpha_complement_list[0]=0;
    4073                         alpha_complement_list[1]=0;
    4074                         alpha_complement_list[2]=0;
    4075                 }
    4076        
    4077                 /*Recover alpha_complement and k: */
    4078                 GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
    4079                 inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);
    4080 
    4081                 /*recover lambda and mu: */
    4082                 inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);
    4083                 inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);
    4084                        
    4085                 /*recover vx and vy: */
    4086                 inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum);
    4087                 inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum);
    4088 
    4089                 /* Get Jacobian determinant: */
    4090                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    4091                
    4092                 /* Get nodal functions value at gaussian point:*/
    4093                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    4094 
    4095                 /*Get nodal functions derivatives*/
    4096                 GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss_l1l2l3);
    4097 
    4098                 /*Get k derivative: dk/dx */
    4099                 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);
    4100 
    4101                 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
    4102                 for (i=0;i<numgrids;i++){
    4103 
    4104                         //standard term dJ/dki
    4105                         grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i];
    4106 
    4107                         //noise dampening d/dki(1/2*(dk/dx)^2)
    4108                         grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
    4109                        
    4110                         //min dampening
    4111                         if(drag<cm_mindmp_value){
    4112                                 grade_g_gaussian[i]+=cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4113                         }
    4114 
    4115                         //max dampening
    4116                         if(drag>cm_maxdmp_value){
    4117                                 grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i];
    4118                         }
    4119                 }
    4120                
    4121                 /*Add gradje_g_gaussian vector to gradje_g: */
    4122                 for( i=0; i<numgrids; i++)grade_g[i]+=grade_g_gaussian[i];
    4123         }
    4124 
    4125         /*Add grade_g to global vector grad_g: */
    4126         VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    4127        
    4128         cleanup_and_return:
    4129         xfree((void**)&first_gauss_area_coord);
    4130         xfree((void**)&second_gauss_area_coord);
    4131         xfree((void**)&third_gauss_area_coord);
    4132         xfree((void**)&gauss_weights);
    4133 
    4134 }
    4135 /*}}}*/
    4136 /*FUNCTION Tria::GradjDragStokes {{{1*/
    4137 void  Tria::GradjDragStokes(Vec grad_g,int analysis_type,int sub_analysis_type){
    4138 
    4139         int i;
    4140 
    4141         /* node data: */
    4142         const int    numgrids=3;
    4143         const int    NDOF2=2;
    4144         double       xyz_list[numgrids][3];
    4145         int          doflist1[numgrids];
    4146         double       dh1dh3[NDOF2][numgrids];
    4147 
    4148         /* grid data: */
    4149         double vx_list[numgrids];
    4150         double vy_list[numgrids];
    4151         double vz_list[numgrids];
    4152         double drag;
    4153         double  thickness_list[numgrids];
    4154         double  bed_list[numgrids];
    4155         double  dragcoefficient_list[numgrids];
    4156         double  drag_p,drag_q;
    4157         double alpha_complement_list[numgrids];
    4158         double alpha_complement;
    4159 
    4160         /* gaussian points: */
    4161         int     num_gauss,ig;
    4162         double* first_gauss_area_coord  =  NULL;
    4163         double* second_gauss_area_coord =  NULL;
    4164         double* third_gauss_area_coord  =  NULL;
    4165         double* gauss_weights           =  NULL;
    4166         double  gauss_weight;
    4167         double  gauss_l1l2l3[3];
    4168         double  gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}};
    4169 
    4170         /* parameters: */
    4171         double  vx,vy,vz;
    4172         double  lambda,mu,xi;
    4173         double  bed,thickness,Neff;
    4174         double  surface_normal[3];
    4175         double  bed_normal[3];
    4176         double  dk[NDOF2];
    4177 
    4178         /*element vector at the gaussian points: */
    4179         double  grade_g[numgrids]={0.0};
    4180         double  grade_g_gaussian[numgrids];
    4181 
    4182         /* Jacobian: */
    4183         double Jdet;
    4184 
    4185         /*nodal functions: */
    4186         double l1l2l3[3];
    4187 
    4188         /* strain rate: */
    4189         double epsilon[3]; /* epsilon=[exx,eyy,exy];*/
    4190 
    4191         /*inputs: */
    4192         bool shelf;
    4193         int  drag_type;
    4194 
    4195         /*parameters: */
    4196         double  cm_noisedmp;
    4197         double  cm_mindmp_slope;
    4198         double  cm_mindmp_value;
    4199         double  cm_maxdmp_value;
    4200         double  cm_maxdmp_slope;
    4201 
    4202         /*retrieve inputs :*/
    4203         inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
    4204         inputs->GetParameterValue(&drag_type,DragTypeEnum);
    4205 
    4206         /*retrieve some parameters: */
    4207         this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum);
    4208         this->parameters->FindParam(&cm_mindmp_value,CmMinDmpValueEnum);
    4209         this->parameters->FindParam(&cm_mindmp_slope,CmMinDmpSlopeEnum);
    4210         this->parameters->FindParam(&cm_maxdmp_value,CmMaxDmpValueEnum);
    4211         this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum);
    4212 
    4213         /*Get out if shelf*/
    4214         if(shelf)return;
    4215 
    4216         /* Get node coordinates and dof list: */
    4217         GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
    4218         GetDofList1(&doflist1[0]);
    4219 
    4220         /*Build alpha_complement_list: */
    4221         if (drag_type==2){
    4222 
    4223                 /*Allocate friction object: */
    4224                 Friction* friction=NewFriction();
    4225 
    4226                 inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum);
    4227                 inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum);
    4228                 inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum);
    4229                 inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum);
    4230                 inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum);
    4231                 inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum);
    4232                 inputs->GetParameterValue(&drag_p,DragPEnum);
    4233                 inputs->GetParameterValue(&drag_q,DragQEnum);
    4234 
    4235 
    4236                 /*Initialize all fields: */
    4237                 friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    4238                 strcpy(friction->element_type,"2d");
    4239 
    4240                 friction->gravity=matpar->GetG();
    4241                 friction->rho_ice=matpar->GetRhoIce();
    4242                 friction->rho_water=matpar->GetRhoWater();
    4243                 friction->K=&dragcoefficient_list[0];
    4244                 friction->bed=&bed_list[0];
    4245                 friction->thickness=&thickness_list[0];
    4246                 friction->vx=&vx_list[0];
    4247                 friction->vy=&vy_list[0];
    4248                 friction->p=drag_p;
    4249                 friction->q=drag_q;
    4250 
    4251 
    4252                 if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!");
    4253 
    4254                 /*Compute alpha complement: */
    4255                 FrictionGetAlphaComplement(&alpha_complement_list[0],friction);
    4256 
    4257                 /*Erase friction object: */
    4258                 DeleteFriction(&friction);
    4259         }
    4260         else{
    4261                 alpha_complement_list[0]=0;
    4262                 alpha_complement_list[1]=0;
    4263                 alpha_complement_list[2]=0;
    4264         }
    4265 
    4266         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    4267         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
    4268 
    4269         /* Start  looping on the number of gaussian points: */
    4270         for (ig=0; ig<num_gauss; ig++){
    4271                 /*Pick up the gaussian point: */
    4272                 gauss_weight=*(gauss_weights+ig);
    4273                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    4274                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    4275                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    4276 
    4277                 /*Recover alpha_complement and k: */
    4278                 GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
     4178                /*Recover alpha_complement and drag: */
     4179                if (drag_type==2) friction->GetAlphaComplement(&alpha_complement, gauss_l1l2l3,VxAverageEnum,VyAverageEnum);
     4180                else alpha_complement=0;
    42794181                inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum);
    42804182
    … …  
    43384240        VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
    43394241
    4340 cleanup_and_return:
     4242        cleanup_and_return:
    43414243        xfree((void**)&first_gauss_area_coord);
    43424244        xfree((void**)&second_gauss_area_coord);
    43434245        xfree((void**)&third_gauss_area_coord);
    43444246        xfree((void**)&gauss_weights);
     4247        delete friction;
    43454248
    43464249}

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

    r3775 r3833  
    1 /*!\file: Friction.cpp
    2  * \brief: wrapper for all friction parameters
    3  */ 
     1/*!\file Friction.c
     2 * \brief: implementation of the Friction object
     3 */
    44
     5/*Headers:*/
     6/*{{{1*/
     7#ifdef HAVE_CONFIG_H
     8        #include "config.h"
     9#else
     10#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
     11#endif
     12
     13#include "stdio.h"
     14#include <string.h>
    515#include "../objects.h"
     16#include "../../EnumDefinitions/EnumDefinitions.h"
    617#include "../../shared/shared.h"
    718#include "../../include/include.h"
     19/*}}}*/
    820
    9 /*--------------------------------------------------
    10         NewFriction
    11   --------------------------------------------------*/
    12 
    13 /*FUNCTION NewFriction{{{1*/
    14 Friction* NewFriction(void)
    15 {
    16         /* create a new Friction object */
    17 
    18         return (Friction*)xmalloc(sizeof(Friction));
     21/*methods: */
     22/*FUNCTION Friction::Friction() {{{1*/
     23Friction::Friction(){
     24        this->element_type=NULL;
     25        this->inputs=NULL;
     26        this->matpar=NULL;
    1927}
    2028/*}}}*/
    21 /*FUNCTION FrictionEcho{{{1*/
    22 void  FrictionEcho(Friction* friction){
     29/*FUNCTION Friction::Friction(char* element_type, Inputs* inputs,Matpar* matpar){{{1*/
     30Friction::Friction(char* element_type_in,Inputs* inputs_in,Matpar* matpar_in){
    2331
    24         int i;
    25        
    26         printf("Friction echo: \n");
    27         printf("   element_type: %s\n",friction->element_type);
    28         printf("   gravity: %g\n",friction->gravity);
    29         printf("   rho_ice: %g\n",friction->rho_ice);
    30         printf("   rho_water: %g\n",friction->rho_water);
    31         printf("   p: %g\n",friction->p);
    32         printf("   q: %g\n",friction->q);
    33         printf("   analysis_type: %i\n",friction->analysis_type);
    34 
    35         printf("K: ");
    36         for(i=0;i<3;i++)printf("%g ",friction->K[i]);
    37         printf("\n");
    38 
    39         printf("bed: ");
    40         for(i=0;i<3;i++)printf("%g ",friction->bed[i]);
    41         printf("\n");
    42 
    43         printf("thickness: ");
    44         for(i=0;i<3;i++)printf("%g ",friction->thickness[i]);
    45         printf("\n");
    46 
    47         printf("vx: ");
    48         for(i=0;i<3;i++)printf("%g ",friction->vx[i]);
    49         printf("\n            ");
    50         printf("vy: ");
    51         for(i=0;i<3;i++)printf("%g ",friction->vy[i]);
    52         printf("\n");
    53         if(friction->vz){
    54                 printf("vz: ");
    55                 for(i=0;i<3;i++)printf("%g ",friction->vz[i]);
    56                 printf("\n");
    57         }
     32        this->inputs=inputs_in;
     33        this->element_type=(char*)xmalloc((strlen(element_type_in)+1)*sizeof(char));
     34        strcpy(this->element_type,element_type);
     35        this->matpar=matpar_in;
    5836
    5937}
    6038/*}}}*/
    61 /*FUNCTION FrictionInit {{{1*/
    62 /*--------------------------------------------------
    63         FrictionInit
    64   --------------------------------------------------*/
    65 
    66 int FrictionInit(Friction* friction)
    67 {
    68         friction->element_type=NULL;
    69         friction->gravity=UNDEF;
    70         friction->rho_ice=UNDEF;
    71         friction->rho_water=UNDEF;
    72         friction->K=NULL;
    73         friction->bed=NULL;
    74         friction->thickness=NULL;
    75         friction->vx=NULL;
    76         friction->vy=NULL;
    77         friction->vz=NULL;
    78         friction->p=UNDEF;
    79         friction->q=UNDEF;
    80         friction->analysis_type=UNDEF;
    81        
    82         return 1;
     39/*FUNCTION Friction::~Friction() {{{1*/
     40Friction::~Friction(){
     41        xfree((void**)&element_type);
    8342}
    8443/*}}}*/
    85 /*FUNCTION DeleteFriction {{{1*/
    86 /*--------------------------------------------------
    87         DeleteFriction
    88   --------------------------------------------------*/
    89 
    90 void DeleteFriction(Friction** pfriction)
    91 {
    92         Friction* friction = *pfriction;
    93 
    94         /*Just  erase element_type: */
    95         xfree((void**)&friction->element_type);
    96 
    97         /*Erase entire structure: */
    98         xfree((void**)pfriction);
     44/*FUNCTION Friction::Echo {{{1*/
     45void Friction::Echo(void){
     46        printf("Friction:\n");
     47        printf("   element_type: %s\n",this->element_type);
     48        inputs->Echo();
     49        matpar->Echo();
    9950}
    10051/*}}}*/
    101 /*FUNCTION FrictionGetAlpha2 {{{1*/
    102 /*--------------------------------------------------
    103         FrictionGetAlpha2
    104   --------------------------------------------------*/
    105 void  FrictionGetAlpha2(double* alpha2, Friction* friction){
     52/*FUNCTION Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){{{1*/
     53void Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){
    10654
    10755        /*This routine calculates the basal friction coefficient
    108         alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p*
    109         alpha2 is assumed to be double[3]: */
     56        alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/
    11057
    11158        /*diverse: */
    112         int     i;
    113         const int     numgrids=3;
    114         double  Neff[numgrids];
    11559        double  r,s;
    116         double  velocity_x[numgrids];
    117         double  velocity_y[numgrids];
    118         double  velocity_z[numgrids];
    119         double  velocity_mag[numgrids];
     60        double  drag_p, drag_q;
     61        double  gravity,rho_ice,rho_water;
     62        double  Neff;
     63        double  thickness,bed;
     64        double  vx,vy,vz,vmag;
     65        double  drag_coefficient;
     66        double  alpha2;
     67
     68
     69        /*Recover parameters: */
     70        inputs->GetParameterValue(&drag_p,DragPEnum);
     71        inputs->GetParameterValue(&drag_q,DragQEnum);
     72        inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
     73        inputs->GetParameterValue(&bed, gauss,BedEnum);
     74        inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     75
     76        /*Get material parameters: */
     77        gravity=matpar->GetG();
     78        rho_ice=matpar->GetRhoIce();
     79        rho_water=matpar->GetRhoWater();
    12080
    12181        //compute r and q coefficients: */
    122         r=friction->q/friction->p;
    123         s=1./friction->p;
     82        r=drag_q/drag_p;
     83        s=1./drag_p;
    12484               
    12585        //From bed and thickness, compute effective pressure when drag is viscous:
    126         for(i=0;i<numgrids;i++){
     86        Neff=gravity*(rho_ice*thickness+rho_water*bed);
     87       
     88        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     89        the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     90        for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     91        replace it by Neff=0 (ie, equival it to an ice shelf)*/
     92        if (Neff<0)Neff=0;
    12793
    128                 Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);
     94        if(strcmp(element_type,"2d")){
     95                inputs->GetParameterValue(&vx, gauss,vxenum);
     96                inputs->GetParameterValue(&vy, gauss,vyenum);
     97                vmag=sqrt(pow(vx,2)+pow(vy,2));
     98        }
     99        else{
     100                inputs->GetParameterValue(&vx, gauss,vxenum);
     101                inputs->GetParameterValue(&vy, gauss,vyenum);
     102                inputs->GetParameterValue(&vz, gauss,vzenum);
     103                vmag=sqrt(pow(vx,2)+pow(vy,2)+pow(vz,2));
     104        }
     105       
     106        alpha2=pow(drag_coefficient,2)*pow(Neff,r)*pow(vmag,(s-1));
    129107
    130                 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    131                 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    132                 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    133                 replace it by Neff=0 (ie, equival it to an ice shelf)*/
    134                 if (Neff[i]<0)Neff[i]=0;
    135 
    136                 //We need the velocity magnitude to evaluate the basal stress:
    137                 if(strcmp(friction->element_type,"2d")){
    138                         velocity_x[i]=friction->vx[i];//velocities of size numgridsx2
    139                         velocity_y[i]=friction->vy[i];
    140                         velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
    141                 }
    142                 else{
    143                         velocity_x[i]=friction->vx[i]; //velocities of size numgridsx3
    144                         velocity_y[i]=friction->vy[i];
    145                         velocity_z[i]=friction->vz[i];
    146                         velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)+pow(velocity_z[i],2));
    147                 }
    148        
    149                 alpha2[i]=pow(friction->K[i],2)*pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
    150         }
     108        /*Assign output pointers:*/
     109        *palpha2=alpha2;
    151110}
    152111/*}}}*/
    153 /*FUNCTION FrictionGetAlphaComplement {{{1*/
    154 /*--------------------------------------------------
    155         FrictionGetAlphaComplement
    156   --------------------------------------------------*/
    157 void  FrictionGetAlphaComplement(double* alpha_complement, Friction* friction){
    158 
     112/*FUNCTION Friction::GetAlphaComplement {{{1*/
     113void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){
     114       
    159115        /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p.
    160116         * FrictionGetAlphaComplement is used in control methods on drag, and it computes:
    … …  
    163119        /*diverse: */
    164120        int     i;
    165         const int     numgrids=3;
    166         double  Neff[numgrids];
     121        double  Neff;
    167122        double  r,s;
    168         double  velocity_x[numgrids];
    169         double  velocity_y[numgrids];
    170         double  velocity_mag[numgrids];
     123        double  vx;
     124        double  vy;
     125        double  vmag;
     126        double  drag_p,drag_q;
     127        double  drag_coefficient;
     128        double  bed,thickness;
     129        double  gravity,rho_ice,rho_water;
     130        double  alpha_complement;
     131
     132        /*Recover parameters: */
     133        inputs->GetParameterValue(&drag_p,DragPEnum);
     134        inputs->GetParameterValue(&drag_q,DragQEnum);
     135        inputs->GetParameterValue(&thickness, gauss,ThicknessEnum);
     136        inputs->GetParameterValue(&bed, gauss,BedEnum);
     137        inputs->GetParameterValue(&drag_coefficient, gauss,DragCoefficientEnum);
     138
     139        /*Get material parameters: */
     140        gravity=matpar->GetG();
     141        rho_ice=matpar->GetRhoIce();
     142        rho_water=matpar->GetRhoWater();
     143
    171144
    172145        //compute r and q coefficients: */
    173         r=friction->q/friction->p;
    174         s=1./friction->p;
     146        r=drag_q/drag_p;
     147        s=1./drag_p;
    175148               
    176149        //From bed and thickness, compute effective pressure when drag is viscous:
    177         for(i=0;i<numgrids;i++){
     150        Neff=gravity*(rho_ice*thickness+rho_water*bed);
    178151
    179                 Neff[i]=friction->gravity*(friction->rho_ice*friction->thickness[i]+friction->rho_water*friction->bed[i]);
     152        /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
     153          the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
     154          for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
     155          replace it by Neff=0 (ie, equival it to an ice shelf)*/
     156        if (Neff<0)Neff=0;
    180157
    181                 /*If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
    182                 the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
    183                 for friction should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
    184                 replace it by Neff=0 (ie, equival it to an ice shelf)*/
    185                 if (Neff[i]<0)Neff[i]=0;
     158        //We need the velocity magnitude to evaluate the basal stress:
     159        inputs->GetParameterValue(&vx, gauss,vxenum);
     160        inputs->GetParameterValue(&vy, gauss,vyenum);
     161        vmag=sqrt(pow(vx,2)+pow(vy,2));
     162       
     163        alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    186164
    187                 //We need the velocity magnitude to evaluate the basal stress:
    188                 velocity_x[i]=friction->vx[i]; //velocities of size numgridsx2
    189                 velocity_y[i]=friction->vy[i];
    190                 velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2));
    191        
    192                 alpha_complement[i]=pow(Neff[i],r)*pow(velocity_mag[i],(s-1));
    193         }
     165        /*Assign output pointers:*/
     166        *palpha_complement=alpha_complement;
     167
    194168}
    195169/*}}}*/

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

    r3683 r3833  
    1 /*! \file Friction.h
    2  *  \brief Friction object, interface, declaration
     1/*!\file Friction.h
     2 * \brief: header file for friction object
    33 */
    44
    5 #ifndef _FRICTION_H
    6 #define _FRICTION_H
     5#ifndef _FRICTION_H_
     6#define _FRICTION_H_
    77
    8 /*!Friction declaration: */
    9 struct Friction {
     8/*Headers:*/
     9/*{{{1*/
     10class Inputs;
     11class Matpar;
     12/*}}}*/
    1013
    11         char*   element_type;
    12         double  gravity;
    13         double  rho_ice;
    14         double  rho_water;
    15         double* K;
    16         double* bed;
    17         double* thickness;
    18         double* vx;
    19         double* vy;
    20         double* vz;
    21         double  p;
    22         double  q;
    23         int     analysis_type;
     14class Friction{
     15
     16        public:
     17
     18                char* element_type;
     19                Inputs* inputs;
     20                Matpar* matpar;
     21
     22                /*methods: */
     23                Friction();
     24                Friction(char* element_type, Inputs* inputs,Matpar* matpar);
     25                ~Friction();
     26       
     27                void  Echo(void);
     28                void  GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum);
     29                void  GetAlphaComplement(double* alpha_complement, double* gauss,int vxenum,int vyenum);
    2430
    2531};
    2632
    27 
    28 /* creation, initialisation: */
    29 
    30         Friction*       NewFriction(void);
    31         int     FrictionInit(Friction* friction);
    32         void  DeleteFriction(Friction** pfriction);
    33         void  FrictionGetAlpha2(double* alpha2, Friction* friction);
    34         void  FrictionGetAlphaComplement(double* alpha_complement, Friction* friction);
    35         void  FrictionEcho(Friction* friction);
    36 
    37 #endif  /* _FRICTION_H */
    38 
     33#endif  /* _FRICTION_H_ */

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

    r3832 r3833  
    2626/*Loads: */
    2727#include "./Loads/Friction.h"
    28 #include "./Loads/Friction2.h"
    2928#include "./Loads/Icefront.h"
    3029#include "./Loads/Numericalflux.h"
Note: See TracChangeset for help on using the changeset viewer.