Changeset 4838
- Timestamp:
- 07/27/10 15:20:20 (15 years ago)
- Location:
- issm/trunk/src/c/objects/Elements
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/objects/Elements/Penta.cpp
r4835 r4838 3024 3024 int numberofdofspernode; 3025 3025 3026 3027 3026 /* 3d gaussian points: */ 3028 3027 int num_gauss,ig; … … 3097 3096 bool shelf; 3098 3097 3098 /*inputs: */ 3099 Input* vx_input=NULL; 3100 Input* vy_input=NULL; 3101 Input* vxold_input=NULL; 3102 Input* vyold_input=NULL; 3103 3099 3104 /*retrieve inputs :*/ 3100 3105 inputs->GetParameterValue(&onwater,ElementOnWaterEnum); … … 3112 3117 bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build 3113 3118 the stiffness matrix. */ 3114 3115 3119 3116 3120 if ((collapse==1) && (onbed==0)){ … … 3136 3140 GetDofList(&doflist[0],&numberofdofspernode); 3137 3141 3138 3139 3142 /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore 3140 3143 get tria gaussian points as well as segment gaussian points. For tria gaussian … … 3148 3151 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); 3149 3152 3153 /*Retrieve all inputs we will be needing: */ 3154 vx_input=inputs->GetInput(VxEnum); 3155 vy_input=inputs->GetInput(VyEnum); 3156 vxold_input=inputs->GetInput(VxOldEnum); 3157 vyold_input=inputs->GetInput(VyOldEnum); 3158 3150 3159 /* Start looping on the number of gaussian points: */ 3151 3160 for (ig1=0; ig1<num_area_gauss; ig1++){ … … 3165 3174 3166 3175 /*Get strain rate from velocity: */ 3167 inputs->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum);3168 inputs->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum);3176 this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); 3177 this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); 3169 3178 3170 3179 /*Get viscosity: */ … … 3344 3353 bool isstokes; 3345 3354 3355 /*inputs: */ 3356 Input* vx_input=NULL; 3357 Input* vy_input=NULL; 3358 Input* vz_input=NULL; 3359 3346 3360 /*retrive parameters: */ 3347 3361 parameters->FindParam(&analysis_type,AnalysisTypeEnum); … … 3379 3393 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 3380 3394 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); 3395 3396 /*Retrieve all inputs we will be needing: */ 3397 vx_input=inputs->GetInput(VxEnum); 3398 vy_input=inputs->GetInput(VyEnum); 3399 vz_input=inputs->GetInput(VzEnum); 3381 3400 3382 3401 /* Start looping on the number of gaussian points: */ … … 3393 3412 3394 3413 /*Compute strain rate: */ 3395 inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);3414 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 3396 3415 3397 3416 /*Get viscosity: */ … … 3462 3481 3463 3482 /*Compute strain rate: */ 3464 inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);3483 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 3465 3484 3466 3485 /*Get viscosity at last iteration: */ … … 4094 4113 bool collapse; 4095 4114 bool onbed; 4115 Input* surface_input=NULL; 4116 Input* thickness_input=NULL; 4096 4117 4097 4118 /*retrieve inputs :*/ … … 4135 4156 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); 4136 4157 4158 /*Retrieve all inputs we will be needing: */ 4159 thickness_input=inputs->GetInput(ThicknessEnum); 4160 surface_input=inputs->GetInput(SurfaceEnum); 4161 4137 4162 /* Start looping on the number of gaussian points: */ 4138 4163 for (ig1=0; ig1<num_area_gauss; ig1++){ … … 4150 4175 4151 4176 /*Compute thickness at gaussian point: */ 4152 inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum);4177 thickness_input->GetParameterValue(&thickness, gauss_coord); 4153 4178 4154 4179 /*Compute slope at gaussian point: */ 4155 inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum);4180 surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord); 4156 4181 4157 4182 /* Get Jacobian determinant: */ … … 4293 4318 bool shelf; 4294 4319 bool isstokes; 4320 Input* vx_input=NULL; 4321 Input* vy_input=NULL; 4322 Input* vz_input=NULL; 4323 Input* bed_input=NULL; 4295 4324 4296 4325 /*retrieve inputs :*/ … … 4326 4355 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 4327 4356 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); 4357 4358 /*Retrieve all inputs we will be needing: */ 4359 vx_input=inputs->GetInput(VxEnum); 4360 vy_input=inputs->GetInput(VyEnum); 4361 vz_input=inputs->GetInput(VzEnum); 4362 bed_input=inputs->GetInput(BedEnum); 4328 4363 4329 4364 /* Start looping on the number of gaussian points: */ … … 4340 4375 4341 4376 /*Compute strain rate and viscosity: */ 4342 inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);4377 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 4343 4378 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 4344 4379 … … 4421 4456 4422 4457 /* Get bed at gaussian point */ 4423 inputs->GetParameterValue(&bed, gauss_coord,BedEnum);4458 bed_input->GetParameterValue(&bed, gauss_coord); 4424 4459 4425 4460 /*Get L matrix : */ … … 4514 4549 bool onwater; 4515 4550 bool onbed; 4551 Input* vx_input=NULL; 4552 Input* vy_input=NULL; 4516 4553 4517 4554 /*retrieve inputs :*/ … … 4541 4578 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss); 4542 4579 4580 /*Retrieve all inputs we will be needing: */ 4581 vx_input=inputs->GetInput(VxEnum); 4582 vy_input=inputs->GetInput(VyEnum); 4583 4543 4584 /* Start looping on the number of gaussian points: */ 4544 4585 for (ig1=0; ig1<num_area_gauss; ig1++){ … … 4557 4598 /*Get velocity derivative, with respect to x and y: */ 4558 4599 4559 inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum);4560 inputs->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord,VyEnum);4600 vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord); 4601 vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord); 4561 4602 dudx=du[0]; 4562 4603 dvdy=dv[1]; … … 4725 4766 bool onbed; 4726 4767 bool shelf; 4768 Input* vx_input=NULL; 4769 Input* vy_input=NULL; 4770 Input* vz_input=NULL; 4727 4771 4728 4772 /*retrieve inputs :*/ … … 4757 4801 /*Get gaussian points: */ 4758 4802 GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss); 4803 4804 /*Retrieve all inputs we will be needing: */ 4805 vx_input=inputs->GetInput(VxEnum); 4806 vy_input=inputs->GetInput(VyEnum); 4807 vz_input=inputs->GetInput(VzEnum); 4759 4808 4760 4809 /* Start looping on the number of gaussian points: */ … … 4771 4820 4772 4821 /*Compute strain rate and viscosity: */ 4773 inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);4822 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 4774 4823 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 4775 4824 … … 6014 6063 } 6015 6064 /*}}}*/ 6065 /*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ 6066 void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){ 6067 /*Compute the 3d Blatter/PattynStrain Rate (5 components): 6068 * 6069 * epsilon=[exx eyy exy exz eyz] 6070 * 6071 * with exz=1/2 du/dz 6072 * eyz=1/2 dv/dz 6073 * 6074 * the contribution of vz is neglected 6075 */ 6076 6077 int i; 6078 6079 double epsilonvx[5]; 6080 double epsilonvy[5]; 6081 6082 /*Check that both inputs have been found*/ 6083 if (!vx_input || !vy_input){ 6084 ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input); 6085 } 6086 6087 /*Get strain rate assuming that epsilon has been allocated*/ 6088 vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); 6089 vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); 6090 6091 /*Sum all contributions*/ 6092 for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; 6093 6094 } 6095 /*}}}*/ 6096 /*FUNCTION Penta::GetStrainRate3d{{{1*/ 6097 void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ 6098 /*Compute the 3d Strain Rate (6 components): 6099 * 6100 * epsilon=[exx eyy ezz exy exz eyz] 6101 */ 6102 6103 int i; 6104 6105 double epsilonvx[6]; 6106 double epsilonvy[6]; 6107 double epsilonvz[6]; 6108 6109 /*Check that both inputs have been found*/ 6110 if (!vx_input || !vy_input || !vz_input){ 6111 ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input); 6112 } 6113 6114 /*Get strain rate assuming that epsilon has been allocated*/ 6115 vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); 6116 vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); 6117 vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); 6118 6119 /*Sum all contributions*/ 6120 for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; 6121 6122 } 6123 /*}}}*/ 6016 6124 /*FUNCTION Penta::GradjB {{{1*/ 6017 6125 void Penta::GradjB(Vec gradient){ -
issm/trunk/src/c/objects/Elements/Penta.h
r4835 r4838 168 168 void GetSolutionFromInputsDiagnosticVert(Vec solutiong); 169 169 void GetSolutionFromInputsThermal(Vec solutiong); 170 void GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);171 void GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord);170 void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input); 171 void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input); 172 172 Penta* GetUpperElement(void); 173 173 void InputExtrude(int enum_type,bool only_if_collapsed);
Note:
See TracChangeset
for help on using the changeset viewer.