Changeset 3833
- Timestamp:
- 05/18/10 21:22:39 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 2 deleted
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Makefile.am
r3832 r3833 50 50 ./objects/Loads/Friction.h\ 51 51 ./objects/Loads/Friction.cpp\ 52 ./objects/Loads/Friction2.h\53 ./objects/Loads/Friction2.cpp\54 52 ./objects/DakotaPlugin.h\ 55 53 ./objects/DakotaPlugin.cpp\ … … 499 497 ./objects/Loads/Friction.h\ 500 498 ./objects/Loads/Friction.cpp\ 501 ./objects/Loads/Friction2.h\502 ./objects/Loads/Friction2.cpp\503 499 ./objects/DakotaPlugin.h\ 504 500 ./objects/DakotaPlugin.cpp\ -
issm/trunk/src/c/objects/Elements/Penta.cpp
r3824 r3833 1180 1180 double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ 1181 1181 double viscosity; 1182 double alpha2_list[numgrids2d];1183 1182 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; 1192 1184 1193 1185 /*parameters: */ … … 1281 1273 if((onbed==1) && (shelf==0)){ 1282 1274 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); 1316 1277 1317 1278 for(i=0;i<numgrids2d;i++){ … … 1357 1318 1358 1319 /*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); 1360 1321 1361 1322 DLStokes[0][0]=alpha2_gauss*gauss_weight*Jdet2d; … … 1384 1345 } 1385 1346 } 1347 1348 /*Free ressources:*/ 1349 delete friction; 1350 1386 1351 } //if ( (onbed==1) && (shelf==0)) 1387 1352 -
issm/trunk/src/c/objects/Elements/Tria.cpp
r3832 r3833 1299 1299 1300 1300 /*friction: */ 1301 Friction 2* friction2=NULL;1301 Friction* friction=NULL; 1302 1302 double alpha2; 1303 1303 … … 1324 1324 /*build friction object, used later on: */ 1325 1325 if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); 1326 friction 2=new Friction2("2d",inputs,matpar);1326 friction=new Friction("2d",inputs,matpar); 1327 1327 1328 1328 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 1339 1339 1340 1340 /*Friction: */ 1341 friction 2->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum);1341 friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxAverageEnum,VyAverageEnum,VzAverageEnum); 1342 1342 1343 1343 // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, … … 1380 1380 xfree((void**)&third_gauss_area_coord); 1381 1381 xfree((void**)&gauss_weights); 1382 delete friction 2;1382 delete friction; 1383 1383 1384 1384 } … … 2902 2902 int drag_type; 2903 2903 double basalfriction; 2904 Friction 2* friction2=NULL;2904 Friction* friction=NULL; 2905 2905 double alpha2,vx,vy; 2906 2906 double geothermalflux_value; … … 2936 2936 inputs->GetParameterValue(&drag_type,DragTypeEnum); 2937 2937 if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); 2938 friction 2=new Friction2("3d",inputs,matpar);2938 friction=new Friction("3d",inputs,matpar); 2939 2939 2940 2940 /* Ice/ocean heat exchange flux on ice shelf base */ … … 2957 2957 inputs->GetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); 2958 2958 2959 friction 2->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum);2959 friction->GetAlpha2(&alpha2,&gauss_coord[0],VxAverageEnum,VyAverageEnum,VzAverageEnum); 2960 2960 inputs->GetParameterValue(&vx, &gauss_coord[0],VxAverageEnum); 2961 2961 inputs->GetParameterValue(&vy, &gauss_coord[0],VyAverageEnum); … … 2981 2981 xfree((void**)&third_gauss_area_coord); 2982 2982 xfree((void**)&gauss_weights); 2983 delete friction 2;2983 delete friction; 2984 2984 2985 2985 } … … 3943 3943 3944 3944 /* grid data: */ 3945 double vx_list[numgrids];3946 double vy_list[numgrids];3947 3945 double adjx_list[numgrids]; 3948 3946 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*/ 4087 void 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: */ 3956 4099 double drag; 4100 double alpha_complement; 4101 Friction* friction=NULL; 3957 4102 3958 4103 /* gaussian points: */ … … 3967 4112 3968 4113 /* 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]; 3969 4119 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;3977 4120 3978 4121 /*element vector at the gaussian points: */ … … 3991 4134 /*inputs: */ 3992 4135 bool shelf; 4136 int drag_type; 3993 4137 3994 4138 /*parameters: */ … … 4001 4145 /*retrieve inputs :*/ 4002 4146 inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); 4147 inputs->GetParameterValue(&drag_type,DragTypeEnum); 4003 4148 4004 4149 /*retrieve some parameters: */ … … 4009 4154 this->parameters->FindParam(&cm_maxdmp_slope,CmMaxDmpSlopeEnum); 4010 4155 4011 4012 4156 /*Get out if shelf*/ 4013 4157 if(shelf)return; … … 4017 4161 GetDofList1(&doflist1[0]); 4018 4162 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: */ 4027 4164 inputs->GetParameterValue(&drag_type,DragTypeEnum); 4165 friction=new Friction("2d",inputs,matpar); 4028 4166 4029 4167 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ … … 4038 4176 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); 4039 4177 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; 4279 4181 inputs->GetParameterValue(&drag, &gauss_l1l2l3[0],DragCoefficientEnum); 4280 4182 … … 4338 4240 VecSetValues(grad_g,numgrids,doflist1,(const double*)grade_g,ADD_VALUES); 4339 4241 4340 cleanup_and_return:4242 cleanup_and_return: 4341 4243 xfree((void**)&first_gauss_area_coord); 4342 4244 xfree((void**)&second_gauss_area_coord); 4343 4245 xfree((void**)&third_gauss_area_coord); 4344 4246 xfree((void**)&gauss_weights); 4247 delete friction; 4345 4248 4346 4249 } -
issm/trunk/src/c/objects/Loads/Friction.cpp
r3775 r3833 1 /*!\file : Friction.cpp2 * \brief: wrapper for all friction parameters3 */ 1 /*!\file Friction.c 2 * \brief: implementation of the Friction object 3 */ 4 4 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> 5 15 #include "../objects.h" 16 #include "../../EnumDefinitions/EnumDefinitions.h" 6 17 #include "../../shared/shared.h" 7 18 #include "../../include/include.h" 19 /*}}}*/ 8 20 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*/ 23 Friction::Friction(){ 24 this->element_type=NULL; 25 this->inputs=NULL; 26 this->matpar=NULL; 19 27 } 20 28 /*}}}*/ 21 /*FUNCTION Friction Echo{{{1*/22 void FrictionEcho(Friction* friction){29 /*FUNCTION Friction::Friction(char* element_type, Inputs* inputs,Matpar* matpar){{{1*/ 30 Friction::Friction(char* element_type_in,Inputs* inputs_in,Matpar* matpar_in){ 23 31 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; 58 36 59 37 } 60 38 /*}}}*/ 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*/ 40 Friction::~Friction(){ 41 xfree((void**)&element_type); 83 42 } 84 43 /*}}}*/ 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*/ 45 void Friction::Echo(void){ 46 printf("Friction:\n"); 47 printf(" element_type: %s\n",this->element_type); 48 inputs->Echo(); 49 matpar->Echo(); 99 50 } 100 51 /*}}}*/ 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*/ 53 void Friction::GetAlpha2(double* palpha2, double* gauss,int vxenum,int vyenum,int vzenum){ 106 54 107 55 /*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**/ 110 57 111 58 /*diverse: */ 112 int i;113 const int numgrids=3;114 double Neff[numgrids];115 59 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(); 120 80 121 81 //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; 124 84 125 85 //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; 127 93 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)); 129 107 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; 151 110 } 152 111 /*}}}*/ 153 /*FUNCTION FrictionGetAlphaComplement {{{1*/ 154 /*-------------------------------------------------- 155 FrictionGetAlphaComplement 156 --------------------------------------------------*/ 157 void FrictionGetAlphaComplement(double* alpha_complement, Friction* friction){ 158 112 /*FUNCTION Friction::GetAlphaComplement {{{1*/ 113 void Friction::GetAlphaComplement(double* palpha_complement, double* gauss,int vxenum,int vyenum){ 114 159 115 /* 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. 160 116 * FrictionGetAlphaComplement is used in control methods on drag, and it computes: … … 163 119 /*diverse: */ 164 120 int i; 165 const int numgrids=3; 166 double Neff[numgrids]; 121 double Neff; 167 122 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 171 144 172 145 //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; 175 148 176 149 //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); 178 151 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; 180 157 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)); 186 164 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 194 168 } 195 169 /*}}}*/ -
issm/trunk/src/c/objects/Loads/Friction.h
r3683 r3833 1 /*! 2 * \brief Friction object, interface, declaration1 /*!\file Friction.h 2 * \brief: header file for friction object 3 3 */ 4 4 5 #ifndef _FRICTION_H 6 #define _FRICTION_H 5 #ifndef _FRICTION_H_ 6 #define _FRICTION_H_ 7 7 8 /*!Friction declaration: */ 9 struct Friction { 8 /*Headers:*/ 9 /*{{{1*/ 10 class Inputs; 11 class Matpar; 12 /*}}}*/ 10 13 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; 14 class 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); 24 30 25 31 }; 26 32 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 26 26 /*Loads: */ 27 27 #include "./Loads/Friction.h" 28 #include "./Loads/Friction2.h"29 28 #include "./Loads/Icefront.h" 30 29 #include "./Loads/Numericalflux.h"
Note:
See TracChangeset
for help on using the changeset viewer.