Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4837) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4838) @@ -3024,5 +3024,4 @@ int numberofdofspernode; - /* 3d gaussian points: */ int num_gauss,ig; @@ -3097,4 +3096,10 @@ bool shelf; + /*inputs: */ + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vxold_input=NULL; + Input* vyold_input=NULL; + /*retrieve inputs :*/ inputs->GetParameterValue(&onwater,ElementOnWaterEnum); @@ -3112,5 +3117,4 @@ bedrock, in which case we spawn a tria element using the 3 first grids, and use it to build the stiffness matrix. */ - if ((collapse==1) && (onbed==0)){ @@ -3136,5 +3140,4 @@ GetDofList(&doflist[0],&numberofdofspernode); - /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore get tria gaussian points as well as segment gaussian points. For tria gaussian @@ -3148,4 +3151,10 @@ 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); + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vxold_input=inputs->GetInput(VxOldEnum); + vyold_input=inputs->GetInput(VyOldEnum); + /* Start looping on the number of gaussian points: */ for (ig1=0; ig1GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum); - inputs->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum); + this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input); + this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss_coord,vxold_input,vyold_input); /*Get viscosity: */ @@ -3344,4 +3353,9 @@ bool isstokes; + /*inputs: */ + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; + /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -3379,4 +3393,9 @@ /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 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); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); /* Start looping on the number of gaussian points: */ @@ -3393,5 +3412,5 @@ /*Compute strain rate: */ - inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); /*Get viscosity: */ @@ -3462,5 +3481,5 @@ /*Compute strain rate: */ - inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); /*Get viscosity at last iteration: */ @@ -4094,4 +4113,6 @@ bool collapse; bool onbed; + Input* surface_input=NULL; + Input* thickness_input=NULL; /*retrieve inputs :*/ @@ -4135,4 +4156,8 @@ 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); + /*Retrieve all inputs we will be needing: */ + thickness_input=inputs->GetInput(ThicknessEnum); + surface_input=inputs->GetInput(SurfaceEnum); + /* Start looping on the number of gaussian points: */ for (ig1=0; ig1GetParameterValue(&thickness, gauss_coord,ThicknessEnum); + thickness_input->GetParameterValue(&thickness, gauss_coord); /*Compute slope at gaussian point: */ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum); + surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord); /* Get Jacobian determinant: */ @@ -4293,4 +4318,8 @@ bool shelf; bool isstokes; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; + Input* bed_input=NULL; /*retrieve inputs :*/ @@ -4326,4 +4355,10 @@ /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ 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); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + bed_input=inputs->GetInput(BedEnum); /* Start looping on the number of gaussian points: */ @@ -4340,5 +4375,5 @@ /*Compute strain rate and viscosity: */ - inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -4421,5 +4456,5 @@ /* Get bed at gaussian point */ - inputs->GetParameterValue(&bed, gauss_coord,BedEnum); + bed_input->GetParameterValue(&bed, gauss_coord); /*Get L matrix : */ @@ -4514,4 +4549,6 @@ bool onwater; bool onbed; + Input* vx_input=NULL; + Input* vy_input=NULL; /*retrieve inputs :*/ @@ -4541,4 +4578,8 @@ 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); + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + /* Start looping on the number of gaussian points: */ for (ig1=0; ig1GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum); - inputs->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord,VyEnum); + vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord); + vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss_coord); dudx=du[0]; dvdy=dv[1]; @@ -4725,4 +4766,7 @@ bool onbed; bool shelf; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; /*retrieve inputs :*/ @@ -4757,4 +4801,9 @@ /*Get gaussian points: */ 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); + + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); /* Start looping on the number of gaussian points: */ @@ -4771,5 +4820,5 @@ /*Compute strain rate and viscosity: */ - inputs->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss_coord,vx_input,vy_input,vz_input); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -6014,4 +6063,63 @@ } /*}}}*/ +/*FUNCTION Penta::GetStrainRate3dPattyn{{{1*/ +void Penta::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input){ + /*Compute the 3d Blatter/PattynStrain Rate (5 components): + * + * epsilon=[exx eyy exy exz eyz] + * + * with exz=1/2 du/dz + * eyz=1/2 dv/dz + * + * the contribution of vz is neglected + */ + + int i; + + double epsilonvx[5]; + double epsilonvy[5]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input){ + ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + +} +/*}}}*/ +/*FUNCTION Penta::GetStrainRate3d{{{1*/ +void Penta::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input){ + /*Compute the 3d Strain Rate (6 components): + * + * epsilon=[exx eyy ezz exy exz eyz] + */ + + int i; + + double epsilonvx[6]; + double epsilonvy[6]; + double epsilonvz[6]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !vz_input){ + ISSMERROR("Input missing. Here are the input pointers we have for vx: %p, vy: %p, vz: %p\n",vx_input,vy_input,vz_input); + } + + /*Get strain rate assuming that epsilon has been allocated*/ + vx_input->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); + vz_input->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); + + /*Sum all contributions*/ + for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; + +} +/*}}}*/ /*FUNCTION Penta::GradjB {{{1*/ void Penta::GradjB(Vec gradient){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 4837) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 4838) @@ -168,6 +168,6 @@ void GetSolutionFromInputsDiagnosticVert(Vec solutiong); void GetSolutionFromInputsThermal(Vec solutiong); - void GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord); - void GetStrainRateStokes(double* epsilon, double* velocity, double* xyz_list, double* gauss_coord); + void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input); + void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input); Penta* GetUpperElement(void); void InputExtrude(int enum_type,bool only_if_collapsed);