Index: /issm/trunk/src/c/Container/Inputs.cpp =================================================================== --- /issm/trunk/src/c/Container/Inputs.cpp (revision 4839) +++ /issm/trunk/src/c/Container/Inputs.cpp (revision 4840) @@ -292,162 +292,4 @@ } /*}}}*/ -/*FUNCTION Inputs::GetStrainRate2d{{{1*/ -void Inputs::GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){ - /*Compute the 2d Strain Rate (3 components): - * - * epsilon=[exx eyy exy] - */ - - vector::iterator object; - int i; - Input* vxinput=NULL; - Input* vyinput=NULL; - double epsilonvx[3]; - double epsilonvy[3]; - bool foundvx=false; - bool foundvy=false; - - /*Go through inputs and find data for vxenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vxinput=(Input*)(*object); - if (vxinput->EnumType()==vxenum){ - foundvx=true; - break; - } - } - /*Go through inputs and find data for vyenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vyinput=(Input*)(*object); - if (vyinput->EnumType()==vyenum){ - foundvy=true; - break; - } - } - - /*Check that both inputs have been found*/ - if (!foundvx || !foundvy){ - ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum)); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vxinput->GetVxStrainRate2d(epsilonvx,xyz_list,gauss); - vyinput->GetVyStrainRate2d(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - -} -/*}}}*/ -/*FUNCTION Inputs::GetStrainRate3d{{{1*/ -void Inputs::GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum){ - /*Compute the 3d Strain Rate (6 components): - * - * epsilon=[exx eyy ezz exy exz eyz] - */ - - int i; - vector::iterator object; - Input* vxinput=NULL; - Input* vyinput=NULL; - Input* vzinput=NULL; - bool foundvx=false; - bool foundvy=false; - bool foundvz=false; - double epsilonvx[6]; - double epsilonvy[6]; - double epsilonvz[6]; - - /*Go through inputs and find data for vxenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vxinput=(Input*)(*object); - if (vxinput->EnumType()==vxenum){ - foundvx=true; - break; - } - } - /*Go through inputs and find data for vyenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vyinput=(Input*)(*object); - if (vyinput->EnumType()==vyenum){ - foundvy=true; - break; - } - } - /*Go through inputs and find data for vzenum, not for Pattyn*/ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vzinput=(Input*)(*object); - if (vzinput->EnumType()==vzenum){ - foundvz=true; - break; - } - } - - /*Check that all inputs have been found*/ - if (!foundvx || !foundvy || !foundvz){ - 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)); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vxinput->GetVxStrainRate3d(epsilonvx,xyz_list,gauss); - vyinput->GetVyStrainRate3d(epsilonvy,xyz_list,gauss); - vzinput->GetVzStrainRate3d(epsilonvz,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<6;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]+epsilonvz[i]; - -} -/*}}}*/ -/*FUNCTION Inputs::GetStrainRate3dPattyn{{{1*/ -void Inputs::GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum){ - /*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; - vector::iterator object; - Input* vxinput=NULL; - Input* vyinput=NULL; - bool foundvx=false; - bool foundvy=false; - double epsilonvx[5]; - double epsilonvy[5]; - - /*Go through inputs and find data for vxenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vxinput=(Input*)(*object); - if (vxinput->EnumType()==vxenum){ - foundvx=true; - break; - } - } - /*Go through inputs and find data for vyenum: */ - for ( object=objects.begin() ; object < objects.end(); object++ ){ - vyinput=(Input*)(*object); - if (vyinput->EnumType()==vyenum){ - foundvy=true; - break; - } - } - - /*Check that all inputs have been found*/ - if (!foundvx || !foundvy){ - ISSMERROR("Could not find input with enum %i (%s) or enum %i (%s)",vxenum,EnumAsString(vxenum),vyenum,EnumAsString(vyenum)); - } - - /*Get strain rate assuming that epsilon has been allocated*/ - vxinput->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); - vyinput->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); - - /*Sum all contributions*/ - for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; - -} -/*}}}*/ /*FUNCTION Inputs::AddInput{{{1*/ int Inputs::AddInput(Input* in_input){ Index: /issm/trunk/src/c/Container/Inputs.h =================================================================== --- /issm/trunk/src/c/Container/Inputs.h (revision 4839) +++ /issm/trunk/src/c/Container/Inputs.h (revision 4840) @@ -45,7 +45,4 @@ void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); - void GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum); - void GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum,int vzenum); - void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, double* gauss, int vxenum, int vyenum); void ChangeEnum(int enumtype,int new_enumtype); @@ -55,5 +52,3 @@ }; - - #endif //ifndef _INPUTS_H_ Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4839) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 4840) @@ -490,8 +490,11 @@ /*flags: */ bool onbed; + Input* pressure_input=NULL; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; /*parameters: */ double stokesreconditioning; - /*retrive parameters: */ @@ -528,4 +531,10 @@ GaussTria(&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights,2); + /*Retrieve all inputs we will be needing: */ + pressure_input=inputs->GetInput(PressureEnum); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + /* Start looping on the number of gaussian points: */ for (ig=0; igGetStrainRate3d(&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]); - inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); + pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); /*Compute Stress*/ @@ -587,4 +596,5 @@ /*inputs: */ bool onwater; + Input* surface_input=NULL; /*retrieve inputs :*/ @@ -603,6 +613,6 @@ GetDofList1(&doflist[0]); - /*recover value of surface at grids: */ - inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); + /*Retrieve all inputs we will be needing: */ + surface_input->GetParameterValues(&surface[0],&gauss[0][0],6); /*pressure is lithostatic: */ @@ -3055,6 +3065,4 @@ bool onbed; bool shelf; - - /*inputs: */ Input* vx_input=NULL; Input* vy_input=NULL; Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4839) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4840) @@ -930,4 +930,12 @@ double cm_maxdmp_slope; + /*inputs: */ + Input* thickness_input=NULL; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* adjointx_input=NULL; + Input* adjointy_input=NULL; + Input* rheologyb_input=NULL; + /*retrieve some parameters: */ this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); @@ -943,4 +951,12 @@ /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); + + /*Retrieve all inputs we will be needing: */ + thickness_input=inputs->GetInput(ThicknessEnum); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + adjointx_input=inputs->GetInput(AdjointxEnum); + adjointy_input=inputs->GetInput(AdjointyEnum); + rheologyb_input=inputs->GetInput(RheologyBEnum); /* Start looping on the number of gaussian points: */ @@ -953,8 +969,8 @@ /*Get thickness: */ - inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum); + thickness_input->GetParameterValue(&thickness, gauss_l1l2l3); /*Get strain rate, if velocity has been supplied: */ - inputs->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum); + this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,vx_input,vy_input); /*Get viscosity complement: */ @@ -962,8 +978,8 @@ /*Get dvx, dvy, dadjx and dadjx: */ - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyEnum); - inputs->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],&gauss_l1l2l3[0],AdjointxEnum); - inputs->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],&gauss_l1l2l3[0],AdjointyEnum); + vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + adjointx_input->GetParameterDerivativeValue(&dadjx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + adjointy_input->GetParameterDerivativeValue(&dadjy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); /* Get Jacobian determinant: */ @@ -977,6 +993,6 @@ /*Get B derivative: dB/dx */ - inputs->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],&gauss_l1l2l3[0],RheologyBEnum); - inputs->GetParameterValue(&B_gauss, gauss_l1l2l3,RheologyBEnum); + rheologyb_input->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + rheologyb_input->GetParameterValue(&B_gauss, gauss_l1l2l3); /*Build gradje_g_gaussian vector (actually -dJ/dB): */ @@ -1063,4 +1079,9 @@ /*inputs: */ bool shelf; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* adjointx_input=NULL; + Input* adjointy_input=NULL; + Input* dragcoefficient_input=NULL; /*parameters: */ @@ -1104,4 +1125,11 @@ /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4); + + /*Retrieve all inputs we will be needing: */ + adjointx_input=inputs->GetInput(AdjointxEnum); + adjointy_input=inputs->GetInput(AdjointyEnum); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + dragcoefficient_input=inputs->GetInput(DragCoefficientEnum); /* Start looping on the number of gaussian points: */ @@ -1119,9 +1147,9 @@ /*Recover alpha_complement and k: */ - inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); + dragcoefficient_input->GetParameterValue(&drag, gauss_l1l2l3); /*recover lambda and mu: */ - inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); - inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); + adjointx_input->GetParameterValue(&lambda, gauss_l1l2l3); + adjointy_input->GetParameterValue(&mu, gauss_l1l2l3); /*recover vx and vy: */ @@ -1139,5 +1167,5 @@ /*Get k derivative: dk/dx */ - inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); + dragcoefficient_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0]); /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */