Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4835) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 4836) @@ -2534,5 +2534,7 @@ int found=0; - /*parameters: */ + /*inputs: */ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; bool artdiff; @@ -2543,4 +2545,8 @@ /*retrieve some parameters: */ this->parameters->FindParam(&artdiff,ArtDiffEnum); + + /*Retrieve all inputs we will be needed*/ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); //Create Artificial diffusivity once for all if requested @@ -2551,6 +2557,6 @@ //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + vxaverage_input->GetParameterAverage(&v_gauss[0]); + vyaverage_input->GetParameterAverage(&v_gauss[1]); K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); @@ -2577,9 +2583,9 @@ //Get vx, vy and their derivatives at gauss point - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + vxaverage_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); + vyaverage_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); + + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); dvxdx=dvx[0]; @@ -2678,4 +2684,8 @@ int found; + /*inputs: */ + Input* vx_input=NULL; + Input* vy_input=NULL; + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -2684,4 +2694,8 @@ /* 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, 2); + + /*Retrieve all inputs we will be needed*/ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); /* Start looping on the number of gaussian points: */ @@ -2702,6 +2716,6 @@ //Get vx, vy and their derivatives at gauss point - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); + vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); DL_scalar=-gauss_weight*Jdettria; @@ -2783,5 +2797,7 @@ int found=0; - /*parameters: */ + /*inputs: */ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; bool artdiff; @@ -2813,4 +2829,8 @@ ny=ny/norm; + /*Retrieve all inputs we will be needed*/ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); + //Create Artificial diffusivity once for all if requested if(artdiff){ @@ -2820,6 +2840,6 @@ //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + vxaverage_input->GetParameterAverage(&v_gauss[0]); + vyaverage_input->GetParameterAverage(&v_gauss[1]); K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! @@ -2846,9 +2866,9 @@ //Get vx, vy and their derivatives at gauss point - inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]); + vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]); + + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); dvxdx=dvx[0]; @@ -3089,5 +3109,4 @@ double gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED - double MAXSLOPE=.06; // 6 % double MOUNTAINKEXPONENT=10; @@ -3097,4 +3116,7 @@ int drag_type; Input* surface_input=NULL; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* vz_input=NULL; /*retrive parameters: */ @@ -3105,4 +3127,7 @@ inputs->GetParameterValue(&drag_type,DragTypeEnum); surface_input=inputs->GetInput(SurfaceEnum); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); /* Get node coordinates and dof list: */ @@ -3460,5 +3485,7 @@ int found; - /*parameters: */ + /*inputs: */ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; double dt; bool artdiff; @@ -3471,4 +3498,8 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); + + /*Retrieve all inputs we will be needing: */ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); //Create Artificial diffusivity once for all if requested @@ -3479,6 +3510,6 @@ //Build K matrix (artificial diffusivity matrix) - inputs->GetParameterAverage(&v_gauss[0],VxAverageEnum); - inputs->GetParameterAverage(&v_gauss[1],VyAverageEnum); + vxaverage_input->GetParameterAverage(&v_gauss[0]); + vyaverage_input->GetParameterAverage(&v_gauss[1]); K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); @@ -3517,9 +3548,9 @@ //Get vx, vy and their derivatives at gauss point - inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); - - inputs->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0],VyAverageEnum); + vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]); + vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]); + + vxaverage_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],&gauss_l1l2l3[0]); + vyaverage_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],&gauss_l1l2l3[0]); dvxdx=dvx[0]; @@ -3623,5 +3654,7 @@ int found; - /*parameters: */ + /*inputs: */ + Input* vxaverage_input=NULL; + Input* vyaverage_input=NULL; double dt; @@ -3635,4 +3668,8 @@ /* 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, 2); + + /*Retrieve all inputs we will be needed*/ + vxaverage_input=inputs->GetInput(VxAverageEnum); + vyaverage_input=inputs->GetInput(VyAverageEnum); /* Start looping on the number of gaussian points: */ @@ -3665,6 +3702,6 @@ //Get vx, vy and their derivatives at gauss point - inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); - inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); + vxaverage_input->GetParameterValue(&vx,&gauss_l1l2l3[0]); + vyaverage_input->GetParameterValue(&vy,&gauss_l1l2l3[0]); DL_scalar=-dt*gauss_weight*Jdettria; @@ -3826,5 +3863,4 @@ rho_ice=matpar->GetRhoIce(); heatcapacity=matpar->GetHeatCapacity(); - GaussTria (&num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); @@ -3904,4 +3940,8 @@ double melting_g; + /*inputs: */ + Input* accumulation_input=NULL; + Input* melting_input=NULL; + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -3911,4 +3951,8 @@ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + /*retrieve inputs :*/ + accumulation_input=inputs->GetInput(AccumulationRateEnum); + melting_input=inputs->GetInput(MeltingRateEnum); + /* Start looping on the number of gaussian points: */ for (ig=0; igGetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); /* Add value into pe_g: */ @@ -3947,5 +3991,4 @@ /*FUNCTION Tria::CreatePVectorBalancedthickness_DG {{{1*/ void Tria::CreatePVectorBalancedthickness_DG(Vec pg){ - /* local declarations */ @@ -3979,4 +4022,9 @@ double dhdt_g; + /*inputs: */ + Input* accumulation_input=NULL; + Input* melting_input=NULL; + Input* dhdt_input=NULL; + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -3985,4 +4033,9 @@ /* 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, 2); + + /*Retrieve all inputs we will be needing: */ + accumulation_input=inputs->GetInput(AccumulationRateEnum); + melting_input=inputs->GetInput(MeltingRateEnum); + dhdt_input=inputs->GetInput(DhDtEnum); /* Start looping on the number of gaussian points: */ @@ -4001,7 +4054,7 @@ /* Get accumulation, melting and thickness at gauss point */ - inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0],DhDtEnum); + accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); + dhdt_input->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0]); /* Add value into pe_g: */ @@ -4054,4 +4107,8 @@ double melting_g; + /*inputs: */ + Input* accumulation_input=NULL; + Input* melting_input=NULL; + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -4060,4 +4117,8 @@ /* 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, 2); + + /*Retrieve all inputs we will be needing: */ + accumulation_input=inputs->GetInput(AccumulationRateEnum); + melting_input=inputs->GetInput(MeltingRateEnum); /* Start looping on the number of gaussian points: */ @@ -4076,6 +4137,6 @@ /* Get accumulation, melting at gauss point */ - inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); /* Add value into pe_g: */ @@ -4136,4 +4197,10 @@ double dbdx,dbdy; + /*inputs: */ + Input* melting_input=NULL; + Input* vx_input=NULL; + Input* vy_input=NULL; + Input* bed_input=NULL; + /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -4142,4 +4209,10 @@ /* 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, 2); + + /*Retrieve all inputs we will be needing: */ + melting_input=inputs->GetInput(MeltingRateEnum); + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + bed_input=inputs->GetInput(BedEnum); /*For icesheets: */ @@ -4154,12 +4227,12 @@ /*Get melting at gaussian point: */ - inputs->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0],MeltingRateEnum); + melting_input->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0]); /*Get velocity at gaussian point: */ - inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); - inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + vx_input->GetParameterValue(&vx, &gauss_l1l2l3[0]); + vy_input->GetParameterValue(&vy, &gauss_l1l2l3[0]); /*Get bed slope: */ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); + bed_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]); dbdx=slope[0]; dbdy=slope[1]; @@ -4170,5 +4243,4 @@ //Get L matrix if viscous basal drag present: GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,NDOF1); - /*Build gaussian vector: */ @@ -4380,6 +4452,9 @@ double thickness_g; - /*parameters: */ + /*inputs: */ double dt; + Input* accumulation_input=NULL; + Input* melting_input=NULL; + Input* thickness_input=NULL; /*retrieve some parameters: */ @@ -4392,4 +4467,9 @@ /* 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, 2); + + /*Retrieve all inputs we will be needing: */ + accumulation_input=inputs->GetInput(AccumulationRateEnum); + melting_input=inputs->GetInput(MeltingRateEnum); + thickness_input=inputs->GetInput(ThicknessEnum); /* Start looping on the number of gaussian points: */ @@ -4408,7 +4488,7 @@ /* Get accumulation, melting and thickness at gauss point */ - inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); + accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); + thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]); /* Add value into pe_g: */ @@ -4461,4 +4541,7 @@ double thickness_g; double dt; + Input* accumulation_input=NULL; + Input* melting_input=NULL; + Input* thickness_input=NULL; /*retrieve some parameters: */ @@ -4471,4 +4554,9 @@ /* 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, 2); + + /*Retrieve all inputs we will be needing: */ + accumulation_input=inputs->GetInput(AccumulationRateEnum); + melting_input=inputs->GetInput(MeltingRateEnum); + thickness_input=inputs->GetInput(ThicknessEnum); /* Start looping on the number of gaussian points: */ @@ -4487,7 +4575,7 @@ /* Get accumulation, melting and thickness at gauss point */ - inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); - inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); - inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); + accumulation_input->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0]); + melting_input->GetParameterValue(&melting_g, &gauss_l1l2l3[0]); + thickness_input->GetParameterValue(&thickness_g, &gauss_l1l2l3[0]); /* Add value into pe_g: */ @@ -4541,4 +4629,7 @@ int analysis_type; + /*inputs :*/ + Input* slope_input=NULL; + /*retrive parameters: */ parameters->FindParam(&analysis_type,AnalysisTypeEnum); @@ -4548,9 +4639,15 @@ GetDofList(&doflist[0],&numberofdofspernode); - /* Get gaussian points and weights: */ GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/ - + /*Retrieve all inputs we will be needing: */ + if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){ + slope_input=inputs->GetInput(SurfaceEnum); + } + if ( (analysis_type==BedSlopeXAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ + slope_input=inputs->GetInput(BedEnum); + } + /* Start looping on the number of gaussian points: */ for (ig=0; igGetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); - } - if ( (analysis_type==BedSlopeXAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){ - inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); - } + slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0]); /* Get Jacobian determinant: */ @@ -4639,4 +4731,7 @@ double t_pmp; double scalar_ocean; + + /*inputs: */ + Input* pressure_input=NULL; /* Get node coordinates and dof list: */ @@ -4660,4 +4755,7 @@ GaussTria (&num_area_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); + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ for (ig=0; igGetParameterValue(&pressure, &gauss_coord[0],PressureEnum); + pressure_input->GetParameterValue(&pressure, &gauss_coord[0]); t_pmp=meltingpoint-beta*pressure; @@ -4749,4 +4847,9 @@ parameters->FindParam(&analysis_type,AnalysisTypeEnum); + /*inputs: */ + Input* vx_input=NULL; //TO BE DELETED + Input* vy_input=NULL; //TO BE DELETED + Input* vz_input=NULL; //TO BE DELETED + Input* geothermalflux_input=NULL; /* Get node coordinates and dof list: */ @@ -4766,10 +4869,16 @@ friction=new Friction("3d",inputs,matpar,analysis_type); + /*Retrieve all inputs we will be needing: */ + vx_input=inputs->GetInput(VxEnum); + vy_input=inputs->GetInput(VyEnum); + vz_input=inputs->GetInput(VzEnum); + geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); + /*COMPUT alpha2_list and basalfriction_list (TO BE DELETED)*/ for(i=0;iGetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum); //TO BE DELETED } - inputs->GetParameterValues(&vx_list[0],&gauss[0][0],3,VxEnum); //TO BE DELETED - inputs->GetParameterValues(&vy_list[0],&gauss[0][0],3,VyEnum); //TO BE DELETED + vx_input->GetParameterValues(&vx_list[0],&gauss[0][0],3); //TO BE DELETED + vy_input->GetParameterValues(&vy_list[0],&gauss[0][0],3); //TO BE DELETED for(i=0;iGetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); + geothermalflux_input->GetParameterValue(&geothermalflux_value, &gauss_coord[0]); /*Friction: */