Index: /issm/trunk/src/c/DataSet/Inputs.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3646) +++ /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3647) @@ -102,6 +102,6 @@ } /*}}}*/ -/*FUNCTION Inputs::GetParameterValueAtNode(double* pvalue, Node* node, int enum_type){{{1*/ -void Inputs::GetParameterValueAtNode(double* pvalue,Node* node,int enum_type){ +/*FUNCTION Inputs::GetParameterValue(double* pvalue, Node* node, int enum_type){{{1*/ +void Inputs::GetParameterValue(double* pvalue,Node* node,int enum_type){ /*given a node, instead of a gauss point, we want to recover a value: probably in an element!: */ @@ -124,5 +124,30 @@ /*Ok, we have an input if we made it here, request the input to return the values: */ - input->GetParameterValueAtNode(pvalue,node); + input->GetParameterValue(pvalue,node); +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValue(double* pvalue, Node* node1, Node* node2,int enum_type){{{1*/ +void Inputs::GetParameterValue(double* pvalue,Node* node1, Node* node2,double gauss_coord,int enum_type){ + + /*given a node, instead of a gauss point, we want to recover a value: probably in an element!: */ + + vector::iterator object; + Input* input=NULL; + + /*Go through inputs and check whether any input with the same name is already in: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + input=(Input*)(*object); + if (input->EnumType()==enum_type)break; + } + + if (!input){ + /*we could not find an input with the correct enum type. No defaults values were provided, + * error out: */ + ISSMERROR("%s%i"," could not find input with enum type ",enum_type); + } + + /*Ok, we have an input if we made it here, request the input to return the values: */ + input->GetParameterValue(pvalue,node1,node2,gauss_coord); } /*}}}*/ Index: /issm/trunk/src/c/DataSet/Inputs.h =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.h (revision 3646) +++ /issm/trunk/src/c/DataSet/Inputs.h (revision 3647) @@ -21,13 +21,17 @@ /*}}}*/ /*numerics: {{{1*/ + int AddInput(Input* in_input); + void GetParameterValue(bool* pvalue,int enum_type); void GetParameterValue(int* pvalue,int enum_type); void GetParameterValue(double* pvalue,int enum_type); - void GetParameterValueAtNode(double* pvalue,Node* node,int enum_type); + void GetParameterValue(double* pvalue,Node* node,int enum_type); + void GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_coord,int enum_type); void GetParameterValue(double* pvalue,double* gauss,int enum_type); void GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue); + void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type); void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues); - int AddInput(Input* in_input); + void GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type); void GetStrainRate(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3646) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3647) @@ -185,5 +185,7 @@ ConvergedEnum, SegmentOnIceShelfEnum, - FillEnum + FillEnum, + InternalEnum, + BoundaryEnum /*}}}*/ Index: /issm/trunk/src/c/objects/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3646) +++ /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3647) @@ -54,6 +54,5 @@ this->parameters=numericalflux_parameters; } -/*}}}*/ -/*FUNCTION Numericalflux::Numericalflux(int id, int i, IoModel* iomodel) {{{1*/ +/*}}}*//*FUNCTION Numericalflux::Numericalflux(int id, int i, IoModel* iomodel) {{{1*/ Numericalflux::Numericalflux(int numericalflux_id,int i, IoModel* iomodel){ @@ -148,4 +147,5 @@ } /*}}}*/ + /*FUNCTION Numericalflux::~Numericalflux(){{{1*/ Numericalflux::~Numericalflux(){ @@ -288,11 +288,16 @@ void Numericalflux::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ - if (strcmp(type,"internal")==0){ - - CreateKMatrixInternal(Kgg,inputs,analysis_type,sub_analysis_type); - } - else if (strcmp(type,"boundary")==0){ - - CreateKMatrixBoundary(Kgg,inputs,analysis_type,sub_analysis_type); + int type; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + if (type==InternalEnum){ + + CreateKMatrixInternal(Kgg,analysis_type,sub_analysis_type); + } + else if (type==BoundaryEnum){ + + CreateKMatrixBoundary(Kgg,analysis_type,sub_analysis_type); } else ISSMERROR("type not supported yet"); @@ -332,24 +337,20 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; double vx,vy; double UdotN; double dt; - int dofs[1]={0}; int found; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vx_average in inputs!"); - found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vy_average in inputs!"); + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Tria** trias=NULL; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); + trias=(Tria**)helements.deliverp(); + + /*recover parameters: */ if (analysis_type==Prognostic2AnalysisEnum){ - if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + parameters->FindParam(&dt,"dt"); } else if (analysis_type==Balancedthickness2AnalysisEnum){ @@ -381,6 +382,7 @@ //Get vx, vy and v.n - GetParameterValue(&vx, &vx_list[0],gauss_coord); - GetParameterValue(&vy, &vy_list[0],gauss_coord); + trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); + trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); + UdotN=vx*normal[0]+vy*normal[1]; if (fabs(UdotN)<1.0e-9 && analysis_type==Balancedthickness2AnalysisEnum) printf("Edge number %i has a flux very small (u.n = %g ), which could lead to unaccurate results\n",id,UdotN); @@ -447,6 +449,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; double vx,vy; double mean_vx; @@ -457,16 +457,15 @@ int found; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vx_average in inputs!"); - found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vy_average in inputs!"); + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Tria** trias=NULL; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); + trias=(Tria**)helements.deliverp(); + + /*recover parameters: */ if (analysis_type==Prognostic2AnalysisEnum){ - if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + parameters->FindParam(&dt,"dt"); } else if (analysis_type==Balancedthickness2AnalysisEnum){ @@ -477,5 +476,5 @@ ISSMERROR("analysis_type %i not supported yet"); } - + /* Get node coordinates, dof list and normal vector: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -484,6 +483,6 @@ /*Check wether it is an inflow or outflow BC*/ - mean_vx=(vx_list[0]+vx_list[1])/2.0; - mean_vy=(vy_list[0]+vy_list[1])/2.0; + trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum); + trias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN<=0){ @@ -507,6 +506,7 @@ //Get vx, vy and v.n - GetParameterValue(&vx, &vx_list[0],gauss_coord); - GetParameterValue(&vy, &vy_list[0],gauss_coord); + trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); + trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); + UdotN=vx*normal[0]+vy*normal[1]; @@ -539,11 +539,16 @@ void Numericalflux::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){ - if (strcmp(type,"internal")==0){ - - CreatePVectorInternal(pg,inputs,analysis_type,sub_analysis_type); - } - else if (strcmp(type,"boundary")==0){ - - CreatePVectorBoundary(pg,inputs,analysis_type,sub_analysis_type); + int type; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + if (type==InternalEnum){ + + CreatePVectorInternal(pg,analysis_type,sub_analysis_type); + } + else if (type==BoundaryEnum){ + + CreatePVectorBoundary(pg,analysis_type,sub_analysis_type); } else ISSMERROR("type not supported yet"); @@ -588,6 +593,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; double vx,vy; double mean_vx; @@ -599,16 +602,15 @@ int found; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("vx_average",&vx_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vx_average in inputs!"); - found=inputs->Recover("vy_average",&vy_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find vy_average in inputs!"); + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + Tria** trias=NULL; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); + trias=(Tria**)helements.deliverp(); + + /*recover parameters: */ if (analysis_type==Prognostic2AnalysisEnum){ - if(!inputs->Recover("dt",&dt))ISSMERROR(" could not find dt in inputs!"); + parameters->FindParam(&dt,"dt"); } else if (analysis_type==Balancedthickness2AnalysisEnum){ @@ -619,5 +621,5 @@ ISSMERROR("analysis_type %i not supported yet"); } - + /* Get node coordinates, dof list and normal vector: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -626,6 +628,6 @@ /*Check wether it is an inflow or outflow BC*/ - mean_vx=(vx_list[0]+vx_list[1])/2.0; - mean_vy=(vy_list[0]+vy_list[1])/2.0; + trias[0]->inputs->GetParameterValue(&mean_vx, nodes[0],nodes[1],.5,VxAverageEnum); + trias[0]->inputs->GetParameterValue(&mean_vy, nodes[0],nodes[1],.5,VyAverageEnum); UdotN=mean_vx*normal[0]+mean_vy*normal[1]; if (UdotN>0){ @@ -649,7 +651,8 @@ //Get vx, vy and v.n - GetParameterValue(&vx, &vx_list[0],gauss_coord); - GetParameterValue(&vy, &vy_list[0],gauss_coord); - GetParameterValue(&thickness,&h[0], gauss_coord); + trias[0]->inputs->GetParameterValue(&vx, nodes[0],nodes[1],gauss_coord,VxAverageEnum); + trias[0]->inputs->GetParameterValue(&vy, nodes[0],nodes[1],gauss_coord,VyAverageEnum); + trias[0]->inputs->GetParameterValue(&thickness, nodes[0],nodes[1],gauss_coord,ThicknessEnum); + UdotN=vx*normal[0]+vy*normal[1]; @@ -701,6 +704,16 @@ int doflist_per_node[MAXDOFSPERNODE]; int numberofdofspernode; + + /*dynamic objects pointed to by hooks: */ + Node** nodes=NULL; + int type; + + /*recover objects from hooks: */ + nodes=(Node**)hnodes.deliverp(); - if (strcmp(type,"internal")==0){ + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + if (type==InternalEnum){ for(i=0;i<4;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); @@ -710,5 +723,5 @@ } } - else if (strcmp(type,"boundary")==0){ + else if (type==BoundaryEnum){ for(i=0;i<2;i++){ nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode); Index: /issm/trunk/src/c/objects/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.h (revision 3646) +++ /issm/trunk/src/c/objects/Numericalflux.h (revision 3647) @@ -28,4 +28,5 @@ Numericalflux(int numericalflux_id,int* numericalflux_node_ids,int numnodes, int* numericalflux_element_ids, int numelements); Numericalflux(int numericalflux_id,Hook* numericalflux_hnodes, Hook* numericalflux_helements, Parameters* parameters, Inputs* numericalflux_inputs); + Numericalflux::Numericalflux(int numericalflux_id,int i, IoModel* iomodel); ~Numericalflux(); /*}}}*/ Index: /issm/trunk/src/c/objects/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Pengrid.cpp (revision 3646) +++ /issm/trunk/src/c/objects/Pengrid.cpp (revision 3647) @@ -364,6 +364,6 @@ //First recover pressure and temperature values, using the element: */ - penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum); - penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum); + penta->inputs->GetParameterValue(&pressure,node,PressureEnum); + penta->inputs->GetParameterValue(&temperature,node,TemperatureEnum); //Recover our data: @@ -463,6 +463,6 @@ //recover slope: */ - penta->inputs->GetParameterValueAtNode(&slope[0],node,BedSlopexEnum); - penta->inputs->GetParameterValueAtNode(&slope[1],node,BedSlopeyEnum); + penta->inputs->GetParameterValue(&slope[0],node,BedSlopexEnum); + penta->inputs->GetParameterValue(&slope[1],node,BedSlopeyEnum); /*recover parameters: */ @@ -511,6 +511,6 @@ //First recover pressure and temperature values, using the element: */ - penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum); - penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum); + penta->inputs->GetParameterValue(&pressure,node,PressureEnum); + penta->inputs->GetParameterValue(&temperature,node,TemperatureEnum); /*recover parameters: */ @@ -622,6 +622,6 @@ //First recover pressure and temperature values, using the element: */ - penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum); - penta->inputs->GetParameterValueAtNode(&temperature,node,TemperatureEnum); + penta->inputs->GetParameterValue(&pressure,node,PressureEnum); + penta->inputs->GetParameterValue(&temperature,node,TemperatureEnum); inputs->GetParameterValue(&melting_offset,MeltingOffsetEnum); parameters->FindParam(&dt,"dt"); @@ -689,5 +689,5 @@ //First recover pressure and penalty_offset - penta->inputs->GetParameterValueAtNode(&pressure,node,PressureEnum); + penta->inputs->GetParameterValue(&pressure,node,PressureEnum); parameters->FindParam(&penalty_offset,"penalty_offset"); Index: /issm/trunk/src/c/objects/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Riftfront.cpp (revision 3646) +++ /issm/trunk/src/c/objects/Riftfront.cpp (revision 3647) @@ -338,8 +338,8 @@ /*First recover velocity: */ - tria1->inputs->GetParameterValueAtNode(&vx1,nodes[0],VxEnum); - tria2->inputs->GetParameterValueAtNode(&vx2,nodes[1],VxEnum); - tria1->inputs->GetParameterValueAtNode(&vy1,nodes[0],VyEnum); - tria2->inputs->GetParameterValueAtNode(&vy2,nodes[1],VyEnum); + tria1->inputs->GetParameterValue(&vx1,nodes[0],VxEnum); + tria2->inputs->GetParameterValue(&vx2,nodes[1],VxEnum); + tria1->inputs->GetParameterValue(&vy1,nodes[0],VyEnum); + tria2->inputs->GetParameterValue(&vy2,nodes[1],VyEnum); /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ @@ -478,8 +478,8 @@ /*recover velocity: */ - tria1->inputs->GetParameterValueAtNode(&vx1,nodes[0],VxEnum); - tria2->inputs->GetParameterValueAtNode(&vx2,nodes[1],VxEnum); - tria1->inputs->GetParameterValueAtNode(&vy1,nodes[0],VyEnum); - tria2->inputs->GetParameterValueAtNode(&vy2,nodes[1],VyEnum); + tria1->inputs->GetParameterValue(&vx1,nodes[0],VxEnum); + tria2->inputs->GetParameterValue(&vx2,nodes[1],VxEnum); + tria1->inputs->GetParameterValue(&vy1,nodes[0],VyEnum); + tria2->inputs->GetParameterValue(&vy2,nodes[1],VyEnum); /*Grid 1 faces grid2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ @@ -561,6 +561,6 @@ /*Recover thickness: */ - tria1->inputs->GetParameterValueAtNode(&h[0],nodes[0],ThicknessEnum); - tria2->inputs->GetParameterValueAtNode(&h[1],nodes[1],ThicknessEnum); + tria1->inputs->GetParameterValue(&h[0],nodes[0],ThicknessEnum); + tria2->inputs->GetParameterValue(&h[1],nodes[1],ThicknessEnum); if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); @@ -688,12 +688,12 @@ /*get thickness: */ - tria1->inputs->GetParameterValueAtNode(&h[0],nodes[0],ThicknessEnum); - tria2->inputs->GetParameterValueAtNode(&h[1],nodes[1],ThicknessEnum); + tria1->inputs->GetParameterValue(&h[0],nodes[0],ThicknessEnum); + tria2->inputs->GetParameterValue(&h[1],nodes[1],ThicknessEnum); if (h[0]!=h[1])ISSMERROR(" different thicknesses not supported for rift fronts"); thickness=h[0]; - tria1->inputs->GetParameterValueAtNode(&b[0],nodes[0],BedEnum); - tria2->inputs->GetParameterValueAtNode(&b[1],nodes[1],BedEnum); + tria1->inputs->GetParameterValue(&b[0],nodes[0],BedEnum); + tria2->inputs->GetParameterValue(&b[1],nodes[1],BedEnum); if (b[0]!=b[1])ISSMERROR(" different beds not supported for rift fronts"); @@ -779,8 +779,8 @@ /*First recover velocity: */ - tria1->inputs->GetParameterValueAtNode(&vx1,nodes[0],VxEnum); - tria2->inputs->GetParameterValueAtNode(&vx2,nodes[1],VxEnum); - tria1->inputs->GetParameterValueAtNode(&vy1,nodes[0],VyEnum); - tria2->inputs->GetParameterValueAtNode(&vy2,nodes[1],VyEnum); + tria1->inputs->GetParameterValue(&vx1,nodes[0],VxEnum); + tria2->inputs->GetParameterValue(&vx2,nodes[1],VxEnum); + tria1->inputs->GetParameterValue(&vy1,nodes[0],VyEnum); + tria2->inputs->GetParameterValue(&vy2,nodes[1],VyEnum); /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ @@ -829,8 +829,8 @@ /*First recover velocity: */ - tria1->inputs->GetParameterValueAtNode(&vx1,nodes[0],VxEnum); - tria2->inputs->GetParameterValueAtNode(&vx2,nodes[1],VxEnum); - tria1->inputs->GetParameterValueAtNode(&vy1,nodes[0],VyEnum); - tria2->inputs->GetParameterValueAtNode(&vy2,nodes[1],VyEnum); + tria1->inputs->GetParameterValue(&vx1,nodes[0],VxEnum); + tria2->inputs->GetParameterValue(&vx2,nodes[1],VxEnum); + tria1->inputs->GetParameterValue(&vy1,nodes[0],VyEnum); + tria2->inputs->GetParameterValue(&vy2,nodes[1],VyEnum); /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */ @@ -885,8 +885,8 @@ /*First recover velocity: */ - tria1->inputs->GetParameterValueAtNode(&vx1,nodes[0],VxEnum); - tria2->inputs->GetParameterValueAtNode(&vx2,nodes[1],VxEnum); - tria1->inputs->GetParameterValueAtNode(&vy1,nodes[0],VyEnum); - tria2->inputs->GetParameterValueAtNode(&vy2,nodes[1],VyEnum); + tria1->inputs->GetParameterValue(&vx1,nodes[0],VxEnum); + tria2->inputs->GetParameterValue(&vx2,nodes[1],VxEnum); + tria1->inputs->GetParameterValue(&vy1,nodes[0],VyEnum); + tria2->inputs->GetParameterValue(&vy2,nodes[1],VyEnum); /*Node 1 faces node 2, compute penetration of 2 into 1 (V2-V1).N (with N normal vector, and V velocity vector: */