Index: /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h =================================================================== --- /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h (revision 3611) +++ /issm/trunk/src/c/ComputeBasalStressx/ComputeBasalStressx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void ComputeBasalStressx( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _COMPUTEBASALSTRESSX_H */ Index: /issm/trunk/src/c/ComputePressurex/ComputePressurex.h =================================================================== --- /issm/trunk/src/c/ComputePressurex/ComputePressurex.h (revision 3611) +++ /issm/trunk/src/c/ComputePressurex/ComputePressurex.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void ComputePressurex( Vec* pp_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _COMPUTEPRESSUREX_H */ Index: /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h =================================================================== --- /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h (revision 3611) +++ /issm/trunk/src/c/ComputeStrainRatex/ComputeStrainRatex.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void ComputeStrainRatex(Vec* eps_g,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _COMPUTESTRAINRATEX_H */ Index: /issm/trunk/src/c/CostFunctionx/CostFunctionx.h =================================================================== --- /issm/trunk/src/c/CostFunctionx/CostFunctionx.h (revision 3611) +++ /issm/trunk/src/c/CostFunctionx/CostFunctionx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _MISFITX_H */ Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3611) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 3612) @@ -225,4 +225,10 @@ dataset->AddObject(tria); } + else if(enum_type==TriaVertexInputEnum){ + TriaVertexInput* triavertexinput=NULL; + triavertexinput=new TriaVertexInput(); + triavertexinput->Demarshall(&marshalled_dataset); + dataset->AddObject(triavertexinput); + } else if(enum_type==SingEnum){ Sing* sing=NULL; @@ -316,25 +322,4 @@ /*Specific methods*/ -/*FUNCTION DataSet::AddEinput{{{1*/ -int DataSet::AddEinput(Einput* in_einput){ - - /*First, go through dataset of inputs and check whether any input - * with the same name is already in. If so, erase the corresponding - * object before adding this new one: */ - vector::iterator object; - Einput* einput=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - einput=(Einput*)(*object); //assume this is an inputs dataset - - if (einput->EnumType()==in_einput->EnumType()){ - this->DeleteObject(einput); - break; - } - } - this->AddObject(in_einput); -} -/*}}}*/ /*FUNCTION DataSet::AddObject{{{1*/ int DataSet::AddObject(Object* object){ Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 3611) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 3612) @@ -12,5 +12,4 @@ #include #include "../toolkits/toolkits.h" -#include "../objects/Einput.h" #include "../objects/Object.h" @@ -47,5 +46,4 @@ int MarshallSize(); int AddObject(Object* object); - int AddEinput(Einput* in_einput); int DeleteObject(int id); int Size(); Index: /issm/trunk/src/c/DataSet/Inputs.cpp =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3612) +++ /issm/trunk/src/c/DataSet/Inputs.cpp (revision 3612) @@ -0,0 +1,273 @@ +/*!\file Inputs.c + * \brief: implementation of the Inputs object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./DataSet.h" +#include "./Inputs.h" + +/*Object constructors and destructor*/ +/*FUNCTION Inputs::Inputs(){{{1*/ +Inputs::Inputs(){ + return; +} +/*}}}*/ +/*FUNCTION Inputs::Inputs(int in_enum){{{1*/ +Inputs::Inputs(int in_enum): DataSet(in_enum){ +{ + //do nothing; + return; +} +/*}}}*/ +/*FUNCTION Inputs::~Inputs(){{{1*/ +Inputs::~Inputs(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type){{{1*/ +void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type){ + + 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->GetParameterValues(pvalue,gauss); + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValue(double* pvalue,double* gauss,int enum_type,double defaultvalue){{{1*/ +void Inputs::GetParameterValue(double* pvalue,double* gauss, int enum_type,double defaultvalue){ + + 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. Return the default value: */ + *pvalue=defaultvalue; + } + else{ + input->GetParameterValue(pvalue,gauss); + } +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type){{{1*/ +void Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type){ + + 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->GetParameterValues(values,gauss_pointers,numgauss); + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues){{{1*/ +void GetParameterValues(double* values,double* gauss_pointers, int numgauss,int enum_type,double* defaultvalues){ + + 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. Return the default values: */ + for(i=0;iGetParameterValues(values,gauss_pointers,numgauss); + } + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValue(bool* pvalue,int enum-type){{{1*/ +void GetParameterValue(bool* pvalue,int enum_type){ + + 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 value: */ + input->GetParameterValue(pvalue); + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValue(int* pvalue,int enum-type){{{1*/ +void GetParameterValue(int* pvalue,int enum_type){ + + 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 value: */ + input->GetParameterValue(pvalue); + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterValue(double* pvalue,int enum-type){{{1*/ +void GetParameterValue(double* pvalue,int enum_type){ + + 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 value: */ + input->GetParameterValue(pvalue); + +} +/*}}}*/ +/*FUNCTION Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){{{1*/ +void Inputs::GetParameterDerivativeValue(double* derivativevalues, double* xyz_list, double* gauss,int enum_type){ + + 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 value: */ + input->GetParameterValue(derivativevalues,xyz_list,gauss); +} +/*}}}*/ +/*FUNCTION Inputs::GetStrainRate(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum){{{1*/ +void Inputs::GetStrainRate(double* epsilon,double* xyz_list, double* gauss, int xenum, int yenum){ + + vector::iterator object; + Input* xinput=NULL; + Input* yinput=NULL; + + /*Go through inputs and find data for xenum: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + xinput=(Input*)(*object); + if (xinput->EnumType()==xenum)break; + } + /*Go through inputs and find data for yenum: */ + for ( object=objects.begin() ; object < objects.end(); object++ ){ + yinput=(Input*)(*object); + if (yinput->EnumType()==yenum)break; + } + + if (!xinput | !yinput){ + /*we could not find one input with the correct enum type. No defaults values were provided, + * error out: */ + ISSMERROR("%s%i%s%i\n"," could not find input with enum type ",xenum," or enum type ",yenum); + } + + /*Ok, we have the inputs, call bilinear operator: */ + xinput->GetStrainRate(epsilon,yinput,xyz_list,gauss); + +} +/*}}}*/ +/*FUNCTION Inputs::AddInput{{{1*/ +int Inputs::AddInput(Input* in_input){ + + /*First, go through dataset of inputs and check whether any input + * with the same name is already in. If so, erase the corresponding + * object before adding this new one: */ + vector::iterator object; + Input* input=NULL; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + input=(Input*)(*object); + + if (input->EnumType()==in_input->EnumType()){ + this->DeleteObject(input); + break; + } + } + this->AddObject(in_input); +} +/*}}}*/ + +/*Object functions*/ Index: /issm/trunk/src/c/DataSet/Inputs.h =================================================================== --- /issm/trunk/src/c/DataSet/Inputs.h (revision 3612) +++ /issm/trunk/src/c/DataSet/Inputs.h (revision 3612) @@ -0,0 +1,33 @@ +/*! \file Inputs.h + * \brief: header file for inputs object + */ + +#ifndef _INPUTS_H_ +#define _INPUTS_H_ + +#include "./DataSet.h" + +class Inputs: public DataSet{ + + public: + + /*constructors, destructors: {{{1*/ + Inputs(); + Inputs(int enum_type); + ~Inputs(); + /*}}}*/ + /*numerics: {{{1*/ + void GetParameterValue(bool* pvalue,int enum_type); + void GetParameterValue(int* pvalue,int enum_type); + void GetParameterValue(double* pvalue,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); + /*}}}*/ + +}; +#endif /* _INPUTS_H */ Index: /issm/trunk/src/c/Dux/Dux.h =================================================================== --- /issm/trunk/src/c/Dux/Dux.h (revision 3611) +++ /issm/trunk/src/c/Dux/Dux.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _DUX_H */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3611) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3612) @@ -99,4 +99,7 @@ InputEnum, TriaVertexInputEnum, + BoolInputEnum, + IntInputEnum, + DoubleInputEnum, /*Params: */ ParamEnum, @@ -124,4 +127,37 @@ VxEnum, VyEnum, + VzEnum, + VxAverageEnum, + VyAverageEnum, + VzAverageEnum, + VxObsEnum, + VyObsEnum, + VzObsEnum, + VxOldEnum, + VyOldEnum, + VzOldEnum, + DhDtEnum, + ThicknessEnum, + SurfaceEnum, + BedEnum, + DragCoefficientEnum, + DragPEnum, + DragQEnum, + DragTypeEnum, + RheologyBEnum, + RheologyNEnum, + MeltingRateEnum, + AccumulationRateEnum, + GeothermalFluxEnum, + ElementOnIceShelfEnum, + ElementOnBedEnum, + ElementOnWaterEnum, + ElementOnSurfaceEnum, + SurfaceAreaEnum, + WeightsEnum, + FitEnum, + AdjointxEnum, + AdjointyEnum, + PressureEnum /*}}}*/ Index: /issm/trunk/src/c/Gradjx/Gradjx.h =================================================================== --- /issm/trunk/src/c/Gradjx/Gradjx.h (revision 3611) +++ /issm/trunk/src/c/Gradjx/Gradjx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void Gradjx( Vec* pgrad_g, int numberofnodes, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type,char* control_type); + int analysis_type,int sub_analysis_type,char* control_type); #endif /* _GRADJX_H */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3611) +++ /issm/trunk/src/c/Makefile.am (revision 3612) @@ -30,6 +30,4 @@ ./objects/BamgOpts.cpp\ ./objects/Element.h\ - ./objects/ElementProperties.h\ - ./objects/ElementProperties.cpp\ ./objects/Model.h\ ./objects/Model.cpp\ @@ -67,4 +65,10 @@ ./objects/TriaVertexInput.h\ ./objects/TriaVertexInput.cpp\ + ./objects/BoolInput.h\ + ./objects/BoolInput.cpp\ + ./objects/IntInput.h\ + ./objects/IntInput.cpp\ + ./objects/DoubleInput.h\ + ./objects/DoubleInput.cpp\ ./objects/Sing.h\ ./objects/Sing.cpp\ @@ -82,6 +86,4 @@ ./objects/Input.cpp\ ./objects/Einput.h\ - ./objects/ParameterInputs.h\ - ./objects/ParameterInputs.cpp\ ./objects/Spc.cpp\ ./objects/Spc.h\ @@ -104,4 +106,6 @@ ./DataSet/DataSet.cpp\ ./DataSet/DataSet.h\ + ./DataSet/Inputs.cpp\ + ./DataSet/Inputs.h\ ./shared/shared.h\ ./shared/Alloc/alloc.h\ @@ -200,5 +204,4 @@ ./io/FetchParams.cpp\ ./io/FetchNodeSets.cpp\ - ./io/ParameterInputsInit.cpp\ ./io/pfopen.cpp\ ./io/pfclose.cpp\ @@ -272,6 +275,4 @@ ./Gradjx/Gradjx.h\ ./Gradjx/Gradjx.cpp\ - ./UpdateFromInputsx/UpdateFromInputsx.h\ - ./UpdateFromInputsx/UpdateFromInputsx.cpp\ ./UpdateInputsx/UpdateInputsx.h\ ./UpdateInputsx/UpdateInputsx.cpp\ @@ -441,6 +442,4 @@ ./objects/BamgOpts.cpp\ ./objects/Element.h\ - ./objects/ElementProperties.h\ - ./objects/ElementProperties.cpp\ ./objects/Model.h\ ./objects/Model.cpp\ @@ -478,4 +477,10 @@ ./objects/TriaVertexInput.h\ ./objects/TriaVertexInput.cpp\ + ./objects/BoolInput.h\ + ./objects/BoolInput.cpp\ + ./objects/IntInput.h\ + ./objects/IntInput.cpp\ + ./objects/DoubleInput.h\ + ./objects/DoubleInput.cpp\ ./objects/Sing.h\ ./objects/Sing.cpp\ @@ -490,9 +495,5 @@ ./objects/Numpar.h\ ./objects/Numpar.cpp\ - ./objects/Input.h\ - ./objects/Input.cpp\ ./objects/Einput.h\ - ./objects/ParameterInputs.h\ - ./objects/ParameterInputs.cpp\ ./objects/Spc.cpp\ ./objects/Spc.h\ @@ -515,4 +516,6 @@ ./DataSet/DataSet.cpp\ ./DataSet/DataSet.h\ + ./DataSet/Inputs.cpp\ + ./DataSet/Inputs.h\ ./shared/shared.h\ ./shared/Threads/issm_threads.h\ @@ -608,5 +611,4 @@ ./io/WriteParams.cpp\ ./io/FetchNodeSets.cpp\ - ./io/ParameterInputsInit.cpp\ ./io/pfopen.cpp\ ./io/pfclose.cpp\ @@ -680,6 +682,4 @@ ./Gradjx/Gradjx.h\ ./Gradjx/Gradjx.cpp\ - ./UpdateFromInputsx/UpdateFromInputsx.h\ - ./UpdateFromInputsx/UpdateFromInputsx.cpp\ ./UpdateInputsx/UpdateInputsx.h\ ./UpdateInputsx/UpdateInputsx.cpp\ Index: /issm/trunk/src/c/Misfitx/Misfitx.h =================================================================== --- /issm/trunk/src/c/Misfitx/Misfitx.h (revision 3611) +++ /issm/trunk/src/c/Misfitx/Misfitx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void Misfitx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _MISFITX_H */ Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3611) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3612) @@ -43,14 +43,12 @@ IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); - IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag"); - IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p"); - IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q"); + IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); + IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); - IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B"); - IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n"); - IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); - IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); for (i=0;inumberofelements;i++){ @@ -71,16 +69,16 @@ /*Free data : */ xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements_type); xfree((void**)&iomodel->thickness); xfree((void**)&iomodel->surface); xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->drag); - xfree((void**)&iomodel->p); - xfree((void**)&iomodel->q); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); + xfree((void**)&iomodel->rheology_B); + xfree((void**)&iomodel->rheology_n); xfree((void**)&iomodel->elementoniceshelf); xfree((void**)&iomodel->elementonwater); - xfree((void**)&iomodel->B); - xfree((void**)&iomodel->n); - xfree((void**)&iomodel->elements_type); - + } else{ // if (strcmp(meshtype,"2d")==0) @@ -88,20 +86,18 @@ /*Fetch data needed: */ IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); - IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag"); - IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p"); - IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q"); + IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); + IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); - IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); - IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B"); - IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n"); - IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); - IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); - IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); - + for (i=0;inumberofelements;i++){ if(iomodel->my_elements[i]){ @@ -120,19 +116,18 @@ /*Free data: */ xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements_type); xfree((void**)&iomodel->thickness); xfree((void**)&iomodel->surface); xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->drag); - xfree((void**)&iomodel->p); - xfree((void**)&iomodel->q); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); + xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->rheology_B); xfree((void**)&iomodel->elementoniceshelf); xfree((void**)&iomodel->elementonbed); xfree((void**)&iomodel->elementonsurface); - xfree((void**)&iomodel->elements_type); - xfree((void**)&iomodel->n); - xfree((void**)&iomodel->B); - xfree((void**)&iomodel->accumulation); - xfree((void**)&iomodel->melting); xfree((void**)&iomodel->elementonwater); + } //if (strcmp(meshtype,"2d")==0) Index: /issm/trunk/src/c/ModelProcessorx/IoModel.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3611) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.cpp (revision 3612) @@ -85,7 +85,7 @@ iomodel->drag_type=0; - iomodel->drag=NULL; - iomodel->p=NULL; - iomodel->q=NULL; + iomodel->drag_cofficient=NULL; + iomodel->drag_p=NULL; + iomodel->drag_q=NULL; @@ -102,6 +102,6 @@ iomodel->rho_ice=0; iomodel->g=0; - iomodel->n=NULL; - iomodel->B=NULL; + iomodel->rheology_n=NULL; + iomodel->rheology_B=NULL; /*!control methods: */ @@ -170,5 +170,5 @@ /*!basal: */ - iomodel->accumulation=NULL; + iomodel->accumulation_rate=NULL; iomodel->dhdt=NULL; @@ -240,7 +240,7 @@ xfree((void**)&iomodel->pressure); xfree((void**)&iomodel->temperature); - xfree((void**)&iomodel->drag); - xfree((void**)&iomodel->p); - xfree((void**)&iomodel->q); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); xfree((void**)&iomodel->elementoniceshelf); xfree((void**)&iomodel->elementonwater); @@ -253,9 +253,9 @@ xfree((void**)&iomodel->edges); xfree((void**)&iomodel->geothermalflux); - xfree((void**)&iomodel->melting); - xfree((void**)&iomodel->accumulation); + xfree((void**)&iomodel->melting_rate); + xfree((void**)&iomodel->accumulation_rate); xfree((void**)&iomodel->dhdt); - xfree((void**)&iomodel->B); - xfree((void**)&iomodel->n); + xfree((void**)&iomodel->rheology_B); + xfree((void**)&iomodel->rheology_n); xfree((void**)&iomodel->fit); xfree((void**)&iomodel->weights); Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3611) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3612) @@ -84,7 +84,7 @@ /*friction: */ int drag_type; - double* drag; - double* p; - double* q; + double* drag_coefficient; + double* drag_p; + double* drag_q; /*boundary conditions: */ @@ -101,6 +101,6 @@ double rho_water,rho_ice; double g; - double* B; - double* n; + double* rheology_B; + double* rheology_n; /*numerical parameters: */ @@ -170,6 +170,6 @@ /*basal: */ - double* melting; - double* accumulation; + double* melting_rate; + double* accumulation_rate; double* dhdt; Index: /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h =================================================================== --- /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 3611) +++ /issm/trunk/src/c/PenaltyConstraintsx/PenaltyConstraintsx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void PenaltyConstraintsx(int* pconverged, int* pnum_unstable_constraints, DataSet* elements,DataSet* nodes, DataSet* vertices, - DataSet* loads,DataSet* materials, DataSet* parameters,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + DataSet* loads,DataSet* materials, DataSet* parameters,int analysis_type,int sub_analysis_type); #endif /* _PENALTYCONSTRAINTSX_H */ Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 3611) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.cpp (revision 3612) @@ -11,5 +11,5 @@ void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, DataSet* parameters, - int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + int kflag,int pflag,int analysis_type,int sub_analysis_type){ int i; Index: /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h =================================================================== --- /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 3611) +++ /issm/trunk/src/c/PenaltySystemMatricesx/PenaltySystemMatricesx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void PenaltySystemMatricesx(Mat Kgg, Vec pg,double* pkmax, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, DataSet* parameters, - int kflag,int pflag,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int kflag,int pflag,int analysis_type,int sub_analysis_type); #endif /* _PENALTYSYSTEMMATRICESX_H */ Index: /issm/trunk/src/c/Qmux/Qmux.h =================================================================== --- /issm/trunk/src/c/Qmux/Qmux.h (revision 3611) +++ /issm/trunk/src/c/Qmux/Qmux.h (revision 3612) @@ -15,7 +15,7 @@ void SpawnCoreSerial(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, mxArray* model,mxArray* inputs,int analysis_type,int sub_analysis_type,int counter); #else -void Qmux(Model* model,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); -void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,ParameterInputs* inputs,int analysis_type,int sub_analysis_type,int counter); -void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,DataSet* processed_results,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void Qmux(Model* model,int analysis_type,int sub_analysis_type); +void SpawnCoreParallel(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, Model* model,int analysis_type,int sub_analysis_type,int counter); +void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Model* model, DataSet* results,DataSet* processed_results,int analysis_type,int sub_analysis_type); #endif Index: /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h =================================================================== --- /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h (revision 3611) +++ /issm/trunk/src/c/SurfaceAreax/SurfaceAreax.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void SurfaceAreax( double* pS, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, DataSet* parameters, - ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int analysis_type,int sub_analysis_type); #endif /* _SURFACEAREAX_H */ Index: /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp (revision 3611) +++ /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.cpp (revision 3612) @@ -11,5 +11,5 @@ void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, DataSet* parameters, - int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type){ + int kflag,int pflag,int connectivity,int numberofdofspernode,int analysis_type,int sub_analysis_type){ int i; @@ -30,4 +30,5 @@ loads->Configure(elements, loads, nodes,vertices, materials,parameters); parameters->Configure(elements,loads, nodes,vertices, materials,parameters); + /*Get size of matrix: */ Index: /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h =================================================================== --- /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h (revision 3611) +++ /issm/trunk/src/c/SystemMatricesx/SystemMatricesx.h (revision 3612) @@ -11,5 +11,5 @@ /* local prototypes: */ void SystemMatricesx(Mat* pKgg, Vec* ppg,DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads,DataSet* materials, DataSet* parameters, - int kflag,int pflag,int connectivity,int numberofdofspernode,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); + int kflag,int pflag,int connectivity,int numberofdofspernode,int analysis_type,int sub_analysis_type); #endif /* _SYSTEMMATRICESX_H */ Index: /issm/trunk/src/c/UpdateInputsx/UpdateInputsx.cpp =================================================================== --- /issm/trunk/src/c/UpdateInputsx/UpdateInputsx.cpp (revision 3611) +++ /issm/trunk/src/c/UpdateInputsx/UpdateInputsx.cpp (revision 3612) @@ -31,6 +31,4 @@ elements->UpdateInputs(serial_solution,analysis_type,sub_analysis_type); - elements->Echo(); - /*Free ressources:*/ xfree((void**)&serial_solution); Index: /issm/trunk/src/c/UpdateVertexPositionsx/UpdateVertexPositionsx.h =================================================================== --- /issm/trunk/src/c/UpdateVertexPositionsx/UpdateVertexPositionsx.h (revision 3611) +++ /issm/trunk/src/c/UpdateVertexPositionsx/UpdateVertexPositionsx.h (revision 3612) @@ -7,5 +7,4 @@ #include "../DataSet/DataSet.h" -#include "../objects/ParameterInputs.h" /* local prototypes: */ Index: sm/trunk/src/c/io/ParameterInputsInit.cpp =================================================================== --- /issm/trunk/src/c/io/ParameterInputsInit.cpp (revision 3611) +++ (revision ) @@ -1,94 +1,0 @@ -/*!\file: ParameterInputsInits.cpp - * \brief inputs routine is in the io for initializing a ParameterInputs structure out of the matlab workspace. - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - - -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "../include/macros.h" -#include "./io.h" - -#ifdef _SERIAL_ -void ParameterInputs::Init( void* vinput_handle){ - - int i; - int nelems; - mxArray* input=NULL; - const char* field_name=NULL; - ConstDataHandle input_handle=NULL; - int type; - - char* string=NULL; - char* name=NULL; - double scalar; - int integer; - int ndof; - int numberofnodes; - double* vector=NULL; - - input_handle=(ConstDataHandle)vinput_handle; - - /*number of inputs: */ - nelems=mxGetNumberOfElements(input_handle); - - /*Go through ParameterInputs matlab structure, and recover each input object. Use it for fill - * the ParameterInputs* "C" class object: */ - - for(i=0;iAdd(name,string); - xfree((void**)&string); - break; - case 1: - //integer - integer=(int)*(mxGetPr(mxGetField(input,0,"integer"))); - this->Add(name,integer); - break; - case 2: - //double - scalar=(double)*(mxGetPr(mxGetField(input,0,"double"))); - this->Add(name,scalar); - break; - case 3: - //double vector - vector=mxGetPr(mxGetField(input,0,"vector")); - ndof=(int)*(mxGetPr(mxGetField(input,0,"ndof"))); - numberofnodes=(int)*(mxGetPr(mxGetField(input,0,"numberofnodes"))); - this->Add(name,vector,ndof,numberofnodes); - break; - default: - ISSMERROR("%s%i"," unknow parameter input type ",type); - } - - } - - -} -#else -void ParameterInputs::Init( void* vinput_handle){ -} -#endif //#ifdef _SERIAL_ Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 3611) +++ /issm/trunk/src/c/issm.h (revision 3612) @@ -39,5 +39,4 @@ #include "./ConfigureObjectsx/ConfigureObjectsx.h" #include "./SystemMatricesx/SystemMatricesx.h" -#include "./UpdateFromInputsx/UpdateFromInputsx.h" #include "./UpdateInputsx/UpdateInputsx.h" #include "./UpdateGeometryx/UpdateGeometryx.h" Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 3612) @@ -243,78 +243,8 @@ } /*}}}*/ -/*FUNCTION Beam::UpdateFromInputs{{{1*/ -void Beam::UpdateFromInputs(void* vinputs){ - - int dofs[1]={0}; - double temperature_list[2]; - double temperature_average; - double B_list[2]; - double B_average; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matpar* matpar=NULL; - Matice* matice=NULL; - Numpar* numpar=NULL; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*Update internal data if inputs holds new values: */ - //if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n"); - //For now - if(this->properties.h) inputs->Recover("thickness",&this->properties.h[0],1,dofs,2,(void**)nodes); - //Later - /* - if(inputs->Recover("thickness",&new_h[0],1,dofs,2,(void**)nodes)){ - //density, needed later: - double di=(this->matpar->GetRhoIce()/this->matpar->GetRhoWater()); - //Go through grids: - for (i=0;i<2;i++){ - if(nodes[i]->IsOnShelf()){ - this->b[i]=this->b[i]-di*(new_h[i]-h[i]); //hydrostatic equilibrium; - } - this->s[i]=this->b[i]+new_h[i]; - this->h[i]=new_h[i]; - } - } - */ - if (this->properties.s) inputs->Recover("surface",&this->properties.s[0],1,dofs,2,(void**)nodes); - if (this->properties.b) inputs->Recover("bed",&this->properties.b[0],1,dofs,2,(void**)nodes); - if (this->properties.k) inputs->Recover("drag",&this->properties.k[0],1,dofs,2,(void**)nodes); - if (this->properties.melting) inputs->Recover("melting",&this->properties.melting[0],1,dofs,2,(void**)nodes); - if (this->properties.accumulation) inputs->Recover("accumulation",&this->properties.accumulation[0],1,dofs,2,(void**)nodes); - if (this->properties.geothermalflux) inputs->Recover("geothermalflux",&this->properties.geothermalflux[0],1,dofs,2,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,2,(void**)nodes)){ - temperature_average=(temperature_list[0]+temperature_list[1])/2.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - if(inputs->Recover("temperature",&temperature_list[0],1,dofs,2,(void**)nodes)){ - temperature_average=(temperature_list[0]+temperature_list[1])/2.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - if(inputs->Recover("B",&B_list[0],1,dofs,2,(void**)nodes)){ - B_average=(B_list[0]+B_list[1])/2.0; - matice->SetB(B_average); - } - -} -/*}}}*/ /*Object functions*/ /*FUNCTION Beam::ComputeBasalStress{{{1*/ -void Beam::ComputeBasalStress(Vec eps,void* inputs,int analysis_type,int sub_analysis_type){ +void Beam::ComputeBasalStress(Vec eps,int analysis_type,int sub_analysis_type){ ISSMERROR("Not implemented yet"); @@ -323,5 +253,5 @@ /*}}}*/ /*FUNCTION Beam::ComputePressure{{{1*/ -void Beam::ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){ +void Beam::ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type){ int i; @@ -366,5 +296,5 @@ /*}}}*/ /*FUNCTION Beam::ComputeStrainRate{{{1*/ -void Beam::ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type){ +void Beam::ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type){ ISSMERROR("Not implemented yet"); @@ -379,5 +309,5 @@ /*FUNCTION Beam::CreateKMatrix{{{1*/ -void Beam::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Beam::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -399,5 +329,5 @@ /*FUNCTION Beam::CreateKMatrixDiagnosticHutter{{{1*/ -void Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Beam::CreateKMatrixDiagnosticHutter(Mat Kgg,int analysis_type,int sub_analysis_type){ @@ -432,5 +362,5 @@ /*}}}*/ /*FUNCTION Beam::CreatePVector{{{1*/ -void Beam::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ +void Beam::CreatePVector(Vec pg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ @@ -450,5 +380,5 @@ /*FUNCTION Beam::CreatePVectorDiagnosticHutter{{{1*/ -void Beam::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Beam::CreatePVectorDiagnosticHutter( Vec pg, int analysis_type,int sub_analysis_type){ int i,j,k; @@ -699,5 +629,5 @@ /*}}}*/ /*FUNCTION Beam::GetOnBed{{{1*/ -int Beam::GetOnBed(){ +bool Beam::GetOnBed(){ ISSMERROR(" not supported yet!"); } @@ -714,5 +644,5 @@ /*}}}*/ /*FUNCTION Beam::GetShelf{{{1*/ -int Beam::GetShelf(){ +bool Beam::GetShelf(){ ISSMERROR(" not supported yet!"); } Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 3611) +++ /issm/trunk/src/c/objects/Beam.h (revision 3612) @@ -11,6 +11,4 @@ #include "./Matice.h" #include "./Matpar.h" -#include "./ParameterInputs.h" -#include "./ElementProperties.h" #include "../ModelProcessorx/IoModel.h" #include "./Hook.h" @@ -20,5 +18,4 @@ class Element; class Hook; -class ElementProperties; class Beam: public Element{ @@ -34,5 +31,5 @@ Hook hnumpar; //hook to 1 numpar - ElementProperties properties; + Inputs* inputs; public: @@ -40,6 +37,6 @@ /*constructors, destructors: {{{1*/ Beam(); - Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id, ElementProperties* beam_properties); - Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, ElementProperties* beam_properties); + Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id); + Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs); Beam(int i, IoModel* iomodel); ~Beam(); @@ -61,23 +58,22 @@ /*}}}*/ /*numerics: {{{1*/ - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); - void UpdateFromInputs(void* inputs); + void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHutter(Mat Kgg,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); - void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHutter(Vec pg,int analysis_type,int sub_analysis_type); void* GetMatPar(); - void ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type); - void ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type); - void ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type); + void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); + void ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type); + void ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type); void GetNodes(void** vpnodes); /*}}}*/ /*not implemented: {{{1*/ - int GetShelf(); - int GetOnBed(); + bool GetShelf(); + bool GetOnBed(); void GetBedList(double*); void GetThicknessList(double* thickness_list); Index: /issm/trunk/src/c/objects/BoolInput.cpp =================================================================== --- /issm/trunk/src/c/objects/BoolInput.cpp (revision 3612) +++ /issm/trunk/src/c/objects/BoolInput.cpp (revision 3612) @@ -0,0 +1,140 @@ +/*!\file BoolInput.c + * \brief: implementation of the BoolInput object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "./BoolInput.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../include/typedefs.h" +#include "../include/types.h" +#include "../include/macros.h" + +/*Object constructors and destructor*/ +/*FUNCTION BoolInput::BoolInput(){{{1*/ +BoolInput::BoolInput(){ + return; +} +/*}}}*/ +/*FUNCTION BoolInput::BoolInput(double* values){{{1*/ +BoolInput::BoolInput(int in_enum_type,IssmBool in_value){ + + enum_type=in_enum_type; + value=in_value; +} +/*}}}*/ +/*FUNCTION BoolInput::~BoolInput(){{{1*/ +BoolInput::~BoolInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION BoolInput::copy{{{1*/ +Object* BoolInput::copy() { + + return new BoolInput(this->enum_type,this->value); + +} +/*}}}*/ +/*FUNCTION BoolInput::DeepEcho{{{1*/ +void BoolInput::DeepEcho(void){ + + printf("BoolInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %s\n",value?"true":"false"); +} +/*}}}*/ +/*FUNCTION BoolInput::Demarshall{{{1*/ +void BoolInput::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(&value,marshalled_dataset,sizeof(value));marshalled_dataset+=sizeof(value); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION BoolInput::Echo {{{1*/ +void BoolInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION BoolInput::Enum{{{1*/ +int BoolInput::Enum(void){ + + return BoolInputEnum; + +} +/*}}}*/ +/*FUNCTION BoolInput::EnumType{{{1*/ +int BoolInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION BoolInput::GetId{{{1*/ +int BoolInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION BoolInput::GetName{{{1*/ +char* BoolInput::GetName(void){ + return "boolinput"; +} +/*}}}*/ +/*FUNCTION BoolInput::Marshall{{{1*/ +void BoolInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of BoolInput: */ + enum_type=BoolInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall BoolInput data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&value,sizeof(value));marshalled_dataset+=sizeof(value); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION BoolInput::MarshallSize{{{1*/ +int BoolInput::MarshallSize(){ + + return sizeof(value)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION BoolInput::MyRank{{{1*/ +int BoolInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/BoolInput.h =================================================================== --- /issm/trunk/src/c/objects/BoolInput.h (revision 3612) +++ /issm/trunk/src/c/objects/BoolInput.h (revision 3612) @@ -0,0 +1,43 @@ +/*! \file BoolInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" +#include "../include/types.h" + +#ifndef _BOOLINPUT_H_ +#define _BOOLINPUT_H_ + +class BoolInput: public Input{ + + private: + /*just hold 3 values for 3 vertices: */ + int enum_type; + IssmBool value; + + public: + + /*constructors, destructors: {{{1*/ + BoolInput(); + BoolInput(int enum_type,IssmBool value); + ~BoolInput(); + /*}}}*/ + /*object management: {{{1*/ + void DeepEcho(); + void Echo(); + int GetId(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + char* GetName(); + void Demarshall(char** pmarshalled_dataset); + int Enum(); + int MyRank(); + Object* copy(); + int EnumType(); + + /*}}}*/ + /*numerics: {{{1*/ + /*}}}*/ + +}; +#endif /* _BOOLINPUT_H */ Index: /issm/trunk/src/c/objects/DoubleInput.cpp =================================================================== --- /issm/trunk/src/c/objects/DoubleInput.cpp (revision 3612) +++ /issm/trunk/src/c/objects/DoubleInput.cpp (revision 3612) @@ -0,0 +1,140 @@ +/*!\file DoubleInput.c + * \brief: implementation of the DoubleInput object + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include "./DoubleInput.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../include/typedefs.h" +#include "../include/types.h" +#include "../include/macros.h" + +/*Object constructors and destructor*/ +/*FUNCTION DoubleInput::DoubleInput(){{{1*/ +DoubleInput::DoubleInput(){ + return; +} +/*}}}*/ +/*FUNCTION DoubleInput::DoubleInput(double* values){{{1*/ +DoubleInput::DoubleInput(int in_enum_type,IssmDouble in_value){ + + enum_type=in_enum_type; + value=in_value; +} +/*}}}*/ +/*FUNCTION DoubleInput::~DoubleInput(){{{1*/ +DoubleInput::~DoubleInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION DoubleInput::copy{{{1*/ +Object* DoubleInput::copy() { + + return new DoubleInput(this->enum_type,this->value); + +} +/*}}}*/ +/*FUNCTION DoubleInput::DeepEcho{{{1*/ +void DoubleInput::DeepEcho(void){ + + printf("DoubleInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %g\n",this->value); +} +/*}}}*/ +/*FUNCTION DoubleInput::Demarshall{{{1*/ +void DoubleInput::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(&value,marshalled_dataset,sizeof(value));marshalled_dataset+=sizeof(value); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION DoubleInput::Echo {{{1*/ +void DoubleInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION DoubleInput::Enum{{{1*/ +int DoubleInput::Enum(void){ + + return DoubleInputEnum; + +} +/*}}}*/ +/*FUNCTION DoubleInput::EnumType{{{1*/ +int DoubleInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION DoubleInput::GetId{{{1*/ +int DoubleInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION DoubleInput::GetName{{{1*/ +char* DoubleInput::GetName(void){ + return "doubleinput"; +} +/*}}}*/ +/*FUNCTION DoubleInput::Marshall{{{1*/ +void DoubleInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of DoubleInput: */ + enum_type=DoubleInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall DoubleInput data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&value,sizeof(value));marshalled_dataset+=sizeof(value); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION DoubleInput::MarshallSize{{{1*/ +int DoubleInput::MarshallSize(){ + + return sizeof(value)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION DoubleInput::MyRank{{{1*/ +int DoubleInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/DoubleInput.h =================================================================== --- /issm/trunk/src/c/objects/DoubleInput.h (revision 3612) +++ /issm/trunk/src/c/objects/DoubleInput.h (revision 3612) @@ -0,0 +1,43 @@ +/*! \file DoubleInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" +#include "../include/types.h" + +#ifndef _DOUBLEINPUT_H_ +#define _DOUBLEINPUT_H_ + +class DoubleInput: public Input{ + + private: + /*just hold 3 values for 3 vertices: */ + int enum_type; + IssmDouble value; + + public: + + /*constructors, destructors: {{{1*/ + DoubleInput(); + DoubleInput(int enum_type,IssmDouble value); + ~DoubleInput(); + /*}}}*/ + /*object management: {{{1*/ + void DeepEcho(); + void Echo(); + int GetId(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + char* GetName(); + void Demarshall(char** pmarshalled_dataset); + int Enum(); + int MyRank(); + Object* copy(); + int EnumType(); + + /*}}}*/ + /*numerics: {{{1*/ + /*}}}*/ + +}; +#endif /* _DOUBLEINPUT_H */ Index: sm/trunk/src/c/objects/Einput.h =================================================================== --- /issm/trunk/src/c/objects/Einput.h (revision 3611) +++ (revision ) @@ -1,20 +1,0 @@ -/*!\file: Einput.h - * \brief abstract class for Einput object - */ - - -#ifndef _EINPUT_H_ -#define _EINPUT_H_ - -#include "./Object.h" - -class Einput: public Object{ - - public: - - virtual ~Einput(){}; - virtual int Enum()=0; //object Enum - virtual int EnumType()=0; //type of input (vx,vy?) - -}; -#endif Index: sm/trunk/src/c/objects/Element.cpp =================================================================== --- /issm/trunk/src/c/objects/Element.cpp (revision 3611) +++ (revision ) @@ -1,12 +1,0 @@ -/*!\file: Element.cpp - * \brief: implementation of the Enum routine, which is the only non-virtual routine from - * the Object derived class Element. - */ - -#include "./Element.h" -#include "../EnumDefinitions/EnumDefinitions.h" - -int Element::Enum(void){ - - return ElementEnum; -} Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 3611) +++ /issm/trunk/src/c/objects/Element.h (revision 3612) @@ -17,26 +17,26 @@ virtual ~Element(){}; + virtual int Enum()=0; virtual void Configure(void* loads,void* nodes,void* materials,void* parameters)=0; - virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type)=0; - virtual int Enum()=0; - virtual void UpdateFromInputs(void* inputs)=0; + + virtual void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type)=0; + virtual void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type)=0; virtual void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type)=0; virtual void GetNodes(void** nodes)=0; virtual void* GetMatPar()=0; - virtual int GetShelf()=0; - virtual int GetOnBed()=0; + virtual bool GetShelf()=0; + virtual bool GetOnBed()=0; virtual void GetThicknessList(double* thickness_list)=0; virtual void GetBedList(double* bed_list)=0; - virtual void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type)=0; - virtual void GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual double Misfit(void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual double CostFunction(void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual double SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void ComputePressure(Vec p_g, void* inputs,int analysis_type,int sub_analysis_type)=0; - virtual void ComputeStrainRate(Vec eps, void* inputs,int analysis_type,int sub_analysis_type)=0; + virtual void Du(Vec du_g,int analysis_type,int sub_analysis_type)=0; + virtual void Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type)=0; + virtual void GradjDrag(Vec grad_g,int analysis_type,int sub_analysis_type)=0; + virtual void GradjB(Vec grad_g,int analysis_type,int sub_analysis_type)=0; + virtual double Misfit(int analysis_type,int sub_analysis_type)=0; + virtual double CostFunction(int analysis_type,int sub_analysis_type)=0; + virtual double SurfaceArea(int analysis_type,int sub_analysis_type)=0; + virtual void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type)=0; + virtual void ComputePressure(Vec p_g, int analysis_type,int sub_analysis_type)=0; + virtual void ComputeStrainRate(Vec eps, int analysis_type,int sub_analysis_type)=0; virtual double MassFlux(double* segment,double* ug)=0; Index: sm/trunk/src/c/objects/ElementProperties.cpp =================================================================== --- /issm/trunk/src/c/objects/ElementProperties.cpp (revision 3611) +++ (revision ) @@ -1,526 +1,0 @@ -/*!\file ElementProperties.c - * \brief: implementation of the ElementProperties object - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "stdio.h" -#include "./ElementProperties.h" -#include -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../shared/shared.h" -#include "../DataSet/DataSet.h" -#include "../include/typedefs.h" -#include "../include/macros.h" - -/*Object constructors and destructor*/ -/*FUNCTION ElementProperties default constructor {{{1*/ -ElementProperties::ElementProperties(){ - - int i; - - /*initialize to 0 every value: */ - this->numnodes=0; - this->h=NULL; - this->s=NULL; - this->b=NULL; - this->k=NULL; - this->melting=NULL; - this->accumulation=NULL; - this->geothermalflux=NULL; - - this->friction_type=UNDEF; - this->p=UNDEF; - this->q=UNDEF; - this->shelf=UNDEF; - this->onbed=UNDEF; - this->onwater=UNDEF; - this->onsurface=UNDEF; - this->collapse=UNDEF; - this->thermal_steadystate=UNDEF; - return; -} -/*}}}*/ -/*FUNCTION ElementProperties copy constructor {{{1*/ -ElementProperties::ElementProperties(ElementProperties* in){ - - int i; - - this->numnodes=in->numnodes; - - if(in->h){ - this->h=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->h[i]=in->h[i]; - } - else this->h=NULL; - - this->numnodes=in->numnodes; - if(in->s){ - this->s=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->s[i]=in->s[i]; - } - else this->s=NULL; - - - this->numnodes=in->numnodes; - if(in->b){ - this->b=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->b[i]=in->b[i]; - } - else this->b=NULL; - - - this->numnodes=in->numnodes; - if(in->k){ - this->k=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->k[i]=in->k[i]; - } - else this->k=NULL; - - - this->numnodes=in->numnodes; - if(in->melting){ - this->melting=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->melting[i]=in->melting[i]; - } - else this->melting=NULL; - - - this->numnodes=in->numnodes; - if(in->accumulation){ - this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->accumulation[i]=in->accumulation[i]; - } - else this->accumulation=NULL; - - - this->numnodes=in->numnodes; - if(in->geothermalflux){ - this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double)); - for(i=0;inumnodes;i++)this->geothermalflux[i]=in->geothermalflux[i]; - } - else this->geothermalflux=NULL; - - - this->friction_type=in->friction_type; - this->p=in->p; - this->q=in->q; - this->shelf=in->shelf; - this->onbed=in->onbed; - this->onwater=in->onwater; - this->onsurface=in->onsurface; - this->collapse=in->collapse; - this->thermal_steadystate=in->thermal_steadystate; - -} -/*}}}*/ -/*FUNCTION ElementProperties constructor {{{1*/ - -ElementProperties::ElementProperties(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b, - double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation, - double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p, - double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, - int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){ - - this->Init(elementproperties_numnodes, elementproperties_h,elementproperties_s,elementproperties_b, - elementproperties_k,elementproperties_melting,elementproperties_accumulation, - elementproperties_geothermalflux, elementproperties_friction_type,elementproperties_p, - elementproperties_q, elementproperties_shelf, elementproperties_onbed, elementproperties_onwater, - elementproperties_onsurface, elementproperties_collapse, elementproperties_thermal_steadystate); -} -/*}}}*/ -/*FUNCTION ElementProperties Initialize propreties, used by constructor{{{1*/ -void ElementProperties::Init(int elementproperties_numnodes, double* elementproperties_h,double* elementproperties_s,double* elementproperties_b, - double* elementproperties_k,double* elementproperties_melting,double* elementproperties_accumulation, - double* elementproperties_geothermalflux, int elementproperties_friction_type,double elementproperties_p, - double elementproperties_q, int elementproperties_shelf, int elementproperties_onbed, bool elementproperties_onwater, - int elementproperties_onsurface, int elementproperties_collapse, int elementproperties_thermal_steadystate){ - - int i; - - this->numnodes=elementproperties_numnodes; - - this->h=(double*)xmalloc(this->numnodes*sizeof(double)); - this->s=(double*)xmalloc(this->numnodes*sizeof(double)); - this->b=(double*)xmalloc(this->numnodes*sizeof(double)); - this->k=(double*)xmalloc(this->numnodes*sizeof(double)); - this->melting=(double*)xmalloc(this->numnodes*sizeof(double)); - this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double)); - this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double)); - - for(i=0;inumnodes;i++){ - if(elementproperties_h) this->h[i]=elementproperties_h[i]; - else this->h[i]=UNDEF; - if(elementproperties_s) this->s[i]=elementproperties_s[i]; - else this->s[i]=UNDEF; - if(elementproperties_b) this->b[i]=elementproperties_b[i]; - else this->b[i]=UNDEF; - if(elementproperties_k) this->k[i]=elementproperties_k[i]; - else this->k[i]=UNDEF; - if(elementproperties_melting) this->melting[i]=elementproperties_melting[i]; - else this->melting[i]=UNDEF; - if(elementproperties_accumulation) this->accumulation[i]=elementproperties_accumulation[i]; - else this->accumulation[i]=UNDEF; - if(elementproperties_geothermalflux)this->geothermalflux[i]=elementproperties_geothermalflux[i]; - else this->geothermalflux[i]=UNDEF; - } - - this->friction_type=elementproperties_friction_type; - this->p=elementproperties_p; - this->q=elementproperties_q; - this->shelf=elementproperties_shelf; - this->onbed=elementproperties_onbed; - this->onwater=elementproperties_onwater; - this->onsurface=elementproperties_onsurface; - this->collapse=elementproperties_collapse; - this->thermal_steadystate=elementproperties_thermal_steadystate; - - return; -} -/*}}}*/ -/*FUNCTION ElementProperties Init propreties from iomodel{{{1*/ -void ElementProperties::Init(int numberofnodes, int* vertex_ids, int element_id, IoModel* iomodel){ //i is the node index - - int i,k; - - /*Initialize number of nodes*/ - this->numnodes=numberofnodes; - - /*Allocate node properties*/ - if (iomodel->thickness) this->h=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->surface) this->s=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->bed) this->b=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->drag) this->k=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->melting) this->melting=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->accumulation) this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double)); - if (iomodel->geothermalflux)this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double)); - - /*Get node properties*/ - for(i=0;inumnodes;i++){ - - //get index of the vertex in C indexing - k=vertex_ids[i]-1; - - //Get properties from iomodel - if(iomodel->thickness) this->h[i]=iomodel->thickness[k]; - if(iomodel->surface) this->s[i]=iomodel->surface[k]; - if(iomodel->bed) this->b[i]=iomodel->bed[k]; - if(iomodel->drag) this->k[i]=iomodel->drag[k]; - if(iomodel->melting) this->melting[i]=iomodel->melting[k]; - if(iomodel->accumulation) this->accumulation[i]=iomodel->accumulation[k]; - if(iomodel->geothermalflux)this->geothermalflux[i]=iomodel->geothermalflux[k]; - } - - /*Get element properties from iomodel*/ - if (iomodel->p) this->p=iomodel->p[element_id-1]; - if (iomodel->q) this->q=iomodel->q[element_id-1]; - if (iomodel->elementoniceshelf)this->shelf=(int)iomodel->elementoniceshelf[element_id-1]; - if (iomodel->elementonbed) this->onbed=(int)iomodel->elementonbed[element_id-1]; - if (iomodel->elementonwater) this->onwater=(bool)iomodel->elementonwater[element_id-1]; - if (iomodel->elementonsurface) this->onsurface=(int)iomodel->elementonsurface[element_id-1]; - - /*Get other properties from iomodel*/ - this->friction_type=iomodel->drag_type; - - /*WARNING: collapse and thermal_steadystate are left undefined*/ - - return; -} -/*}}}*/ -/*FUNCTION ElementProperties destructor {{{1*/ -ElementProperties::~ElementProperties(){ //destructor - - if(this->h)xfree((void**)&this->h); - if(this->s)xfree((void**)&this->s); - if(this->b)xfree((void**)&this->b); - if(this->k)xfree((void**)&this->k); - if(this->melting)xfree((void**)&this->melting); - if(this->accumulation)xfree((void**)&this->accumulation); - if(this->geothermalflux)xfree((void**)&this->geothermalflux); - - return; -} -/*}}}*/ -/*FUNCTION ElementProperties Marshall{{{1*/ -void ElementProperties::Marshall(char** pmarshalled_dataset){ - - char* marshalled_dataset=NULL; - int enum_type=0; - int nill=0; - - /*recover marshalled_dataset: */ - marshalled_dataset=*pmarshalled_dataset; - - /*get enum type of ElementProperties: */ - enum_type=ElementPropertiesEnum; - - /*marshall enum: */ - memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); - - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - - if(this->h){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->h,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - if(this->s){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->s,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - if(this->b){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->b,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - if(this->k){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->k,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - if(this->melting){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->melting,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - if(this->accumulation){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->accumulation,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - - if(this->geothermalflux){ - memcpy(marshalled_dataset,&this->numnodes,sizeof(this->numnodes));marshalled_dataset+=sizeof(this->numnodes); - memcpy(marshalled_dataset,this->geothermalflux,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - else{ memcpy(marshalled_dataset,&nill,sizeof(nill));marshalled_dataset+=sizeof(nill); } - - memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); - memcpy(marshalled_dataset,&p,sizeof(p));marshalled_dataset+=sizeof(p); - memcpy(marshalled_dataset,&q,sizeof(q));marshalled_dataset+=sizeof(q); - memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf); - memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed); - memcpy(marshalled_dataset,&onwater,sizeof(onwater));marshalled_dataset+=sizeof(onwater); - memcpy(marshalled_dataset,&onsurface,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(marshalled_dataset,&collapse,sizeof(collapse));marshalled_dataset+=sizeof(collapse); - memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate); - - *pmarshalled_dataset=marshalled_dataset; - return; -} -/*}}}*/ -/*FUNCTION ElementProperties MarshallSize{{{1*/ -int ElementProperties::MarshallSize(){ - - int size; - - size=sizeof(numnodes); - - if(this->h)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->s)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->b)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->k)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->melting)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->accumulation)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - if(this->geothermalflux)size+=sizeof(int)+numnodes*sizeof(double); else size+=sizeof(int); - - size+= sizeof(friction_type)+ - sizeof(p)+ - sizeof(q)+ - sizeof(shelf)+ - sizeof(onbed)+ - sizeof(onwater)+ - sizeof(onsurface)+ - sizeof(collapse)+ - sizeof(thermal_steadystate)+ - sizeof(int); //sizeof(int) for enum type - - return size; - -} -/*}}}*/ -/*FUNCTION ElementProperties Demarshall{{{1*/ -void ElementProperties::Demarshall(char** pmarshalled_dataset){ - - char* marshalled_dataset=NULL; - int nonnill=0; - int enum_type; - - /*recover marshalled_dataset: */ - marshalled_dataset=*pmarshalled_dataset; - - /* Recover enum_type as ElementProperties is not directly called by dataset */ - memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int); - - memcpy(&numnodes,marshalled_dataset,sizeof(numnodes));marshalled_dataset+=sizeof(numnodes); - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->h=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->h,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->s=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->s,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->b=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->b,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->k=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->k,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->melting=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->melting,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->accumulation=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->accumulation,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - memcpy(&nonnill,marshalled_dataset,sizeof(nonnill));marshalled_dataset+=sizeof(nonnill); - if(nonnill){ - this->geothermalflux=(double*)xmalloc(this->numnodes*sizeof(double)); - memcpy(this->geothermalflux,marshalled_dataset,this->numnodes*sizeof(double));marshalled_dataset+=this->numnodes*sizeof(double); - } - - - memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type); - memcpy(&p,marshalled_dataset,sizeof(p));marshalled_dataset+=sizeof(p); - memcpy(&q,marshalled_dataset,sizeof(q));marshalled_dataset+=sizeof(q); - memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf); - memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed); - memcpy(&onwater,marshalled_dataset,sizeof(onwater));marshalled_dataset+=sizeof(onwater); - memcpy(&onsurface,marshalled_dataset,sizeof(onsurface));marshalled_dataset+=sizeof(onsurface); - memcpy(&collapse,marshalled_dataset,sizeof(collapse));marshalled_dataset+=sizeof(collapse); - memcpy(&thermal_steadystate,marshalled_dataset,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate); - - /*return: */ - *pmarshalled_dataset=marshalled_dataset; - return; -} -/*}}}*/ -/*FUNCTION ElementProperties Echo{{{1*/ -void ElementProperties::Echo(void){ - - int i; - - printf("Element properties:\n"); - if (h) {printf(" h=[ "); for(i=0;ih[i]);printf("]\n");} - if (s) {printf(" s=[ "); for(i=0;is[i]);printf("]\n");} - if (b) {printf(" b=[ "); for(i=0;ib[i]);printf("]\n");} - if (k) {printf(" k=[ "); for(i=0;ik[i]);printf("]\n");} - if (melting) {printf(" melting=[ "); for(i=0;imelting[i]);printf("]\n");} - if (accumulation) {printf(" accumulation=[ "); for(i=0;iaccumulation[i]);printf("]\n");} - if (geothermalflux) {printf(" geothermalflux=[ "); for(i=0;igeothermalflux[i]);printf("]\n");} - printf(" friction_type: %i\n",friction_type); - printf(" p: %g\n",p); - printf(" q: %g\n",q); - printf(" shelf: %i\n",shelf); - printf(" onbed: %i\n",onbed); - printf(" onwater: %i\n",onwater); - printf(" onsurface: %i\n",onsurface); - printf(" collapse: %i\n",collapse); - printf(" thermal_steadystate: %i\n",thermal_steadystate); -} -/*}}}*/ -/*FUNCTION ElementProperties DeepEcho{{{1*/ -void ElementProperties::DeepEcho(void){ - this->Echo(); -} -/*}}}*/ -/*FUNCTION ElementProperties Spawn{{{1*/ -ElementProperties* ElementProperties::Spawn(int* indices, int numindices){ - - int i; - - int sub_numnodes=0; - double* sub_h=NULL; - double* sub_s=NULL; - double* sub_b=NULL; - double* sub_k=NULL; - double* sub_melting=NULL; - double* sub_accumulation=NULL; - double* sub_geothermalflux=NULL; - - int sub_friction_type; - double sub_p; - double sub_q; - int sub_shelf; - int sub_onbed; - bool sub_onwater; - int sub_onsurface; - int sub_collapse; - int sub_thermal_steadystate; - - ElementProperties* sub_properties = NULL; - - sub_numnodes=numindices; - - /*go pickup the correct objects, for the correct ids: */ - if(this->h)sub_h=(double*)xmalloc(numindices*sizeof(double)); - if(this->s)sub_s=(double*)xmalloc(numindices*sizeof(double)); - if(this->b)sub_b=(double*)xmalloc(numindices*sizeof(double)); - if(this->k)sub_k=(double*)xmalloc(numindices*sizeof(double)); - if(this->melting)sub_melting=(double*)xmalloc(numindices*sizeof(double)); - if(this->accumulation)sub_accumulation=(double*)xmalloc(numindices*sizeof(double)); - if(this->geothermalflux)sub_geothermalflux=(double*)xmalloc(numindices*sizeof(double)); - - for(i=0;ih)sub_h[i]=this->h[indices[i]]; - if(this->s)sub_s[i]=this->s[indices[i]]; - if(this->b)sub_b[i]=this->b[indices[i]]; - if(this->k)sub_k[i]=this->k[indices[i]]; - if(this->melting)sub_melting[i]=this->melting[indices[i]]; - if(this->accumulation)sub_accumulation[i]=this->accumulation[indices[i]]; - if(this->geothermalflux)sub_geothermalflux[i]=this->geothermalflux[indices[i]]; - } - - sub_friction_type=this->friction_type; - sub_p=this->p; - sub_q=this->q; - sub_shelf=this->shelf; - sub_onbed=this->onbed; - sub_onwater=this->onwater; - sub_onsurface=this->onsurface; - sub_collapse=this->collapse; - sub_thermal_steadystate=this->thermal_steadystate; - - sub_properties=new ElementProperties(numindices,sub_h, sub_s, sub_b, sub_k, sub_melting, sub_accumulation, sub_geothermalflux, sub_friction_type, sub_p, sub_q, sub_shelf, sub_onbed, sub_onwater, sub_onsurface, sub_collapse, sub_thermal_steadystate); - - - xfree((void**)&sub_h); - xfree((void**)&sub_s); - xfree((void**)&sub_b); - xfree((void**)&sub_k); - xfree((void**)&sub_melting); - xfree((void**)&sub_accumulation); - xfree((void**)&sub_geothermalflux); - - return sub_properties; - -} -/*}}}*/ Index: sm/trunk/src/c/objects/ElementProperties.h =================================================================== --- /issm/trunk/src/c/objects/ElementProperties.h (revision 3611) +++ (revision ) @@ -1,55 +1,0 @@ -/*!\file: ElementProperties.h - * \brief prototype for ElementProperties.h - */ - -#ifndef _ELEMENTPROPERTIES_H_ -#define _ELEMENTPROPERTIES_H_ - -/*indefinitions: */ -struct IoModel; - -class ElementProperties{ - - public: - - /*nodal property values: */ - int numnodes; - double* h; - double* s; - double* b; - double* k; - double* melting; - double* accumulation; - double* geothermalflux; - - /*element property values: */ - int friction_type; - double p; - double q; - int shelf; - int onbed; - bool onwater; - int onsurface; - int collapse; - int thermal_steadystate; - - ElementProperties(); - ElementProperties(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate); - ElementProperties(ElementProperties* properties); - ~ElementProperties(); - - void Init(int numnodes, int* vertex_ids, int element_id, IoModel* iomodel); - void Init(int numnodes, double* h, double* s, double* b, double* k, double* melting, double* accumulation, double* geothermalflux, int friction_type, double p, double q, int shelf, int onbed, bool onwater, int onsurface, int collapse, int thermal_steadystate); - - void Marshall(char** pmarshalled_dataset); - int MarshallSize(); - void Demarshall(char** pmarshalled_dataset); - void copy(ElementProperties* properties); - void Echo(void); - void DeepEcho(void); - ElementProperties* Spawn(int* indices, int numindices); - -}; -#endif //ifndef _ELEMENTPROPERTIES_H_ - - Index: /issm/trunk/src/c/objects/Friction.cpp =================================================================== --- /issm/trunk/src/c/objects/Friction.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Friction.cpp (revision 3612) @@ -45,9 +45,16 @@ printf("\n"); - printf("velocities: "); - for(i=0;i<3;i++)printf("%g ",friction->velocities[i*2+0]); + printf("vx: "); + for(i=0;i<3;i++)printf("%g ",friction->vx[i]); printf("\n "); - for(i=0;i<3;i++)printf("%g ",friction->velocities[i*2+1]); + printf("vy: "); + for(i=0;i<3;i++)printf("%g ",friction->vy[i]); printf("\n"); + if(friction->vz){ + printf("vz: "); + for(i=0;i<3;i++)printf("%g ",friction->vz[i]); + printf("\n"); + } + } /*}}}*/ @@ -66,5 +73,7 @@ friction->bed=NULL; friction->thickness=NULL; - friction->velocities=NULL; + friction->vx=NULL; + friction->vy=NULL; + friction->vz=NULL; friction->p=UNDEF; friction->q=UNDEF; @@ -127,12 +136,12 @@ //We need the velocity magnitude to evaluate the basal stress: if(strcmp(friction->element_type,"2d")){ - velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2 - velocity_y[i]=*(friction->velocities+2*i+1); + velocity_x[i]=friction->vx[i];//velocities of size numgridsx2 + velocity_y[i]=friction->vy[i]; velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)); } else{ - velocity_x[i]=*(friction->velocities+3*i+0); //velocities of size numgridsx3 - velocity_y[i]=*(friction->velocities+3*i+1); - velocity_z[i]=*(friction->velocities+3*i+2); + velocity_x[i]=friction->vx[i]; //velocities of size numgridsx3 + velocity_y[i]=friction->vy[i]; + velocity_z[i]=friction->vz[i]; velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)+pow(velocity_z[i],2)); } @@ -177,6 +186,6 @@ //We need the velocity magnitude to evaluate the basal stress: - velocity_x[i]=*(friction->velocities+2*i+0); //velocities of size numgridsx2 - velocity_y[i]=*(friction->velocities+2*i+1); + velocity_x[i]=friction->vx[i]; //velocities of size numgridsx2 + velocity_y[i]=friction->vy[i]; velocity_mag[i]=sqrt(pow(velocity_x[i],2)+pow(velocity_y[i],2)); Index: /issm/trunk/src/c/objects/Friction.h =================================================================== --- /issm/trunk/src/c/objects/Friction.h (revision 3611) +++ /issm/trunk/src/c/objects/Friction.h (revision 3612) @@ -16,5 +16,7 @@ double* bed; double* thickness; - double* velocities; + double* vx; + double* vy; + double* vz; double p; double q; Index: /issm/trunk/src/c/objects/Hook.cpp =================================================================== --- /issm/trunk/src/c/objects/Hook.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Hook.cpp (revision 3612) @@ -16,5 +16,4 @@ #include "./Hook.h" #include "../EnumDefinitions/EnumDefinitions.h" -#include "./ParameterInputs.h" #include "../shared/shared.h" #include "../include/typedefs.h" Index: sm/trunk/src/c/objects/Input.cpp =================================================================== --- /issm/trunk/src/c/objects/Input.cpp (revision 3611) +++ (revision ) @@ -1,302 +1,0 @@ -/*!\file Input.c - * \brief: implementation of the Input object - */ - - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "stdio.h" -#include -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../shared/shared.h" -#include "../include/typedefs.h" -#include "../include/macros.h" -#include "./objects.h" - - -/*Object constructors and destructor*/ -/*FUNCTION Input::constructor {{{1*/ -Input::Input(){ - return; -} -/*}}}*/ -/*FUNCTION Input::destructor {{{1*/ -Input::~Input(){ - - xfree((void**)&vector); - return; -} -/*}}}*/ -/*FUNCTION Input::creation(char* input_name,int input_integer) {{{1*/ -Input::Input(char* input_name,int input_integer){ - - strcpy(name,input_name); - type=INTEGER; - integer=input_integer; - - ndof=0; - numberofnodes=0; - vector=NULL; -} -/*}}}*/ -/*FUNCTION Input::creation(char* input_name,double input_scalar) {{{1*/ -Input::Input(char* input_name,double input_scalar){ - - strcpy(name,input_name); - type=DOUBLE; - scalar=input_scalar; - - ndof=0; - numberofnodes=0; - vector=NULL; - -} -/*}}}*/ -/*FUNCTION Input::creation(char* input_name,char* input_string) {{{1*/ -Input::Input(char* input_name,char* input_string){ - - strcpy(name,input_name); - type=STRING; - strcpy(string,input_string); - ndof=0; - numberofnodes=0; - vector=NULL; -} -/*}}}*/ -/*FUNCTION Input::creation(char* input_name,double* input_vector,int input_ndof,int input_numberofnodes) {{{1*/ -Input::Input(char* input_name,double* input_vector,int input_ndof,int input_numberofnodes){ - - strcpy(name,input_name); - type=DOUBLEVEC; - vector=(double*)xmalloc(input_ndof*input_numberofnodes*sizeof(double)); - memcpy(vector,input_vector,input_ndof*input_numberofnodes*sizeof(double)); - ndof=input_ndof; - numberofnodes=input_numberofnodes; - -} -/*}}}*/ -/*FUNCTION Input::creation(char* input_name,Vec input_vector,int input_ndof, int input_numberofnode) {{{1*/ -Input::Input(char* input_name,Vec input_vector,int input_ndof, int input_numberofnodes){ - - strcpy(name,input_name); - type=DOUBLEVEC; - VecToMPISerial(&vector,input_vector); - ndof=input_ndof; - numberofnodes=input_numberofnodes; - -} -/*}}}*/ - -/*Object marshall*/ -/*FUNCTION Input::Demarshall {{{1*/ -void Input::Demarshall(char** pmarshalled_dataset){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ -/*FUNCTION Input::Marshall {{{1*/ -void Input::Marshall(char** pmarshalled_dataset){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ -/*FUNCTION Input::MarshallSize {{{1*/ -int Input::MarshallSize(){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ - -/*Object functions*/ -/*FUNCTION Input::copy {{{1*/ -Object* Input::copy() { - - return new Input(*this); - -} -/*}}}*/ -/*FUNCTION Input::DeepEcho {{{1*/ -void Input::DeepEcho(void){ - - int i,j; - - printf("Input:\n"); - printf(" name: %s\n",name); - printf(" type: %i\n",type); - if(type==STRING) printf(" string: %s\n",string); - if(type==INTEGER) printf(" integer: %i\n",integer); - if(type==DOUBLE) printf(" double: %g\n",scalar); - if(type==DOUBLEVEC) printf(" doublevec:\n"); - for(i=0;iRecover(&values[0],ndof,dofs,3,this->nodes); //ie: recover values for dofs 0 and 2 of all nodes of the tria - - int* doflist=NULL; - int dof1; - int dof; - int i,j; - - Node** nodes=NULL; - Node* node=NULL; - - /*recover pointer to nodes: */ - nodes=(Node**)vpnodesin; - - if (!valuesin) ISSMERROR("valuesin has not been allocated properly (pointer is NULL)"); - if (type!=DOUBLEVEC) ISSMERROR("%s%i%s"," cannot recover values from a ",type," input type"); - - /*Ok, we are trying to fill values. values is of size ndof*numnodes. The dofs - * for picking up the values are found in the numnodes nodes. We firt build the dof list: */ - doflist=(int*)xmalloc(ndofin*numnodesin*sizeof(int)); - - for(i=0;iGetDofList1(); - - for(j=0;j +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../include/typedefs.h" +#include "../include/types.h" +#include "../include/macros.h" + +/*Object constructors and destructor*/ +/*FUNCTION IntInput::IntInput(){{{1*/ +IntInput::IntInput(){ + return; +} +/*}}}*/ +/*FUNCTION IntInput::IntInput(double* values){{{1*/ +IntInput::IntInput(int in_enum_type,IssmInt in_value){ + + enum_type=in_enum_type; + value=in_value; +} +/*}}}*/ +/*FUNCTION IntInput::~IntInput(){{{1*/ +IntInput::~IntInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION IntInput::copy{{{1*/ +Object* IntInput::copy() { + + return new IntInput(this->enum_type,this->value); + +} +/*}}}*/ +/*FUNCTION IntInput::DeepEcho{{{1*/ +void IntInput::DeepEcho(void){ + + printf("IntInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %i\n",this->value); +} +/*}}}*/ +/*FUNCTION IntInput::Demarshall{{{1*/ +void IntInput::Demarshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int i; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*this time, no need to get enum type, the pointer directly points to the beginning of the + *object data (thanks to DataSet::Demarshall):*/ + memcpy(&enum_type,marshalled_dataset,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(&value,marshalled_dataset,sizeof(value));marshalled_dataset+=sizeof(value); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION IntInput::Echo {{{1*/ +void IntInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION IntInput::Enum{{{1*/ +int IntInput::Enum(void){ + + return IntInputEnum; + +} +/*}}}*/ +/*FUNCTION IntInput::EnumType{{{1*/ +int IntInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION IntInput::GetId{{{1*/ +int IntInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION IntInput::GetName{{{1*/ +char* IntInput::GetName(void){ + return "intinput"; +} +/*}}}*/ +/*FUNCTION IntInput::Marshall{{{1*/ +void IntInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of IntInput: */ + enum_type=IntInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall IntInput data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&value,sizeof(value));marshalled_dataset+=sizeof(value); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION IntInput::MarshallSize{{{1*/ +int IntInput::MarshallSize(){ + + return sizeof(value)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION IntInput::MyRank{{{1*/ +int IntInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/IntInput.h =================================================================== --- /issm/trunk/src/c/objects/IntInput.h (revision 3612) +++ /issm/trunk/src/c/objects/IntInput.h (revision 3612) @@ -0,0 +1,43 @@ +/*! \file IntInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" +#include "../include/types.h" + +#ifndef _INTINPUT_H_ +#define _INTINPUT_H_ + +class IntInput: public Input{ + + private: + /*just hold 3 values for 3 vertices: */ + int enum_type; + IssmInt value; + + public: + + /*constructors, destructors: {{{1*/ + IntInput(); + IntInput(int enum_type,IssmInt value); + ~IntInput(); + /*}}}*/ + /*object management: {{{1*/ + void DeepEcho(); + void Echo(); + int GetId(); + void Marshall(char** pmarshalled_dataset); + int MarshallSize(); + char* GetName(); + void Demarshall(char** pmarshalled_dataset); + int Enum(); + int MyRank(); + Object* copy(); + int EnumType(); + + /*}}}*/ + /*numerics: {{{1*/ + /*}}}*/ + +}; +#endif /* _INTINPUT_H */ Index: /issm/trunk/src/c/objects/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Matice.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Matice.cpp (revision 3612) @@ -53,8 +53,8 @@ /*Compute B on the element if provided*/ if (num_vertices==3 || num_vertices==6){ - if (iomodel->B){ + if (iomodel->rheology_B){ B_avg=0; for(j=0;jB+((int)*(iomodel->elements+num_vertices*i+j)-1)); + B_avg+=*(iomodel->rheology_B+((int)*(iomodel->elements+num_vertices*i+j)-1)); } B_avg=B_avg/num_vertices; @@ -62,19 +62,19 @@ } else if (num_vertices==1 || num_vertices==2){ - if (iomodel->B){ - B_avg=*(iomodel->B+i); + if (iomodel->rheology_B){ + B_avg=*(iomodel->rheology_B+i); } } else ISSMERROR("num_vertices = %i not supported yet",num_vertices); - if (iomodel->B) matice_B=B_avg; + if (iomodel->rheology_B) matice_B=B_avg; else matice_B=UNDEF; - if (iomodel->n){ + if (iomodel->rheology_n){ if (num_vertices==3 || num_vertices==6){ - matice_n=(double)*(iomodel->n+i); + matice_n=(double)*(iomodel->rheology_n+i); } else if (num_vertices==1 || num_vertices==2){ /*n is on the elements for now, so just take the first one*/ - matice_n=(double)*(iomodel->n); + matice_n=(double)*(iomodel->rheology_n); } else ISSMERROR("num_vertices = %i not supported yet",num_vertices); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3612) @@ -17,5 +17,4 @@ #include #include "../EnumDefinitions/EnumDefinitions.h" -#include "./ParameterInputs.h" #include "../shared/shared.h" #include "../include/typedefs.h" Index: /issm/trunk/src/c/objects/Numpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Numpar.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Numpar.cpp (revision 3612) @@ -39,4 +39,5 @@ cm_maxdmp_value=UNDEF; cm_maxdmp_slope=UNDEF; + dt=UNDEF; return; @@ -82,4 +83,5 @@ memcpy(marshalled_dataset,&cm_maxdmp_value,sizeof(cm_maxdmp_value));marshalled_dataset+=sizeof(cm_maxdmp_value); memcpy(marshalled_dataset,&cm_maxdmp_slope,sizeof(cm_maxdmp_slope));marshalled_dataset+=sizeof(cm_maxdmp_slope); + memcpy(marshalled_dataset,&dt,sizeof(dt));marshalled_dataset+=sizeof(dt); *pmarshalled_dataset=marshalled_dataset; @@ -101,4 +103,5 @@ +sizeof(cm_maxdmp_value) +sizeof(cm_maxdmp_slope) + +sizeof(dt) +sizeof(int); //sizeof(int) for enum type } @@ -128,4 +131,5 @@ memcpy(&cm_maxdmp_value,marshalled_dataset,sizeof(cm_maxdmp_value));marshalled_dataset+=sizeof(cm_maxdmp_value); memcpy(&cm_maxdmp_slope,marshalled_dataset,sizeof(cm_maxdmp_slope));marshalled_dataset+=sizeof(cm_maxdmp_slope); + memcpy(&dt,marshalled_dataset,sizeof(dt));marshalled_dataset+=sizeof(dt); /*return: */ @@ -158,4 +162,5 @@ parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); + parameters->FindParam(&dt,"dt"); return; @@ -185,4 +190,5 @@ printf(" cm_maxdmp_value: %g\n",cm_maxdmp_value); printf(" cm_maxdmp_slope: %g\n",cm_maxdmp_slope); + printf(" dt: %g\n",dt); } /*}}}*/ @@ -203,4 +209,5 @@ printf(" cm_maxdmp_value: %g\n",cm_maxdmp_value); printf(" cm_maxdmp_slope: %g\n",cm_maxdmp_slope); + printf(" dt: %g\n",dt); } /*}}}*/ @@ -246,5 +253,6 @@ inputs->Recover("cm_maxdmp_value",&cm_maxdmp_value); inputs->Recover("cm_maxdmp_slope",&cm_maxdmp_slope); - -} -/*}}}*/ + inputs->Recover("dt",&dt); + +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Numpar.h =================================================================== --- /issm/trunk/src/c/objects/Numpar.h (revision 3611) +++ /issm/trunk/src/c/objects/Numpar.h (revision 3612) @@ -27,4 +27,5 @@ double cm_maxdmp_value; double cm_maxdmp_slope; + double dt; Numpar(); Index: /issm/trunk/src/c/objects/OptArgs.h =================================================================== --- /issm/trunk/src/c/objects/OptArgs.h (revision 3611) +++ /issm/trunk/src/c/objects/OptArgs.h (revision 3612) @@ -24,12 +24,8 @@ #else -#include "./ParameterInputs.h" -class ParameterInputs; - struct OptArgs{ Model* model; double* param_g; double* grad_g; - ParameterInputs* inputs; int n; }; Index: sm/trunk/src/c/objects/ParameterInputs.cpp =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.cpp (revision 3611) +++ (revision ) @@ -1,386 +1,0 @@ -/*!\file ParameterInputs.c - * \brief: implementation of the ParameterInputs object - */ - - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "stdio.h" -#include -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../shared/shared.h" -#include "../include/typedefs.h" -#include "../include/macros.h" -#include "./objects.h" - -/*Object constructors and destructor*/ -/*FUNCTION ParameterInputs::constructor {{{1*/ -ParameterInputs::ParameterInputs(){ - - dataset=new DataSet(); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::destructor {{{1*/ -ParameterInputs::~ParameterInputs(){ - - delete dataset; - -} -/*}}}*/ -/*FUNCTION ParameterInputs::purge {{{1*/ -void ParameterInputs::purge(char* name){ - - int i; - Input* input=NULL; - - /*Go through dataset, and figure out if an Input - * already has the name "name". If so, delete it : */ - - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - /*delete object: */ - dataset->DeleteObject(input); - } - } -} -/*}}}*/ - -/*Object functions*/ -/*FUNCTION ParameterInputs::Add(char* name,char* string) {{{1*/ -void ParameterInputs::Add(char* name,char* string){ - - Input* input=NULL; - - /*First, purge object with same name: */ - this->purge(name); - - /*We are sure an input of the same name does not exist. Create new - * input: */ - input=new Input(name,string); - - /*Add input to dataset: */ - dataset->AddObject(input); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Add(char* name,int integer) {{{1*/ -void ParameterInputs::Add(char* name,int integer){ - - Input* input=NULL; - - /*First, purge object with same name: */ - this->purge(name); - - /*We are sure an input of the same name does not exist. Create new - * input: */ - input=new Input(name,integer); - - /*Add input to dataset: */ - dataset->AddObject(input); -} -/*}}}*/ -/*FUNCTION ParameterInputs::Add(char* name,double scalar) {{{1*/ -void ParameterInputs::Add(char* name,double scalar){ - - Input* input=NULL; - - /*First, purge object with same name: */ - this->purge(name); - - /*We are sure an input of the same name does not exist. Create new - * input: */ - input=new Input(name,scalar); - - /*Add input to dataset: */ - dataset->AddObject(input); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes) {{{1*/ -void ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes){ - - Input* input=NULL; - - /*First, purge object with same name: */ - this->purge(name); - - /*We are sure an input of the same name does not exist. Create new - * input: */ - input=new Input(name,vector,ndof,numberofnodes); - - /*Add input to dataset: */ - dataset->AddObject(input); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes) {{{1*/ -void ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes){ - - Input* input=NULL; - - /*First, purge object with same name: */ - this->purge(name); - - /*We are sure an input of the same name does not exist. Create new - * input: */ - input=new Input(name,petscvector,ndof,numberofnodes); - - /*Add input to dataset: */ - dataset->AddObject(input); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::DeepEcho {{{1*/ -void ParameterInputs::DeepEcho(){ - - printf("ParameterInputs echo: \n"); - dataset->DeepEcho(); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Echo {{{1*/ -void ParameterInputs::Echo(){ - - printf("ParameterInputs echo: \n"); - dataset->Echo(); - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Get {{{1*/ -Vec ParameterInputs::Get(char* name,int* dofs, int numdofs){ - - int i; - Input* input=NULL; - int found=0; - - /*output: */ - Vec ug=NULL; - - /*Go through dataset, and figure out if an Input - * has the name "name": */ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - found=1; - break; - } - } - - if(found==0)return NULL; - - /*call submethod: */ - ug=input->Get(dofs,numdofs); - - return ug; -} -/*}}}*/ -/*FUNCTION ParameterInputs::IsPresent(char* name,char* string) {{{1*/ -int ParameterInputs::IsPresent(char* name){ - - /*Intermediary*/ - int i; - Input* input=NULL; - - /*Go through dataset, and figure out if the input is present*/ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - if (input->IsPresent(name)) return 1; - } - - /*Input not found... return 0*/ - return 0; -} -/*}}}*/ -/*FUNCTION ParameterInputs::Recover(char* name, char** pstring) {{{1*/ -int ParameterInputs::Recover(char* name, char** pstring){ - - int found=0; - int i; - Input* input=NULL; - char* string=NULL; - - /*Go through dataset, and figure out if an Input - * has the name "name": */ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - - /*Get string out of this input: */ - input->Recover(&string); - found=1;break; - } - } - /*Assign output pointers:*/ - *pstring=string; - return found; - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Recover(char* name, int* pinteger) {{{1*/ -int ParameterInputs::Recover(char* name, int* pinteger){ - - int found=0; - int i; - Input* input=NULL; - int integer; - - /*Go through dataset, and figure out if an Input - * has the name "name": */ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - - /*Get string out of this input: */ - input->Recover(&integer); - found=1; break; - - } - } - /*Assign output pointers:*/ - *pinteger=integer; - - return found; -} -/*}}}*/ -/*FUNCTION ParameterInputs::Recover(char* name, double* pscalar) {{{1*/ -int ParameterInputs::Recover(char* name, double* pscalar){ - - int found=0; - int i; - Input* input=NULL; - double scalar; - - /*Go through dataset, and figure out if an Input - * has the name "name": */ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - - /*Get string out of this input: */ - input->Recover(&scalar); - found=1; break; - - } - } - /*Assign output pointers:*/ - *pscalar=scalar; - - return found; - -} -/*}}}*/ -/*FUNCTION ParameterInputs::Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** pnodes) {{{1*/ -int ParameterInputs::Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** pnodes){ - - int i; - Node** nodes=NULL; - Input* input=NULL; - int found=0; - - /*First recover nodes pointer :*/ - nodes=(Node**)pnodes; - - /*Go through dataset, and figure out if an Input - * has the name "name": */ - for(i=0;iSize();i++){ - input=(Input*)dataset->GetObjectByOffset(i); - - if (input->IsSameName(name)){ - found=1; - break; - } - } - - if(found==0)return 0; - - /*call submethod: */ - input->Recover(values,ndof,dofs,numnodes,(void**)nodes); - return 1; -} -/*}}}*/ -/*FUNCTION ParameterInputs::UpdateFromDakota {{{1*/ -void ParameterInputs::UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters,double* qmu_part,int qmu_npart){ - - /*Go through all dakota descriptors, ex: "rho_ice","thermal_conductivity","thickness1","thickness2", etc ..., and - * for each descriptor, take the variable value and plug it into the inputs: */ - - int i,j,k; - int found; - - double* distributed_values=NULL; - char* root=NULL; //root name of distributed variable, ex: thickness, drag, etc ... - double* parameter=NULL; - int numberofnodes; - - char* descriptor=NULL; - - /*retrive some necessary parameters first: */ - found=parameters->FindParam(&numberofnodes,"numberofnodes"); - if(!found)ISSMERROR("coud not find numberofnodes in parameters dataset!"); - - for(i=0;iAdd(descriptor,variables[i]); - } - else{ - - /*Ok, variable is distributed. Root name of variable is also known. Now, allocate distributed_values and fill the - * distributed_values with the next qmu_npart variables: */ - - distributed_values=(double*)xmalloc(qmu_npart*sizeof(double)); - for(j=0;jFindParam(¶meter,NULL,NULL,root); - if(!found)ISSMERROR("%s%s%s"," could not find parameter ",root," for Dakota input update"); - - /*We've got the parameter, we need to update it using qmu_part (a partitioning vector), and the distributed_values: */ - for(k=0;kAdd(root,parameter,1,numberofnodes); - - /*increment i to skip the distributed values just collected: */ - i+=qmu_npart-1; //careful, the for loop will add 1. - - /*Free allocations: */ - xfree((void**)¶meter); - xfree((void**)&distributed_values); - xfree((void**)&root); - - } - } - - ISSMERROR("create UpdateFromDakota routine, to call UpdateFromDakota in the elements, so it does not conflict with UpdateFromInputs! in the elements"); - - - /*Free ressources:*/ - xfree((void**)&distributed_values); - xfree((void**)&root); - xfree((void**)¶meter); - -} -/*}}}*/ Index: sm/trunk/src/c/objects/ParameterInputs.h =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.h (revision 3611) +++ (revision ) @@ -1,46 +1,0 @@ -/*!\file: ParameterInputs.h - * \brief dataset filled with input objects. Used by CreateKMatrix and CreatePVector to hold - * parameters coming from the workspace. - */ - -#ifndef _PARAMETERINPUTS_H_ -#define _PARAMETERINPUTS_H_ - -#include "../toolkits/toolkits.h" -#include "../DataSet/DataSet.h" - -class DataSet; - -class ParameterInputs{ - - private: - DataSet* dataset; //cannot declare DataSet* dataset, to avoid circular referencing. - - public: - ParameterInputs(); - ~ParameterInputs(); - void Echo(); - void DeepEcho(); - - void Add(char* name,char* string); - void Add(char* name,int integer); - void Add(char* name,double scalar); - void Add(char* name,double* vector,int ndof,int numberofnodes); - void Add(char* name,Vec petscvector,int ndof, int numberofnodes); - - int Recover(char* name, char** pstring); - int Recover(char* name, int* pinteger); - int Recover(char* name, double* pscalar); - int Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** nodes); //void** because otherwise we get a circular dependency with Node - - Vec Get(char* field_name,int* dofs, int numdofs); - int IsPresent(char* name); - - /*declaration found in io: */ - void Init( void* data_handle); - void purge(char* name); - void UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters, double* qmu_part,int qmu_npart); - -}; - -#endif Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3612) @@ -266,68 +266,4 @@ } /*}}}*/ -/*FUNCTION UpdateFromInputs {{{1*/ -void Penta::UpdateFromInputs(void* vinputs){ - - int dofs[1]={0}; - double temperature_list[6]; - double temperature_average; - double B_list[6]; - double B_average; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matice* matice=NULL; - - ParameterInputs* inputs=NULL; - - /*If on water, skip: */ - if(this->properties.onwater)return; - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matice=(Matice*)hmatice.delivers(); - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*Update internal data if inputs holds new values: */ - if (this->properties.h) inputs->Recover("thickness",&this->properties.h[0],1,dofs,6,(void**)nodes); - if (this->properties.s) inputs->Recover("surface",&this->properties.s[0],1,dofs,6,(void**)nodes); - if (this->properties.b) inputs->Recover("bed",&this->properties.b[0],1,dofs,6,(void**)nodes); - if (this->properties.k) inputs->Recover("drag",&this->properties.k[0],1,dofs,6,(void**)nodes); - if (this->properties.melting) inputs->Recover("melting",&this->properties.melting[0],1,dofs,6,(void**)nodes); - if (this->properties.accumulation) inputs->Recover("accumulation",&this->properties.accumulation[0],1,dofs,6,(void**)nodes); - if (this->properties.geothermalflux) inputs->Recover("geothermalflux",&this->properties.geothermalflux[0],1,dofs,6,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature",&temperature_list[0],1,dofs,6,(void**)nodes)){ - if(matice && !this->properties.collapse){ - //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) - // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - } - - if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ - if(matice && this->properties.collapse){ - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; - B_average=Paterson(temperature_average); - //B_average=(Paterson(temperature_list[0])+Paterson(temperature_list[1])+Paterson(temperature_list[2]) - // +Paterson(temperature_list[3])+Paterson(temperature_list[4])+Paterson(temperature_list[5]))/6.0; - matice->SetB(B_average); - } - } - - if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ - if(matice){ - B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0; - matice->SetB(B_average); - } - } - -} -/*}}}*/ /*FUNCTION UpdateFromDakota {{{1*/ void Penta::UpdateFromDakota(void* vinputs){ @@ -345,5 +281,5 @@ /*Object functions*/ /*FUNCTION ComputeBasalStress {{{1*/ -void Penta::ComputeBasalStress(Vec sigma_b,void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type){ int i,j; @@ -477,5 +413,5 @@ /*}}}*/ /*FUNCTION ComputePressure {{{1*/ -void Penta::ComputePressure(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::ComputePressure(Vec pg,int analysis_type,int sub_analysis_type){ int i; @@ -519,5 +455,5 @@ /*}}}*/ /*FUNCTION ComputeStrainRate {{{1*/ -void Penta::ComputeStrainRate(Vec eps,void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type){ ISSMERROR("Not implemented yet"); @@ -526,5 +462,5 @@ /*}}}*/ /*FUNCTION CostFunction {{{1*/ -double Penta::CostFunction(void* inputs,int analysis_type,int sub_analysis_type){ +double Penta::CostFunction(int analysis_type,int sub_analysis_type){ double J; @@ -560,5 +496,5 @@ /*FUNCTION CreateKMatrix {{{1*/ -void Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -608,5 +544,5 @@ /*}}}*/ /*FUNCTION CreateKMatrixDiagnosticHoriz {{{1*/ -void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, int analysis_type,int sub_analysis_type){ @@ -844,5 +780,5 @@ /*}}}*/ /*FUNCTION CreateKMatrixDiagnosticStokes {{{1*/ -void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixDiagnosticStokes( Mat Kgg, int analysis_type,int sub_analysis_type){ int i,j; @@ -1139,5 +1075,5 @@ /*}}}*/ /*FUNCTION CreateKMatrixDiagnosticVert {{{1*/ -void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -1275,5 +1211,5 @@ /*}}}*/ /*FUNCTION CreateKMatrixMelting {{{1*/ -void Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixMelting(Mat Kgg,int analysis_type,int sub_analysis_type){ Tria* tria=NULL; @@ -1296,5 +1232,5 @@ /*FUNCTION CreateKMatrixPrognostic {{{1*/ -void Penta::CreateKMatrixPrognostic(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -1317,5 +1253,5 @@ /*FUNCTION CreateKMatrixSlopeCompute {{{1*/ -void Penta::CreateKMatrixSlopeCompute(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixSlopeCompute(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -1337,5 +1273,5 @@ /*}}}*/ /*FUNCTION CreateKMatrixThermal {{{1*/ -void Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreateKMatrixThermal(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -1600,5 +1536,5 @@ /*}}}*/ /*FUNCTION CreatePVector {{{1*/ -void Penta::CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Penta::CreatePVector(Vec pg, int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -1646,5 +1582,5 @@ /*}}}*/ /*FUNCTION CreatePVectorDiagnosticHoriz {{{1*/ -void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorDiagnosticHoriz( Vec pg, int analysis_type,int sub_analysis_type){ int i,j; @@ -1801,5 +1737,5 @@ /*}}}*/ /*FUNCTION CreatePVectorDiagnosticStokes {{{1*/ -void Penta::CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorDiagnosticStokes( Vec pg, int analysis_type,int sub_analysis_type){ /*indexing: */ @@ -2051,5 +1987,5 @@ /*}}}*/ /*FUNCTION CreatePVectorDiagnosticVert {{{1*/ -void Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorDiagnosticVert( Vec pg, int analysis_type,int sub_analysis_type){ int i; @@ -2190,5 +2126,5 @@ /*}}}*/ /*FUNCTION CreatePVectorMelting {{{1*/ -void Penta::CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorMelting( Vec pg, int analysis_type,int sub_analysis_type){ return; } @@ -2196,5 +2132,5 @@ /*FUNCTION CreatePVectorPrognostic {{{1*/ -void Penta::CreatePVectorPrognostic( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorPrognostic( Vec pg, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2216,5 +2152,5 @@ /*FUNCTION CreatePVectorSlopeCompute {{{1*/ -void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorSlopeCompute( Vec pg, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -2235,5 +2171,5 @@ /*}}}*/ /*FUNCTION CreatePVectorThermal {{{1*/ -void Penta::CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type){ +void Penta::CreatePVectorThermal( Vec pg, int analysis_type,int sub_analysis_type){ @@ -2438,5 +2374,5 @@ /*}}}*/ /*FUNCTION Du {{{1*/ -void Penta::Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::Du(Vec du_g,int analysis_type,int sub_analysis_type){ int i; @@ -3647,5 +3583,5 @@ /*}}}*/ /*FUNCTION GetOnBed {{{1*/ -int Penta::GetOnBed(){ +bool Penta::GetOnBed(){ return this->properties.onbed; } @@ -3792,5 +3728,5 @@ /*}}}*/ /*FUNCTION Gradj {{{1*/ -void Penta::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ +void Penta::Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type){ /*If on water, skip grad (=0): */ @@ -3807,5 +3743,5 @@ /*}}}*/ /*FUNCTION GradjDrag {{{1*/ -void Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::GradjDrag(Vec grad_g,int analysis_type,int sub_analysis_type){ Tria* tria=NULL; @@ -3840,5 +3776,5 @@ /*}}}*/ /*FUNCTION GradjB {{{1*/ -void Penta::GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){ +void Penta::GradjB(Vec grad_g,int analysis_type,int sub_analysis_type){ Tria* tria=NULL; @@ -3873,5 +3809,5 @@ /*}}}*/ /*FUNCTION Misfit {{{1*/ -double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){ +double Penta::Misfit(int analysis_type,int sub_analysis_type){ double J; @@ -4006,5 +3942,5 @@ /*}}}1*/ /*FUNCTION SurfaceArea {{{1*/ -double Penta::SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type){ +double Penta::SurfaceArea(int analysis_type,int sub_analysis_type){ double S; Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3611) +++ /issm/trunk/src/c/objects/Penta.h (revision 3612) @@ -9,5 +9,4 @@ class Node; class Hook; -class ElementProperties; class DataSet; struct IoModel; @@ -21,7 +20,5 @@ #include "./Tria.h" #include "./Hook.h" -#include "./ElementProperties.h" #include "../ModelProcessorx/IoModel.h" -#include "./ParameterInputs.h" #include "./Node.h" @@ -36,6 +33,5 @@ Hook hnumpar; //hook to 1 numpar - ElementProperties properties; - DataSet* inputs; + Inputs* inputs; public: @@ -43,6 +39,6 @@ /*FUNCTION constructors, destructors {{{1*/ Penta(); - Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* properties); - Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties); + Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id); + Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* inputs); Penta(int i, IoModel* iomodel); ~Penta(); @@ -62,5 +58,4 @@ void* SpawnTria(int g0, int g1, int g2); void UpdateFromDakota(void* inputs); - void UpdateFromInputs(void* inputs); void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); void SetClone(int* minranks); @@ -68,21 +63,21 @@ /*}}}*/ /*FUNCTION element numerical routines {{{1*/ - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type,int sub_analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); + void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHoriz( Mat Kgg, int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticVert( Mat Kgg, int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); void* GetMatPar(); - int GetShelf(); + bool GetShelf(); void GetNodes(void** nodes); - int GetOnBed(); - void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type); - void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); - void GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - double Misfit(void* inputs,int analysis_type,int sub_analysis_type); - double SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type); - double CostFunction(void* inputs,int analysis_type,int sub_analysis_type); + bool GetOnBed(); + void Du(Vec du_g,int analysis_type,int sub_analysis_type); + void Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type); + void GradjDrag(Vec grad_g,int analysis_type,int sub_analysis_type); + void GradjB(Vec grad_g,int analysis_type,int sub_analysis_type); + double Misfit(int analysis_type,int sub_analysis_type); + double SurfaceArea(int analysis_type,int sub_analysis_type); + double CostFunction(int analysis_type,int sub_analysis_type); void GetThicknessList(double* thickness_list); @@ -99,20 +94,20 @@ void GetNodalFunctionsDerivativesReference(double* dl1dl6,double* gauss_coord); void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_coord); - void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHoriz( Vec pg, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticVert( Vec pg, int analysis_type,int sub_analysis_type); void GetParameterValue(double* pvalue, double* v_list,double* gauss_coord); void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord); void GetNodalFunctions(double* l1l6, double* gauss_coord); void FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed); - void ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type); - void ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type); - void ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorPrognostic( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); + void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); + void ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type); + void ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type); + void CreateKMatrixSlopeCompute(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorSlopeCompute( Vec pg, int analysis_type,int sub_analysis_type); + void CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorPrognostic( Vec pg, int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreatePVectorDiagnosticStokes( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticStokes( Mat Kgg, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticStokes( Vec pg, int analysis_type,int sub_analysis_type); void ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp); void GetMatrixInvert(double* Ke_invert, double* Ke); @@ -127,12 +122,12 @@ void ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp); void GetNodalFunctionsStokes(double* l1l7, double* gauss_coord); - void CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixThermal(Mat Kgg,int analysis_type,int sub_analysis_type); void GetB_conduct(double* B_conduct, double* xyz_list, double* gauss_coord); void GetB_advec(double* B_advec, double* xyz_list, double* gauss_coord); void GetBprime_advec(double* Bprime_advec, double* xyz_list, double* gauss_coord); void GetB_artdiff(double* B_artdiff, double* xyz_list, double* gauss_coord); - void CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVectorThermal( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorMelting( Vec pg, void* vinputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixMelting(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorThermal( Vec pg, int analysis_type,int sub_analysis_type); + void CreatePVectorMelting( Vec pg, int analysis_type,int sub_analysis_type); void GetPhi(double* phi, double* epsilon, double viscosity); double MassFlux(double* segment,double* ug); Index: /issm/trunk/src/c/objects/Result.cpp =================================================================== --- /issm/trunk/src/c/objects/Result.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Result.cpp (revision 3612) @@ -15,5 +15,4 @@ #include "../EnumDefinitions/EnumDefinitions.h" #include "../include/macros.h" -#include "./ParameterInputs.h" #include "../shared/shared.h" #include "../include/typedefs.h" Index: /issm/trunk/src/c/objects/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Riftfront.h (revision 3611) +++ /issm/trunk/src/c/objects/Riftfront.h (revision 3612) @@ -10,5 +10,4 @@ #include "./Element.h" #include "./Node.h" -#include "./ParameterInputs.h" #define MAX_RIFTFRONT_GRIDS 2 //max number of grids on a rift flank, only 2 because 2d for now. Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 3612) @@ -226,11 +226,26 @@ } /*}}}*/ -/*FUNCTION Sing::UpdateFromInputs {{{1*/ -void Sing::UpdateFromInputs(void* vinputs){ - - int dofs[1]={0}; - double temperature; - double B; - double B_average; +/*FUNCTION Sing::UpdateInputs {{{1*/ +void Sing::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ + +/*Object functions*/ +/*FUNCTION Sing::ComputeBasalStress {{{1*/ +void Sing::ComputeBasalStress(Vec p_g,int analysis_type,int sub_analysis_type){ + + ISSMERROR("Not implemented yet"); + +} +/*}}}*/ +/*FUNCTION Sing::ComputePressure {{{1*/ +void Sing::ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type){ + + int i; + const int numgrids=1; + int doflist[numgrids]; + double pressure[numgrids]; + double rho_ice,g; /*dynamic objects pointed to by hooks: */ @@ -240,8 +255,6 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&doflist[0]); /*recover objects from hooks: */ @@ -251,79 +264,4 @@ numpar=(Numpar*)hnumpar.delivers(); - /*Update internal data if inputs holds new values: */ - //if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n"); - //For now - if(this->properties.h) inputs->Recover("thickness",&this->properties.h[0],1,dofs,1,(void**)nodes); - //Later - /* - if(inputs->Recover("thickness",&new_h[0],1,dofs,1,(void**)nodes)){ - //density, needed later: - double di=(this->matpar->GetRhoIce()/this->matpar->GetRhoWater()); - //Go through grids: - for (i=0;i<1;i++){ - if(nodes[i]->IsOnShelf()){ - this->b[i]=this->b[i]-di*(new_h[i]-h[i]); //hydrostatic equilibrium; - } - this->s[i]=this->b[i]+new_h[i]; - this->h[i]=new_h[i]; - } - } - */ - if (this->properties.k) inputs->Recover("drag",&this->properties.k[0],1,dofs,1,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature_average",&temperature,1,dofs,1,(void**)nodes)){ - B_average=Paterson(temperature); - matice->SetB(B_average); - } - else if(inputs->Recover("temperature",&temperature,1,dofs,1,(void**)nodes)){ - B=Paterson(temperature); - matice->SetB(B); - } - - if(inputs->Recover("B",&B,1,dofs,1,(void**)nodes)){ - matice->SetB(B); - } - -} -/*}}}*/ -/*FUNCTION Sing::UpdateInputs {{{1*/ -void Sing::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){ - ISSMERROR(" not supported yet!"); -} -/*}}}*/ - -/*Object functions*/ -/*FUNCTION Sing::ComputeBasalStress {{{1*/ -void Sing::ComputeBasalStress(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){ - - ISSMERROR("Not implemented yet"); - -} -/*}}}*/ -/*FUNCTION Sing::ComputePressure {{{1*/ -void Sing::ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){ - - int i; - const int numgrids=1; - int doflist[numgrids]; - double pressure[numgrids]; - double rho_ice,g; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matpar* matpar=NULL; - Matice* matice=NULL; - Numpar* numpar=NULL; - - /*Get dof list on which we will plug the pressure values: */ - GetDofList1(&doflist[0]); - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - /*pressure is lithostatic: */ rho_ice=matpar->GetRhoIce(); @@ -337,5 +275,5 @@ /*}}}*/ /*FUNCTION Sing::ComputeStrainRate {{{1*/ -void Sing::ComputeStrainRate(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type){ +void Sing::ComputeStrainRate(Vec p_g,int analysis_type,int sub_analysis_type){ ISSMERROR("Not implemented yet"); @@ -350,5 +288,5 @@ /*FUNCTION Sing::CreateKMatrix {{{1*/ -void Sing::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Sing::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -366,5 +304,5 @@ /*FUNCTION Sing::CreateKMatrixDiagnosticHutter {{{1*/ -void Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Sing::CreateKMatrixDiagnosticHutter(Mat Kgg,int analysis_type,int sub_analysis_type){ const int numgrids=1; @@ -387,5 +325,5 @@ /*}}}*/ /*FUNCTION Sing::CreatePVector {{{1*/ -void Sing::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ +void Sing::CreatePVector(Vec pg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ @@ -403,5 +341,5 @@ /*FUNCTION Sing::CreatePVectorDiagnosticHutter {{{1*/ -void Sing::CreatePVectorDiagnosticHutter( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Sing::CreatePVectorDiagnosticHutter( Vec pg, int analysis_type,int sub_analysis_type){ @@ -561,5 +499,5 @@ /*}}}*/ /*FUNCTION Sing::GetOnBed {{{1*/ -int Sing::GetOnBed(){ +bool Sing::GetOnBed(){ ISSMERROR(" not supported yet!"); } Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 3611) +++ /issm/trunk/src/c/objects/Sing.h (revision 3612) @@ -10,11 +10,8 @@ #include "./Matice.h" #include "./Matpar.h" -#include "./ParameterInputs.h" -#include "./ElementProperties.h" #include "../ModelProcessorx/IoModel.h" #include "./Hook.h" class Hook; -class ElementProperties; class Sing: public Element{ @@ -30,5 +27,5 @@ Hook hnumpar; //hook to 1 numpar - ElementProperties properties; + Inputs* inputs; public: @@ -36,6 +33,6 @@ /*constructors, destructors: {{{1*/ Sing(); - Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id, ElementProperties* sing_properties); - Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, ElementProperties* sing_properties); + Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id); + Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs); Sing(int i, IoModel* iomodel); ~Sing(); @@ -56,22 +53,21 @@ /*}}}*/ /*numerics: {{{1*/ - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); - void UpdateFromInputs(void* inputs); + void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void CreateKMatrixDiagnosticHutter(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHutter(Mat Kgg,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); - void CreatePVectorDiagnosticHutter(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHutter(Vec pg,int analysis_type,int sub_analysis_type); void* GetMatPar(); - void ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type); - void ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type); - void ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type); + void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); + void ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type); + void ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type); void GetNodes(void** vpnodes); /*}}}*/ /*not implemented: {{{1*/ - int GetShelf(); - int GetOnBed(); + bool GetShelf(); + bool GetOnBed(); void GetBedList(double*); void GetThicknessList(double* thickness_list); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3612) @@ -19,4 +19,5 @@ #include "../shared/shared.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Inputs.h" #include "../include/typedefs.h" #include "../include/macros.h" @@ -29,27 +30,25 @@ } /*}}}*/ -/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/ -Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id, ElementProperties* triaproperties): +/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id){{{1*/ +Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id): hnodes(tria_node_ids,3), hmatice(&tria_matice_id,1), hmatpar(&tria_matpar_id,1), - hnumpar(&tria_numpar_id,1), - properties(triaproperties) + hnumpar(&tria_numpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=tria_id; - this->inputs=new DataSet(); + this->inputs=new Inputs(); return; } /*}}}*/ -/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties,DataSet* tria_inputs) {{{1*/ -Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, ElementProperties* tria_properties,DataSet* tria_inputs): +/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, Inputs* tria_inputs) {{{1*/ +Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs): hnodes(tria_hnodes), hmatice(tria_hmatice), hmatpar(tria_hmatpar), - hnumpar(tria_hnumpar), - properties(tria_properties) + hnumpar(tria_hnumpar) { @@ -57,8 +56,8 @@ this->id=tria_id; if(tria_inputs){ - this->inputs=tria_inputs->Copy(); + this->inputs=(Inputs*)tria_inputs->Copy(); } else{ - this->inputs=new DataSet(); + this->inputs=new Inputs(); } @@ -67,15 +66,15 @@ /*}}}*/ /*FUNCTION Tria::Tria(int i, IoModel* iomodel){{{1*/ -Tria::Tria(int i, IoModel* iomodel){ //i is the element index - - int j; +Tria::Tria(int index, IoModel* iomodel){ //i is the element index + + int i,j; int tria_node_ids[3]; int tria_matice_id; int tria_matpar_id; int tria_numpar_id; + double nodeinputs[3]; /*id: */ - this->id=i+1; - this->inputs=new DataSet(); + this->id=index+1; /*hooks: */ @@ -83,15 +82,15 @@ if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum){ /*Discontinuous Galerkin*/ - tria_node_ids[0]=3*i+1; - tria_node_ids[1]=3*i+2; - tria_node_ids[2]=3*i+3; + tria_node_ids[0]=3*index+1; + tria_node_ids[1]=3*index+2; + tria_node_ids[2]=3*index+3; } else{ /*Continuous Galerkin*/ - for(j=0;j<3;j++){ - tria_node_ids[j]=(int)*(iomodel->elements+3*i+j); //ids for vertices are in the elements array from Matlab + for(i=0;i<3;i++){ + tria_node_ids[i]=(int)*(iomodel->elements+3*index+i); //ids for vertices are in the elements array from Matlab } } - tria_matice_id=i+1; //refers to the corresponding ice material object + tria_matice_id=index+1; //refers to the corresponding ice material object tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object tria_numpar_id=1; //refers to numerical parameters object @@ -101,5 +100,47 @@ this->hmatpar.Init(&tria_matpar_id,1); this->hnumpar.Init(&tria_numpar_id,1); - this->properties.Init(3,tria_node_ids,this->id,iomodel); + + //intialize inputs, and add as many inputs per element as requested: + this->inputs=new Inputs(); + + if (iomodel->thickness) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(DragCoefficientEnum,nodeinputs)); + + if (iomodel->drag_p) this->inputs->AddInput(new DoubleInput(DragPEnum,iomodel->drag_p[index])); + if (iomodel->drag_q) this->inputs->AddInput(new DoubleInput(DragQEnum,iomodel->drag_q[index])); + this->inputs->AddInput(new IntInput(DragTypeEnum,iomodel->drag_type)); + + } + if (iomodel->melting_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_node_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs)); + } + + if (iomodel->elementoniceshelf) this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,(IssmBool)iomodel->elementoniceshelf[index])); + if (iomodel->elementonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->elementonbed[index])); + if (iomodel->elementonwater) this->inputs->AddInput(new BoolInput(ElementOnWaterEnum,(IssmBool)iomodel->elementonwater[index])); + if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); + } @@ -141,5 +182,5 @@ Object* Tria::copy() { - return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties,this->inputs); + return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); } @@ -166,8 +207,6 @@ hnumpar.Demarshall(&marshalled_dataset); - /*demarshall properties: */ - properties.Demarshall(&marshalled_dataset); - - inputs=DataSetDemarshallRaw(&marshalled_dataset); + /*demarshall inputs: */ + inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); /*return: */ @@ -186,6 +225,5 @@ hmatpar.DeepEcho(); hnumpar.DeepEcho(); - properties.DeepEcho(); - printf("Element inputs\n"); + printf(" inputs\n"); inputs->DeepEcho(); @@ -203,6 +241,5 @@ hmatpar.Echo(); hnumpar.Echo(); - properties.Echo(); - printf("Element inputs\n"); + printf(" inputs\n"); inputs->Echo(); @@ -236,7 +273,4 @@ hnumpar.Marshall(&marshalled_dataset); - /*Marshall properties: */ - properties.Marshall(&marshalled_dataset); - /*Marshall inputs: */ marshalled_inputs_size=inputs->MarshallSize(); @@ -257,5 +291,4 @@ +hmatpar.MarshallSize() +hnumpar.MarshallSize() - +properties.MarshallSize() +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type @@ -264,6 +297,6 @@ /*Updates: */ -/*FUNCTION Tria::UpdateFromInputs {{{1*/ -void Tria::UpdateFromInputs(void* vinputs){ +/*FUNCTION Tria::UpdateFromDakota {{{1*/ +void Tria::UpdateFromDakota(void* vinputs){ int i; @@ -280,77 +313,4 @@ Matice* matice=NULL; Numpar* numpar=NULL; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*Update internal data if inputs holds new values: */ - //if (id==1) printf("WARNING if QMU: no hydrostatic equilibrium is applied here (conflict with prognostic, which does not have matpar)\n"); - //For now - if(this->properties.h) inputs->Recover("thickness",&this->properties.h[0],1,dofs,3,(void**)nodes); - //Later - /* - if(inputs->Recover("thickness",&new_h[0],1,dofs,3,(void**)nodes)){ - //density, needed later: - double di=(this->matpar->GetRhoIce()/this->matpar->GetRhoWater()); - //Go through grids: - for (i=0;i<3;i++){ - if(nodes[i]->IsOnShelf()){ - this->b[i]=this->b[i]-di*(new_h[i]-h[i]); //hydrostatic equilibrium; - } - this->s[i]=this->b[i]+new_h[i]; - this->h[i]=new_h[i]; - } - } - */ - if (this->properties.s) inputs->Recover("surface",&this->properties.s[0],1,dofs,3,(void**)nodes); - if (this->properties.b) inputs->Recover("bed",&this->properties.b[0],1,dofs,3,(void**)nodes); - if (this->properties.k) inputs->Recover("drag",&this->properties.k[0],1,dofs,3,(void**)nodes); - if (this->properties.melting) inputs->Recover("melting",&this->properties.melting[0],1,dofs,3,(void**)nodes); - if (this->properties.accumulation) inputs->Recover("accumulation",&this->properties.accumulation[0],1,dofs,3,(void**)nodes); - if (this->properties.geothermalflux) inputs->Recover("geothermalflux",&this->properties.geothermalflux[0],1,dofs,3,(void**)nodes); - - //Update material if necessary - if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ - temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0; - B_average=Paterson(temperature_average); - matice->SetB(B_average); - } - - if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ - B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; - matice->SetB(B_average); - } - -} -/*}}}*/ -/*FUNCTION Tria::UpdateFromDakota {{{1*/ -void Tria::UpdateFromDakota(void* vinputs){ - - int i; - int dofs[1]={0}; - double temperature_list[3]; - double temperature_average; - double B_list[3]; - double B_average; - double new_h[3]; - - /*dynamic objects pointed to by hooks: */ - Node** nodes=NULL; - Matpar* matpar=NULL; - Matice* matice=NULL; - Numpar* numpar=NULL; - - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; /*recover objects from hooks: */ @@ -458,6 +418,6 @@ /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddEinput(new TriaVertexInput(VxEnum,vx)); - this->inputs->AddEinput(new TriaVertexInput(VyEnum,vy)); + this->inputs->AddInput(new TriaVertexInput(VxEnum,vx)); + this->inputs->AddInput(new TriaVertexInput(VyEnum,vy)); } @@ -497,5 +457,5 @@ /*Object functions*/ /*FUNCTION Tria::ComputeBasalStress {{{1*/ -void Tria::ComputeBasalStress(Vec eps,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::ComputeBasalStress(Vec eps,int analysis_type,int sub_analysis_type){ int i; @@ -523,11 +483,13 @@ /*}}}*/ /*FUNCTION Tria::ComputePressure {{{1*/ -void Tria::ComputePressure(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::ComputePressure(Vec pg,int analysis_type,int sub_analysis_type){ int i; - const int numgrids=3; - int doflist[numgrids]; - double pressure[numgrids]; + const int numvertices=3; + int doflist[numvertices]; + double pressure[numvertices]; + double thickness[numvertices]; double rho_ice,g; + double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}}; /*dynamic objects pointed to by hooks: */ @@ -537,8 +499,4 @@ Numpar* numpar=NULL; - - /*Get dof list on which we will plug the pressure values: */ - GetDofList1(&doflist[0]); - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -547,18 +505,25 @@ numpar=(Numpar*)hnumpar.delivers(); + /*Get dof list on which we will plug the pressure values: */ + GetDofList1(&doflist[0]); + /*pressure is lithostatic: */ rho_ice=matpar->GetRhoIce(); g=matpar->GetG(); - for(i=0;iproperties.h[i]; + + /*recover value of thickness at gauss points (0,0,1), (0,1,0),(1,0,0): */ + inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum); + + for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); /*If on water, return 0: */ - if(this->properties.onwater)return 0; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return 0; /*recover objects from hooks: */ @@ -646,5 +614,5 @@ /*First, get Misfit*/ - Jelem=Misfit(inputs,analysis_type,sub_analysis_type); + Jelem=Misfit(analysis_type,sub_analysis_type); /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ @@ -664,7 +632,6 @@ /*Add Tikhonov regularization term to misfit*/ if (strcmp(numpar->control_type,"drag")==0){ - if (!this->properties.shelf){ - - GetParameterDerivativeValue(&dk[0], &this->properties.k[0],&xyz_list[0][0], gauss_l1l2l3); + if (shelf){ + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); Jelem+=numpar->cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; @@ -672,11 +639,6 @@ } else if (strcmp(numpar->control_type,"B")==0){ - - if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ - ISSMERROR("parameter B not found in input"); - } - GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); + inputs->GetParameterDerivativeValue(&dB[0], &xyz_list[0][0], &gauss_l1l2l3[0],RheologyBEnum); Jelem+=numpar->cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; - } else{ @@ -697,10 +659,10 @@ /*FUNCTION Tria::CreateKMatrix {{{1*/ -void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==DiagnosticAnalysisEnum){ @@ -708,5 +670,5 @@ if (sub_analysis_type==HorizAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); @@ -715,25 +677,25 @@ else if (analysis_type==SlopecomputeAnalysisEnum){ - CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixSlopeCompute( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==PrognosticAnalysisEnum){ - CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixPrognostic( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==Prognostic2AnalysisEnum){ - CreateKMatrixPrognostic2(Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixPrognostic2(Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixBalancedthickness( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==Balancedthickness2AnalysisEnum){ - CreateKMatrixBalancedthickness2( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixBalancedthickness2( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixBalancedvelocities( Kgg,analysis_type,sub_analysis_type); } else{ @@ -745,5 +707,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ -void Tria::CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixBalancedthickness(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -766,4 +728,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* matrices: */ @@ -782,7 +745,6 @@ /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={0.0}; - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; + double vx_list[numgrids]; + double vy_list[numgrids]; double dvx[2]; double dvy[2]; @@ -790,4 +752,6 @@ double dvxdx,dvydy; double v_gauss[2]={0.0}; + + double K[2][2]={0.0}; double KDL[2][2]={0.0}; @@ -801,8 +765,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; /*recover objects from hooks: */ @@ -812,16 +772,11 @@ numpar=(Numpar*)hnumpar.delivers(); - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity_average in inputs!"); - - for(i=0;iGetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); //Create Artificial diffusivity once for all if requested @@ -832,6 +787,6 @@ //Build K matrix (artificial diffusivity matrix) - v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); - v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + v_gauss[0]=ONETHIRD*(vx_list[0]+vx_list[1]+vx_list[2]); + v_gauss[1]=ONETHIRD*(vy_list[0]+vy_list[1]+vy_list[2]); K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); @@ -858,9 +813,9 @@ //Get vx, vy and their derivatives at gauss point - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); - - GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); + this->GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); + this->GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + + this->GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); + this->GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); dvxdx=dvx[0]; @@ -922,5 +877,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedthickness2 {{{1*/ -void Tria::CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixBalancedthickness2(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -955,6 +910,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; double vx,vy; int dofs[1]={0}; @@ -967,8 +920,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; /*recover objects from hooks: */ @@ -977,10 +926,4 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*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!"); /* Get node coordinates and dof list: */ @@ -1008,6 +951,6 @@ //Get vx, vy and their derivatives at gauss point - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); DL_scalar=-gauss_weight*Jdettria; @@ -1039,5 +982,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedvelocities {{{1*/ -void Tria::CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixBalancedvelocities(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -1060,4 +1003,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* matrices: */ @@ -1076,6 +1020,6 @@ /*input parameters for structural analysis (diagnostic): */ double surface_normal[3]; + double surface_list[3]; double nx,ny,norm; - double vxvy_list[numgrids][2]={0.0}; double vx_list[numgrids]={0.0}; double vy_list[numgrids]={0.0}; @@ -1096,9 +1040,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -1107,12 +1046,7 @@ numpar=(Numpar*)hnumpar.delivers(); - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity_average in inputs!"); - - for(i=0;iGetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); /* Get node coordinates and dof list: */ @@ -1121,5 +1055,6 @@ /*Modify z so that it reflects the surface*/ - for(i=0;iproperties.s[i]; + inputs->GetParameterValues(&surface_list[0],&gaussgrids[0][0],3,SurfaceEnum); + for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); /*First, if we are on water, return empty matrix: */ - if(this->properties.onwater) return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater) return; /*recover objects from hooks: */ @@ -1294,8 +1229,4 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - inputs->Recover("old_velocity",&oldvxvy_list[0][0],2,dofs,numgrids,(void**)nodes); /* Get node coordinates and dof list: */ @@ -1319,9 +1250,9 @@ /*Compute thickness at gaussian point: */ - GetParameterValue(&thickness, &this->properties.h[0],gauss_l1l2l3); + inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum); /*Get strain rate from velocity: */ - GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); - GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); + inputs->GetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum); + inputs->GetStrainRate(&oldepsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxOldEnum,VyOldEnum); /*Get viscosity: */ @@ -1360,6 +1291,6 @@ /*Do not forget to include friction: */ - if(!this->properties.shelf){ - CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); + if(!shelf){ + CreateKMatrixDiagnosticHorizFriction(Kgg,analysis_type,sub_analysis_type); } @@ -1373,5 +1304,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/ -void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,int analysis_type,int sub_analysis_type){ @@ -1394,4 +1325,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* matrices: */ @@ -1411,6 +1343,10 @@ /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - int dofs[2]={0,1}; + double vx_list[numgrids]; + double vy_list[numgrids]; + double thickness_list[numgrids]; + double bed_list[numgrids]; + double dragcoefficient_list[numgrids]; + double drag_p,drag_q; /*friction: */ @@ -1420,4 +1356,8 @@ double MAXSLOPE=.06; // 6 % double MOUNTAINKEXPONENT=10; + + /*inputs: */ + bool shelf; + int drag_type; /*dynamic objects pointed to by hooks: */ @@ -1427,8 +1367,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; /*recover objects from hooks: */ @@ -1437,4 +1373,8 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); /* Get node coordinates and dof list: */ @@ -1445,13 +1385,19 @@ for(i=0;iproperties.shelf){ + if (shelf){ /*no friction, do nothing*/ return; } - if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + + /*Recover inputs: */ + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); + inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); + inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); + inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); /*Build alpha2_list used by drag stiffness matrix*/ @@ -1465,10 +1411,11 @@ friction->rho_ice=matpar->GetRhoIce(); friction->rho_water=matpar->GetRhoWater(); - friction->K=&this->properties.k[0]; - friction->bed=&this->properties.b[0]; - friction->thickness=&this->properties.h[0]; - friction->velocities=&vxvy_list[0][0]; - friction->p=this->properties.p; - friction->q=this->properties.q; + friction->K=&dragcoefficient_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->vx=&vx_list[0]; + friction->vy=&vy_list[0]; + friction->p=drag_p; + friction->q=drag_q; /*Compute alpha2_list: */ @@ -1492,5 +1439,5 @@ // If we have a slope > 6% for this element, it means we are on a mountain. In this particular case, //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */ - GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3); + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2)); @@ -1537,5 +1484,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticSurfaceVert {{{1*/ -void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,int analysis_type,int sub_analysis_type){ int i,j; @@ -1584,9 +1531,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -1675,5 +1617,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixMelting {{{1*/ -void Tria::CreateKMatrixMelting(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixMelting(Mat Kgg,int analysis_type,int sub_analysis_type){ /*indexing: */ @@ -1771,5 +1713,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/ -void Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -1792,4 +1734,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* matrices: */ @@ -1807,5 +1750,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={0.0}; double vx_list[numgrids]={0.0}; double vy_list[numgrids]={0.0}; @@ -1827,9 +1769,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -1838,19 +1775,14 @@ numpar=(Numpar*)hnumpar.delivers(); - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity_average in inputs!"); - - for(i=0;iRecover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + /*recover dt: */ + dt=numpar->dt; /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); + + /*Recover velocity: */ + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); //Create Artificial diffusivity once for all if requested @@ -1861,6 +1793,6 @@ //Build K matrix (artificial diffusivity matrix) - v_gauss[0]=ONETHIRD*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); - v_gauss[1]=ONETHIRD*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + v_gauss[0]=ONETHIRD*(vx_list[0]+vx_list[1]+vx_list[2]); + v_gauss[1]=ONETHIRD*(vy_list[0]+vy_list[1]+vy_list[2]); K[0][0]=pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); @@ -1966,5 +1898,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixPrognostic2 {{{1*/ -void Tria::CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixPrognostic2(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -2001,6 +1933,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0.0}; - double vy_list[numgrids]={0.0}; double vx,vy; double dt; @@ -2014,9 +1944,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**) hnodes.deliverp(); @@ -2025,11 +1950,6 @@ numpar=(Numpar*)hnumpar.delivers(); - /*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!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + /*recover dt: */ + dt=numpar->dt; /* Get node coordinates and dof list: */ @@ -2069,6 +1989,6 @@ //Get vx, vy and their derivatives at gauss point - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&vx,&gauss_l1l2l3[0],VxAverageEnum); + inputs->GetParameterValue(&vy,&gauss_l1l2l3[0],VyAverageEnum); DL_scalar=-dt*gauss_weight*Jdettria; @@ -2103,5 +2023,5 @@ /*FUNCTION Tria::CreateKMatrixSlopeCompute {{{1*/ -void Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixSlopeCompute(Mat Kgg,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -2195,5 +2115,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixThermal {{{1*/ -void Tria::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreateKMatrixThermal(Mat Kgg,int analysis_type,int sub_analysis_type){ int i,j; @@ -2230,5 +2150,4 @@ double tl1l2l3D[3]; double D_scalar; - ParameterInputs* inputs=NULL; /*dynamic objects pointed to by hooks: */ @@ -2238,7 +2157,4 @@ Numpar* numpar=NULL; - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -2248,6 +2164,5 @@ /*recover extra inputs from users, dt: */ - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + dt=numpar->dt; /* Get node coordinates and dof list: */ @@ -2307,10 +2222,10 @@ /*}}}*/ /*FUNCTION Tria::CreatePVector {{{1*/ -void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVector(Vec pg,int analysis_type,int sub_analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ if (analysis_type==ControlAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); } @@ -2318,5 +2233,5 @@ if (sub_analysis_type==HorizAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); @@ -2324,25 +2239,25 @@ else if (analysis_type==SlopecomputeAnalysisEnum){ - CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorSlopeCompute( pg,analysis_type,sub_analysis_type); } else if (analysis_type==PrognosticAnalysisEnum){ - CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorPrognostic( pg,analysis_type,sub_analysis_type); } else if (analysis_type==Prognostic2AnalysisEnum){ - CreatePVectorPrognostic2( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorPrognostic2( pg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedthicknessAnalysisEnum){ - CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorBalancedthickness( pg,analysis_type,sub_analysis_type); } else if (analysis_type==Balancedthickness2AnalysisEnum){ - CreatePVectorBalancedthickness2( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorBalancedthickness2( pg,analysis_type,sub_analysis_type); } else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorBalancedvelocities( pg,analysis_type,sub_analysis_type); } else{ @@ -2353,5 +2268,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedthickness {{{1*/ -void Tria::CreatePVectorBalancedthickness(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorBalancedthickness(Vec pg ,int analysis_type,int sub_analysis_type){ @@ -2382,12 +2297,6 @@ /*input parameters for structural analysis (diagnostic): */ - double accumulation_list[numgrids]={0.0}; double accumulation_g; - double melting_list[numgrids]={0.0}; double melting_g; - double thickness_list[numgrids]={0.0}; - double thickness_g; - int dofs[1]={0}; - int found=0; /*dynamic objects pointed to by hooks: */ @@ -2397,9 +2306,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -2407,10 +2311,4 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find melting in inputs!"); /* Get node coordinates and dof list: */ @@ -2435,7 +2333,7 @@ GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode); - /* Get accumulation, melting and thickness at gauss point */ - GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); - GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); + /* Get accumulation, melting at gauss point */ + inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); /* Add value into pe_g: */ @@ -2456,5 +2354,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedthickness2 {{{1*/ -void Tria::CreatePVectorBalancedthickness2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorBalancedthickness2(Vec pg ,int analysis_type,int sub_analysis_type){ @@ -2485,12 +2383,7 @@ /*input parameters for structural analysis (diagnostic): */ - double accumulation_list[numgrids]={0.0}; double accumulation_g; - double melting_list[numgrids]={0.0}; double melting_g; - double dhdt_list[numgrids]={0.0}; double dhdt_g; - int dofs[1]={0}; - int found=0; /*dynamic objects pointed to by hooks: */ @@ -2500,9 +2393,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**) hnodes.deliverp(); @@ -2510,12 +2398,4 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find melting in inputs!"); - found=inputs->Recover("dhdt",&dhdt_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find dhdt in inputs!"); /* Get node coordinates and dof list: */ @@ -2541,7 +2421,7 @@ /* Get accumulation, melting and thickness at gauss point */ - GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); - GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); - GetParameterValue(&dhdt_g, &dhdt_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&dhdt_g, &gauss_l1l2l3[0],DhDtEnum); /* Add value into pe_g: */ @@ -2562,5 +2442,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorBalancedvelocities {{{1*/ -void Tria::CreatePVectorBalancedvelocities(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorBalancedvelocities(Vec pg ,int analysis_type,int sub_analysis_type){ @@ -2591,10 +2471,6 @@ /*input parameters for structural analysis (diagnostic): */ - double accumulation_list[numgrids]={0.0}; double accumulation_g; - double melting_list[numgrids]={0.0}; double melting_g; - int dofs[1]={0}; - int found=0; /*dynamic objects pointed to by hooks: */ @@ -2604,9 +2480,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -2620,10 +2491,4 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find melting in inputs!"); /* Get node coordinates and dof list: */ @@ -2649,6 +2514,6 @@ /* Get accumulation, melting at gauss point */ - GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); - GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); /* Add value into pe_g: */ @@ -2669,5 +2534,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorDiagnosticBaseVert {{{1*/ -void Tria::CreatePVectorDiagnosticBaseVert(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorDiagnosticBaseVert(Vec pg,int analysis_type,int sub_analysis_type){ int i,j; @@ -2704,13 +2569,8 @@ /*input parameters for structural analysis (diagnostic): */ - double* velocity_param=NULL; - double vx_list[numgrids]={0,0,0}; - double vy_list[numgrids]={0,0,0}; double vx,vy; double meltingvalue; double slope[2]; double dbdx,dbdy; - int dofs1[1]={0}; - int dofs2[1]={1}; /*dynamic objects pointed to by hooks: */ @@ -2720,9 +2580,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -2731,8 +2586,4 @@ numpar=(Numpar*)hnumpar.delivers(); - /* recover input parameters: */ - if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))ISSMERROR(" cannot compute vertical velocity without horizontal velocity"); - inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); - /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); @@ -2756,12 +2607,12 @@ /*Get melting at gaussian point: */ - GetParameterValue(&meltingvalue, &this->properties.melting[0],gauss_l1l2l3); + inputs->GetParameterValue(&meltingvalue, &gauss_l1l2l3[0],MeltingRateEnum); /*Get velocity at gaussian point: */ - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); /*Get bed slope: */ - GetParameterDerivativeValue(&slope[0], &this->properties.b[0],&xyz_list[0][0], gauss_l1l2l3); + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); dbdx=slope[0]; dbdy=slope[1]; @@ -2796,5 +2647,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorDiagnosticHoriz {{{1*/ -void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorDiagnosticHoriz( Vec pg, int analysis_type,int sub_analysis_type){ int i,j; @@ -2841,11 +2692,14 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; + /*inputs: */ + bool onwater; + int drag_type; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); /*First, if we are on water, return empty vector: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return; /*recover objects from hooks: */ @@ -2875,12 +2729,11 @@ /*Compute thickness at gaussian point: */ - GetParameterValue(&thickness, &this->properties.h[0],gauss_l1l2l3); - - GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_l1l2l3); + inputs->GetParameterValue(&thickness, &gauss_l1l2l3[0],ThicknessEnum); + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the * element itself: */ - if(this->properties.friction_type==1){ - GetParameterValue(&plastic_stress, &this->properties.k[0],gauss_l1l2l3); + if(drag_type==1){ + inputs->GetParameterValue(&plastic_stress, &gauss_l1l2l3[0],DragCoefficientEnum); } @@ -2895,5 +2748,5 @@ /*Build pe_g_gaussian vector: */ - if(this->properties.friction_type==1){ + if(drag_type==1){ for (i=0;iRecover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find melting in inputs!"); - found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find thickness in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + /*recover dt: */ + dt=numpar->dt; /* Get node coordinates and dof list: */ @@ -3015,7 +2851,7 @@ /* Get accumulation, melting and thickness at gauss point */ - GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); - GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); - GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); /* Add value into pe_g: */ @@ -3036,5 +2872,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorPrognostic2 {{{1*/ -void Tria::CreatePVectorPrognostic2(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorPrognostic2(Vec pg ,int analysis_type,int sub_analysis_type){ /* local declarations */ @@ -3064,13 +2900,8 @@ /*input parameters for structural analysis (diagnostic): */ - double accumulation_list[numgrids]={0.0}; double accumulation_g; - double melting_list[numgrids]={0.0}; double melting_g; - double thickness_list[numgrids]={0.0}; double thickness_g; double dt; - int dofs[1]={0}; - int found=0; /*dynamic objects pointed to by hooks: */ @@ -3080,9 +2911,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -3091,13 +2917,7 @@ numpar=(Numpar*)hnumpar.delivers(); - /*recover extra inputs from users, at current convergence iteration: */ - found=inputs->Recover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find accumulation in inputs!"); - found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find melting in inputs!"); - found=inputs->Recover("thickness",&thickness_list[0],1,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find thickness in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); + /*recover dt: */ + dt=numpar->dt; + /* Get node coordinates and dof list: */ @@ -3123,7 +2943,7 @@ /* Get accumulation, melting and thickness at gauss point */ - GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); - GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); - GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&accumulation_g, &gauss_l1l2l3[0],AccumulationRateEnum); + inputs->GetParameterValue(&melting_g, &gauss_l1l2l3[0],MeltingRateEnum); + inputs->GetParameterValue(&thickness_g, &gauss_l1l2l3[0],ThicknessEnum); /* Add value into pe_g: */ @@ -3145,5 +2965,5 @@ /*FUNCTION Tria::CreatePVectorSlopeCompute {{{1*/ -void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorSlopeCompute( Vec pg, int analysis_type,int sub_analysis_type){ int i,j; @@ -3175,5 +2995,4 @@ double pe_g[numdof]; double pe_g_gaussian[numdof]; - double param[3]; double slope[2]; @@ -3197,10 +3016,4 @@ for(i=0;iproperties.s[i]; - } - if ( (sub_analysis_type==BedXAnalysisEnum) || (sub_analysis_type==BedYAnalysisEnum)){ - for(i=0;iproperties.b[i]; - } /* Get gaussian points and weights: */ @@ -3216,6 +3029,11 @@ gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); - GetParameterDerivativeValue(&slope[0], ¶m[0],&xyz_list[0][0], gauss_l1l2l3); - + if ( (sub_analysis_type==SurfaceXAnalysisEnum) || (sub_analysis_type==SurfaceYAnalysisEnum)){ + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],SurfaceEnum); + } + if ( (sub_analysis_type==BedXAnalysisEnum) || (sub_analysis_type==BedYAnalysisEnum)){ + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],&gauss_l1l2l3[0],BedEnum); + } + /* Get Jacobian determinant: */ GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3); @@ -3249,5 +3067,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorThermalShelf {{{1*/ -void Tria::CreatePVectorThermalShelf( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorThermalShelf( Vec pg, int analysis_type,int sub_analysis_type){ int i,found; @@ -3270,5 +3088,4 @@ /*inputs: */ double dt; - double pressure_list[3]; double pressure; @@ -3297,9 +3114,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -3320,11 +3132,6 @@ beta=matpar->GetBeta(); meltingpoint=matpar->GetMeltingPoint(); - - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); - found=inputs->Recover("pressure",&pressure_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find pressure in inputs!"); + dt=numpar->dt; + /* Ice/ocean heat exchange flux on ice shelf base */ @@ -3346,5 +3153,5 @@ /*Get geothermal flux and basal friction */ - GetParameterValue(&pressure,&pressure_list[0],gauss_coord); + inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); t_pmp=meltingpoint-beta*pressure; @@ -3372,5 +3179,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVectorThermalSheet {{{1*/ -void Tria::CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){ +void Tria::CreatePVectorThermalSheet( Vec pg, int analysis_type,int sub_analysis_type){ int i,found; @@ -3382,7 +3189,4 @@ int numberofdofspernode; double xyz_list[numgrids][3]; - double vxvyvz_list[numgrids][3]; - double vx_list[numgrids]; - double vy_list[numgrids]; double rho_ice; @@ -3398,4 +3202,12 @@ double geothermalflux_value; + double vx_list[numgrids]; + double vy_list[numgrids]; + double thickness_list[numgrids]; + double bed_list[numgrids]; + double dragcoefficient_list[numgrids]; + double drag_p,drag_q; + int drag_type; + /* gaussian points: */ int num_area_gauss,ig; @@ -3406,5 +3218,5 @@ double gauss_weight; double gauss_coord[3]; - int dofs1[1]={0}; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /*matrices: */ @@ -3414,7 +3226,4 @@ double scalar; - int dofs[3]={0,1,2}; - - ParameterInputs* inputs=NULL; /*dynamic objects pointed to by hooks: */ @@ -3424,7 +3233,4 @@ Numpar* numpar=NULL; - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -3432,4 +3238,7 @@ matice=(Matice*)hmatice.delivers(); numpar=(Numpar*)hnumpar.delivers(); + + /*retrieve inputs :*/ + inputs->GetParameterValue(&drag_type,DragTypeEnum); /* Get node coordinates and dof list: */ @@ -3440,17 +3249,15 @@ rho_ice=matpar->GetRhoIce(); heatcapacity=matpar->GetHeatCapacity(); - - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); - - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity in inputs!"); - - for(i=0;idt; + + + /*Recover inputs: */ + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); + inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); + inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); + inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); /*Build alpha2_list used by drag stiffness matrix*/ @@ -3458,18 +3265,19 @@ /*Initialize all fields: */ - if (this->properties.friction_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); + if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!"); friction->element_type=(char*)xmalloc((strlen("3d")+1)*sizeof(char)); strcpy(friction->element_type,"3d"); - + friction->gravity=matpar->GetG(); friction->rho_ice=matpar->GetRhoIce(); friction->rho_water=matpar->GetRhoWater(); - friction->K=&this->properties.k[0]; - friction->bed=&this->properties.b[0]; - friction->thickness=&this->properties.h[0]; - friction->velocities=&vxvyvz_list[0][0]; - friction->p=this->properties.p; - friction->q=this->properties.q; + friction->K=&dragcoefficient_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->vx=&vx_list[0]; + friction->vy=&vy_list[0]; + friction->p=drag_p; + friction->q=drag_q; /*Compute alpha2_list: */ @@ -3501,5 +3309,5 @@ /*Get geothermal flux and basal friction */ - GetParameterValue(&geothermalflux_value,&this->properties.geothermalflux[0],gauss_coord); + inputs->GetParameterValue(&geothermalflux_value, &gauss_coord[0],GeothermalFluxEnum); GetParameterValue(&basalfriction,&basalfriction_list[0],gauss_coord); @@ -3527,5 +3335,5 @@ /*}}}*/ /*FUNCTION Tria::Du {{{1*/ -void Tria::Du(Vec du_g,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::Du(Vec du_g,int analysis_type,int sub_analysis_type){ int i; @@ -3538,12 +3346,9 @@ int doflist[numdof]; int numberofdofspernode; - int dofs2[2]={0,1}; - int dofs1[1]={0}; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* grid data: */ - double vxvy_list[numgrids][2]; double vx_list[numgrids]; double vy_list[numgrids]; - double obs_vxvy_list[numgrids][2]; double obs_vx_list[numgrids]; double obs_vy_list[numgrids]; @@ -3566,5 +3371,5 @@ /*element vector : */ - double due_g[numdof]; + double due_g[numdof]={0,0,0,0,0,0}; double due_g_gaussian[numdof]; @@ -3588,9 +3393,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -3603,27 +3403,16 @@ GetDofList(&doflist[0],&numberofdofspernode); - /* Set due_g to 0: */ - for(i=0;iRecover("fit",&fit)) ISSMERROR(" missing fit input parameter"); - if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity_obs input parameter"); - } - if(fit==3 && !inputs->Recover("surfacearea",&S)){ - ISSMERROR("missing surface area input parameter"); - } - if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity input parameter"); - } - if(!inputs->Recover("weights",&weights_list[0],1,dofs1,numgrids,(void**)nodes)){ - ISSMERROR("missing weights input parameter"); - } - - for(i=0;iGetParameterValues(&obs_vx_list[0],&gaussgrids[0][0],3,VxObsEnum); + inputs->GetParameterValues(&obs_vy_list[0],&gaussgrids[0][0],3,VyObsEnum); + + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyEnum); + + inputs->GetParameterValues(&weights_list[0],&gaussgrids[0][0],3,WeightsEnum); + + inputs->GetParameterValue(&fit,FitEnum); + if(fit==3){ + inputs->GetParameterValue(&S,SurfaceAreaEnum); } @@ -3927,8 +3716,10 @@ /*}}}*/ /*FUNCTION Tria::GetBedList {{{1*/ -void Tria::GetBedList(double* bed_list){ - - int i; - for(i=0;i<3;i++)bed_list[i]=this->properties.b[i]; +void Tria::GetBedList(double* bedlist){ + + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + + inputs->GetParameterValues(bedlist,&gaussgrids[0][0],3,BedEnum); } @@ -4296,6 +4087,10 @@ /*}}}*/ /*FUNCTION Tria::GetOnBed {{{1*/ -int Tria::GetOnBed(){ - return this->properties.onbed; +bool Tria::GetOnBed(){ + + bool onbed; + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + return onbed; } /*}}}*/ @@ -4344,6 +4139,12 @@ /*}}}*/ /*FUNCTION Tria::GetShelf {{{1*/ -int Tria::GetShelf(){ - return this->properties.shelf; +bool Tria::GetShelf(){ + /*inputs: */ + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + return shelf; } /*}}}*/ @@ -4370,19 +4171,27 @@ void Tria::GetThicknessList(double* thickness_list){ - int i; - for(i=0;i<3;i++)thickness_list[i]=this->properties.h[i]; + const int numgrids=3; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; + inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],3,ThicknessEnum); + } /*}}}*/ /*FUNCTION Tria::Gradj {{{1*/ -void Tria::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){ +void Tria::Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type){ + + /*inputs: */ + bool onwater; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); /*If on water, grad = 0: */ - if(this->properties.onwater)return; + if(onwater)return; if (strcmp(control_type,"drag")==0){ - GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); + GradjDrag( grad_g,analysis_type,sub_analysis_type); } else if (strcmp(control_type,"B")==0){ - GradjB( grad_g,inputs,analysis_type,sub_analysis_type); + GradjB( grad_g,analysis_type,sub_analysis_type); } else ISSMERROR("%s%s","control type not supported yet: ",control_type); @@ -4390,5 +4199,5 @@ /*}}}*/ /*FUNCTION Tria::GradjB{{{1*/ -void Tria::GradjB(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::GradjB(Vec grad_g,int analysis_type,int sub_analysis_type){ int i; @@ -4404,10 +4213,4 @@ /* grid data: */ - double vx_list[numgrids]; - double vy_list[numgrids]; - double vxvy_list[numgrids][2]; - double adjx_list[numgrids]; - double adjy_list[numgrids]; - double adjxadjy_list[numgrids][2]; double B[numgrids]; @@ -4456,9 +4259,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -4473,24 +4271,4 @@ /* Set grade_g to 0: */ for(i=0;iRecover("thickness",&this->properties.h[0],1,dofs,numgrids,(void**)nodes); - if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity input parameter"); - } - if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing adjoint input parameter"); - } - if(!inputs->Recover("B",&B[0],1,dofs1,numgrids,(void**)nodes)){ - ISSMERROR("parameter B not found in input"); - } - - /*Initialize parameter lists: */ - for(i=0;iproperties.h[0],gauss_l1l2l3); + inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum); /*Get strain rate, if velocity has been supplied: */ - GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3); + inputs->GetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_l1l2l3,VxEnum,VyEnum); /*Get viscosity complement: */ @@ -4515,8 +4293,8 @@ /*Get dvx, dvy, dadjx and dadjx: */ - GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dadjx[0], &adjx_list[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterDerivativeValue(&dadjy[0], &adjy_list[0],&xyz_list[0][0], gauss_l1l2l3); + 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); /* Get Jacobian determinant: */ @@ -4530,6 +4308,6 @@ /*Get B derivative: dB/dx */ - GetParameterDerivativeValue(&dB[0], &B[0],&xyz_list[0][0], gauss_l1l2l3); - GetParameterValue(&B_gauss, &B[0],gauss_l1l2l3); + inputs->GetParameterDerivativeValue(&dB[0],&xyz_list[0][0],&gauss_l1l2l3[0],RheologyBEnum); + inputs->GetParameterValue(&B_gauss, gauss_l1l2l3,RheologyBEnum); /*Build gradje_g_gaussian vector (actually -dJ/dB): */ @@ -4568,5 +4346,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDrag {{{1*/ -void Tria::GradjDrag(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::GradjDrag(Vec grad_g,int analysis_type,int sub_analysis_type){ @@ -4584,12 +4362,14 @@ double vx_list[numgrids]; double vy_list[numgrids]; - double vxvy_list[numgrids][2]; double adjx_list[numgrids]; double adjy_list[numgrids]; - double adjxadjy_list[numgrids][2]; + double thickness_list[numgrids]; + double bed_list[numgrids]; + double dragcoefficient_list[numgrids]; + double drag_p; + double drag_q; + int drag_type; double drag; - int dofs1[1]={0}; - int dofs2[2]={0,1}; /* gaussian points: */ @@ -4601,4 +4381,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* parameters: */ @@ -4609,10 +4390,9 @@ /*drag: */ - double pcoeff,qcoeff; double alpha_complement_list[numgrids]; double alpha_complement; /*element vector at the gaussian points: */ - double grade_g[numgrids]; + double grade_g[numgrids]={0,0,0}; double grade_g_gaussian[numgrids]; @@ -4625,4 +4405,7 @@ /* strain rate: */ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + + /*inputs: */ + bool shelf; /*dynamic objects pointed to by hooks: */ @@ -4632,9 +4415,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -4643,6 +4421,9 @@ numpar=(Numpar*)hnumpar.delivers(); + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*Get out if shelf*/ - if(this->properties.shelf) return; + if(shelf)return; /* Get node coordinates and dof list: */ @@ -4650,25 +4431,13 @@ GetDofList1(&doflist1[0]); - /* Set grade_g to 0: */ - for(i=0;iRecover("drag",&this->properties.k[0],1,dofs1,numgrids,(void**)nodes); - inputs->Recover("bed",&this->properties.b[0],1,dofs1,numgrids,(void**)nodes); - inputs->Recover("thickness",&this->properties.h[0],1,dofs1,numgrids,(void**)nodes); - if(!inputs->Recover("velocity",&vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity input parameter"); - } - if(!inputs->Recover("adjoint",&adjxadjy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing adjoint input parameter"); - } - - /*Initialize parameter lists: */ - for(i=0;iGetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); + inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],3,ThicknessEnum); + inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],3,BedEnum); + inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],3,DragCoefficientEnum); + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ @@ -4684,5 +4453,5 @@ /*Build alpha_complement_list: */ - if (!this->properties.shelf && (this->properties.friction_type==2)){ + if (drag_type==2){ /*Allocate friction object: */ @@ -4693,14 +4462,17 @@ strcpy(friction->element_type,"2d"); + friction->gravity=matpar->GetG(); friction->rho_ice=matpar->GetRhoIce(); friction->rho_water=matpar->GetRhoWater(); - friction->K=&this->properties.k[0]; - friction->bed=&this->properties.b[0]; - friction->thickness=&this->properties.h[0]; - friction->velocities=&vxvy_list[0][0]; - friction->p=this->properties.p; - friction->q=this->properties.q; - + friction->K=&dragcoefficient_list[0]; + friction->bed=&bed_list[0]; + friction->thickness=&thickness_list[0]; + friction->vx=&vx_list[0]; + friction->vy=&vy_list[0]; + friction->p=drag_p; + friction->q=drag_q; + + if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!"); @@ -4719,13 +4491,13 @@ /*Recover alpha_complement and k: */ GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3); - GetParameterValue(&drag, &this->properties.k[0],gauss_l1l2l3); + inputs->GetParameterValue(&drag, gauss_l1l2l3,DragCoefficientEnum); /*recover lambda and mu: */ - GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3); - GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&lambda, gauss_l1l2l3,AdjointxEnum); + inputs->GetParameterValue(&mu, gauss_l1l2l3,AdjointyEnum); /*recover vx and vy: */ - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&vx, gauss_l1l2l3,VxEnum); + inputs->GetParameterValue(&vy, gauss_l1l2l3,VyEnum); /* Get Jacobian determinant: */ @@ -4739,5 +4511,5 @@ /*Get k derivative: dk/dx */ - GetParameterDerivativeValue(&dk[0], &this->properties.k[0],&xyz_list[0][0], gauss_l1l2l3); + inputs->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],&gauss_l1l2l3[0],DragCoefficientEnum); /*Build gradje_g_gaussian vector (actually -dJ/ddrag): */ @@ -4777,5 +4549,5 @@ /*}}}*/ /*FUNCTION Tria::GradjDragStokes {{{1*/ -void Tria::GradjDragStokes(Vec grad_g,void* vinputs,int analysis_type,int sub_analysis_type){ +void Tria::GradjDragStokes(Vec grad_g,int analysis_type,int sub_analysis_type){ int i; @@ -4837,5 +4609,7 @@ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ - ParameterInputs* inputs=NULL; + /*inputs: */ + bool shelf; + int drag_type; /*dynamic objects pointed to by hooks: */ @@ -4845,9 +4619,10 @@ Numpar* numpar=NULL; + /*retrieve inputs :*/ + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + inputs->GetParameterValue(&drag_type,DragTypeEnum); + /*Get out if shelf*/ - if(this->properties.shelf) return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(shelf)return; /*recover objects from hooks: */ @@ -4897,5 +4672,5 @@ /*Build alpha_complement_list: */ - if (!this->properties.shelf && (this->properties.friction_type==2)){ + if (drag_type==2){ /*Allocate friction object: */ @@ -5080,5 +4855,5 @@ /*}}}*/ /*FUNCTION Tria::Misfit {{{1*/ -double Tria::Misfit(void* vinputs,int analysis_type,int sub_analysis_type){ +double Tria::Misfit(int analysis_type,int sub_analysis_type){ int i; @@ -5127,6 +4902,4 @@ double fit=-1; - ParameterInputs* inputs=NULL; - /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; @@ -5135,9 +4908,12 @@ Numpar* numpar=NULL; + /*inputs: */ + bool onwater; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*If on water, return 0: */ - if(this->properties.onwater)return 0; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return 0; /*recover objects from hooks: */ @@ -5329,5 +5105,5 @@ /*}}}*/ /*FUNCTION Tria::SurfaceArea {{{1*/ -double Tria::SurfaceArea(void* vinputs,int analysis_type,int sub_analysis_type){ +double Tria::SurfaceArea(int analysis_type,int sub_analysis_type){ int i; @@ -5346,9 +5122,15 @@ Node** nodes=NULL; + /*inputs: */ + bool onwater; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); /*If on water, return 0: */ - if(this->properties.onwater)return 0; + if(onwater)return 0; /* Get node coordinates and dof list: */ Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 3611) +++ /issm/trunk/src/c/objects/Tria.h (revision 3612) @@ -9,6 +9,6 @@ class Element; class Hook; -class ElementProperties; class DataSet; +class Inputs; class Node; struct IoModel; @@ -18,7 +18,7 @@ #include "./Hook.h" #include "./Node.h" -#include "./ElementProperties.h" #include "../ModelProcessorx/IoModel.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Inputs.h" class Tria: public Element{ @@ -33,6 +33,5 @@ Hook hnumpar; //hook to 1 numpar - ElementProperties properties; - DataSet* inputs; + Inputs* inputs; public: @@ -40,6 +39,6 @@ /*FUNCTION constructors, destructors {{{1*/ Tria(); - Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id, ElementProperties* tria_properties); - Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, ElementProperties* tria_properties,DataSet* tria_inputs); + Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id, int tria_numpar_id); + Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs); Tria(int i, IoModel* iomodel); ~Tria(); @@ -60,12 +59,12 @@ /*}}}*/ /*FUNCTION element numerical routines {{{1*/ - void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVector(Vec pg, void* inputs, int analysis_type,int sub_analysis_type); + void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVector(Vec pg, int analysis_type,int sub_analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void GetDofList1(int* doflist); - void CreateKMatrixDiagnosticHoriz(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHoriz(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticHorizFriction(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixDiagnosticSurfaceVert(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixSlopeCompute(Mat Kgg,int analysis_type,int sub_analysis_type); void GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3); void GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3); @@ -83,41 +82,40 @@ void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3); void GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3); - void Du(Vec du_g,void* inputs,int analysis_type,int sub_analysis_type); - void Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type); - void GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - void GradjDragStokes(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); + void Du(Vec du_g,int analysis_type,int sub_analysis_type); + void Gradj(Vec grad_g,int analysis_type,int sub_analysis_type,char* control_type); + void GradjDrag(Vec grad_g,int analysis_type,int sub_analysis_type); + void GradjDragStokes(Vec grad_g,int analysis_type,int sub_analysis_type); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); - void GradjB(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type); - double Misfit(void* inputs,int analysis_type,int sub_analysis_type); - double SurfaceArea(void* inputs,int analysis_type,int sub_analysis_type); - double CostFunction(void* inputs,int analysis_type,int sub_analysis_type); + void GradjB(Vec grad_g,int analysis_type,int sub_analysis_type); + double Misfit(int analysis_type,int sub_analysis_type); + double SurfaceArea(int analysis_type,int sub_analysis_type); + double CostFunction(int analysis_type,int sub_analysis_type); - void CreatePVectorDiagnosticHoriz(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVectorDiagnosticBaseVert(Vec pg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticHoriz(Vec pg,int analysis_type,int sub_analysis_type); + void CreatePVectorDiagnosticBaseVert(Vec pg,int analysis_type,int sub_analysis_type); + void CreatePVectorSlopeCompute( Vec pg, int analysis_type,int sub_analysis_type); void* GetMatPar(); - int GetShelf(); + bool GetShelf(); void GetNodes(void** nodes); - int GetOnBed(); - + bool GetOnBed(); void GetThicknessList(double* thickness_list); void GetBedList(double* bed_list); - void ComputeBasalStress(Vec sigma_b,void* inputs,int analysis_type,int sub_analysis_type); - void ComputePressure(Vec p_g,void* inputs,int analysis_type,int sub_analysis_type); - void ComputeStrainRate(Vec eps,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixThermal(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type); - void CreatePVectorThermalShelf( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); - void CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorPrognostic(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixPrognostic2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorPrognostic2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorBalancedthickness(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixBalancedthickness2(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorBalancedthickness2(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); - void CreatePVectorBalancedvelocities(Vec pg,void* vinputs,int analysis_type,int sub_analysis_type); + void ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type); + void ComputePressure(Vec p_g,int analysis_type,int sub_analysis_type); + void ComputeStrainRate(Vec eps,int analysis_type,int sub_analysis_type); + void CreateKMatrixThermal(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreateKMatrixMelting(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorThermalShelf( Vec pg, int analysis_type,int sub_analysis_type); + void CreatePVectorThermalSheet( Vec pg, int analysis_type,int sub_analysis_type); + void CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorPrognostic(Vec pg,int analysis_type,int sub_analysis_type); + void CreateKMatrixPrognostic2(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorPrognostic2(Vec pg,int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedthickness(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedthickness(Vec pg,int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedthickness2(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedthickness2(Vec pg,int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedvelocities(Mat Kgg,int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedvelocities(Vec pg,int analysis_type,int sub_analysis_type); double MassFlux(double* segment,double* ug); double GetArea(void); @@ -126,5 +124,4 @@ /*updates:*/ void UpdateFromDakota(void* inputs); - void UpdateFromInputs(void* inputs); void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); void UpdateInputsDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/TriaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/TriaVertexInput.h (revision 3611) +++ /issm/trunk/src/c/objects/TriaVertexInput.h (revision 3612) @@ -3,10 +3,10 @@ */ -#include "./Einput.h" +#include "./Input.h" #ifndef _TRIAVERTEXINPUT_H_ #define _TRIAVERTEXINPUT_H_ -class TriaVertexInput: public Einput{ +class TriaVertexInput: public Input{ private: Index: /issm/trunk/src/c/objects/Vertex.cpp =================================================================== --- /issm/trunk/src/c/objects/Vertex.cpp (revision 3611) +++ /issm/trunk/src/c/objects/Vertex.cpp (revision 3612) @@ -182,23 +182,4 @@ } /*}}}*/ -/*FUNCTION UpdateFromInputs {{{2*/ -void Vertex::UpdateFromInputs(void* vinputs){ - - ParameterInputs* inputs=NULL; - Vertex* vertex=NULL; - int dof[1]={0}; - - vertex=this; - - /*Recover parameter inputs: */ - inputs=(ParameterInputs*)vinputs; - - /*Update internal data if inputs holds new values: */ - inputs->Recover("x",&x,1,dof,1,(void**)&vertex); - inputs->Recover("y",&y,1,dof,1,(void**)&vertex); - inputs->Recover("z",&z,1,dof,1,(void**)&vertex); - -} -/*}}}*/ /*FUNCTION UpdateVertexPosition {{{2*/ void Vertex::UpdatePosition(double* thickness,double* bed){ Index: /issm/trunk/src/c/objects/Vertex.h =================================================================== --- /issm/trunk/src/c/objects/Vertex.h (revision 3611) +++ /issm/trunk/src/c/objects/Vertex.h (revision 3612) @@ -46,5 +46,4 @@ int MyRank(); void UpdateFromDakota(void* vinputs); - void UpdateFromInputs(void* vinputs); void UpdatePosition(double* thickness,double* bed); Index: /issm/trunk/src/c/objects/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 3611) +++ /issm/trunk/src/c/objects/objects.h (revision 3612) @@ -25,4 +25,9 @@ class NodeSets; class Model; +class TriaVertexInput; +class DoubleInput; +class IntInput; +class BoolInput; +class Input; /*Abstract class: */ @@ -50,9 +55,12 @@ #include "./NodeSets.h" #include "./Model.h" +#include "./Input.h" +#include "./TriaVertexInput.h" +#include "./BoolInput.h" +#include "./IntInput.h" +#include "./DoubleInput.h" /*C objects: */ #include "./Contour.h" -#include "./ParameterInputs.h" -#include "./Input.h" #include "./Friction.h" #include "./SolverEnum.h" Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 3611) +++ /issm/trunk/src/c/parallel/parallel.h (revision 3612) @@ -10,36 +10,35 @@ struct OptArgs; -class ParameterInputs; class FemModel; -void gradjcompute_core(DataSet* results,Model* model, ParameterInputs* inputs); +void gradjcompute_core(DataSet* results,Model* model); -void diagnostic_core(DataSet* results,Model* model, ParameterInputs* inputs); -void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs); -void prognostic2_core(DataSet* results,Model* model, ParameterInputs* inputs); -void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs); -void balancedthickness2_core(DataSet* results,Model* model, ParameterInputs* inputs); -void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs); -void slopecompute_core(DataSet* results,Model* model, ParameterInputs* inputs); -void control_core(DataSet* results,Model* model, ParameterInputs* inputs); +void diagnostic_core(DataSet* results,Model* model); +void prognostic_core(DataSet* results,Model* model); +void prognostic2_core(DataSet* results,Model* model); +void balancedthickness_core(DataSet* results,Model* model); +void balancedthickness2_core(DataSet* results,Model* model); +void balancedvelocities_core(DataSet* results,Model* model); +void slopecompute_core(DataSet* results,Model* model); +void control_core(DataSet* results,Model* model); -void thermal_core(DataSet* results,Model* model, ParameterInputs* inputs); -void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void thermal_core(DataSet* results,Model* model); +void thermal_core_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type); -void steadystate_core(DataSet* results,Model* model, ParameterInputs* inputs); +void steadystate_core(DataSet* results,Model* model); -void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); -void diagnostic_core_linear(Vec* ppg,FemModel* fem,ParameterInputs* inputs,int analysis_type,int sub_analysis_type); +void diagnostic_core_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, DataSet* loads, FemModel* fem,int analysis_type,int sub_analysis_type); +void diagnostic_core_linear(Vec* ppg,FemModel* fem,int analysis_type,int sub_analysis_type); void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,DataSet* parameters); -void transient_core(DataSet* results,Model* model, ParameterInputs* inputs); -void transient_core_2d(DataSet* results,Model* model, ParameterInputs* inputs); -void transient_core_3d(DataSet* results,Model* model, ParameterInputs* inputs); +void transient_core(DataSet* results,Model* model); +void transient_core_2d(DataSet* results,Model* model); +void transient_core_3d(DataSet* results,Model* model); //int GradJOrth(WorkspaceParams* workspaceparams); -int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel); +int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); -int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*,ParameterInputs*),FemModel* femmodel); +int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel); double objectivefunctionC(double search_scalar,OptArgs* optargs); @@ -52,5 +51,5 @@ void WriteLockFile(char* filename); -void ControlInitialization(Model* model, ParameterInputs* inputs); +void ControlInitialization(Model* model); void ControlRestart(Model* model,double* param_g); Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 3611) +++ /issm/trunk/src/m/classes/@model/model.m (revision 3612) @@ -124,10 +124,10 @@ md.g=0; md.yts=0; - md.drag=NaN; md.drag_type=0; - md.p=NaN; - md.q=NaN; - md.B=NaN; - md.n=NaN; + md.drag_coefficient=NaN; + md.drag_p=NaN; + md.drag_q=NaN; + md.rheology_B=NaN; + md.rheology_n=NaN; %Geometrical parameters @@ -160,5 +160,5 @@ md.vel_bal=NaN; md.vel_obs_raw=NaN; - md.accumulation=NaN; + md.accumulation_rate=NaN; md.dhdt=NaN; md.geothermalflux=NaN; @@ -237,5 +237,5 @@ md.vel=NaN; md.temperature=NaN; %temperature solution vector - md.melting=NaN; + md.melting_rate=NaN; md.pressure=NaN; Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 3611) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 3612) @@ -64,11 +64,11 @@ WriteData(fid,md.pressure,'Mat','pressure'); WriteData(fid,md.temperature,'Mat','temperature'); -WriteData(fid,md.melting,'Mat','melting'); +WriteData(fid,md.melting_rate,'Mat','melting_rate'); WriteData(fid,md.drag_type,'Integer','drag_type'); -WriteData(fid,md.drag,'Mat','drag'); -WriteData(fid,md.p,'Mat','p'); -WriteData(fid,md.q,'Mat','q'); +WriteData(fid,md.drag_coefficient,'Mat','drag_coefficient'); +WriteData(fid,md.drag_p,'Mat','drag_p'); +WriteData(fid,md.drag_q,'Mat','drag_q'); WriteData(fid,md.elementoniceshelf,'Mat','elementoniceshelf'); @@ -86,5 +86,5 @@ WriteData(fid,md.geothermalflux,'Mat','geothermalflux'); WriteData(fid,md.melting,'Mat','melting'); -WriteData(fid,md.accumulation,'Mat','accumulation'); +WriteData(fid,md.accumulation_rate,'Mat','accumulation_rate'); WriteData(fid,md.dhdt,'Mat','dhdt'); @@ -93,6 +93,6 @@ WriteData(fid,md.rho_ice,'Scalar','rho_ice'); WriteData(fid,md.g,'Scalar','g'); -WriteData(fid,md.B,'Mat','B'); -WriteData(fid,md.n,'Mat','n'); +WriteData(fid,md.rheology_B,'Mat','rheology_B'); +WriteData(fid,md.rheology_n,'Mat','rheology_n'); %Control methods Index: /issm/trunk/src/m/enum/AccumulationRateEnum.m =================================================================== --- /issm/trunk/src/m/enum/AccumulationRateEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/AccumulationRateEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=AccumulationRateEnum() +%ACCUMULATIONRATEENUM - Enum of AccumulationRate +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=AccumulationRateEnum() + +macro=105; Index: /issm/trunk/src/m/enum/AdjointxEnum.m =================================================================== --- /issm/trunk/src/m/enum/AdjointxEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/AdjointxEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=AdjointxEnum() +%ADJOINTXENUM - Enum of Adjointx +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=AdjointxEnum() + +macro=114; Index: /issm/trunk/src/m/enum/AdjointyEnum.m =================================================================== --- /issm/trunk/src/m/enum/AdjointyEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/AdjointyEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=AdjointyEnum() +%ADJOINTYENUM - Enum of Adjointy +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=AdjointyEnum() + +macro=115; Index: /issm/trunk/src/m/enum/AirEnum.m =================================================================== --- /issm/trunk/src/m/enum/AirEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/AirEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=AirEnum() -macro=77; +macro=80; Index: /issm/trunk/src/m/enum/BedEnum.m =================================================================== --- /issm/trunk/src/m/enum/BedEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/BedEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=BedEnum() +%BEDENUM - Enum of Bed +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=BedEnum() + +macro=97; Index: /issm/trunk/src/m/enum/BoolInputEnum.m =================================================================== --- /issm/trunk/src/m/enum/BoolInputEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/BoolInputEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=BoolInputEnum() +%BOOLINPUTENUM - Enum of BoolInput +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=BoolInputEnum() + +macro=67; Index: /issm/trunk/src/m/enum/DhDtEnum.m =================================================================== --- /issm/trunk/src/m/enum/DhDtEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DhDtEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DhDtEnum() +%DHDTENUM - Enum of DhDt +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DhDtEnum() + +macro=94; Index: /issm/trunk/src/m/enum/DofVecEnum.m =================================================================== --- /issm/trunk/src/m/enum/DofVecEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/DofVecEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=DofVecEnum() -macro=71; +macro=74; Index: /issm/trunk/src/m/enum/DoubleInputEnum.m =================================================================== --- /issm/trunk/src/m/enum/DoubleInputEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DoubleInputEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DoubleInputEnum() +%DOUBLEINPUTENUM - Enum of DoubleInput +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DoubleInputEnum() + +macro=69; Index: /issm/trunk/src/m/enum/DragCoefficientEnum.m =================================================================== --- /issm/trunk/src/m/enum/DragCoefficientEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DragCoefficientEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DragCoefficientEnum() +%DRAGCOEFFICIENTENUM - Enum of DragCoefficient +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DragCoefficientEnum() + +macro=98; Index: /issm/trunk/src/m/enum/DragPEnum.m =================================================================== --- /issm/trunk/src/m/enum/DragPEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DragPEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DragPEnum() +%DRAGPENUM - Enum of DragP +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DragPEnum() + +macro=99; Index: /issm/trunk/src/m/enum/DragQEnum.m =================================================================== --- /issm/trunk/src/m/enum/DragQEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DragQEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DragQEnum() +%DRAGQENUM - Enum of DragQ +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DragQEnum() + +macro=100; Index: /issm/trunk/src/m/enum/DragTypeEnum.m =================================================================== --- /issm/trunk/src/m/enum/DragTypeEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/DragTypeEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=DragTypeEnum() +%DRAGTYPEENUM - Enum of DragType +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=DragTypeEnum() + +macro=101; Index: /issm/trunk/src/m/enum/ElementOnBedEnum.m =================================================================== --- /issm/trunk/src/m/enum/ElementOnBedEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/ElementOnBedEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=ElementOnBedEnum() +%ELEMENTONBEDENUM - Enum of ElementOnBed +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=ElementOnBedEnum() + +macro=108; Index: /issm/trunk/src/m/enum/ElementOnIceShelfEnum.m =================================================================== --- /issm/trunk/src/m/enum/ElementOnIceShelfEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/ElementOnIceShelfEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=ElementOnIceShelfEnum() +%ELEMENTONICESHELFENUM - Enum of ElementOnIceShelf +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=ElementOnIceShelfEnum() + +macro=107; Index: /issm/trunk/src/m/enum/ElementOnSurfaceEnum.m =================================================================== --- /issm/trunk/src/m/enum/ElementOnSurfaceEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/ElementOnSurfaceEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=ElementOnSurfaceEnum() +%ELEMENTONSURFACEENUM - Enum of ElementOnSurface +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=ElementOnSurfaceEnum() + +macro=110; Index: /issm/trunk/src/m/enum/ElementOnWaterEnum.m =================================================================== --- /issm/trunk/src/m/enum/ElementOnWaterEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/ElementOnWaterEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=ElementOnWaterEnum() +%ELEMENTONWATERENUM - Enum of ElementOnWater +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=ElementOnWaterEnum() + +macro=109; Index: /issm/trunk/src/m/enum/FitEnum.m =================================================================== --- /issm/trunk/src/m/enum/FitEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/FitEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=FitEnum() +%FITENUM - Enum of Fit +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=FitEnum() + +macro=113; Index: /issm/trunk/src/m/enum/GeographyEnum.m =================================================================== --- /issm/trunk/src/m/enum/GeographyEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/GeographyEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=GeographyEnum() -macro=72; +macro=75; Index: /issm/trunk/src/m/enum/GeothermalFluxEnum.m =================================================================== --- /issm/trunk/src/m/enum/GeothermalFluxEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/GeothermalFluxEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=GeothermalFluxEnum() +%GEOTHERMALFLUXENUM - Enum of GeothermalFlux +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=GeothermalFluxEnum() + +macro=106; Index: /issm/trunk/src/m/enum/IceEnum.m =================================================================== --- /issm/trunk/src/m/enum/IceEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/IceEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=IceEnum() -macro=76; +macro=79; Index: /issm/trunk/src/m/enum/IceSheetEnum.m =================================================================== --- /issm/trunk/src/m/enum/IceSheetEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/IceSheetEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=IceSheetEnum() -macro=73; +macro=76; Index: /issm/trunk/src/m/enum/IceShelfEnum.m =================================================================== --- /issm/trunk/src/m/enum/IceShelfEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/IceShelfEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=IceShelfEnum() -macro=74; +macro=77; Index: /issm/trunk/src/m/enum/IntInputEnum.m =================================================================== --- /issm/trunk/src/m/enum/IntInputEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/IntInputEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=IntInputEnum() +%INTINPUTENUM - Enum of IntInput +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=IntInputEnum() + +macro=68; Index: /issm/trunk/src/m/enum/MelangeEnum.m =================================================================== --- /issm/trunk/src/m/enum/MelangeEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/MelangeEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=MelangeEnum() -macro=78; +macro=81; Index: /issm/trunk/src/m/enum/MeltingRateEnum.m =================================================================== --- /issm/trunk/src/m/enum/MeltingRateEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/MeltingRateEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=MeltingRateEnum() +%MELTINGRATEENUM - Enum of MeltingRate +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=MeltingRateEnum() + +macro=104; Index: /issm/trunk/src/m/enum/ParamEnum.m =================================================================== --- /issm/trunk/src/m/enum/ParamEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/ParamEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=ParamEnum() -macro=67; +macro=70; Index: /issm/trunk/src/m/enum/ResultEnum.m =================================================================== --- /issm/trunk/src/m/enum/ResultEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/ResultEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=ResultEnum() -macro=68; +macro=71; Index: /issm/trunk/src/m/enum/RgbEnum.m =================================================================== --- /issm/trunk/src/m/enum/RgbEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/RgbEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=RgbEnum() -macro=69; +macro=72; Index: /issm/trunk/src/m/enum/RheologyBEnum.m =================================================================== --- /issm/trunk/src/m/enum/RheologyBEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/RheologyBEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=RheologyBEnum() +%RHEOLOGYBENUM - Enum of RheologyB +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=RheologyBEnum() + +macro=102; Index: /issm/trunk/src/m/enum/RheologyNEnum.m =================================================================== --- /issm/trunk/src/m/enum/RheologyNEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/RheologyNEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=RheologyNEnum() +%RHEOLOGYNENUM - Enum of RheologyN +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=RheologyNEnum() + +macro=103; Index: /issm/trunk/src/m/enum/SpcEnum.m =================================================================== --- /issm/trunk/src/m/enum/SpcEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/SpcEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=SpcEnum() -macro=70; +macro=73; Index: /issm/trunk/src/m/enum/SurfaceAreaEnum.m =================================================================== --- /issm/trunk/src/m/enum/SurfaceAreaEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/SurfaceAreaEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=SurfaceAreaEnum() +%SURFACEAREAENUM - Enum of SurfaceArea +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=SurfaceAreaEnum() + +macro=111; Index: /issm/trunk/src/m/enum/SurfaceEnum.m =================================================================== --- /issm/trunk/src/m/enum/SurfaceEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/SurfaceEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=SurfaceEnum() +%SURFACEENUM - Enum of Surface +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=SurfaceEnum() + +macro=96; Index: /issm/trunk/src/m/enum/ThicknessEnum.m =================================================================== --- /issm/trunk/src/m/enum/ThicknessEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/ThicknessEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=ThicknessEnum() +%THICKNESSENUM - Enum of Thickness +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=ThicknessEnum() + +macro=95; Index: /issm/trunk/src/m/enum/VxAverageEnum.m =================================================================== --- /issm/trunk/src/m/enum/VxAverageEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VxAverageEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VxAverageEnum() +%VXAVERAGEENUM - Enum of VxAverage +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VxAverageEnum() + +macro=85; Index: /issm/trunk/src/m/enum/VxEnum.m =================================================================== --- /issm/trunk/src/m/enum/VxEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/VxEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=VxEnum() -macro=79; +macro=82; Index: /issm/trunk/src/m/enum/VxObsEnum.m =================================================================== --- /issm/trunk/src/m/enum/VxObsEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VxObsEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VxObsEnum() +%VXOBSENUM - Enum of VxObs +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VxObsEnum() + +macro=88; Index: /issm/trunk/src/m/enum/VxOldEnum.m =================================================================== --- /issm/trunk/src/m/enum/VxOldEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VxOldEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VxOldEnum() +%VXOLDENUM - Enum of VxOld +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VxOldEnum() + +macro=91; Index: /issm/trunk/src/m/enum/VyAverageEnum.m =================================================================== --- /issm/trunk/src/m/enum/VyAverageEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VyAverageEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VyAverageEnum() +%VYAVERAGEENUM - Enum of VyAverage +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VyAverageEnum() + +macro=86; Index: /issm/trunk/src/m/enum/VyEnum.m =================================================================== --- /issm/trunk/src/m/enum/VyEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/VyEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=VyEnum() -macro=80; +macro=83; Index: /issm/trunk/src/m/enum/VyObsEnum.m =================================================================== --- /issm/trunk/src/m/enum/VyObsEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VyObsEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VyObsEnum() +%VYOBSENUM - Enum of VyObs +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VyObsEnum() + +macro=89; Index: /issm/trunk/src/m/enum/VyOldEnum.m =================================================================== --- /issm/trunk/src/m/enum/VyOldEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VyOldEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VyOldEnum() +%VYOLDENUM - Enum of VyOld +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VyOldEnum() + +macro=92; Index: /issm/trunk/src/m/enum/VzAverageEnum.m =================================================================== --- /issm/trunk/src/m/enum/VzAverageEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VzAverageEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VzAverageEnum() +%VZAVERAGEENUM - Enum of VzAverage +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VzAverageEnum() + +macro=87; Index: /issm/trunk/src/m/enum/VzEnum.m =================================================================== --- /issm/trunk/src/m/enum/VzEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VzEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VzEnum() +%VZENUM - Enum of Vz +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VzEnum() + +macro=84; Index: /issm/trunk/src/m/enum/VzObsEnum.m =================================================================== --- /issm/trunk/src/m/enum/VzObsEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VzObsEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VzObsEnum() +%VZOBSENUM - Enum of VzObs +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VzObsEnum() + +macro=90; Index: /issm/trunk/src/m/enum/VzOldEnum.m =================================================================== --- /issm/trunk/src/m/enum/VzOldEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/VzOldEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=VzOldEnum() +%VZOLDENUM - Enum of VzOld +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=VzOldEnum() + +macro=93; Index: /issm/trunk/src/m/enum/WaterEnum.m =================================================================== --- /issm/trunk/src/m/enum/WaterEnum.m (revision 3611) +++ /issm/trunk/src/m/enum/WaterEnum.m (revision 3612) @@ -7,3 +7,3 @@ % macro=WaterEnum() -macro=75; +macro=78; Index: /issm/trunk/src/m/enum/WeightsEnum.m =================================================================== --- /issm/trunk/src/m/enum/WeightsEnum.m (revision 3612) +++ /issm/trunk/src/m/enum/WeightsEnum.m (revision 3612) @@ -0,0 +1,9 @@ +function macro=WeightsEnum() +%WEIGHTSENUM - Enum of Weights +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=WeightsEnum() + +macro=112; Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m (revision 3611) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m (revision 3612) @@ -14,4 +14,5 @@ converged=0; count=1; soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3) (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node + soln(count).u_f= Reducevectorgtof(soln(count).u_g,m.nodesets); @@ -56,5 +57,4 @@ %Update elements with new solution [m.elements]=UpdateInputs(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,soln(count).u_g,analysis_type,sub_analysis_type); - error; %Deal with penalty loads @@ -63,4 +63,5 @@ %penalty constraints [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes,m.vertices,loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + error; displaystring(m.parameters.verbose,'%s%i',' number of unstable constraints: ',num_unstable_constraints); Index: /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp =================================================================== --- /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 3611) +++ /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 3612) @@ -50,4 +50,7 @@ PenaltyConstraintsx(&converged, &num_unstable_constraints, elements,nodes,vertices, loads,materials,parameters,inputs,analysis_type,sub_analysis_type); + elements->Echo(); + mexErrMsgTxt(" "); + /*write output datasets: */ WriteData(LOADS,loads);