Index: /issm/trunk/src/c/DataSet/DataSet.cpp =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.cpp (revision 245) +++ /issm/trunk/src/c/DataSet/DataSet.cpp (revision 246) @@ -881,5 +881,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::CreateKMatrix" -void DataSet::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void DataSet::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ vector::iterator object; @@ -905,5 +905,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::CreatePVector" -void DataSet::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){ +void DataSet::CreatePVector(Vec pg,void* inputs,int analysis_type){ vector::iterator object; @@ -929,5 +929,5 @@ #undef __FUNCT__ #define __FUNCT__ "UpdateFromInputs" -void DataSet::UpdateFromInputs(ParameterInputs* inputs){ +void DataSet::UpdateFromInputs(void* inputs){ vector::iterator object; @@ -965,5 +965,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::PenaltyCreateKMatrix" -void DataSet::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){ +void DataSet::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ vector::iterator object; @@ -983,5 +983,5 @@ #undef __FUNCT__ #define __FUNCT__ "DataSet::PenaltyCreatePVector" -void DataSet::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){ +void DataSet::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ vector::iterator object; @@ -1009,5 +1009,5 @@ #undef __FUNCT__ #define __FUNCT__ "RiftConstraints" -void DataSet::RiftConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){ +void DataSet::RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){ throw ErrorException(__FUNCT__," not implemented yet!"); @@ -1026,5 +1026,5 @@ #undef __FUNCT__ #define __FUNCT__ "MeltingConstraints" -void DataSet::MeltingConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type){ +void DataSet::MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type){ throw ErrorException(__FUNCT__," not implemented yet!"); @@ -1065,5 +1065,5 @@ -void DataSet::Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){ +void DataSet::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){ @@ -1083,5 +1083,5 @@ } -void DataSet::Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type){ +void DataSet::Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type){ @@ -1101,5 +1101,5 @@ } -void DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){ +void DataSet::Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type){ double J=0;; Index: /issm/trunk/src/c/DataSet/DataSet.h =================================================================== --- /issm/trunk/src/c/DataSet/DataSet.h (revision 245) +++ /issm/trunk/src/c/DataSet/DataSet.h (revision 246) @@ -10,16 +10,14 @@ #define DATASET_H_ -#include "../objects/Object.h" -#include "../objects/ParameterInputs.h" #include #include "../toolkits/toolkits.h" +#include "../objects/objects.h" -using namespace std; - -class DataSet { +class DataSet{ private: + /*internals: */ - vector objects; + std::vector objects; /*type of dataset: */ @@ -63,24 +61,21 @@ void SetSorting(int* in_sorted_ids,int* in_id_offsets); void Sort(); - void CreateKMatrix(Mat Kgg,ParameterInputs* inputs, int analysis_type); - void CreatePVector(Vec pg,ParameterInputs* inputs, int analysis_type); - void UpdateFromInputs(ParameterInputs* inputs); - void PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs, int analysis_type); + void CreatePVector(Vec pg,void* inputs, int analysis_type); + void UpdateFromInputs(void* inputs); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); int RiftIsPresent(); - void RiftConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type); + void RiftConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type); int MeltingIsPresent(); - void MeltingConstraints(int* pnum_unstable_constraints,ParameterInputs* inputs,int analysis_type); + void MeltingConstraints(int* pnum_unstable_constraints,void* inputs,int analysis_type); DataSet* Copy(void); - void Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type); - void Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type); - void Misfit(double* pJ, double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type); + void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type); + void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type); + void Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type); void VelocityExtrude(Vec ug,double* ug_serial); int DeleteObject(Object* object); - - }; - /*This routine cannot be object oriented, but need for demarshalling: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 245) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 246) @@ -57,4 +57,8 @@ int StokesEnum(void){ return 4; } +/*inputs: */ +int InputEnum(void){ return 310; } + + /*functions on enums: */ Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 245) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 246) @@ -21,4 +21,5 @@ int LoadEnum(void); int RgbEnum(void); +int InputEnum(void); /*formulations: */ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 245) +++ /issm/trunk/src/c/Makefile.am (revision 246) @@ -42,4 +42,6 @@ ./objects/Matpar.h\ ./objects/Matpar.cpp\ + ./objects/Input.h\ + ./objects/Input.cpp\ ./objects/ParameterInputs.h\ ./objects/ParameterInputs.cpp\ @@ -262,4 +264,6 @@ ./objects/Matpar.h\ ./objects/Matpar.cpp\ + ./objects/Input.h\ + ./objects/Input.cpp\ ./objects/ParameterInputs.h\ ./objects/ParameterInputs.cpp\ Index: /issm/trunk/src/c/UpdateFromInputsx/UpdateFromInputsx.cpp =================================================================== --- /issm/trunk/src/c/UpdateFromInputsx/UpdateFromInputsx.cpp (revision 245) +++ /issm/trunk/src/c/UpdateFromInputsx/UpdateFromInputsx.cpp (revision 246) @@ -26,5 +26,5 @@ loads->Configure(elements, loads, nodes, materials); nodes->Configure(elements, loads, nodes, materials); - + /*Update elements, nodes, loads and materials from inputs: */ elements->UpdateFromInputs(inputs); Index: /issm/trunk/src/c/io/ParameterInputsInit.cpp =================================================================== --- /issm/trunk/src/c/io/ParameterInputsInit.cpp (revision 245) +++ /issm/trunk/src/c/io/ParameterInputsInit.cpp (revision 246) @@ -9,5 +9,4 @@ #endif -#ifdef _SERIAL_ #include "../objects/objects.h" @@ -18,59 +17,81 @@ #undef __FUNCT__ #define __FUNCT__ "ParameterInputsInit" -void ParameterInputsInit(ParameterInputs** pinputs, const mxArray* inputhandle){ - - int dummy; + +#ifdef _SERIAL_ +void ParameterInputs::Init( void* vinput_handle){ + int i; int nelems; - int field_name_length; + mxArray* input=NULL; + const char* field_name=NULL; + ConstDataHandle input_handle=NULL; + int type; + + char* string=NULL; + char* name=NULL; + double scalar; + int integer; + int ndof; + int numberofnodes; + double* vector=NULL; + + input_handle=(ConstDataHandle)vinput_handle; - /*output: */ - ParameterInputs* inputs=NULL; + /*number of inputs: */ + nelems=mxGetNumberOfElements(input_handle); - /*First allocate ParameterInputs: */ - inputs=NewParameterInputs(); + /*Go through ParameterInputs matlab structure, and recover each input object. Use it for fill + * the ParameterInputs* "C" class object: */ - /*Now populate the structure: */ + for(i=0;infields = mxGetNumberOfFields(inputhandle); //tells us how many parameters we have - - /*Check that we don't have a multi-dimensional structure array: */ - nelems = mxGetNumberOfElements(inputhandle); - - if(nelems!=1){ - _printf_("%s%s\n",__FUNCT__," error message: input parameter structure must be of size 1!"); - mexErrMsgTxt(" "); - } - - /*Allocate fields and field names: */ - inputs->fields=(double**)xmalloc(inputs->nfields*sizeof(double*)); - inputs->field_names=(char**)xmalloc(inputs->nfields*sizeof(char*)); - inputs->fields_length=(int*)xmalloc(inputs->nfields*sizeof(int)); + /*Get i'th matlab input: */ + input=mxGetField(input_handle,i,"input"); - /*Initialize, we never know: */ - for(i=0;infields;i++){ - inputs->fields[i]=NULL; - inputs->field_names[i]=NULL; - inputs->fields_length[i]=0; + /*This input has the following fields: name, type, string, integer, double, and vector. Plus ndof and numberofnodes: */ + //first get type + type=(int)*mxGetPr(mxGetField(input,0,"type")); + + /*Get name: */ + name=mxArrayToString(mxGetField(input,0,"name")); + + if (strcmp(name,"")==0)continue; //this input has no name! can't register it. keep going. + + //string + switch(type){ + case 0: + //we have a string, recover it. + string=mxArrayToString(mxGetField(input,0,"string")); + //add to inputs. + this->Add(name,string); + xfree((void**)&string); + break; + case 1: + //integer + integer=(int)*(mxGetPr(mxGetField(input,0,"integer"))); + this->Add(name,integer); + break; + case 2: + //double + scalar=(double)*(mxGetPr(mxGetField(input,0,"double"))); + this->Add(name,scalar); + break; + case 3: + //double vector + vector=mxGetPr(mxGetField(input,0,"vector")); + ndof=(int)*(mxGetPr(mxGetField(input,0,"ndof"))); + numberofnodes=(int)*(mxGetPr(mxGetField(input,0,"numberofnodes"))); + this->Add(name,vector,ndof,numberofnodes); + break; + default: + throw ErrorException(__FUNCT__,exprintf("%s%i"," unknow parameter input type ",type)); + } + } - /*Go through input structure fields and populate our ParameterInputs object:*/ - for (i=0;infields;i++){ - - /*Recover field value: */ - FetchData((void**)&inputs->fields[i],&inputs->fields_length[i],&dummy,mxGetFieldByNumber( inputhandle,0,i),"Matrix","Mat"); - - /*Recover field name: */ - const char* field_name=mxGetFieldNameByNumber(inputhandle,i); - field_name_length=strlen(field_name)+1; - inputs->field_names[i]=(char*)xmalloc(field_name_length*sizeof(char)); - strcpy(inputs->field_names[i],field_name); - } - - /*Assign output poiters:*/ - *pinputs=inputs; - } +#else +void ParameterInputs::Init( void* vinput_handle){ +} #endif //#ifdef _SERIAL- Index: /issm/trunk/src/c/objects/Element.cpp =================================================================== --- /issm/trunk/src/c/objects/Element.cpp (revision 245) +++ /issm/trunk/src/c/objects/Element.cpp (revision 246) @@ -5,9 +5,6 @@ -#include "./Object.h" #include "./Element.h" #include "../EnumDefinitions/EnumDefinitions.h" -#include "stdio.h" -#include "../DataSet/DataSet.h" int Element::Enum(void){ Index: /issm/trunk/src/c/objects/Element.h =================================================================== --- /issm/trunk/src/c/objects/Element.h (revision 245) +++ /issm/trunk/src/c/objects/Element.h (revision 246) @@ -10,8 +10,5 @@ #include "./Object.h" -#include "./ParameterInputs.h" #include "../toolkits/toolkits.h" -#include "./Node.h" -#include "./Matpar.h" class Element: public Object{ @@ -27,18 +24,18 @@ virtual void Demarshall(char** pmarshalled_dataset)=0; virtual void Configure(void* loads,void* nodes,void* materials)=0; - virtual void CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type)=0; - virtual void CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type)=0; - virtual void UpdateFromInputs(ParameterInputs* inputs)=0; - virtual void GetNodes(Node** nodes)=0; - virtual Matpar* GetMatPar()=0; + virtual void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type)=0; + virtual void CreatePVector(Vec pg, void* inputs, int analysis_type)=0; + virtual void UpdateFromInputs(void* inputs)=0; + virtual void GetNodes(void** nodes)=0; + virtual void* GetMatPar()=0; virtual int GetShelf()=0; virtual int GetOnBed()=0; - virtual void GetThicknessList(double* thickness_list)=0; - virtual void GetBedList(double* bed_list)=0; - virtual void Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type)=0; - virtual void Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type)=0; - virtual void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type)=0; - virtual void GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type)=0; - virtual double Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type)=0; + virtual void GetThicknessList(double* thickness_list)=0; + virtual void GetBedList(double* bed_list)=0; + virtual void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0; + virtual void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type)=0; + virtual void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0; + virtual void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type)=0; + virtual double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type)=0; int Enum(); Index: /issm/trunk/src/c/objects/FemModel.h =================================================================== --- /issm/trunk/src/c/objects/FemModel.h (revision 245) +++ /issm/trunk/src/c/objects/FemModel.h (revision 246) @@ -9,6 +9,8 @@ #include "../DataSet/DataSet.h" -struct FemModel { - +class DataSet; + +struct FemModel{ + DataSet* elements; DataSet* nodes; Index: /issm/trunk/src/c/objects/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Icefront.cpp (revision 245) +++ /issm/trunk/src/c/objects/Icefront.cpp (revision 246) @@ -240,5 +240,5 @@ #define __FUNCT__ "Icefront::CreateKMatrix" -void Icefront::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void Icefront::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ /*No stiffness loads applied, do nothing: */ @@ -250,5 +250,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::CreatePVector" -void Icefront::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){ +void Icefront::CreatePVector(Vec pg, void* inputs, int analysis_type){ /*Just branch to the correct element icefront vector generator, according to the type of analysis we are carrying out: */ @@ -264,5 +264,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHoriz" -void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type){ /*Branck on the type of icefront: */ @@ -278,5 +278,5 @@ #define __FUNCT__ "Icefront::CreatePVectorDiagnosticHorizSegment" -void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,ParameterInputs* inputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHorizSegment( Vec pg,void* vinputs, int analysis_type){ int i,j; @@ -306,7 +306,10 @@ Matpar* matpar=NULL; Node** element_nodes=NULL; + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; /*Recover material and fill parameters: */ - matpar=element->GetMatPar(); + matpar=(Matpar*)element->GetMatPar(); fill=element->GetShelf(); @@ -323,5 +326,5 @@ if(element_type==TriaEnum())element_nodes=(Node**)xmalloc(3*sizeof(Node*)); if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*)); - element->GetNodes(element_nodes); + element->GetNodes((void**)element_nodes); /*go through first 3 grids (all grids for tria, bed grids for penta) and identify 1st and 2nd grid: */ @@ -389,5 +392,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefont::CreatePVectorDiagnosticHorizQuad" -void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type){ +void Icefront::CreatePVectorDiagnosticHorizQuad( Vec pg,void* vinputs, int analysis_type){ int i,j; @@ -413,4 +416,5 @@ Matpar* matpar=NULL; Node** element_nodes=NULL; + ParameterInputs* inputs=NULL; /*quad grids: */ @@ -427,5 +431,6 @@ double v35[3]; double v45[3]; - + + inputs=(ParameterInputs*)vinputs; /*check icefront is associated to a pentaelem: */ @@ -434,5 +439,5 @@ } /*Recover material and fill parameters: */ - matpar=element->GetMatPar(); + matpar=(Matpar*)element->GetMatPar(); fill=element->GetShelf(); @@ -446,5 +451,5 @@ //Identify which grids are comprised in the quad: if(element_type==PentaEnum())element_nodes=(Node**)xmalloc(6*sizeof(Node*)); - element->GetNodes(element_nodes); + element->GetNodes((void**)element_nodes); grid1=-1; grid2=-1; grid3=-1; grid4=-1; @@ -549,5 +554,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::UpdateFromInputs" -void Icefront::UpdateFromInputs(ParameterInputs* inputs){ +void Icefront::UpdateFromInputs(void* vinputs){ int i; @@ -561,34 +566,19 @@ double* h_param=NULL; double* b_param=NULL; - - /*Get dof list: */ - nodes[0]->GetDofList(&doflist0[0],&numberofdofspernode); - nodes[1]->GetDofList(&doflist1[0],&numberofdofspernode); + int dofs[1]={0}; + + ParameterInputs* inputs=NULL; + + inputs=(ParameterInputs*)vinputs; + if(strcmp(type,"quad")==0){ - nodes[2]->GetDofList(&doflist2[0],&numberofdofspernode); - nodes[3]->GetDofList(&doflist3[0],&numberofdofspernode); - } - - /*Update internal data if inputs holds new values: */ - h_param=ParameterInputsRecover(inputs,"thickness"); - b_param=ParameterInputsRecover(inputs,"bed"); - - if(h_param){ - h[0]=h_param[doflist0[0]]; - h[1]=h_param[doflist1[0]]; - if(strcmp(type,"quad")==0){ - h[2]=h_param[doflist2[0]]; - h[3]=h_param[doflist3[0]]; - } - } - - if(b_param){ - b[0]=b_param[doflist0[0]]; - b[1]=b_param[doflist1[0]]; - if(strcmp(type,"quad")==0){ - b[2]=b_param[doflist2[0]]; - b[3]=b_param[doflist3[0]]; - } - } + inputs->Recover("thickness",&h[0],1,dofs,4,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,4,(void**)nodes); + } + else{ + inputs->Recover("thickness",&h[0],1,dofs,2,(void**)nodes); + inputs->Recover("bed",&h[0],1,dofs,2,(void**)nodes); + } + } @@ -699,5 +689,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::PenaltyCreateKMatrix" -void Icefront::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){ +void Icefront::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ /*do nothing: */ } @@ -705,5 +695,5 @@ #undef __FUNCT__ #define __FUNCT__ "Icefront::PenaltyCreatePVector" -void Icefront::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){ +void Icefront::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ /*do nothing: */ } Index: /issm/trunk/src/c/objects/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Icefront.h (revision 245) +++ /issm/trunk/src/c/objects/Icefront.h (revision 246) @@ -14,4 +14,5 @@ #define ICEFRONTSTRING 20 //max string length +class Element; class Icefront: public Load { @@ -56,17 +57,17 @@ void DistributeNumDofs(int* numdofspernode,int analysis_type); void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type); - void UpdateFromInputs(ParameterInputs* inputs); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); + void UpdateFromInputs(void* inputs); - void CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type); - void CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type); - void CreatePVectorDiagnosticHorizSegment( Vec pg,ParameterInputs* inputs, int analysis_type); - void CreatePVectorDiagnosticHorizQuad( Vec pg,ParameterInputs* inputs, int analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type); + void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs, int analysis_type); + void CreatePVectorDiagnosticHorizSegment( Vec pg,void* inputs, int analysis_type); + void CreatePVectorDiagnosticHorizQuad( Vec pg,void* inputs, int analysis_type); void GetDofList(int* doflist,int* pnumberofdofs); void SegmentPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal,double length,int fill); void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list, int fill); - void PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); Object* copy(); }; Index: /issm/trunk/src/c/objects/Input.cpp =================================================================== --- /issm/trunk/src/c/objects/Input.cpp (revision 246) +++ /issm/trunk/src/c/objects/Input.cpp (revision 246) @@ -0,0 +1,216 @@ +/*!\file Input.c + * \brief: implementation of the Input object + */ + + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../shared/shared.h" +#include "../include/typedefs.h" +#include "../include/macros.h" +#include "./objects.h" + + +Input::Input(){ + return; +} + +Input::~Input(){ + + xfree((void**)&vector); + return; +} + +Input::Input(char* input_name,char* input_string){ + + strcpy(name,input_name); + type=STRING; + strcpy(string,input_string); + ndof=0; + numberofnodes=0; + vector=NULL; +} + +Input::Input(char* input_name,int input_integer){ + + strcpy(name,input_name); + type=INTEGER; + integer=input_integer; + + ndof=0; + numberofnodes=0; + vector=NULL; + +} + +Input::Input(char* input_name,double input_scalar){ + + strcpy(name,input_name); + type=DOUBLE; + scalar=input_scalar; + + ndof=0; + numberofnodes=0; + vector=NULL; + +} + +Input::Input(char* input_name,double* input_vector,int input_ndof,int input_numberofnodes){ + + strcpy(name,input_name); + type=DOUBLEVEC; + vector=(double*)xmalloc(input_ndof*input_numberofnodes*sizeof(double)); + memcpy(vector,input_vector,input_ndof*input_numberofnodes*sizeof(double)); + ndof=input_ndof; + numberofnodes=input_numberofnodes; + +} + +Input::Input(char* input_name,Vec input_vector,int input_ndof, int input_numberofnodes){ + + strcpy(name,input_name); + type=DOUBLEVEC; + VecToMPISerial(&vector,input_vector); + ndof=input_ndof; + numberofnodes=input_numberofnodes; + +} + +void Input::Echo(void){ + + int i,j; + + printf("Input:\n"); + printf(" name: %s\n",name); + printf(" type: %i\n",type); + if(type==STRING) printf(" string: %s\n",string); + if(type==INTEGER) printf(" integer: %i\n",integer); + if(type==DOUBLE) printf(" double: %g\n",scalar); + if(type==DOUBLEVEC) printf(" doublevec:\n"); + for(i=0;iRecover(&values[0],ndof,dofs,3,this->nodes); //ie: recover values for dofs 0 and 2 of all nodes of the tria + + int* doflist=NULL; + int dof1; + int dof; + int i,j; + Node** nodes=NULL; + Node* node=NULL; + + /*recover pointer to nodes: */ + nodes=(Node**)vpnodes; + + if (type!=DOUBLEVEC) throw ErrorException(__FUNCT__,exprintf("%s%i%s"," cannot recover values from a ",type," input type")); + + /*Ok, we are trying to fill values. values is of size ndof*numnodes. The dofs + * for picking up the values are found in the numnodes nodes. We firt build the dof list: */ + doflist=(int*)xmalloc(ndof*numnodes*sizeof(int)); + + for(i=0;iGetDoflist1(); + + for(j=0;jRecover("x",&x[0],1,dof,1,(void**)&node); + inputs->Recover("y",&x[1],1,dof,1,(void**)&node); + inputs->Recover("z",&x[2],1,dof,1,(void**)&node); + + } @@ -494,2 +487,6 @@ } + +int Node::GetDoflist1(void){ + return doflist1[0]; +} Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 245) +++ /issm/trunk/src/c/objects/Node.h (revision 246) @@ -8,5 +8,4 @@ #include "./Object.h" #include "../toolkits/toolkits.h" -#include "./ParameterInputs.h" #define MAXDOFSPERNODE 4 @@ -67,9 +66,10 @@ void CreateVecSets(Vec pv_g,Vec pv_m,Vec pv_n,Vec pv_f,Vec pv_s); void GetDofList(int* outdoflist,int* pnumberofdofspernode); + int GetDoflist1(void); double GetX(); double GetY(); double GetZ(); Object* copy(); - void UpdateFromInputs(ParameterInputs* inputs); + void UpdateFromInputs(void* inputs); void Configure(void* pnodes); Node* GetUpperNode(); Index: /issm/trunk/src/c/objects/OptArgs.h =================================================================== --- /issm/trunk/src/c/objects/OptArgs.h (revision 245) +++ /issm/trunk/src/c/objects/OptArgs.h (revision 246) @@ -24,4 +24,5 @@ #include "./FemModel.h" #include "./ParameterInputs.h" +class ParameterInputs; struct OptArgs{ Index: /issm/trunk/src/c/objects/ParameterInputs.cpp =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.cpp (revision 245) +++ /issm/trunk/src/c/objects/ParameterInputs.cpp (revision 246) @@ -1,380 +1,241 @@ -/*!\file: ParameterInputs.cpp - * \brief routines for loading from matlab the parameter inputs to our non-linear iterations. - */ - - -#include "./ParameterInputs.h" +/*!\file ParameterInputs.c + * \brief: implementation of the ParameterInputs object + */ + + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include +#include "../EnumDefinitions/EnumDefinitions.h" #include "../shared/shared.h" +#include "../include/typedefs.h" #include "../include/macros.h" - -/*-------------------------------------------------- - NewParameterInputs - --------------------------------------------------*/ - -ParameterInputs* NewParameterInputs(void) { - - ParameterInputs* inputs=NULL; - - inputs=(ParameterInputs*)xmalloc(sizeof(ParameterInputs)); - - inputs->nfields=0; - inputs->fields=NULL; - inputs->field_names=NULL; - inputs->fields_length=NULL; - - return inputs; -} - - -/*-------------------------------------------------- - ParameterInputsEcho - --------------------------------------------------*/ -int ParameterInputsEcho(ParameterInputs* inputs,int longecho){ - - int noerr=1; - int i,j; - - /*short echo: */ - if(!longecho){ - _printf_("field names, field pointers: \n"); - for (i=0;infields;i++){ - _printf_("%s: %p\n",inputs->field_names[i],inputs->fields[i]); - } - } - /*long echo: */ - else{ - for (i=0;infields;i++){ - double* values=inputs->fields[i]; - _printf_("field names: %s\n",inputs->field_names[i]); - _printf_("field values: \n"); - for(j=0;jfields_length[i];j++){ - printf(" %g\n",values[j]); - } - } - } - return noerr; -} - -/*-------------------------------------------------- - ParameterInputsRecover - --------------------------------------------------*/ -double* ParameterInputsRecover(ParameterInputs* inputs,char* field_name){ - - double* param=NULL; - int i; - - /*Find field_name input field, and return the value: */ - for(i=0;infields;i++){ - if (strcmp(inputs->field_names[i],field_name)==0){ - param=inputs->fields[i]; - } - } - return param; -} - -/*-------------------------------------------------- - DeleteParameterInputs - --------------------------------------------------*/ -void DeleteParameterInputs(ParameterInputs** pinputs){ - - ParameterInputs* inputs=NULL; - int i; - char* field_name=NULL; - double* field=NULL; - - #ifdef _PARALLEL_ //inputs operation is only done on the cluster, on the serial side, we leave the matlab memory manager handle inputs. - inputs=*pinputs; - - if(inputs){ - - for(i=0;infields;i++){ - field_name=inputs->field_names[i]; - xfree((void**)&field_name); - - /*In parallel, erase: */ - field=inputs->fields[i]; - xfree((void**)&field); - } - - xfree((void**)&inputs->field_names); - xfree((void**)&inputs->fields); - xfree((void**)&inputs->fields_length); - } - xfree((void**)pinputs); - #else - #endif - return; -} - -/*-------------------------------------------------- - ParameterInputsAddFromVec - --------------------------------------------------*/ -int ParameterInputsAddFromVec(ParameterInputs* inputs,Vec vector,char* pfield_name){ - - int noerr=1; - double* double_vector=NULL; - double* double_vector2=NULL; - int vector_size; - char* field_name; - char** field_names; - int field_name_size; - double** fields=NULL; - int* fields_length; - int i; - - /*Take vector, serialize it, and plug it into inputs with field name "field_name": */ - if (vector){ - VecGetSize(vector,&vector_size); - VecToMPISerial(&double_vector,vector); +#include "./objects.h" + +ParameterInputs::ParameterInputs(){ + + dataset=new DataSet(); + +} - if(inputs->nfields==0){ - /*Ok, we need to allocate fields, field_names and fields_length: */ - inputs->nfields=1; - - inputs->field_names=(char**)xmalloc(sizeof(char*)); - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[0]=field_name; - - inputs->fields=(double**)xmalloc(sizeof(double*)); - inputs->fields[0]=double_vector; - - inputs->fields_length=(int*)xmalloc(sizeof(int)); - inputs->fields_length[0]=vector_size; - } - else{ - /*First, figure out if there is already an existing field with name pfield_name: */ - for(i=0;infields;i++){ - if (strcmp(inputs->field_names[i],pfield_name)==0){ - /*Erase old field and plug new field: */ - double_vector2=inputs->fields[i]; - xfree((void**)&double_vector2); //no leak! - inputs->fields[i]=double_vector; - inputs->fields_length[i]=vector_size; - return noerr; - } - } - - /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate - * the inputs: */ - field_names=inputs->field_names; //save field_names pointer. - inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*)); - for(i=0;infields;i++){ - inputs->field_names[i]=field_names[i]; - } - xfree((void**)&field_names); - - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[inputs->nfields]=field_name; - - - - fields=inputs->fields; //save fields pointer - inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*)); - for(i=0;infields;i++){ - inputs->fields[i]=fields[i]; - } - xfree((void**)&fields); - inputs->fields[inputs->nfields]=double_vector; - - fields_length=inputs->fields_length; - inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int)); - for(i=0;infields;i++){ - inputs->fields_length[i]=fields_length[i]; - } - inputs->fields_length[inputs->nfields]=vector_size; - xfree((void**)&fields_length); - - inputs->nfields++; - } - } - else{ - /*Do nothing, there is no input to add!: */ - ; - } - return noerr; -} - -/*-------------------------------------------------- - ParameterInputsAddFromDouble - --------------------------------------------------*/ -int ParameterInputsAddFromDouble(ParameterInputs* inputs,double scalar, char* pfield_name){ - - int noerr=1; - double* double_vector2=NULL; - double* double_vector=NULL; - int vector_size; - char* field_name; - char** field_names; - int field_name_size; - double** fields=NULL; - int* fields_length; - int i; - - /*Create double_vector pointing to scalar: */ - vector_size=1; - double_vector=(double*)xmalloc(sizeof(double)); - double_vector[0]=scalar; - - - if(inputs->nfields==0){ - /*Ok, we need to allocate fields, field_names and fields_length: */ - inputs->nfields=1; +ParameterInputs::~ParameterInputs(){ + + delete dataset; + +} + +void ParameterInputs::Echo(){ + + printf("ParameterInputs echo: \n"); + dataset->Echo(); + +} + +void ParameterInputs::purge(char* name){ + + int i; + Input* input=NULL; + + /*Go through dataset, and figure out if an Input + * already has the name "name". If so, delete it : */ + + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + + if (input->IsSameName(name)){ + /*delete object: */ + dataset->DeleteObject(input); + } + } +} - inputs->field_names=(char**)xmalloc(sizeof(char*)); - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[0]=field_name; - - inputs->fields=(double**)xmalloc(sizeof(double*)); - inputs->fields[0]=double_vector; - - inputs->fields_length=(int*)xmalloc(sizeof(int)); - inputs->fields_length[0]=vector_size; - } - else{ - /*First, figure out if there is already an existing field with name pfield_name: */ - for(i=0;infields;i++){ - if (strcmp(inputs->field_names[i],pfield_name)==0){ - /*Erase old field and plug new field: */ - double_vector2=inputs->fields[i]; - xfree((void**)&double_vector2); //no leak! - inputs->fields[i]=double_vector; - inputs->fields_length[i]=vector_size; - return noerr; - } - } - - /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate - * the inputs: */ - field_names=inputs->field_names; //save field_names pointer. - inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*)); - for(i=0;infields;i++){ - inputs->field_names[i]=field_names[i]; - } - xfree((void**)&field_names); - - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[inputs->nfields]=field_name; - - - - fields=inputs->fields; //save fields pointer - inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*)); - for(i=0;infields;i++){ - inputs->fields[i]=fields[i]; - } - xfree((void**)&fields); - inputs->fields[inputs->nfields]=double_vector; - - fields_length=inputs->fields_length; - inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int)); - for(i=0;infields;i++){ - inputs->fields_length[i]=fields_length[i]; - } - inputs->fields_length[inputs->nfields]=vector_size; - xfree((void**)&fields_length); - - inputs->nfields++; - } - return noerr; -} - -/*-------------------------------------------------- - ParameterInputsAddFromMat - --------------------------------------------------*/ -int ParameterInputsAddFromMat(ParameterInputs* inputs,double* vector,int vector_size,char* pfield_name){ - - int noerr=1; - char* field_name; - char** field_names; - int field_name_size; - double** fields=NULL; - int* fields_length; - double* double_vector=NULL; - double* double_vector2=NULL; - int i; - - /*Make a copy of input vector, so as to not create any leaks: */ - double_vector=(double*)xmalloc(vector_size*sizeof(double)); - memcpy(double_vector,vector,vector_size*sizeof(double)); - - /*Take vector, serialize it, and plug it into inputs with field name "field_name": */ - if (vector){ - - if(inputs->nfields==0){ - /*Ok, we need to allocate fields, field_names and fields_length: */ - inputs->nfields=1; - - inputs->field_names=(char**)xmalloc(sizeof(char*)); - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[0]=field_name; - - inputs->fields=(double**)xmalloc(sizeof(double*)); - inputs->fields[0]=double_vector; - - inputs->fields_length=(int*)xmalloc(sizeof(int)); - inputs->fields_length[0]=vector_size; - } - else{ - - /*First, figure out if there is already an existing field with name pfield_name: */ - for(i=0;infields;i++){ - if (strcmp(inputs->field_names[i],pfield_name)==0){ - /*Just xfree old field and plug new field: */ - double_vector2=inputs->fields[i]; - xfree((void**)&double_vector2); - inputs->fields[i]=double_vector; - inputs->fields_length[i]=vector_size; - return noerr; - } - } - - /*Ok, if we are here, the field did not already exist, we need to add it. First reallocate - * the inputs: */ - field_names=inputs->field_names; //save field_names pointer. - inputs->field_names=(char**)xmalloc((inputs->nfields+1)*sizeof(char*)); - for(i=0;infields;i++){ - inputs->field_names[i]=field_names[i]; - } - xfree((void**)&field_names); - - field_name_size=strlen(pfield_name)+1; - field_name=(char*)xmalloc(field_name_size*sizeof(char)); - strcpy(field_name,pfield_name); - inputs->field_names[inputs->nfields]=field_name; - - - fields=inputs->fields; //save fields pointer - inputs->fields=(double**)xmalloc((inputs->nfields+1)*sizeof(double*)); - for(i=0;infields;i++){ - inputs->fields[i]=fields[i]; - } - xfree((void**)&fields); - inputs->fields[inputs->nfields]=double_vector; - - fields_length=inputs->fields_length; - inputs->fields_length=(int*)xmalloc((inputs->nfields+1)*sizeof(int)); - for(i=0;infields;i++){ - inputs->fields_length[i]=fields_length[i]; - } - inputs->fields_length[inputs->nfields]=vector_size; - xfree((void**)&fields_length); - - inputs->nfields++; - } - } - else{ - /*Do nothing, there is no input to add!: */ - ; - } - return noerr; -} + +void ParameterInputs::Add(char* name,char* string){ + + Input* input=NULL; + + /*First, purge object with same name: */ + this->purge(name); + + /*We are sure an input of the same name does not exist. Create new + * input: */ + input=new Input(name,string); + + /*Add input to dataset: */ + dataset->AddObject(input); + +} + +void ParameterInputs::Add(char* name,int integer){ + + Input* input=NULL; + + /*First, purge object with same name: */ + this->purge(name); + + /*We are sure an input of the same name does not exist. Create new + * input: */ + input=new Input(name,integer); + + /*Add input to dataset: */ + dataset->AddObject(input); +} + + +void ParameterInputs::Add(char* name,double scalar){ + + Input* input=NULL; + + /*First, purge object with same name: */ + this->purge(name); + + /*We are sure an input of the same name does not exist. Create new + * input: */ + input=new Input(name,scalar); + + /*Add input to dataset: */ + dataset->AddObject(input); + +} + +void ParameterInputs::Add(char* name,double* vector,int ndof,int numberofnodes){ + + Input* input=NULL; + + /*First, purge object with same name: */ + this->purge(name); + + /*We are sure an input of the same name does not exist. Create new + * input: */ + input=new Input(name,vector,ndof,numberofnodes); + + /*Add input to dataset: */ + dataset->AddObject(input); + +} + +void ParameterInputs::Add(char* name,Vec petscvector,int ndof, int numberofnodes){ + + Input* input=NULL; + + /*First, purge object with same name: */ + this->purge(name); + + /*We are sure an input of the same name does not exist. Create new + * input: */ + input=new Input(name,petscvector,ndof,numberofnodes); + + /*Add input to dataset: */ + dataset->AddObject(input); + +} +int ParameterInputs::Recover(char* name, char** pstring){ + + int found=0; + int i; + Input* input=NULL; + char* string=NULL; + + /*Go through dataset, and figure out if an Input + * has the name "name": */ + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + + if (input->IsSameName(name)){ + + /*Get string out of this input: */ + input->Recover(&string); + found=1;break; + } + } + /*Assign output pointers:*/ + *pstring=string; + return found; + +} +int ParameterInputs::Recover(char* name, int* pinteger){ + + int found=0; + int i; + Input* input=NULL; + int integer; + + /*Go through dataset, and figure out if an Input + * has the name "name": */ + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + + if (input->IsSameName(name)){ + + /*Get string out of this input: */ + input->Recover(&integer); + found=1; break; + + } + } + /*Assign output pointers:*/ + *pinteger=integer; + + return found; +} +int ParameterInputs::Recover(char* name, double* pscalar){ + + int found=0; + int i; + Input* input=NULL; + double scalar; + + /*Go through dataset, and figure out if an Input + * has the name "name": */ + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + + if (input->IsSameName(name)){ + + /*Get string out of this input: */ + input->Recover(&scalar); + found=1; break; + + } + } + /*Assign output pointers:*/ + *pscalar=scalar; + + return found; + +} + +#undef __FUNCT__ +#define __FUNCT__ "ParameterInputs::Recover" +int ParameterInputs::Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** pnodes){ + + int i; + Node** nodes=NULL; + Input* input=NULL; + int found=0; + + /*First recover nodes pointer :*/ + nodes=(Node**)pnodes; + + /*Go through dataset, and figure out if an Input + * has the name "name": */ + for(i=0;iSize();i++){ + input=(Input*)dataset->GetObjectByOffset(i); + + if (input->IsSameName(name)){ + found=1; + break; + } + } + + if(found==0)return 0; + + /*call submethod: */ + input->Recover(values,ndof,dofs,numnodes,(void**)nodes); + return 1; +} Index: /issm/trunk/src/c/objects/ParameterInputs.h =================================================================== --- /issm/trunk/src/c/objects/ParameterInputs.h (revision 245) +++ /issm/trunk/src/c/objects/ParameterInputs.h (revision 246) @@ -1,32 +1,40 @@ /*!\file: ParameterInputs.h - * \brief prototype for parameter inputs, from matlab workspace. + * \brief dataset filled with input objects. Used by CreateKMatrix and CreatePVector to hold + * parameters coming from the workspace. */ -#ifndef PARAMETER_INPUTS_H_ -#define PARAMETER_INPUTS_H_ + +#ifndef _PARAMETERINPUTS_H_ +#define _PARAMETERINPUTS_H_ #include "../toolkits/toolkits.h" +#include "../DataSet/DataSet.h" -struct ParameterInputs{ +class DataSet; - int nfields; //number of inputs - char** field_names; //characterization of each field. - double** fields; //value of each field. - int* fields_length; //length of each values field +class ParameterInputs{ + private: + DataSet* dataset; //cannot declare DataSet* dataset, to avoid circular referencing. + + public: + ParameterInputs(); + ~ParameterInputs(); + void Echo(); + + void Add(char* name,char* string); + void Add(char* name,int integer); + void Add(char* name,double scalar); + void Add(char* name,double* vector,int ndof,int numberofnodes); + void Add(char* name,Vec petscvector,int ndof, int numberofnodes); + + int Recover(char* name, char** pstring); + int Recover(char* name, int* pinteger); + int Recover(char* name, double* pscalar); + int Recover(char* name,double* values, int ndof, int* dofs,int numnodes,void** nodes); //void** because otherwise we get a circular dependency with Node + + /*declaration found in io: */ + void Init( void* data_handle); + void purge(char* name); }; -ParameterInputs* NewParameterInputs(); -void DeleteParameterInputs(ParameterInputs** pinputs); -int ParameterInputsEcho(ParameterInputs* inputs,int longecho); -double* ParameterInputsRecover(ParameterInputs* inputs,char* field_name); - -int ParameterInputsAddFromVec(ParameterInputs* inputs,Vec vector,char* field_name); -int ParameterInputsAddFromMat(ParameterInputs* inputs,double* vector,int vector_size,char* field_name); -int ParameterInputsAddFromDouble(ParameterInputs* inputs,double scalar, char* pfield_name); - -#ifdef _SERIAL_ -#include "mex.h" -void ParameterInputsInit(ParameterInputs** pthis, const mxArray* inputhandle); #endif - -#endif Index: /issm/trunk/src/c/objects/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Penpair.cpp (revision 245) +++ /issm/trunk/src/c/objects/Penpair.cpp (revision 246) @@ -208,5 +208,5 @@ #define __FUNCT__ "Penpair::CreateKMatrix" -void Penpair::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void Penpair::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ /*No loads applied, do nothing: */ @@ -217,5 +217,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::CreatePVector" -void Penpair::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){ +void Penpair::CreatePVector(Vec pg, void* inputs, int analysis_type){ /*No loads applied, do nothing: */ @@ -225,5 +225,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::UpdateFromInputs" -void Penpair::UpdateFromInputs(ParameterInputs* inputs){ +void Penpair::UpdateFromInputs(void* inputs){ } @@ -231,5 +231,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::PenaltyCreateKMatrix" -void Penpair::PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type){ +void Penpair::PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type){ if(numdofs==1){ @@ -245,5 +245,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penpair::PenaltyCreatePVector" -void Penpair::PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type){ +void Penpair::PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type){ if(numdofs==1){ Index: /issm/trunk/src/c/objects/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Penpair.h (revision 245) +++ /issm/trunk/src/c/objects/Penpair.h (revision 246) @@ -9,4 +9,5 @@ #include "./Element.h" +class Element; class Penpair: public Load{ @@ -52,9 +53,9 @@ void DistributeNumDofs(int* numdofspernode,int analysis_type); void Configure(void* elements,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type); - void CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type); - void UpdateFromInputs(ParameterInputs* inputs); - void PenaltyCreateKMatrix(Mat Kgg,ParameterInputs* inputs,double kmax,int analysis_type); - void PenaltyCreatePVector(Vec pg,ParameterInputs* inputs,double kmax,int analysis_type); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type); + void UpdateFromInputs(void* inputs); + void PenaltyCreateKMatrix(Mat Kgg,void* inputs,double kmax,int analysis_type); + void PenaltyCreatePVector(Vec pg,void* inputs,double kmax,int analysis_type); Object* copy(); }; Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 245) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 246) @@ -279,5 +279,5 @@ #define __FUNCT__ "Penta::CreateKMatrix" -void Penta::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void Penta::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -301,5 +301,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz" -void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type){ +void Penta::CreateKMatrixDiagnosticHoriz( Mat Kgg, void* vinputs, int analysis_type){ @@ -368,8 +368,7 @@ /*input parameters for structural analysis (diagnostic): */ - double* velocity_param=NULL; - double vxvy_list[numgrids][2]; - double* oldvelocity_param=NULL; - double oldvxvy_list[numgrids][2]; + 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; @@ -381,6 +380,11 @@ double MOUNTAINKEXPONENT=10; + ParameterInputs* inputs=NULL; + /*Collapsed formulation: */ Tria* tria=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /*Figure out if this pentaelem is collapsed. If so, then bailout, except if it is at the @@ -397,5 +401,5 @@ /*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=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->CreateKMatrix(Kgg,inputs, analysis_type); delete tria; @@ -407,6 +411,6 @@ /*recover extra inputs from users, at current convergence iteration: */ - velocity_param=ParameterInputsRecover(inputs,"velocity"); - oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity"); + 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: */ @@ -418,28 +422,4 @@ for(j=0;jCreateKMatrixDiagnosticHorizFriction(Kgg,inputs,analysis_type); delete tria; @@ -574,5 +554,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert" -void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type){ +void Penta::CreateKMatrixDiagnosticVert( Mat Kgg, void* vinputs, int analysis_type){ /* local declarations */ @@ -610,12 +590,18 @@ double Bprime[NDOF1][numgrids]; double DL_scalar; - + + ParameterInputs* inputs=NULL; + /*Collapsed formulation: */ Tria* tria=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + /*If this element is on the surface, we have a dynamic boundary condition that applies, as a stiffness * matrix: */ if(onsurface){ - tria=SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface + tria=(Tria*)SpawnTria(3,4,5); //nodes 3,4 and 5 are on the surface tria->CreateKMatrixDiagnosticSurfaceVert(Kgg,inputs, analysis_type); delete tria; @@ -706,5 +692,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVector" -void Penta::CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type){ +void Penta::CreatePVector(Vec pg, void* inputs, int analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -724,60 +710,34 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::UpdateFromInputs" -void Penta::UpdateFromInputs(ParameterInputs* inputs){ - - int i; - int doflist[MAXDOFSPERNODE*3]; - int numberofdofspernode; - - /*element: */ - double* h_param=NULL; - double* s_param=NULL; - double* b_param=NULL; - double* k_param=NULL; - double* melting_param=NULL; - double* accumulation_param=NULL; - - /*material: */ - double* temperature_average_param=NULL; +void Penta::UpdateFromInputs(void* vinputs){ + + int dofs[1]={0}; + double temperature_list[6]; double temperature_average; - double* B_param=NULL; + double B_list[6]; double B_average; - - /*Get dof list: */ - GetDofList(&doflist[0],&numberofdofspernode); + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /*Update internal data if inputs holds new values: */ - h_param=ParameterInputsRecover(inputs,"thickness"); - s_param=ParameterInputsRecover(inputs,"surface"); - b_param=ParameterInputsRecover(inputs,"bed"); - k_param=ParameterInputsRecover(inputs,"drag"); - melting_param=ParameterInputsRecover(inputs,"melting"); - accumulation_param=ParameterInputsRecover(inputs,"accumulation"); - - for(i=0;i<6;i++){ - if(h_param)h[i]=h_param[doflist[i*numberofdofspernode+0]]; - if(s_param)s[i]=s_param[doflist[i*numberofdofspernode+0]]; - if(b_param)b[i]=b_param[doflist[i*numberofdofspernode+0]]; - if(k_param)k[i]=k_param[doflist[i*numberofdofspernode+0]]; - if(melting_param)melting[i]=melting_param[doflist[i*numberofdofspernode+0]]; - if(accumulation_param)accumulation[i]=accumulation_param[doflist[i*numberofdofspernode+0]]; - } - + inputs->Recover("thickness",&h[0],1,dofs,6,(void**)nodes); + inputs->Recover("surface",&s[0],1,dofs,6,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,6,(void**)nodes); + inputs->Recover("drag",&k[0],1,dofs,6,(void**)nodes); + inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes); + inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes); + //Update material if necessary - temperature_average_param=ParameterInputsRecover(inputs,"temperature_average"); - B_param=ParameterInputsRecover(inputs,"B"); - if(temperature_average_param){ - temperature_average=0; - for(i=0;i<6;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]]; - temperature_average=temperature_average/6.0; + if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){ + temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2]+temperature_list[3]+temperature_list[4]+temperature_list[5])/6.0; B_average=Paterson(temperature_average); matice->SetB(B_average); } - - if(B_param){ - B_average=0; - for(i=0;i<6;i++) B_average+=B_param[doflist[i*numberofdofspernode+0]]; - B_average=B_average/6.0; + + if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){ + B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0; matice->SetB(B_average); } @@ -785,5 +745,5 @@ } -Matpar* Penta::GetMatPar(){ +void* Penta::GetMatPar(){ return matpar; } @@ -794,6 +754,8 @@ -void Penta::GetNodes(Node** pnodes){ +void Penta::GetNodes(void** vpnodes){ int i; + Node** pnodes=(Node**)vpnodes; + for(i=0;i<6;i++){ pnodes[i]=nodes[i]; @@ -823,5 +785,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Du" -void Penta::Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type){ +void Penta::Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type){ int i; @@ -834,5 +796,5 @@ else{ - tria=SpawnTria(3,4,5); //grids 0, 1 and 2 make the new tria. + tria=(Tria*)SpawnTria(3,4,5); //grids 0, 1 and 2 make the new tria. tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type); delete tria; @@ -843,5 +805,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Gradj" -void Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){ +void Penta::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){ if (strcmp(control_type,"drag")==0){ @@ -856,5 +818,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjDrag" -void Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type){ +void Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){ @@ -867,5 +829,5 @@ else{ - tria=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type); delete tria; @@ -876,5 +838,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::GradjB" -void Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type){ +void Penta::GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type){ throw ErrorException(__FUNCT__," not supported yet!"); } @@ -882,5 +844,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::Misfit" -double Penta::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){ +double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type){ double J; @@ -894,5 +856,5 @@ /*use the Tria::CreateMisfit routine, on a Tria Element which is made of the 3 upper grids at the surface of the ice sheet. : */ - tria=SpawnTria(3,4,5); + tria=(Tria*)SpawnTria(3,4,5); J=tria->Misfit( velocity,obs_velocity,inputs,analysis_type); delete tria; @@ -903,5 +865,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::SpawnTria" -Tria* Penta::SpawnTria(int g0, int g1, int g2){ +void* Penta::SpawnTria(int g0, int g1, int g2){ /*out of grids g0,g1 and g2 from Penta, build a tria element: */ @@ -1345,5 +1307,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz" -void Penta::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type){ +void Penta::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs,int analysis_type){ int i,j; @@ -1388,6 +1350,11 @@ double pe_g_gaussian[numdofs]; + ParameterInputs* inputs=NULL; + /*Spawning: */ Tria* tria=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; @@ -1405,5 +1372,5 @@ /*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=SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. tria->CreatePVector(pg,inputs, analysis_type); delete tria; @@ -1644,5 +1611,5 @@ #undef __FUNCT__ #define __FUNCT__ "Penta:CreatePVectorDiagnosticVert" -void Penta::CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type){ +void Penta::CreatePVectorDiagnosticVert( Vec pg, void* vinputs,int analysis_type){ int i; @@ -1685,15 +1652,21 @@ /*input parameters for structural analysis (diagnostic): */ - double* velocity_param=NULL; - double vx_list[numgrids]; - double vy_list[numgrids]; + double vx_list[numgrids]={0,0,0}; + double vy_list[numgrids]={0,0,0}; double du[3]; double dv[3]; double dudx,dvdy; + int dofs1[1]={0}; + int dofs2[1]={1}; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the *diagnostic base vertical stifness: */ if(onbed){ - tria=SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock + tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 are on the bedrock tria->CreatePVectorDiagnosticBaseVert(pg,inputs, analysis_type); delete tria; @@ -1709,15 +1682,6 @@ /* recover input parameters: */ - velocity_param=ParameterInputsRecover(inputs,"velocity"); - - if(!velocity_param)throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!"); - - /*Initialize vx_list and vy_list: */ - for(i=0;iRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!"); + inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); /*Get gaussian points and weights :*/ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 245) +++ /issm/trunk/src/c/objects/Penta.h (revision 246) @@ -12,4 +12,5 @@ #include "./Node.h" #include "./Tria.h" +#include "./ParameterInputs.h" class Penta: public Element{ @@ -71,24 +72,24 @@ int MyRank(); void Configure(void* loads,void* nodes,void* materials); - void CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type); - void CreateKMatrixDiagnosticHoriz( Mat Kgg, ParameterInputs* inputs, int analysis_type); - void CreateKMatrixDiagnosticVert( Mat Kgg, ParameterInputs* inputs, int analysis_type); - void CreatePVector(Vec pg, ParameterInputs* inputs, int analysis_type); - void UpdateFromInputs(ParameterInputs* inputs); + void CreateKMatrix(Mat Kgg,void* inputs,int analysis_type); + void CreateKMatrixDiagnosticHoriz( Mat Kgg, void* inputs, int analysis_type); + void CreateKMatrixDiagnosticVert( Mat Kgg, void* inputs, int analysis_type); + void CreatePVector(Vec pg, void* inputs, int analysis_type); + void UpdateFromInputs(void* inputs); void GetDofList(int* doflist,int* pnumberofdofs); - Matpar* GetMatPar(); + void* GetMatPar(); int GetShelf(); - void GetNodes(Node** nodes); + void GetNodes(void** nodes); int GetOnBed(); - void Du(Vec du_g,double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type); - void Gradj(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type,char* control_type); - void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type); - void GradjB(Vec grad_g,double* u_g,double* lambda_g,ParameterInputs* inputs,int analysis_type); - double Misfit(double* u_g,double* u_g_obs,ParameterInputs* inputs,int analysis_type); + void Du(Vec du_g,double* u_g,double* u_g_obs,void* inputs,int analysis_type); + void Gradj(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,char* control_type); + void GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type); + void GradjB(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type); + double Misfit(double* u_g,double* u_g_obs,void* inputs,int analysis_type); void GetThicknessList(double* thickness_list); void GetBedList(double* bed_list); Object* copy(); - Tria* SpawnTria(int g0, int g1, int g2); + void* SpawnTria(int g0, int g1, int g2); void GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3l4); @@ -102,6 +103,6 @@ void GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4); void GetJacobianInvert(double* Jinv, double* xyz_list,double* gauss_l1l2l3l4); - void CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs,int analysis_type); - void CreatePVectorDiagnosticVert( Vec pg, ParameterInputs* inputs,int analysis_type); + void CreatePVectorDiagnosticHoriz( Vec pg, void* inputs,int analysis_type); + void CreatePVectorDiagnosticVert( Vec pg, void* inputs,int analysis_type); void GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4); void GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 245) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 246) @@ -266,5 +266,5 @@ #define __FUNCT__ "Tria::CreateKMatrix" -void Tria::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void Tria::CreateKMatrix(Mat Kgg,void* inputs,int analysis_type){ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ @@ -284,5 +284,5 @@ #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz" -void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type){ +void Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type){ @@ -328,13 +328,17 @@ /*input parameters for structural analysis (diagnostic): */ - double* velocity_param=NULL; - double* oldvelocity_param=NULL; - double vxvy_list[numgrids][2]; - double oldvxvy_list[numgrids][2]; + double vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; + double oldvxvy_list[numgrids][2]={{0,0},{0,0},{0,0}}; double thickness; + int dofs[2]={0,1}; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /*recover extra inputs from users, at current convergence iteration: */ - velocity_param=ParameterInputsRecover(inputs,"velocity"); - oldvelocity_param=ParameterInputsRecover(inputs,"old_velocity"); + 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: */ @@ -344,25 +348,4 @@ /* Set Ke_gg to 0: */ for(i=0;iRecover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); /*Build alpha2_list used by drag stiffness matrix*/ @@ -691,5 +667,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::CreatePVector" -void Tria::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){ +void Tria::CreatePVector(Vec pg,void* inputs,int analysis_type){ /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */ @@ -708,5 +684,5 @@ #define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz" -void Tria::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){ +void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type){ int i,j; @@ -746,4 +722,9 @@ /*input parameters for structural analysis (diagnostic): */ double thickness; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /* Get node coordinates and dof list: */ @@ -871,60 +852,34 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::UpdateFromInputs" -void Tria::UpdateFromInputs(ParameterInputs* inputs){ - - int i; - int doflist[MAXDOFSPERNODE*3]; - int numberofdofspernode; - - /*element: */ - double* h_param=NULL; - double* s_param=NULL; - double* b_param=NULL; - double* k_param=NULL; - double* melting_param=NULL; - double* accumulation_param=NULL; - - /*material: */ - double* temperature_average_param=NULL; +void Tria::UpdateFromInputs(void* vinputs){ + + int dofs[1]={0}; + double temperature_list[3]; double temperature_average; - double* B_param=NULL; + double B_list[3]; double B_average; - - /*Get dof list: */ - GetDofList(&doflist[0],&numberofdofspernode); + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /*Update internal data if inputs holds new values: */ - h_param=ParameterInputsRecover(inputs,"thickness"); - s_param=ParameterInputsRecover(inputs,"surface"); - b_param=ParameterInputsRecover(inputs,"bed"); - k_param=ParameterInputsRecover(inputs,"drag"); - melting_param=ParameterInputsRecover(inputs,"melting"); - accumulation_param=ParameterInputsRecover(inputs,"accumulation"); - - for(i=0;i<3;i++){ - if(h_param)h[i]=h_param[doflist[i*numberofdofspernode+0]]; - if(s_param)s[i]=s_param[doflist[i*numberofdofspernode+0]]; - if(b_param)b[i]=b_param[doflist[i*numberofdofspernode+0]]; - if(k_param)k[i]=k_param[doflist[i*numberofdofspernode+0]]; - if(melting_param)melting[i]=melting_param[doflist[i*numberofdofspernode+0]]; - if(accumulation_param)accumulation[i]=accumulation_param[doflist[i*numberofdofspernode+0]]; - } - + inputs->Recover("thickness",&h[0],1,dofs,3,(void**)nodes); + inputs->Recover("surface",&s[0],1,dofs,3,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,3,(void**)nodes); + inputs->Recover("drag",&k[0],1,dofs,3,(void**)nodes); + inputs->Recover("melting",&melting[0],1,dofs,3,(void**)nodes); + inputs->Recover("accumulation_param",&accumulation[0],1,dofs,3,(void**)nodes); + //Update material if necessary - temperature_average_param=ParameterInputsRecover(inputs,"temperature_average"); - B_param=ParameterInputsRecover(inputs,"B"); - if(temperature_average_param){ - temperature_average-0; - for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]]; - temperature_average=temperature_average/3.0; + if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,3,(void**)nodes)){ + temperature_average=(temperature_list[0]+temperature_list[1]+temperature_list[2])/3.0; B_average=Paterson(temperature_average); matice->SetB(B_average); } - - if(B_param){ - B_average=0; - for(i=0;i<3;i++) B_average+=B_param[doflist[i*numberofdofspernode+0]]; - B_average=B_average/3.0; + + if(inputs->Recover("B",&B_list[0],1,dofs,3,(void**)nodes)){ + B_average=(B_list[0]+B_list[1]+B_list[2])/3.0; matice->SetB(B_average); } @@ -1311,5 +1266,5 @@ } -Matpar* Tria::GetMatPar(){ +void* Tria::GetMatPar(){ return matpar; } @@ -1320,6 +1275,8 @@ -void Tria::GetNodes(Node** pnodes){ +void Tria::GetNodes(void** vpnodes){ int i; + Node** pnodes=(Node**)vpnodes; + for(i=0;i<3;i++){ pnodes[i]=nodes[i]; @@ -1346,5 +1303,7 @@ Object* Tria::copy() { + return new Tria(*this); + } @@ -1352,5 +1311,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::Du" -void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){ +void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,void* vinputs,int analysis_type){ int i; @@ -1404,6 +1363,10 @@ double scaley=0; double scale=0; - double* fit_param=NULL; double fit=-1; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /* Get node coordinates and dof list: */ @@ -1415,9 +1378,5 @@ /* Recover input data: */ - fit_param=ParameterInputsRecover(inputs,"fit"); - if(fit_param){ - fit=*fit_param; - } - else throw ErrorException(__FUNCT__," missing fit input parameter"); + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); /*Initialize velocities: */ @@ -1570,5 +1529,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::Gradj" -void Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){ +void Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,void* inputs,int analysis_type,char* control_type){ if (strcmp(control_type,"drag")==0){ @@ -1583,5 +1542,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::GradjDrag" -void Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){ +void Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){ @@ -1603,12 +1562,6 @@ double adjy_list[numgrids]; - double* basal_drag_param=NULL; - double K_list[numgrids]; double drag; - - double* bed_param=NULL; - double bed_list[numgrids]; - double* thickness_param=NULL; - double thickness_list[numgrids]; + int dofs[1]={0}; /* gaussian points: */ @@ -1650,4 +1603,9 @@ double epsilon[3]; /* epsilon=[exx,eyy,exy];*/ + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + /* Get node coordinates and dof list: */ GetElementNodeData( &xyz_list[0][0], nodes, numgrids); @@ -1658,7 +1616,7 @@ /* recover input parameters: */ - basal_drag_param=ParameterInputsRecover(inputs,"drag"); - bed_param=ParameterInputsRecover(inputs,"bed"); - thickness_param=ParameterInputsRecover(inputs,"thickness"); + inputs->Recover("drag",&k[0],1,dofs,numgrids,(void**)nodes); + inputs->Recover("bed",&b[0],1,dofs,numgrids,(void**)nodes); + inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes); /*Initialize parameter lists: */ @@ -1670,22 +1628,4 @@ adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]]; adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]]; - if(basal_drag_param){ - K_list[i]=basal_drag_param[doflist[i*numberofdofspernode+0]]; - } - else{ - K_list[i]=k[i]; - } - if(bed_param){ - bed_list[i]=bed_param[doflist[i*numberofdofspernode+0]]; - } - else{ - bed_list[i]=b[i]; - } - if(thickness_param){ - thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]]; - } - else{ - thickness_list[i]=h[i]; - } } @@ -1723,7 +1663,7 @@ friction->rho_ice=matpar->GetRhoIce(); friction->rho_water=matpar->GetRhoWater(); - friction->K=&K_list[0]; - friction->bed=&bed_list[0]; - friction->thickness=&thickness_list[0]; + friction->K=&k[0]; + friction->bed=&b[0]; + friction->thickness=&h[0]; friction->velocities=&vxvy_list[0][0]; friction->p=p; @@ -1746,5 +1686,5 @@ /*Recover alpha_complement and k: */ GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3); - GetParameterValue(&drag, &K_list[0],gauss_l1l2l3); + GetParameterValue(&drag, &k[0],gauss_l1l2l3); #ifdef _DEBUG_ printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag); @@ -1796,5 +1736,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::GradjB" -void Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){ +void Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,void* vinputs,int analysis_type){ int i; @@ -1802,6 +1742,7 @@ /* node data: */ const int numgrids=3; - const int numdofs=2*numgrids; - const int NDOF2=2; + const int NDOF1=1; + const int NDOF2=1; + const int numdofs=NDOF1*numgrids; double xyz_list[numgrids][3]; int doflist[numdofs]; @@ -1814,7 +1755,4 @@ double adjx_list[numgrids]; double adjy_list[numgrids]; - - double* thickness_param=NULL; - double thickness_list[numgrids]; /* gaussian points: */ @@ -1850,5 +1788,10 @@ double lambda,mu; double thickness; - + int dofs[1]={0}; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /* Get node coordinates and dof list: */ @@ -1860,5 +1803,5 @@ /* recover input parameters: */ - thickness_param=ParameterInputsRecover(inputs,"thickness"); + inputs->Recover("thickness",&h[0],1,dofs,numgrids,(void**)nodes); /*Initialize parameter lists: */ @@ -1870,10 +1813,4 @@ adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]]; adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]]; - if(thickness_param){ - thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]]; - } - else{ - thickness_list[i]=h[i]; - } } @@ -1940,5 +1877,5 @@ #undef __FUNCT__ #define __FUNCT__ "Tria::Misfit" -double Tria::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){ +double Tria::Misfit(double* velocity,double* obs_velocity,void* vinputs,int analysis_type){ int i; @@ -1979,6 +1916,10 @@ double scalex=1; double scaley=1; - double* fit_param=NULL; double fit=-1; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; /* Get node coordinates and dof list: */ @@ -1987,9 +1928,5 @@ /* Recover input data: */ - fit_param=ParameterInputsRecover(inputs,"fit"); - if(fit_param){ - fit=*fit_param; - } - else ErrorException(__FUNCT__," missing fit input parameter"); + if(!inputs->Recover("fit",&fit)) throw ErrorException(__FUNCT__," missing fit input parameter"); /*Initialize velocities: */ @@ -2102,7 +2039,92 @@ + +#undef __FUNCT__ +#define __FUNCT__ "Tria::CreateKMatrixDiagnosticBaseVert" +void Tria::CreateKMatrixDiagnosticBaseVert(Mat Kgg,void* vinputs,int analysis_type){ + + int i,j; + + /* node data: */ + const int numgrids=3; + const int NDOF1=1; + const int numdofs=NDOF1*numgrids; + double xyz_list[numgrids][3]; + int doflist[numdofs]; + int numberofdofspernode; + + /* gaussian points: */ + int num_gauss,ig; + double* first_gauss_area_coord = NULL; + double* second_gauss_area_coord = NULL; + double* third_gauss_area_coord = NULL; + double* gauss_weights = NULL; + double gauss_weight; + double gauss_l1l2l3[3]; + + + /* matrices: */ + double L[1][3]; + double DL_scalar; + + double Ke_gg[3][3]; //stiffness matrix + double Ke_gg_gaussian[3][3]; //stiffness matrix evaluated at the gaussian point. + double Jdet; + + ParameterInputs* inputs=NULL; + + /*recover pointers: */ + inputs=(ParameterInputs*)vinputs; + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Set Ke_gg to 0: */ + for(i=0;iRecover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity"); + inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes); /* Get node coordinates and dof list: */ @@ -2281,18 +2315,4 @@ for(i=0;iFindParam((void*)&mincontrolconstraint,"mincontrolconstraint"); femmodel.parameters->FindParam((void*)&maxcontrolconstraint,"maxcontrolconstraint"); - gsize=femmodel.nodes->NumberOfDofs(); - + femmodel.parameters->FindParam((void*)&gsize,"gsize"); femmodel.parameters->FindParam((void*)&p_g,"p_g"); femmodel.parameters->FindParam((void*)&u_g_obs,"u_g_obs"); @@ -96,5 +98,21 @@ /*Initialize inputs:*/ - inputs=NewParameterInputs(); + + femmodel.parameters->FindParam((void*)&u_g_initial,"u_g"); + femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + + inputs=new ParameterInputs; + inputs->Add("velocity",u_g_initial,3,numberofnodes); + + /*erase useless parameters: */ + param=(Param*)femmodel.parameters->FindParamObject("p_g"); + femmodel.parameters->DeleteObject((Object*)param); + + param=(Param*)femmodel.parameters->FindParamObject("u_g_obs"); + femmodel.parameters->DeleteObject((Object*)param); + + param=(Param*)femmodel.parameters->FindParamObject("u_g"); + femmodel.parameters->DeleteObject((Object*)param); + /*Start looping: */ @@ -102,6 +120,6 @@ _printf_("\n%s%i%s%i\n"," control method step ",n+1,"/",nsteps); - ParameterInputsAddFromMat(inputs,p_g,gsize,control_type); - ParameterInputsAddFromDouble(inputs,fit[n],"fit"); + inputs->Add(control_type,p_g,2,numberofnodes); + inputs->Add("fit",fit[n]); _printf_("%s\n"," computing gradJ..."); @@ -135,6 +153,6 @@ if (((n+1)%5)==0){ _printf_("%s\n"," saving temporary results..."); - ParameterInputsAddFromMat(inputs,p_g,gsize,control_type); - ParameterInputsAddFromDouble(inputs,fit[n],"fit"); + inputs->Add(control_type,p_g,2,numberofnodes); + inputs->Add("fit",fit[n]); diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); OutputControl(u_g,p_g,J,nsteps,femmodel.partition,outputfilename,femmodel.nodesets); @@ -150,6 +168,6 @@ /*Write results to disk: */ _printf_("%s\n"," preparing final velocity solution..."); - ParameterInputsAddFromMat(inputs,p_g,gsize,control_type); - ParameterInputsAddFromDouble(inputs,fit[n],"fit"); + inputs->Add(control_type,p_g,2,numberofnodes); + inputs->Add("fit",fit[n]); diagnostic_core_nonlinear(&u_g,NULL,NULL,inputs,&femmodel); Index: /issm/trunk/src/c/parallel/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic.cpp (revision 245) +++ /issm/trunk/src/c/parallel/diagnostic.cpp (revision 246) @@ -24,4 +24,5 @@ char* lockname=NULL; char* analysis_type="diagnostic_horiz"; + int numberofnodes; /*Intermediary: */ @@ -32,5 +33,4 @@ double* u_g_initial=NULL; - int gsize; Param* param=NULL; @@ -60,9 +60,10 @@ _printf_("initialize inputs:\n"); - inputs=NewParameterInputs(); - femmodel.parameters->FindParam((void*)&u_g,"u_g"); - femmodel.parameters->FindParam((void*)&gsize,"gsize"); - ParameterInputsAddFromMat(inputs,u_g_initial,gsize,"velocity"); + femmodel.parameters->FindParam((void*)&u_g_initial,"u_g"); + femmodel.parameters->FindParam((void*)&numberofnodes,"numberofnodes"); + inputs=new ParameterInputs; + inputs->Add("velocity",u_g_initial,3,numberofnodes); + param=(Param*)femmodel.parameters->FindParamObject("u_g"); femmodel.parameters->DeleteObject((Object*)param); Index: /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 245) +++ /issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp (revision 246) @@ -31,4 +31,5 @@ int num_unstable_constraints; int count; + int numberofnodes; Vec dug; @@ -47,4 +48,5 @@ fem->parameters->FindParam((void*)&connectivity,"connectivity"); fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode"); + fem->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); fem->parameters->FindParam((void*)&solver_string,"solverstring"); fem->parameters->FindParam((void*)&eps_rel,"eps_rel"); @@ -64,6 +66,6 @@ /*Set input parameters: */ - if(old_ug)ParameterInputsAddFromVec(inputs,old_ug,"old_velocity"); - if(ug)ParameterInputsAddFromVec(inputs,ug,"velocity"); + if(old_ug)inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes); + if(ug)inputs->Add("velocity",ug,numberofdofspernode,numberofnodes); //save pointer to old velocity @@ -112,6 +114,6 @@ //Deal with penalty loads if (debug) _printf_(" penalty constraints\n"); - ParameterInputsAddFromVec(inputs,ug,"velocity"); - ParameterInputsAddFromVec(inputs,old_ug,"old_velocity"); + inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes); + inputs->Add("velocity",ug,numberofdofspernode,numberofnodes); PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,loads,fem->materials,inputs,analysis_type); @@ -157,6 +159,6 @@ /*Set input parameters: */ - ParameterInputsAddFromVec(inputs,ug,"velocity"); - ParameterInputsAddFromVec(inputs,old_ug,"old_velocity"); + inputs->Add("old_velocity",old_ug,numberofdofspernode,numberofnodes); + inputs->Add("velocity",ug,numberofdofspernode,numberofnodes); kflag=1; pflag=0; //stiffness generation only Index: /issm/trunk/src/c/parallel/objectivefunctionC.cpp =================================================================== --- /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 245) +++ /issm/trunk/src/c/parallel/objectivefunctionC.cpp (revision 246) @@ -41,4 +41,5 @@ Vec u_g=NULL; double* u_g_double=NULL; + int numberofnodes; @@ -58,4 +59,5 @@ femmodel->parameters->FindParam((void*)&fit,"fit"); femmodel->parameters->FindParam((void*)&analysis_type,"analysis_type"); + femmodel->parameters->FindParam((void*)&numberofnodes,"numberofnodes"); /*First copy p_g so we don't modify it: */ @@ -70,5 +72,5 @@ /*Add new parameter to inputs: */ - ParameterInputsAddFromMat(inputs,p_g_copy,gsize,control_type); + inputs->Add(control_type,p_g_copy,2,numberofnodes); //Run diagnostic with updated parameters. @@ -77,5 +79,5 @@ //Compute misfit for this velocity field. - ParameterInputsAddFromDouble(inputs,fit[n],"fit"); + inputs->Add("fit",fit[n]); Misfitx( &J, femmodel->elements,femmodel->nodes, femmodel->loads, femmodel->materials, u_g_double,u_g_obs, inputs,analysis_type); Index: /issm/trunk/src/m/classes/@input/display.m =================================================================== --- /issm/trunk/src/m/classes/@input/display.m (revision 246) +++ /issm/trunk/src/m/classes/@input/display.m (revision 246) @@ -0,0 +1,31 @@ +function display(input) +%DISPLAY - displays the fields of an input object +% +% echo function for 'input' class + +STRING=0; +INTEGER=1; +DOUBLE=2; +DOUBLEVEC=3; + +if isempty(input.name), + %do nothing + return; +end + +%disp(sprintf('\n%s = \n',inputname(1))); + +switch input.type, + case STRING, + disp(sprintf(' %s(string): %s',input.name,input.string)); + case INTEGER, + disp(sprintf(' %s(integer): %i',input.name,input.integer)); + case DOUBLE, + disp(sprintf(' %s(double): %i',input.name,input.double)); + case DOUBLEVEC, + disp(sprintf(' %s(vector): ndof %i numberofnodes %i ',input.name,input.ndof,input.numberofnodes)); + for i=1:length(input.vector), + disp(sprintf(' %g',input.vector(i))); + end +end + Index: /issm/trunk/src/m/classes/@input/get.m =================================================================== --- /issm/trunk/src/m/classes/@input/get.m (revision 246) +++ /issm/trunk/src/m/classes/@input/get.m (revision 246) @@ -0,0 +1,47 @@ +function field=get(input,varargin) +%GET - returns field from input +% usage: +% field=get(input); +% field=get(input,doflist); %dof flag list for doublevec vectors. +% +% ex: velocity=get(input) +% or: velocity=get(input,[1 1 0]); %get only vx and vy +% or: velocity=get(input,[1 1 1]); %get vx,vy and vz +% + + +STRING=0; +INTEGER=1; +DOUBLE=2; +DOUBLEVEC=3; + +field=[]; %in case nothing is returned +if input.type==STRING, + field=input.string; +end + +if input.type==INTEGER, + field=input.integer; +end + +if input.type==DOUBLE, + field=input.double; +end + +if input.type==DOUBLEVEC, + if nargin==1, + field=input.double; + else + doflist=varargin{1}; + numdofs=numel(find(doflist)); + field=zeros(numdofs*input.numberofnodes,1); + + count=0; + for i=1:length(doflist), + if doflist(i), + count=count+1; + field(count:numdofs:end)=input.vector(i:input.ndof:end); + end + end + end +end Index: /issm/trunk/src/m/classes/@input/getname.m =================================================================== --- /issm/trunk/src/m/classes/@input/getname.m (revision 246) +++ /issm/trunk/src/m/classes/@input/getname.m (revision 246) @@ -0,0 +1,6 @@ +function name=getname(input) +%GETNAME - returns name of an input +% + +name=input.name; + Index: /issm/trunk/src/m/classes/@input/input.m =================================================================== --- /issm/trunk/src/m/classes/@input/input.m (revision 246) +++ /issm/trunk/src/m/classes/@input/input.m (revision 246) @@ -0,0 +1,135 @@ +function input = input(varargin) +%INPUTS - contructor for input objects +% +% Usage: +% input = input() +% input = input(name,field,field_type) + +STRING=0; +INTEGER=1; +DOUBLE=2; +DOUBLEVEC=3; + +switch nargin +case 0 + input.name=''; + input.type=DOUBLE; + + %placeholders: + input.string=''; + input.integer=NaN; + input.double=NaN; + input.vector=[]; + input.ndof=NaN; + input.numberofnodes=NaN; + + input=class(input,'input'); + + +case 3 + name=varargin{1}; + field=varargin{2}; + field_type=varargin{3}; + + if ~ischar(name), + error('in input(name,field,field_type) constructor call, name is a string'); + end + + if ~ischar(field_type), + error('in input(name,field,field_type) constructor call, field_type is a string: ''string'',''integer'',''double'' or ''doublevec'' '); + end + + input.name=''; + input.type=DOUBLE; + + %placeholders: + input.string=''; + input.integer=NaN; + input.double=NaN; + input.vector=[]; + input.ndof=NaN; + input.numberofnodes=NaN; + + %fill fields + input.name=name; + if strcmpi(field_type,'string'), + input.type=STRING; + input.string=field; + end + + if strcmpi(field_type,'integer'), + input.type=INTEGER; + input.integer=field; + end + + if strcmpi(field_type,'double'), + input.type=DOUBLE; + input.double=field; + end + + if strcmpi(field_type,'doublevec'), + input.type=DOUBLEVEC; + input.vector=field; + input.ndof=1; + input.numberofnodes=length(input.vector); + end + input=class(input,'input'); + +case 5 + name=varargin{1}; + field=varargin{2}; + field_type=varargin{3}; + ndof=varargin{4}; + numberofnodes=varargin{5}; + + if ~ischar(name), + error('in input(name,field,field_type) constructor call, name is a string'); + end + + if ~ischar(field_type), + error('in input(name,field,field_type) constructor call, field_type is a string: ''string'',''integer'',''double'' or ''doublevec'' '); + end + + input.name=''; + input.type=DOUBLE; + + %placeholders: + input.string=''; + input.integer=NaN; + input.double=NaN; + input.vector=[]; + input.ndof=NaN; + input.numberofnodes=NaN; + + %fill fields + input.name=name; + if strcmpi(field_type,'string'), + input.type=STRING; + input.string=field; + end + + if strcmpi(field_type,'integer'), + input.type=INTEGER; + input.integer=field; + end + + if strcmpi(field_type,'double'), + input.type=DOUBLE; + input.double=field; + end + + if strcmpi(field_type,'doublevec'), + input.type=DOUBLEVEC; + input.vector=field; + input.ndof=ndof; + input.numberofnodes=numberofnodes; + if (input.ndof*input.numberofnodes)~=length(input.vector), + error(['length of vector ''',input.name,''' should be ndof*numberofnodes=' num2str(ndof*numberofnodes)]); + end + end + input=class(input,'input'); + + +otherwise + error('input constructor error message: 0 of 1 argument only in input.'); +end Index: /issm/trunk/src/m/classes/@parameterinputs/add.m =================================================================== --- /issm/trunk/src/m/classes/@parameterinputs/add.m (revision 246) +++ /issm/trunk/src/m/classes/@parameterinputs/add.m (revision 246) @@ -0,0 +1,34 @@ +function parameterinputs = add(parameterinputs,field_name,field,field_type,varargin) +%ADD - add an input to parameter inputs +% +% Usage: +% parameterinputs = add(parameterinputs,field_name,field,field_type) + +found=0; + +if length(parameterinputs)==0, + parameterinputs + parameterinputs(end+1).input=input(field_name,field,field_type,varargin); +else + for i=1:length(parameterinputs), + + if strcmpi(field_name,getname(parameterinputs(i).input)), + %replace this field by new field + parameterinputs(i).input=input(field_name,field,field_type,varargin{:}); + found=1; + break; + end + end + + + %did we find the input name? + if found, + %we are done + else + %add at the end + parameterinputs(end+1).input=input(field_name,field,field_type,varargin{:}); + end +end +if ~isa(parameterinputs,'parameterinputs'), + parameterinputs=class(parameterinputs,'parameterinputs'); +end Index: /issm/trunk/src/m/classes/@parameterinputs/display.m =================================================================== --- /issm/trunk/src/m/classes/@parameterinputs/display.m (revision 246) +++ /issm/trunk/src/m/classes/@parameterinputs/display.m (revision 246) @@ -0,0 +1,13 @@ +function display(parameterinputs) +%DISPLAY - displays the fields of an parameterinputs object +% +% echo function for 'parameterinputs' class +disp('parameter parameterinputs echo'); + +if length(parameterinputs)==0, + disp(' empty'); +else + for i=1:length(parameterinputs), + parameterinputs(i).input + end +end Index: /issm/trunk/src/m/classes/@parameterinputs/get.m =================================================================== --- /issm/trunk/src/m/classes/@parameterinputs/get.m (revision 246) +++ /issm/trunk/src/m/classes/@parameterinputs/get.m (revision 246) @@ -0,0 +1,21 @@ +function field = get(parameterinputs,field_name,varargin) +%ADD - get a parameter input with name "field_name". +% +% Usage: +% parameterinputs = get(parameterinputs,field_name,varargin) +% optional argument: list of dofs being kept. +% ex: +% parameterinputs=get(parameterinputs,'velocity',[1 1 0]); %only keep vx,vy +% parameterinputs=get(parameterinputs,'velocity',[1 1 1]); %keep vx,vy,vz + +field=[]; %in case we find nothing + +for i=1:length(parameterinputs), + + if strcmpi(field_name,getname(parameterinputs(i).input)), + + %found the field, call the get method + field=get(parameterinputs(i).input,varargin{:}); + return; + end +end Index: /issm/trunk/src/m/classes/@parameterinputs/parameterinputs.m =================================================================== --- /issm/trunk/src/m/classes/@parameterinputs/parameterinputs.m (revision 246) +++ /issm/trunk/src/m/classes/@parameterinputs/parameterinputs.m (revision 246) @@ -0,0 +1,21 @@ +function parameterinputs = parameterinputs(varargin) +%INPUTS - contructor for parameterinputs objects +% +% Usage: +% parameterinputs = parameterinputs(varargin) + +switch nargin +case 0 + parameterinputs.input=input; + parameterinputs=class(parameterinputs,'parameterinputs'); + +case 1 + %If single argument of class parameterinputs, we have a copy constructor. + if (isa(varargin{1},'parameterinputs')) + parameterinputs = varargin{1}; + else + error('parameterinputs constructor error message: copy constructor called on a non ''parameterinputs'' class object'); + end +otherwise + error('parameterinputs constructor error message: 0 of 1 argument only in input.'); +end Index: /issm/trunk/src/m/solutions/cielo/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 245) +++ /issm/trunk/src/m/solutions/cielo/diagnostic.m (revision 246) @@ -19,11 +19,6 @@ %initialize inputs - if strcmp(md.analysis_type,'diagnostic_stokes'), - inputs.velocity=m_dh.parameters.u_g; - else - inputs.velocity=zeros(2*m_dh.parameters.numberofnodes,1); - inputs.velocity(1:2:end)=m_dh.parameters.u_g(1:3:end); - inputs.velocity(2:2:end)=m_dh.parameters.u_g(2:3:end); - end + inputs=parameterinputs; + inputs=add(inputs,'velocity',m_dh.parameters.u_g,'doublevec',3,m_dh.parameters.numberofnodes); %Get horizontal solution. Index: /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 245) +++ /issm/trunk/src/m/solutions/cielo/diagnostic_core_nonlinear.m (revision 246) @@ -12,10 +12,5 @@ converged=0; - [velocity_param velocity_is_present]=recover_input(inputs,'velocity'); - if velocity_is_present - soln(count).u_g=velocity_param; - else - soln(count).u_g=[]; - end + soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node soln(count).u_f=[]; @@ -27,14 +22,14 @@ %add velocity to inputs - inputs.velocity=full(soln(count).u_g); + inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); if count>1, - inputs.old_velocity=full(soln(count-1).u_g); + inputs=add(inputs,'old_velocity',soln(count-1).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); else - inputs.old_velocity=full(soln(count).u_g); + inputs=add(inputs,'old_velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); end %Update inputs in datasets [m.elements,m.nodes, loads,m.materials]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,inputs); - + %system matrices [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs); @@ -67,5 +62,5 @@ %Deal with penalty loads - inputs.velocity=soln(count).u_g; + inputs=add(inputs,'velocity',soln(count).u_g,'vector',m.parameters.numberofdofspernode,m.parameters.numberofnodes); %penalty constraints @@ -73,33 +68,31 @@ % Figure out if convergence is reached. - if(count>=3 | velocity_is_present), - dug=soln(count).u_g-soln(count-1).u_g; - nduinf=norm(dug,inf)*m.parameters.yts; - ndu=norm(dug,2); - nu=norm(soln(count-1).u_g,2); + dug=soln(count).u_g-soln(count-1).u_g; + nduinf=norm(dug,inf)*m.parameters.yts; + ndu=norm(dug,2); + nu=norm(soln(count-1).u_g,2); - %Relative criterion - if (ndu/nu<=m.parameters.eps_rel), - if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: norm(du)/norm(u)=',ndu/nu,' < ',m.parameters.eps_rel)); end - converged=1; + %Relative criterion + if (ndu/nu<=m.parameters.eps_rel), + if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: norm(du)/norm(u)=',ndu/nu,' < ',m.parameters.eps_rel)); end + converged=1; + else + if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel)); end + converged=0; + end + + %Absolute criterion + if ~isnan(m.parameters.eps_abs), + if (nduinf<=m.parameters.eps_abs), + if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: max(du)=',nduinf,' < ',m.parameters.eps_abs)); end else - if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: norm(du)/norm(u)=',ndu/nu,' > ',m.parameters.eps_rel)); end - converged=0; + if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: max(du)=',nduinf,' > ',m.parameters.eps_abs)); end + converged=0; end + end - %Absolute criterion - if ~isnan(m.parameters.eps_abs), - if (nduinf<=m.parameters.eps_abs), - if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: max(du)=',nduinf,' < ',m.parameters.eps_abs)); end - else - if m.parameters.debug, disp(sprintf('%s %g %s %g',' convergence criterion: max(du)=',nduinf,' > ',m.parameters.eps_abs)); end - converged=0; - end - end - - %rift convergence criterion - if ~constraints_converged, - converged=0; - end + %rift convergence criterion + if ~constraints_converged, + converged=0; end @@ -118,5 +111,5 @@ if nout==2, - inputs.velocity=full(soln(count).u_g); + inputs=add(inputs,'velocity',soln(count).u_g,'vector',m.parameters.numberofdofspernode,m.parameters.numberofnodes); m.parameters.kflag=1; m.parameters.pflag=0; [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs); Index: /issm/trunk/src/mex/Du/Du.cpp =================================================================== --- /issm/trunk/src/mex/Du/Du.cpp (revision 245) +++ /issm/trunk/src/mex/Du/Du.cpp (revision 246) @@ -40,5 +40,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Call core code: */ @@ -56,5 +57,5 @@ xfree((void**)&u_g_obs); VecFree(&du_g); - DeleteParameterInputs(&inputs); + delete inputs; /*end module: */ Index: /issm/trunk/src/mex/Gradj/Gradj.cpp =================================================================== --- /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 245) +++ /issm/trunk/src/mex/Gradj/Gradj.cpp (revision 246) @@ -42,5 +42,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Call core code: */ @@ -58,5 +59,5 @@ xfree((void**)&lambda_g); xfree((void**)&control_type); - DeleteParameterInputs(&inputs); + delete inputs; VecFree(&grad_g); Index: /issm/trunk/src/mex/Misfit/Misfit.cpp =================================================================== --- /issm/trunk/src/mex/Misfit/Misfit.cpp (revision 245) +++ /issm/trunk/src/mex/Misfit/Misfit.cpp (revision 246) @@ -39,5 +39,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Call core code: */ @@ -54,5 +55,5 @@ xfree((void**)&u_g); xfree((void**)&u_g_obs); - DeleteParameterInputs(&inputs); + delete inputs /*end module: */ Index: /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp =================================================================== --- /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 245) +++ /issm/trunk/src/mex/PenaltyConstraints/PenaltyConstraints.cpp (revision 246) @@ -38,5 +38,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Generate internal degree of freedom numbers: */ @@ -54,5 +55,5 @@ delete loads; delete materials; - DeleteParameterInputs(&inputs); + delete inputs; /*end module: */ Index: /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp =================================================================== --- /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (revision 245) +++ /issm/trunk/src/mex/PenaltySystemMatrices/PenaltySystemMatrices.cpp (revision 246) @@ -40,5 +40,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Generate stiffnesses from penalties: */ @@ -54,5 +55,5 @@ delete loads; delete materials; - DeleteParameterInputs(&inputs); + delete inputs; MatFree(&Kgg); VecFree(&pg); Index: /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp =================================================================== --- /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp (revision 245) +++ /issm/trunk/src/mex/SystemMatrices/SystemMatrices.cpp (revision 246) @@ -44,5 +44,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Generate internal degree of freedom numbers: */ @@ -59,5 +60,5 @@ delete loads; delete materials; - DeleteParameterInputs(&inputs); + delete inputs; MatFree(&Kgg); VecFree(&pg); Index: /issm/trunk/src/mex/UpdateFromInputs/UpdateFromInputs.cpp =================================================================== --- /issm/trunk/src/mex/UpdateFromInputs/UpdateFromInputs.cpp (revision 245) +++ /issm/trunk/src/mex/UpdateFromInputs/UpdateFromInputs.cpp (revision 246) @@ -30,5 +30,6 @@ /*Fetch inputs: */ - ParameterInputsInit(&inputs,INPUTS); + inputs=new ParameterInputs; + inputs->Init(INPUTS); /*!Generate internal degree of freedom numbers: */ @@ -46,5 +47,5 @@ delete loads; delete materials; - DeleteParameterInputs(&inputs); + delete inputs; /*end module: */