Index: /issm/trunk/src/c/objects/Beam.cpp =================================================================== --- /issm/trunk/src/c/objects/Beam.cpp (revision 3619) +++ /issm/trunk/src/c/objects/Beam.cpp (revision 3620) @@ -11,4 +11,5 @@ #include "stdio.h" #include "./Beam.h" +#include "./BeamVertexInput.h" #include #include "../EnumDefinitions/EnumDefinitions.h" @@ -21,40 +22,45 @@ /*FUNCTION Beam::Beam(){{{1*/ Beam::Beam(){ + this->inputs=NULL; return; } /*}}}*/ /*FUNCTION Beam::Beam(int id, int* node_ids, int matice_id, int matpar_id, int numpar_id, ElementProperties* properties){{{1*/ -Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id, ElementProperties* beamproperties): +Beam::Beam(int beam_id,int* beam_node_ids, int beam_matice_id, int beam_matpar_id, int beam_numpar_id): hnodes(beam_node_ids,2), hmatice(&beam_matice_id,1), hmatpar(&beam_matpar_id,1), - hnumpar(&beam_numpar_id,1), - properties(beamproperties) + hnumpar(&beam_numpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=beam_id; - - return; + this->inputs=new Inputs(); } /*}}}*/ /*FUNCTION Beam::Beam(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, Hook* hnumpar, ElementProperties* properties){{{1*/ -Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, ElementProperties* beam_properties): +Beam::Beam(int beam_id,Hook* beam_hnodes, Hook* beam_hmatice, Hook* beam_hmatpar, Hook* beam_hnumpar, Inputs* beam_inputs): hnodes(beam_hnodes), hmatice(beam_hmatice), hmatpar(beam_hmatpar), - hnumpar(beam_hnumpar), - properties(beam_properties) + hnumpar(beam_hnumpar) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=beam_id; - + if(beam_inputs){ + this->inputs=(Inputs*)beam_inputs->Copy(); + } + else{ + this->inputs=new Inputs(); + } return; } /*}}}*/ /*FUNCTION Beam::Beam(int i, IoModel* iomodel){{{1*/ -Beam::Beam(int i,IoModel* iomodel){ - +Beam::Beam(int index,IoModel* iomodel){ + + + int i; /*beam constructor input: */ @@ -64,21 +70,16 @@ int beam_numpar_id; int beam_node_ids[2]; - double beam_h[2]; - double beam_s[2]; - double beam_b[2]; - double beam_k[2]; - int beam_onbed; - - + double nodeinputs[2]; + /*id: */ - beam_id=i+1; + beam_id=index+1; this->id=beam_id; /*hooks: */ - beam_matice_id=i+1; //refers to the corresponding material property card + beam_matice_id=index+1; //refers to the corresponding material property card beam_matpar_id=iomodel->numberofvertices2d*(iomodel->numlayers-1)+1;//refers to the corresponding matpar property card beam_numpar_id=1; - beam_node_ids[0]=i+1; - beam_node_ids[1]=(int)iomodel->uppernodes[i]; //grid that lays right on top + beam_node_ids[0]=index+1; + beam_node_ids[1]=(int)iomodel->uppernodes[index]; //grid that lays right on top this->hnodes.Init(beam_node_ids,2); @@ -87,20 +88,34 @@ this->hnumpar.Init(&beam_numpar_id,1); - /*properties: */ - beam_h[0]=iomodel->thickness[i]; - beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)]; - beam_s[0]=iomodel->surface[i]; - beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)]; - beam_b[0]=iomodel->bed[i]; - beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)]; - beam_k[0]=iomodel->drag[i]; - beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)]; - beam_onbed=(bool)iomodel->gridonbed[i]; - - this->properties.Init(2,beam_h, beam_s, beam_b, beam_k,NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, beam_onbed, UNDEF, UNDEF, UNDEF, UNDEF); + //intialize inputs, and add as many inputs per element as requested: + this->inputs=new Inputs(); + + if (iomodel->thickness) { + nodeinputs[0]=iomodel->thickness[index]; + nodeinputs[1]=iomodel->thickness[(int)(iomodel->uppernodes[index]-1)]; + this->inputs->AddInput(new BeamVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + nodeinputs[0]=iomodel->surface[index]; + nodeinputs[1]=iomodel->surface[(int)(iomodel->uppernodes[index]-1)]; + this->inputs->AddInput(new BeamVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + nodeinputs[0]=iomodel->bed[index]; + nodeinputs[1]=iomodel->bed[(int)(iomodel->uppernodes[index]-1)]; + this->inputs->AddInput(new BeamVertexInput(BedEnum,nodeinputs)); + } + if (iomodel->drag_coefficient) { + nodeinputs[0]=iomodel->drag_coefficient[index]; + nodeinputs[1]=iomodel->drag_coefficient[(int)(iomodel->uppernodes[index]-1)]; + this->inputs->AddInput(new BeamVertexInput(DragCoefficientEnum,nodeinputs)); + } + if (iomodel->gridonbed) this->inputs->AddInput(new BoolInput(ElementOnBedEnum,(IssmBool)iomodel->gridonbed[index])); + } /*}}}*/ /*FUNCTION Beam::~Beam(){{{1*/ Beam::~Beam(){ + delete inputs; return; } @@ -136,5 +151,5 @@ Object* Beam::copy() { - return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties); + return new Beam(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); } @@ -149,5 +164,6 @@ hmatpar.DeepEcho(); hnumpar.DeepEcho(); - properties.DeepEcho(); + printf(" inputs\n"); + inputs->DeepEcho(); return; @@ -173,6 +189,6 @@ hnumpar.Demarshall(&marshalled_dataset); - /*demarshall properties: */ - properties.Demarshall(&marshalled_dataset); + /*demarshall inputs: */ + inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); /*return: */ @@ -190,5 +206,6 @@ hmatpar.Echo(); hnumpar.Echo(); - properties.Echo(); + printf(" inputs\n"); + inputs->Echo(); return; @@ -200,4 +217,6 @@ char* marshalled_dataset=NULL; int enum_type=0; + char* marshalled_inputs=NULL; + int marshalled_inputs_size; /*recover marshalled_dataset: */ @@ -219,6 +238,11 @@ hnumpar.Marshall(&marshalled_dataset); - /*Marshall properties: */ - properties.Marshall(&marshalled_dataset); + /*Marshall inputs: */ + marshalled_inputs_size=inputs->MarshallSize(); + marshalled_inputs=inputs->Marshall(); + memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); + marshalled_dataset+=marshalled_inputs_size; + + xfree((void**)&marshalled_inputs); *pmarshalled_dataset=marshalled_dataset; @@ -234,5 +258,5 @@ +hmatpar.MarshallSize() +hnumpar.MarshallSize() - +properties.MarshallSize() + +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type } @@ -259,6 +283,8 @@ int doflist[numgrids]; double pressure[numgrids]; + double surface[numgrids]; double rho_ice,g; double xyz_list[numgrids][3]; + double gauss[numgrids][numgrids]={{1,0},{0,1}}; /*dynamic objects pointed to by hooks: */ @@ -286,8 +312,11 @@ rho_ice=matpar->GetRhoIce(); g=matpar->GetG(); + + /*recover value of surface at gauss points (0,1) and (1,0): */ + inputs->GetParameterValues(&surface[0],&gauss[0][0],2,SurfaceEnum); for(i=0;iproperties.s[i]-xyz_list[i][2]); - } - + pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); + } + /*plug local pressure values into global pressure vector: */ VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES); @@ -303,5 +332,5 @@ /*}}}*/ /*FUNCTION Beam::CostFunction{{{1*/ -double Beam::CostFunction(void*,int,int){ +double Beam::CostFunction(int,int){ ISSMERROR(" not supported yet!"); } @@ -316,5 +345,5 @@ if (sub_analysis_type==HutterAnalysisEnum) { - CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type); } else @@ -338,8 +367,12 @@ double Ke_gg[numdofs][numdofs]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}}; int numberofdofspernode; - + bool onbed; + + + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + GetDofList(&doflist[0],&numberofdofspernode); - if (this->properties.onbed){ + if (onbed){ Ke_gg[0][0]=1; Ke_gg[1][1]=1; @@ -367,5 +400,5 @@ if (analysis_type==DiagnosticAnalysisEnum) { if (sub_analysis_type==HutterAnalysisEnum) { - CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type); } else @@ -403,4 +436,5 @@ double* gauss_weights=NULL; double gauss_weight; + double gauss1[2]={1,0}; int ig; double l1l2[2]; @@ -412,4 +446,7 @@ double Jdet; double ub,vb; + double surface,thickness; + + bool onbed; /*dynamic objects pointed to by hooks: */ @@ -419,9 +456,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -433,4 +465,7 @@ GetDofList(&doflist[0],&numberofdofspernode); + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*recover parameters: */ rho_ice=matpar->GetRhoIce(); @@ -440,9 +475,8 @@ //recover extra inputs - found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,1,(void**)&nodes[0]); //only recover from fist node, second node will have the same slope - if(!found)ISSMERROR(" surfaceslopex missing from inputs!"); - - found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,1,(void**)&nodes[0]);//only recover from fist node, second node will have the same slope - if(!found)ISSMERROR(" surfaceslopey missing from inputs!"); + inputs->GetParameterValue(&slope[0],&gauss1[0],SurfaceSlopexEnum); + inputs->GetParameterValue(&slope[1],&gauss1[0],SurfaceSlopeyEnum); + inputs->GetParameterValue(&surface,&gauss1[0],SurfaceEnum); + inputs->GetParameterValue(&thickness,&gauss1[0],ThicknessEnum); //Get all element grid data: @@ -474,5 +508,5 @@ for(j=0;jproperties.s[0]-z_g)/B,n)*slope[j]*Jdet*gauss_weight; + pe_g_gaussian[NDOF2+j]=constant_part*pow((surface-z_g)/B,n)*slope[j]*Jdet*gauss_weight; } @@ -484,8 +518,8 @@ //Deal with lower surface - if (this->properties.onbed){ + if (onbed){ //compute ub - constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*this->properties.h[0]; + constant_part=-1.58*pow((double)10.0,-(double)10.0)*rho_ice*gravity*thickness; ub=constant_part*slope[0]; vb=constant_part*slope[1]; @@ -505,5 +539,5 @@ /*}}}*/ /*FUNCTION Beam::Du{{{1*/ -void Beam::Du(_p_Vec*,void*,int,int){ +void Beam::Du(Vec,int,int){ ISSMERROR(" not supported yet!"); } @@ -654,15 +688,15 @@ /*}}}*/ /*FUNCTION Beam::Gradj{{{1*/ -void Beam::Gradj(_p_Vec*, void*, int, int,char*){ +void Beam::Gradj(Vec, int, int,char*){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjB{{{1*/ -void Beam::GradjB(_p_Vec*, void*, int, int){ +void Beam::GradjB(Vec, int, int){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Beam::GradjDrag{{{1*/ -void Beam::GradjDrag(_p_Vec*, void*, int,int ){ +void Beam::GradjDrag(Vec, int,int ){ ISSMERROR(" not supported yet!"); } @@ -674,5 +708,5 @@ /*}}}*/ /*FUNCTION Beam::Misfit{{{1*/ -double Beam::Misfit(void*,int,int){ +double Beam::Misfit(int,int){ ISSMERROR(" not supported yet!"); } @@ -685,5 +719,5 @@ /*}}}*/ /*FUNCTION Beam::SurfaceArea{{{1*/ -double Beam::SurfaceArea(void*,int,int){ +double Beam::SurfaceArea(int,int){ ISSMERROR(" not supported yet!"); } Index: /issm/trunk/src/c/objects/Beam.h =================================================================== --- /issm/trunk/src/c/objects/Beam.h (revision 3619) +++ /issm/trunk/src/c/objects/Beam.h (revision 3620) @@ -12,4 +12,5 @@ #include "./Matpar.h" #include "../ModelProcessorx/IoModel.h" +#include "../DataSet/Inputs.h" #include "./Hook.h" @@ -30,5 +31,4 @@ Hook hmatpar; //hook to 1 matpar Hook hnumpar; //hook to 1 numpar - Inputs* inputs; @@ -52,5 +52,5 @@ int GetId(); int MyRank(); - void Configure(void* loads,void* nodes,void* materials,void* parameters); + void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); Object* copy(); void SetClone(int* minranks); @@ -78,11 +78,11 @@ void GetBedList(double*); void GetThicknessList(double* thickness_list); - void Du(_p_Vec*,void*, int,int); - void Gradj(_p_Vec*, void*, int, int,char*); - void GradjDrag(_p_Vec*, void*, int,int ); - void GradjB(_p_Vec*, void*, int,int ); - double Misfit(void*,int,int); - double SurfaceArea(void*,int,int); - double CostFunction(void*,int,int); + void Du(Vec, int,int); + void Gradj(Vec, int, int,char*); + void GradjDrag(Vec, int,int ); + void GradjB(Vec, int,int ); + double Misfit(int,int); + double SurfaceArea(int,int); + double CostFunction(int,int); void GetNodalFunctions(double* l1l2, double gauss_coord); void GetParameterValue(double* pvalue, double* value_list,double gauss_coord); Index: /issm/trunk/src/c/objects/BeamVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/BeamVertexInput.cpp (revision 3620) +++ /issm/trunk/src/c/objects/BeamVertexInput.cpp (revision 3620) @@ -0,0 +1,140 @@ +/*!\file BeamVertexInput.c + * \brief: implementation of the BeamVertexInput 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 "./BeamVertexInput.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 BeamVertexInput::BeamVertexInput(){{{1*/ +BeamVertexInput::BeamVertexInput(){ + return; +} +/*}}}*/ +/*FUNCTION BeamVertexInput::BeamVertexInput(double* values){{{1*/ +BeamVertexInput::BeamVertexInput(int in_enum_type,double* in_values){ + + enum_type=in_enum_type; + values[0]=in_values[0]; + values[1]=in_values[1]; +} +/*}}}*/ +/*FUNCTION BeamVertexInput::~BeamVertexInput(){{{1*/ +BeamVertexInput::~BeamVertexInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION BeamVertexInput::copy{{{1*/ +Object* BeamVertexInput::copy() { + + return new BeamVertexInput(this->enum_type,this->values); + +} +/*}}}*/ +/*FUNCTION BeamVertexInput::DeepEcho{{{1*/ +void BeamVertexInput::DeepEcho(void){ + + printf("BeamVertexInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %g|%g\n",this->values[0],this->values[1]); +} +/*}}}*/ +/*FUNCTION BeamVertexInput::Demarshall{{{1*/ +void BeamVertexInput::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(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION BeamVertexInput::Echo {{{1*/ +void BeamVertexInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION BeamVertexInput::Enum{{{1*/ +int BeamVertexInput::Enum(void){ + + return BeamVertexInputEnum; + +} +/*}}}*/ +/*FUNCTION BeamVertexInput::EnumType{{{1*/ +int BeamVertexInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION BeamVertexInput::GetId{{{1*/ +int BeamVertexInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION BeamVertexInput::GetName{{{1*/ +char* BeamVertexInput::GetName(void){ + return "triavertexinput"; +} +/*}}}*/ +/*FUNCTION BeamVertexInput::Marshall{{{1*/ +void BeamVertexInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of BeamVertexInput: */ + enum_type=BeamVertexInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall BeamVertexInput data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION BeamVertexInput::MarshallSize{{{1*/ +int BeamVertexInput::MarshallSize(){ + + return sizeof(values)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION BeamVertexInput::MyRank{{{1*/ +int BeamVertexInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/BeamVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/BeamVertexInput.h (revision 3620) +++ /issm/trunk/src/c/objects/BeamVertexInput.h (revision 3620) @@ -0,0 +1,42 @@ +/*! \file BeamVertexInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" + +#ifndef _BEAMVERTEXINPUT_H_ +#define _BEAMVERTEXINPUT_H_ + +class BeamVertexInput: public Input{ + + private: + /*just hold 2 values for 2 vertices: */ + int enum_type; + double values[2]; + + public: + + /*constructors, destructors: {{{1*/ + BeamVertexInput(); + BeamVertexInput(int enum_type,double* values); + ~BeamVertexInput(); + /*}}}*/ + /*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 /* _BEAMVERTEXINPUT_H */ Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 3619) +++ /issm/trunk/src/c/objects/Element.h (revision 3620) @@ -10,4 +10,6 @@ #include "./Object.h" +#include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "../toolkits/toolkits.h" @@ -18,5 +20,5 @@ virtual ~Element(){}; virtual int Enum()=0; - virtual void Configure(void* loads,void* nodes,void* materials,void* parameters)=0; + virtual void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters)=0; virtual void CreateKMatrix(Mat Kgg,int analysis_type,int sub_analysis_type)=0; Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 3619) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 3620) @@ -39,5 +39,5 @@ /*FUNCTION FemModel::FemModel {{{1*/ FemModel::FemModel(DataSet* femmodel_elements,DataSet* femmodel_nodes,DataSet* femmodel_vertices, DataSet* femmodel_constraints,DataSet* femmodel_loads, - DataSet* femmodel_materials,DataSet* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg, + DataSet* femmodel_materials,Parameters* femmodel_parameters, DofVec* femmodel_partition,DofVec* femmodel_tpartition,DofVec* femmodel_yg, Mat femmodel_Rmg,Mat femmodel_Gmn,NodeSets* femmodel_nodesets,Vec femmodel_ys,Vec femmodel_ys0){ @@ -266,5 +266,5 @@ /*}}}*/ /*FUNCTION FemModel::get_parameters {{{1*/ -DataSet* FemModel::get_parameters(void){ +Parameters* FemModel::get_parameters(void){ return parameters; } Index: /issm/trunk/src/c/objects/FemModel.h =================================================================== --- /issm/trunk/src/c/objects/FemModel.h (revision 3619) +++ /issm/trunk/src/c/objects/FemModel.h (revision 3620) @@ -11,4 +11,5 @@ #include "./Object.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Parameters.h" #include "./DofVec.h" #include "../toolkits/toolkits.h" @@ -27,5 +28,5 @@ DataSet* loads; DataSet* materials; - DataSet* parameters; + Parameters* parameters; DofVec* partition; @@ -41,5 +42,5 @@ FemModel(); ~FemModel(); - FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,DataSet* parameters, + FemModel(DataSet* elements,DataSet* nodes,DataSet* vertices, DataSet* constraints,DataSet* loads,DataSet* materials,Parameters* parameters, DofVec* partition,DofVec* tpartition,DofVec* yg,Mat Rmg,Mat Gmn,NodeSets* nodesets,Vec ys,Vec ys0); @@ -69,5 +70,5 @@ DataSet* get_loads(void); DataSet* get_materials(void); - DataSet* get_parameters(void); + Parameters* get_parameters(void); DofVec* get_partition(void); DofVec* get_tpartition(void); Index: /issm/trunk/src/c/objects/Model.cpp =================================================================== --- /issm/trunk/src/c/objects/Model.cpp (revision 3619) +++ /issm/trunk/src/c/objects/Model.cpp (revision 3620) @@ -47,5 +47,5 @@ DataSet* loads=NULL; DataSet* materials=NULL; - DataSet* parameters=NULL; + Parameters* parameters=NULL; DofVec* partition=NULL; DofVec* tpartition=NULL; @@ -121,5 +121,5 @@ DataSet* loads=NULL; DataSet* materials=NULL; - DataSet* parameters=NULL; + Parameters* parameters=NULL; DofVec* partition=NULL; DofVec* tpartition=NULL; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 3619) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 3620) @@ -10,4 +10,5 @@ #include "stdio.h" +#include "./PentaVertexInput.h" #include "./Penta.h" #include @@ -20,72 +21,119 @@ /*FUNCTION Penta default constructor {{{1*/ Penta::Penta(){ + this->inputs=NULL; return; } /*}}}*/ /*FUNCTION Penta constructor {{{1*/ -Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id, ElementProperties* pentaproperties): +Penta::Penta(int penta_id,int* penta_node_ids, int penta_matice_id, int penta_matpar_id, int penta_numpar_id): hnodes(penta_node_ids,6), hmatice(&penta_matice_id,1), hmatpar(&penta_matpar_id,1), - hnumpar(&penta_numpar_id,1), - properties(pentaproperties) + hnumpar(&penta_numpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=penta_id; - - return; + this->inputs=new Inputs(); + } /*}}}*/ /*FUNCTION Penta other constructor {{{1*/ -Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, ElementProperties* penta_properties): +Penta::Penta(int penta_id,Hook* penta_hnodes, Hook* penta_hmatice, Hook* penta_hmatpar, Hook* penta_hnumpar, Inputs* penta_inputs): hnodes(penta_hnodes), hmatice(penta_hmatice), hmatpar(penta_hmatpar), - hnumpar(penta_hnumpar), - properties(penta_properties) + hnumpar(penta_hnumpar) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=penta_id; - - return; + if(penta_inputs){ + this->inputs=(Inputs*)penta_inputs->Copy(); + } + else{ + this->inputs=new Inputs(); + } } /*}}}*/ /*FUNCTION Penta other constructor {{{1*/ -Penta::Penta(int i, IoModel* iomodel){ //i is the element index - - int j; +Penta::Penta(int index, IoModel* iomodel){ //i is the element index + + IssmInt i; int penta_node_ids[6]; int penta_matice_id; int penta_matpar_id; int penta_numpar_id; + double nodeinputs[6]; + bool collapse; /*id: */ - this->id=i+1; + this->id=index+1; /*hooks: */ - for(j=0;j<6;j++){ //go recover node ids, needed to initialize the node hook. - penta_node_ids[j]=(int)*(iomodel->elements+6*i+j); //ids for vertices are in the elements array from Matlab - } - penta_matice_id=i+1; //refers to the corresponding ice material object + for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. + penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab + } + penta_matice_id=index+1; //refers to the corresponding ice material object penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object penta_numpar_id=1; //refers to numerical parameters object - this->hnodes.Init(penta_node_ids,6); + this->hnodes.Init(&penta_node_ids[0],6); this->hmatice.Init(&penta_matice_id,1); this->hmatpar.Init(&penta_matpar_id,1); this->hnumpar.Init(&penta_numpar_id,1); - this->properties.Init(6,penta_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<6;i++)nodeinputs[i]=iomodel->thickness[penta_node_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(ThicknessEnum,nodeinputs)); + } + if (iomodel->surface) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->surface[penta_node_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(SurfaceEnum,nodeinputs)); + } + if (iomodel->bed) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->bed[penta_node_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_node_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_node_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(MeltingRateEnum,nodeinputs)); + } + if (iomodel->accumulation_rate) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->accumulation_rate[penta_node_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(AccumulationRateEnum,nodeinputs)); + } + if (iomodel->geothermalflux) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->geothermalflux[penta_node_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(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])); //elements of type 3 are macayeal type penta. we collapse the formulation on their base. if(iomodel->elements_type){ if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ - this->properties.collapse=1; + collapse=1; } else{ - this->properties.collapse=0; + collapse=0; } } + this->inputs->AddInput(new BoolInput(CollapseEnum,(IssmBool)collapse)); } @@ -93,4 +141,5 @@ /*FUNCTION Penta destructor {{{1*/ Penta::~Penta(){ + delete inputs; return; } @@ -125,5 +174,5 @@ /*FUNCTION copy {{{1*/ Object* Penta::copy() { - return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,&this->properties); + return new Penta(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); } /*}}}*/ @@ -148,6 +197,6 @@ hnumpar.Demarshall(&marshalled_dataset); - /*demarshall properties: */ - properties.Demarshall(&marshalled_dataset); + /*demarshall inputs: */ + inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); /*return: */ @@ -166,5 +215,6 @@ hmatpar.DeepEcho(); hnumpar.DeepEcho(); - properties.DeepEcho(); + printf(" inputs\n"); + inputs->DeepEcho(); return; @@ -181,7 +231,14 @@ hmatpar.Echo(); hnumpar.Echo(); - properties.Echo(); - - return; + printf(" inputs\n"); + inputs->Echo(); + +} +/*}}}*/ +/*FUNCTION Enum {{{1*/ +int Penta::Enum(void){ + + return PentaEnum; + } /*}}}*/ @@ -191,4 +248,6 @@ char* marshalled_dataset=NULL; int enum_type=0; + char* marshalled_inputs=NULL; + int marshalled_inputs_size; /*recover marshalled_dataset: */ @@ -210,6 +269,12 @@ hnumpar.Marshall(&marshalled_dataset); - /*Marshall properties: */ - properties.Marshall(&marshalled_dataset); + /*Marshall inputs: */ + marshalled_inputs_size=inputs->MarshallSize(); + marshalled_inputs=inputs->Marshall(); + memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); + marshalled_dataset+=marshalled_inputs_size; + + xfree((void**)&marshalled_inputs); + *pmarshalled_dataset=marshalled_dataset; @@ -225,5 +290,5 @@ +hmatpar.MarshallSize() +hnumpar.MarshallSize() - +properties.MarshallSize() + +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type } @@ -240,6 +305,5 @@ Hook* tria_hmatpar=NULL; Hook* tria_hnumpar=NULL; - ElementProperties* tria_properties=NULL; - DataSet* tria_inputs=NULL; + Inputs* tria_inputs=NULL; indices[0]=g0; @@ -251,8 +315,7 @@ tria_hmatpar=this->hmatpar.Spawn(&zero,1); tria_hnumpar=this->hnumpar.Spawn(&zero,1); - tria_properties=(ElementProperties*)this->properties.Spawn(indices,3); - tria_inputs=(DataSet*)this->inputs->Spawn(indices,3); - - tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_properties,tria_inputs); + tria_inputs=(Inputs*)this->inputs->Spawn(indices,3); + + tria=new Tria(this->id,tria_hnodes,tria_hmatice,tria_hmatpar,tria_hnumpar,tria_inputs); delete tria_hnodes; @@ -260,5 +323,4 @@ delete tria_hmatpar; delete tria_hnumpar; - delete tria_properties; delete tria_inputs; @@ -266,4 +328,6 @@ } /*}}}*/ + +/*updates: */ /*FUNCTION UpdateFromDakota {{{1*/ void Penta::UpdateFromDakota(void* vinputs){ @@ -275,4 +339,113 @@ /*FUNCTION Penta::UpdateInputs {{{1*/ void Penta::UpdateInputs(double* solution, int analysis_type, int sub_analysis_type){ + + /*Just branch to the correct UpdateInputs generator, according to the type of analysis we are carrying out: */ + if (analysis_type==ControlAnalysisEnum){ + + UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==DiagnosticAnalysisEnum){ + + if (sub_analysis_type==HorizAnalysisEnum){ + + UpdateInputsDiagnosticHoriz( solution,analysis_type,sub_analysis_type); + } + else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); + + } + else if (analysis_type==SlopecomputeAnalysisEnum){ + + UpdateInputsSlopeCompute( solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==PrognosticAnalysisEnum){ + + UpdateInputsPrognostic( solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==Prognostic2AnalysisEnum){ + + UpdateInputsPrognostic2(solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + + UpdateInputsBalancedthickness( solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==Balancedthickness2AnalysisEnum){ + + UpdateInputsBalancedthickness2( solution,analysis_type,sub_analysis_type); + } + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ + + UpdateInputsBalancedvelocities( solution,analysis_type,sub_analysis_type); + } + else{ + + ISSMERROR("%s%i%s\n","analysis: ",analysis_type," not supported yet"); + } +} +/*Object functions*/ +/*FUNCTION Penta::UpdateInputsDiagnosticHoriz {{{1*/ +void Penta::UpdateInputsDiagnosticHoriz(double* solution, int analysis_type, int sub_analysis_type){ + + + int i; + + const int numvertices=6; + const int numdofpervertex=2; + const int numdof=numdofpervertex*numvertices; + + int doflist[numdof]; + double values[numdof]; + double vx[numvertices]; + double vy[numvertices]; + + int dummy; + + /*Get dof list: */ + GetDofList(&doflist[0],&dummy); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iinputs->AddInput(new PentaVertexInput(VxEnum,vx)); + this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); +} + +/*}}}*/ +/*FUNCTION Penta::UpdateInputsSlopeCompute {{{1*/ +void Penta::UpdateInputsSlopeCompute(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::UpdateInputsPrognostic {{{1*/ +void Penta::UpdateInputsPrognostic(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::UpdateInputsPrognostic2 {{{1*/ +void Penta::UpdateInputsPrognostic2(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::UpdateInputsBalancedthickness {{{1*/ +void Penta::UpdateInputsBalancedthickness(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::UpdateInputsBalancedthickness2 {{{1*/ +void Penta::UpdateInputsBalancedthickness2(double* solution, int analysis_type, int sub_analysis_type){ + ISSMERROR(" not supported yet!"); +} +/*}}}*/ +/*FUNCTION Penta::UpdateInputsBalancedvelocities {{{1*/ +void Penta::UpdateInputsBalancedvelocities(double* solution, int analysis_type, int sub_analysis_type){ ISSMERROR(" not supported yet!"); } @@ -295,6 +468,4 @@ double bed; double basalforce[3]; - double vxvyvz_list[numgrids][3]; - double pressure_list[numgrids]; double epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ double devstresstensor[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/ @@ -304,5 +475,4 @@ int dofv[3]={0,1,2}; int dofp[1]={3}; - ParameterInputs* inputs=NULL; double Jdet2d; Tria* tria=NULL; @@ -323,4 +493,7 @@ double surface=0; double value=0; + + /*flags: */ + bool onbed; /*dynamic objects pointed to by hooks: */ @@ -333,7 +506,4 @@ if (analysis_type!=DiagnosticAnalysisEnum || sub_analysis_type!=StokesAnalysisEnum) ISSMERROR("Not supported yet!"); - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); @@ -342,5 +512,8 @@ numpar=(Numpar*)hnumpar.delivers(); - if(!this->properties.onbed){ + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + if(!onbed){ //put zero VecSetValue(sigma_b,id-1,0.0,INSERT_VALUES); @@ -351,8 +524,4 @@ rho_ice=matpar->GetRhoIce(); gravity=matpar->GetG(); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofv,numgrids,(void**)nodes); - inputs->Recover("velocity",&pressure_list[0] ,1,dofp,numgrids,(void**)nodes); /* Get node coordinates and dof list: */ @@ -378,7 +547,7 @@ /*Compute strain rate viscosity and pressure: */ - GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); - GetParameterValue(&pressure,&pressure_list[0],&gauss_coord[0]); + inputs->GetParameterValue(&pressure, &gauss_coord[0],PressureEnum); /*Compute Stress*/ @@ -417,8 +586,13 @@ int i; const int numgrids=6; - int doflist[numgrids]; + int doflist[numgrids]; double pressure[numgrids]; double rho_ice,g; - double xyz_list[numgrids][3]; + double surface[numgrids]; + double xyz_list[numgrids][3]; + double gauss[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + + /*inputs: */ + bool onwater; /*dynamic objects pointed to by hooks: */ @@ -430,6 +604,9 @@ matpar=(Matpar*)hmatpar.delivers(); + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Get node data: */ @@ -441,4 +618,7 @@ /*Get dof list on which we will plug the pressure values: */ GetDofList1(&doflist[0]); + + /*recover value of surface at grids: */ + inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum); /*pressure is lithostatic: */ @@ -446,5 +626,5 @@ g=matpar->GetG(); for(i=0;iproperties.s[i]-xyz_list[i][2]); + pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]); } @@ -467,19 +647,31 @@ Tria* tria=NULL; + /*flags: */ + bool onbed; + bool onwater; + bool collapse; + bool onsurface; + + /*recover some inputs: */ + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + /*If on water, return 0: */ - if(this->properties.onwater)return 0; + if(onwater)return 0; /*Bail out if this element if: * -> Non collapsed and not on the surface * -> collapsed (2d model) and not on bed) */ - if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){ + if ((!collapse && !onsurface) || (collapse && !onbed)){ return 0; } - else if (this->properties.collapse){ + else if (collapse){ /*This element should be collapsed into a tria element at its base. Create this tria element, * and compute CostFunction*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->CostFunction(inputs,analysis_type,sub_analysis_type); + J=tria->CostFunction(analysis_type,sub_analysis_type); delete tria; return J; @@ -488,5 +680,5 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->CostFunction(inputs,analysis_type,sub_analysis_type); + J=tria->CostFunction(analysis_type,sub_analysis_type); delete tria; return J; @@ -501,5 +693,5 @@ if (analysis_type==ControlAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); } @@ -508,13 +700,13 @@ if (sub_analysis_type==HorizAnalysisEnum){ - CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHoriz( Kgg,analysis_type,sub_analysis_type); } else if (sub_analysis_type==VertAnalysisEnum){ - CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticVert( Kgg,analysis_type,sub_analysis_type); } else if (sub_analysis_type==StokesAnalysisEnum){ - CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticStokes( Kgg,analysis_type,sub_analysis_type); } @@ -523,9 +715,9 @@ 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==BalancedthicknessAnalysisEnum){ @@ -539,9 +731,9 @@ else if (analysis_type==ThermalAnalysisEnum){ - CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixThermal( Kgg,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum){ - CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixMelting( Kgg,analysis_type,sub_analysis_type); } else{ @@ -645,11 +837,12 @@ double Bprime[5][numdof]; double L[2][numdof]; - double D[5][5]={{ 0,0,0,0,0 },{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0},{0,0,0,0,0}}; // material matrix, simple scalar matrix. + double D[5][5]={0.0}; // material matrix, simple scalar matrix. double D_scalar; - double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag + double DL[2][2]={0.0}; //for basal drag double DL_scalar; /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix + double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag @@ -657,12 +850,6 @@ /*slope: */ - double slope[2]={0.0,0.0}; + double slope[2]={0.0}; double slope_magnitude; - - /*input parameters for structural analysis (diagnostic): */ - double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; - int dofs[2]={0,1}; - double thickness; /*friction: */ @@ -672,6 +859,4 @@ double MAXSLOPE=.06; // 6 % double MOUNTAINKEXPONENT=10; - - ParameterInputs* inputs=NULL; /*Collapsed formulation: */ @@ -688,9 +873,18 @@ numpar=(Numpar*)hnumpar.delivers(); + /*inputs: */ + bool onwater; + bool collapse; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*If on water, skip stiffness: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return; /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the @@ -699,15 +893,15 @@ - if ((this->properties.collapse==1) && (this->properties.onbed==0)){ + if ((collapse==1) && (onbed==0)){ /*This element should be collapsed, but this element is not on the bedrock, therefore all its * dofs have already been frozen! Do nothing: */ return; } - else if ((this->properties.collapse==1) && (this->properties.onbed==1)){ + else if ((collapse==1) && (onbed==1)){ /*This element should be collapsed into a tria element at its base. Create this tria element, *and use its CreateKMatrix functionality to fill the global stiffness matrix: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); delete tria; return; @@ -716,8 +910,4 @@ /*Implement standard penta element: */ - - /*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: */ @@ -725,10 +915,4 @@ GetDofList(&doflist[0],&numberofdofspernode); - /* Set Ke_gg to 0: */ - for(i=0;iGetStrainRate(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum); + inputs->GetStrainRate(&oldepsilon[0],&xyz_list[0][0],gauss_coord,VxOldEnum,VyOldEnum); /*Get viscosity: */ @@ -774,5 +958,5 @@ GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord); - /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant + /*Build the D matrix: we plug the gaussian weight, the viscosity, and the jacobian determinant onto this scalar matrix, so that we win some computational time: */ @@ -803,5 +987,5 @@ //Deal with 2d friction at the bedrock interface - if((this->properties.onbed && !this->properties.shelf)){ + if((onbed && !shelf)){ /*Build a tria element using the 3 grids of the base of the penta. Then use @@ -810,5 +994,5 @@ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type,sub_analysis_type); + tria->CreateKMatrixDiagnosticHorizFriction(Kgg,analysis_type,sub_analysis_type); delete tria; } @@ -845,5 +1029,5 @@ int dofs[3]={0,1,2}; - double K_terms[numdof][numdof]; + double K_terms[numdof][numdof]={0.0}; /*Material properties: */ @@ -860,5 +1044,4 @@ double surface_normal[3]; double bed_normal[3]; - double vxvyvz_list[numgrids][3]; double thickness; @@ -889,4 +1072,5 @@ double* area_gauss_weights = NULL; double* vert_gauss_weights = NULL; + double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; /* specific gaussian point: */ @@ -902,5 +1086,11 @@ double alpha2_gauss; - ParameterInputs* inputs=NULL; + double vx_list[numgrids]; + double vy_list[numgrids]; + double vz_list[numgrids]; + double thickness_list[numgrids]; + double bed_list[numgrids]; + double dragcoefficient_list[numgrids]; + double drag_p,drag_q; /*dynamic objects pointed to by hooks: */ @@ -910,6 +1100,16 @@ Numpar* numpar=NULL; + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*If on water, skip stiffness: */ - if(this->properties.onwater)return; + if(onwater)return; /*recover objects from hooks: */ @@ -919,13 +1119,4 @@ numpar=(Numpar*)hnumpar.delivers(); - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - - /* Set K_terms to 0: */ - for(i=0;iRecover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore @@ -965,8 +1153,6 @@ gauss_coord[3]=*(vert_gauss_coord+igvert); - /*Compute thickness: */ - /*Compute strain rate: */ - GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); /*Get viscosity: */ @@ -1003,5 +1189,5 @@ } - if((this->properties.onbed==1) && (this->properties.shelf==0)){ + if((onbed==1) && (shelf==0)){ /*Build alpha2_list used by drag stiffness matrix*/ @@ -1012,14 +1198,24 @@ strcpy(friction->element_type,"2d"); + inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],6,VxEnum); + inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],6,VyEnum); + inputs->GetParameterValues(&vz_list[0],&gaussgrids[0][0],6,VzEnum); + inputs->GetParameterValues(&dragcoefficient_list[0],&gaussgrids[0][0],6,DragCoefficientEnum); + inputs->GetParameterValues(&bed_list[0],&gaussgrids[0][0],6,BedEnum); + inputs->GetParameterValues(&thickness_list[0],&gaussgrids[0][0],6,ThicknessEnum); + inputs->GetParameterValue(&drag_p,DragPEnum); + inputs->GetParameterValue(&drag_q,DragQEnum); + 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->analysis_type=analysis_type; + 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->vz=&vz_list[0]; + friction->p=drag_p; + friction->q=drag_q; /*Compute alpha2_list: */ @@ -1058,5 +1254,5 @@ /*Compute strain rate: */ - GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); /*Get viscosity at last iteration: */ @@ -1098,5 +1294,5 @@ } } - } //if ( (this->properties.onbed==1) && (shelf==0)) + } //if ( (onbed==1) && (shelf==0)) /*Reduce the matrix */ @@ -1155,5 +1351,5 @@ /* matrices: */ - double Ke_gg[numdof][numdof]; + double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; double Jdet; @@ -1168,21 +1364,23 @@ nodes=(Node**)hnodes.deliverp(); - ParameterInputs* inputs=NULL; - /*Collapsed formulation: */ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onsurface; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + /*If on water, skip stiffness: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - + if(onwater)return; /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness * matrix: */ - if(this->properties.onsurface){ + if(onsurface){ tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface - tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrixDiagnosticSurfaceVert(Kgg, analysis_type,sub_analysis_type); delete tria; } @@ -1194,10 +1392,4 @@ GetDofList(&doflist[0],&numberofdofspernode); - /* Set Ke_gg to 0: */ - for(i=0;iGetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; - - if (!this->properties.onbed){ + if(onwater)return; + + if (!onbed){ return; } @@ -1274,5 +1474,5 @@ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrixMelting(Kgg, analysis_type,sub_analysis_type); delete tria; return; @@ -1287,13 +1487,21 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); delete tria; return; @@ -1308,13 +1516,22 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrix(Kgg, analysis_type,sub_analysis_type); delete tria; return; @@ -1356,14 +1573,5 @@ double K[2][2]={0.0}; - double vxvyvz_list[numgrids][3]; - double vx_list[numgrids]; - int vx_list_exists; - double u; - double vy_list[numgrids]; - int vy_list_exists; - double v; - double vz_list[numgrids]; - int vz_list_exists; - double w; + double u,v,w; /*matrices: */ @@ -1404,8 +1612,18 @@ /*Collapsed formulation: */ Tria* tria=NULL; - ParameterInputs* inputs=NULL; + + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*recover objects from hooks: */ @@ -1413,7 +1631,4 @@ matpar=(Matpar*)hmatpar.delivers(); numpar=(Numpar*)hnumpar.delivers(); - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; /* Get node coordinates and dof list: */ @@ -1427,17 +1642,5 @@ heatcapacity=matpar->GetHeatCapacity(); thermalconductivity=matpar->GetThermalConductivity(); - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); - - for(i=0;idt; /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore @@ -1495,7 +1698,7 @@ //Build the D matrix - GetParameterValue(&u, &vx_list[0],gauss_coord); - GetParameterValue(&v, &vy_list[0],gauss_coord); - GetParameterValue(&w, &vz_list[0],gauss_coord); + inputs->GetParameterValue(&u, gauss_coord,VxEnum); + inputs->GetParameterValue(&v, gauss_coord,VyEnum); + inputs->GetParameterValue(&w, gauss_coord,VzEnum); D_scalar=gauss_weight*Jdet; @@ -1577,8 +1780,8 @@ //Ice/ocean heat exchange flux on ice shelf base - if(this->properties.onbed && this->properties.shelf){ + if(onbed && shelf){ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrixThermal(Kgg,inputs, analysis_type,sub_analysis_type); + tria->CreateKMatrixThermal(Kgg, analysis_type,sub_analysis_type); delete tria; } @@ -1591,5 +1794,5 @@ if (analysis_type==ControlAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); } else if (analysis_type==DiagnosticAnalysisEnum){ @@ -1597,13 +1800,13 @@ if (sub_analysis_type==HorizAnalysisEnum){ - CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHoriz( pg,analysis_type,sub_analysis_type); } else if (sub_analysis_type==VertAnalysisEnum){ - CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticVert( pg,analysis_type,sub_analysis_type); } else if (sub_analysis_type==StokesAnalysisEnum){ - CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticStokes( pg,analysis_type,sub_analysis_type); } else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"); @@ -1611,25 +1814,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,analysis_type,sub_analysis_type); + } + else if (analysis_type==BalancedthicknessAnalysisEnum){ + CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); } - else if (analysis_type==BalancedthicknessAnalysisEnum){ + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - - CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type); - } else if (analysis_type==ThermalAnalysisEnum){ - CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorThermal( pg,analysis_type,sub_analysis_type); } else if (analysis_type==MeltingAnalysisEnum){ - CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorMelting( pg,analysis_type,sub_analysis_type); } else{ @@ -1720,8 +1923,6 @@ /*element vector at the gaussian points: */ - double pe_g[numdof]; + double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; - - ParameterInputs* inputs=NULL; /*dynamic objects pointed to by hooks: */ @@ -1732,9 +1933,16 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool collapse; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip load: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return; /*recover objects from hooks: */ @@ -1747,15 +1955,15 @@ the load vector. */ - if ((this->properties.collapse==1) && (this->properties.onbed==0)){ + if ((collapse==1) && (onbed==0)){ /*This element should be collapsed, but this element is not on the bedrock, therefore all its * dofs have already been frozen! Do nothing: */ return; } - else if ((this->properties.collapse==1) && (this->properties.onbed==1)){ + else if ((collapse==1) && (onbed==1)){ /*This element should be collapsed into a tria element at its base. Create this tria element, *and use its CreatePVector functionality to return an elementary load vector: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVector(pg, analysis_type,sub_analysis_type); delete tria; return; @@ -1768,7 +1976,4 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set pe_g to 0: */ - for(i=0;iproperties.h[0],gauss_coord); + inputs->GetParameterValue(&thickness, gauss_coord,ThicknessEnum); /*Compute slope at gaussian point: */ - GetParameterDerivativeValue(&slope[0], &this->properties.s[0],&xyz_list[0][0], gauss_coord); + inputs->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss_coord,SurfaceEnum); /* Get Jacobian determinant: */ @@ -1820,5 +2025,5 @@ } //for (ig1=0; ig1properties.onbed==1)) + } //else if ((collapse==1) && (onbed==1)) /*Add pe_g to global vector pg: */ @@ -1857,5 +2062,4 @@ double bed_normal[3]; double bed; - double vxvyvz_list[numgrids][3]; /* gaussian points: */ @@ -1898,7 +2102,6 @@ double D_scalar; double tBD[27][8]; - double P_terms[numdof]; - - ParameterInputs* inputs=NULL; + double P_terms[numdof]={0.0}; + Tria* tria=NULL; @@ -1909,9 +2112,16 @@ Numpar* numpar=NULL; + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*If on water, skip load: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return; /*recover objects from hooks: */ @@ -1921,8 +2131,4 @@ numpar=(Numpar*)hnumpar.delivers(); - /* Set P_terms to 0: */ - for(i=0;iGetRhoIce(); gravity=matpar->GetG(); - - /*recover extra inputs from users, at current convergence iteration: */ - inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); /* Get node coordinates and dof list: */ @@ -1960,5 +2163,5 @@ /*Compute strain rate and viscosity: */ - GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -2014,5 +2217,5 @@ /*Deal with 2d friction at the bedrock interface: */ - if ( (this->properties.onbed==1) && (this->properties.shelf==1)){ + if ( (onbed==1) && (shelf==1)){ for(i=0;iproperties.b[0],gauss_coord); + inputs->GetParameterValue(&bed, gauss_coord,BedEnum); /*Get L matrix : */ @@ -2062,5 +2265,5 @@ } } - } //if ( (this->properties.onbed==1) && (this->properties.shelf==1)) + } //if ( (onbed==1) && (shelf==1)) /*Reduce the matrix */ @@ -2118,5 +2321,5 @@ /*element vector at the gaussian points: */ - double pe_g[numdof]; + double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; double l1l6[6]; @@ -2126,6 +2329,4 @@ /*input parameters for structural analysis (diagnostic): */ - double vx_list[numgrids]={0,0,0,0,0,0}; - double vy_list[numgrids]={0,0,0,0,0,0}; double du[3]; double dv[3]; @@ -2137,20 +2338,24 @@ Node** nodes=NULL; - ParameterInputs* inputs=NULL; + /*inputs: */ + bool onwater; + bool onbed; /*recover objects from hooks: */ nodes=(Node**)hnodes.deliverp(); - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the *diagnostic base vertical stifness: */ - if(this->properties.onbed){ + if(onbed){ tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock - tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVectorDiagnosticBaseVert(pg, analysis_type,sub_analysis_type); delete tria; } @@ -2161,11 +2366,4 @@ GetDofList(&doflist[0],&numberofdofspernode); - /* Set pe_g to 0: */ - for(i=0;iRecover("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 gaussian points and weights :*/ order_area_gauss=2; @@ -2189,6 +2387,7 @@ /*Get velocity derivative, with respect to x and y: */ - GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_coord); - GetParameterDerivativeValue(&dv[0], &vy_list[0],&xyz_list[0][0], gauss_coord); + + inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VxEnum); + inputs->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss_coord,VyEnum); dudx=du[0]; dvdy=dv[1]; @@ -2235,13 +2434,21 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVector(pg, analysis_type,sub_analysis_type); delete tria; return; @@ -2255,13 +2462,21 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Is this element on the bed? :*/ - if(!this->properties.onbed)return; + if(!onbed)return; /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVector(pg, analysis_type,sub_analysis_type); delete tria; return; @@ -2299,8 +2514,4 @@ double dt; - double vx_list[numgrids]; - double vy_list[numgrids]; - double vz_list[numgrids]; - double vxvyvz_list[numgrids][3]; double temperature_list[numgrids]; double temperature; @@ -2340,5 +2551,4 @@ /*Collapsed formulation: */ Tria* tria=NULL; - ParameterInputs* inputs=NULL; /*dynamic objects pointed to by hooks: */ @@ -2346,10 +2556,18 @@ Matpar* matpar=NULL; Matice* matice=NULL; + Numpar* numpar=NULL; + + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); /*If on water, skip: */ - if(this->properties.onwater)return; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; + if(onwater)return; /*recover objects from hooks: */ @@ -2357,4 +2575,5 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); + numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -2369,21 +2588,5 @@ beta=matpar->GetBeta(); meltingpoint=matpar->GetMeltingPoint(); - - /*recover extra inputs from users, dt and velocity: */ - found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find velocity in inputs!"); - found=inputs->Recover("dt",&dt); - if(!found)ISSMERROR(" could not find dt in inputs!"); - - if(dt){ - found=inputs->Recover("temperature",&temperature_list[0],1,dofs1,numgrids,(void**)nodes); - if(!found)ISSMERROR(" could not find temperature in inputs!"); - } - - for(i=0;idt; /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore @@ -2409,5 +2612,5 @@ /*Compute strain rate and viscosity: */ - GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord); + inputs->GetStrainRateStokes(&epsilon[0],&xyz_list[0][0],gauss_coord,VxEnum,VyEnum,VzEnum); matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); @@ -2433,5 +2636,5 @@ /* Build transient now */ if(dt){ - GetParameterValue(&temperature, &temperature_list[0],gauss_coord); + inputs->GetParameterValue(&temperature, gauss_coord,TemperatureEnum); scalar_transient=temperature*Jdet*gauss_weight; for(i=0;iproperties.onbed && this->properties.shelf){ + if(onbed && shelf){ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalShelf(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVectorThermalShelf(pg, analysis_type,sub_analysis_type); delete tria; } /* Geothermal flux on ice sheet base and basal friction */ - if(this->properties.onbed && !this->properties.shelf){ + if(onbed && !shelf){ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVectorThermalSheet(pg,inputs, analysis_type,sub_analysis_type); + tria->CreatePVectorThermalSheet(pg, analysis_type,sub_analysis_type); delete tria; } @@ -2477,19 +2680,31 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool collapse; + bool onsurface; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*Bail out if this element if: * -> Non collapsed and not on the surface * -> collapsed (2d model) and not on bed) */ - if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){ + if ((!collapse && !onsurface) || (collapse && !onbed)){ return; } - else if (this->properties.collapse){ + else if (collapse){ /*This element should be collapsed into a tria element at its base. Create this tria element, * and compute Du*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->Du(du_g,inputs,analysis_type,sub_analysis_type); + tria->Du(du_g,analysis_type,sub_analysis_type); delete tria; return; @@ -2498,15 +2713,8 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - tria->Du(du_g,inputs,analysis_type,sub_analysis_type); + tria->Du(du_g,analysis_type,sub_analysis_type); delete tria; return; } -} -/*}}}*/ -/*FUNCTION Enum {{{1*/ -int Penta::Enum(void){ - - return PentaEnum; - } /*}}}*/ @@ -2527,8 +2735,16 @@ nodes=(Node**)hnodes.deliverp(); + /*inputs: */ + bool collapse; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*Figure out if we should extrude for this element: */ if (iscollapsed){ /*From higher level, we are told to extrude only elements that have the collapse flag on: */ - if (this->properties.collapse)extrude=1; + if (collapse)extrude=1; else extrude=0; } @@ -2540,5 +2756,5 @@ /*Now, extrusion starts from the bed on, so double check this element is on * the bedrock: */ - if(this->properties.onbed==0)extrude=0; + if(onbed==0)extrude=0; /*Go on and extrude field: */ @@ -2822,8 +3038,10 @@ /*}}}*/ /*FUNCTION GetBedList {{{1*/ -void Penta::GetBedList(double* bed_list){ - - int i; - for(i=0;i<6;i++)bed_list[i]=this->properties.b[i]; +void Penta::GetBedList(double* bedlist){ + + const int numgrids=6; + double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + + inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum); } @@ -3682,5 +3900,10 @@ /*FUNCTION GetOnBed {{{1*/ bool Penta::GetOnBed(){ - return this->properties.onbed; + + bool onbed; + + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + + return onbed; } /*}}}*/ @@ -3757,6 +3980,11 @@ /*}}}*/ /*FUNCTION GetShelf {{{1*/ -int Penta::GetShelf(){ - return this->properties.shelf; +bool Penta::GetShelf(){ + bool onshelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onshelf,ElementOnIceShelfEnum); + + return onshelf; } /*}}}*/ @@ -3821,6 +4049,8 @@ void Penta::GetThicknessList(double* thickness_list){ - int i; - for(i=0;i<6;i++)thickness_list[i]=this->properties.h[i]; + const int numgrids=6; + double gaussgrids[numgrids][numgrids]={{1,0,0,0},{0,1,0,0},{0,0,1,0},{1,0,0,1},{0,1,0,1},{0,0,1,1}}; + inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum); + } /*}}}*/ @@ -3828,12 +4058,18 @@ void Penta::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, skip 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); @@ -3845,12 +4081,22 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool onbed; + bool shelf; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; + if(onwater)return; /*If on shelf, skip: */ - if(this->properties.shelf)return; + if(shelf)return; /*Bail out if this element does not touch the bedrock: */ - if (!this->properties.onbed) return; + if (!onbed) return; if (sub_analysis_type==HorizAnalysisEnum){ @@ -3858,5 +4104,5 @@ /*MacAyeal or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDrag( grad_g,inputs,analysis_type,sub_analysis_type); + tria->GradjDrag( grad_g,analysis_type,sub_analysis_type); delete tria; return; @@ -3866,5 +4112,5 @@ /*Stokes*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->GradjDragStokes( grad_g,inputs,analysis_type,sub_analysis_type); + tria->GradjDragStokes( grad_g,analysis_type,sub_analysis_type); delete tria; return; @@ -3878,15 +4124,25 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool collapse; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + /*If on water, skip: */ - if(this->properties.onwater)return; - - if (this->properties.collapse){ + if(onwater)return; + + if (collapse){ /*Bail out element if collapsed (2d) and not on bed*/ - if (!this->properties.onbed) return; + if (!onbed) return; /*This element should be collapsed into a tria element at its base. Create this tria element, * and compute gardj*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); + tria->GradjB(grad_g,analysis_type,sub_analysis_type); delete tria; return; @@ -3895,5 +4151,5 @@ /*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - tria->GradjB(grad_g,inputs,analysis_type,sub_analysis_type); + tria->GradjB(grad_g,analysis_type,sub_analysis_type); delete tria; return; @@ -3912,19 +4168,31 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool collapse; + bool onsurface; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + /*If on water, return 0: */ - if(this->properties.onwater)return 0; + if(onwater)return 0; /*Bail out if this element if: * -> Non collapsed and not on the surface * -> collapsed (2d model) and not on bed) */ - if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){ + if ((!collapse && !onsurface) || (collapse && !onbed)){ return 0; } - else if (this->properties.collapse){ + else if (collapse){ /*This element should be collapsed into a tria element at its base. Create this tria element, * and compute Misfit*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type); + J=tria->Misfit(analysis_type,sub_analysis_type); delete tria; return J; @@ -3933,5 +4201,5 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - J=tria->Misfit(inputs,analysis_type,sub_analysis_type); + J=tria->Misfit(analysis_type,sub_analysis_type); delete tria; return J; @@ -4045,19 +4313,31 @@ Tria* tria=NULL; + /*inputs: */ + bool onwater; + bool collapse; + bool onsurface; + bool onbed; + + /*retrieve inputs :*/ + inputs->GetParameterValue(&onwater,ElementOnWaterEnum); + inputs->GetParameterValue(&collapse,CollapseEnum); + inputs->GetParameterValue(&onbed,ElementOnBedEnum); + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + /*If on water, return 0: */ - if(this->properties.onwater)return 0; + if(onwater)return 0; /*Bail out if this element if: * -> Non collapsed and not on the surface * -> collapsed (2d model) and not on bed) */ - if ((!this->properties.collapse && !this->properties.onsurface) || (this->properties.collapse && !this->properties.onbed)){ + if ((!collapse && !onsurface) || (collapse && !onbed)){ return 0; } - else if (this->properties.collapse){ + else if (collapse){ /*This element should be collapsed into a tria element at its base. Create this tria element, * and compute SurfaceArea*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face). - S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type); + S=tria->SurfaceArea(analysis_type,sub_analysis_type); delete tria; return S; @@ -4066,5 +4346,5 @@ tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face). - S=tria->SurfaceArea(inputs,analysis_type,sub_analysis_type); + S=tria->SurfaceArea(analysis_type,sub_analysis_type); delete tria; return S; Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 3619) +++ /issm/trunk/src/c/objects/Penta.h (revision 3620) @@ -16,5 +16,4 @@ #include "./Element.h" #include "./Matpar.h" -#include "./Numpar.h" #include "./Matice.h" #include "./Tria.h" @@ -45,5 +44,5 @@ /*}}}*/ /*FUNCTION object management {{{1*/ - void Configure(void* loads,void* nodes,void* materials,void* parameters); + void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); Object* copy(); void DeepEcho(); @@ -57,6 +56,4 @@ int MyRank(); void* SpawnTria(int g0, int g1, int g2); - void UpdateFromDakota(void* inputs); - void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); void SetClone(int* minranks); @@ -107,8 +104,4 @@ void CreateKMatrixPrognostic(Mat Kgg,int analysis_type,int sub_analysis_type); void CreatePVectorPrognostic( Vec pg, 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 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 CreateKMatrixDiagnosticStokes( Mat Kgg, int analysis_type,int sub_analysis_type); @@ -136,4 +129,16 @@ void GetPhi(double* phi, double* epsilon, double viscosity); double MassFlux(double* segment,double* ug); + + /*updates: */ + void UpdateFromDakota(void* inputs); + void UpdateInputs(double* solution, int analysis_type, int sub_analysis_type); + void UpdateInputsDiagnosticHoriz( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsSlopeCompute( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsPrognostic( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsPrognostic2(double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsBalancedthickness( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsBalancedthickness2( double* solution,int analysis_type,int sub_analysis_type); + void UpdateInputsBalancedvelocities( double* solution,int analysis_type,int sub_analysis_type); + /*}}}*/ Index: /issm/trunk/src/c/objects/PentaVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/PentaVertexInput.cpp (revision 3620) +++ /issm/trunk/src/c/objects/PentaVertexInput.cpp (revision 3620) @@ -0,0 +1,144 @@ +/*!\file PentaVertexInput.c + * \brief: implementation of the PentaVertexInput 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 "./PentaVertexInput.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 PentaVertexInput::PentaVertexInput(){{{1*/ +PentaVertexInput::PentaVertexInput(){ + return; +} +/*}}}*/ +/*FUNCTION PentaVertexInput::PentaVertexInput(double* values){{{1*/ +PentaVertexInput::PentaVertexInput(int in_enum_type,double* in_values){ + + enum_type=in_enum_type; + values[0]=in_values[0]; + values[1]=in_values[1]; + values[2]=in_values[2]; + values[3]=in_values[3]; + values[4]=in_values[4]; + values[5]=in_values[5]; +} +/*}}}*/ +/*FUNCTION PentaVertexInput::~PentaVertexInput(){{{1*/ +PentaVertexInput::~PentaVertexInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION PentaVertexInput::copy{{{1*/ +Object* PentaVertexInput::copy() { + + return new PentaVertexInput(this->enum_type,this->values); + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::DeepEcho{{{1*/ +void PentaVertexInput::DeepEcho(void){ + + printf("PentaVertexInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %g|%g|%g|%g|%g|%g\n",this->values[0],this->values[1],this->values[2],this->values[3],this->values[4],this->values[5]); +} +/*}}}*/ +/*FUNCTION PentaVertexInput::Demarshall{{{1*/ +void PentaVertexInput::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(&values,marshalled_dataset,sizeof(values));marshalled_dataset+=sizeof(values); + + /*return: */ + *pmarshalled_dataset=marshalled_dataset; + return; +} +/*}}}*/ +/*FUNCTION PentaVertexInput::Echo {{{1*/ +void PentaVertexInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION PentaVertexInput::Enum{{{1*/ +int PentaVertexInput::Enum(void){ + + return PentaVertexInputEnum; + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::EnumType{{{1*/ +int PentaVertexInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION PentaVertexInput::GetId{{{1*/ +int PentaVertexInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION PentaVertexInput::GetName{{{1*/ +char* PentaVertexInput::GetName(void){ + return "pentavertexinput"; +} +/*}}}*/ +/*FUNCTION PentaVertexInput::Marshall{{{1*/ +void PentaVertexInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of PentaVertexInput: */ + enum_type=PentaVertexInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall PentaVertexInput data: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + memcpy(marshalled_dataset,&values,sizeof(values));marshalled_dataset+=sizeof(values); + + *pmarshalled_dataset=marshalled_dataset; +} +/*}}}*/ +/*FUNCTION PentaVertexInput::MarshallSize{{{1*/ +int PentaVertexInput::MarshallSize(){ + + return sizeof(values)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION PentaVertexInput::MyRank{{{1*/ +int PentaVertexInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/PentaVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/PentaVertexInput.h (revision 3620) +++ /issm/trunk/src/c/objects/PentaVertexInput.h (revision 3620) @@ -0,0 +1,42 @@ +/*! \file PentaVertexInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" + +#ifndef _PENTAVERTEXINPUT_H_ +#define _PENTAVERTEXINPUT_H_ + +class PentaVertexInput: public Input{ + + private: + /*just hold 3 values for 3 vertices: */ + int enum_type; + double values[6]; + + public: + + /*constructors, destructors: {{{1*/ + PentaVertexInput(); + PentaVertexInput(int enum_type,double* values); + ~PentaVertexInput(); + /*}}}*/ + /*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 /* _PENTAVERTEXINPUT_H */ Index: /issm/trunk/src/c/objects/Sing.cpp =================================================================== --- /issm/trunk/src/c/objects/Sing.cpp (revision 3619) +++ /issm/trunk/src/c/objects/Sing.cpp (revision 3620) @@ -11,8 +11,10 @@ #include "stdio.h" #include "./Sing.h" +#include "./SingVertexInput.h" #include #include "../EnumDefinitions/EnumDefinitions.h" #include "../shared/shared.h" #include "../DataSet/DataSet.h" +#include "../DataSet/Inputs.h" #include "../include/typedefs.h" #include "../include/macros.h" @@ -21,18 +23,19 @@ /*FUNCTION Sing::constructor {{{1*/ Sing::Sing(){ + this->inputs=NULL; return; } /*}}}*/ /*FUNCTION Sing constructor {{{1*/ -Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id, ElementProperties* singproperties): +Sing::Sing(int sing_id,int* sing_node_ids, int sing_matice_id, int sing_matpar_id, int sing_numpar_id): hnodes(sing_node_ids,1), hmatice(&sing_matice_id,1), hmatpar(&sing_matpar_id,1), - hnumpar(&sing_numpar_id,1), - properties(singproperties) + hnumpar(&sing_numpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=sing_id; + this->inputs=new Inputs(); return; @@ -40,15 +43,19 @@ /*}}}*/ /*FUNCTION Sing other constructor {{{1*/ -Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, ElementProperties* sing_properties): +Sing::Sing(int sing_id,Hook* sing_hnodes, Hook* sing_hmatice, Hook* sing_hmatpar, Hook* sing_hnumpar, Inputs* sing_inputs): hnodes(sing_hnodes), hmatice(sing_hmatice), hmatpar(sing_hmatpar), - hnumpar(sing_hnumpar), - properties(sing_properties) + hnumpar(sing_hnumpar) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=sing_id; - + if(sing_inputs){ + this->inputs=(Inputs*)sing_inputs->Copy(); + } + else{ + this->inputs=new Inputs(); + } return; } @@ -80,9 +87,9 @@ this->hnumpar.Init(&sing_numpar_id,1); - /*properties: */ - sing_h=iomodel->thickness[i]; - sing_k=iomodel->drag[i]; - - this->properties.Init(1,&sing_h, NULL, NULL, &sing_k, NULL, NULL, NULL, UNDEF, UNDEF, UNDEF, UNDEF, UNDEF,UNDEF, UNDEF,UNDEF,UNDEF); + //intialize inputs, and add as many inputs per element as requested: + this->inputs=new Inputs(); + + if (iomodel->thickness) this->inputs->AddInput(new SingVertexInput(ThicknessEnum,iomodel->thickness[i])); + if (iomodel->drag_coefficient) this->inputs->AddInput(new SingVertexInput(DragCoefficientEnum,iomodel->drag_coefficient[i])); } @@ -90,4 +97,5 @@ /*FUNCTION Sing::destructor {{{1*/ Sing::~Sing(){ + delete inputs; return; } @@ -136,5 +144,6 @@ hmatpar.DeepEcho(); hnumpar.DeepEcho(); - properties.DeepEcho(); + printf(" inputs\n"); + inputs->DeepEcho(); return; @@ -160,6 +169,6 @@ hnumpar.Demarshall(&marshalled_dataset); - /*demarshall properties: */ - properties.Demarshall(&marshalled_dataset); + /*demarshall inputs: */ + inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); /*return: */ @@ -178,7 +187,6 @@ hmatpar.Echo(); hnumpar.Echo(); - properties.Echo(); - - return; + printf(" inputs\n"); + inputs->Echo(); } /*}}}*/ @@ -188,4 +196,6 @@ char* marshalled_dataset=NULL; int enum_type=0; + char* marshalled_inputs=NULL; + int marshalled_inputs_size; /*recover marshalled_dataset: */ @@ -207,6 +217,11 @@ hnumpar.Marshall(&marshalled_dataset); - /*Marshall properties: */ - properties.Marshall(&marshalled_dataset); + /*Marshall inputs: */ + marshalled_inputs_size=inputs->MarshallSize(); + marshalled_inputs=inputs->Marshall(); + memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); + marshalled_dataset+=marshalled_inputs_size; + + xfree((void**)&marshalled_inputs); *pmarshalled_dataset=marshalled_dataset; @@ -222,5 +237,5 @@ +hmatpar.MarshallSize() +hnumpar.MarshallSize() - +properties.MarshallSize() + +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type } @@ -243,32 +258,28 @@ 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]; + int dof; + double pressure; + double thickness; 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]); + GetDofList1(&dof); /*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(); g=matpar->GetG(); - pressure[0]=rho_ice*g*this->properties.h[0]; + inputs->GetParameterValue(&thickness,ThicknessEnum); + pressure=rho_ice*g*thickness; /*plug local pressure values into global pressure vector: */ - VecSetValues(p_g,numgrids,doflist,(const double*)pressure,INSERT_VALUES); + VecSetValue(p_g,dof,pressure,INSERT_VALUES); } @@ -282,5 +293,5 @@ /*}}}*/ /*FUNCTION Sing::CostFunction {{{1*/ -double Sing::CostFunction(void*, int,int){ +double Sing::CostFunction( int,int){ ISSMERROR(" not supported yet!"); } @@ -293,5 +304,5 @@ if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){ - CreateKMatrixDiagnosticHutter( Kgg,inputs,analysis_type,sub_analysis_type); + CreateKMatrixDiagnosticHutter( Kgg,analysis_type,sub_analysis_type); } @@ -313,9 +324,4 @@ int numberofdofspernode; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - GetDofList(&doflist[0],&numberofdofspernode); @@ -330,5 +336,5 @@ if ((analysis_type==DiagnosticAnalysisEnum) && (sub_analysis_type==HutterAnalysisEnum)){ - CreatePVectorDiagnosticHutter( pg,inputs,analysis_type,sub_analysis_type); + CreatePVectorDiagnosticHutter( pg,analysis_type,sub_analysis_type); } @@ -365,9 +371,4 @@ Numpar* numpar=NULL; - ParameterInputs* inputs=NULL; - - /*recover pointers: */ - inputs=(ParameterInputs*)vinputs; - /*recover objects from hooks: */ node=(Node**)hnodes.deliverp(); @@ -376,8 +377,6 @@ numpar=(Numpar*)hnumpar.delivers(); - found=inputs->Recover("surfaceslopex",&slope[0],1,dofs,numgrids,(void**)&node[0]); - if(!found)ISSMERROR(" surfaceslopex missing from inputs!"); - found=inputs->Recover("surfaceslopey",&slope[1],1,dofs,numgrids,(void**)&node[0]); - if(!found)ISSMERROR(" surfaceslopey missing from inputs!"); + inputs->GetParameterValue(&slope[0],SurfaceSlopexEnum); + inputs->GetParameterValue(&slope[1],SurfaceSlopeyEnum); GetDofList(&doflist[0],&numberofdofspernode); @@ -391,5 +390,5 @@ n=matice->GetN(); B=matice->GetB(); - thickness=this->properties.h[0]; + inputs->GetParameterValue(&thickness,ThicknessEnum); ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0]; @@ -408,5 +407,5 @@ /*}}}*/ /*FUNCTION Sing::Du {{{1*/ -void Sing::Du(_p_Vec*,void*,int,int){ +void Sing::Du(Vec,int,int){ ISSMERROR(" not supported yet!"); } @@ -504,5 +503,5 @@ /*}}}*/ /*FUNCTION Sing::GetShelf {{{1*/ -int Sing::GetShelf(){ +bool Sing::GetShelf(){ ISSMERROR(" not supported yet!"); } @@ -514,15 +513,15 @@ /*}}}*/ /*FUNCTION Sing::Gradj {{{1*/ -void Sing::Gradj(_p_Vec*, void*, int, int ,char*){ +void Sing::Gradj(Vec, int, int ,char*){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradB {{{1*/ -void Sing::GradjB(_p_Vec*, void*, int,int){ +void Sing::GradjB(Vec, int,int){ ISSMERROR(" not supported yet!"); } /*}}}*/ /*FUNCTION Sing::GradjDrag {{{1*/ -void Sing::GradjDrag(_p_Vec*, void*, int,int){ +void Sing::GradjDrag(Vec, int,int){ ISSMERROR(" not supported yet!"); } @@ -534,5 +533,5 @@ /*}}}*/ /*FUNCTION Sing::Misfit {{{1*/ -double Sing::Misfit(void*, int,int){ +double Sing::Misfit( int,int){ ISSMERROR(" not supported yet!"); } @@ -545,5 +544,5 @@ /*}}}*/ /*FUNCTION Sing::SurfaceArea {{{1*/ -double Sing::SurfaceArea(void*, int,int){ +double Sing::SurfaceArea( int,int){ ISSMERROR(" not supported yet!"); } Index: /issm/trunk/src/c/objects/Sing.h =================================================================== --- /issm/trunk/src/c/objects/Sing.h (revision 3619) +++ /issm/trunk/src/c/objects/Sing.h (revision 3620) @@ -12,6 +12,5 @@ #include "../ModelProcessorx/IoModel.h" #include "./Hook.h" - -class Hook; +#include "../DataSet/Inputs.h" class Sing: public Element{ @@ -26,5 +25,4 @@ Hook hmatpar; //hook to 1 matpar Hook hnumpar; //hook to 1 numpar - Inputs* inputs; @@ -39,5 +37,5 @@ /*}}}*/ /*object management: {{{1*/ - void Configure(void* loads,void* nodes,void* materials,void* parameters); + void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); Object* copy(); void DeepEcho(); @@ -72,11 +70,11 @@ void GetBedList(double*); void GetThicknessList(double* thickness_list); - void Du(_p_Vec*,void*,int,int); - void Gradj(_p_Vec*, void*, int, int,char*); - void GradjDrag(_p_Vec*, void*, int,int); - void GradjB(_p_Vec*, void*, int,int); - double Misfit(void*,int,int); - double SurfaceArea(void*,int,int); - double CostFunction(void*,int,int); + void Du(Vec,int,int); + void Gradj(Vec, int, int,char*); + void GradjDrag(Vec , int,int); + void GradjB(Vec, int,int); + double Misfit(int,int); + double SurfaceArea(int,int); + double CostFunction(int,int); double MassFlux(double* segment,double* ug); /*}}}*/ Index: /issm/trunk/src/c/objects/SingVertexInput.cpp =================================================================== --- /issm/trunk/src/c/objects/SingVertexInput.cpp (revision 3620) +++ /issm/trunk/src/c/objects/SingVertexInput.cpp (revision 3620) @@ -0,0 +1,139 @@ +/*!\file SingVertexInput.c + * \brief: implementation of the SingVertexInput 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 "./SingVertexInput.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 SingVertexInput::SingVertexInput(){{{1*/ +SingVertexInput::SingVertexInput(){ + return; +} +/*}}}*/ +/*FUNCTION SingVertexInput::SingVertexInput(double* values){{{1*/ +SingVertexInput::SingVertexInput(int in_enum_type,double in_value){ + + enum_type=in_enum_type; + value=in_value; +} +/*}}}*/ +/*FUNCTION SingVertexInput::~SingVertexInput(){{{1*/ +SingVertexInput::~SingVertexInput(){ + return; +} +/*}}}*/ + +/*Object management*/ +/*FUNCTION SingVertexInput::copy{{{1*/ +Object* SingVertexInput::copy() { + + return new SingVertexInput(this->enum_type,this->value); + +} +/*}}}*/ +/*FUNCTION SingVertexInput::DeepEcho{{{1*/ +void SingVertexInput::DeepEcho(void){ + + printf("SingVertexInput:\n"); + printf(" enum: %i\n",this->enum_type); + printf(" %g\n",this->value); +} +/*}}}*/ +/*FUNCTION SingVertexInput::Demarshall{{{1*/ +void SingVertexInput::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 SingVertexInput::Echo {{{1*/ +void SingVertexInput::Echo(void){ + this->DeepEcho(); +} +/*}}}*/ +/*FUNCTION SingVertexInput::Enum{{{1*/ +int SingVertexInput::Enum(void){ + + return SingVertexInputEnum; + +} +/*}}}*/ +/*FUNCTION SingVertexInput::EnumType{{{1*/ +int SingVertexInput::EnumType(void){ + + return this->enum_type; + +} +/*}}}*/ +/*FUNCTION SingVertexInput::GetId{{{1*/ +int SingVertexInput::GetId(void){ return -1; } +/*}}}*/ +/*FUNCTION SingVertexInput::GetName{{{1*/ +char* SingVertexInput::GetName(void){ + return "triavertexinput"; +} +/*}}}*/ +/*FUNCTION SingVertexInput::Marshall{{{1*/ +void SingVertexInput::Marshall(char** pmarshalled_dataset){ + + char* marshalled_dataset=NULL; + int enum_type=0; + + /*recover marshalled_dataset: */ + marshalled_dataset=*pmarshalled_dataset; + + /*get enum type of SingVertexInput: */ + enum_type=SingVertexInputEnum; + + /*marshall enum: */ + memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type); + + /*marshall SingVertexInput 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 SingVertexInput::MarshallSize{{{1*/ +int SingVertexInput::MarshallSize(){ + + return sizeof(value)+ + sizeof(enum_type)+ + +sizeof(int); //sizeof(int) for enum type +} +/*}}}*/ +/*FUNCTION SingVertexInput::MyRank{{{1*/ +int SingVertexInput::MyRank(void){ + extern int my_rank; + return my_rank; +} +/*}}}*/ + +/*Object functions*/ + Index: /issm/trunk/src/c/objects/SingVertexInput.h =================================================================== --- /issm/trunk/src/c/objects/SingVertexInput.h (revision 3620) +++ /issm/trunk/src/c/objects/SingVertexInput.h (revision 3620) @@ -0,0 +1,42 @@ +/*! \file SingVertexInput.h + * \brief: header file for triavertexinput object + */ + +#include "./Input.h" + +#ifndef _SINGVERTEXINPUT_H_ +#define _SINGVERTEXINPUT_H_ + +class SingVertexInput: public Input{ + + private: + /*just hold 1 value for 1 vertex: */ + int enum_type; + double value; + + public: + + /*constructors, destructors: {{{1*/ + SingVertexInput(); + SingVertexInput(int enum_type,double value); + ~SingVertexInput(); + /*}}}*/ + /*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 /* _SINGVERTEXINPUT_H */ Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 3619) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 3620) @@ -30,25 +30,23 @@ } /*}}}*/ -/*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): +/*FUNCTION Tria::Tria(int id, int* node_ids, int matice_id, int matpar_id){{{1*/ +Tria::Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id): hnodes(tria_node_ids,3), hmatice(&tria_matice_id,1), - hmatpar(&tria_matpar_id,1), - hnumpar(&tria_numpar_id,1) + hmatpar(&tria_matpar_id,1) { /*all the initialization has been done by the initializer, just fill in the id: */ this->id=tria_id; + this->parameters=NULL; this->inputs=new Inputs(); - return; -} -/*}}}*/ -/*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): +} +/*}}}*/ +/*FUNCTION Tria::Tria(int id, Hook* hnodes, Hook* hmatice, Hook* hmatpar, DataSet* parameters, Inputs* tria_inputs) {{{1*/ +Tria::Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* tria_parameters, Inputs* tria_inputs): hnodes(tria_hnodes), hmatice(tria_hmatice), - hmatpar(tria_hmatpar), - hnumpar(tria_hnumpar) + hmatpar(tria_hmatpar) { @@ -61,6 +59,6 @@ this->inputs=new Inputs(); } - - return; + /*point parameters: */ + this->parameters=tria_parameters; } /*}}}*/ @@ -72,5 +70,4 @@ int tria_matice_id; int tria_matpar_id; - int tria_numpar_id; double nodeinputs[3]; @@ -94,10 +91,8 @@ 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 this->hnodes.Init(tria_node_ids,3); this->hmatice.Init(&tria_matice_id,1); this->hmatpar.Init(&tria_matpar_id,1); - this->hnumpar.Init(&tria_numpar_id,1); //intialize inputs, and add as many inputs per element as requested: @@ -143,4 +138,7 @@ if (iomodel->elementonsurface) this->inputs->AddInput(new BoolInput(ElementOnSurfaceEnum,(IssmBool)iomodel->elementonsurface[index])); + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + } @@ -149,4 +147,5 @@ Tria::~Tria(){ delete inputs; + this->parameters=NULL; return; } @@ -155,18 +154,7 @@ /*Object management: */ /*FUNCTION Tria::Configure {{{1*/ -void Tria::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){ +void Tria::Configure(DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ int i; - - DataSet* loadsin=NULL; - DataSet* nodesin=NULL; - DataSet* materialsin=NULL; - DataSet* parametersin=NULL; - - /*Recover pointers :*/ - loadsin=(DataSet*)ploadsin; - nodesin=(DataSet*)pnodesin; - materialsin=(DataSet*)pmaterialsin; - parametersin=(DataSet*)pparametersin; /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective @@ -175,5 +163,7 @@ hmatice.configure(materialsin); hmatpar.configure(materialsin); - hnumpar.configure(parametersin); + + /*point parameters to real dataset: */ + this->parameters=parametersin; } @@ -182,5 +172,5 @@ Object* Tria::copy() { - return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,&this->hnumpar,this->inputs); + return new Tria(this->id,&this->hnodes,&this->hmatice,&this->hmatpar,this->parameters,this->inputs); } @@ -205,8 +195,10 @@ hmatice.Demarshall(&marshalled_dataset); hmatpar.Demarshall(&marshalled_dataset); - hnumpar.Demarshall(&marshalled_dataset); /*demarshall inputs: */ inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); + + /*parameters: may not exist even yet, so let Configure handle it: */ + this->parameters=NULL; /*return: */ @@ -224,8 +216,9 @@ hmatice.DeepEcho(); hmatpar.DeepEcho(); - hnumpar.DeepEcho(); + printf(" parameters\n"); + parameters->DeepEcho(); printf(" inputs\n"); inputs->DeepEcho(); - + return; } @@ -240,9 +233,8 @@ hmatice.Echo(); hmatpar.Echo(); - hnumpar.Echo(); + printf(" parameters\n"); + parameters->Echo(); printf(" inputs\n"); inputs->Echo(); - - return; } /*}}}*/ @@ -271,5 +263,4 @@ hmatice.Marshall(&marshalled_dataset); hmatpar.Marshall(&marshalled_dataset); - hnumpar.Marshall(&marshalled_dataset); /*Marshall inputs: */ @@ -279,4 +270,8 @@ marshalled_dataset+=marshalled_inputs_size; + /*parameters: don't do anything about it. parameters are marshalled somewhere else!*/ + + xfree((void**)&marshalled_inputs); + *pmarshalled_dataset=marshalled_dataset; return; @@ -290,5 +285,4 @@ +hmatice.MarshallSize() +hmatpar.MarshallSize() - +hnumpar.MarshallSize() +inputs->MarshallSize() +sizeof(int); //sizeof(int) for enum type @@ -312,5 +306,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -318,5 +311,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*Update internal data if inputs holds new values: */ @@ -468,5 +460,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -474,5 +465,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*plug local pressure values into global pressure vector: */ @@ -497,5 +487,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -503,5 +492,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*Get dof list on which we will plug the pressure values: */ @@ -536,5 +524,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -542,5 +529,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*plug local pressure values into global pressure vector: */ @@ -583,4 +569,6 @@ double dk[NDOF2]; double dB[NDOF2]; + char* control_type=NULL; + double cm_noisedmp; /* Jacobian: */ @@ -591,5 +579,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*inputs: */ @@ -608,5 +595,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&control_type,"control_type"); + this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); /* Get node coordinates and dof list: */ @@ -631,17 +621,17 @@ /*Add Tikhonov regularization term to misfit*/ - if (strcmp(numpar->control_type,"drag")==0){ + if (strcmp(control_type,"drag")==0){ 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; + Jelem+=cm_noisedmp*1/2*(pow(dk[0],2)+pow(dk[1],2))*Jdet*gauss_weight; } } - else if (strcmp(numpar->control_type,"B")==0){ + else if (strcmp(control_type,"B")==0){ 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; + Jelem+=cm_noisedmp*1/2*(pow(dB[0],2)+pow(dB[1],2))*Jdet*gauss_weight; } else{ - ISSMERROR("%s%s","unsupported control type: ",numpar->control_type); + ISSMERROR("%s%s","unsupported control type: ",control_type); } @@ -652,4 +642,5 @@ xfree((void**)&third_gauss_area_coord); xfree((void**)&gauss_weights); + xfree((void**)&control_type); /*Return: */ @@ -763,5 +754,7 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; + + /*parameters: */ + double artdiff; @@ -770,5 +763,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -780,6 +772,9 @@ inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); + /*retrieve some parameters: */ + this->parameters->FindParam(&artdiff,"artdiff"); + //Create Artificial diffusivity once for all if requested - if(numpar->artdiff){ + if(artdiff){ //Get the Jacobian determinant gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; @@ -848,5 +843,5 @@ for( i=0; iartdiff){ + if(artdiff){ /* Compute artificial diffusivity */ @@ -918,5 +913,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; @@ -925,5 +919,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -1034,9 +1027,11 @@ int found=0; + /*parameters: */ + double artdiff; + /*dynamic objects pointed to by hooks: */ Node** nodes=NULL; Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -1044,9 +1039,11 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*Recover velocity: */ inputs->GetParameterValues(&vx_list[0],&gaussgrids[0][0],3,VxAverageEnum); inputs->GetParameterValues(&vy_list[0],&gaussgrids[0][0],3,VyAverageEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&artdiff,"artdiff"); /* Get node coordinates and dof list: */ @@ -1075,5 +1072,5 @@ //Create Artificial diffusivity once for all if requested - if(numpar->artdiff){ + if(artdiff){ //Get the Jacobian determinant gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; @@ -1131,5 +1128,5 @@ for( i=0; iartdiff){ + if(artdiff){ /* Compute artificial diffusivity */ @@ -1195,9 +1192,12 @@ double B[3][numdof]; double Bprime[3][numdof]; - double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}}; // material matrix, simple scalar matrix. + double D[3][3]={0.0}; // material matrix, simple scalar matrix. double D_scalar; + /*parameters: */ + double viscosity_overshoot; + /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. @@ -1212,5 +1212,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*inputs: */ @@ -1220,4 +1219,7 @@ inputs->GetParameterValue(&onwater,ElementOnWaterEnum); inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&viscosity_overshoot,"viscosity_overshoot"); /*First, if we are on water, return empty matrix: */ @@ -1228,12 +1230,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;iviscosity_overshoot*(viscosity-oldviscosity); + newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity); D_scalar=newviscosity*thickness*gauss_weight*Jdet; @@ -1333,5 +1331,5 @@ /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag @@ -1365,5 +1363,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; @@ -1372,5 +1369,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*retrieve inputs :*/ @@ -1381,7 +1377,4 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;idt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); + this->parameters->FindParam(&artdiff,"artdiff"); /* Get node coordinates and dof list: */ @@ -1787,5 +1775,5 @@ //Create Artificial diffusivity once for all if requested - if(numpar->artdiff==1){ + if(artdiff==1){ //Get the Jacobian determinant gauss_l1l2l3[0]=ONETHIRD; gauss_l1l2l3[1]=ONETHIRD; gauss_l1l2l3[2]=ONETHIRD; @@ -1867,5 +1855,5 @@ for( i=0; iartdiff==1){ + if(artdiff==1){ /* Compute artificial diffusivity */ @@ -1934,7 +1922,9 @@ /*input parameters for structural analysis (diagnostic): */ double vx,vy; - double dt; int dofs[1]={0}; int found; + + /*parameters: */ + double dt; /*dynamic objects pointed to by hooks: */ @@ -1942,5 +1932,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -1948,8 +1937,7 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*recover dt: */ - dt=numpar->dt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); /* Get node coordinates and dof list: */ @@ -2050,5 +2038,5 @@ /* local element matrices: */ - double Ke_gg[numdof][numdof]; //local element stiffness matrix + double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. @@ -2059,5 +2047,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2065,12 +2052,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set Ke_gg to 0: */ - for(i=0;idt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); /* Get node coordinates and dof list: */ @@ -2304,5 +2287,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2310,5 +2292,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -2391,5 +2372,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2397,5 +2377,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -2478,5 +2457,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2484,11 +2462,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*recover objects from hooks: */ - nodes=(Node**)hnodes.deliverp(); - matpar=(Matpar*)hmatpar.delivers(); - matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -2562,5 +2533,5 @@ /*element vector at the gaussian points: */ - double pe_g[numdof]; + double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; @@ -2578,5 +2549,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2584,12 +2554,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set pe_g to 0: */ - for(i=0;idt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); /* Get node coordinates and dof list: */ @@ -2909,5 +2870,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -2915,9 +2875,7 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); - - /*recover dt: */ - dt=numpar->dt; - + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); /* Get node coordinates and dof list: */ @@ -2993,5 +2951,5 @@ /*element vector at the gaussian points: */ - double pe_g[numdof]; + double pe_g[numdof]={0.0}; double pe_g_gaussian[numdof]; double slope[2]; @@ -3001,5 +2959,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -3007,12 +2964,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList(&doflist[0],&numberofdofspernode); - - /* Set pe_g to 0: */ - for(i=0;iGetBeta(); meltingpoint=matpar->GetMeltingPoint(); - dt=numpar->dt; - + + /*retrieve some solution parameters: */ + this->parameters->FindParam(&dt,"dt"); /* Ice/ocean heat exchange flux on ice shelf base */ @@ -3231,5 +3183,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -3237,5 +3188,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*retrieve inputs :*/ @@ -3249,5 +3199,7 @@ rho_ice=matpar->GetRhoIce(); heatcapacity=matpar->GetHeatCapacity(); - dt=numpar->dt; + + /*retrieve some parameters: */ + this->parameters->FindParam(&dt,"dt"); @@ -3369,7 +3321,8 @@ double obs_velocity_mag,velocity_mag; double dux,duy; + double meanvel, epsvel; /*element vector : */ - double due_g[numdof]={0,0,0,0,0,0}; + double due_g[numdof]={0.0}; double due_g_gaussian[numdof]; @@ -3391,5 +3344,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -3397,5 +3349,8 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,"meanvel"); + this->parameters->FindParam(&epsvel,"epsvel"); /* Get node coordinates and dof list: */ @@ -3460,6 +3415,6 @@ */ for (i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),2); - scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),2); + scalex=pow(meanvel/(obs_vx_list[i]+epsvel),2); + scaley=pow(meanvel/(obs_vy_list[i]+epsvel),2); if(obs_vx_list[i]==0)scalex=0; if(obs_vy_list[i]==0)scaley=0; @@ -3482,7 +3437,7 @@ */ for (i=0;iepsvel; //epsvel to avoid velocity being nil. - obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. - scale=-8*pow(numpar->meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); + velocity_mag=sqrt(pow(vx_list[i],2)+pow(vy_list[i],2))+epsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_vx_list[i],2)+pow(obs_vy_list[i],2))+epsvel; //epsvel to avoid observed velocity being nil. + scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag); dux_list[i]=scale*vx_list[i]; duy_list[i]=scale*vy_list[i]; @@ -3501,5 +3456,5 @@ */ for (i=0;iepsvel); + scale=1.0/(S*sqrt(pow(vx_list[i]-obs_vx_list[i],2)+pow(vy_list[i]-obs_vx_list[i],2))+epsvel); dux_list[i]=scale*(obs_vx_list[i]-vx_list[i]); duy_list[i]=scale*(obs_vy_list[i]-vy_list[i]); @@ -3518,8 +3473,8 @@ */ for (i=0;imeanvel,(double)2)*( - log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)) * 1/(vx_list[i]+numpar->epsvel)); - duy_list[i] = - pow(numpar->meanvel,(double)2)*( - log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)) * 1/(vy_list[i]+numpar->epsvel)); + dux_list[i] = - pow(meanvel,(double)2)*( + log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)) * 1/(vx_list[i]+epsvel)); + duy_list[i] = - pow(meanvel,(double)2)*( + log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)) * 1/(vy_list[i]+epsvel)); } } @@ -4228,5 +4183,5 @@ /*element vector at the gaussian points: */ - double grade_g[numgrids]; + double grade_g[numgrids]={0.0}; double grade_g_gaussian[numgrids]; @@ -4252,4 +4207,11 @@ double dB[NDOF2]; double B_gauss; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; /*dynamic objects pointed to by hooks: */ @@ -4257,5 +4219,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -4263,12 +4224,15 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); + this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); + this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); + this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); + this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids); GetDofList1(&doflist1[0]); - - /* Set grade_g to 0: */ - for(i=0;icm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]); + grade_g_gaussian[i]-=cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dB[0]+dh1dh3[1][i]*dB[1]); //min dampening - if(B_gausscm_mindmp_value){ - grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(B_gaussnumpar->cm_maxdmp_value){ - grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(B_gauss>cm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; } @@ -4394,5 +4358,5 @@ /*element vector at the gaussian points: */ - double grade_g[numgrids]={0,0,0}; + double grade_g[numgrids]={0.0}; double grade_g_gaussian[numgrids]; @@ -4408,4 +4372,11 @@ /*inputs: */ bool shelf; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; /*dynamic objects pointed to by hooks: */ @@ -4413,5 +4384,4 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*recover objects from hooks: */ @@ -4419,8 +4389,14 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /*retrieve inputs :*/ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); + this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); + this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); + this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); + this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); /*Get out if shelf*/ @@ -4520,14 +4496,14 @@ //noise dampening d/dki(1/2*(dk/dx)^2) - grade_g_gaussian[i]+=-numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + grade_g_gaussian[i]+=-cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); //min dampening - if(dragcm_mindmp_value){ - grade_g_gaussian[i]+=numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(dragnumpar->cm_maxdmp_value){ - grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(drag>cm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; } } @@ -4564,13 +4540,11 @@ double vy_list[numgrids]; double vz_list[numgrids]; - double vxvyvz_list[numgrids][3]; - double adjx_list[numgrids]; - double adjy_list[numgrids]; - double adjz_list[numgrids]; - double adjxyz_list[numgrids][3]; - double drag; - int dofs1[1]={0}; - int dofs3[3]={0,1,2}; + double thickness_list[numgrids]; + double bed_list[numgrids]; + double dragcoefficient_list[numgrids]; + double drag_p,drag_q; + double alpha_complement_list[numgrids]; + double alpha_complement; /* gaussian points: */ @@ -4582,4 +4556,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* parameters: */ @@ -4591,11 +4566,6 @@ double dk[NDOF2]; - /*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}; double grade_g_gaussian[numgrids]; @@ -4612,4 +4582,11 @@ bool shelf; int drag_type; + + /*parameters: */ + double cm_noisedmp; + double cm_mindmp_slope; + double cm_mindmp_value; + double cm_maxdmp_value; + double cm_maxdmp_slope; /*dynamic objects pointed to by hooks: */ @@ -4617,9 +4594,15 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*retrieve inputs :*/ inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum); inputs->GetParameterValue(&drag_type,DragTypeEnum); + + /*retrieve some parameters: */ + this->parameters->FindParam(&cm_noisedmp,"cm_noisedmp"); + this->parameters->FindParam(&cm_mindmp_value,"cm_mindmp_value"); + this->parameters->FindParam(&cm_mindmp_slope,"cm_mindmp_slope"); + this->parameters->FindParam(&cm_maxdmp_value,"cm_maxdmp_value"); + this->parameters->FindParam(&cm_maxdmp_slope,"cm_maxdmp_slope"); /*Get out if shelf*/ @@ -4630,5 +4613,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -4636,26 +4618,48 @@ 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",&vxvyvz_list[0][0],3,dofs3,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity input parameter"); - } - if(!inputs->Recover("adjoint",&adjxyz_list[0][0],3,dofs3,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(&vz_list[0],&gaussgrids[0][0],3,VzAverageEnum); + 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); + + + /*Initialize all fields: */ + friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); + strcpy(friction->element_type,"2d"); + + friction->gravity=matpar->GetG(); + friction->rho_ice=matpar->GetRhoIce(); + friction->rho_water=matpar->GetRhoWater(); + 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!"); + + /*Compute alpha complement: */ + FrictionGetAlphaComplement(&alpha_complement_list[0],friction); + + /*Erase friction object: */ + DeleteFriction(&friction); + } + else{ + alpha_complement_list[0]=0; + alpha_complement_list[1]=0; + alpha_complement_list[2]=0; } @@ -4671,51 +4675,17 @@ gauss_l1l2l3[2]=*(third_gauss_area_coord+ig); - /*Build alpha_complement_list: */ - if (drag_type==2){ - - /*Allocate friction object: */ - Friction* friction=NewFriction(); - - /*Initialize all fields: */ - friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char)); - 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=&vxvyvz_list[0][0]; - friction->p=this->properties.p; - friction->q=this->properties.q; - - if(friction->p!=1) ISSMERROR("non-linear friction not supported yet in control methods!"); - - /*Compute alpha complement: */ - FrictionGetAlphaComplement(&alpha_complement_list[0],friction); - - /*Erase friction object: */ - DeleteFriction(&friction); - } - else{ - alpha_complement_list[0]=0; - alpha_complement_list[1]=0; - alpha_complement_list[2]=0; - } - /*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[0],DragCoefficientEnum); /*recover lambda mu and xi: */ - GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3); - GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3); - GetParameterValue(&xi, &adjz_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&lambda, &gauss_l1l2l3[0],AdjointxEnum); + inputs->GetParameterValue(&mu, &gauss_l1l2l3[0],AdjointyEnum); + inputs->GetParameterValue(&xi, &gauss_l1l2l3[0],AdjointzEnum); /*recover vx vy and vz: */ - GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3); - GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3); - GetParameterValue(&vz, &vz_list[0],gauss_l1l2l3); + inputs->GetParameterValue(&vx, &gauss_l1l2l3[0],VxEnum); + inputs->GetParameterValue(&vy, &gauss_l1l2l3[0],VyEnum); + inputs->GetParameterValue(&vz, &gauss_l1l2l3[0],VzEnum); /*Get normal vecyor to the bed */ @@ -4736,5 +4706,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): */ @@ -4748,14 +4718,14 @@ //Add regularization term - grade_g_gaussian[i]+= - numpar->cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); + grade_g_gaussian[i]+= - cm_noisedmp*Jdet*gauss_weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]); //min dampening - if(dragcm_mindmp_value){ - grade_g_gaussian[i]+= numpar->cm_mindmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(dragnumpar->cm_maxdmp_value){ - grade_g_gaussian[i]+= - numpar->cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; + if(drag>cm_maxdmp_value){ + grade_g_gaussian[i]+= - cm_maxdmp_slope*Jdet*gauss_weight*l1l2l3[i]; } } @@ -4839,11 +4809,10 @@ /*get thickness and velocity at two segment extremities: */ - GetParameterValue(&h1, &this->properties.h[0],gauss_1); - GetParameterValue(&h2, &this->properties.h[0],gauss_2); - GetParameterValue(&vx1, &vx_list[0],gauss_1); - GetParameterValue(&vy1, &vy_list[0],gauss_1); - GetParameterValue(&vx2, &vx_list[0],gauss_2); - GetParameterValue(&vy2, &vy_list[0],gauss_2); - + inputs->GetParameterValue(&h1, &gauss_1[0],ThicknessEnum); + inputs->GetParameterValue(&h2, &gauss_2[0],ThicknessEnum); + inputs->GetParameterValue(&vx1, &gauss_1[0],VxEnum); + inputs->GetParameterValue(&vx2, &gauss_2[0],VxEnum); + inputs->GetParameterValue(&vy1, &gauss_1[0],VyEnum); + inputs->GetParameterValue(&vy2, &gauss_2[0],VyEnum); mass_flux= rho_ice*length*( @@ -4871,8 +4840,6 @@ /* 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]; @@ -4888,4 +4855,5 @@ double gauss_weight; double gauss_l1l2l3[3]; + double gaussgrids[numgrids][numgrids]={{1,0,0},{0,1,0},{0,0,1}}; /* parameters: */ @@ -4906,8 +4874,9 @@ Matpar* matpar=NULL; Matice* matice=NULL; - Numpar* numpar=NULL; /*inputs: */ bool onwater; + + double meanvel, epsvel; /*retrieve inputs :*/ @@ -4921,5 +4890,4 @@ matpar=(Matpar*)hmatpar.delivers(); matice=(Matice*)hmatice.delivers(); - numpar=(Numpar*)hnumpar.delivers(); /* Get node coordinates and dof list: */ @@ -4927,26 +4895,18 @@ /* Recover input data: */ - if(!inputs->Recover("fit",&fit)) ISSMERROR(" missing fit input parameter"); - if(fit==3 && !inputs->Recover("surfacearea",&S)){ - ISSMERROR("missing surface area input parameter"); - } - if(!inputs->Recover("velocity_obs",&obs_vxvy_list[0][0],2,dofs2,numgrids,(void**)nodes)){ - ISSMERROR("missing velocity_obs 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"); - } - - /*Initialize velocities: */ - for(i=0;iGetParameterValue(&fit,FitEnum); + if(fit==3){ + inputs->GetParameterValue(&S,SurfaceAreaEnum); + } + inputs->GetParameterValues(&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); + + /*retrieve some parameters: */ + this->parameters->FindParam(&meanvel,"meanvel"); + this->parameters->FindParam(&epsvel,"epsvel"); + /* Compute Misfit at the 3 nodes * Here we integrate linearized functions: @@ -4979,6 +4939,6 @@ */ for (i=0;imeanvel/(obs_vx_list[i]+numpar->epsvel),(double)2); - scaley=pow(numpar->meanvel/(obs_vy_list[i]+numpar->epsvel),(double)2); + scalex=pow(meanvel/(obs_vx_list[i]+epsvel),(double)2); + scaley=pow(meanvel/(obs_vy_list[i]+epsvel),(double)2); if(obs_vx_list[i]==0)scalex=0; if(obs_vy_list[i]==0)scaley=0; @@ -4995,7 +4955,7 @@ */ for (i=0;iepsvel; //epsvel to avoid velocity being nil. - obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+numpar->epsvel; //epsvel to avoid observed velocity being nil. - misfit_list[i]=4*pow(numpar->meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); + velocity_mag=sqrt(pow(vx_list[i],(double)2)+pow(vy_list[i],(double)2))+epsvel; //epsvel to avoid velocity being nil. + obs_velocity_mag=sqrt(pow(obs_vx_list[i],(double)2)+pow(obs_vy_list[i],(double)2))+epsvel; //epsvel to avoid observed velocity being nil. + misfit_list[i]=4*pow(meanvel,(double)2)*pow(log(velocity_mag/obs_velocity_mag),(double)2); } } @@ -5020,7 +4980,7 @@ */ for (i=0;imeanvel,(double)2)*( - pow(log((fabs(vx_list[i])+numpar->epsvel)/(fabs(obs_vx_list[i])+numpar->epsvel)),(double)2) + - pow(log((fabs(vy_list[i])+numpar->epsvel)/(fabs(obs_vy_list[i])+numpar->epsvel)),(double)2) ); + misfit_list[i]=0.5*pow(meanvel,(double)2)*( + pow(log((fabs(vx_list[i])+epsvel)/(fabs(obs_vx_list[i])+epsvel)),(double)2) + + pow(log((fabs(vy_list[i])+epsvel)/(fabs(obs_vy_list[i])+epsvel)),(double)2) ); } } Index: /issm/trunk/src/c/objects/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Tria.h (revision 3619) +++ /issm/trunk/src/c/objects/Tria.h (revision 3620) @@ -31,7 +31,7 @@ Hook hmatice; //hook to 1 matice Hook hmatpar; //hook to 1 matpar - Hook hnumpar; //hook to 1 numpar - Inputs* inputs; + Parameters* parameters; //pointer to solution parameters + Inputs* inputs; public: @@ -39,11 +39,11 @@ /*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); - Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Hook* tria_hnumpar, Inputs* tria_inputs); + Tria(int tria_id,int* tria_node_ids, int tria_matice_id, int tria_matpar_id); + Tria(int tria_id,Hook* tria_hnodes, Hook* tria_hmatice, Hook* tria_hmatpar, Parameters* parameters, Inputs* tria_inputs); Tria(int i, IoModel* iomodel); ~Tria(); /*}}}*/ /*FUNCTION object management {{{1*/ - void Configure(void* loads,void* nodes,void* materials,void* parameters); + void Configure(DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); Object* copy(); void DeepEcho();