Index: /issm/trunk/src/c/Container/Parameters.cpp =================================================================== --- /issm/trunk/src/c/Container/Parameters.cpp (revision 6212) +++ /issm/trunk/src/c/Container/Parameters.cpp (revision 6213) @@ -175,4 +175,31 @@ } /*}}}*/ +/*FUNCTION Parameters::FindParam(int** pintarray,int* pM,int enum_type){{{1*/ +int Parameters::FindParam(int** pintarray,int* pM, int enum_type){ + + /*Go through a dataset, and find a Param* object + *which parameter name is "name" : */ + + vector::iterator object; + Param* param=NULL; + + int found=0; + + for ( object=objects.begin() ; object < objects.end(); object++ ){ + + /*Ok, this object is a parameter, recover it and ask which name it has: */ + param=(Param*)(*object); + + if(param->EnumType()==enum_type){ + /*Ok, this is the one! Recover the value of this parameter: */ + param->GetParameterValue(pintarray,pM); + found=1; + break; + } + } + return found; + +} +/*}}}*/ /*FUNCTION Parameters::FindParam(double** pdoublearray,int* pM,int enum_type){{{1*/ int Parameters::FindParam(double** pdoublearray,int* pM, int enum_type){ Index: /issm/trunk/src/c/Container/Parameters.h =================================================================== --- /issm/trunk/src/c/Container/Parameters.h (revision 6212) +++ /issm/trunk/src/c/Container/Parameters.h (revision 6213) @@ -31,4 +31,5 @@ int FindParam(char** pstring,int enum_type); int FindParam(char*** pstringarray,int* pM,int enum_type); + int FindParam(int** pintarray,int* pM,int enum_type); int FindParam(double** pdoublearray,int* pM,int enum_type); int FindParam(double** pdoublearray,int* pM,int* pN,int enum_type); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 6212) +++ /issm/trunk/src/c/Makefile.am (revision 6213) @@ -383,4 +383,5 @@ ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\ ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\ + ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\ ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\ @@ -952,4 +953,5 @@ ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\ ./modules/ModelProcessorx/Control/CreateParametersControl.cpp\ + ./modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp\ ./modules/ModelProcessorx/Thermal/UpdateElementsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateNodesThermal.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/UpdateElementsBalancedthickness.cpp (revision 6213) @@ -37,4 +37,5 @@ } if(iomodel->control_analysis){ + IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type"); IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs"); IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); @@ -66,3 +67,4 @@ xfree((void**)&iomodel->thickness_obs); xfree((void**)&iomodel->weights); + xfree((void**)&iomodel->control_type); } Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp (revision 6213) @@ -27,12 +27,7 @@ /*What control type?*/ - if (iomodel->control_type!=DragCoefficientEnum && - iomodel->control_type!=RheologyBbarEnum && - iomodel->control_type!=DhDtEnum && - iomodel->control_type!=VxEnum && - iomodel->control_type!=VyEnum - ) ISSMERROR("control_type %s not supported yet!",EnumToString(iomodel->control_type)); - - parameters->AddObject(new IntParam(ControlTypeEnum,iomodel->control_type)); + IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type"); + parameters->AddObject(new IntVecParam(ControlTypeEnum,iomodel->control_type,iomodel->num_control_type)); + xfree((void**)&iomodel->control_type); /*What solution type?*/ Index: /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 6213) +++ /issm/trunk/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp (revision 6213) @@ -0,0 +1,80 @@ +/* + * UpdateElementsAndMaterialsControl: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + /*Intermediary*/ + int i; + int counter; + Element *element = NULL; + Material *material = NULL; + + /*Now, return if no control*/ + if (!iomodel->control_analysis); + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type"); + IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); + IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs"); + IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); + IoModelFetchData(&iomodel->thickness_obs,NULL,NULL,iomodel_handle,"thickness_obs"); + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case DhDtEnum: + IoModelFetchData(&iomodel->dhdt,NULL,NULL,iomodel_handle,"dhdt"); + break; + case VxEnum: + IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); + break; + case VyEnum: + IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); + break; + case DragCoefficientEnum: + IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + break; + case RheologyBbarEnum: + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + break; + default: + ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); + } + } + + /*Update elements: */ + counter=0; + for (i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel); //we need i to index into elements. + + material=(Material*)materials->GetObjectByOffset(counter); + material->Update(i,iomodel); //we need i to index into elements. + counter++; + } + } + + /*Free data: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->weights); + xfree((void**)&iomodel->control_type); + xfree((void**)&iomodel->vx_obs); + xfree((void**)&iomodel->vy_obs); + xfree((void**)&iomodel->thickness_obs); + xfree((void**)&iomodel->dhdt); + xfree((void**)&iomodel->vx); + xfree((void**)&iomodel->vy); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->rheology_B); +} Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 6213) @@ -18,12 +18,14 @@ void CreateDataSets(Elements** pelements,Nodes** pnodes, Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,const int solution_type,const int analysis_type,const int nummodels,int analysis_counter){ - bool continuous=true; - Elements* elements=NULL; + bool continuous = true; + Elements *elements = NULL; + Materials *materials = NULL; /*Create elements, vertices and materials, independent of analysis_type: */ CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,iomodel_handle,nummodels); - /*Recover elements, for future update: */ + /*Recover elements and materials, for future update: */ elements=*pelements; + materials=*pmaterials; /*Now, branch onto analysis dependent model generation: */ @@ -102,4 +104,7 @@ } + /*Update Elements and Materials For Control methods*/ + UpdateElementsAndMaterialsControl(elements,materials,iomodel,iomodel_handle); + /*Generate objects that are not dependent on any analysis_type: */ CreateParameters(pparameters,iomodel,iomodel_handle,solution_type,analysis_type,analysis_counter); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 6213) @@ -40,4 +40,5 @@ IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type"); /*Create elements and materials: */ @@ -61,4 +62,5 @@ xfree((void**)&iomodel->rheology_B); xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->control_type); /*Add new constrant material property tgo materials, at the end: */ @@ -95,4 +97,3 @@ *pvertices=vertices; *pmaterials=materials; - } Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 6213) @@ -47,4 +47,5 @@ } if(iomodel->control_analysis){ + IoModelFetchData(&iomodel->control_type,NULL,NULL,iomodel_handle,"control_type"); IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs"); IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); @@ -87,3 +88,4 @@ xfree((void**)&iomodel->vy_obs); xfree((void**)&iomodel->weights); + xfree((void**)&iomodel->control_type); } Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 6212) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 6213) @@ -21,4 +21,5 @@ void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int solution_type,int analysis_type); void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int solution_type,int analysis_type); +void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel,ConstDataHandle iomodel_handle); /*Creation of fem datasets: specialised drivers: */ Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 6212) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 6213) @@ -56,4 +56,5 @@ virtual void DeleteResults(void)=0; virtual void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type)=0; + virtual void Update(int index, IoModel* iomodel)=0; virtual void UpdateGeometry(void)=0; virtual void InputToResult(int enum_type,int step,double time)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6212) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 6213) @@ -641,4 +641,5 @@ } if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); this->GetDofList1(&doflist1[0]); @@ -661,4 +662,5 @@ } if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); this->GetDofList1(&doflist1[0]); @@ -959,5 +961,6 @@ /*Intermediary*/ - int control_type; + int num_controls; + int* control_type=NULL; Input* input=NULL; double cm_min,cm_max; @@ -966,23 +969,30 @@ this->parameters->FindParam(&cm_min,CmMinEnum); this->parameters->FindParam(&cm_max,CmMaxEnum); - this->parameters->FindParam(&control_type,ControlTypeEnum); - - if(control_type==RheologyBbarEnum){ - if (!IsOnBed()) return; - input=(Input*)matice->inputs->GetInput(RheologyBEnum); ISSMASSERT(input); - } - else{ - input=(Input*)this->inputs->GetInput(control_type); ISSMASSERT(input); - } - - if (input->Enum()!=ControlInputEnum) ISSMERROR("input %s is not a ControlInput",EnumToString(control_type)); - - ((ControlInput*)input)->UpdateValue(scalar); - input->Constrain(cm_min,cm_max); - if (save_parameter) ((ControlInput*)input)->SaveValue(); - - if(control_type==RheologyBbarEnum){ - this->InputExtrude(RheologyBEnum,MaterialsEnum); - } + this->parameters->FindParam(&num_controls,NumControlsEnum); + this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + + for(int i=0;iinputs->GetInput(RheologyBEnum); ISSMASSERT(input); + } + else{ + input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); + } + + if (input->Enum()!=ControlInputEnum) ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); + + ((ControlInput*)input)->UpdateValue(scalar); + input->Constrain(cm_min,cm_max); + if (save_parameter) ((ControlInput*)input)->SaveValue(); + + if(control_type[i]==RheologyBbarEnum){ + this->InputExtrude(RheologyBEnum,MaterialsEnum); + } + } + + /*Clean up and return*/ + xfree((void**)&control_type); } /*}}}*/ @@ -1688,9 +1698,9 @@ } /*}}}*/ -/*FUNCTION Penta::Update {{{1*/ +/*FUNCTION Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type) {{{1*/ void Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type){ /*Intermediaries*/ - IssmInt i; + IssmInt i,j; int penta_type; int penta_node_ids[6]; @@ -1729,89 +1739,6 @@ this->SetHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type - //add as many inputs per element as requested: - if (iomodel->thickness) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); - } - if (iomodel->surface) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); - } - if (iomodel->bed) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); - } - if (iomodel->drag_coefficient) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(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<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); - } - if (iomodel->accumulation_rate) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); - } - if (iomodel->geothermalflux) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); - } - if (iomodel->pressure) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); - } - if (iomodel->temperature) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); - } - if (iomodel->dhdt) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); - } - /*vx,vy and vz: */ - if (iomodel->vx) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs)); - } - if (iomodel->vy) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs)); - } - if (iomodel->vz) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs)); - } - if (iomodel->vx_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs)); - } - if (iomodel->vy_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs)); - } - if (iomodel->vz_obs) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs)); - } - if (iomodel->weights) { - for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new PentaVertexInput(WeightsEnum,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])); + /*Fill with IoModel*/ + this->Update(index,iomodel); /*Defaults if not provided in iomodel*/ @@ -1888,35 +1815,143 @@ } +} +/*}}}*/ +/*FUNCTION Penta::Update(int index,IoModel* iomodel) {{{1*/ +void Penta::Update(int index,IoModel* iomodel){ + + /*Intermediaries*/ + IssmInt i,j; + int penta_vertex_ids[6]; + double nodeinputs[6]; + + /*Checks if debuging*/ + /*{{{2*/ + ISSMASSERT(iomodel->elements); + /*}}}*/ + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<6;i++){ + penta_vertex_ids[i]=(int)iomodel->elements[6*index+i]; //ids for vertices are in the elements array from Matlab + } + + //add as many inputs per element as requested: + if (iomodel->thickness) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->thickness[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(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<6;i++)nodeinputs[i]=iomodel->melting_rate[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(GeothermalFluxEnum,nodeinputs)); + } + if (iomodel->pressure) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->pressure[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(PressureEnum,nodeinputs)); + } + if (iomodel->temperature) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); + } + if (iomodel->dhdt) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(DhDtEnum,nodeinputs)); + } + /*vx,vy and vz: */ + if (iomodel->vx) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VxOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVxEnum,nodeinputs)); + } + if (iomodel->vy) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VyOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVyEnum,nodeinputs)); + } + if (iomodel->vz) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VzOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new PentaVertexInput(QmuVzEnum,nodeinputs)); + } + if (iomodel->vx_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VxObsEnum,nodeinputs)); + } + if (iomodel->vy_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VyObsEnum,nodeinputs)); + } + if (iomodel->vz_obs) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->vz_obs[penta_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new PentaVertexInput(VzObsEnum,nodeinputs)); + } + if (iomodel->weights) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->weights[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(WeightsEnum,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])); + /*Control Inputs*/ - if (iomodel->control_analysis){ - switch(iomodel->control_type){ - case DhDtEnum: - if (iomodel->dhdt){ - for(i=0;i<6;i++)nodeinputs[i]=iomodel->dhdt[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case VxEnum: - if (iomodel->vx){ - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vx[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case VyEnum: - if (iomodel->vy){ - for(i=0;i<6;i++)nodeinputs[i]=iomodel->vy[penta_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case DragCoefficientEnum: - if (iomodel->drag_coefficient){ - for(i=0;i<6;i++)nodeinputs[i]=iomodel->drag_coefficient[penta_vertex_ids[i]-1]; - this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs)); - } - break; - case RheologyBbarEnum: - /*Matice will take care of it*/ break; - default: - ISSMERROR("Control %s not implemented yet",EnumToString(iomodel->control_type)); + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case DhDtEnum: + if (iomodel->dhdt){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->dhdt[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(DhDtEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case VxEnum: + if (iomodel->vx){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->vx[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VxEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case VyEnum: + if (iomodel->vy){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->vy[penta_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VyEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case DragCoefficientEnum: + if (iomodel->drag_coefficient){ + for(j=0;j<6;j++)nodeinputs[j]=iomodel->drag_coefficient[penta_vertex_ids[j]-1]; + this->inputs->AddInput(new ControlInput(DragCoefficientEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + case RheologyBbarEnum: + /*Matice will take care of it*/ break; + default: + ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); + } } } Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 6212) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 6213) @@ -119,4 +119,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); + void Update(int index, IoModel* iomodel); void UpdateGeometry(void); double TimeAdapt(); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6212) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6213) @@ -641,4 +641,5 @@ } if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); this->GetDofList1(&doflist1[0]); @@ -661,4 +662,5 @@ } if (!input) ISSMERROR("Input %s not found",EnumToString(enum_type)); + if (input->Enum()!=ControlInputEnum) ISSMERROR("Input %s is not a ControlInput",EnumToString(enum_type)); this->GetDofList1(&doflist1[0]); @@ -677,5 +679,6 @@ /* Intermediaries */ int i,j,ig; - int control_type; + int num_controls; + int *control_type=NULL; double Jelem = 0; double cm_noisedmp; @@ -686,5 +689,6 @@ /*retrieve parameters and inputs*/ - this->parameters->FindParam(&control_type,ControlTypeEnum); + this->parameters->FindParam(&num_controls,NumControlsEnum); + this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); this->parameters->FindParam(&cm_noisedmp,CmNoiseDmpEnum); @@ -694,42 +698,46 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(2); - for (ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); - - /*Add Tikhonov regularization term to misfit: - * - * WARNING: for now, the regularization is only taken into account by the gradient - * the misfit reflects the response only! - * - * */ - if (control_type==DragCoefficientEnum){ - Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); - drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); - //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss); + + /*Add Tikhonov regularization term to misfit: + * + * WARNING: for now, the regularization is only taken into account by the gradient + * the misfit reflects the response only! + * + * */ + if (control_type[i]==DragCoefficientEnum){ + Input* drag_input=inputs->GetInput(DragCoefficientEnum); ISSMASSERT(drag_input); + drag_input->GetParameterDerivativeValue(&dk[0],&xyz_list[0][0],gauss); + //Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss->weight; + } + else if (control_type[i]==RheologyBbarEnum){ + //nothing + } + else if (control_type[i]==DhDtEnum){ + //nothing + } + else if (control_type[i]==VxEnum){ + //nothing + } + else if (control_type[i]==VyEnum){ + //nothing + } + else{ + ISSMERROR("unsupported control type: %s",EnumToString(control_type[i])); + } } - else if (control_type==RheologyBbarEnum){ - //nothing - } - else if (control_type==DhDtEnum){ - //nothing - } - else if (control_type==VxEnum){ - //nothing - } - else if (control_type==VyEnum){ - //nothing - } - else{ - ISSMERROR("unsupported control type: %s",EnumToString(control_type)); - } - } - - /*clean-up and return: */ - delete gauss; + + delete gauss; + } + + /*Clean up and return*/ + xfree((void**)&control_type); return Jelem; } @@ -1187,5 +1195,6 @@ /*Intermediary*/ - int control_type; + int num_controls; + int* control_type=NULL; Input* input=NULL; double cm_min,cm_max; @@ -1194,20 +1203,28 @@ this->parameters->FindParam(&cm_min,CmMinEnum); this->parameters->FindParam(&cm_max,CmMaxEnum); - this->parameters->FindParam(&control_type,ControlTypeEnum); - - if(control_type==RheologyBbarEnum){ - input=(Input*)matice->inputs->GetInput(control_type); ISSMASSERT(input); - } - else{ - input=(Input*)this->inputs->GetInput(control_type); ISSMASSERT(input); - } - - if (input->Enum()!=ControlInputEnum){ - ISSMERROR("input %s is not a ControlInput",EnumToString(control_type)); - } - - ((ControlInput*)input)->UpdateValue(scalar); - input->Constrain(cm_min,cm_max); - if (save_parameter) ((ControlInput*)input)->SaveValue(); + this->parameters->FindParam(&num_controls,NumControlsEnum); + this->parameters->FindParam(&control_type,NULL,ControlTypeEnum); + + for(int i=0;iinputs->GetInput(control_type[i]); ISSMASSERT(input); + } + else{ + input=(Input*)this->inputs->GetInput(control_type[i]); ISSMASSERT(input); + } + + if (input->Enum()!=ControlInputEnum){ + ISSMERROR("input %s is not a ControlInput",EnumToString(control_type[i])); + } + + ((ControlInput*)input)->UpdateValue(scalar); + input->Constrain(cm_min,cm_max); + if (save_parameter) ((ControlInput*)input)->SaveValue(); + + } + + /*Clean up and return*/ + xfree((void**)&control_type); } /*}}}*/ @@ -2117,9 +2134,9 @@ } /*}}}*/ -/*FUNCTION Tria::Update{{{1*/ +/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index /*Intermediaries*/ - int i; + int i,j; int tria_node_ids[3]; int tria_vertex_ids[3]; @@ -2164,93 +2181,7 @@ /*hooks: */ this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type - - /*add as many inputs per element as requested:*/ - if (iomodel->thickness) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs)); - } - if (iomodel->surface) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs)); - } - if (iomodel->bed) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_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_vertex_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->thickness_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_obs[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(ThicknessObsEnum,nodeinputs)); - } - if (iomodel->melting_rate) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs)); - } - if (iomodel->accumulation_rate) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs)); - } - if (iomodel->geothermalflux) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs)); - } - if (iomodel->dhdt){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs)); - } - if (iomodel->pressure){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); - } - if (iomodel->temperature) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs)); - } - /*vx,vy and vz: */ - if (iomodel->vx) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); - } - if (iomodel->vy) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); - } - if (iomodel->vz) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); - this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); - if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); - } - if (iomodel->vx_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs)); - } - if (iomodel->vy_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs)); - } - if (iomodel->vz_obs) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs)); - } - if (iomodel->weights) { - for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new TriaVertexInput(WeightsEnum,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])); + + /*Fill with IoModel*/ + this->Update(index,iomodel); /*Defaults if not provided in iomodel*/ @@ -2296,42 +2227,150 @@ } - /*Control Inputs*/ - if (iomodel->control_analysis){ - switch(iomodel->control_type){ - case DhDtEnum: - if (iomodel->dhdt){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(DhDtEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - case VxEnum: - if (iomodel->vx){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - case VyEnum: - if (iomodel->vy){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts; - this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - case DragCoefficientEnum: - if (iomodel->drag_coefficient){ - for(i=0;i<3;i++)nodeinputs[i]=iomodel->drag_coefficient[tria_vertex_ids[i]-1]; - this->inputs->AddInput(new ControlInput(DragCoefficientEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - case RheologyBbarEnum: - /*Matice will take care of it*/ break; - default: - ISSMERROR("Control %s not implemented yet",EnumToString(iomodel->control_type)); - } - - } - //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. this->parameters=NULL; - +} +/*}}}*/ +/*FUNCTION Tria::Update(int index, IoModel* iomodel){{{1*/ +void Tria::Update(int index, IoModel* iomodel){ //i is the element index + + /*Intermediaries*/ + int i,j; + int tria_vertex_ids[3]; + double nodeinputs[3]; + + /*Checks if debuging*/ + /*{{{2*/ + ISSMASSERT(iomodel->elements); + /*}}}*/ + + /*Recover vertices ids needed to initialize inputs*/ + for(i=0;i<3;i++){ + tria_vertex_ids[i]=(int)iomodel->elements[3*index+i]; //ids for vertices are in the elements array from Matlab + } + + /*add as many inputs per element as requested:*/ + if (iomodel->thickness) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->surface[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_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_vertex_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->thickness_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->thickness_obs[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(ThicknessObsEnum,nodeinputs)); + } + if (iomodel->melting_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->melting_rate[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->geothermalflux[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(GeothermalFluxEnum,nodeinputs)); + } + if (iomodel->dhdt){ + for(i=0;i<3;i++)nodeinputs[i]=iomodel->dhdt[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(DhDtEnum,nodeinputs)); + } + if (iomodel->pressure){ + for(i=0;i<3;i++)nodeinputs[i]=iomodel->pressure[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(PressureEnum,nodeinputs)); + } + if (iomodel->temperature) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->temperature[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(TemperatureEnum,nodeinputs)); + } + /*vx,vy and vz: */ + if (iomodel->vx) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VxEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VxOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVxEnum,nodeinputs)); + } + if (iomodel->vy) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VyEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VyOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVyEnum,nodeinputs)); + } + if (iomodel->vz) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VzEnum,nodeinputs)); + this->inputs->AddInput(new TriaVertexInput(VzOldEnum,nodeinputs)); + if(iomodel->qmu_analysis)this->inputs->AddInput(new TriaVertexInput(QmuVzEnum,nodeinputs)); + } + if (iomodel->vx_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vx_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VxObsEnum,nodeinputs)); + } + if (iomodel->vy_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vy_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VyObsEnum,nodeinputs)); + } + if (iomodel->vz_obs) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->vz_obs[tria_vertex_ids[i]-1]/iomodel->yts; + this->inputs->AddInput(new TriaVertexInput(VzObsEnum,nodeinputs)); + } + if (iomodel->weights) { + for(i=0;i<3;i++)nodeinputs[i]=iomodel->weights[tria_vertex_ids[i]-1]; + this->inputs->AddInput(new TriaVertexInput(WeightsEnum,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])); + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case DhDtEnum: + if (iomodel->dhdt){ + for(j=0;j<3;j++)nodeinputs[j]=iomodel->dhdt[tria_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(DhDtEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + case VxEnum: + if (iomodel->vx){ + for(j=0;j<3;j++)nodeinputs[j]=iomodel->vx[tria_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VxEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + case VyEnum: + if (iomodel->vy){ + for(j=0;j<3;j++)nodeinputs[j]=iomodel->vy[tria_vertex_ids[j]-1]/iomodel->yts; + this->inputs->AddInput(new ControlInput(VyEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + case DragCoefficientEnum: + if (iomodel->drag_coefficient){ + for(j=0;j<3;j++)nodeinputs[j]=iomodel->drag_coefficient[tria_vertex_ids[j]-1]; + this->inputs->AddInput(new ControlInput(DragCoefficientEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + case RheologyBbarEnum: + /*Matice will take care of it*/ break; + default: + ISSMERROR("Control %s not implemented yet",EnumToString((int)iomodel->control_type[i])); + } + } + } } /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6212) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6213) @@ -122,4 +122,5 @@ double SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); + void Update(int index, IoModel* iomodel); void UpdateGeometry(void); double TimeAdapt(); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 6212) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 6213) @@ -84,4 +84,5 @@ xfree((void**)&this->rheology_B); xfree((void**)&this->rheology_n); + xfree((void**)&this->control_type); xfree((void**)&this->cm_responses); xfree((void**)&this->weights); @@ -152,5 +153,4 @@ /*Get control parameters: */ IoModelFetchData(&this->num_control_type,iomodel_handle,"num_control_type"); - IoModelFetchData(&this->control_type,iomodel_handle,"control_type"); /*!Get solution parameters: */ @@ -295,4 +295,5 @@ /*!solution parameters: */ + this->control_type=NULL; this->cm_responses=NULL; this->weights=NULL; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 6212) +++ /issm/trunk/src/c/objects/IoModel.h (revision 6213) @@ -122,5 +122,5 @@ /*control methods: */ int num_control_type; - int control_type; + double* control_type; /*solution parameters: */ Index: /issm/trunk/src/c/objects/Materials/Material.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Material.h (revision 6212) +++ /issm/trunk/src/c/objects/Materials/Material.h (revision 6213) @@ -22,4 +22,5 @@ virtual void InputDuplicate(int original_enum,int new_enum)=0; virtual void Configure(Elements* elements)=0; + virtual void Update(int index, IoModel* iomodel)=0; }; Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 6212) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 6213) @@ -57,14 +57,14 @@ /*Control Inputs*/ - if (iomodel->control_analysis){ - switch(iomodel->control_type){ - case RheologyBbarEnum: - if (iomodel->rheology_B){ - for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; - this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); - } - break; - default: - /*Nothing*/; + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + } } } @@ -91,12 +91,14 @@ /*Control Inputs*/ - if (iomodel->control_analysis){ - switch(iomodel->control_type){ - case RheologyBbarEnum: - if (iomodel->rheology_B){ - for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; - this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); - } - break; + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(i=0;irheology_B[int(iomodel->elements[num_vertices*index+i]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + } } } @@ -246,4 +248,59 @@ return matice; +} +/*}}}*/ +/*FUNCTION Matice::Update{{{1*/ +void Matice::Update(int index, IoModel* iomodel){ + + int i,j; + + /*if 2d*/ + if(iomodel->dim==2){ + + /*Intermediaries*/ + const int num_vertices = 3; //Tria has 3 vertices + double nodeinputs[num_vertices]; + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBbarEnum,TriaVertexInputEnum,nodeinputs)); + } + break; + } + } + } + } + + /*if 3d*/ + else if(iomodel->dim==3){ + + /*Intermediaries*/ + const int num_vertices = 6; //Penta has 6 vertices + double nodeinputs[num_vertices]; + + /*Control Inputs*/ + if (iomodel->control_analysis && iomodel->control_type){ + for(i=0;inum_control_type;i++){ + switch((int)iomodel->control_type[i]){ + case RheologyBbarEnum: + if (iomodel->rheology_B){ + for(j=0;jrheology_B[int(iomodel->elements[num_vertices*index+j]-1)]; + this->inputs->AddInput(new ControlInput(RheologyBEnum,PentaVertexInputEnum,nodeinputs)); + } + break; + } + } + } + } + else{ + ISSMERROR(" Mesh type not supported yet!"); + } + + return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.h (revision 6212) +++ /issm/trunk/src/c/objects/Materials/Matice.h (revision 6213) @@ -52,4 +52,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); + void Update(int index, IoModel* iomodel); /*}}}*/ /*Material virtual functions resolution: {{{1*/ Index: /issm/trunk/src/c/objects/Materials/Matpar.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.h (revision 6212) +++ /issm/trunk/src/c/objects/Materials/Matpar.h (revision 6213) @@ -56,4 +56,5 @@ void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromSolution(double* solution); + void Update(int index, IoModel* iomodel){ISSMERROR("not implemented yet");}; /*}}}*/ /*Material virtual functions resolution: {{{1*/ Index: /issm/trunk/src/c/objects/Params/DoubleParam.cpp =================================================================== --- /issm/trunk/src/c/objects/Params/DoubleParam.cpp (revision 6212) +++ /issm/trunk/src/c/objects/Params/DoubleParam.cpp (revision 6213) @@ -129,5 +129,5 @@ } /*}}}*/ -/*FUNCTION DoubleParam::GetParameterValue{{{1*/ +/*FUNCTION DoubleParam::GetParameterValue(int* pinteger){{{1*/ void DoubleParam::GetParameterValue(int* pinteger){ #ifdef _SERIAL_ @@ -138,5 +138,5 @@ } /*}}}*/ -/*FUNCTION DoubleParam::GetParameterValue{{{1*/ +/*FUNCTION DoubleParam::GetParameterValue(bool* pbool){{{1*/ void DoubleParam::GetParameterValue(bool* pbool){ #ifdef _SERIAL_ @@ -148,4 +148,20 @@ #else ISSMERROR("Double param of enum %i (%s) cannot return an bool",enum_type,EnumToString(enum_type)); +#endif +} +/*}}}*/ +/*FUNCTION DoubleParam::GetParameterValue(int** pintarray,int* pM){{{1*/ +void DoubleParam::GetParameterValue(int** pintarray,int* pM){ +#ifdef _SERIAL_ + int* output=NULL; + + output=(int*)xmalloc(1*sizeof(int)); + *output=(int)value; + + /*Assign output pointers:*/ + if(pM) *pM=1; + *pintarray=output; +#else + ISSMERROR("Double param of enum %i (%s) cannot return an array of integers",enum_type,EnumToString(enum_type)); #endif } Index: /issm/trunk/src/c/objects/Params/DoubleParam.h =================================================================== --- /issm/trunk/src/c/objects/Params/DoubleParam.h (revision 6212) +++ /issm/trunk/src/c/objects/Params/DoubleParam.h (revision 6213) @@ -52,5 +52,5 @@ void GetParameterValue(bool* pbool); void GetParameterValue(int* pinteger); - void GetParameterValue(int** pintarray,int* pM){ISSMERROR("Double param of enum %i (%s) cannot return an array of integers",enum_type,EnumToString(enum_type));} + void GetParameterValue(int** pintarray,int* pM); void GetParameterValue(double* pdouble){*pdouble=value;} void GetParameterValue(char** pstring){ISSMERROR("Double param of enum %i (%s) cannot return a string",enum_type,EnumToString(enum_type));} Index: /issm/trunk/src/c/objects/Params/IntVecParam.cpp =================================================================== --- /issm/trunk/src/c/objects/Params/IntVecParam.cpp (revision 6212) +++ /issm/trunk/src/c/objects/Params/IntVecParam.cpp (revision 6213) @@ -36,4 +36,14 @@ } /*}}}*/ +/*FUNCTION IntVecParam::IntVecParam(int enum_type,double* values,int M){{{1*/ +IntVecParam::IntVecParam(int in_enum_type,double* in_values, int in_M){ + + enum_type=in_enum_type; + M=in_M; + + values=(int*)xmalloc(M*sizeof(int)); + for(int i=0;iparameters->FindParam(&nsteps,NStepsEnum); - femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&num_controls,NumControlsEnum); + femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); femmodel->parameters->FindParam(&optscal,NULL,OptScalEnum); @@ -122,11 +124,12 @@ /*some results not computed by steadystate_core or diagnostic_core: */ if(!qmu_analysis){ //do not save this if we are running the control core from a qmu run! - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself! + for(i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]); femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0)); - femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0)); + //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0)); } cleanup_and_return: /*Free ressources: */ + xfree((void**)&control_type); xfree((void**)&responses); xfree((void**)&optscal); Index: /issm/trunk/src/c/solutions/controlrestart.cpp =================================================================== --- /issm/trunk/src/c/solutions/controlrestart.cpp (revision 6212) +++ /issm/trunk/src/c/solutions/controlrestart.cpp (revision 6213) @@ -9,10 +9,12 @@ void controlrestart(FemModel* femmodel,double* J){ - int control_type; + int num_controls; + int* control_type = NULL; int nsteps; bool qmu_analysis=true; /*retrieve output file name: */ - femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&num_controls,NumControlsEnum); + femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); femmodel->parameters->FindParam(&nsteps,NStepsEnum); femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); @@ -22,7 +24,7 @@ if(!qmu_analysis){ /*we essentially want J and the parameter: */ - InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself! + for(int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]); femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0)); - femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0)); + //femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumToString(control_type),1,0)); /*write to disk: */ @@ -30,3 +32,5 @@ } + /*Clean up and return*/ + xfree((void**)&control_type); } Index: /issm/trunk/src/c/solutions/gradient_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6212) +++ /issm/trunk/src/c/solutions/gradient_core.cpp (revision 6213) @@ -17,5 +17,6 @@ /*parameters: */ bool control_steady; - int control_type; + int num_controls; + int* control_type=NULL; int verbose; @@ -27,29 +28,35 @@ /*retrieve parameters:*/ femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum); - femmodel->parameters->FindParam(&control_type,ControlTypeEnum); + femmodel->parameters->FindParam(&num_controls,NumControlsEnum); + femmodel->parameters->FindParam(&control_type,NULL,ControlTypeEnum); femmodel->parameters->FindParam(&verbose,VerboseEnum); - if(verbose)_printf_("%s\n"," compute gradient:"); - Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type); + for (int i=0;i0 && search_scalar==0){ - _printf_("%s"," orthogonalization...\n"); - if(verbose)_printf_("%s\n"," retrieve old gradient:"); - ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type); - Orthx(&new_gradient,gradient,old_gradient); - VecFree(&old_gradient); + if(verbose)_printf_(" compute gradient of J with respect to %s\n",EnumToString(control_type[i])); + Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type[i]); + + if(control_steady)diagnostic_core(femmodel); + + if (step>0 && search_scalar==0){ + _printf_("%s"," orthogonalization...\n"); + ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,control_type[i]); + Orthx(&new_gradient,gradient,old_gradient); + VecFree(&old_gradient); + } + else{ + _printf_("%s"," normalizing directions...\n"); + Orthx(&new_gradient,gradient,NULL); + } + VecFree(&gradient); + + /*plug back into inputs: */ + ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type[i],new_gradient); + + /*Free ressources and return:*/ + VecFree(&new_gradient); } - else{ - _printf_("%s"," normalizing directions...\n"); - Orthx(&new_gradient,gradient,NULL); - } - VecFree(&gradient); - /*plug back into inputs: */ - ControlInputSetGradientx(femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials, femmodel->parameters,control_type,new_gradient); - - /*Free ressources and return:*/ - VecFree(&new_gradient); + /*Clean up and return*/ + xfree((void**)&control_type); } Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 6212) +++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp (revision 6213) @@ -31,6 +31,5 @@ int n; double *optscal = NULL; - double *responses =NULL; - int control_type; + double *responses = NULL; int solution_type; int analysis_type; @@ -47,5 +46,4 @@ femmodel->parameters->FindParam(&responses,NULL,CmResponsesEnum); femmodel->parameters->FindParam(&isstokes,IsStokesEnum); - femmodel->parameters->FindParam(&control_type,ControlTypeEnum); femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);