Changeset 4840
- Timestamp:
- 07/27/10 16:04:10 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Inputs.cpp
r4698 r4840 292 292 } 293 293 /*}}}*/ 294 /*FUNCTION Inputs::GetStrainRate2d{{{1*/295 void Inputs::GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){296 /*Compute the 2d Strain Rate (3 components):297 *298 * epsilon=[exx eyy exy]299 */300 301 vector<Object*>::iterator object;302 int i;303 Input* vxinput=NULL;304 Input* vyinput=NULL;305 double epsilonvx[3];306 double epsilonvy[3];307 bool foundvx=false;308 bool foundvy=false;309 310 /*Go through inputs and find data for vxenum: */311 for ( object=objects.begin() ; object < objects.end(); object++ ){312 vxinput=(Input*)(*object);313 if (vxinput->EnumType()==vxenum){314 foundvx=true;315 break;316 }317 }318 /*Go through inputs and find data for vyenum: */319 for ( object=objects.begin() ; object < objects.end(); object++ ){320 vyinput=(Input*)(*object);321 if (vyinput->EnumType()==vyenum){322 foundvy=true;323 break;324 }325 }326 327 /*Check that both inputs have been found*/328 if (!foundvx || !foundvy){329 ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum));330 }331 332 /*Get strain rate assuming that epsilon has been allocated*/333 vxinput->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);334 vyinput->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);335 336 /*Sum all contributions*/337 for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];338 339 }340 /*}}}*/341 /*FUNCTION Inputs::GetStrainRate3d{{{1*/342 void Inputs::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum){343 /*Compute the 3d Strain Rate (6 components):344 *345 * epsilon=[exx eyy ezz exy exz eyz]346 */347 348 int i;349 vector<Object*>::iterator object;350 Input* vxinput=NULL;351 Input* vyinput=NULL;352 Input* vzinput=NULL;353 bool foundvx=false;354 bool foundvy=false;355 bool foundvz=false;356 double epsilonvx[6];357 double epsilonvy[6];358 double epsilonvz[6];359 360 /*Go through inputs and find data for vxenum: */361 for ( object=objects.begin() ; object < objects.end(); object++ ){362 vxinput=(Input*)(*object);363 if (vxinput->EnumType()==vxenum){364 foundvx=true;365 break;366 }367 }368 /*Go through inputs and find data for vyenum: */369 for ( object=objects.begin() ; object < objects.end(); object++ ){370 vyinput=(Input*)(*object);371 if (vyinput->EnumType()==vyenum){372 foundvy=true;373 break;374 }375 }376 /*Go through inputs and find data for vzenum, not for Pattyn*/377 for ( object=objects.begin() ; object < objects.end(); object++ ){378 vzinput=(Input*)(*object);379 if (vzinput->EnumType()==vzenum){380 foundvz=true;381 break;382 }383 }384 385 /*Check that all inputs have been found*/386 if (!foundvx || !foundvy || !foundvz){387 ISSMERROR("Could not find input with enum %i (%s), enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum),vzenum,EnumAsString(vzenum));388 }389 390 /*Get strain rate assuming that epsilon has been allocated*/391 vxinput->GetVxStrainRate3d(epsilonvx,xyz_list,gauss);392 vyinput->GetVyStrainRate3d(epsilonvy,xyz_list,gauss);393 vzinput->GetVzStrainRate3d(epsilonvz,xyz_list,gauss);394 395 /*Sum all contributions*/396 for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i];397 398 }399 /*}}}*/400 /*FUNCTION Inputs::GetStrainRate3dPattyn{{{1*/401 void Inputs::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){402 /*Compute the 3d Blatter/PattynStrain Rate (5 components):403 *404 * epsilon=[exx eyy exy exz eyz]405 *406 * with exz=1/2 du/dz407 * eyz=1/2 dv/dz408 *409 * the contribution of vz is neglected410 */411 412 int i;413 vector<Object*>::iterator object;414 Input* vxinput=NULL;415 Input* vyinput=NULL;416 bool foundvx=false;417 bool foundvy=false;418 double epsilonvx[5];419 double epsilonvy[5];420 421 /*Go through inputs and find data for vxenum: */422 for ( object=objects.begin() ; object < objects.end(); object++ ){423 vxinput=(Input*)(*object);424 if (vxinput->EnumType()==vxenum){425 foundvx=true;426 break;427 }428 }429 /*Go through inputs and find data for vyenum: */430 for ( object=objects.begin() ; object < objects.end(); object++ ){431 vyinput=(Input*)(*object);432 if (vyinput->EnumType()==vyenum){433 foundvy=true;434 break;435 }436 }437 438 /*Check that all inputs have been found*/439 if (!foundvx || !foundvy){440 ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum));441 }442 443 /*Get strain rate assuming that epsilon has been allocated*/444 vxinput->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);445 vyinput->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);446 447 /*Sum all contributions*/448 for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];449 450 }451 /*}}}*/452 294 /*FUNCTION Inputs::AddInput{{{1*/ 453 295 int Inputs::AddInput(Input* in_input){ -
issm/trunk/src/c/Container/Inputs.h
r4698 r4840 45 45 46 46 void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); 47 void GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum);48 void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum);49 void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum);50 47 51 48 void ChangeEnum(int enumtype,int new_enumtype); … … 55 52 }; 56 53 57 58 59 54 #endif //ifndef _INPUTS_H_ -
issm/trunk/src/c/objects/Elements/Penta.cpp
r4839 r4840 490 490 /*flags: */ 491 491 bool onbed; 492 Input* pressure_input=NULL; 493 Input* vx_input=NULL; 494 Input* vy_input=NULL; 495 Input* vz_input=NULL; 492 496 493 497 /*parameters: */ 494 498 double stokesreconditioning; 495 496 499 497 500 /*retrive parameters: */ … … 528 531 GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2); 529 532 533 /*Retrieve all inputs we will be needing: */ 534 pressure_input=inputs->GetInput(PressureEnum); 535 vx_input=inputs->GetInput(VxEnum); 536 vy_input=inputs->GetInput(VyEnum); 537 vz_input=inputs->GetInput(VzEnum); 538 530 539 /* Start looping on the number of gaussian points: */ 531 540 for (ig=0; ig<num_gauss; ig++){ … … 539 548 540 549 /*Compute strain rate viscosity and pressure: */ 541 inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum);550 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); 542 551 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); 543 inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum);552 pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); 544 553 545 554 /*Compute Stress*/ … … 587 596 /*inputs: */ 588 597 bool onwater; 598 Input* surface_input=NULL; 589 599 590 600 /*retrieve inputs :*/ … … 603 613 GetDofList1(&doflist[0]); 604 614 605 /* recover value of surface at grids: */606 inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);615 /*Retrieve all inputs we will be needing: */ 616 surface_input->GetParameterValues(&surface[0],&gauss[0][0],6); 607 617 608 618 /*pressure is lithostatic: */ … … 3055 3065 bool onbed; 3056 3066 bool shelf; 3057 3058 /*inputs: */3059 3067 Input* vx_input=NULL; 3060 3068 Input* vy_input=NULL; -
issm/trunk/src/c/objects/Elements/Tria.cpp
r4839 r4840 930 930 double cm_maxdmp_slope; 931 931 932 /*inputs: */ 933 Input* thickness_input=NULL; 934 Input* vx_input=NULL; 935 Input* vy_input=NULL; 936 Input* adjointx_input=NULL; 937 Input* adjointy_input=NULL; 938 Input* rheologyb_input=NULL; 939 932 940 /*retrieve some parameters: */ 933 941 this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); … … 943 951 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 944 952 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); 953 954 /*Retrieve all inputs we will be needing: */ 955 thickness_input=inputs->GetInput(ThicknessEnum); 956 vx_input=inputs->GetInput(VxEnum); 957 vy_input=inputs->GetInput(VyEnum); 958 adjointx_input=inputs->GetInput(AdjointxEnum); 959 adjointy_input=inputs->GetInput(AdjointyEnum); 960 rheologyb_input=inputs->GetInput(RheologyBEnum); 945 961 946 962 /* Start looping on the number of gaussian points: */ … … 953 969 954 970 /*Get thickness: */ 955 inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum);971 thickness_input->GetParameterValue(&thickness, gauss_l1l2l3); 956 972 957 973 /*Get strain rate, if velocity has been supplied: */ 958 inputs->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum);974 this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input); 959 975 960 976 /*Get viscosity complement: */ … … 962 978 963 979 /*Get dvx, dvy, dadjx and dadjx: */ 964 inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxEnum);965 inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyEnum);966 inputs->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],&gauss_l1l2l3[0],AdjointxEnum);967 inputs->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],&gauss_l1l2l3[0],AdjointyEnum);980 vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 981 vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 982 adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 983 adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 968 984 969 985 /* Get Jacobian determinant: */ … … 977 993 978 994 /*Get B derivative: dB/dx */ 979 inputs->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],&gauss_l1l2l3[0],RheologyBEnum);980 inputs->GetParameterValue(&B_gauss, gauss_l1l2l3,RheologyBEnum);995 rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 996 rheologyb_input->GetParameterValue(&B_gauss, gauss_l1l2l3); 981 997 982 998 /*Build gradje_g_gaussian vector (actually -dJ/dB): */ … … 1063 1079 /*inputs: */ 1064 1080 bool shelf; 1081 Input* vx_input=NULL; 1082 Input* vy_input=NULL; 1083 Input* adjointx_input=NULL; 1084 Input* adjointy_input=NULL; 1085 Input* dragcoefficient_input=NULL; 1065 1086 1066 1087 /*parameters: */ … … 1104 1125 /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 1105 1126 GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); 1127 1128 /*Retrieve all inputs we will be needing: */ 1129 adjointx_input=inputs->GetInput(AdjointxEnum); 1130 adjointy_input=inputs->GetInput(AdjointyEnum); 1131 vx_input=inputs->GetInput(VxEnum); 1132 vy_input=inputs->GetInput(VyEnum); 1133 dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); 1106 1134 1107 1135 /* Start looping on the number of gaussian points: */ … … 1119 1147 1120 1148 /*Recover alpha_complement and k: */ 1121 inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum);1149 dragcoefficient_input->GetParameterValue(&drag, gauss_l1l2l3); 1122 1150 1123 1151 /*recover lambda and mu: */ 1124 inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum);1125 inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum);1152 adjointx_input->GetParameterValue(&lambda, gauss_l1l2l3); 1153 adjointy_input->GetParameterValue(&mu, gauss_l1l2l3); 1126 1154 1127 1155 /*recover vx and vy: */ … … 1139 1167 1140 1168 /*Get k derivative: dk/dx */ 1141 inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum);1169 dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0]); 1142 1170 1143 1171 /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
Note:
See TracChangeset
for help on using the changeset viewer.